diff --git a/AMDiS/src/AMDiS.cc b/AMDiS/src/AMDiS.cc index 8d65eff5f2a590dc29eddd116a959ec62017bd2b..998ea2210ccb343986c2735d443dd02e7359852c 100644 --- a/AMDiS/src/AMDiS.cc +++ b/AMDiS/src/AMDiS.cc @@ -125,6 +125,7 @@ namespace AMDiS { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS MeshDistributor::globalMeshDistributor->exitParallelization(); + delete MeshDistributor::globalMeshDistributor; #ifdef HAVE_PARALLEL_MTL4 if (mtl_environment) delete mtl_environment; diff --git a/AMDiS/src/ComponentTraverseInfo.h b/AMDiS/src/ComponentTraverseInfo.h index 5e55fccf23de619ea3a3ea609f8b6422705aae4d..63fc773daee41710febb38a1f0e7e375e375fc46 100644 --- a/AMDiS/src/ComponentTraverseInfo.h +++ b/AMDiS/src/ComponentTraverseInfo.h @@ -177,11 +177,9 @@ namespace AMDiS { return vectorComponents[row].getStatus(); } - /** \brief - * Returns true, if for the given matrix component the row and the col FE spaces - * are equal. Note that this is also the case, if there is another aux FE space, - * that may be different from both, the row and the col FE spaces. - */ + /// Returns true, if for the given matrix component the row and the col FE spaces + /// are equal. Note that this is also the case, if there is another aux FE space, + /// that may be different from both, the row and the col FE spaces. inline bool eqSpaces(int row, int col) const { int status = matrixComponents[row][col].getStatus(); diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh index 48aeaf7a1b9aa54cac4e38b9b87a3b5c44524815..39ecb551407ed2488f7969748c491d35ebc57ef2 100644 --- a/AMDiS/src/DOFVector.hh +++ b/AMDiS/src/DOFVector.hh @@ -1495,11 +1495,12 @@ namespace AMDiS { for (int iq = 0; iq < nPoints; iq++) { nullify(grd1); - for (int j = 0; j < nBasFcts; j++) // #BasisFunctions + for (int j = 0; j < nBasFcts; j++) { // #BasisFunctions for (int k = 0; k < parts; k++) // #edges (2d) or #faces (3d) grd1[k] += quadFast->getGradient(iq, j, k) * localVec[j]; + } - for (int l=0; l < dow; l++) { + for (int l = 0; l < dow; l++) { nullify(grdAtQPs[iq][l]); for (int k = 0; k < parts; k++) grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k]; diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h index f09ebc85942dcdbc7531a9ba4b7731e365643a0b..241e765ade7f2cd01b0016932c9ee346d098a9a5 100644 --- a/AMDiS/src/ElInfo.h +++ b/AMDiS/src/ElInfo.h @@ -51,20 +51,16 @@ namespace AMDiS { /// Protected constructor. Avoids instatiation of the basis class ElInfo(); - /** \brief - * Protected constructor. Avoids instatiation of the basis class. - * \param mesh pointer to the corresponding mesh. - */ + /// Protected constructor. Avoids instatiation of the basis class. + /// \param mesh pointer to the corresponding mesh. ElInfo(Mesh *mesh); public: /// Virtual destructor because ElInfo is pure virtual. virtual ~ElInfo(); - /** \brief - * Assignement operator. - * \param rhs right hand side. - */ + /// Assignement operator. + /// \param rhs right hand side. ElInfo& operator=(const ElInfo& rhs) { mesh = rhs.mesh; @@ -137,37 +133,29 @@ namespace AMDiS { return iChild; } - /** \brief - * Get ElInfo's \ref coord[i]. This is a WorldVector filled with the - * coordinates of the i-th vertex of element \ref el. - */ + /// Get ElInfo's \ref coord[i]. This is a WorldVector filled with + /// the coordinates of the i-th vertex of element \ref el. inline WorldVector& getCoord(int i) { return coord[i]; } - /** \brief - * Get ElInfo's \ref coord[i]. This is a WorldVector filled with the - * coordinates of the i-th vertex of element \ref el. - */ + /// Get ElInfo's \ref coord[i]. This is a WorldVector filled with the + /// coordinates of the i-th vertex of element \ref el. inline const WorldVector& getCoord(int i) const { return coord[i]; } - /** \brief - * Get ElInfo's \ref coord. This is a FixVec > filled with the - * coordinates of the all vertice of element \ref el. - */ + /// Get ElInfo's \ref coord. This is a FixVec > filled + /// with the coordinates of the all vertice of element \ref el. inline FixVec, VERTEX>& getCoords() { return coord; } - /** \brief - * Get ElInfo's \ref coord. This is a FixVec > filled with the - * coordinates of the all vertice of element \ref el. - */ + /// Get ElInfo's \ref coord. This is a FixVec > filled + /// with the coordinates of the all vertice of element \ref el. inline const FixVec, VERTEX>& getCoords() const { return coord; @@ -370,83 +358,67 @@ namespace AMDiS { /** \} */ - /** \brief - * Returns the absolute value of the determinant of the affine linear - * parametrization's Jacobian - */ + /// Returns the absolute value of the determinant of the affine linear + /// parametrization's Jacobian virtual double calcDet() const; - /** \brief - * Used by non static method \ref calcDet(). Calculates the determinant - * for a given vector of vertex coordinates. - */ + /// Used by non static method \ref calcDet(). Calculates the determinant + /// for a given vector of vertex coordinates. double calcDet(const FixVec, VERTEX> &coords) const; /// from CFE_Integration double calcSurfaceDet(VectorOfFixVecs > &surfVert) const; - /// Checks whether flag is set in ElInfo's \ref fillFlag. If not, the program exits. + /// Checks whether flag is set in ElInfo's \ref fillFlag. If not, the + /// program exits. void testFlag(const Flag& flag) const; - /** \brief - * Transforms local barycentric coordinates of a point defined on this element - * to global world coordinates. - */ + /// Transforms local barycentric coordinates of a point defined on this + /// element to global world coordinates. void coordToWorld(const DimVec& lambda, WorldVector& world) const; /// Fills ElInfo's \ref det and \ref grdLambda entries. virtual void fillDetGrdLambda(); - /** \brief - * Returns a pointer to a vector, which contains the barycentric coordinates - * with respect to \ref element of a point with world coordinates world. - * The barycentric coordinates are stored in lambda. - * pure virtual => must be overriden in sub-class. - */ + /// Returns a pointer to a vector, which contains the barycentric coordinates + /// with respect to \ref element of a point with world coordinates world. + /// The barycentric coordinates are stored in lambda. + /// pure virtual => must be overriden in sub-class. virtual const int worldToCoord(const WorldVector& world, DimVec* lambda) const = 0; - /** \brief - * Fills this ElInfo with macro element information of mel. - * pure virtual => must be overriden in sub-class. - */ + /// Fills this ElInfo with macro element information of mel. + /// pure virtual => must be overriden in sub-class. virtual void fillMacroInfo(const MacroElement *mel) = 0; - /** \brief - * Fills this ElInfo for the child ichild using hierarchy information and - * parent data parentInfo. - * pure virtual => must be overriden in sub-class. - */ + /// Fills this ElInfo for the child ichild using hierarchy information and + /// parent data parentInfo. + /// pure virtual => must be overriden in sub-class. virtual void fillElInfo(int ichild, const ElInfo *parentInfo) = 0; - void fillElInfo(const MacroElement *mel, int refinementPathLength, unsigned long refinementPath); + void fillElInfo(const MacroElement *mel, + int refinementPathLength, + unsigned long refinementPath); - /** \brief - * Calculates the Jacobian of the barycentric coordinates on \element and - * stores the matrix in grd_lam. The return value of the function is the - * absolute value of the determinant of the affine linear paraetrization's - * Jacobian. - * pure virtual => must be overriden in sub-class. - */ + /// Calculates the Jacobian of the barycentric coordinates on \element and + /// stores the matrix in grd_lam. The return value of the function is the + /// absolute value of the determinant of the affine linear paraetrization's + /// Jacobian. + /// pure virtual => must be overriden in sub-class. virtual double calcGrdLambda(DimVec >& grd_lam) = 0; - /** \brief - * calculates a normal of the given side (1d, 2d: edge, 3d: face) of \ref element. - * Returns the absolute value of the determinant of the - * transformation to the reference element. - * pure virtual => must be overriden in sub-class. - */ + /// calculates a normal of the given side (1d, 2d: edge, 3d: face) of \ref element. + /// Returns the absolute value of the determinant of the + /// transformation to the reference element. + /// pure virtual => must be overriden in sub-class. virtual double getNormal(int side, WorldVector &normal) const = 0; - - /** \brief - * calculates a normal of the element in dim of world = dim + 1. - * Returns the absolute value of the determinant of the - * transformation to the reference element. - * pure virtual => must be overriden in sub-class. - */ + /// calculates a normal of the element in dim of world = dim + 1. + /// Returns the absolute value of the determinant of the + /// transformation to the reference element. + /// pure virtual => must be overriden in sub-class. virtual double getElementNormal(WorldVector &elementNormal) const { FUNCNAME("ElInfo::getElementNormal()"); diff --git a/AMDiS/src/FiniteElemSpace.cc b/AMDiS/src/FiniteElemSpace.cc index 5964a88e5dae7fbed0822d38ffe0b99980cd96f5..62c9c057ca19e59fff304b3dd308a1e337beb56a 100644 --- a/AMDiS/src/FiniteElemSpace.cc +++ b/AMDiS/src/FiniteElemSpace.cc @@ -63,7 +63,6 @@ namespace AMDiS { } - FiniteElemSpace::FiniteElemSpace() {} @@ -104,6 +103,14 @@ namespace AMDiS { } + void FiniteElemSpace::destroyFeSpaces() + { + for (unsigned int i = 0; i < feSpaces.size(); i++) + delete feSpaces[i]; + + feSpaces.resize(0); + } + #if DEBUG FiniteElemSpace* FiniteElemSpace::provideFeSpace(Mesh *mesh) { diff --git a/AMDiS/src/FiniteElemSpace.h b/AMDiS/src/FiniteElemSpace.h index 177a9fd9cf5a9f20e31dc08249add5826ce59067..ddcb9c9dd2030a3585131bc1973b4002979d7898 100644 --- a/AMDiS/src/FiniteElemSpace.h +++ b/AMDiS/src/FiniteElemSpace.h @@ -51,6 +51,8 @@ namespace AMDiS { Mesh *mesh, string name = ""); + static void destroyFeSpaces(); + #if DEBUG /// For debugging it may be useful to get some FE space for a given mesh at a /// position in code where it is not possible to access the FE space directly. The @@ -104,10 +106,8 @@ namespace AMDiS { getHighest(vector& feSpaces); protected: - /** \brief - * Constructs a FiniteElemSpace with name name_ and the given DOFAdmin, - * BasisFunction and Mesh. - */ + /// Constructs a FiniteElemSpace with name name_ and the given DOFAdmin, + /// BasisFunction and Mesh. FiniteElemSpace(DOFAdmin* admin, const BasisFunction* basisFcts, Mesh* mesh, diff --git a/AMDiS/src/FirstOrderTerm.cc b/AMDiS/src/FirstOrderTerm.cc index cacc0986c61978c24bf94b076f8f5c0d48f7fc43..ce0c3944ace2bc030582e20e14ba80dc8d581b06 100644 --- a/AMDiS/src/FirstOrderTerm.cc +++ b/AMDiS/src/FirstOrderTerm.cc @@ -99,7 +99,7 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, @@ -140,7 +140,7 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo, @@ -281,7 +281,7 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } @@ -477,7 +477,7 @@ namespace AMDiS { Quadrature *quad) { getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs); - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void FctVecAtQP_FOT::getLb(const ElInfo *elInfo, @@ -686,7 +686,7 @@ namespace AMDiS { int nVecs = static_cast(vecs_.size()); int nGrads = static_cast(grads_.size()); - coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); for (int i = 0; i < nVecs; i++) getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]); @@ -713,7 +713,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - lb(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), Lb[iq], 1.0); + lb(grdLambda, (*f_)(coordsAtQPs[iq], vecsArg, gradsArg), Lb[iq], 1.0); } } @@ -741,7 +741,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - const WorldVector &b = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg); + const WorldVector &b = (*f_)(coordsAtQPs[iq], vecsArg, gradsArg); for (int i = 0; i < dow; i++) resultQP += grdUhAtQP[iq][i] * b[i]; @@ -788,7 +788,7 @@ namespace AMDiS { int nGrads = static_cast(grads_.size()); elInfo->getElementNormal(elementNormal_); - coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); for (int i = 0; i < nVecs; i++) getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]); @@ -815,18 +815,18 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - lb(grdLambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg), + lb(grdLambda, (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg), Lb[iq], 1.0); } } - void GeneralParametric_FOT::eval( int nPoints, - const mtl::dense_vector& uhAtQP, - const mtl::dense_vector >& grdUhAtQP, - const mtl::dense_vector >& D2UhAtQP, - mtl::dense_vector& result, - double factor) + void GeneralParametric_FOT::eval(int nPoints, + const mtl::dense_vector& uhAtQP, + const mtl::dense_vector >& grdUhAtQP, + const mtl::dense_vector >& D2UhAtQP, + mtl::dense_vector& result, + double factor) { int dow = Global::getGeo(WORLD); int nVecs = static_cast(vecs_.size()); @@ -845,7 +845,7 @@ namespace AMDiS { gradsArg[i] = gradsAtQPs_[i][iq]; const WorldVector &b = - (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg); + (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg); for (int i = 0; i < dow; i++) resultQP += grdUhAtQP[iq][i] * b[i]; diff --git a/AMDiS/src/FirstOrderTerm.h b/AMDiS/src/FirstOrderTerm.h index 01bafb942671515a3ece3ec1fdd6eba1b5e98f2e..773bf488dcb32f2991ce5c4441037e53dada7323 100644 --- a/AMDiS/src/FirstOrderTerm.h +++ b/AMDiS/src/FirstOrderTerm.h @@ -336,7 +336,7 @@ namespace AMDiS { protected: /// Stores coordinates at quadrature points. Set in \ref initElement(). - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function avaluated at world coordinates. AbstractFunction > *g; @@ -375,7 +375,7 @@ namespace AMDiS { protected: /// Stores coordinates at quadrature points. Set in \ref initElement(). - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function evaluated at world coordinates. AbstractFunction > *g; @@ -495,7 +495,7 @@ namespace AMDiS { protected: /// Stores coordinates at quadrature points. Set in \ref initElement(). - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function avaluated at world coordinates. AbstractFunction, WorldVector > *g; @@ -657,9 +657,13 @@ namespace AMDiS { protected: DOFVectorBase* vec; + mtl::dense_vector vecAtQPs; - WorldVector *coordsAtQPs; + + mtl::dense_vector > coordsAtQPs; + BinaryAbstractFunction, double> *f; + WorldVector *b; }; @@ -772,7 +776,7 @@ namespace AMDiS { vector, vector > > *f_; - WorldVector *coordsAtQPs_; + mtl::dense_vector > coordsAtQPs; vector > vecsAtQPs; @@ -818,7 +822,8 @@ namespace AMDiS { vector, vector > > *f_; - WorldVector *coordsAtQPs_; + mtl::dense_vector > coordsAtQPs; + WorldVector elementNormal_; vector > vecsAtQPs; diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h index dfa76cae772575e8dd726670a469b74a5e1c9cf0..c89c563dafa06b719d0840a4879be2a4bf54da35 100644 --- a/AMDiS/src/FixVec.h +++ b/AMDiS/src/FixVec.h @@ -53,20 +53,16 @@ namespace AMDiS { { public: - /** \brief - * constructor without initialisation. initType must be NO_INIT. If dim is - * not spezified, a FixVec for DIM_OF_WORLD is created. - */ + /// Constructor without initialisation. initType must be NO_INIT. If dim is + /// not spezified, a FixVec for DIM_OF_WORLD is created. FixVec(int dim = -1, InitType initType = NO_INIT) : Vector(calcSize(dim)) { TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or missing initializer\n"); } - /** \brief - * constructor with value list initialisation. initType must be VALUE_LIST. - * ini is an array which contains the initialisation values. - */ + /// constructor with value list initialisation. initType must be VALUE_LIST. + /// ini is an array which contains the initialisation values. FixVec(int dim, InitType initType, const T* ini) : Vector(calcSize(dim)) { @@ -74,10 +70,8 @@ namespace AMDiS { setValues(ini); } - /** \brief - * constructor with default value initialisation. initType must be - * DEFAULT_VALUE. All vector entries are set to ini. - */ + /// constructor with default value initialisation. initType must be + /// DEFAULT_VALUE. All vector entries are set to ini. FixVec(int dim, InitType initType, const T& ini) : Vector(calcSize(dim)) { @@ -131,11 +125,9 @@ namespace AMDiS { class VectorOfFixVecs { public: - /** \brief - * constructs a VectorOfFixVecs without initialisation. dim is passed to - * FixVec's constructors. size_ is the number of contained FixVecs. initType - * must be NO_INIT. - */ + /// constructs a VectorOfFixVecs without initialisation. dim is passed to + /// FixVec's constructors. size_ is the number of contained FixVecs. initType + /// must be NO_INIT. VectorOfFixVecs(int d, int s, InitType initType) : size(s), dim(d) @@ -147,11 +139,9 @@ namespace AMDiS { vec[i] = new FixVecType(dim, NO_INIT); } - /** \brief - * constructs a VectorOfFixVecs via an value list. dim is passed to - * FixVec's constructors. size_ is the number of contained FixVecs. initType - * must be VALUE_LIST. ini contains the initialisation values. - */ + /// constructs a VectorOfFixVecs via an value list. dim is passed to + /// FixVec's constructors. size_ is the number of contained FixVecs. initType + /// must be VALUE_LIST. ini contains the initialisation values. VectorOfFixVecs(int d, int s, InitType initType, const FixVecType* ini) : size(s), dim(d) @@ -163,11 +153,9 @@ namespace AMDiS { vec[i] = new FixVecType(ini[i]); } - /** \brief - * constructs a VectorOfFixVecs with an default value. dim is passed to - * FixVec's constructors. size_ is the number of contained FixVecs. initType - * must be DEFAULT_VALUE. All entries are set to ini. - */ + /// constructs a VectorOfFixVecs with an default value. dim is passed to + /// FixVec's constructors. size_ is the number of contained FixVecs. initType + /// must be DEFAULT_VALUE. All entries are set to ini. VectorOfFixVecs(int d, int s, InitType initType, const FixVecType& ini) : size(s), dim(d) @@ -273,11 +261,9 @@ namespace AMDiS { class MatrixOfFixVecs { public: - /** \brief - * Constructs the matrix without initialisation. r is the number of rows, - * c is the number of columns. The other parameters are analog to the - * VectorOfFixVecs constructors. - */ + /// Constructs the matrix without initialisation. r is the number of rows, + /// c is the number of columns. The other parameters are analog to the + /// VectorOfFixVecs constructors. MatrixOfFixVecs(int dim, int r, int c, InitType initType) : rows(r), columns(c) diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc index 49d309125c98ad7098eb3cc2bc21e78d8e53e655..21d597780e3f4f1f6049849ae69e1f67006beaf5 100644 --- a/AMDiS/src/Lagrange.cc +++ b/AMDiS/src/Lagrange.cc @@ -55,6 +55,7 @@ namespace AMDiS { setFunctionPointer(); } + Lagrange::~Lagrange() { for (int i = 0; i < static_cast(bary->size()); i++) @@ -64,6 +65,7 @@ namespace AMDiS { } } + Lagrange* Lagrange::getLagrange(int dim, int degree) { std::list::iterator it; @@ -76,6 +78,7 @@ namespace AMDiS { return newLagrange; } + void Lagrange::clear() { for (std::list::iterator it = allBasFcts.begin(); @@ -86,6 +89,7 @@ namespace AMDiS { } } + void Lagrange::setFunctionPointer() { if (static_cast(phiDimDegree[dim][degree].size()) == 0) { @@ -193,6 +197,7 @@ namespace AMDiS { } } + void Lagrange::setNDOF() { if (static_cast(baryDimDegree[dim][degree].size()) == 0) { @@ -781,9 +786,9 @@ namespace AMDiS { return NULL; } - void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const + void Lagrange::getBound(const ElInfo* elInfo, BoundaryType* bound) const { - el_info->testFlag(Mesh::FILL_BOUND); + elInfo->testFlag(Mesh::FILL_BOUND); // boundaries int index = 0; @@ -795,7 +800,7 @@ namespace AMDiS { for (int i = 0; i < dim; i++) { int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim); for (int j = offset; j < jto; j++) { - boundaryType = el_info->getBoundary(j); + boundaryType = elInfo->getBoundary(j); int kto = (*nDOF)[INDEX_OF_DIM(i, dim)]; for (int k = 0; k < kto; k++) bound[index++] = boundaryType; @@ -811,7 +816,7 @@ namespace AMDiS { } - void Lagrange::interpol(const ElInfo *el_info, + void Lagrange::interpol(const ElInfo *elInfo, int no, const int *b_no, AbstractFunction > *f, @@ -821,32 +826,32 @@ namespace AMDiS { WorldVector x; - el_info->testFlag(Mesh::FILL_COORDS); + elInfo->testFlag(Mesh::FILL_COORDS); if (b_no) { TEST_EXIT_DBG(no >= 0 && no < getNumber())("Something is wrong!\n"); for (int i = 0; i < no; i++) { if (b_no[i] < Global::getGeo(VERTEX, dim)) { - rvec[i] = (*f)(el_info->getCoord(b_no[i])); + rvec[i] = (*f)(elInfo->getCoord(b_no[i])); } else { - el_info->coordToWorld(*(*bary)[b_no[i]], x); + elInfo->coordToWorld(*(*bary)[b_no[i]], x); rvec[i] = (*f)(x); } } } else { int vertices = Global::getGeo(VERTEX, dim); for (int i = 0; i < vertices; i++) - rvec[i] = (*f)(el_info->getCoord(i)); + rvec[i] = (*f)(elInfo->getCoord(i)); for (int i = vertices; i < nBasFcts; i++) { - el_info->coordToWorld(*(*bary)[i], x); + elInfo->coordToWorld(*(*bary)[i], x); rvec[i] = (*f)(x); } } } - void Lagrange::interpol(const ElInfo *el_info, + void Lagrange::interpol(const ElInfo *elInfo, int no, const int *b_no, AbstractFunction, WorldVector > *f, @@ -856,7 +861,7 @@ namespace AMDiS { WorldVector x; - el_info->testFlag(Mesh::FILL_COORDS); + elInfo->testFlag(Mesh::FILL_COORDS); int vertices = Global::getGeo(VERTEX, dim); @@ -865,17 +870,17 @@ namespace AMDiS { for (int i = 0; i < no; i++) { if (b_no[i] < Global::getGeo(VERTEX, dim)) { - rvec[i] = (*f)(el_info->getCoord(b_no[i])); + rvec[i] = (*f)(elInfo->getCoord(b_no[i])); } else { - el_info->coordToWorld(*(*bary)[b_no[i]], x); + elInfo->coordToWorld(*(*bary)[b_no[i]], x); rvec[i] = (*f)(x); } } } else { for (int i = 0; i < vertices; i++) - rvec[i] = (*f)(el_info->getCoord(i)); + rvec[i] = (*f)(elInfo->getCoord(i)); for (int i = vertices; i < nBasFcts; i++) { - el_info->coordToWorld(*(*bary)[i], x); + elInfo->coordToWorld(*(*bary)[i], x); rvec[i] = (*f)(x); } } @@ -992,14 +997,14 @@ namespace AMDiS { WorldVector x; TraverseStack stack; - ElInfo *el_info = stack.traverseFirst(fh->getFeSpace()->getMesh(), -1, + ElInfo *elInfo = stack.traverseFirst(fh->getFeSpace()->getMesh(), -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); - while (el_info) { - const DegreeOfFreedom *dof = getLocalIndices(el_info->getElement(), admin, NULL); - double det = el_info->getDet(); + while (elInfo) { + const DegreeOfFreedom *dof = getLocalIndices(elInfo->getElement(), admin, NULL); + double det = elInfo->getDet(); for (int iq = 0; iq < nPoints; iq++) { - el_info->coordToWorld(quad->getLambda(iq), x); + elInfo->coordToWorld(quad->getLambda(iq), x); wdetf_qp[iq] = quad->getWeight(iq) * det * ((*f)(x)); } @@ -1011,7 +1016,7 @@ namespace AMDiS { (*fh)[dof[j]] += val; } - el_info = stack.traverseNext(el_info); + elInfo = stack.traverseNext(elInfo); } delete [] wdetf_qp; @@ -1020,7 +1025,10 @@ namespace AMDiS { // ===== refineInter functions ===== - void Lagrange::refineInter0(DOFIndexed *drv, RCNeighbourList* list, int n, BasisFunction* basFct) + void Lagrange::refineInter0(DOFIndexed *drv, + RCNeighbourList* list, + int n, + BasisFunction* basFct) { FUNCNAME("Lagrange::refineInter0()"); @@ -1041,7 +1049,9 @@ namespace AMDiS { (*drv)[dof_new] = (*drv)[dof0]; } - void Lagrange::refineInter1(DOFIndexed *drv, RCNeighbourList* list, int n, + void Lagrange::refineInter1(DOFIndexed *drv, + RCNeighbourList* list, + int n, BasisFunction* basFct) { FUNCNAME("Lagrange::refineInter1_1d()"); @@ -1061,7 +1071,9 @@ namespace AMDiS { (*drv)[dof_new] = 0.5 * ((*drv)[dof0] + (*drv)[dof1]); } - void Lagrange::refineInter2_1d(DOFIndexed *drv, RCNeighbourList* list, int n, + void Lagrange::refineInter2_1d(DOFIndexed *drv, + RCNeighbourList* list, + int n, BasisFunction* basFct) { FUNCNAME("Lagrange::refineInter2_1d()"); diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h index 2e0dbf79031aea2aefda2bd53c2725ae62b2cce7..52b56780cd5e778b082f47aa98d072a0f6f8ec13 100644 --- a/AMDiS/src/Lagrange.h +++ b/AMDiS/src/Lagrange.h @@ -43,11 +43,9 @@ namespace AMDiS { class Lagrange : public BasisFunction { protected: - /** \brief - * Constructs lagrange basis functions with the given dim and degree. - * Constructor is protected to avoid multiple instantiation of identical - * basis functions. Use \ref getLagrange instead. - */ + /// Constructs lagrange basis functions with the given dim and degree. + /// Constructor is protected to avoid multiple instantiation of identical + /// basis functions. Use \ref getLagrange instead. Lagrange(int dim_, int degree_); /** \brief @@ -56,11 +54,9 @@ namespace AMDiS { virtual ~Lagrange(); public: - /** \brief - * Returns a pointer to lagrange basis functions with the given dim and - * degree. Multiple instantiation of identical basis functions is avoided - * by rembering once created basis functions in \ref allBasFcts. - */ + /// Returns a pointer to lagrange basis functions with the given dim and + /// degree. Multiple instantiation of identical basis functions is avoided + /// by rembering once created basis functions in \ref allBasFcts. static Lagrange* getLagrange(int dim, int degree); /// Implements BasisFunction::interpol @@ -141,10 +137,8 @@ namespace AMDiS { static void clear(); protected: - /** \brief - * sets the barycentric coordinates (stored in \ref bary) of the local - * basis functions. - */ + /// sets the barycentric coordinates (stored in \ref bary) of the local + /// basis functions. void setBary(); /// Recursive calculation of coordinates. Used by \ref setBary @@ -183,10 +177,8 @@ namespace AMDiS { static std::vector D2PhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1]; /** \} */ - /** \brief - * List of all used BasisFunctions in the whole program. Avoids duplicate - * instantiation of identical BasisFunctions. - */ + /// List of all used BasisFunctions in the whole program. Avoids duplicate + /// instantiation of identical BasisFunctions. static std::list allBasFcts; @@ -286,11 +278,9 @@ namespace AMDiS { class Phi : public BasFctType { public: - /** \brief - * Constructs the local lagrange basis function for the given position, - * positionIndex and nodeIndex. owner_ is a pointer to the Lagrange object - * this basis function belongs to. - */ + /// Constructs the local lagrange basis function for the given position, + /// positionIndex and nodeIndex. owner_ is a pointer to the Lagrange + /// object this basis function belongs to. Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex); /// Destructor @@ -403,10 +393,8 @@ namespace AMDiS { - /** \brief - * AbstractFunction which implements gradients of lagrange basis functions. - * See \ref Phi - */ + /// AbstractFunction which implements gradients of lagrange basis functions. + /// See \ref Phi class GrdPhi : public GrdBasFctType { public: @@ -568,10 +556,8 @@ namespace AMDiS { - /** \brief - * AbstractFunction which implements second derivatives of Lagrange basis - * functions. See \ref Phi - */ + /// AbstractFunction which implements second derivatives of Lagrange basis + /// functions. See \ref Phi class D2Phi : public D2BasFctType { public: diff --git a/AMDiS/src/MatrixVector.h b/AMDiS/src/MatrixVector.h index b8d06ce809b692b8111835e99dd2c86367c2c740..3947983b6af6a80c10baa744bc409354b61cb8ec 100644 --- a/AMDiS/src/MatrixVector.h +++ b/AMDiS/src/MatrixVector.h @@ -67,7 +67,8 @@ namespace AMDiS { /// Copy constructor. Vector(const Vector& rhs) - : Serializable(),size(rhs.size) + : Serializable(), + size(rhs.size) { valArray = new T[rhs.size]; *this = rhs; // uses operator=() @@ -281,13 +282,13 @@ namespace AMDiS { return !(*this == rhs); } - /// Acces to i-th matrix row. + /// Access to i-th matrix row. inline T *operator[](int i) { return this->valArray + cols * i; } - /// Acces to i-th matrix row for constant matrices. + /// Access to i-th matrix row for constant matrices. inline const T *operator[](int i) const { return this->valArray + cols * i; diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index e29b0b0735effdd199a295b469b15392cfb04779..10f067efbfdf9f412995d5bb46729af5ba5fc273 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -64,7 +64,7 @@ namespace AMDiS { Mesh(string name, int dim); /// Destructor - virtual ~Mesh(); + ~Mesh(); /// Reads macro triangulation. void initialize(); @@ -76,10 +76,8 @@ namespace AMDiS { * \{ */ - /** \brief - * Returns geometric information about this mesh. With GeoIndex p it is - * specified which information is requested. - */ + /// Returns geometric information about this mesh. With GeoIndex p it is + /// specified which information is requested. inline int getGeo(GeoIndex p) const { return Global::getGeo(p, dim); @@ -388,7 +386,7 @@ namespace AMDiS { void dofCompress(); /// Adds a DOFAdmin to the mesh - virtual void addDOFAdmin(DOFAdmin *admin); + void addDOFAdmin(DOFAdmin *admin); /// Recalculates the number of leave elements. void updateNumberOfLeaves(); diff --git a/AMDiS/src/OperatorTerm.h b/AMDiS/src/OperatorTerm.h index 1bfa57f8de70d6a8f095501c13c9188fe069e9da..25b09bf7d81b613722b31d259855d6b06ee8efb0 100644 --- a/AMDiS/src/OperatorTerm.h +++ b/AMDiS/src/OperatorTerm.h @@ -41,11 +41,9 @@ namespace AMDiS { class OperatorTerm { public: - /** \brief - * Constructs an OperatorTerm with initially no properties. - * degree_ is used to determine the degree of the needed quadrature - * for the assemblage. - */ + /// Constructs an OperatorTerm with initially no properties. + /// degree_ is used to determine the degree of the needed quadrature + /// for the assemblage. OperatorTerm(int deg) : properties(0), degree(deg), @@ -56,20 +54,22 @@ namespace AMDiS { /// Destructor. virtual ~OperatorTerm() {} - /** \brief - * Virtual method. It's called by SubAssembler::initElement() for - * each OperatorTerm belonging to this SubAssembler. E.g., vectors - * and coordinates at quadrature points can be calculated here. - */ - virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL) + /// Virtual method. It's called by SubAssembler::initElement() for + /// each OperatorTerm belonging to this SubAssembler. E.g., vectors + /// and coordinates at quadrature points can be calculated here. + virtual void initElement(const ElInfo*, SubAssembler*, + Quadrature *quad = NULL) {} - virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo, - SubAssembler*, Quadrature *quad = NULL) + virtual void initElement(const ElInfo* largeElInfo, + const ElInfo* smallElInfo, + SubAssembler*, + Quadrature *quad = NULL) {} - /// Returs \auxFeSpaces, the list of all aux fe spaces the operator makes use off. + /// Returs \auxFeSpaces, the list of all aux fe spaces the operator makes + /// use off. inline std::set& getAuxFeSpaces() { return auxFeSpaces; @@ -78,7 +78,7 @@ namespace AMDiS { /// Specifies whether the matrix of the term is symmetric void setSymmetric(bool symm); - /// Returns true, if the term is piecewise constant, returns false otherwise + /// Returns true, if the term is piecewise constant, returns false otherwise. inline bool isPWConst() { return (degree == 0); diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc index c0e66f3e9fa9f363fc52cc0225aa5bb6cca384db..cb5833b441d465b370c3b6d6210047e1b2d79c04 100644 --- a/AMDiS/src/ProblemStat.cc +++ b/AMDiS/src/ProblemStat.cc @@ -99,6 +99,35 @@ namespace AMDiS { Parameters::get("debug->write asm info", writeAsmInfo); } + + ProblemStatSeq::~ProblemStatSeq() + { + delete rhs; + rhs = NULL; + delete solution; + solution = NULL; + + for (int i = 0; i < nComponents; i++) + for (int j = 0; j < nComponents; j++) + if ((*systemMatrix)[i][j]) { + delete (*systemMatrix)[i][j]; + (*systemMatrix)[i][j] = NULL; + } + + delete systemMatrix; + systemMatrix = NULL; + + for (unsigned int i = 0; i < meshes.size(); i++) + delete meshes[i]; + + for (unsigned int i = 0; i < estimator.size(); i++) + delete estimator[i]; + + for (unsigned int i = 0; i < marker.size(); i++) + delete marker[i]; + } + + void ProblemStatSeq::initialize(Flag initFlag, ProblemStatSeq *adoptProblem, Flag adoptFlag) @@ -395,7 +424,8 @@ namespace AMDiS { int degree = 1; Parameters::get(name + "->polynomial degree[" + boost::lexical_cast(i) + "]", degree); - + TEST_EXIT(degree > 0) + ("Poynomial degree in component %d must be larger than zero!\n", i); TEST_EXIT(componentSpaces[i] == NULL)("feSpace already created\n"); if (feSpaceMap[pair(componentMeshes[i], degree)] == NULL) { diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h index 7fd02cb10f270ce5514f7506a5645db4e6a4686c..4fd293e238cca5284d6d18df6392a54240f6614d 100644 --- a/AMDiS/src/ProblemStat.h +++ b/AMDiS/src/ProblemStat.h @@ -68,7 +68,7 @@ namespace AMDiS { ProblemIterationInterface *problemIteration = NULL); /// Destructor - virtual ~ProblemStatSeq() {} + virtual ~ProblemStatSeq(); /// Initialisation of the problem. virtual void initialize(Flag initFlag, diff --git a/AMDiS/src/SecondOrderTerm.cc b/AMDiS/src/SecondOrderTerm.cc index 2c3bd0c42a1056b1ecbe32c25d5cee76b44cf169..cf5ca67e64efc7cdf70990b7626a4c61d160d420 100644 --- a/AMDiS/src/SecondOrderTerm.cc +++ b/AMDiS/src/SecondOrderTerm.cc @@ -346,7 +346,7 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, @@ -382,6 +382,11 @@ namespace AMDiS { void CoordsAtQP_SOT::weakEval(const std::vector > &grdUhAtQP, std::vector > &result) { + FUNCNAME("CoordsAtQP_SOT::weakEval()"); + TEST_EXIT_DBG(grdUhAtQP.size() == num_rows(coordsAtQPs)) + ("Wrong sizes! grdUhAtQP = %d coordsAtQPs = %d\n", + grdUhAtQP.size(), num_rows(coordsAtQPs)); + int nPoints = grdUhAtQP.size(); for (int iq = 0; iq < nPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); @@ -627,7 +632,7 @@ namespace AMDiS { { getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs); getGradientsAtQPs(vec, elInfo, subAssembler, quad, gradAtQPs); - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, @@ -689,7 +694,7 @@ namespace AMDiS { Quadrature *quad) { getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs); - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, @@ -755,7 +760,7 @@ namespace AMDiS { Quadrature *quad) { getGradientsAtQPs(vec, elInfo, subAssembler, quad, gradAtQPs); - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, @@ -919,7 +924,7 @@ namespace AMDiS { int nVecs = static_cast(vecs_.size()); int nGrads = static_cast(grads_.size()); - coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); for (int i = 0; i < nVecs; i++) getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]); @@ -945,7 +950,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - lalt(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), + lalt(grdLambda, (*f_)(coordsAtQPs[iq], vecsArg, gradsArg), LALt[iq], symmetric_, 1.0); } } @@ -972,7 +977,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - WorldMatrix A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg); + WorldMatrix A = (*f_)(coordsAtQPs[iq], vecsArg, gradsArg); if (num_rows(D2UhAtQP) > 0) { for (int i = 0; i < dow; i++) @@ -1004,7 +1009,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg); + A = (*f_)(coordsAtQPs[iq], vecsArg, gradsArg); result[iq] += A * grdUhAtQP[iq]; } } @@ -1053,7 +1058,7 @@ namespace AMDiS { int nGrads = static_cast(grads_.size()); elInfo->getElementNormal(elementNormal_); - coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); for (int i = 0; i < nVecs; i++) getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]); @@ -1078,7 +1083,7 @@ namespace AMDiS { vecsArg[i] = vecsAtQPs[i][iq]; for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - lalt(grdLambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg), + lalt(grdLambda, (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg), LALt[iq], symmetric_, 1.0); } } @@ -1106,7 +1111,7 @@ namespace AMDiS { gradsArg[i] = gradsAtQPs_[i][iq]; WorldMatrix A = - (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg); + (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg); if (num_rows(D2UhAtQP) > 0) { for (int i = 0; i < dow; i++) @@ -1138,7 +1143,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - A = (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg); + A = (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg); result[iq] += A * grdUhAtQP[iq]; } } @@ -1262,7 +1267,7 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void CoordsAtQP_IJ_SOT::getLALt(const ElInfo *elInfo, diff --git a/AMDiS/src/SecondOrderTerm.h b/AMDiS/src/SecondOrderTerm.h index b2f67ceaf82f887b17d03fd7afb1cd171d094f39..879eb8d797df9185a23cefa01c9647f6c81852b8 100644 --- a/AMDiS/src/SecondOrderTerm.h +++ b/AMDiS/src/SecondOrderTerm.h @@ -485,7 +485,7 @@ namespace AMDiS { protected: /// Stores coordinates at quadrature points. Set in \ref initElement(). - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function evaluated at quadrature points. AbstractFunction > *g; @@ -704,15 +704,13 @@ namespace AMDiS { TertiaryAbstractFunction, WorldVector > *f; - /** \brief - * Pointer to a WorldVector array containing the gradients of the DOFVector - * at each quadrature point. - */ + /// Pointer to a WorldVector array containing the gradients of the + /// DOFVector at each quadrature point. mtl::dense_vector vecAtQPs; mtl::dense_vector > gradAtQPs; - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; }; @@ -754,7 +752,7 @@ namespace AMDiS { /// Pointer to an array containing the DOFVector evaluated at quadrature points. mtl::dense_vector vecAtQPs; - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function evaluated at \ref vecAtQPs. BinaryAbstractFunction > *f; @@ -809,13 +807,11 @@ namespace AMDiS { AbstractFunction, WorldMatrix > *divFct; - /** \brief - * Pointer to a WorldVector array containing the gradients of the DOFVector - * at each quadrature point. - */ + /// Pointer to a WorldVector array containing the gradients of the + /// DOFVector at each quadrature point. mtl::dense_vector > gradAtQPs; - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// True, if \ref f provides always a symmetric WorldMatrix. bool symmetric; @@ -913,7 +909,7 @@ namespace AMDiS { AbstractFunction, WorldMatrix > *divFct_; - WorldVector *coordsAtQPs_; + mtl::dense_vector > coordsAtQPs; std::vector > vecsAtQPs; @@ -972,9 +968,10 @@ namespace AMDiS { AbstractFunction, WorldMatrix > *divFct_; - WorldVector *coordsAtQPs_; WorldVector elementNormal_; + mtl::dense_vector > coordsAtQPs; + std::vector > vecsAtQPs; std::vector > > gradsAtQPs_; @@ -1124,7 +1121,7 @@ namespace AMDiS { private: /// Stores coordinates at quadrature points. Set in \ref initElement(). - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function evaluated at quadrature points. AbstractFunction > *g; diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc index 3c3abf252c0f5bfaaba33bf36629329f1311a003..92ab6acfb7e7a95e5ebb309a503d3b6dc2f1ef8e 100644 --- a/AMDiS/src/SubAssembler.cc +++ b/AMDiS/src/SubAssembler.cc @@ -31,7 +31,6 @@ namespace AMDiS { colFeSpace(assembler->colFeSpace), nRow(0), nCol(0), - coordsAtQPs(NULL), coordsNumAllocated(0), quadrature(quadrat), psiFast(NULL), @@ -112,17 +111,17 @@ namespace AMDiS { coordsValid = false; // set values at QPs invalid - std::map::iterator itAny; + map::iterator itAny; for (itAny = cachedValuesAtQPs.begin(); itAny != cachedValuesAtQPs.end(); ++itAny) (*itAny).second->valid = false; // set gradients at QPs invalid - std::map::iterator itAnyGrd; + map::iterator itAnyGrd; for (itAnyGrd = cachedGradientsAtQPs.begin(); itAnyGrd != cachedGradientsAtQPs.end(); ++itAnyGrd) (*itAnyGrd).second->valid = false; // calls initElement of each term - for (std::vector::iterator it = terms.begin(); + for (vector::iterator it = terms.begin(); it != terms.end(); ++it) { if (largeElInfo == NULL) (*it)->initElement(smallElInfo, this, quad); @@ -132,39 +131,27 @@ namespace AMDiS { } - WorldVector* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo, - Quadrature *quad) + void SubAssembler::getCoordsAtQPs(const ElInfo* elInfo, + Quadrature *quad, + mtl::dense_vector > &coordsAtQPs) { - Quadrature *localQuad = quad ? quad : quadrature; - - const int nPoints = localQuad->getNumPoints(); - // already calculated for this element ? - if (coordsValid) - return coordsAtQPs; - - if (coordsAtQPs) { - if (coordsNumAllocated != nPoints) { - delete [] coordsAtQPs; - coordsAtQPs = new WorldVector[nPoints]; - coordsNumAllocated = nPoints; - } - } else { - coordsAtQPs = new WorldVector[nPoints]; - coordsNumAllocated = nPoints; + if (coordsValid) { + coordsAtQPs = cacheCoordsAtQPs; + return; } - // set new values - WorldVector* k = &(coordsAtQPs[0]); - for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l) - elInfo->coordToWorld(localQuad->getLambda(l), *k); + Quadrature *localQuad = quad ? quad : quadrature; + const int nPoints = localQuad->getNumPoints(); + + coordsAtQPs.change_dim(nPoints); + for (int i = 0; i < nPoints; i++) + elInfo->coordToWorld(localQuad->getLambda(i), coordsAtQPs[i]); // mark values as valid coordsValid = true; - - return coordsAtQPs; + cacheCoordsAtQPs.change_dim(num_rows(coordsAtQPs)); + cacheCoordsAtQPs = coordsAtQPs; } - - } diff --git a/AMDiS/src/SubAssembler.h b/AMDiS/src/SubAssembler.h index 9903841c718de42ad79a21cd7c50629ca396a06b..f59bb44b0266d0b580cd54db7c77bfa34a516b30 100644 --- a/AMDiS/src/SubAssembler.h +++ b/AMDiS/src/SubAssembler.h @@ -31,6 +31,8 @@ #include namespace AMDiS { + + using namespace std; /** * \ingroup Assembler @@ -48,12 +50,10 @@ namespace AMDiS { class SubAssembler { public: - /** \brief - * Creates a SubAssembler belonging to assembler for the terms of given - * order of Operator op. If order is equal to one, type spezifies what kind - * of FirstOrderType are to assemble. During construction of a SubAssembler - * the needs and properties of the terms are considered. - */ + /// Creates a SubAssembler belonging to assembler for the terms of given + /// order of Operator op. If order is equal to one, type spezifies what kind + /// of FirstOrderType are to assemble. During construction of a SubAssembler + /// the needs and properties of the terms are considered. SubAssembler(Operator *op, Assembler *assembler, Quadrature *quadrat, @@ -75,7 +75,7 @@ namespace AMDiS { ElementVector& vec) = 0; /// Returns \ref terms - inline std::vector *getTerms() + inline vector *getTerms() { return &terms; } @@ -92,13 +92,12 @@ namespace AMDiS { quadrature = q; } - /** \brief - * Returns a vector with the world coordinates of the quadrature points - * of \ref quadrature on the element of elInfo. - * Used by \ref OperatorTerm::initElement(). - */ - WorldVector* getCoordsAtQPs(const ElInfo* elInfo, - Quadrature *quad = NULL); + /// Creates a vector with the world coordinates of the quadrature points + /// of \ref quadrature on the element of elInfo. + /// Used by \ref OperatorTerm::initElement(). + void getCoordsAtQPs(const ElInfo* elInfo, + Quadrature *quad, + mtl::dense_vector >& coordsAtQPs); /// DOFVector dv evaluated at quadrature points. /// Used by \ref OperatorTerm::initElement(). @@ -119,17 +118,17 @@ namespace AMDiS { /// Gradients of DOFVector dv evaluated at quadrature points. /// Used by \ref OperatorTerm::initElement(). template - void getGradientsAtQPs( DOFVectorBase* dv, - const ElInfo* elInfo, - Quadrature *quad, - mtl::dense_vector::type>& grdAtQPs); + void getGradientsAtQPs(DOFVectorBase* dv, + const ElInfo* elInfo, + Quadrature *quad, + mtl::dense_vector::type>& grdAtQPs); template - void getGradientsAtQPs( DOFVectorBase* dv, - const ElInfo* smallElInfo, - const ElInfo* largeElInfo, - Quadrature *quad, - mtl::dense_vector::type>& grdAtQPs); + void getGradientsAtQPs(DOFVectorBase* dv, + const ElInfo* smallElInfo, + const ElInfo* largeElInfo, + Quadrature *quad, + mtl::dense_vector::type>& grdAtQPs); /// The comp'th component of the derivative of DOFVector dv evaluated at /// quadrature points. Used by \ref OperatorTerm::initElement(). @@ -170,7 +169,7 @@ namespace AMDiS { } /// Returns \ref name. - std::string getName() const + string getName() const { return name; } @@ -213,12 +212,12 @@ namespace AMDiS { Quadrature* quad; }; - std::map cachedValuesAtQPs; - std::map cachedGradientsAtQPs; + map cachedValuesAtQPs; + map cachedGradientsAtQPs; /// Set and updated by \ref initElement() for each ElInfo. /// coordsAtQPs[i] points to the coordinates of the i-th quadrature point. - WorldVector *coordsAtQPs; + mtl::dense_vector > cacheCoordsAtQPs; /// Used for \ref getCoordsAtQPs(). bool coordsValid; @@ -240,7 +239,7 @@ namespace AMDiS { bool symmetric; /// List of all terms with a contribution to this SubAssembler - std::vector terms; + vector terms; /// bool opt; @@ -248,8 +247,9 @@ namespace AMDiS { /// bool firstCall; - /// Name of the assembler. Is used to print information about used assembler. - std::string name; + /// Name of the assembler. Is used to print information about + /// used assembler. + string name; friend class Assembler; }; diff --git a/AMDiS/src/SubAssembler.hh b/AMDiS/src/SubAssembler.hh index 7dd43945270948ea2231634d5640da5dcfc82094..0768aa7f6b5226f4ada8dc0f9c9f5927819db4b4 100644 --- a/AMDiS/src/SubAssembler.hh +++ b/AMDiS/src/SubAssembler.hh @@ -111,11 +111,10 @@ namespace AMDiS { template - void SubAssembler::getGradientsAtQPs( DOFVectorBase* vec, - const ElInfo* elInfo, - Quadrature *quad, - mtl::dense_vector::type>& grdAtQPs - ) + void SubAssembler::getGradientsAtQPs(DOFVectorBase* vec, + const ElInfo* elInfo, + Quadrature *quad, + mtl::dense_vector::type>& grdAtQPs) { FUNCNAME("SubAssembler::getGradientsAtQPs()"); @@ -124,7 +123,8 @@ namespace AMDiS { grdAtQPs.change_dim(localQuad->getNumPoints()); if (cachedGradientsAtQPs[vec] && cachedGradientsAtQPs[vec]->valid) { - grdAtQPs = boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); + grdAtQPs = + boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); return; } @@ -136,7 +136,8 @@ namespace AMDiS { boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values).change_dim(localQuad->getNumPoints()); - mtl::dense_vector::type>& values = boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); + mtl::dense_vector::type>& values = + boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); const BasisFunction *psi = rowFeSpace->getBasisFcts(); const BasisFunction *phi = colFeSpace->getBasisFcts(); @@ -171,11 +172,11 @@ namespace AMDiS { template - void SubAssembler::getGradientsAtQPs( DOFVectorBase* vec, - const ElInfo* smallElInfo, - const ElInfo* largeElInfo, - Quadrature *quad, - mtl::dense_vector::type>& grdAtQPs) + void SubAssembler::getGradientsAtQPs(DOFVectorBase* vec, + const ElInfo* smallElInfo, + const ElInfo* largeElInfo, + Quadrature *quad, + mtl::dense_vector::type>& grdAtQPs) { FUNCNAME("SubAssembler::getGradientsAtQPs()"); @@ -185,7 +186,8 @@ namespace AMDiS { grdAtQPs.change_dim(localQuad->getNumPoints()); if (cachedGradientsAtQPs[vec] && cachedGradientsAtQPs[vec]->valid) { - grdAtQPs = boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); + grdAtQPs = + boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); return; } @@ -197,13 +199,15 @@ namespace AMDiS { boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values).change_dim(localQuad->getNumPoints()); - mtl::dense_vector::type>& values = boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); + mtl::dense_vector::type>& values = + boost::any_cast::type>& >(cachedGradientsAtQPs[vec]->values); cachedGradientsAtQPs[vec]->valid = true; vec->getGrdAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values); grdAtQPs = values; } + template void SubAssembler::getDerivativeAtQPs(DOFVectorBase* vec, const ElInfo* elInfo, diff --git a/AMDiS/src/ZeroOrderTerm.cc b/AMDiS/src/ZeroOrderTerm.cc index 57d0276d7bed59417fec05a39e2a3cda8d534136..85e28dde8b71f8cb5a8d879b534b7838f646c3e1 100644 --- a/AMDiS/src/ZeroOrderTerm.cc +++ b/AMDiS/src/ZeroOrderTerm.cc @@ -264,7 +264,7 @@ namespace AMDiS { Quadrature *quad) { getGradientsAtQPs(vec, elInfo, subAssembler, quad, gradAtQPs); - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) @@ -303,7 +303,7 @@ namespace AMDiS { { getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs); getGradientsAtQPs(vec, elInfo, subAssembler, quad, gradAtQPs); - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) @@ -343,7 +343,7 @@ namespace AMDiS { Quadrature *quad) { getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs); - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) @@ -707,10 +707,9 @@ namespace AMDiS { { int size = static_cast(vecs.size()); for (int i = 0; i < size; i++) - getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad, gradsAtQPs[i]); + getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad, gradsAtQPs[i]); } - void VecDivergence_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) { int size = static_cast(vecs.size()); @@ -720,12 +719,12 @@ namespace AMDiS { C[iq] += gradsAtQPs[i][iq][i]; } - void VecDivergence_ZOT::eval( int nPoints, - const mtl::dense_vector& uhAtQP, - const mtl::dense_vector >& grdUhAtQP, - const mtl::dense_vector >& D2UhAtQP, - mtl::dense_vector& result, - double fac) + void VecDivergence_ZOT::eval(int nPoints, + const mtl::dense_vector& uhAtQP, + const mtl::dense_vector >& grdUhAtQP, + const mtl::dense_vector >& D2UhAtQP, + mtl::dense_vector& result, + double fac) { int size = static_cast(vecs.size()); @@ -924,7 +923,7 @@ namespace AMDiS { int nVecs = static_cast(vecs_.size()); int nGrads = static_cast(grads_.size()); - coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); for (int i = 0; i < nVecs; i++) getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]); @@ -943,7 +942,7 @@ namespace AMDiS { for (unsigned int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - C[iq] += (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg); + C[iq] += (*f_)(coordsAtQPs[iq], vecsArg, gradsArg); } } @@ -963,7 +962,7 @@ namespace AMDiS { for (unsigned int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs_[i][iq]; - result[iq] += fac * (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg) * uhAtQP[iq]; + result[iq] += fac * (*f_)(coordsAtQPs[iq], vecsArg, gradsArg) * uhAtQP[iq]; } } @@ -1002,7 +1001,7 @@ namespace AMDiS { int nGrads = static_cast(grads_.size()); elInfo->getElementNormal(elementNormal_); - coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); for (int i = 0; i < nVecs; i++) getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]); @@ -1024,7 +1023,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs[i][iq]; - C[iq] += (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg); + C[iq] += (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg); } } @@ -1047,7 +1046,7 @@ namespace AMDiS { for (int i = 0; i < nGrads; i++) gradsArg[i] = gradsAtQPs[i][iq]; result[iq] += - fac * (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg) * uhAtQP[iq]; + fac * (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg) * uhAtQP[iq]; } } @@ -1058,7 +1057,7 @@ namespace AMDiS { SubAssembler* subAssembler, Quadrature *quad) { - coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad); + subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs); } void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector& C) diff --git a/AMDiS/src/ZeroOrderTerm.h b/AMDiS/src/ZeroOrderTerm.h index 986cfd535d156cc498548b2197d7e00041ec0310..7c6016a5bfbef0fb652bf90190ff89d794e75d03 100644 --- a/AMDiS/src/ZeroOrderTerm.h +++ b/AMDiS/src/ZeroOrderTerm.h @@ -305,14 +305,11 @@ namespace AMDiS { /// Function wich maps \ref gradAtQPs to a double. BinaryAbstractFunction, WorldVector > *f; - /** \brief - * Pointer to a WorldVector array containing the gradients of the DOFVector - * at each quadrature point. - */ + /// Pointer to a WorldVector array containing the gradients of the + /// DOFVector at each quadrature point. mtl::dense_vector > gradAtQPs; - WorldVector* coordsAtQPs; - + mtl::dense_vector > coordsAtQPs; }; /** @@ -352,7 +349,7 @@ namespace AMDiS { /// Gradient at quadrature points. mtl::dense_vector > gradAtQPs; - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function for c. TertiaryAbstractFunction, WorldVector > *f; @@ -392,8 +389,7 @@ namespace AMDiS { /// Vector v at quadrature points. mtl::dense_vector vecAtQPs; - /// Gradient at quadrature points. - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function for c. BinaryAbstractFunction > *f; @@ -854,7 +850,7 @@ namespace AMDiS { std::vector, std::vector > > *f_; - WorldVector *coordsAtQPs_; + mtl::dense_vector > coordsAtQPs; std::vector > vecsAtQPs; @@ -904,7 +900,7 @@ namespace AMDiS { std::vector, std::vector > > *f_; - WorldVector *coordsAtQPs_; + mtl::dense_vector > coordsAtQPs; WorldVector elementNormal_; @@ -946,7 +942,7 @@ namespace AMDiS { protected: /// Stores coordinates at quadrature points. Set in \ref initElement(). - WorldVector* coordsAtQPs; + mtl::dense_vector > coordsAtQPs; /// Function for c. AbstractFunction > *g; diff --git a/AMDiS/src/est/ResidualEstimator.cc b/AMDiS/src/est/ResidualEstimator.cc index 0c1921708be0326899fd7686d8febf059707cb05..5e6f4a8314d2b220e99251942054332eb9b2e99c 100644 --- a/AMDiS/src/est/ResidualEstimator.cc +++ b/AMDiS/src/est/ResidualEstimator.cc @@ -268,6 +268,8 @@ namespace AMDiS { elBoundGrdUhNeigh[subBound[i]] = stdMpiGrdUh.getRecvData(it->first)[i]; } } + + delete elInfo; } #endif diff --git a/AMDiS/src/io/MacroReader.cc b/AMDiS/src/io/MacroReader.cc index b87e3f8097f26798133b482e0e6772cbf06618a2..6e6e4b7a594deee7cfdfa16839aecfc3ea69d13d 100644 --- a/AMDiS/src/io/MacroReader.cc +++ b/AMDiS/src/io/MacroReader.cc @@ -325,14 +325,8 @@ namespace AMDiS { if (check) { checkMesh(mesh); - if (mesh->getDim() > 1) { - char filenew[128]; - strncpy(filenew, filename.c_str(), 128); - filenew[127] = 0; - strncat(filenew, ".new", 128); - filenew[127] = 0; - macroTest(mesh, filenew); - } + if (mesh->getDim() > 1) + macroTest(mesh); } return macroInfo; @@ -558,19 +552,16 @@ namespace AMDiS { (wenn nameneu=NULL wird keine MAcro-Datei erzeugt) */ - void MacroReader::macroTest(Mesh *mesh, const char *nameneu) + void MacroReader::macroTest(Mesh *mesh) { FUNCNAME("MacroReader::macroTest()"); int i = macrotest(mesh); if (i >= 0) { - ERROR("There is a cycle beginning in macro element %d\n", i); - ERROR("Entries in MacroElement structures get reordered\n"); + WARNING("There is a cycle beginning in macro element %d\n", i); + WARNING("Entries in MacroElement structures get reordered\n"); umb(NULL, mesh, umbVkantMacro); - - if (nameneu) - ERROR_EXIT("mesh->feSpace\n"); } } diff --git a/AMDiS/src/io/MacroReader.h b/AMDiS/src/io/MacroReader.h index c8bcc36087f3e64b040f1d7072708500172fada8..b407638b5615b12fa4845c614dccfa9adaa53405 100644 --- a/AMDiS/src/io/MacroReader.h +++ b/AMDiS/src/io/MacroReader.h @@ -105,7 +105,7 @@ namespace AMDiS { static int macrotest(Mesh *mesh); - static void macroTest(Mesh *mesh, const char *nameneu); + static void macroTest(Mesh *mesh); static bool newEdge(Mesh *mesh, MacroElement *mel, int mel_edge_no, int *n_neigh); diff --git a/AMDiS/src/io/MacroWriter.h b/AMDiS/src/io/MacroWriter.h index 26d4cf1c3f9a408b5431b0a60eeff4f20b093fe4..1fe2621936623e14781ac0fd3d13f59aaef312a2 100644 --- a/AMDiS/src/io/MacroWriter.h +++ b/AMDiS/src/io/MacroWriter.h @@ -94,10 +94,8 @@ namespace AMDiS { /// Maps internal element indices to global output indices. static std::map outputIndices; - /** \brief - * periodicConnections[i][j] stores whether the connection at side j of - * the element with output index i has already been written. - */ + /// periodicConnections[i][j] stores whether the connection at side j of + /// the element with output index i has already been written. static std::vector > periodicConnections; /// diff --git a/AMDiS/src/io/VtkWriter.cc b/AMDiS/src/io/VtkWriter.cc index fa56f4ae9737c6f2bf732adfe086d3fa6d9eac86..82c4a6be44472237d76fcd7f75ad0b9a457bf9ef 100644 --- a/AMDiS/src/io/VtkWriter.cc +++ b/AMDiS/src/io/VtkWriter.cc @@ -72,6 +72,8 @@ namespace AMDiS { { FUNCNAME("VtkWriter::writeParallelFile()"); + using boost::lexical_cast; + boost::iostreams::filtering_ostream file; { ofstream swapfile(name.c_str(), ios::out | ios::trunc); @@ -102,13 +104,11 @@ namespace AMDiS { file << " \n"; for (int i = 0; i < nRanks; i++) { - stringstream oss; - oss << fnPrefix << "-p" << i << "-" << fnPostfix; - boost::filesystem::path filepath(oss.str()); + string pname(fnPrefix + "-p" + lexical_cast(i) + "-" + fnPostfix); + boost::filesystem::path filepath(pname); file << " \n"; - + << boost::filesystem::extension(filepath) << "\"/>\n"; } file << " \n"; diff --git a/AMDiS/src/parallel/MatrixNnzStructure.cc b/AMDiS/src/parallel/MatrixNnzStructure.cc index d978f5672fec4c5d647343680d9e2cb4e7e447e6..c7d4c977f60237f3bd7aff902e965a9cdfdc654f 100644 --- a/AMDiS/src/parallel/MatrixNnzStructure.cc +++ b/AMDiS/src/parallel/MatrixNnzStructure.cc @@ -18,6 +18,12 @@ namespace AMDiS { + MatrixNnzStructure::~MatrixNnzStructure() + { + clear(); + } + + void MatrixNnzStructure::clear() { if (dnnz) { @@ -342,6 +348,8 @@ namespace AMDiS { return; TEST_EXIT_DBG(nRows0 > 0)("Should not happen!\n"); + + clear(); dnnz = new int[nRows0]; for (int i = 0; i < nRows0; i++) diff --git a/AMDiS/src/parallel/MatrixNnzStructure.h b/AMDiS/src/parallel/MatrixNnzStructure.h index 57d9a84befd70eeaf51ba074e1a65de5310e912e..cb9441e5ddcc01dfd87323b36b9b065c9ff52984 100644 --- a/AMDiS/src/parallel/MatrixNnzStructure.h +++ b/AMDiS/src/parallel/MatrixNnzStructure.h @@ -38,6 +38,8 @@ namespace AMDiS { onnz(NULL) {} + ~MatrixNnzStructure(); + void clear(); void create(Matrix &mat, diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index f95c6dc48c044b51075ee03564b25f0a7bea9ed5..13fc20291689bc0023326966815aed1e8c3c4bbe 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -130,6 +130,13 @@ namespace AMDiS { } + MeshDistributor::~MeshDistributor() + { + if (partitioner) + delete partitioner; + } + + void MeshDistributor::initParallelization() { FUNCNAME("MeshDistributor::initParallelization()"); @@ -931,6 +938,13 @@ namespace AMDiS { { FUNCNAME("MeshDistributor::checkMeshChange()"); + double vm, rss; + processMemUsage(vm, rss); + mpi::globalAdd(vm); + mpi::globalAdd(rss); + MSG("MEM-RSS 0 = %f\n", rss); + + MPI::COMM_WORLD.Barrier(); double first = MPI::Wtime(); @@ -1053,6 +1067,12 @@ namespace AMDiS { // === Print imbalance factor. === printImbalanceFactor(); + + processMemUsage(vm, rss); + mpi::globalAdd(vm); + mpi::globalAdd(rss); + MSG("MEM-RSS 1 = %f\n", rss); + } diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index c9a7328c51179087b9c25383f9c94b5e04964ebc..397eeef8b11b9510a9e3745061c4ef4b953aac21 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -60,9 +60,9 @@ namespace AMDiS { private: MeshDistributor(); - virtual ~MeshDistributor() {} - public: + ~MeshDistributor(); + void initParallelization(); void exitParallelization(); @@ -122,7 +122,7 @@ namespace AMDiS { /// leaf elements. void setInitialElementWeights(); - inline virtual string getName() + inline string getName() { return name; } diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc index e58d33f1fc381097ad7257ccfb810dd0a13a2d86..d08f82332930787df190d6667b9561aaa114ea8a 100644 --- a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc +++ b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc @@ -76,10 +76,10 @@ namespace AMDiS { TEST_EXIT(uniqueCoarseMap.size() <= 2) ("Not yet implemented for more than two coarse spaces!\n"); if (uniqueCoarseMap.size() == 2) { -// if (uniqueCoarseMap[0]->getBasisFcts()->getDegree() < -// uniqueCoarseMap[1]->getBasisFcts()->getDegree()) { + if (uniqueCoarseMap[0]->getFeSpace()->getBasisFcts()->getDegree() < + uniqueCoarseMap[1]->getFeSpace()->getBasisFcts()->getDegree()) { swap(uniqueCoarseMap[0], uniqueCoarseMap[1]); - // } + } } // === Create pointers to PETSc matrix and vector objects. === @@ -128,6 +128,7 @@ namespace AMDiS { // === If required, recompute non zero structure of the matrix. === bool localMatrix = (coarseSpaceMap.size() && subdomainLevel == 0); + if (checkMeshChange()) { // Mesh has been changed, recompute interior DOF mapping. vector feSpaces = getComponentFeSpaces(seqMat); @@ -178,7 +179,7 @@ namespace AMDiS { } else { MatCreateAIJ(mpiCommGlobal, nRankInteriorRows, nRankInteriorRows, nOverallInteriorRows, nOverallInteriorRows, - 0, nnz[0][0].dnnz, 0, nnz[0][0].onnz, + 0, nnz[0][0].dnnz, 0, nnz[0][0].onnz, &mat[0][0]); MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); } diff --git a/AMDiS/src/parallel/ParallelDofMapping.h b/AMDiS/src/parallel/ParallelDofMapping.h index e609dea867ea91f265983b8bb88f014232423c36..da79892c0f98f58e2a4e3d577e484bc512f698e2 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.h +++ b/AMDiS/src/parallel/ParallelDofMapping.h @@ -461,6 +461,12 @@ namespace AMDiS { return dofToMatIndex.get(ithComponent, d) - rStartDofs; } + inline const FiniteElemSpace* getFeSpace(int i = 0) + { + TEST_EXIT_DBG(i < feSpacesUnique.size())("Wrong component number!\n"); + return feSpacesUnique[i]; + } + // Writes all data of this object to an output stream. void serialize(ostream &out) { @@ -495,6 +501,11 @@ namespace AMDiS { VecCreateMPI(mpiComm, getRankDofs(), nGlobalRows, &vec); } + inline void createLocalVec(Vec &vec) + { + VecCreateSeq(PETSC_COMM_SELF, getRankDofs(), &vec); + } + protected: /// Insert a new FE space DOF mapping for a given FE space. void addFeSpace(const FiniteElemSpace* feSpace); diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index f5ce440ea1069fd18a4892ca95606bbaa2de68d8..7a8200c1c76c6858d94a0ff748adfc2885a9368d 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -34,7 +34,8 @@ namespace AMDiS { rStartInterior(0), nGlobalOverallInterior(0), printTimings(false), - enableStokesMode(false) + stokesMode(false), + pressureComponent(-1) { FUNCNAME("PetscSolverFeti::PetscSolverFeti()"); @@ -78,12 +79,16 @@ namespace AMDiS { vector& uniqueFe = meshDistributor->getFeSpaces(); - enableStokesMode = false; - Parameters::get("parallel->stokes mode", enableStokesMode); - if (enableStokesMode) { + stokesMode = false; + Parameters::get("parallel->feti->stokes mode", stokesMode); + if (stokesMode) { + Parameters::get("parallel->feti->pressure component", pressureComponent); + TEST_EXIT(pressureComponent >= 0) + ("FETI-DP in Stokes mode, no pressure component defined!\n"); + MSG("========= !!!! SUPER WARNING !!! ========\n"); MSG("Use make use of stokes mode which is work in progress! We guarantee for nothing and even less than this!\n"); - fullInterface = feSpaces[0]; + pressureFeSpace = feSpaces[pressureComponent]; } if (subdomain == NULL) { @@ -105,7 +110,7 @@ namespace AMDiS { localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0); lagrangeMap.init(levelData, feSpaces, uniqueFe); - if (fullInterface != NULL) + if (stokesMode) interfaceDofMap.init(levelData, feSpaces, uniqueFe); if (fetiPreconditioner == FETI_DIRICHLET) { @@ -151,7 +156,7 @@ namespace AMDiS { else localDofMap.setDofComm(meshDistributor->getDofCommSd()); - if (fullInterface != NULL) { + if (stokesMode) { interfaceDofMap.clear(); interfaceDofMap.setDofComm(meshDistributor->getDofComm()); interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0); @@ -175,7 +180,7 @@ namespace AMDiS { if (fetiPreconditioner == FETI_DIRICHLET) interiorDofMap.update(); - if (fullInterface != NULL) + if (stokesMode) interfaceDofMap.update(); for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { @@ -211,9 +216,9 @@ namespace AMDiS { MSG("FETI-DP data created on mesh level %d\n", meshLevel); for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); - MSG("FETI-DP data for %d-ith FE space:\n", i); + MSG("FETI-DP data for %d-ith FE space: %p\n", i, feSpace); - if (feSpace == fullInterface) { + if (feSpace == pressureFeSpace) { MSG(" nRankInterface = %d nOverallInterface = %d\n", interfaceDofMap[feSpace].nRankDofs, interfaceDofMap[feSpace].nOverallDofs); @@ -237,19 +242,9 @@ namespace AMDiS { } subdomain->setDofMapping(&localDofMap); - - if (enableStokesMode) { - MSG("WARNING: MAKE THIS MORE GENERAL!!!!!\n"); - subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, 0); - // subdomain->setCoarseSpaceDofMapping(&primalDofMap, 0); - subdomain->setCoarseSpaceDofMapping(&primalDofMap, 1); - subdomain->setCoarseSpaceDofMapping(&primalDofMap, 2); - subdomain->setCoarseSpaceDofMapping(&primalDofMap, 3); - subdomain->setCoarseSpaceDofMapping(&primalDofMap, 4); - - } else { - subdomain->setCoarseSpaceDofMapping(&primalDofMap); - } + subdomain->setCoarseSpaceDofMapping(&primalDofMap); + if (stokesMode) + subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent); if (printTimings) { MPI::COMM_WORLD.Barrier(); @@ -264,7 +259,7 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::createPrimals()"); - if (feSpace == fullInterface) + if (feSpace == pressureFeSpace) return; // === Define all vertices on the interior boundaries of the macro mesh === @@ -309,7 +304,7 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::createDuals()"); - if (feSpace == fullInterface) + if (feSpace == pressureFeSpace) return; // === Create global index of the dual nodes on each rank. === @@ -333,7 +328,7 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::createInterfaceNodes()"); - if (feSpace != fullInterface) + if (feSpace != pressureFeSpace) return; DofContainer allBoundaryDofs; @@ -352,7 +347,7 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::createLagrange()"); - if (feSpace == fullInterface) + if (feSpace == pressureFeSpace) return; boundaryDofRanks[feSpace].clear(); @@ -776,7 +771,7 @@ namespace AMDiS { lagrangeMap.createVec(fetiData.tmp_vec_lagrange); primalDofMap.createVec(fetiData.tmp_vec_primal0); - if (enableStokesMode == false) { + if (stokesMode == false) { MatCreateShell(mpiCommGlobal, lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(), @@ -807,42 +802,26 @@ namespace AMDiS { KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000); KSPSetFromOptions(ksp_feti); - if (enableStokesMode) { - Vec nullSpaceInterface; - Vec nullSpaceLagrange; - Vec nullSpaceBasis; - - interfaceDofMap.createVec(nullSpaceInterface); - VecSet(nullSpaceInterface, 1.0); - - lagrangeMap.createVec(nullSpaceLagrange); - VecSet(nullSpaceLagrange, 0.0); - - Vec vecArray[2] = {nullSpaceInterface, nullSpaceLagrange}; - VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecArray, &nullSpaceBasis); - - - MatNullSpace matNullspace; - MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullspace); - MatSetNullSpace(mat_feti, matNullspace); - KSPSetNullSpace(ksp_feti, matNullspace); - } - // === Create FETI-DP preconditioner object. === - if (fetiPreconditioner != FETI_NONE) { + if (fetiPreconditioner != FETI_NONE || stokesMode) { MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled); MatScale(mat_lagrange_scaled, 0.5); } + // === Create null space objects. === + + createNullSpace(); + + switch (fetiPreconditioner) { case FETI_DIRICHLET: TEST_EXIT(meshLevel == 0) ("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n"); - TEST_EXIT(!enableStokesMode) + TEST_EXIT(!stokesMode) ("Stokes mode does not yet support the Dirichlet precondition!\n"); KSPCreate(PETSC_COMM_SELF, &ksp_interior); @@ -905,34 +884,48 @@ namespace AMDiS { break; case FETI_LUMPED: - fetiLumpedPreconData.enableStokesMode = enableStokesMode; - fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled; - fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals; - - for (unsigned int i = 0; i < feSpaces.size(); i++) { - if (enableStokesMode && feSpaces[i] == fullInterface) - continue; + { + FetiLumpedPreconData *lumpedData = + (stokesMode ? &fetiInterfaceLumpedPreconData : &fetiLumpedPreconData); - DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(); - for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { - DegreeOfFreedom d = it->first; - int matIndexLocal = localDofMap.getLocalMatIndex(i, d); - int matIndexDual = dualDofMap.getLocalMatIndex(i, d); - fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual; + lumpedData->mat_lagrange_scaled = &mat_lagrange_scaled; + lumpedData->mat_duals_duals = &mat_duals_duals; + + localDofMap.createVec(lumpedData->tmp_vec_b0); + MatGetVecs(mat_duals_duals, PETSC_NULL, + &(lumpedData->tmp_vec_duals0)); + MatGetVecs(mat_duals_duals, PETSC_NULL, + &(lumpedData->tmp_vec_duals1)); + + for (unsigned int i = 0; i < feSpaces.size(); i++) { + if (stokesMode && i == pressureComponent) + continue; + + DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(); + for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { + DegreeOfFreedom d = it->first; + int matIndexLocal = localDofMap.getLocalMatIndex(i, d); + int matIndexDual = dualDofMap.getLocalMatIndex(i, d); + lumpedData->localToDualMap[matIndexLocal] = matIndexDual; + } + } + + if (stokesMode) { + localDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_vec_b1); + fetiInterfaceLumpedPreconData.subSolver = subdomain; + } + + KSPGetPC(ksp_feti, &precon_feti); + PCSetType(precon_feti, PCSHELL); + if (stokesMode) { + PCShellSetContext(precon_feti, static_cast(&fetiInterfaceLumpedPreconData)); + PCShellSetApply(precon_feti, petscApplyFetiInterfaceLumpedPrecon); + } else { + PCShellSetContext(precon_feti, static_cast(&fetiLumpedPreconData)); + PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon); } } - localDofMap.createVec(fetiLumpedPreconData.tmp_vec_b); - MatGetVecs(mat_duals_duals, PETSC_NULL, - &(fetiLumpedPreconData.tmp_vec_duals0)); - MatGetVecs(mat_duals_duals, PETSC_NULL, - &(fetiLumpedPreconData.tmp_vec_duals1)); - - KSPGetPC(ksp_feti, &precon_feti); - PCSetType(precon_feti, PCSHELL); - PCShellSetContext(precon_feti, static_cast(&fetiLumpedPreconData)); - PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon); - break; default: break; @@ -954,7 +947,7 @@ namespace AMDiS { VecDestroy(&fetiData.tmp_vec_lagrange); VecDestroy(&fetiData.tmp_vec_primal0); - if (enableStokesMode) { + if (stokesMode) { VecDestroy(&fetiData.tmp_vec_b1); VecDestroy(&fetiData.tmp_vec_primal1); VecDestroy(&fetiData.tmp_vec_interface); @@ -985,12 +978,17 @@ namespace AMDiS { break; case FETI_LUMPED: - fetiLumpedPreconData.mat_lagrange_scaled = NULL; - fetiLumpedPreconData.mat_duals_duals = NULL; - - VecDestroy(&fetiLumpedPreconData.tmp_vec_b); - VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0); - VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1); + { + FetiLumpedPreconData &lumpedData = + (stokesMode ? fetiInterfaceLumpedPreconData : fetiLumpedPreconData); + + lumpedData.mat_lagrange_scaled = NULL; + lumpedData.mat_duals_duals = NULL; + + VecDestroy(&lumpedData.tmp_vec_b0); + VecDestroy(&lumpedData.tmp_vec_duals0); + VecDestroy(&lumpedData.tmp_vec_duals1); + } break; default: break; @@ -998,6 +996,76 @@ namespace AMDiS { } + void PetscSolverFeti::createNullSpace() + { + FUNCNAME("PetscSolverFeti::createNullSpace()"); + + if (!stokesMode) + return; + + Vec ktest0, ktest1; + localDofMap.createLocalVec(ktest0); + localDofMap.createLocalVec(ktest1); + DofMap& m = localDofMap[pressureFeSpace].getMap(); + for (DofMap::iterator it = m.begin(); it != m.end(); ++it) { + if (meshDistributor->getDofMap()[pressureFeSpace].isRankDof(it->first)) { + int index = localDofMap.getLocalMatIndex(pressureComponent, it->first); + VecSetValue(ktest0, index, 1.0, INSERT_VALUES); + } + } + VecAssemblyBegin(ktest0); + VecAssemblyEnd(ktest0); + MatMult(subdomain->getMatInterior(), ktest0, ktest1); + + PetscScalar *valarray; + VecGetArray(ktest1, &valarray); + + Vec ktest2, ktest3; + VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, + localDofMap.getRankDofs(), nGlobalOverallInterior, + valarray, &ktest2); + localDofMap.createVec(ktest3, nGlobalOverallInterior); + + Vec vecArray[2]; + interfaceDofMap.createVec(vecArray[0]); + lagrangeMap.createVec(vecArray[1]); + + VecSet(vecArray[0], 1.0); + MatMult(subdomain->getMatInteriorCoarse(1), vecArray[0], ktest3); + VecAXPY(ktest3, 1.0, ktest2); + MatMult(mat_lagrange_scaled, ktest3, vecArray[1]); + VecScale(vecArray[1], -1.0); + + Vec nullSpaceBasis; + VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecArray, &nullSpaceBasis); + + MatNullSpace matNullSpace; + MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, + &matNullSpace); + MatSetNullSpace(mat_feti, matNullSpace); + KSPSetNullSpace(ksp_feti, matNullSpace); + MatNullSpaceDestroy(&matNullSpace); + +#if (DEBUG != 0) + Vec vecSol; + VecDuplicate(nullSpaceBasis, &vecSol); + MatMult(mat_feti, nullSpaceBasis, vecSol); + PetscReal norm; + VecNorm(vecSol, NORM_2, &norm); + MSG("NORM: %e\n", norm); +#endif + + VecDestroy(&ktest0); + VecDestroy(&ktest1); + VecDestroy(&ktest2); + VecDestroy(&ktest3); + + VecDestroy(&(vecArray[0])); + VecDestroy(&(vecArray[1])); + VecDestroy(&nullSpaceBasis); + } + + void PetscSolverFeti::dbgMatrix() { FUNCNAME("PetscSolverFeti::dbgMatrix()"); @@ -1029,7 +1097,7 @@ namespace AMDiS { MSG("Interior-Primal matrix [%d x %d] nnz = %d\n", nRow, nCol, static_cast(minfo.nz_used)); - if (enableStokesMode) { + if (stokesMode) { MatGetSize(subdomain->getMatCoarse(1, 1), &nRow, &nCol); MatGetInfo(subdomain->getMatCoarse(1, 1), MAT_GLOBAL_SUM, &minfo); MSG("Interface matrix [%d x %d] nnz = %d\n", nRow, nCol, @@ -1240,7 +1308,7 @@ namespace AMDiS { VecGetArray(local_sol_primal, &localSolPrimal); // === And copy from PETSc local vectors to the DOF vectors. === - + int cnt = 0; for (int i = 0; i < vec.getSize(); i++) { DOFVector& dofVec = *(vec.getDOFVector(i)); @@ -1274,18 +1342,20 @@ namespace AMDiS { { FUNCNAME("PetscSolverFeti::recoverInterfaceSolution()"); - if (!enableStokesMode) + if (!stokesMode) return; vector globalIsIndex, localIsIndex; globalIsIndex.reserve(interfaceDofMap.getLocalDofs()); localIsIndex.reserve(interfaceDofMap.getLocalDofs()); - MSG("USDED FIXED PRESSURE VARIABLE 0!!!\n"); + const FiniteElemSpace *pressureFeSpace = + vec.getDOFVector(pressureComponent)->getFeSpace(); int cnt = 0; - DofMap& dofMap = interfaceDofMap[fullInterface].getMap(); + DofMap& dofMap = interfaceDofMap[pressureFeSpace].getMap(); for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) { - globalIsIndex.push_back(interfaceDofMap.getMatIndex(0, it->first)); + globalIsIndex.push_back(interfaceDofMap.getMatIndex(pressureComponent, + it->first)); localIsIndex.push_back(cnt++); } @@ -1323,10 +1393,11 @@ namespace AMDiS { // === And copy from PETSc local vectors to the DOF vectors. === cnt = 0; - DOFVector& dofVec = *(vec.getDOFVector(2)); - for (DofMap::iterator it = interfaceDofMap[fullInterface].getMap().begin(); - it != interfaceDofMap[fullInterface].getMap().end(); ++it) - dofVec[it->first] = localSolInterface[cnt++]; + DOFVector& dofVec = *(vec.getDOFVector(pressureComponent)); + for (DofMap::iterator it = interfaceDofMap[pressureFeSpace].getMap().begin(); + it != interfaceDofMap[pressureFeSpace].getMap().end(); ++it) { + dofVec[it->first] = localSolInterface[cnt++]; + } VecRestoreArray(local_sol_interface, &localSolInterface); VecDestroy(&local_sol_interface); @@ -1462,12 +1533,14 @@ namespace AMDiS { if (!dofMat) continue; + if (stokesMode && + (rowComponent == pressureComponent || + colComponent == pressureComponent)) + continue; + const FiniteElemSpace *rowFeSpace = feSpaces[rowComponent]; const FiniteElemSpace *colFeSpace = feSpaces[colComponent]; - if (enableStokesMode && (rowFeSpace == fullInterface || colFeSpace == fullInterface)) - continue; - traits::col::type col(dofMat->getBaseMatrix()); traits::const_value::type value(dofMat->getBaseMatrix()); @@ -1636,7 +1709,7 @@ namespace AMDiS { Vec vecRhs, vecSol; Vec vecRhsLagrange, vecSolLagrange, vecRhsInterface, vecSolInterface; - if (!enableStokesMode) { + if (!stokesMode) { MatGetVecs(mat_lagrange, PETSC_NULL, &vecRhsLagrange); MatGetVecs(mat_lagrange, PETSC_NULL, &vecSolLagrange); @@ -1669,7 +1742,7 @@ namespace AMDiS { double wtime = MPI::Wtime(); - if (!enableStokesMode) { + if (!stokesMode) { // d = L inv(K_BB) f_B - L inv(K_BB) K_BPi inv(S_PiPi) [f_Pi - K_PiB inv(K_BB) f_B] // vecRhs = L * inv(K_BB) * f_B subdomain->solveGlobal(subdomain->getVecRhsInterior(), tmp_b0); @@ -1746,7 +1819,7 @@ namespace AMDiS { // === Solve for u_primals. === - if (!enableStokesMode) { + if (!stokesMode) { MatMultTranspose(mat_lagrange, vecSol, tmp_b0); VecAYPX(tmp_b0, -1.0, subdomain->getVecRhsInterior()); } else { @@ -1761,14 +1834,13 @@ namespace AMDiS { MatMult(subdomain->getMatCoarseInterior(), tmp_b1, tmp_primal0); VecAYPX(tmp_primal0, -1.0, subdomain->getVecRhsCoarse()); - if (enableStokesMode) { + if (stokesMode) { MatMult(subdomain->getMatCoarse(0, 1), vecSolInterface, tmp_primal1); VecAXPY(tmp_primal0, -1.0, tmp_primal1); } KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0); - // === Solve for u_b. === MatMult(subdomain->getMatInteriorCoarse(), tmp_primal0, tmp_b1); @@ -1779,7 +1851,7 @@ namespace AMDiS { // === And recover AMDiS solution vectors. === recoverSolution(tmp_b0, tmp_primal0, vec); - if (enableStokesMode) + if (stokesMode) recoverInterfaceSolution(vecSolInterface, vec); @@ -1799,7 +1871,7 @@ namespace AMDiS { VecDestroy(&tmp_primal0); VecDestroy(&tmp_primal1); - if (enableStokesMode) { + if (stokesMode) { VecDestroy(&tmp_interface); VecDestroy(&vecRhsInterface); VecDestroy(&vecSolInterface); @@ -1855,4 +1927,94 @@ namespace AMDiS { MeshDistributor::globalMeshDistributor->synchVector(vec); } + + void PetscSolverFeti::createInteriorMat(Mat &mat) + { + FUNCNAME("PetscSolverFeti::createInteriorMat()"); + + Mat &localMat = subdomain->getMatInterior(); + int nRow, nCol; + MatGetSize(localMat, &nRow, &nCol); + + TEST_EXIT_DBG(nRow == nCol)("Should not happen!\n"); + + int *dnnz = new int[nRow]; + for (int i = 0; i < nRow; i++) { + MatGetRow(localMat, i, &nCol, PETSC_NULL, PETSC_NULL); + dnnz[i] = nCol; + MatRestoreRow(localMat, i, &nCol, PETSC_NULL, PETSC_NULL); + } + + MatCreateAIJ(PETSC_COMM_WORLD, + nRow, nRow, + nGlobalOverallInterior, nGlobalOverallInterior, + 0, dnnz, 0, PETSC_NULL, &mat); + + PetscInt rStart, rEnd; + MatGetOwnershipRange(mat, &rStart, &rEnd); + + for (int i = 0; i < nRow; i++) { + const PetscInt *cols; + const PetscScalar *vals; + + MatGetRow(localMat, i, &nCol, &cols, &vals); + for (int j = 0; j < nCol; j++) + MatSetValue(mat, i + rStart, cols[j] + rStart, vals[j], INSERT_VALUES); + } + + MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY); + + delete [] dnnz; + } + + + void PetscSolverFeti::createNestedFetiMat(Mat &mat) + { + FUNCNAME("PetscSolverFeti::createNestedFetiMat()"); + + if (stokesMode) { + vector nestMat; + nestMat.resize(16); + + createInteriorMat(nestMat[0]); + nestMat[1] = subdomain->getMatInteriorCoarse(0); + nestMat[2] = subdomain->getMatInteriorCoarse(1); + MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &nestMat[3]); + nestMat[4] = subdomain->getMatCoarseInterior(0); + + nestMat[5] = subdomain->getMatCoarse(0, 0); + nestMat[6] = subdomain->getMatCoarse(0, 1); + nestMat[7] = PETSC_NULL; + nestMat[8] = subdomain->getMatCoarseInterior(1); + nestMat[9] = subdomain->getMatCoarse(1, 0); + nestMat[10] = PETSC_NULL; + nestMat[11] = PETSC_NULL; + nestMat[12] = mat_lagrange; + nestMat[13] = PETSC_NULL; + nestMat[14] = PETSC_NULL; + nestMat[15] = PETSC_NULL; + + Mat nestFetiMat; + MatCreateNest(mpiCommGlobal, 4, PETSC_NULL, 4, PETSC_NULL, + &(nestMat[0]), &mat); + } else { + vector nestMat; + nestMat.resize(9); + + createInteriorMat(nestMat[0]); + nestMat[1] = subdomain->getMatInteriorCoarse(0); + MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &nestMat[2]); + nestMat[3] = subdomain->getMatCoarseInterior(0); + nestMat[4] = subdomain->getMatCoarse(0, 0); + nestMat[5] = PETSC_NULL; + nestMat[6] = mat_lagrange; + nestMat[7] = PETSC_NULL; + nestMat[8] = PETSC_NULL; + + Mat nestFetiMat; + MatCreateNest(mpiCommGlobal, 3, PETSC_NULL, 3, PETSC_NULL, + &(nestMat[0]), &mat); + } + } } diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h index 16f3c0d5e70d0d5ae0fb4dbfaec733807aab39fc..f32f281d034b615831542d557449898fd71e6062 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.h +++ b/AMDiS/src/parallel/PetscSolverFeti.h @@ -141,6 +141,10 @@ namespace AMDiS { /// Destroys FETI-DP operator, \ref ksp_feti void destroyFetiKsp(); + /// Create the null space of the FETI-DP operator (if there is one) and + /// attachets it to the corresponding matrices and KSP objects. + void createNullSpace(); + /// In debug modes, this function runs some debug tests on the FETI /// matrices. In optimized mode, nothing is done here. void dbgMatrix(); @@ -205,12 +209,18 @@ namespace AMDiS { inline bool isInterface(const FiniteElemSpace *feSpace, DegreeOfFreedom dof) { - if (feSpace == fullInterface) + if (feSpace == pressureFeSpace) return interfaceDofMap[feSpace].isSet(dof); return false; } + /// For debugging only! + void createInteriorMat(Mat &mat); + + /// For debugging only! + void createNestedFetiMat(Mat &mat); + protected: /// Mapping from primal DOF indices to a global index of primals. ParallelDofMapping primalDofMap; @@ -281,6 +291,8 @@ namespace AMDiS { FetiLumpedPreconData fetiLumpedPreconData; + FetiInterfaceLumpedPreconData fetiInterfaceLumpedPreconData; + /// Matrices for Dirichlet preconditioner. Mat mat_interior_interior, mat_duals_duals, mat_interior_duals, mat_duals_interior; @@ -296,11 +308,13 @@ namespace AMDiS { int nGlobalOverallInterior; - const FiniteElemSpace* fullInterface; - bool printTimings; - bool enableStokesMode; + bool stokesMode; + + const FiniteElemSpace* pressureFeSpace; + + int pressureComponent; }; } diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc index 664e9f84d7d17313f7f1a6a04f90a5dee380d1e3..61f86cdcad2aec885139af2e7901033893046271 100644 --- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc +++ b/AMDiS/src/parallel/PetscSolverFetiOperators.cc @@ -231,51 +231,102 @@ namespace AMDiS { } - PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y) + PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec xvec, Vec yvec) { // Get data for the preconditioner void *ctx; PCShellGetContext(pc, &ctx); FetiLumpedPreconData* data = static_cast(ctx); - Vec xvec, yvec; - if (data->enableStokesMode) { - Vec x_interface, x_lagrange, y_interface, y_lagrange; - VecNestGetSubVec(x, 0, &x_interface); - VecNestGetSubVec(x, 1, &x_lagrange); - VecNestGetSubVec(y, 0, &y_interface); - VecNestGetSubVec(y, 1, &y_lagrange); - - xvec = x_lagrange; - yvec = y_lagrange; - - VecCopy(x_interface, y_interface); - // VecScale(y_interface, 0.2519); - } else { - xvec = x; - yvec = y; - } + // Multiply with scaled Lagrange constraint matrix. + MatMultTranspose(*(data->mat_lagrange_scaled), xvec, data->tmp_vec_b0); + + + // === Restriction of the B nodes to the boundary nodes. === + + int nLocalB; + int nLocalDuals; + VecGetLocalSize(data->tmp_vec_b0, &nLocalB); + VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals); + + PetscScalar *local_b, *local_duals; + VecGetArray(data->tmp_vec_b0, &local_b); + VecGetArray(data->tmp_vec_duals0, &local_duals); + + for (map::iterator it = data->localToDualMap.begin(); + it != data->localToDualMap.end(); ++it) + local_duals[it->second] = local_b[it->first]; + + VecRestoreArray(data->tmp_vec_b0, &local_b); + VecRestoreArray(data->tmp_vec_duals0, &local_duals); + + + // === K_DD === + + MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1); + + + // === Prolongation from local dual nodes to B nodes. + + VecGetArray(data->tmp_vec_b0, &local_b); + VecGetArray(data->tmp_vec_duals1, &local_duals); + + for (map::iterator it = data->localToDualMap.begin(); + it != data->localToDualMap.end(); ++it) + local_b[it->first] = local_duals[it->second]; + + VecRestoreArray(data->tmp_vec_b0, &local_b); + VecRestoreArray(data->tmp_vec_duals0, &local_duals); + + + // Multiply with scaled Lagrange constraint matrix. + MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b0, yvec); + + return 0; + } + + + PetscErrorCode petscApplyFetiInterfaceLumpedPrecon(PC pc, Vec xvec, Vec yvec) + { + // Get data for the preconditioner + void *ctx; + PCShellGetContext(pc, &ctx); + FetiInterfaceLumpedPreconData* data = + static_cast(ctx); + + Vec x_interface, x_lagrange, y_interface, y_lagrange; + VecNestGetSubVec(xvec, 0, &x_interface); + VecNestGetSubVec(xvec, 1, &x_lagrange); + VecNestGetSubVec(yvec, 0, &y_interface); + VecNestGetSubVec(yvec, 1, &y_lagrange); + + + VecCopy(x_interface, y_interface); + // VecScale(y_interface, 0.2519); + // Multiply with scaled Lagrange constraint matrix. - MatMultTranspose(*(data->mat_lagrange_scaled), xvec, data->tmp_vec_b); + MatMult(data->subSolver->getMatInteriorCoarse(1), x_interface, data->tmp_vec_b0); + MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b1); + VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1); // === Restriction of the B nodes to the boundary nodes. === int nLocalB; int nLocalDuals; - VecGetLocalSize(data->tmp_vec_b, &nLocalB); + VecGetLocalSize(data->tmp_vec_b0, &nLocalB); VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals); PetscScalar *local_b, *local_duals; - VecGetArray(data->tmp_vec_b, &local_b); + VecGetArray(data->tmp_vec_b0, &local_b); VecGetArray(data->tmp_vec_duals0, &local_duals); for (map::iterator it = data->localToDualMap.begin(); it != data->localToDualMap.end(); ++it) local_duals[it->second] = local_b[it->first]; - VecRestoreArray(data->tmp_vec_b, &local_b); + VecRestoreArray(data->tmp_vec_b0, &local_b); VecRestoreArray(data->tmp_vec_duals0, &local_duals); @@ -286,19 +337,21 @@ namespace AMDiS { // === Prolongation from local dual nodes to B nodes. - VecGetArray(data->tmp_vec_b, &local_b); + VecGetArray(data->tmp_vec_b0, &local_b); VecGetArray(data->tmp_vec_duals1, &local_duals); for (map::iterator it = data->localToDualMap.begin(); it != data->localToDualMap.end(); ++it) local_b[it->first] = local_duals[it->second]; - VecRestoreArray(data->tmp_vec_b, &local_b); + VecRestoreArray(data->tmp_vec_b0, &local_b); VecRestoreArray(data->tmp_vec_duals0, &local_duals); // Multiply with scaled Lagrange constraint matrix. - MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, yvec); + MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b0, y_lagrange); + + // MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface); return 0; } diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.h b/AMDiS/src/parallel/PetscSolverFetiOperators.h index 6b0fa999a83a30e5769cda7cba60856c3d211217..e305b25c40e08299ef2eef2a9353c1881179c143 100644 --- a/AMDiS/src/parallel/PetscSolverFetiOperators.h +++ b/AMDiS/src/parallel/PetscSolverFetiOperators.h @@ -42,6 +42,8 @@ namespace AMDiS { /// y = PC * x PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y); + + PetscErrorCode petscApplyFetiInterfaceLumpedPrecon(PC pc, Vec x, Vec y); } #endif diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h index 2de8c3056fde4b0c6cceec9b6a47a399ad22832e..b8dd365fe2f04d5805aa3a2d308c9afcb18aa235 100644 --- a/AMDiS/src/parallel/PetscSolverFetiStructs.h +++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h @@ -99,17 +99,13 @@ namespace AMDiS { struct FetiLumpedPreconData { - /// Defines whether the precondition is applied for a Stokes like FETI-DP - /// problem with interface variables which handled in a different way. - bool enableStokesMode; - /// Matrix of scaled Lagrange variables. Mat *mat_lagrange_scaled; Mat *mat_duals_duals; /// Temporal vector on the B variables. - Vec tmp_vec_b; + Vec tmp_vec_b0; /// Temporal vector on the dual variables. Vec tmp_vec_duals0, tmp_vec_duals1; @@ -118,6 +114,14 @@ namespace AMDiS { }; + struct FetiInterfaceLumpedPreconData : public FetiLumpedPreconData { + /// Temporal vectors on the B variables. + Vec tmp_vec_b1; + + PetscSolver* subSolver; + }; + + typedef enum { FETI_NONE = 0, FETI_DIRICHLET = 1, diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc index 756a7a5bc5f1d388c40fe905a0fe733c20b9801a..55e96dacf0309364e68eda993f30ea60b5efcca6 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc @@ -252,7 +252,7 @@ namespace AMDiS { if (coarseSpaceMap.size()) { for (int i = 0; i < vec->getSize(); i++) setDofVector(getVecRhsInterior(), - getVecSolCoarseByComponent(i), vec->getDOFVector(i), i); + getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i); } else { for (int i = 0; i < vec->getSize(); i++) setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i); @@ -379,12 +379,9 @@ namespace AMDiS { Vec tmp; if (mpiCommLocal.Get_size() == 1) - VecCreateSeq(mpiCommLocal, interiorMap->getRankDofs(), &tmp); + interiorMap->createLocalVec(tmp); else - VecCreateMPI(mpiCommLocal, - interiorMap->getRankDofs(), - interiorMap->getOverallDofs(), - &tmp); + interiorMap->createVec(tmp); PetscScalar *tmpValues, *rhsValues; VecGetArray(tmp, &tmpValues);