Commit a6bacc8a by Thomas Witkowski

### Work on repartitioning with higher order finite elements.

parent 79bcb0ed
 // ============================================================================ // == == // == AMDiS - Adaptive multidimensional simulations == // == == // == http://www.amdis-fem.org == // == == // ============================================================================ // // Software License for AMDiS // // Copyright (c) 2010 Dresden University of Technology // All rights reserved. // Authors: Simon Vey, Thomas Witkowski et al. // // This file is part of AMDiS // // See also license.opensource.txt in the distribution. /** \file Containers.h */ #ifndef AMDIS_CONTAINERS_H #define AMDIS_CONTAINERS_H #include #include namespace AMDiS { /* template class amdis_map : public std::map { public: inline void insert(Key &key, T &v) { std::map::insert(std::pair(key, v)); } inline T& operator[](const Key& key) { return std::map::operator[](key); } }; */ } #endif
 ... @@ -627,12 +627,13 @@ namespace AMDiS { ... @@ -627,12 +627,13 @@ namespace AMDiS { void Element::getAllDofs(const FiniteElemSpace* feSpace, void Element::getAllDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, BoundaryObject bound, DofContainer& dofs) DofContainer& dofs, bool baseDofPtr) { { getNodeDofs(feSpace, bound, dofs); getNodeDofs(feSpace, bound, dofs, baseDofPtr); if (feSpace->getBasisFcts()->getDegree() > 1) if (feSpace->getBasisFcts()->getDegree() > 1) getHigherOrderDofs(feSpace, bound, dofs); getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr); } } } }
 ... @@ -82,10 +82,8 @@ namespace AMDiS { ... @@ -82,10 +82,8 @@ namespace AMDiS { /// /// void deleteElementDOFs(); void deleteElementDOFs(); /** \brief /// Clone this Element and return a reference to it. Because also the DOFs * Clone this Element and return a reference to it. Because also the DOFs /// are cloned, \ref Mesh::serializedDOfs must be used. * are cloned, \ref Mesh::serializedDOfs must be used. */ Element* cloneWithDOFs(); Element* cloneWithDOFs(); /** \name getting methods /** \name getting methods ... @@ -111,10 +109,8 @@ namespace AMDiS { ... @@ -111,10 +109,8 @@ namespace AMDiS { return child[i]; return child[i]; } } /** \brief /// Returns true if Element is a leaf element (\ref child[0] == NULL), returns * Returns true if Element is a leaf element (\ref child[0] == NULL), returns /// false otherwise. * false otherwise. */ inline const bool isLeaf() const inline const bool isLeaf() const { { return (child[0] == NULL); return (child[0] == NULL); ... @@ -155,10 +151,8 @@ namespace AMDiS { ... @@ -155,10 +151,8 @@ namespace AMDiS { dofValid = b; dofValid = b; } } /** \brief /// Returns \ref elementData's error estimation, if Element is a leaf element * Returns \ref elementData's error estimation, if Element is a leaf element /// and has leaf data. * and has leaf data. */ inline double getEstimation(int row) const inline double getEstimation(int row) const { { if (isLeaf()) { if (isLeaf()) { ... @@ -172,10 +166,8 @@ namespace AMDiS { ... @@ -172,10 +166,8 @@ namespace AMDiS { return 0.0; return 0.0; } } /** \brief /// Returns Element's coarsening error estimation, if Element is a leaf * Returns Element's coarsening error estimation, if Element is a leaf /// element and if it has leaf data and if this leaf data are coarsenable. * element and if it has leaf data and if this leaf data are coarsenable. */ inline double getCoarseningEstimation(int row) inline double getCoarseningEstimation(int row) { { if (isLeaf()) { if (isLeaf()) { ... @@ -195,10 +187,8 @@ namespace AMDiS { ... @@ -195,10 +187,8 @@ namespace AMDiS { /// Returns local vertex number of the j-th vertex of the i-th edge /// Returns local vertex number of the j-th vertex of the i-th edge virtual int getVertexOfEdge(int i, int j) const = 0; virtual int getVertexOfEdge(int i, int j) const = 0; /** \brief /// Returns local vertex number of the vertexIndex-th vertex of the * Returns local vertex number of the vertexIndex-th vertex of the /// positionIndex-th part of type position (vertex, edge, face) * positionIndex-th part of type position (vertex, edge, face) */ virtual int getVertexOfPosition(GeoIndex position, virtual int getVertexOfPosition(GeoIndex position, int positionIndex, int positionIndex, int vertexIndex) const = 0; int vertexIndex) const = 0; ... @@ -215,10 +205,6 @@ namespace AMDiS { ... @@ -215,10 +205,6 @@ namespace AMDiS { /// /// virtual DofFace getFace(int localFaceIndex) const = 0; virtual DofFace getFace(int localFaceIndex) const = 0; /// Returns either vertex, edge or face of the element. /* template */ /* T getObject(int index); */ /// Returns the number of parts of type i in this element /// Returns the number of parts of type i in this element virtual int getGeo(GeoIndex i) const = 0; virtual int getGeo(GeoIndex i) const = 0; ... @@ -271,10 +257,8 @@ namespace AMDiS { ... @@ -271,10 +257,8 @@ namespace AMDiS { elementData = ed; elementData = ed; } } /** \brief /// Sets \ref newCoord of Element. Needed by refinement, if Element has a * Sets \ref newCoord of Element. Needed by refinement, if Element has a /// boundary edge on a curved boundary. * boundary edge on a curved boundary. */ inline void setNewCoord(WorldVector* coord) inline void setNewCoord(WorldVector* coord) { { newCoord = coord; newCoord = coord; ... @@ -292,10 +276,8 @@ namespace AMDiS { ... @@ -292,10 +276,8 @@ namespace AMDiS { dof[pos] = p; dof[pos] = p; } } /** \brief /// Checks whether Element is a leaf element and whether it has leaf data. * Checks whether Element is a leaf element and whether it has leaf data. /// If the checks don't fail, leaf data's error estimation is set to est. * If the checks don't fail, leaf data's error estimation is set to est. */ inline void setEstimation(double est, int row) inline void setEstimation(double est, int row) { { FUNCNAME("Element::setEstimation()"); FUNCNAME("Element::setEstimation()"); ... @@ -312,10 +294,8 @@ namespace AMDiS { ... @@ -312,10 +294,8 @@ namespace AMDiS { } } } } /** \brief /// Sets Element's coarsening error estimation, if Element is a leaf element * Sets Element's coarsening error estimation, if Element is a leaf element /// and if it has leaf data and if this leaf data are coarsenable. * and if it has leaf data and if this leaf data are coarsenable. */ inline void setCoarseningEstimation(double est, int row) inline void setCoarseningEstimation(double est, int row) { { if (isLeaf()) { if (isLeaf()) { ... @@ -363,18 +343,19 @@ namespace AMDiS { ... @@ -363,18 +343,19 @@ namespace AMDiS { virtual const FixVec& virtual const FixVec& sortFaceIndices(int face, FixVec *vec) const = 0; sortFaceIndices(int face, FixVec *vec) const = 0; /// Returns a copy of itself. Needed by Mesh to create Elements by a prototype. /// Returns a copy of itself. Needed by Mesh to create Elements by /// a prototype. virtual Element *clone() = 0; virtual Element *clone() = 0; /** \brief /// Returns which side of child[childnr] corresponds to side sidenr of * Returns which side of child[childnr] corresponds to side sidenr of this /// this Element. If the child has no corresponding side, the return value * Element. If the child has no corresponding side, the return value is negative. /// is negative. */ virtual int getSideOfChild(int childnr, int sidenr, int elType = 0) const = 0; virtual int getSideOfChild(int childnr, int sidenr, int elType = 0) const = 0; /** \brief /** \brief * Generalization of \ref getSideOfChild to arbitrary subObject. Thus, e.g., in 3d * Generalization of \ref getSideOfChild to arbitrary subObject. Thus, * we can ask for the local id of a verte, edge or face on the elements children. * e.g., in 3d we can ask for the local id of a verte, edge or face * on the elements children. * * * \param[in] childnr Either 0 or 1 for the left or right children. * \param[in] childnr Either 0 or 1 for the left or right children. * \param[in] subObj Defines whether we ask for VERTEX, EDGE or FACE. * \param[in] subObj Defines whether we ask for VERTEX, EDGE or FACE. ... @@ -384,12 +365,12 @@ namespace AMDiS { ... @@ -384,12 +365,12 @@ namespace AMDiS { virtual int getSubObjOfChild(int childnr, GeoIndex subObj, int ithObj, virtual int getSubObjOfChild(int childnr, GeoIndex subObj, int ithObj, int elType = 0) const = 0; int elType = 0) const = 0; /** \brief /// Returns which vertex of elements parent corresponds to the vertexnr of * Returns which vertex of elements parent corresponds to the vertexnr of /// the element, if the element is the childnr-th child of the parent. * the element, if the element is the childnr-th child of the parent. /// If the vertex is the ner vertex at the refinement edge, -1 is returned. * If the vertex is the ner vertex at the refinement edge, -1 is returned. virtual int getVertexOfParent(int childnr, */ int vertexnr, virtual int getVertexOfParent(int childnr, int vertexnr, int elType = 0) const = 0; int elType = 0) const = 0; /// Returns whether Element is a Line /// Returns whether Element is a Line virtual bool isLine() const = 0; virtual bool isLine() const = 0; ... @@ -416,14 +397,19 @@ namespace AMDiS { ... @@ -416,14 +397,19 @@ namespace AMDiS { * nodes alonge this vertex/edge/face are assembled and put together to * nodes alonge this vertex/edge/face are assembled and put together to * a list. * a list. * * * \param[in] feSpace FE space which is used to get the dofs. * \param[in] feSpace FE space which is used to get the dofs. * \param[in] bound Defines the vertex/edge/face of the element on * \param[in] bound Defines the vertex/edge/face of the element on * which all vertex dofs are assembled. * which all vertex dofs are assembled. * \param[out] dofs List of dofs, where the result is stored. * \param[out] dofs List of dofs, where the result is stored. * \param[in] baseDofPtr If true, the base DOF pointes are stored. Thus, * dof* [\ref dof] of the element is inserted. If * false, &(dof[.][n0]) is put to the result vector, * with n0 beging the number of predofs. */ */ virtual void getNodeDofs(const FiniteElemSpace* feSpace, virtual void getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, BoundaryObject bound, DofContainer& dofs) const = 0; DofContainer& dofs, bool baseDofPtr = false) const = 0; /** \brief /** \brief * Traverses a vertex/edge/face of a given element (this includes also all * Traverses a vertex/edge/face of a given element (this includes also all ... @@ -431,18 +417,26 @@ namespace AMDiS { ... @@ -431,18 +417,26 @@ namespace AMDiS { * to higher order basis functions alonge this vertex/edge/face are * to higher order basis functions alonge this vertex/edge/face are * assembled and put together to a list. * assembled and put together to a list. * * * \param[in] feSpace FE space which is used to get the dofs. * \param[in] feSpace FE space which is used to get the dofs. * \param[in] bound Defines the edge/face of the element on which * \param[in] bound Defines the edge/face of the element on which * all non vertex dofs are assembled. * all non vertex dofs are assembled. * \param[out] dofs All dofs are put to this dof list. * \param[out] dofs All dofs are put to this dof list. * \param[in] baseDofPtr If true, the base DOF pointes are stored. Thus, * dof* [\ref dof] of the element is inserted. If * false, &(dof[.][n0]) is put to the result vector, * with n0 beging the number of predofs. */ */ virtual void getHigherOrderDofs(const FiniteElemSpace* feSpace, virtual void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, BoundaryObject bound, DofContainer& dofs) const = 0; DofContainer& dofs, bool baseDofPtr = false) const = 0; /// Combines \ref getNodeDofs and \ref getHigherOrderDofs to one function. /// See parameter description there. void getAllDofs(const FiniteElemSpace* feSpace, void getAllDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, BoundaryObject bound, DofContainer& dofs); DofContainer& dofs, bool baseDofPtr = false); /** \} */ /** \} */ ... ...
 ... @@ -82,7 +82,7 @@ namespace AMDiS { ... @@ -82,7 +82,7 @@ namespace AMDiS { /// Returns a pointer to the starting position of the current DOF /// Returns a pointer to the starting position of the current DOF /// array. Makes only sence, if \ref nextStrict() is used for traverse. /// array. Makes only sence, if \ref nextStrict() is used for traverse. inline const DegreeOfFreedom* getStartDof() inline const DegreeOfFreedom* getBaseDof() { { return dofs[node0 + elementPos]; return dofs[node0 + elementPos]; } } ... ...
 ... @@ -137,4 +137,18 @@ namespace AMDiS { ... @@ -137,4 +137,18 @@ namespace AMDiS { } } } } } } const FiniteElemSpace* FiniteElemSpace::getHighest(vector& feSpaces) { const FiniteElemSpace *feSpace = feSpaces[0]; for (unsigned int i = 1; i < feSpaces.size(); i++) if (feSpaces[i]->getBasisFcts()->getDegree() > feSpace->getBasisFcts()->getDegree()) feSpace = feSpaces[i]; return feSpace; } } }
 ... @@ -97,6 +97,11 @@ namespace AMDiS { ... @@ -97,6 +97,11 @@ namespace AMDiS { int calcMemoryUsage(); int calcMemoryUsage(); static void clear(); static void clear(); /// Returns for a set of FE spaces that FE space having basis functions with /// the highest degree. static const FiniteElemSpace* getHighest(vector& feSpaces); protected: protected: /** \brief /** \brief ... ...
 ... @@ -173,23 +173,23 @@ namespace AMDiS { ... @@ -173,23 +173,23 @@ namespace AMDiS { return "Line"; return "Line"; } } void getNodeDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const void getNodeDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&, bool) const { { FUNCNAME("Line::getNodeDofs()"); FUNCNAME("Line::getNodeDofs()"); ERROR_EXIT("Not yet implemented!\n"); ERROR_EXIT("Not yet implemented!\n"); } } void getHigherOrderDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const void getHigherOrderDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&, bool) const { { FUNCNAME("Line::getHigherOrderDofs()"); FUNCNAME("Line::getHigherOrderDofs()"); ERROR_EXIT("Not yet implemented!\n"); ERROR_EXIT("Not yet implemented!\n"); } } protected: protected: /** \brief /// vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th * vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th /// edge of this element. * edge of this element. */ static const int vertexOfEdge[1][2]; static const int vertexOfEdge[1][2]; static const int sideOfChild[2][2]; static const int sideOfChild[2][2]; ... ...
 ... @@ -279,16 +279,10 @@ namespace AMDiS { ... @@ -279,16 +279,10 @@ namespace AMDiS { TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n"); TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n"); // Search for the FE space with the highest degree of polynomials. Using this // === Search for the FE space with the highest degree of polynomials. === // FE space ensures that deleting DOFs defined on it, also DOFs of lower // === Using this FE space ensures that deleting DOFs defined on it, === // order FE spaces will be deleted correctly. // === also DOFs of lower order FE spaces will be deleted correctly. === const FiniteElemSpace *feSpace = FiniteElemSpace::getHighest(feSpaces); const FiniteElemSpace *feSpace = feSpaces[0]; for (unsigned int i = 1; i < feSpaces.size(); i++) if (feSpaces[i]->getBasisFcts()->getDegree() > feSpace->getBasisFcts()->getDegree()) feSpace = feSpaces[i]; // === Determine to all DOFs in mesh the macro elements where the DOF === // === Determine to all DOFs in mesh the macro elements where the DOF === ... @@ -305,8 +299,8 @@ namespace AMDiS { ... @@ -305,8 +299,8 @@ namespace AMDiS { while (elInfo) { while (elInfo) { elDofIter.reset(elInfo->getElement()); elDofIter.reset(elInfo->getElement()); do { do { dofsOwner[elDofIter.getStartDof()].insert(elInfo->getMacroElement()); dofsOwner[elDofIter.getBaseDof()].insert(elInfo->getMacroElement()); dofsPosIndex[elDofIter.getStartDof()] = elDofIter.getPosIndex(); dofsPosIndex[elDofIter.getBaseDof()] = elDofIter.getPosIndex(); } while (elDofIter.nextStrict()); } while (elDofIter.nextStrict()); elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo); ... @@ -322,9 +316,9 @@ namespace AMDiS { ... @@ -322,9 +316,9 @@ namespace AMDiS { deque newMacroElements; deque newMacroElements; for (deque::iterator elIter = macroElements.begin(); for (deque::iterator elIter = macroElements.begin(); elIter != macroElements.end(); ++elIter) { elIter != macroElements.end(); ++elIter) { // If the current mesh macro element should not be deleted, i.e., it is not a // If the current mesh macro element should not be deleted, i.e., it is not // member of the list of macro elements to be deleted, is is inserted to the new // a member of the list of macro elements to be deleted, is is inserted to // macro element list. // the new macro element list. if (macros.find(*elIter) == macros.end()) if (macros.find(*elIter) == macros.end()) newMacroElements.push_back(*elIter); newMacroElements.push_back(*elIter); } } ... @@ -334,15 +328,15 @@ namespace AMDiS { ... @@ -334,15 +328,15 @@ namespace AMDiS { macroElements = newMacroElements; macroElements = newMacroElements; // === For all macro elements to be deleted, delete them also to be neighbours === // === For all macro elements to be deleted, delete them also to be === // === of some other macro elements. Furtheremore, delete the whole element === // === neighbours of some other macro elements. Furtheremore, delete the === // === hierarchie structure of the macro element. === // === whole element hierarchie structure of the macro element. === for (std::set::iterator macroIt = macros.begin(); for (std::set::iterator macroIt = macros.begin(); macroIt != macros.end(); ++macroIt) { macroIt != macros.end(); ++macroIt) { // Go through all neighbours of the macro element and remove this macro element // Go through all neighbours of the macro element and remove this macro // to be neighbour of some other macro element. // element to be neighbour of some other macro element. for (int i = 0; i < getGeo(NEIGH); i++) for (int i = 0; i < getGeo(NEIGH); i++) if ((*macroIt)->getNeighbour(i)) if ((*macroIt)->getNeighbour(i)) for (int j = 0; j < getGeo(NEIGH); j++) for (int j = 0; j < getGeo(NEIGH); j++) ... @@ -361,7 +355,7 @@ namespace AMDiS { ... @@ -361,7 +355,7 @@ namespace AMDiS { } } // We will delete at least some element DOFs in the next but will keep the // We will delete at least some element DOFs in the next but will keep the // element object still in memory. Therefore the DOF pointer must be set to be // element object still in memory. Therefore the DOF pointer must be set // invalid. // invalid. (*macroIt)->getElement()->setDofValid(false); (*macroIt)->getElement()->setDofValid(false); } } ... ...
 ... @@ -10,15 +10,17 @@ ... @@ -10,15 +10,17 @@ // See also license.opensource.txt in the distribution. // See also license.opensource.txt in the distribution. #include "Debug.h" #include "DOFVector.h" #include "Element.h" #include "ElementDofIterator.h" #include "ElInfo.h" #include "MeshStructure.h" #include "MeshStructure.h" #include "MeshStructure_ED.h" #include "MeshStructure_ED.h" #include "Mesh.h" #include "Mesh.h" #include "Element.h" #include "Traverse.h" #include "ElInfo.h" #include "RefinementManager.h" #include "RefinementManager.h" #include "Debug.h" #include "Traverse.h" #include "DOFVector.h" namespace AMDiS { namespace AMDiS { ... @@ -416,66 +418,104 @@ namespace AMDiS { ... @@ -416,66 +418,104 @@ namespace AMDiS { } } void MeshStructure::getMeshStructureValues(Mesh *mesh, void MeshStructure::getMeshStructureValues(int macroElIndex, int macroElIndex, const DOFVector* vec, const DOFVector* vec, std::vector& values) std::vector& values) { { FUNCNAME("MeshStructure::getMeshStructureValues()"); FUNCNAME("MeshStructure::getMeshStructureValues()"); TEST_EXIT_DBG(mesh)("No mesh defined!\n"); TEST_EXIT_DBG(vec)("No DOFVector defined!\n"); TEST_EXIT_DBG(vec)("No DOFVector defined!\n"); const FiniteElemSpace *feSpace = vec->getFeSpace(); Mesh *mesh = feSpace->getMesh(); bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1); values.clear(); values.clear(); ElementDofIterator elDofIter(feSpace); TraverseStack stack; TraverseStack stack; ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_EVERY_EL_PREORDER); Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { while (elInfo) { // For the macro element the mesh structure code stores all vertex values. if (elInfo->getLevel() == 0) if (elInfo->getLevel() == 0) for (int i = 0; i < mesh->getGeo(VERTEX); i++) for (int i = 0; i < mesh->getGeo(VERTEX); i++) values.push_back((*vec)[elInfo->getElement()->getDof(i, 0)]);