00001 #ifndef _MY_DealBase_
00002
00003 #define _MY_DealBase_
00004
00005
00006
00007 #include "globals.h"
00008 #include <grid/tria.h>
00009 #include <numerics/vectors.h>
00010 #include <numerics/matrices.h>
00011 #include "sfunctions.h"
00012 #include "gnuplotanim.h"
00013 #include <lac/vector.h>
00014 #include <lac/full_matrix.h>
00015 #include <lac/sparse_matrix.h>
00016 #include <lac/solver_cg.h>
00017 #include <lac/vector_memory.h>
00018 #include <lac/precondition.h>
00019 #include <iostream>
00025 class DealBase
00026 {
00027 private:
00028 Triangulation<DIM> *m_pTriangulation;
00029 Cell END_CELL;
00030 protected:
00031 void setTriangulation(Triangulation<DIM>& trig);
00032 int n_cells;
00033
00034 public:
00035 static const int INVALID_INDEX;
00036
00037
00038 DealBase(){m_pTriangulation = NULL;n_cells=-1;}
00039 DealBase(const Triangulation<DIM> &tria){setTriangulation((Triangulation<DIM>&) tria);}
00040 ~DealBase(){}
00041 Triangulation<DIM>& getTriangulation() const {assert(m_pTriangulation);return *m_pTriangulation;}
00042
00047
00048 Cell getCell(int index);
00049 Cell endCell() const {return getTriangulation().end();}
00050 Cell beginCell() const {return getTriangulation().begin_active();}
00051 Face beginFace(){return getTriangulation().begin_active_face();}
00052 Face endFace(){return getTriangulation().end_face();}
00053 void getFaceBarycenter(Face &face,Point3D &p);
00054 Point3D getFaceBarycenter(Face &face);
00055 double getCentralValue(const Cell &cell,const VecDouble &values) const;
00056 double getFaceValue(const Face &face,const VecDouble &values);
00057 double getFaceValue(const Cell &cell,FaceDirection3D dir,const VecDouble &values);
00058 unsigned getFaceValueIndex(const Cell &cell,FaceDirection3D dir);
00059 bool getFaceFlag(const Face &face,const VecBool& values);
00060 Cell getCellAtPoint(const Point<DIM> &p);
00061 double getNeighborValue(const Cell &cell,CellDirection3D dir,const VecDouble &cValues);
00062 double getVerticeValue(const Cell &cell,VertexDirection3D vertex,const VecDouble &values);
00063 double getVerticeValue(unsigned vIndex,const VecDouble &values){return values(vIndex);}
00064 unsigned getCentralPointIndex(Cell &cell);
00065 void getVerticePoint(unsigned vIndex,Point<DIM>& X) const;
00066 Cell getEndCell(){return END_CELL;}
00067 Cell getAdjCell(Cell &cell,CellDirection3D dir);
00068 Cell getAdjCell(Cell &cell,CellDirection3D dir1,CellDirection3D dir2);
00069 bool getCellFlag(const Cell &cell,const VecBool &cVFlags) const;
00070 bool getVerticeFlag(const Cell &cell,VertexDirection3D dir,const VecBool &vVFlags) const;
00071
00072
00073 void setCentralValue(const Cell &cell,double dd,VecDouble &values);
00074 void setCentralValue(DoFCell &cell,double dd,VecDouble &values);
00075 void setFaceValue(const Face &face,double dd,VecDouble &values);
00076 void setFaceValue(const Cell &cell,FaceDirection3D dir,double dd,VecDouble &values);
00077 void setFaceFlag(const Face &face,bool b,VecBool &values);
00078 void setCellFlag(const Cell &cell,const bool flag,VecBool &cVFlags) const;
00079 void setVerticeFlag(const Cell &cell,VertexDirection3D dir,const bool flag,VecBool &vVFlags) const;
00082 double getL2NormAtCells(const VecDouble &cValues) const;
00083 double getIntegralAtCells(const VecDouble &cValues) const;
00084
00085
00086
00087 unsigned numCells() const {return n_cells;}
00088 unsigned numFaces() const {return getTriangulation().n_active_faces();}
00089 unsigned numVertices() const {return getTriangulation().n_vertices();}
00090 double faceArea(Face face);
00091
00096 int getAdjCellIndex(unsigned iCell,CellDirection3D dir);
00097 double getAdjCellValue(unsigned iCell,CellDirection3D dir,const VecDouble &v);
00098
00099 int getFaceIndex(unsigned iCell,FaceDirection3D dir);
00100 double getFaceValue(unsigned iCell,FaceDirection3D dir, const VecDouble &v);
00101
00102 bool isCellAtBoundary(unsigned iCell);
00103 bool hasAdjCell(unsigned iCell, CellDirection3D dir);
00113 void setCentralValuesFromFunction(Function3D &f,VecDouble &cValues,unsigned component=0);
00114 void setCentralValuesInFunctionDomain(Function3D &f,MapIntDouble &values,unsigned component=0);
00115 void setVerticesValuesFromFunction(Function<DIM> &f,VecDouble &vValues,unsigned deg=0);
00116 void setFacesValuesFromFunction(Function3D &f,VecDouble &fValues,unsigned deg=0,bool bOnlyPrescribedBoundary=false);
00117 void setFacesValuesFromFunction(Function3D &f,Matrix &MValues,bool bOnlyPrescribedBoundary=false);
00118 void tagFacesInDomain(Function3D &f,VecTag &tags,char value);
00119 void tagFacesInDomain(Function3D &f,VecBool &tags);
00120
00121 void setVerticeValue(Cell &cell,VertexDirection3D vertex,double dd,VecDouble &values);
00122 void setMaterialIdOfPrescribedBoundaryCells(GeneralFunctionInterface &f,unsigned material_id);
00123 void setMaterialIdFromFunctionDomain(GeneralFunctionInterface &f,unsigned material_id);
00124 void getCellsInFunctionDomain(Function3D &f,VecIndex &vec);
00125
00133 void printMatrixMaple(FullMatrix<double> &M);
00134 void printMatrixMaple(SparseMatrix<double> &M);
00135 void printVertices();
00136 void printCells();
00137 void printFaces();
00138 void printPoint(const Point<DIM> &P,std::ostream &out);
00139 void printVector(const VecDouble &vSol);
00140 void printPoint(const Point<DIM> &p);
00141 void printCentralValues(const VecDouble& cValues);
00142 void printCentralValues(const VecDouble& cValues1,const VecDouble& cValues2);
00143 void printCentralValues(const VecDouble& cValues1,const VecDouble& cValues2,const VecDouble& cValues3,const VecDouble& cValues4);
00144 void printFacesValuesAtBoundary(const VecDouble& fValues,const VecDouble& fValues2);
00145 void printFacesValuesAtBoundary(const VecDouble& fValues1,const VecDouble& fValues2,const VecDouble& fValues3,const VecDouble& fValues4);
00146 void printFacesValuesAtBoundary(const VecDouble& fValues);
00147 static void printTriangulation(Triangulation<3> &tria,std::ostream &out = std::cout);
00149 void printCell(Cell cell,std::ostream &out=std::cout);
00150 void printCells(VecIndex &vec,std::ostream &out=std::cout);
00151
00152 bool isFirstCell(Cell& cell){return cell == getTriangulation().begin_active();}
00153 bool isValid(Cell& cell){return (cell.state() == IteratorState::valid);}
00154 bool hasNeighbor(Cell &cell,CellDirection3D dir){return (cell->neighbor_index(dir) != DealBase::INVALID_INDEX);}
00155 void flagPrescribedFaces(Function3D &f,VecBool &vFlags);
00156 void setPrescribedBoundaryId(Function3D &f,unsigned boundary_id);
00157 double getDX (Cell &cell) {return cell->vertex(VERTEX_100)(0)-cell->vertex(VERTEX_000)(0);}
00158 double getDZ (Cell &cell) {return cell->vertex(VERTEX_001)(2)-cell->vertex(VERTEX_000)(2);}
00159 double getDY (Cell &cell) {return cell->vertex(VERTEX_010)(1)-cell->vertex(VERTEX_000)(1);}
00160
00161 };
00162
00163 #endif
