From a6bacc8a739c92b367549957fe3c474e2f950002 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski Date: Thu, 9 Feb 2012 14:08:33 +0000 Subject: [PATCH] Work on repartitioning with higher order finite elements. --- AMDiS/src/Containers.h | 50 +++++++++++ AMDiS/src/Element.cc | 7 +- AMDiS/src/Element.h | 114 ++++++++++++------------- AMDiS/src/ElementDofIterator.h | 2 +- AMDiS/src/FiniteElemSpace.cc | 14 +++ AMDiS/src/FiniteElemSpace.h | 5 ++ AMDiS/src/Line.h | 12 +-- AMDiS/src/Mesh.cc | 36 ++++---- AMDiS/src/MeshStructure.cc | 72 ++++++++++++---- AMDiS/src/MeshStructure.h | 26 +++++- AMDiS/src/Tetrahedron.cc | 104 +++++++++++++--------- AMDiS/src/Tetrahedron.h | 47 +++++----- AMDiS/src/Triangle.cc | 29 +++++-- AMDiS/src/Triangle.h | 12 +-- AMDiS/src/parallel/ElementObjectData.h | 1 + AMDiS/src/parallel/MeshDistributor.cc | 37 +++----- AMDiS/src/parallel/MeshDistributor.h | 1 + AMDiS/src/parallel/MeshManipulation.cc | 29 +++++-- AMDiS/src/parallel/MeshManipulation.h | 2 +- 19 files changed, 378 insertions(+), 222 deletions(-) create mode 100644 AMDiS/src/Containers.h diff --git a/AMDiS/src/Containers.h b/AMDiS/src/Containers.h new file mode 100644 index 00000000..88782f99 --- /dev/null +++ b/AMDiS/src/Containers.h @@ -0,0 +1,50 @@ +// ============================================================================ +// == == +// == 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 diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc index 21e1e92c..f0c6fb8a 100644 --- a/AMDiS/src/Element.cc +++ b/AMDiS/src/Element.cc @@ -627,12 +627,13 @@ namespace AMDiS { void Element::getAllDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) + DofContainer& dofs, + bool baseDofPtr) { - getNodeDofs(feSpace, bound, dofs); + getNodeDofs(feSpace, bound, dofs, baseDofPtr); if (feSpace->getBasisFcts()->getDegree() > 1) - getHigherOrderDofs(feSpace, bound, dofs); + getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr); } } diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index b69e1606..99aeb660 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -82,10 +82,8 @@ namespace AMDiS { /// void deleteElementDOFs(); - /** \brief - * Clone this Element and return a reference to it. Because also the DOFs - * are cloned, \ref Mesh::serializedDOfs must be used. - */ + /// Clone this Element and return a reference to it. Because also the DOFs + /// are cloned, \ref Mesh::serializedDOfs must be used. Element* cloneWithDOFs(); /** \name getting methods @@ -111,10 +109,8 @@ namespace AMDiS { return child[i]; } - /** \brief - * Returns true if Element is a leaf element (\ref child[0] == NULL), returns - * false otherwise. - */ + /// Returns true if Element is a leaf element (\ref child[0] == NULL), returns + /// false otherwise. inline const bool isLeaf() const { return (child[0] == NULL); @@ -155,10 +151,8 @@ namespace AMDiS { dofValid = b; } - /** \brief - * Returns \ref elementData's error estimation, if Element is a leaf element - * and has leaf data. - */ + /// Returns \ref elementData's error estimation, if Element is a leaf element + /// and has leaf data. inline double getEstimation(int row) const { if (isLeaf()) { @@ -172,10 +166,8 @@ namespace AMDiS { return 0.0; } - /** \brief - * 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. - */ + /// 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. inline double getCoarseningEstimation(int row) { if (isLeaf()) { @@ -195,10 +187,8 @@ namespace AMDiS { /// Returns local vertex number of the j-th vertex of the i-th edge virtual int getVertexOfEdge(int i, int j) const = 0; - /** \brief - * Returns local vertex number of the vertexIndex-th vertex of the - * positionIndex-th part of type position (vertex, edge, face) - */ + /// Returns local vertex number of the vertexIndex-th vertex of the + /// positionIndex-th part of type position (vertex, edge, face) virtual int getVertexOfPosition(GeoIndex position, int positionIndex, int vertexIndex) const = 0; @@ -215,10 +205,6 @@ namespace AMDiS { /// 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 virtual int getGeo(GeoIndex i) const = 0; @@ -271,10 +257,8 @@ namespace AMDiS { elementData = ed; } - /** \brief - * Sets \ref newCoord of Element. Needed by refinement, if Element has a - * boundary edge on a curved boundary. - */ + /// Sets \ref newCoord of Element. Needed by refinement, if Element has a + /// boundary edge on a curved boundary. inline void setNewCoord(WorldVector* coord) { newCoord = coord; @@ -292,10 +276,8 @@ namespace AMDiS { dof[pos] = p; } - /** \brief - * 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. - */ + /// 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. inline void setEstimation(double est, int row) { FUNCNAME("Element::setEstimation()"); @@ -312,10 +294,8 @@ namespace AMDiS { } } - /** \brief - * 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. - */ + /// 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. inline void setCoarseningEstimation(double est, int row) { if (isLeaf()) { @@ -363,18 +343,19 @@ namespace AMDiS { virtual const FixVec& 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; - /** \brief - * Returns which side of child[childnr] corresponds to side sidenr of this - * Element. If the child has no corresponding side, the return value is negative. - */ + /// Returns which side of child[childnr] corresponds to side sidenr of + /// this Element. If the child has no corresponding side, the return value + /// is negative. virtual int getSideOfChild(int childnr, int sidenr, int elType = 0) const = 0; /** \brief - * Generalization of \ref getSideOfChild to arbitrary subObject. Thus, e.g., in 3d - * we can ask for the local id of a verte, edge or face on the elements children. + * Generalization of \ref getSideOfChild to arbitrary subObject. Thus, + * 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] subObj Defines whether we ask for VERTEX, EDGE or FACE. @@ -384,12 +365,12 @@ namespace AMDiS { virtual int getSubObjOfChild(int childnr, GeoIndex subObj, int ithObj, int elType = 0) const = 0; - /** \brief - * Returns which vertex of elements parent corresponds to the vertexnr of - * 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. - */ - virtual int getVertexOfParent(int childnr, int vertexnr, int elType = 0) const = 0; + /// Returns which vertex of elements parent corresponds to the vertexnr of + /// 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. + virtual int getVertexOfParent(int childnr, + int vertexnr, + int elType = 0) const = 0; /// Returns whether Element is a Line virtual bool isLine() const = 0; @@ -416,14 +397,19 @@ namespace AMDiS { * nodes alonge this vertex/edge/face are assembled and put together to * a list. * - * \param[in] feSpace FE space which is used to get the dofs. - * \param[in] bound Defines the vertex/edge/face of the element on - * which all vertex dofs are assembled. - * \param[out] dofs List of dofs, where the result is stored. + * \param[in] feSpace FE space which is used to get the dofs. + * \param[in] bound Defines the vertex/edge/face of the element on + * which all vertex dofs are assembled. + * \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, BoundaryObject bound, - DofContainer& dofs) const = 0; + DofContainer& dofs, + bool baseDofPtr = false) const = 0; /** \brief * Traverses a vertex/edge/face of a given element (this includes also all @@ -431,18 +417,26 @@ namespace AMDiS { * to higher order basis functions alonge this vertex/edge/face are * assembled and put together to a list. * - * \param[in] feSpace FE space which is used to get the dofs. - * \param[in] bound Defines the edge/face of the element on which - * all non vertex dofs are assembled. - * \param[out] dofs All dofs are put to this dof list. + * \param[in] feSpace FE space which is used to get the dofs. + * \param[in] bound Defines the edge/face of the element on which + * all non vertex dofs are assembled. + * \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, 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, BoundaryObject bound, - DofContainer& dofs); + DofContainer& dofs, + bool baseDofPtr = false); /** \} */ diff --git a/AMDiS/src/ElementDofIterator.h b/AMDiS/src/ElementDofIterator.h index 6001b5fa..8b02219b 100644 --- a/AMDiS/src/ElementDofIterator.h +++ b/AMDiS/src/ElementDofIterator.h @@ -82,7 +82,7 @@ namespace AMDiS { /// Returns a pointer to the starting position of the current DOF /// array. Makes only sence, if \ref nextStrict() is used for traverse. - inline const DegreeOfFreedom* getStartDof() + inline const DegreeOfFreedom* getBaseDof() { return dofs[node0 + elementPos]; } diff --git a/AMDiS/src/FiniteElemSpace.cc b/AMDiS/src/FiniteElemSpace.cc index 46c76d0f..5964a88e 100644 --- a/AMDiS/src/FiniteElemSpace.cc +++ b/AMDiS/src/FiniteElemSpace.cc @@ -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; + } + } diff --git a/AMDiS/src/FiniteElemSpace.h b/AMDiS/src/FiniteElemSpace.h index 518da190..177a9fd9 100644 --- a/AMDiS/src/FiniteElemSpace.h +++ b/AMDiS/src/FiniteElemSpace.h @@ -97,6 +97,11 @@ namespace AMDiS { int calcMemoryUsage(); 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: /** \brief diff --git a/AMDiS/src/Line.h b/AMDiS/src/Line.h index 2c8de3ad..9a6a1c98 100644 --- a/AMDiS/src/Line.h +++ b/AMDiS/src/Line.h @@ -173,23 +173,23 @@ namespace AMDiS { return "Line"; } - void getNodeDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const + void getNodeDofs(const FiniteElemSpace*, BoundaryObject, + DofContainer&, bool) const { FUNCNAME("Line::getNodeDofs()"); ERROR_EXIT("Not yet implemented!\n"); } - void getHigherOrderDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const + void getHigherOrderDofs(const FiniteElemSpace*, BoundaryObject, + DofContainer&, bool) const { FUNCNAME("Line::getHigherOrderDofs()"); ERROR_EXIT("Not yet implemented!\n"); } protected: - /** \brief - * vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th - * edge of this element. - */ + /// vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th + /// edge of this element. static const int vertexOfEdge[1][2]; static const int sideOfChild[2][2]; diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index edbc576a..b43b0ca9 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -279,16 +279,10 @@ namespace AMDiS { TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n"); - - // === Search for the FE space with the highest degree of polynomials. === - // === Using this FE space ensures that deleting DOFs defined on it, === - // === also DOFs of lower order FE spaces will be deleted correctly. === - - 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]; + // Search for the FE space with the highest degree of polynomials. Using this + // FE space ensures that deleting DOFs defined on it, also DOFs of lower + // order FE spaces will be deleted correctly. + const FiniteElemSpace *feSpace = FiniteElemSpace::getHighest(feSpaces); // === Determine to all DOFs in mesh the macro elements where the DOF === @@ -305,8 +299,8 @@ namespace AMDiS { while (elInfo) { elDofIter.reset(elInfo->getElement()); do { - dofsOwner[elDofIter.getStartDof()].insert(elInfo->getMacroElement()); - dofsPosIndex[elDofIter.getStartDof()] = elDofIter.getPosIndex(); + dofsOwner[elDofIter.getBaseDof()].insert(elInfo->getMacroElement()); + dofsPosIndex[elDofIter.getBaseDof()] = elDofIter.getPosIndex(); } while (elDofIter.nextStrict()); elInfo = stack.traverseNext(elInfo); @@ -322,9 +316,9 @@ namespace AMDiS { deque newMacroElements; for (deque::iterator elIter = macroElements.begin(); elIter != macroElements.end(); ++elIter) { - // If the current mesh macro element should not be deleted, i.e., it is not a - // member of the list of macro elements to be deleted, is is inserted to the new - // macro element list. + // If the current mesh macro element should not be deleted, i.e., it is not + // a member of the list of macro elements to be deleted, is is inserted to + // the new macro element list. if (macros.find(*elIter) == macros.end()) newMacroElements.push_back(*elIter); } @@ -334,15 +328,15 @@ namespace AMDiS { macroElements = newMacroElements; - // === For all macro elements to be deleted, delete them also to be neighbours === - // === of some other macro elements. Furtheremore, delete the whole element === - // === hierarchie structure of the macro element. === + // === For all macro elements to be deleted, delete them also to be === + // === neighbours of some other macro elements. Furtheremore, delete the === + // === whole element hierarchie structure of the macro element. === for (std::set::iterator macroIt = macros.begin(); macroIt != macros.end(); ++macroIt) { - // Go through all neighbours of the macro element and remove this macro element - // to be neighbour of some other macro element. + // Go through all neighbours of the macro element and remove this macro + // element to be neighbour of some other macro element. for (int i = 0; i < getGeo(NEIGH); i++) if ((*macroIt)->getNeighbour(i)) for (int j = 0; j < getGeo(NEIGH); j++) @@ -361,7 +355,7 @@ namespace AMDiS { } // 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. (*macroIt)->getElement()->setDofValid(false); } diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc index d3650818..8e03083a 100644 --- a/AMDiS/src/MeshStructure.cc +++ b/AMDiS/src/MeshStructure.cc @@ -10,15 +10,17 @@ // 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_ED.h" #include "Mesh.h" -#include "Element.h" -#include "Traverse.h" -#include "ElInfo.h" #include "RefinementManager.h" -#include "Debug.h" -#include "DOFVector.h" +#include "Traverse.h" + namespace AMDiS { @@ -416,66 +418,104 @@ namespace AMDiS { } - void MeshStructure::getMeshStructureValues(Mesh *mesh, - int macroElIndex, + void MeshStructure::getMeshStructureValues(int macroElIndex, const DOFVector* vec, std::vector& values) { FUNCNAME("MeshStructure::getMeshStructureValues()"); - TEST_EXIT_DBG(mesh)("No mesh 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(); + ElementDofIterator elDofIter(feSpace); TraverseStack stack; ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { + // For the macro element the mesh structure code stores all vertex values. if (elInfo->getLevel() == 0) for (int i = 0; i < mesh->getGeo(VERTEX); i++) values.push_back((*vec)[elInfo->getElement()->getDof(i, 0)]); - if (!elInfo->getElement()->isLeaf()) - values.push_back((*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)]); + if (!elInfo->getElement()->isLeaf()) { + // If no leaf element store the value of the "new" DOF that is created + // by bisectioning of this element. + + DegreeOfFreedom dof0 = + elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0); + values.push_back((*vec)[dof0]); + } else { + // If leaf element store all non vertex values of this element, thus + // only relevant for higher order basis functions. + + if (feSpaceHasNonVertexDofs) { + elDofIter.reset(elInfo->getElement()); + do { + if (elDofIter.getPosIndex() != VERTEX) + values.push_back((*vec)[elDofIter.getDof()]); + } while (elDofIter.next()); + } + } elInfo = stack.traverseNext(elInfo); } } - void MeshStructure::setMeshStructureValues(Mesh *mesh, - int macroElIndex, + void MeshStructure::setMeshStructureValues(int macroElIndex, DOFVector* vec, const std::vector& values) { FUNCNAME("MeshStructure::setMeshStructureValues()"); - TEST_EXIT_DBG(mesh)("No mesh defined!\n"); TEST_EXIT_DBG(vec)("No DOFVector defined!\n"); - TEST_EXIT_DBG(static_cast(values.size()) >= mesh->getGeo(VERTEX)) - ("Should not happen!\n"); + const FiniteElemSpace *feSpace = vec->getFeSpace(); + Mesh *mesh = feSpace->getMesh(); + bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1); unsigned int counter = 0; + TEST_EXIT_DBG(static_cast(values.size()) >= mesh->getGeo(VERTEX)) + ("Should not happen!\n"); + + ElementDofIterator elDofIter(feSpace); TraverseStack stack; ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { + // For the macro element all vertex nodes are set first. if (elInfo->getLevel() == 0) for (int i = 0; i < mesh->getGeo(VERTEX); i++) (*vec)[elInfo->getElement()->getDof(i, 0)] = values[counter++]; if (!elInfo->getElement()->isLeaf()) { + // If no leaf element set the value of the "new" DOF that is created + // by bisectioning of this element. TEST_EXIT_DBG(counter < values.size())("Should not happen!\n"); (*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)] = values[counter++]; + } else { + // On leaf elements set all non vertex values (thus DOFs of higher order + // basis functions). + + if (feSpaceHasNonVertexDofs) { + elDofIter.reset(elInfo->getElement()); + do { + if (elDofIter.getPosIndex() != VERTEX) + (*vec)[elDofIter.getDof()] = values[counter++]; + } while (elDofIter.next()); + } } elInfo = stack.traverseNext(elInfo); } + TEST_EXIT_DBG(values.size() == counter)("Should not happen!\n"); } } diff --git a/AMDiS/src/MeshStructure.h b/AMDiS/src/MeshStructure.h index e9f2732d..c9b86f3e 100644 --- a/AMDiS/src/MeshStructure.h +++ b/AMDiS/src/MeshStructure.h @@ -154,14 +154,32 @@ namespace AMDiS { /// Returns true, if the given mesh structure code is equal to this one. bool compare(MeshStructure &other); - void getMeshStructureValues(Mesh *mesh, - int macroElIndex, + /** \brief + * Creates a value array of a given DOFVector. This value array corresponds + * to the mesh structure code of the element and thus can easily be used + * to reconstruct the values of the DOFVector on the same element (e.g., after + * the mesh and the value array has been redistributed in parallel + * computations). + * + * \param[in] macroElIndex Index of the macro element for which the value + * structure code must be created. + * \param[in] vec DOFVector to be used for creating the value code. + * \param[out] values Resulting value structure code. + */ + void getMeshStructureValues(int macroElIndex, const DOFVector* vec, std::vector& values); - void setMeshStructureValues(Mesh *mesh, - int macroElIndex, + /** \brief + * Uses a value structure code, e.g. created by \ref getMeshStructureValues, + * to reconstruct the data of a DOFVector on a given macro element. + * + * \param[in] macroElIndex Macro element index the code is related to. + * \param[out] vec DOFVector that should be reconstructed. + * \param[in] values Value structure code. + */ + void setMeshStructureValues(int macroElIndex, DOFVector* vec, const std::vector& values); diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index c590c795..4da2debf 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -202,22 +202,27 @@ namespace AMDiS { void Tetrahedron::getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const + DofContainer& dofs, + bool baseDofPtr) const { FUNCNAME("Tetrahedron::getNodeDofs()"); switch (bound.subObj) { case VERTEX: { - int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX); - dofs.push_back(&(dof[bound.ithObj][n0])); + if (baseDofPtr) { + dofs.push_back(dof[bound.ithObj]); + } else { + int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX); + dofs.push_back(&(dof[bound.ithObj][n0])); + } } break; case EDGE: - getNodeDofsAtEdge(feSpace, bound, dofs); + getNodeDofsAtEdge(feSpace, bound, dofs, baseDofPtr); break; case FACE: - getNodeDofsAtFace(feSpace, bound, dofs); + getNodeDofsAtFace(feSpace, bound, dofs, baseDofPtr); break; default: ERROR_EXIT("Should not happen!\n"); @@ -227,7 +232,8 @@ namespace AMDiS { void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const + DofContainer& dofs, + bool baseDofPtr) const { FUNCNAME("Tetrahedron::getNodeDofsAtFace()"); @@ -239,21 +245,22 @@ namespace AMDiS { { BoundaryObject nextBound = bound; prepareNextBound(nextBound, 1); - child[1]->getNodeDofs(feSpace, nextBound, dofs); + child[1]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr); } break; case 1: { BoundaryObject nextBound = bound; prepareNextBound(nextBound, 0); - child[0]->getNodeDofs(feSpace, nextBound, dofs); + child[0]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr); } break; case 2: case 3: { - int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX); + int n0 = + (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(FACE)); BoundaryObject nextBound0 = bound, nextBound1 = bound; prepareNextBound(nextBound0, 0); prepareNextBound(nextBound1, 1); @@ -267,19 +274,19 @@ namespace AMDiS { addDof = false; if (bound.reverseMode) { - child[1]->getNodeDofs(feSpace, nextBound1, dofs); + child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr); if (addDof) dofs.push_back(&(child[0]->getDof()[3][n0])); - child[0]->getNodeDofs(feSpace, nextBound0, dofs); + child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr); } else { - child[0]->getNodeDofs(feSpace, nextBound0, dofs); + child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr); if (addDof) dofs.push_back(&(child[0]->getDof()[3][n0])); - child[1]->getNodeDofs(feSpace, nextBound1, dofs); + child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr); } } break; @@ -291,14 +298,15 @@ namespace AMDiS { void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const + DofContainer& dofs, + bool baseDofPtr) const { FUNCNAME("Tetrahedron::getNodeDofsAtEdge()"); if (!child[0]) return; - int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX); + int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(EDGE)); BoundaryObject nextBound0 = bound, nextBound1 = bound; prepareNextBound(nextBound0, 0); prepareNextBound(nextBound1, 1); @@ -309,13 +317,13 @@ namespace AMDiS { nextBound1.reverseMode = false; if (bound.reverseMode) { - child[1]->getNodeDofs(feSpace, nextBound0, dofs); + child[1]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr); dofs.push_back(&(child[0]->getDof()[3][n0])); - child[0]->getNodeDofs(feSpace, nextBound1, dofs); + child[0]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr); } else { - child[0]->getNodeDofs(feSpace, nextBound0, dofs); + child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr); dofs.push_back(&(child[0]->getDof()[3][n0])); - child[1]->getNodeDofs(feSpace, nextBound1, dofs); + child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr); } break; @@ -324,17 +332,17 @@ namespace AMDiS { ("Should not happen!\n"); if (nextBound0.ithObj != -1) - child[0]->getNodeDofs(feSpace, nextBound0, dofs); + child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr); break; default: TEST_EXIT_DBG(nextBound0.ithObj == -1 || nextBound1.ithObj == -1) ("This should not happen!\n"); if (nextBound0.ithObj != -1) - child[0]->getNodeDofs(feSpace, nextBound0, dofs); + child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr); if (nextBound1.ithObj != -1) - child[1]->getNodeDofs(feSpace, nextBound1, dofs); + child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr); } } @@ -367,7 +375,8 @@ namespace AMDiS { void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const + DofContainer& dofs, + bool baseDofPtr) const { FUNCNAME("Tetrahedron::getHigherOrderDofs()"); @@ -394,14 +403,14 @@ namespace AMDiS { if (bound.reverseMode) { if (nextBound1.ithObj >= 0) - child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs); + child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, baseDofPtr); if (nextBound0.ithObj >= 0) - child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs); + child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs, baseDofPtr); } else { if (nextBound0.ithObj >= 0) - child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs); + child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs, baseDofPtr); if (nextBound1.ithObj >= 0) - child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs); + child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, baseDofPtr); } } else { // Either the edge is not contained in further refined children, or @@ -410,11 +419,20 @@ namespace AMDiS { ElementDofIterator elDofIter(feSpace, true); elDofIter.reset(this); - do { - if (elDofIter.getCurrentPos() == 1 && - elDofIter.getCurrentElementPos() == bound.ithObj) - dofs.push_back(elDofIter.getDofPtr()); - } while(elDofIter.next()); + + if (baseDofPtr) { + do { + if (elDofIter.getCurrentPos() == 1 && + elDofIter.getCurrentElementPos() == bound.ithObj) + dofs.push_back(elDofIter.getBaseDof()); + } while (elDofIter.nextStrict()); + } else { + do { + if (elDofIter.getCurrentPos() == 1 && + elDofIter.getCurrentElementPos() == bound.ithObj) + dofs.push_back(elDofIter.getDofPtr()); + } while (elDofIter.next()); + } } } @@ -433,21 +451,29 @@ namespace AMDiS { if (childFace0 != -1) { nextBound.ithObj = childFace0; - child[0]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[0]->getHigherOrderDofs(feSpace, nextBound, dofs, baseDofPtr); } if (childFace1 != -1) { nextBound.ithObj = childFace1; - child[1]->getHigherOrderDofs(feSpace, nextBound, dofs); + child[1]->getHigherOrderDofs(feSpace, nextBound, dofs, baseDofPtr); } } else { ElementDofIterator elDofIter(feSpace, true); elDofIter.reset(this); - do { - if (elDofIter.getCurrentPos() == 2 && - elDofIter.getCurrentElementPos() == bound.ithObj) - dofs.push_back(elDofIter.getDofPtr()); - } while(elDofIter.next()); + if (baseDofPtr) { + do { + if (elDofIter.getCurrentPos() == 2 && + elDofIter.getCurrentElementPos() == bound.ithObj) + dofs.push_back(elDofIter.getBaseDof()); + } while (elDofIter.nextStrict()); + } else { + do { + if (elDofIter.getCurrentPos() == 2 && + elDofIter.getCurrentElementPos() == bound.ithObj) + dofs.push_back(elDofIter.getDofPtr()); + } while (elDofIter.next()); + } } } break; diff --git a/AMDiS/src/Tetrahedron.h b/AMDiS/src/Tetrahedron.h index b801be84..089eb842 100644 --- a/AMDiS/src/Tetrahedron.h +++ b/AMDiS/src/Tetrahedron.h @@ -138,28 +138,29 @@ namespace AMDiS { void getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const; + DofContainer& dofs, + bool baseDofPtr = false) const; void getNodeDofsAtFace(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const; + DofContainer& dofs, + bool baseDofPtr) const; void getNodeDofsAtEdge(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const; + DofContainer& dofs, + bool baseDofPtr) const; void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const; + DofContainer& dofs, + bool baseDofPtr = false) const; void prepareNextBound(BoundaryObject &bound, int ithChild) const; public: - /** \brief - * nChildEdge[el_type][ichild][dir] - * gives local index of new edge on child[ichild] part of face [2+dir] on - * the parent - */ + /// nChildEdge[el_type][ichild][dir] gives local index of new edge on + /// child[ichild] part of face [2 + dir] on the parent. static const unsigned char nChildEdge[3][2][2]; /** \brief @@ -216,27 +217,28 @@ namespace AMDiS { TEST_EXIT_DBG(child == 0 || child == 1)("Child must be in (0,1)!\n"); TEST_EXIT_DBG(side >= 0 && side <= 3)("Side must be between 0 and 3!\n"); - TEST_EXIT_DBG(elType >= 0 && elType <= 2)("ElType must be between 0 and 2!\n"); + TEST_EXIT_DBG(elType >= 0 && elType <= 2) + ("ElType must be between 0 and 2!\n"); return sideOfChild[elType][child][side]; } - /** \brief - * Returns for an edge number its local edge number on a child element. See - * \ref edgeOfChild for mor information. - */ + /// Returns for an edge number its local edge number on a child element. See + /// \ref edgeOfChild for mor information. inline int getEdgeOfChild(int child, int edge, int elType) const { FUNCNAME("Tetrahedron::getEdgeOfChild()"); TEST_EXIT_DBG(child == 0 || child == 1)("Child must be in (0,1)!\n"); TEST_EXIT_DBG(edge >= 0 && edge <= 5)("Side must be between 0 and 3!\n"); - TEST_EXIT_DBG(elType >= 0 && elType <= 2)("ElType must be between 0 and 2!\n"); + TEST_EXIT_DBG(elType >= 0 && elType <= 2) + ("ElType must be between 0 and 2!\n"); return edgeOfChild[elType][child][edge]; } - int getSubObjOfChild(int childnr, GeoIndex subObj, int ithObj, int elType) const + int getSubObjOfChild(int childnr, GeoIndex subObj, + int ithObj, int elType) const { FUNCNAME("Tetrahedron::getSubObjOfChild()"); @@ -255,7 +257,8 @@ namespace AMDiS { TEST_EXIT_DBG(child == 0 || child == 1)("Child must be in (0,1)!\n"); TEST_EXIT_DBG(side >= 0 && side <= 3)("Side must be between 0 and 3!\n"); - TEST_EXIT_DBG(elType >= 0 && elType <= 2)("ElType must be between 0 and 2!\n"); + TEST_EXIT_DBG(elType >= 0 && elType <= 2) + ("ElType must be between 0 and 2!\n"); return vertexOfParent[elType][child][side]; } @@ -323,12 +326,10 @@ namespace AMDiS { /// See \ref Element::getSideOfChild for more information. static const int sideOfChild[3][2][4]; - /** \brief - * edgeOfChild[elType][i][j] is the local edge number of the j-th edge within - * the i-th children of an element of elType. If the value is -1, the edge is - * not included in the element's child. Note that the 0 edge is included in - * both children only by its half. - */ + /// edgeOfChild[elType][i][j] is the local edge number of the j-th edge within + /// the i-th children of an element of elType. If the value is -1, the edge is + /// not included in the element's child. Note that the 0 edge is included in + /// both children only by its half. static const int edgeOfChild[3][2][6]; /// See \ref Element::getVertexOfParent for more information. diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc index 8fd3b780..74c3a9b9 100644 --- a/AMDiS/src/Triangle.cc +++ b/AMDiS/src/Triangle.cc @@ -89,15 +89,17 @@ namespace AMDiS { void Triangle::getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const + DofContainer& dofs, + bool baseDofPtr) const { FUNCNAME("Triangle::getNodeDofs()"); // Get displacement if more than one FE space is defined on mesh. - int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX); + int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX)); if (bound.subObj == VERTEX) { dofs.push_back(&(dof[bound.ithObj][n0])); + return; } @@ -114,7 +116,7 @@ namespace AMDiS { if (bound.reverseMode) { nextBound.ithObj = 1; child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs); - dofs.push_back(&(elDofs[2][n0])); + dofs.push_back(&(elDofs[2][n0])); nextBound.ithObj = 0; child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs); } else { @@ -175,7 +177,8 @@ namespace AMDiS { void Triangle::getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const + DofContainer& dofs, + bool baseDofPtr) const { FUNCNAME("Triange::getHigherOrderDofs()"); @@ -232,12 +235,20 @@ namespace AMDiS { DofContainer addDofs; ElementDofIterator elDofIter(feSpace, true); elDofIter.reset(this); - do { - if (elDofIter.getCurrentPos() == 1 && - elDofIter.getCurrentElementPos() == bound.ithObj) - addDofs.push_back(elDofIter.getDofPtr()); - } while(elDofIter.next()); + if (baseDofPtr) { + do { + if (elDofIter.getCurrentPos() == 1 && + elDofIter.getCurrentElementPos() == bound.ithObj) + addDofs.push_back(elDofIter.getBaseDof()); + } while (elDofIter.nextStrict()); + } else { + do { + if (elDofIter.getCurrentPos() == 1 && + elDofIter.getCurrentElementPos() == bound.ithObj) + addDofs.push_back(elDofIter.getDofPtr()); + } while (elDofIter.next()); + } if (bound.reverseMode) for (int i = addDofs.size() - 1; i >= 0; i--) diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h index fa1658e1..85da4d62 100644 --- a/AMDiS/src/Triangle.h +++ b/AMDiS/src/Triangle.h @@ -196,17 +196,17 @@ namespace AMDiS { void getNodeDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const; + DofContainer& dofs, + bool baseDofPtr = false) const; void getHigherOrderDofs(const FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs) const; + DofContainer& dofs, + bool baseDofPtr = false) const; protected: - /** \brief - * vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th - * edge of this element. - */ + /// vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th + /// edge of this element. static const int vertexOfEdge[3][2]; static const int sideOfChild[2][3]; diff --git a/AMDiS/src/parallel/ElementObjectData.h b/AMDiS/src/parallel/ElementObjectData.h index dd0737f2..25d85e28 100644 --- a/AMDiS/src/parallel/ElementObjectData.h +++ b/AMDiS/src/parallel/ElementObjectData.h @@ -28,6 +28,7 @@ #include #include +#include "Containers.h" #include "Global.h" #include "Boundary.h" #include "Serializer.h" diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index e38fe776..514db842 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -173,8 +173,8 @@ namespace AMDiS { for (vector::iterator it = allMacroElements.begin(); it != allMacroElements.end(); ++it) { - macroElIndexMap[(*it)->getIndex()] = (*it)->getElement(); - macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType(); + macroElIndexMap.insert(make_pair((*it)->getIndex(), (*it)->getElement())); + macroElIndexTypeMap.insert(make_pair((*it)->getIndex(), (*it)->getElType())); } createBoundaryDofs(); @@ -621,12 +621,6 @@ namespace AMDiS { for (std::set::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) { for (std::set::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) { - - TEST_EXIT_DBG(macroElIndexMap.count(*elIt0))("Should not happen!\n"); - TEST_EXIT_DBG(macroElIndexMap.count(*elIt1))("Should not happen!\n"); - TEST_EXIT_DBG(edgeNoInEl.count(*elIt0))("Should not happen!\n"); - TEST_EXIT_DBG(edgeNoInEl.count(*elIt1))("Should not happen!\n"); - pair edge0 = make_pair(macroElIndexMap[*elIt0], edgeNoInEl[*elIt0]); pair edge1 = @@ -1239,7 +1233,7 @@ namespace AMDiS { for (unsigned int i = 0; i < interchangeVectors.size(); i++) { vector valVec; - elCode.getMeshStructureValues(mesh, *elIt, interchangeVectors[i], valVec); + elCode.getMeshStructureValues(*elIt, interchangeVectors[i], valVec); sendValues[it->first].push_back(valVec); } } @@ -1283,8 +1277,7 @@ namespace AMDiS { recvCodes[i].fitMeshToStructure(mesh, refineManager, false, *elIt); for (unsigned int k = 0; k < interchangeVectors.size(); k++) - recvCodes[i].setMeshStructureValues(mesh, - *elIt, + recvCodes[i].setMeshStructureValues(*elIt, interchangeVectors[k], recvValues[j++]); @@ -1312,7 +1305,8 @@ namespace AMDiS { partitioner->getElementInRank()[neighIndex] == false) { (*it)->setNeighbour(i, NULL); } else { - TEST_EXIT_DBG(elIndexMap.count(neighIndex) == 1)("Should not happen!\n"); + TEST_EXIT_DBG(elIndexMap.count(neighIndex) == 1) + ("Should not happen!\n"); (*it)->setNeighbour(i, elIndexMap[neighIndex]); } @@ -1341,11 +1335,8 @@ namespace AMDiS { mesh->removeMacroElements(deleteMacroElements, feSpaces); // === Remove double DOFs. === - - TEST_EXIT(feSpaces.size() == 1) - ("Repartitioning not yet implemented for multiple FE spaces!\n"); MeshManipulation meshManipulation(mesh); - meshManipulation.deleteDoubleDofs(feSpaces[0], newMacroEl, elObjects); + meshManipulation.deleteDoubleDofs(feSpaces, newMacroEl, elObjects); mesh->dofCompress(); partitioner->createPartitionMap(partitionMap); @@ -1432,15 +1423,16 @@ namespace AMDiS { // === Fills macro element data structures. === TraverseStack stack; - ElInfo *elInfo = - stack.traverseFirst(mesh, -1, - Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND); + ElInfo *elInfo = stack.traverseFirst(mesh, -1, + Mesh::CALL_LEAF_EL | + Mesh::FILL_NEIGH | + Mesh::FILL_BOUND); while (elInfo) { TEST_EXIT_DBG(elInfo->getLevel() == 0)("Should not happen!\n"); Element *el = elInfo->getElement(); - macroElIndexMap[el->getIndex()] = el; - macroElIndexTypeMap[el->getIndex()] = elInfo->getType(); + macroElIndexMap.insert(make_pair(el->getIndex(), el)); + macroElIndexTypeMap.insert(make_pair(el->getIndex(), elInfo->getType())); // Add all sub object of the element to the variable elObjects. elObjects.addElement(elInfo); @@ -1485,9 +1477,6 @@ namespace AMDiS { int owner = elObjects.getIterateOwner(); ElementObjectData& rankBoundEl = objData[mpiRank]; - TEST_EXIT_DBG(macroElIndexMap[rankBoundEl.elIndex]) - ("Should not happen!\n"); - AtomicBoundary bound; bound.rankObj.el = macroElIndexMap[rankBoundEl.elIndex]; bound.rankObj.elIndex = rankBoundEl.elIndex; diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index 0c209bb2..a7c91421 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -33,6 +33,7 @@ #include "parallel/PeriodicMap.h" #include "parallel/StdMpi.h" #include "AMDiS_fwd.h" +#include "Containers.h" #include "Global.h" #include "ProblemTimeInterface.h" #include "ProblemIterationInterface.h" diff --git a/AMDiS/src/parallel/MeshManipulation.cc b/AMDiS/src/parallel/MeshManipulation.cc index d992c72a..d678b1f2 100644 --- a/AMDiS/src/parallel/MeshManipulation.cc +++ b/AMDiS/src/parallel/MeshManipulation.cc @@ -41,12 +41,17 @@ namespace AMDiS { } - void MeshManipulation::deleteDoubleDofs(const FiniteElemSpace *feSpace, + void MeshManipulation::deleteDoubleDofs(vector& feSpaces, std::set& newMacroEl, ElementObjects &objects) { FUNCNAME("MeshManipulation::deleteDoubleDofs()"); + // Search for the FE space with the highest degree of polynomials. Using this + // FE space ensures that deleting DOFs defined on it, also DOFs of lower + // order FE spaces will be deleted correctly. + const FiniteElemSpace *feSpace = FiniteElemSpace::getHighest(feSpaces); + // Create data structure that maps macro element indices to the // corresponding pointers. map macroIndexMap; @@ -54,14 +59,11 @@ namespace AMDiS { it != newMacroEl.end(); ++it) macroIndexMap[(*it)->getIndex()] = *it; - std::set macrosProcessed; - map mapDelDofs; - - // === Traverse mesh and put all "old" macro element to macrosProcessed === // === that stores all macro elements which are really connected to the === // === overall mesh structure. === + std::set macrosProcessed; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL); while (elInfo) { @@ -80,6 +82,9 @@ namespace AMDiS { // === structure by deleting all double DOFs on macro element's vertices, === // === edges and faces. === + map mapDelDofs; + map dofPosIndex; + for (std::set::iterator it = newMacroEl.begin(); it != newMacroEl.end(); ++it) { @@ -105,7 +110,8 @@ namespace AMDiS { const DegreeOfFreedom *dof1 = el1->getDof(elIt->ithObject); mapDelDofs[dof0] = dof1; - + dofPosIndex[dof0] = VERTEX; + break; } } @@ -123,7 +129,8 @@ namespace AMDiS { continue; if (macrosProcessed.count(elIt->elIndex) == 1) { - TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))("Should not happen!\n"); + TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex)) + ("Should not happen!\n"); Element *el0 = (*it)->getElement(); Element *el1 = macroIndexMap[elIt->elIndex]->getElement(); @@ -142,8 +149,10 @@ namespace AMDiS { debug::testDofsByCoords(feSpace, dofs0, dofs1); #endif - for (unsigned int i = 0; i < dofs0.size(); i++) + for (unsigned int i = 0; i < dofs0.size(); i++) { mapDelDofs[dofs0[i]] = dofs1[i]; + dofPosIndex[dofs0[i]] = EDGE; + } break; } @@ -182,8 +191,10 @@ namespace AMDiS { debug::testDofsByCoords(feSpace, dofs0, dofs1); #endif - for (unsigned int i = 0; i < dofs0.size(); i++) + for (unsigned int i = 0; i < dofs0.size(); i++) { mapDelDofs[dofs0[i]] = dofs1[i]; + dofPosIndex[dofs0[i]] = FACE; + } break; } diff --git a/AMDiS/src/parallel/MeshManipulation.h b/AMDiS/src/parallel/MeshManipulation.h index cce08a18..c64bd8ae 100644 --- a/AMDiS/src/parallel/MeshManipulation.h +++ b/AMDiS/src/parallel/MeshManipulation.h @@ -37,7 +37,7 @@ namespace AMDiS { ~MeshManipulation(); - void deleteDoubleDofs(const FiniteElemSpace *feSpace, + void deleteDoubleDofs(vector& feSpaces, std::set& newMacroEl, ElementObjects &elObj); -- GitLab