diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am index db4664cc18c0977ecd2ac482871e5ade376e8bd3..332e434310090487128aecbe3587eefa8492f930 100644 --- a/AMDiS/bin/Makefile.am +++ b/AMDiS/bin/Makefile.am @@ -58,7 +58,7 @@ endif libamdis_la_SOURCES = \ $(PARALLEL_AMDIS_SOURCES) \ $(SOURCE_DIR)/MultiGridPreconWrapper.h $(SOURCE_DIR)/MultiGridPreconWrapper.cc \ -$(SOURCE_DIR)/LagrangeInterpolRestrict.h $(SOURCE_DIR)/LagrangeInterpolRestrict.cc \ +$(SOURCE_DIR)/LagrangeInterpolRestrict.h \ $(SOURCE_DIR)/BiCGStab.h $(SOURCE_DIR)/BiCGStab.hh \ $(SOURCE_DIR)/BiCGStab2.h $(SOURCE_DIR)/BiCGStab2.hh \ $(SOURCE_DIR)/InterpolRestrictMatrix.h $(SOURCE_DIR)/InterpolRestrictMatrix.cc \ @@ -84,7 +84,8 @@ $(SOURCE_DIR)/ProblemVec.h $(SOURCE_DIR)/ProblemVec.cc \ $(SOURCE_DIR)/DualTraverse.h $(SOURCE_DIR)/DualTraverse.cc \ $(SOURCE_DIR)/ElementPartition_ED.h $(SOURCE_DIR)/SurfacePartition_ED.h \ $(SOURCE_DIR)/ElementData.h $(SOURCE_DIR)/ElementData.cc \ -$(SOURCE_DIR)/CreatorMap.h $(SOURCE_DIR)/CreatorMap.cc $(SOURCE_DIR)/CreatorInterface.h \ +$(SOURCE_DIR)/CreatorMap.h $(SOURCE_DIR)/CreatorMap.hh $(SOURCE_DIR)/CreatorMap.cc \ +$(SOURCE_DIR)/CreatorInterface.h \ $(SOURCE_DIR)/ElementFunction.h \ $(SOURCE_DIR)/ProblemInterpolScal.h $(SOURCE_DIR)/ProblemInterpolScal.cc \ $(SOURCE_DIR)/ProblemInterpolVec.h $(SOURCE_DIR)/ProblemInterpolVec.cc \ diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in index 36ffced96d669ad8022eb77f9ba141a6fb6fbd67..d80d44a906c9e319d01e8468cbb3e00729a1950f 100644 --- a/AMDiS/bin/Makefile.in +++ b/AMDiS/bin/Makefile.in @@ -81,7 +81,6 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ConditionalEstimator.h \ $(SOURCE_DIR)/MultiGridPreconWrapper.h \ $(SOURCE_DIR)/MultiGridPreconWrapper.cc \ $(SOURCE_DIR)/LagrangeInterpolRestrict.h \ - $(SOURCE_DIR)/LagrangeInterpolRestrict.cc \ $(SOURCE_DIR)/BiCGStab.h $(SOURCE_DIR)/BiCGStab.hh \ $(SOURCE_DIR)/BiCGStab2.h $(SOURCE_DIR)/BiCGStab2.hh \ $(SOURCE_DIR)/InterpolRestrictMatrix.h \ @@ -113,8 +112,8 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ConditionalEstimator.h \ $(SOURCE_DIR)/ElementPartition_ED.h \ $(SOURCE_DIR)/SurfacePartition_ED.h \ $(SOURCE_DIR)/ElementData.h $(SOURCE_DIR)/ElementData.cc \ - $(SOURCE_DIR)/CreatorMap.h $(SOURCE_DIR)/CreatorMap.cc \ - $(SOURCE_DIR)/CreatorInterface.h \ + $(SOURCE_DIR)/CreatorMap.h $(SOURCE_DIR)/CreatorMap.hh \ + $(SOURCE_DIR)/CreatorMap.cc $(SOURCE_DIR)/CreatorInterface.h \ $(SOURCE_DIR)/ElementFunction.h \ $(SOURCE_DIR)/ProblemInterpolScal.h \ $(SOURCE_DIR)/ProblemInterpolScal.cc \ @@ -260,7 +259,6 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ConditionalEstimator.h \ @USE_PARALLEL_AMDIS_TRUE@ libamdis_la-ParMetisPartitioner.lo am_libamdis_la_OBJECTS = $(am__objects_1) \ libamdis_la-MultiGridPreconWrapper.lo \ - libamdis_la-LagrangeInterpolRestrict.lo \ libamdis_la-InterpolRestrictMatrix.lo \ libamdis_la-DOFIndexed.lo libamdis_la-GNUPlotWriter.lo \ libamdis_la-StlVector.lo libamdis_la-V3Vector.lo \ @@ -501,7 +499,7 @@ INCLUDES = $(AMDIS_INCLUDES) $(PARALLEL_INCLUDES) libamdis_la_SOURCES = \ $(PARALLEL_AMDIS_SOURCES) \ $(SOURCE_DIR)/MultiGridPreconWrapper.h $(SOURCE_DIR)/MultiGridPreconWrapper.cc \ -$(SOURCE_DIR)/LagrangeInterpolRestrict.h $(SOURCE_DIR)/LagrangeInterpolRestrict.cc \ +$(SOURCE_DIR)/LagrangeInterpolRestrict.h \ $(SOURCE_DIR)/BiCGStab.h $(SOURCE_DIR)/BiCGStab.hh \ $(SOURCE_DIR)/BiCGStab2.h $(SOURCE_DIR)/BiCGStab2.hh \ $(SOURCE_DIR)/InterpolRestrictMatrix.h $(SOURCE_DIR)/InterpolRestrictMatrix.cc \ @@ -527,7 +525,8 @@ $(SOURCE_DIR)/ProblemVec.h $(SOURCE_DIR)/ProblemVec.cc \ $(SOURCE_DIR)/DualTraverse.h $(SOURCE_DIR)/DualTraverse.cc \ $(SOURCE_DIR)/ElementPartition_ED.h $(SOURCE_DIR)/SurfacePartition_ED.h \ $(SOURCE_DIR)/ElementData.h $(SOURCE_DIR)/ElementData.cc \ -$(SOURCE_DIR)/CreatorMap.h $(SOURCE_DIR)/CreatorMap.cc $(SOURCE_DIR)/CreatorInterface.h \ +$(SOURCE_DIR)/CreatorMap.h $(SOURCE_DIR)/CreatorMap.hh $(SOURCE_DIR)/CreatorMap.cc \ +$(SOURCE_DIR)/CreatorInterface.h \ $(SOURCE_DIR)/ElementFunction.h \ $(SOURCE_DIR)/ProblemInterpolScal.h $(SOURCE_DIR)/ProblemInterpolScal.cc \ $(SOURCE_DIR)/ProblemInterpolVec.h $(SOURCE_DIR)/ProblemInterpolVec.cc \ @@ -816,7 +815,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-InterpolRestrictMatrix.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-JacobiSmoother.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Lagrange.Plo@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-LagrangeInterpolRestrict.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-LeafData.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Line.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MacroElement.Plo@am__quote@ @@ -944,13 +942,6 @@ libamdis_la-MultiGridPreconWrapper.lo: $(SOURCE_DIR)/MultiGridPreconWrapper.cc @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-MultiGridPreconWrapper.lo `test -f '$(SOURCE_DIR)/MultiGridPreconWrapper.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/MultiGridPreconWrapper.cc -libamdis_la-LagrangeInterpolRestrict.lo: $(SOURCE_DIR)/LagrangeInterpolRestrict.cc -@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-LagrangeInterpolRestrict.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-LagrangeInterpolRestrict.Tpo" -c -o libamdis_la-LagrangeInterpolRestrict.lo `test -f '$(SOURCE_DIR)/LagrangeInterpolRestrict.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/LagrangeInterpolRestrict.cc; \ -@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-LagrangeInterpolRestrict.Tpo" "$(DEPDIR)/libamdis_la-LagrangeInterpolRestrict.Plo"; else rm -f "$(DEPDIR)/libamdis_la-LagrangeInterpolRestrict.Tpo"; exit 1; fi -@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/LagrangeInterpolRestrict.cc' object='libamdis_la-LagrangeInterpolRestrict.lo' libtool=yes @AMDEPBACKSLASH@ -@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ -@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-LagrangeInterpolRestrict.lo `test -f '$(SOURCE_DIR)/LagrangeInterpolRestrict.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/LagrangeInterpolRestrict.cc - libamdis_la-InterpolRestrictMatrix.lo: $(SOURCE_DIR)/InterpolRestrictMatrix.cc @am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-InterpolRestrictMatrix.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-InterpolRestrictMatrix.Tpo" -c -o libamdis_la-InterpolRestrictMatrix.lo `test -f '$(SOURCE_DIR)/InterpolRestrictMatrix.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/InterpolRestrictMatrix.cc; \ @am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-InterpolRestrictMatrix.Tpo" "$(DEPDIR)/libamdis_la-InterpolRestrictMatrix.Plo"; else rm -f "$(DEPDIR)/libamdis_la-InterpolRestrictMatrix.Tpo"; exit 1; fi diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h index a8ce1710b927eea83d99e0cf7fe1711df6ad901e..262b7c8b51b7441d3d705ed04a6bc4a0adce8cb7 100644 --- a/AMDiS/src/AMDiS.h +++ b/AMDiS/src/AMDiS.h @@ -85,7 +85,9 @@ #include "RefinementManager2d.h" #include "RefinementManager3d.h" #include "RobinBC.h" +#include "SmootherBase.h" #include "SolutionDataStorage.h" +#include "SparseVector.h" #include "SurfaceOperator.h" #include "SurfaceQuadrature.h" #include "SystemVector.h" diff --git a/AMDiS/src/CoarseningManager2d.cc b/AMDiS/src/CoarseningManager2d.cc index a54e24c6a6fc5970db040e14cd34bffea7997989..250bb4c20b27aec9c023ca7f5e2d71462f3a2df3 100644 --- a/AMDiS/src/CoarseningManager2d.cc +++ b/AMDiS/src/CoarseningManager2d.cc @@ -19,41 +19,38 @@ namespace AMDiS { void CoarseningManager2d::coarsenTriangle(Triangle *el) { - FUNCNAME("CoarseningManager2d::coarseTriangle"); - Triangle *child[2]; + FUNCNAME("CoarseningManager2d::coarseTriangle()"); - child[0] = dynamic_cast<Triangle*>(const_cast<Element*>( el->getChild(0))); - child[1] = dynamic_cast<Triangle*>(const_cast<Element*>( el->getChild(1))); + Triangle *child[2]; + child[0] = dynamic_cast<Triangle*>(const_cast<Element*>(el->getChild(0))); + child[1] = dynamic_cast<Triangle*>(const_cast<Element*>(el->getChild(1))); TEST_EXIT_DBG(child[0]->getMark() < 0 && child[1]->getMark() < 0) ("element %d with children[%d,%d] must not be coarsend!\n", el->getIndex(), child[0]->getIndex(), child[1]->getIndex()); - if (mesh->getNumberOfDOFs(EDGE)) - { - /****************************************************************************/ - /* remove dof from common edge of child[0] and child[1] */ - /****************************************************************************/ - mesh->freeDOF(const_cast<int*>( child[0]->getDOF(4)), EDGE); - } - - if (mesh->getNumberOfDOFs(CENTER)) - { - /****************************************************************************/ - /* remove dof from the barycenters of child[0] and child[1] */ - /****************************************************************************/ - int node = mesh->getNode(CENTER); - - mesh->freeDOF(const_cast<int*>( child[0]->getDOF(node)), CENTER); - mesh->freeDOF(const_cast<int*>( child[1]->getDOF(node)), CENTER); - } + if (mesh->getNumberOfDOFs(EDGE)) { + /****************************************************************************/ + /* remove dof from common edge of child[0] and child[1] */ + /****************************************************************************/ + mesh->freeDOF(const_cast<int*>(child[0]->getDOF(4)), EDGE); + } + if (mesh->getNumberOfDOFs(CENTER)) { + /****************************************************************************/ + /* remove dof from the barycenters of child[0] and child[1] */ + /****************************************************************************/ + int node = mesh->getNode(CENTER); + + mesh->freeDOF(const_cast<int*>(child[0]->getDOF(node)), CENTER); + mesh->freeDOF(const_cast<int*>(child[1]->getDOF(node)), CENTER); + } + el->coarsenElementData(child[0], child[1]); el->setFirstChild(NULL); el->setSecondChild(NULL); - mesh->freeElement(child[0]); mesh->freeElement(child[1]); @@ -72,8 +69,8 @@ namespace AMDiS { /****************************************************************************/ void CoarseningManager2d::coarsenPatch(RCNeighbourList *coarsenList, - int n_neigh, - int bound) + int n_neigh, + int bound) { Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( coarsenList->getElement(0))); Triangle *neigh = dynamic_cast<Triangle*>(const_cast<Element*>( coarsenList->getElement(1))); @@ -195,9 +192,7 @@ namespace AMDiS { if (coarse_list.doCoarsePatch(n_neigh)) { int n_neigh_periodic; - DegreeOfFreedom *next_edge[2]; - RCNeighbourList *periodicList; while (edge[0] != NULL) { diff --git a/AMDiS/src/CreatorInterface.h b/AMDiS/src/CreatorInterface.h index 687827e8046cbf1a6b0c4a1145eefdcc4c429fe0..ea5e3a4674737e377fe1cf5fb344028f509ddbd8 100644 --- a/AMDiS/src/CreatorInterface.h +++ b/AMDiS/src/CreatorInterface.h @@ -40,7 +40,7 @@ namespace AMDiS { class CreatorInterface { public: - virtual ~CreatorInterface() {}; + virtual ~CreatorInterface() {} /** \brief * Must be implemented by sub classes of CreatorInterface. @@ -51,14 +51,14 @@ namespace AMDiS { /** \brief * Can be implemented by sub classes. */ - virtual void free(BaseClass *) {}; + virtual void free(BaseClass *) {} /** \brief * */ virtual bool isNullCreator() { return false; - }; + } }; /** \brief @@ -73,14 +73,14 @@ namespace AMDiS { */ BaseClass* create() { return NULL; - }; + } /** \brief * */ virtual bool isNullCreator() { return true; - }; + } }; } diff --git a/AMDiS/src/CreatorMap.cc b/AMDiS/src/CreatorMap.cc index 5ebfa498fc6079d9c8927764266c2830e261a7c1..d945a5ddb57a669f8d6624713611dbaa55ec872a 100644 --- a/AMDiS/src/CreatorMap.cc +++ b/AMDiS/src/CreatorMap.cc @@ -261,5 +261,4 @@ namespace AMDiS { creator = NEW ElementRegion_ED::Creator; addCreator("ElementRegion_ED", creator); } - } diff --git a/AMDiS/src/CreatorMap.h b/AMDiS/src/CreatorMap.h index 33deb028bbbe12304b840d76aa8561fdb7fb1924..7807ddc8b78dab65bf6058d2273d3444aae97601 100644 --- a/AMDiS/src/CreatorMap.h +++ b/AMDiS/src/CreatorMap.h @@ -39,63 +39,66 @@ namespace AMDiS { * to the map. */ template<typename BaseClass> - class CreatorMap + class CreatorMap + { + public: + /** \brief + * Adds a new creator together with the given key to the map. + */ + static void addCreator(std::string key, CreatorInterface<BaseClass>* creator) { - public: - /** \brief - * Adds a new creator together with the given key to the map. - */ - static void addCreator(std::string key, CreatorInterface<BaseClass>* creator) - { - FUNCNAME("CreatorMap::addCreator()"); - init(); - TEST_EXIT(creatorMap[key] == NULL) - ("there is already a creator for key %s\n",key.c_str()); - creatorMap[key] = creator; - }; - - /** \brief - * Creates a object of the type corresponding to key. - */ - static CreatorInterface<BaseClass>* getCreator(std::string key) { - FUNCNAME("CreatorMap::getCreator()"); - init(); - CreatorInterface<BaseClass> *creator = creatorMap[key]; - TEST_EXIT(creator)("no creator for key %s\n", key.c_str()); - return creator; - }; - - static void addDefaultCreators(); - - protected: - /** \brief - * Constructor is protected because derived maps should be singleton. - */ - static void init() { - if(!initialized) { - initialized = true; - NullCreator<BaseClass> *nullCreator = new NullCreator<BaseClass>; - addCreator("no", nullCreator); - addCreator("0", nullCreator); - addDefaultCreators(); - } - }; - - protected: - /** \brief - * STL map containing the pairs of keys and creators. - */ - static std::map< std::string, CreatorInterface<BaseClass>* > creatorMap; - - static bool initialized; - }; + FUNCNAME("CreatorMap::addCreator()"); + init(); + TEST_EXIT(creatorMap[key] == NULL) + ("there is already a creator for key %s\n",key.c_str()); + creatorMap[key] = creator; + } + + /** \brief + * Creates a object of the type corresponding to key. + */ + static CreatorInterface<BaseClass>* getCreator(std::string key) { + FUNCNAME("CreatorMap::getCreator()"); + init(); + CreatorInterface<BaseClass> *creator = creatorMap[key]; + TEST_EXIT(creator)("no creator for key %s\n", key.c_str()); + return creator; + } + + static void clear(); + + static void addDefaultCreators(); + + protected: + /** \brief + * Constructor is protected because derived maps should be singleton. + */ + static void init() { + if (!initialized) { + initialized = true; + NullCreator<BaseClass> *nullCreator = new NullCreator<BaseClass>; + addCreator("no", nullCreator); + addDefaultCreators(); + } + } + + protected: + /** \brief + * STL map containing the pairs of keys and creators. + */ + static std::map< std::string, CreatorInterface<BaseClass>* > creatorMap; + + static bool initialized; + }; template<typename BaseClass> - std::map< std::string, CreatorInterface<BaseClass>* > CreatorMap<BaseClass>::creatorMap; + std::map< std::string, CreatorInterface<BaseClass>* > CreatorMap<BaseClass>::creatorMap; template<typename BaseClass> - bool CreatorMap<BaseClass>::initialized = false; + bool CreatorMap<BaseClass>::initialized = false; } +#include "CreatorMap.hh" + #endif diff --git a/AMDiS/src/CreatorMap.hh b/AMDiS/src/CreatorMap.hh new file mode 100644 index 0000000000000000000000000000000000000000..b82231017634d1bbcc3c1c17aeaf38041e094264 --- /dev/null +++ b/AMDiS/src/CreatorMap.hh @@ -0,0 +1,18 @@ +#include "MemoryManager.h" +#include <map> + +namespace AMDiS { + + template<typename BaseClass> + void CreatorMap<BaseClass>::clear() + { + typename std::map< std::string, CreatorInterface<BaseClass>* >::iterator it; + + for (it = creatorMap.begin(); + it != creatorMap.end(); + ++it) { + DELETE it->second; + } + } + +} diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc index ae0d0ec85f16ad5b7545993d34097a123d5b5371..28f2427b4ed357e1b68776349b8eb400f5998549 100644 --- a/AMDiS/src/DOFMatrix.cc +++ b/AMDiS/src/DOFMatrix.cc @@ -61,10 +61,15 @@ namespace AMDiS { { FUNCNAME("DOFMatrix::~DOFMatrix()"); + matrix.clear(); + if (rowFESpace && rowFESpace->getAdmin()) { (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this); } + if (boundaryManager) + DELETE boundaryManager; + if (elementMatrix) DELETE elementMatrix; } diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc index 047f2b40df1d155d7b5e308fcbe7a755397dfda2..0b3150cb3b535e71110a27f72865234fb53032f3 100644 --- a/AMDiS/src/Element.cc +++ b/AMDiS/src/Element.cc @@ -7,18 +7,18 @@ namespace AMDiS { + std::map<DegreeOfFreedom*, bool> Element::deletedDOFs; int Element::getRegion() const { - ElementRegion_ED* red_; - if (!elementData) return -1; - red_ = dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION)); + ElementRegion_ED* red = dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION)); - if (red_) - return red_->getRegion(); + if (red) { + return red->getRegion(); + } return -1; } @@ -35,12 +35,9 @@ namespace AMDiS { Element::Element(Mesh *aMesh) { mesh = aMesh; - index = mesh ? mesh->getNextElementIndex() : -1; - child[0] = NULL; child[1] = NULL; - newCoord = NULL; elementData = NULL; @@ -53,7 +50,7 @@ namespace AMDiS { // call destructor through Mesh::freeElement !!! Element::~Element() - { + { if (child[0]) { DELETE child[0]; } @@ -64,6 +61,29 @@ namespace AMDiS { if (newCoord) { DELETE newCoord; } + + if (elementData) { + elementData->deleteDecorated(); + DELETE elementData; + } + } + + bool Element::deleteElementData(int typeID) + { + FUNCNAME("Element::deleteElementData()"); + + if (elementData) { + if (elementData->isOfType(typeID)) { + ElementData *tmp = elementData; + elementData = elementData->getDecorated(); + DELETE tmp; + tmp = NULL; + return true; + } else { + return elementData->deleteDecorated(typeID); + } + } + return false; } void Element::deleteElementDOFs() @@ -81,22 +101,19 @@ namespace AMDiS { if (ndof > 0) { for (int i = 0; i < mesh->getGeo(position); i++) { - if (dof[j] != NULL) { - if (Mesh::serializedDOFs.count(dof[j][0]) == 0) { - // std::cout << "FREE INNER: " << ndof << std::endl; + if (dof[j]) { + if (deletedDOFs.count(dof[j]) == 0) { + deletedDOFs[dof[j]] = true; FREE_MEMORY(dof[j], DegreeOfFreedom, ndof); - Mesh::serializedDOFs[dof[j][0]] = NULL; } - } - + } j++; } } } - // std::cout << "FREE OUTER: " << mesh->getNumberOfNodes() << std::endl; FREE_MEMORY(dof, DegreeOfFreedom*, mesh->getNumberOfNodes()); - + if (child[0]) { child[0]->deleteElementDOFs(); } @@ -108,14 +125,16 @@ namespace AMDiS { Element* Element::cloneWithDOFs() { Element *el; - + if (isLine()) { - el = NEW Line(mesh); + el = NEW Line(NULL); } else if (isTriangle()) { - el = NEW Triangle(mesh); + el = NEW Triangle(NULL); } else { - el = NEW Tetrahedron(mesh); + el = NEW Tetrahedron(NULL); } + + el->mesh = mesh; el->index = index; el->mark = mark; if (newCoord) { @@ -123,9 +142,9 @@ namespace AMDiS { *nc = *newCoord; el->newCoord = nc; } - + /* =========== And here we clone the DOFs =========== */ - + el->dof = GET_MEMORY(DegreeOfFreedom*, mesh->getNumberOfNodes()); int dim = mesh->getDim(); @@ -158,7 +177,7 @@ namespace AMDiS { } } } - + /* =========== And clone the children ============= */ if (child[0]) { diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index ab62a7b4483ae91fa53dc28c8ed392dcb7d833f7..bd7769a2d4d158dbf625ea5cab168c21df30b192 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -464,7 +464,7 @@ namespace AMDiS { /** \brief * Coarsens Element's leaf data */ - inline void coarsenElementData(Element* child1, Element* child2, int elType=0) { + inline void coarsenElementData(Element* child1, Element* child2, int elType = 0) { ElementData *childData; childData = child1->getElementData(); if (childData) { @@ -498,22 +498,9 @@ namespace AMDiS { } /** \brief - * kills \ref elementData + * Deletes the \ref elementData with a specific typeID. */ - bool deleteElementData(int typeID) { - FUNCNAME("Element::deleteElementData()"); - if (elementData) { - if (elementData->isOfType(typeID)) { - ElementData *tmp = elementData; - elementData = elementData->getDecorated(); - DELETE tmp; - return true; - } else { - return elementData->deleteDecorated(typeID); - } - } - return false; - } + bool deleteElementData(int typeID); /** \brief * Returns whether element is refined at side side @@ -620,6 +607,12 @@ namespace AMDiS { ElementData* elementData; + /** \brief + * This map is used for deletion of all DOFs of all elements of a mesh. Once + * a DOF-vector (all DOFS at a node, edge, etc.) is deleted, its address is + * added to this map to note not to delete it a second time. + */ + static std::map<DegreeOfFreedom*, bool> deletedDOFs; friend class Mesh; }; diff --git a/AMDiS/src/ElementData.cc b/AMDiS/src/ElementData.cc index d48d41375d504e8453c07ab4c6e8d569afd7110e..c39fdbb3aeb14da95a2d4d3466b36a4c97542812 100644 --- a/AMDiS/src/ElementData.cc +++ b/AMDiS/src/ElementData.cc @@ -1,8 +1,50 @@ +#include "MemoryManager.h" #include "ElementData.h" namespace AMDiS { - void ElementData::serialize(std::ostream& out) { + void ElementData::coarsenElementData(Element* parent, + Element* thisChild, + Element* otherChild, + int elTypeParent) + { + if (decorated_) { + decorated_->coarsenElementData(parent, thisChild, otherChild, elTypeParent); + delete decorated_; + decorated_ = NULL; + } + } + + bool ElementData::deleteDecorated(int typeID) + { + if (decorated_) { + if (decorated_->isOfType(typeID)) { + ElementData *tmp = decorated_; + decorated_ = decorated_->decorated_; + delete tmp; + tmp = NULL; + return true; + } else { + return decorated_->deleteDecorated(typeID); + } + } + return false; + } + + void ElementData::deleteDecorated() + { + if (decorated_) { + decorated_->deleteDecorated(); + DELETE decorated_; + } + } + + ElementData::~ElementData() + { + } + + void ElementData::serialize(std::ostream& out) + { std::string decoratedType; if (decorated_) { decoratedType = decorated_->getTypeName(); @@ -13,7 +55,8 @@ namespace AMDiS { } } - void ElementData::deserialize(std::istream& in) { + void ElementData::deserialize(std::istream& in) + { TEST_EXIT(decorated_ == NULL) ("there are already decorated element data\n"); std::string decoratedType; @@ -24,7 +67,7 @@ namespace AMDiS { decorated_->deserialize(in); } else { decorated_ = NULL; - }; + } } } diff --git a/AMDiS/src/ElementData.h b/AMDiS/src/ElementData.h index 6eae31603825dfde77d3d0f335dc390bd8c24ac5..a4e23eb249fe385b001ddad198f8b8c26a9796ee 100644 --- a/AMDiS/src/ElementData.h +++ b/AMDiS/src/ElementData.h @@ -56,7 +56,7 @@ namespace AMDiS { /** \brief * destructor */ - virtual ~ElementData() {} + virtual ~ElementData(); /** \brief * Refinement of parent to child1 and child2. @@ -85,23 +85,13 @@ namespace AMDiS { virtual void coarsenElementData(Element* parent, Element* thisChild, Element* otherChild, - int elTypeParent) - { - if(decorated_) { - decorated_->coarsenElementData(parent, - thisChild, - otherChild, - elTypeParent); - delete decorated_; - decorated_ = NULL; - } - } + int elTypeParent); /** \brief * Returns a copy of this ElementData object including all decorated data. */ virtual ElementData *clone() const { - if(decorated_) { + if (decorated_) { return decorated_->clone(); } return NULL; @@ -149,19 +139,16 @@ namespace AMDiS { return NULL; } - inline bool deleteDecorated(int typeID) { - if(decorated_) { - if(decorated_->isOfType(typeID)) { - ElementData *tmp = decorated_; - decorated_ = decorated_->decorated_; - delete tmp; - return true; - } else { - return decorated_->deleteDecorated(typeID); - } - } - return false; - } + /** \ref + * Search the \ref decorated_ chain for a specific type ID, and delets + * this entry. + */ + bool deleteDecorated(int typeID); + + /** \ref + * Delets the whole \ref decorated_ chain. + */ + void deleteDecorated(); inline ElementData *getDecorated() { return decorated_; diff --git a/AMDiS/src/FiniteElemSpace.cc b/AMDiS/src/FiniteElemSpace.cc index b121e53719d756215e9e7d4c8db67aeeb03edbcc..0370144fda3e9b7ef0d07519116097dee1a283d0 100644 --- a/AMDiS/src/FiniteElemSpace.cc +++ b/AMDiS/src/FiniteElemSpace.cc @@ -28,8 +28,8 @@ namespace AMDiS { TEST_EXIT(ndof)("no n_dof or basFcts->n_dof\n"); for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) { admin_local = &(mesh->getDOFAdmin(i)); - int j; - for (j = 0; j <= mesh->getDim(); j++) { + int j = 0; + for (; j <= mesh->getDim(); j++) { if (admin_local->getNumberOfDOFs(j) != (*ndof)[j]) break; } @@ -50,6 +50,10 @@ namespace AMDiS { { } + FiniteElemSpace::~FiniteElemSpace() + { + } + FiniteElemSpace& FiniteElemSpace::operator=(const FiniteElemSpace& feSpace) { if (&feSpace == this) @@ -93,4 +97,10 @@ namespace AMDiS { } } + void FiniteElemSpace::clear() + { + for (int i = 0; i < feSpaces.size(); i++) { + DELETE feSpaces[i]; + } + } } diff --git a/AMDiS/src/FiniteElemSpace.h b/AMDiS/src/FiniteElemSpace.h index 95701ad260e4fecf1d989fec9928cb6ba5dad244..55c061d031e69e699d3599c90fe2d176e6d39cd6 100644 --- a/AMDiS/src/FiniteElemSpace.h +++ b/AMDiS/src/FiniteElemSpace.h @@ -74,7 +74,7 @@ namespace AMDiS { /** \brief * destructor */ - virtual ~FiniteElemSpace() {} + ~FiniteElemSpace(); FiniteElemSpace& operator=(const FiniteElemSpace& feSpace); @@ -105,6 +105,8 @@ namespace AMDiS { inline Mesh* getMesh() const { return mesh; } + + static void clear(); protected: /** \brief diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc index b504aaf33f1fe4bf097fa9073468ac4d9943bdde..62b4a203c176ee12cf3139a9d45a79f2b24d059b 100644 --- a/AMDiS/src/Global.cc +++ b/AMDiS/src/Global.cc @@ -322,7 +322,14 @@ namespace AMDiS { // set msgWait Parameters::getGlobalParameter(0, "WAIT", "%d", &d); Msg::setMsgWait(!(d == 0)); - }; + } + + void Global::clear() + { + DELETE referenceElement[1]; + DELETE referenceElement[2]; + DELETE referenceElement[3]; + } int fac(int i) { diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h index d497dd664aba8e0df7d61030649d33ef42ee12e7..35b3170d1adc4a3c1d690876a58c6f284b8a4096 100644 --- a/AMDiS/src/Global.h +++ b/AMDiS/src/Global.h @@ -474,6 +474,8 @@ namespace AMDiS { */ static void init(); + static void clear(); + private: /** \brief * Global is a pure static class. So the constructor is private to avoid diff --git a/AMDiS/src/InterpolRestrictMatrix.h b/AMDiS/src/InterpolRestrictMatrix.h index 00f9a857ddf385cd783fc06a33f4b8aa50e5c099..d45ec5ee749db4d45bce3a40ab2275c2341e42ac 100644 --- a/AMDiS/src/InterpolRestrictMatrix.h +++ b/AMDiS/src/InterpolRestrictMatrix.h @@ -50,7 +50,7 @@ namespace AMDiS { void print(); protected: - virtual void fillMemory() {}; + virtual void fillMemory() {} protected: int numFcts_; diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc index 872d1437be868bf1d4c919c3c50e9f68855933f7..42b3f42a3d23c5d52f36a55141c50caab85cba1b 100644 --- a/AMDiS/src/Lagrange.cc +++ b/AMDiS/src/Lagrange.cc @@ -49,14 +49,22 @@ namespace AMDiS { Lagrange::~Lagrange() { + std::cout << "DELETE LAGRANGE" << std::endl; + std::cout << "BARY ADR = " << bary << std::endl; + for (int i = 0; i < static_cast<int>(bary->size()); i++) { + if ((*bary)[i]) { + std::cout << "DEL ADR = " << (*bary)[i] << std::endl; + DELETE (*bary)[i]; + (*bary)[i] = NULL; + } + } } Lagrange* Lagrange::getLagrange(int dim, int degree) { std::list<Lagrange*>::iterator it; - - for(it = allBasFcts.begin(); it != allBasFcts.end(); it++) { - if(((*it)->dim == dim) && ((*it)->degree == degree)) { + for (it = allBasFcts.begin(); it != allBasFcts.end(); it++) { + if (((*it)->dim == dim) && ((*it)->degree == degree)) { return (*it); } } @@ -66,6 +74,14 @@ namespace AMDiS { return newLagrange; } + void Lagrange::clear() + { + std::list<Lagrange*>::iterator it; + for (it = allBasFcts.begin(); it != allBasFcts.end(); it++) { + DELETE *it; + } + } + void Lagrange::setFunctionPointer() { if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) { @@ -631,10 +647,12 @@ namespace AMDiS { vec = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0); } - if (dimIndex == numCoords-1) { - (*vec)[coordInd[dimIndex]] = double(rest)/degree; + if (dimIndex == numCoords - 1) { + (*vec)[coordInd[dimIndex]] = double(rest) / degree; DimVec<double>* newCoords = NEW DimVec<double>(*vec); + std::cout << "BARY ADR = " << bary << std::endl; bary->push_back(newCoords); + std::cout << "ADD ADR = " << newCoords << std::endl; } else { for (int i = rest - 1; i >= 1; i--) { (*vec)[coordInd[dimIndex]] = double(i) / degree; @@ -645,14 +663,14 @@ namespace AMDiS { void Lagrange::setBary() { - bary = &baryDimDegree[dim][degree]; + bary = &baryDimDegree[dim][degree]; if (static_cast<int>(bary->size()) == 0) { for (int i = 0; i <= dim; i++) { // for all positions int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim); for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/... int *coordInd = GET_MEMORY(int, i + 1); // indices of relevant coords - for (int k = 0; k < i+1; k++) { + for (int k = 0; k < i + 1; k++) { coordInd[k] = Global::getReferenceElement(dim)-> getVertexOfPosition(INDEX_OF_DIM(i,dim), j, k); } diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h index a977e2fac6f394aa75aa2549631cdb56ae94d7f3..c57b691159fabc6a7d102dd116de2529db24cbbd 100644 --- a/AMDiS/src/Lagrange.h +++ b/AMDiS/src/Lagrange.h @@ -175,6 +175,8 @@ namespace AMDiS { AbstractFunction<WorldVector<double>, WorldVector<double> >* f, DOFVector<WorldVector<double> >* fh); + static void clear(); + protected: /** \brief * sets the barycentric coordinates (stored in \ref bary) of the local diff --git a/AMDiS/src/LagrangeInterpolRestrict.cc b/AMDiS/src/LagrangeInterpolRestrict.cc deleted file mode 100644 index a68321ef2319d2d8fca463d26ddef5a3bb0952a4..0000000000000000000000000000000000000000 --- a/AMDiS/src/LagrangeInterpolRestrict.cc +++ /dev/null @@ -1,26 +0,0 @@ -#include "LagrangeInterpolRestrict.h" - -namespace AMDiS { - - InterpolRestrictMatrix *interpolMatrices[3][4] = { - {new IM_1d_1, new IM_1d_2, NULL, NULL}, - {new IM_2d_1, new IM_2d_2, new IM_2d_3, new IM_2d_4}, - {new IM_3d_1, new IM_3d_2, new IM_3d_3, new IM_3d_4} - }; - - InterpolRestrictMatrix *restrictMatrices[3][4] = { - { new InterpolRestrictMatrix(*(interpolMatrices[0][0]), true), - new InterpolRestrictMatrix(*(interpolMatrices[0][1]), true), - NULL, - NULL}, - { new InterpolRestrictMatrix(*(interpolMatrices[1][0]), true), - new InterpolRestrictMatrix(*(interpolMatrices[1][1]), true), - new InterpolRestrictMatrix(*(interpolMatrices[1][2]), true), - new InterpolRestrictMatrix(*(interpolMatrices[1][3]), true) }, - { new InterpolRestrictMatrix(*(interpolMatrices[2][0]), true), - new InterpolRestrictMatrix(*(interpolMatrices[2][1]), true), - new InterpolRestrictMatrix(*(interpolMatrices[2][2]), true), - new InterpolRestrictMatrix(*(interpolMatrices[2][3]), true) } - }; - -} diff --git a/AMDiS/src/LagrangeInterpolRestrict.h b/AMDiS/src/LagrangeInterpolRestrict.h index a4740d5f6dcdabf1f4930880c037a293e4e0d9c4..fb85fba81095a29552310fc84da64ebdfb3ecda7 100644 --- a/AMDiS/src/LagrangeInterpolRestrict.h +++ b/AMDiS/src/LagrangeInterpolRestrict.h @@ -12,7 +12,9 @@ namespace AMDiS { class IM_1d_1 : public InterpolRestrictMatrix { public: - IM_1d_1() : InterpolRestrictMatrix(2, 3, 2, false) {}; + IM_1d_1() : + InterpolRestrictMatrix(2, 3, 2, false) + {} protected: void fillMemory() @@ -30,7 +32,7 @@ namespace AMDiS { colDOFPtrs_[0] = parentDOFs_ + 0; colDOFPtrs_[1] = parentDOFs_ + 1; - }; + } }; @@ -41,7 +43,9 @@ namespace AMDiS { class IM_1d_2 : public InterpolRestrictMatrix { public: - IM_1d_2() : InterpolRestrictMatrix(3, 5, 3, false) {}; + IM_1d_2() : + InterpolRestrictMatrix(3, 5, 3, false) + {} protected: void fillMemory() @@ -68,7 +72,7 @@ namespace AMDiS { colDOFPtrs_[0] = parentDOFs_ + 0; colDOFPtrs_[1] = parentDOFs_ + 1; colDOFPtrs_[2] = parentDOFs_ + 2; - }; + } }; @@ -79,7 +83,9 @@ namespace AMDiS { class IM_1d_4 : public InterpolRestrictMatrix { public: - IM_1d_4() : InterpolRestrictMatrix(5, 9, 5, false) {}; + IM_1d_4() : + InterpolRestrictMatrix(5, 9, 5, false) + {} protected: void fillMemory() @@ -97,7 +103,7 @@ namespace AMDiS { colDOFPtrs_[0] = parentDOFs_ + 0; colDOFPtrs_[1] = parentDOFs_ + 1; - }; + } }; @@ -108,7 +114,9 @@ namespace AMDiS { class IM_2d_1 : public InterpolRestrictMatrix { public: - IM_2d_1() : InterpolRestrictMatrix(3, 4, 3, false) {}; + IM_2d_1() : + InterpolRestrictMatrix(3, 4, 3, false) + {} void fillMemory() { @@ -128,7 +136,7 @@ namespace AMDiS { colDOFPtrs_[0] = parentDOFs_ + 0; colDOFPtrs_[1] = parentDOFs_ + 1; colDOFPtrs_[2] = parentDOFs_ + 2; - }; + } }; @@ -139,7 +147,9 @@ namespace AMDiS { class IM_2d_2 : public InterpolRestrictMatrix { public: - IM_2d_2() : InterpolRestrictMatrix(6, 9, 6, false) {}; + IM_2d_2() : + InterpolRestrictMatrix(6, 9, 6, false) + {} void fillMemory() { @@ -181,7 +191,7 @@ namespace AMDiS { colDOFPtrs_[3] = parentDOFs_ + 3; colDOFPtrs_[4] = parentDOFs_ + 4; colDOFPtrs_[5] = parentDOFs_ + 5; - }; + } }; @@ -192,7 +202,9 @@ namespace AMDiS { class IM_2d_3 : public InterpolRestrictMatrix { public: - IM_2d_3() : InterpolRestrictMatrix(10, 16, 10, false) {}; + IM_2d_3() : + InterpolRestrictMatrix(10, 16, 10, false) + {} void fillMemory() { @@ -276,7 +288,7 @@ namespace AMDiS { colDOFPtrs_[7] = parentDOFs_ + 7; colDOFPtrs_[8] = parentDOFs_ + 8; colDOFPtrs_[9] = parentDOFs_ + 9; - }; + } }; @@ -287,7 +299,9 @@ namespace AMDiS { class IM_2d_4 : public InterpolRestrictMatrix { public: - IM_2d_4() : InterpolRestrictMatrix(15, 25, 15, false) {}; + IM_2d_4() : + InterpolRestrictMatrix(15, 25, 15, false) + {} void fillMemory() { @@ -444,7 +458,7 @@ namespace AMDiS { colDOFPtrs_[12] = parentDOFs_ + 12; colDOFPtrs_[13] = parentDOFs_ + 13; colDOFPtrs_[14] = parentDOFs_ + 14; - }; + } }; @@ -456,7 +470,9 @@ namespace AMDiS { class IM_3d_1 : public InterpolRestrictMatrix { public: - IM_3d_1() : InterpolRestrictMatrix(4, 5, 4, false) {}; + IM_3d_1() : + InterpolRestrictMatrix(4, 5, 4, false) + {} void fillMemory() { @@ -479,7 +495,7 @@ namespace AMDiS { colDOFPtrs_[1] = parentDOFs_ + 1; colDOFPtrs_[2] = parentDOFs_ + 2; colDOFPtrs_[3] = parentDOFs_ + 3; - }; + } }; @@ -491,7 +507,9 @@ namespace AMDiS { class IM_3d_2 : public InterpolRestrictMatrix { public: - IM_3d_2() : InterpolRestrictMatrix(10, 14, 10, false) {}; + IM_3d_2() : + InterpolRestrictMatrix(10, 14, 10, false) + {} void fillMemory() { @@ -553,7 +571,7 @@ namespace AMDiS { colDOFPtrs_[7] = parentDOFs_ + 7; colDOFPtrs_[8] = parentDOFs_ + 8; colDOFPtrs_[9] = parentDOFs_ + 9; - }; + } }; @@ -565,7 +583,9 @@ namespace AMDiS { class IM_3d_3 : public InterpolRestrictMatrix { public: - IM_3d_3() : InterpolRestrictMatrix(20, 30, 20, true) {}; + IM_3d_3() : + InterpolRestrictMatrix(20, 30, 20, true) + {} void fillMemory() { @@ -774,7 +794,7 @@ namespace AMDiS { type0_colDOFPtrs_[17] = parentDOFs_ + 17; type0_colDOFPtrs_[18] = parentDOFs_ + 18; type0_colDOFPtrs_[19] = parentDOFs_ + 19; - }; + } }; @@ -785,7 +805,9 @@ namespace AMDiS { class IM_3d_4 : public InterpolRestrictMatrix { public: - IM_3d_4() : InterpolRestrictMatrix(35, 55, 35, true) {}; + IM_3d_4() : + InterpolRestrictMatrix(35, 55, 35, true) + {} void fillMemory() { @@ -1254,20 +1276,8 @@ namespace AMDiS { type0_colDOFPtrs_[32] = parentDOFs_ + 32; type0_colDOFPtrs_[33] = parentDOFs_ + 33; type0_colDOFPtrs_[34] = parentDOFs_ + 34; - }; + } }; - - /** \brief - * interpolMatrices[n][m] stores the interpolation matrix for - * n-dimensional elements and m-dimensional basis functions. - */ - extern InterpolRestrictMatrix *interpolMatrices[3][4]; - - /** \brief - * interpolMatrices[n][m] stores the restriction matrix for - * n-dimensional elements and m-dimensional basis functions. - */ - extern InterpolRestrictMatrix *restrictMatrices[3][4]; } #endif diff --git a/AMDiS/src/LeafData.cc b/AMDiS/src/LeafData.cc index b3fc1041f73dd5312d8c57641c7a064c2d6f06a1..25d7e32b819be85377607c5a891b93d903c79f73 100644 --- a/AMDiS/src/LeafData.cc +++ b/AMDiS/src/LeafData.cc @@ -31,7 +31,7 @@ namespace AMDiS { Element* otherChild, int elTypeParent) - { + { bool test = otherChild->deleteElementData(ESTIMATABLE); TEST_EXIT_DBG(test)("couldn't delete LeafDataEstimatable at otherChild"); @@ -60,10 +60,9 @@ namespace AMDiS { { bool test = otherChild->deleteElementData(ESTIMATABLE); TEST_EXIT_DBG(test)("couldn't delete LeafDataEstimatableVec at otherChild"); - parent->setElementData(NEW LeafDataEstimatableVec(parent->getElementData())); ElementData::coarsenElementData(parent, thisChild, otherChild, elTypeParent); - } + } bool LeafDataCoarsenable::refineElementData(Element* parent, Element* child1, @@ -105,7 +104,7 @@ namespace AMDiS { int elTypeParent) { bool test = otherChild->deleteElementData(COARSENABLE); - TEST_EXIT(test)("couldn't delete LeafDataCoarsenableVec at otherChild"); + TEST_EXIT_DBG(test)("couldn't delete LeafDataCoarsenableVec at otherChild"); parent->setElementData(NEW LeafDataCoarsenableVec(parent->getElementData())); ElementData::coarsenElementData(parent, thisChild, otherChild, elTypeParent); } @@ -117,67 +116,59 @@ namespace AMDiS { { ElementData::refineElementData(parent, child1, child2, elType); - int i, j; - Element* child[2] = {child1, child2}; - int dim = parent->getMesh()->getDim(); - LeafDataPeriodic *ld[2] = {NULL, NULL}; - std::list<LeafDataPeriodic::PeriodicInfo>::iterator it; - for(it = periodicInfoList.begin(); - it != periodicInfoList.end(); - ++it) - { - BoundaryType type = it->type; - - int parentSide = it->elementSide; - - int mode = it->periodicMode; - - // for both children - for(i=0; i < 2; i++) { + for (it = periodicInfoList.begin(); + it != periodicInfoList.end(); + ++it) { - // get childs side - int sideOfChild = parent->getSideOfChild(i, parentSide, elType); - if(sideOfChild == -1) continue; + BoundaryType type = it->type; + int parentSide = it->elementSide; + int mode = it->periodicMode; - // create new leaf data if necessary - if(!ld[i]) { - ld[i] = NEW LeafDataPeriodic(child[i]->getElementData()); - child[i]->setElementData(ld[i]); - } + // for both children + for (int i = 0; i < 2; i++) { + // get childs side + int sideOfChild = parent->getSideOfChild(i, parentSide, elType); + if (sideOfChild == -1) + continue; + + // create new leaf data if necessary + if (!ld[i]) { + ld[i] = NEW LeafDataPeriodic(child[i]->getElementData()); + child[i]->setElementData(ld[i]); + } // create new periodic coords - DimVec<WorldVector<double> > coords(dim-1, NO_INIT); - - // for each vertex of childs side - for(j=0; j < dim; j++) { - // get parents vertex nr. - int childVertex = child[i]->getVertexOfPosition(INDEX_OF_DIM(dim-1,dim), - sideOfChild, - j); - int parentVertex = child[i]->getVertexOfParent(i, childVertex, elType); + DimVec<WorldVector<double> > coords(dim - 1, NO_INIT); - if(parentVertex == -1) { - // create coords for new vertex - WorldVector<double> newCoords; - newCoords = (*(it->periodicCoords))[0]; - newCoords += (*(it->periodicCoords))[1]; - newCoords *= 0.5; - coords[j] = newCoords; - } else { - int posAtSide = parent->getPositionOfVertex(parentSide, - parentVertex); - coords[j] = (*(it->periodicCoords))[posAtSide]; - } + // for each vertex of childs side + for (int j = 0; j < dim; j++) { + // get parents vertex nr. + int childVertex = child[i]->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), + sideOfChild, j); + int parentVertex = child[i]->getVertexOfParent(i, childVertex, elType); + + if (parentVertex == -1) { + // create coords for new vertex + WorldVector<double> newCoords; + newCoords = (*(it->periodicCoords))[0]; + newCoords += (*(it->periodicCoords))[1]; + newCoords *= 0.5; + coords[j] = newCoords; + } else { + int posAtSide = parent->getPositionOfVertex(parentSide, + parentVertex); + coords[j] = (*(it->periodicCoords))[posAtSide]; } - - ld[i]->addPeriodicInfo(mode, type, sideOfChild, &coords); } + + ld[i]->addPeriodicInfo(mode, type, sideOfChild, &coords); } + } return false; } @@ -187,14 +178,12 @@ namespace AMDiS { periodicMode = rhs.periodicMode; type = rhs.type; elementSide = rhs.elementSide; - int i, dim; - if(rhs.periodicCoords) { - dim = rhs.periodicCoords->getSize() - 1; + if (rhs.periodicCoords) { + int dim = rhs.periodicCoords->getSize() - 1; periodicCoords = NEW DimVec<WorldVector<double> >(dim, NO_INIT); - for(i = 0; i < dim + 1; i++) { + for (int i = 0; i < dim + 1; i++) { (*periodicCoords)[i] = (*(rhs.periodicCoords))[i]; - } - + } } else { periodicCoords = NULL; } @@ -209,11 +198,10 @@ namespace AMDiS { elementSide(side), periodicCoords(NULL) { - int i, dim; - if(coords) { - dim = coords->getSize() - 1; + if (coords) { + int dim = coords->getSize() - 1; periodicCoords = NEW DimVec<WorldVector<double> >(dim, NO_INIT); - for(i = 0; i < dim + 1; i++) { + for (int i = 0; i < dim + 1; i++) { (*periodicCoords)[i] = (*coords)[i]; } } diff --git a/AMDiS/src/LeafData.h b/AMDiS/src/LeafData.h index 14e8f278f38e4695f5375c7ebd47bdacbb4eba29..77fb546aaeddc607e098b0f1da528fe75a23bff0 100644 --- a/AMDiS/src/LeafData.h +++ b/AMDiS/src/LeafData.h @@ -48,25 +48,28 @@ namespace AMDiS { MEMORY_MANAGED(LeafDataEstimatable); inline bool isOfType(int type) const { - if(type == ESTIMATABLE) + if (type == ESTIMATABLE) { return true; + } + return false; - }; + } class Creator : public CreatorInterface<ElementData> { public: ElementData* create() { return NEW LeafDataEstimatable; - }; + } }; /** \brief * constructor */ LeafDataEstimatable(ElementData *decorated = NULL) - : ElementData(decorated), errorEstimate(0.) - {}; + : ElementData(decorated), + errorEstimate(0.0) + {} /** \brief * Refinement of parent to child1 and child2. @@ -87,12 +90,16 @@ namespace AMDiS { /** \brief * Sets \ref errorEstimate */ - inline void setErrorEstimate(int, double est) { errorEstimate = est; }; + inline void setErrorEstimate(int, double est) { + errorEstimate = est; + } /** \brief * Returns \ref errorEstimate */ - inline double getErrorEstimate(int) { return errorEstimate; }; + inline double getErrorEstimate(int) { + return errorEstimate; + } /** \brief * Implements ElementData::clone(). @@ -108,26 +115,30 @@ namespace AMDiS { // return the clone return newObj; - }; + } /** \brief * Returns the name of element data type. */ - inline std::string getTypeName() const { return "LeafDataEstimatable"; }; + inline std::string getTypeName() const { + return "LeafDataEstimatable"; + } - inline const int getTypeID() const { return ESTIMATABLE; }; + inline const int getTypeID() const { + return ESTIMATABLE; + } void serialize(std::ostream& out) { ElementData::serialize(out); out.write(reinterpret_cast<const char*>(&errorEstimate), sizeof(double)); - }; + } void deserialize(std::istream& in) { ElementData::deserialize(in); in.read(reinterpret_cast<char*>(&errorEstimate), sizeof(double)); - }; + } private: double errorEstimate; @@ -144,14 +155,14 @@ namespace AMDiS { public: ElementData* create() { return NEW LeafDataEstimatableVec; - }; + } }; inline bool isOfType(int type) const { - if(type == ESTIMATABLE) + if (type == ESTIMATABLE) return true; return false; - }; + } /** \brief * constructor @@ -159,8 +170,7 @@ namespace AMDiS { LeafDataEstimatableVec(ElementData *decorated = NULL) : ElementData(decorated) { - //errorEstimate.set(0.0); - }; + } /** \brief * Refinement of parent to child1 and child2. @@ -183,14 +193,14 @@ namespace AMDiS { */ inline void setErrorEstimate(int index, double est) { errorEstimate[index] = est; - }; + } /** \brief * Returns \ref errorEstimate */ inline double getErrorEstimate(int index) { return errorEstimate[index]; - }; + } /** \brief * Implements ElementData::clone(). @@ -207,7 +217,7 @@ namespace AMDiS { // return the clone return newObj; - }; + } void serialize(std::ostream& out) { @@ -219,28 +229,29 @@ namespace AMDiS { out.write(reinterpret_cast<const char*>(&(it->first)), sizeof(int)); out.write(reinterpret_cast<const char*>(&(it->second)), sizeof(double)); } - }; + } void deserialize(std::istream& in) { ElementData::deserialize(in); - unsigned int i, size; + unsigned size; in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int)); - for(i = 0; i < size; i++) { + for(unsigned int i = 0; i < size; i++) { int index; double estimate; in.read(reinterpret_cast<char*>(&index), sizeof(int)); in.read(reinterpret_cast<char*>(&estimate), sizeof(double)); errorEstimate[index] = estimate; } - }; + } - std::string getTypeName() const - { + std::string getTypeName() const { return "LeafDataEstimatableVec"; - }; + } - inline const int getTypeID() const { return ESTIMATABLE; }; + inline const int getTypeID() const { + return ESTIMATABLE; + } private: std::map<int, double> errorEstimate; @@ -249,12 +260,12 @@ namespace AMDiS { class LeafDataCoarsenableInterface { public: - virtual ~LeafDataCoarsenableInterface() {}; + virtual ~LeafDataCoarsenableInterface() {} /** \brief * Sets \ref coarseningError */ - virtual void setCoarseningErrorEstimate(int index, double est)=0; + virtual void setCoarseningErrorEstimate(int index, double est) = 0; /** \brief * Returns \ref coarseningError @@ -271,23 +282,28 @@ namespace AMDiS { class Creator : public CreatorInterface<ElementData> { public: - ElementData* create() { + ElementData* create() + { return NEW LeafDataCoarsenable; - }; + } }; inline bool isOfType(int type) const { if(type == COARSENABLE) return true; return false; - }; + } /** \brief * constructor */ LeafDataCoarsenable(ElementData *decorated = NULL) - : ElementData(decorated), coarseningError(0.0) - {}; + : ElementData(decorated), + coarseningError(0.00) + {} + + ~LeafDataCoarsenable() + {} /** \brief * Refinement of parent to child1 and child2. @@ -317,21 +333,21 @@ namespace AMDiS { // return the clone return newObj; - }; + } /** \brief * Sets \ref coarseningError */ virtual void setCoarseningErrorEstimate(int , double est) { coarseningError = est; - }; + } /** \brief * Returns \ref coarseningError */ virtual double getCoarseningErrorEstimate(int) { return coarseningError; - }; + } // ===== Serializable implementation ===== @@ -339,20 +355,21 @@ namespace AMDiS { { ElementData::serialize(out); out.write(reinterpret_cast<const char*>(&coarseningError), sizeof(double)); - }; + } void deserialize(std::istream& in) { ElementData::deserialize(in); in.read(reinterpret_cast<char*>(&coarseningError), sizeof(double)); - }; + } - std::string getTypeName() const - { + std::string getTypeName() const { return "LeafDataCoarsenable"; - }; + } - inline const int getTypeID() const { return COARSENABLE; }; + inline const int getTypeID() const { + return COARSENABLE; + } private: double coarseningError; /**< \brief coarsening error */ diff --git a/AMDiS/src/MacroElement.cc b/AMDiS/src/MacroElement.cc index d970805cbb66693c4a9decced7a0f0284eee18dc..561189177cc2c0d1904b000cbabe37ea3546417e 100644 --- a/AMDiS/src/MacroElement.cc +++ b/AMDiS/src/MacroElement.cc @@ -26,8 +26,9 @@ namespace AMDiS { MacroElement::~MacroElement() { - if (element) + if (element) { DELETE element; + } } MacroElement& MacroElement::operator=(const MacroElement &el) diff --git a/AMDiS/src/MacroElement.h b/AMDiS/src/MacroElement.h index dbbcc72d005fd4362e71fe3755446a6cdea0a283..584ed99680d02ca61f685831cbd1a621da0abf45 100644 --- a/AMDiS/src/MacroElement.h +++ b/AMDiS/src/MacroElement.h @@ -199,7 +199,7 @@ namespace AMDiS { boundary[i] = b; } - inline void setCoord(int i, const WorldVector<double> &c) { + inline void setCoord(int i, const WorldVector<double> c) { coord[i] = c; } diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc index 230c947f5188c11dfed0a3b01a1ffb2a423014a2..996f6b88e953071b6b19f11c3f2b641c669dc4e2 100644 --- a/AMDiS/src/MacroReader.cc +++ b/AMDiS/src/MacroReader.cc @@ -226,7 +226,7 @@ namespace AMDiS { for (int k = 0; k < mesh->getGeo(VERTEX); k++) { (*(mel + i))->setCoord(k, coords[melVertex[i][k]]); - const_cast<Element*>((*(mel+i))->getElement())-> + const_cast<Element*>((*(mel + i))->getElement())-> setDOF(k, dof[melVertex[i][k]]); } } @@ -245,8 +245,8 @@ namespace AMDiS { MacroElement *neigh = const_cast<MacroElement*>(mel[i]->getNeighbour(k)); if (neigh) { - int j; - for (j = 0; j < mesh->getGeo(NEIGH); j++) + int j = 0; + for (; j < mesh->getGeo(NEIGH); j++) if (neigh->getNeighbour(j) == *(mel + i)) break; @@ -298,7 +298,7 @@ namespace AMDiS { if (mesh->getNumberOfDOFs(CENTER)) { for (int i = 0; i < mesh->getNumberOfMacros(); i++) { const_cast<Element*>(mel[i]->getElement())-> - setDOF(mesh->getNode(CENTER),mesh->getDOF(CENTER)); + setDOF(mesh->getNode(CENTER), mesh->getDOF(CENTER)); } } @@ -336,7 +336,7 @@ namespace AMDiS { } } - return(macroInfo); + return macroInfo; } /****************************************************************************/ @@ -363,15 +363,16 @@ namespace AMDiS { } dof = GET_MEMORY(DegreeOfFreedom*, nv); - coords = NEW WorldVector<double>[nv]; + coords = GET_MEMORY(WorldVector<double>, nv); mel_vertex = GET_MEMORY(int*, ne); for (int i = 0; i < ne; i++) { - mel_vertex[i]=GET_MEMORY(int, mesh->getGeo(VERTEX)); + mel_vertex[i] = GET_MEMORY(int, mesh->getGeo(VERTEX)); } - for (int i = 0; i < nv; i++) + for (int i = 0; i < nv; i++) { dof[i] = mesh->getDOF(VERTEX); + } for (int i = 0; i < ne; i++) { mel[i]->element = mesh->createNewElement(); @@ -391,8 +392,11 @@ namespace AMDiS { FREE_MEMORY(mel_vertex[i], int, mesh->getGeo(VERTEX)); FREE_MEMORY(mel_vertex, int*, ne); + FREE_MEMORY(coords, WorldVector<double>, nv); + coords = NULL; + FREE_MEMORY(dof, DegreeOfFreedom*, nv); + dof = NULL; - coords = NULL; /* must no be freed!!! still used in the mesh!!! */ mesh = NULL; neigh_set = false; } @@ -474,14 +478,15 @@ namespace AMDiS { void MacroInfo::readAMDiSMacro(const char *filename, Mesh* mesh) { - FUNCNAME("MacroInfo::readAMDiSMacro"); - FILE *file; - int dim; - int dow, nv, ne, j, k; - double dbl; - char name[128], line[256]; - int line_no, n_keys, i_key, sort_key[N_KEYS], nv_key, ne_key; - int key_def[N_KEYS] = {0,0,0,0,0,0,0,0,0,0,0,0}; + FUNCNAME("MacroInfo::readAMDiSMacro()"); + + FILE *file; + int dim; + int dow, nv, ne, j, k; + double dbl; + char name[128], line[256]; + int line_no, n_keys, i_key, sort_key[N_KEYS], nv_key, ne_key; + int key_def[N_KEYS] = {0,0,0,0,0,0,0,0,0,0,0,0}; const char *key; DimVec<int> *ind = NULL; @@ -791,6 +796,10 @@ namespace AMDiS { } } + if (ind) { + DELETE ind; + } + fclose(file); } diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index 8d9e8c70067e0a4c3d6050efd8ff68d00d2b6acc..5af9d8dfad381f12a9c39abb839c7fbd6b9e52c3 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -120,7 +120,7 @@ namespace AMDiS { Mesh::~Mesh() { - serializedDOFs.empty(); + Element::deletedDOFs.clear(); for (std::deque<MacroElement*>::const_iterator it = macroElements.begin(); it != macroElements.end(); @@ -129,11 +129,21 @@ namespace AMDiS { DELETE *it; } - serializedDOFs.empty(); + Element::deletedDOFs.clear(); - if (macroFileInfo) { + if (macroFileInfo != NULL) { DELETE macroFileInfo; } + if (elementPrototype) { + DELETE elementPrototype; + } + if (elementDataPrototype) { + DELETE elementDataPrototype; + } + + for (int i = 0; i < static_cast<int>(admin.size()); i++) { + DELETE admin[i]; + } } Mesh& Mesh::operator=(const Mesh& m) @@ -173,9 +183,9 @@ namespace AMDiS { admin[i]->setMesh(this); } + /* ====================== Copy macro elements =================== */ - macroElements.clear(); - + // mapIndex[i] is the index of the MacroElement element in the vector // macroElements, for which holds: element->getIndex() = i std::map<int, int> mapIndex; @@ -392,7 +402,7 @@ namespace AMDiS { int ndof = getNumberOfDOFs(position); if (ndof <= 0) - return(NULL); + return NULL; DegreeOfFreedom *dof = GET_MEMORY(DegreeOfFreedom, ndof); @@ -410,7 +420,7 @@ namespace AMDiS { } } - return(dof); + return dof; } @@ -419,13 +429,13 @@ namespace AMDiS { FUNCNAME("Mesh::createDOFPtrs()"); if (nNodeEl <= 0) - return(NULL); + return NULL; DegreeOfFreedom **ptrs = GET_MEMORY(DegreeOfFreedom*, nNodeEl); for (int i = 0; i < nNodeEl; i++) ptrs[i] = NULL; - return(ptrs); + return ptrs; } void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs) @@ -452,7 +462,7 @@ namespace AMDiS { addDOFAdmin(localAdmin); - return(localAdmin); + return localAdmin; } @@ -469,7 +479,7 @@ namespace AMDiS { } } - return(localAdmin); + return localAdmin; } void Mesh::freeDOF(DegreeOfFreedom* dof, GeoIndex position) @@ -598,7 +608,7 @@ namespace AMDiS { DELETE mel_info; - return(inside); + return inside; } bool Mesh::findElementAtPoint(const WorldVector<double>& xy, @@ -615,7 +625,7 @@ namespace AMDiS { DELETE el_info; - return(val); + return val; } @@ -638,7 +648,7 @@ namespace AMDiS { final_lambda[i] = lambda[i]; } - return(true); + return true; } else { /* outside */ if (g_xy0) { /* find boundary point of [xy0, xy] */ el_info->worldToCoord(*(g_xy0), &c_lambda); @@ -652,9 +662,9 @@ namespace AMDiS { if (dim == 3) MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s); - return(false); /* ??? */ + return false; /* ??? */ } else { - return(false); + return false; } } } @@ -770,7 +780,7 @@ namespace AMDiS { final_el_info); DELETE c_el_info; - return(inside); + return inside; } @@ -1041,13 +1051,13 @@ namespace AMDiS { macroFileInfo->clear(getNumberOfEdges(), getNumberOfVertices()); DELETE macroFileInfo; + macroFileInfo = NULL; } int Mesh::calcMemoryUsage() { - int result = 0; + int result = sizeof(Mesh); - result += sizeof(Mesh); for (int i = 0; i < static_cast<int>(macroElements.size()); i++) { result += macroElements[i]->calcMemoryUsage(); } diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index 8830c7c9c66a495c7a3d35bf3b3063ee229200f6..449e1365a02b73960d3eaa6eb217973cbf18e5e6 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -863,7 +863,7 @@ namespace AMDiS { /** \brief * List of all MacroElements of this Mesh */ - std::deque<MacroElement*> macroElements; + std::deque<MacroElement*> macroElements; /** \brief * Needed during DOF compression (\ref DOFAdmin::compress). diff --git a/AMDiS/src/MultiGridSolverBase.h b/AMDiS/src/MultiGridSolverBase.h index 4a79297f91eb94ff0a5fe0a2b79166a652074892..8480db2f2441e69d8e5e7ddd0ef8a98cfb345145 100755 --- a/AMDiS/src/MultiGridSolverBase.h +++ b/AMDiS/src/MultiGridSolverBase.h @@ -22,6 +22,8 @@ #ifndef AMDIS_MULTIGRIDSOLVERBASE_H #define AMDIS_MULTIGRIDSOLVERBASE_H +#include "InterpolRestrictMatrix.h" +#include "LagrangeInterpolRestrict.h" #include <string> namespace AMDiS { @@ -33,12 +35,12 @@ namespace AMDiS { /** \brief * Constructor, reads all multigrid parameters from the init file. */ - MultiGridSolverBase(const std::string &name) ; + MultiGridSolverBase(const std::string &name); /** \brief * Destructor */ - virtual ~MultiGridSolverBase() {}; + virtual ~MultiGridSolverBase() {} /** \brief * Abstract function for solving an equation using multigrid method. @@ -145,6 +147,18 @@ namespace AMDiS { * zero value. */ int fmg_; + + /** \brief + * interpolMatrices[n][m] stores the interpolation matrix for + * n-dimensional elements and m-dimensional basis functions. + */ + InterpolRestrictMatrix *interpolMatrices[3][4]; + + /** \brief + * interpolMatrices[n][m] stores the restriction matrix for + * n-dimensional elements and m-dimensional basis functions. + */ + InterpolRestrictMatrix *restrictMatrices[3][4]; }; } diff --git a/AMDiS/src/MultiGridSolverBase.hh b/AMDiS/src/MultiGridSolverBase.hh index e5254920bbc8d50febd83eaf31ced1dfa5a3bfb9..c4ff0c2f57cb052c426b390beb23a634cd2f6c1c 100644 --- a/AMDiS/src/MultiGridSolverBase.hh +++ b/AMDiS/src/MultiGridSolverBase.hh @@ -22,6 +22,36 @@ namespace AMDiS { GET_PARAMETER(0, name_ + "->post smoothing steps", "%d", &postSmoothingSteps_); GET_PARAMETER(0, name_ + "->min level", "%d", &minLevel_); GET_PARAMETER(0, name_ + "->full mg", "%d", &fmg_); + + interpolMatrices[0][0] = new IM_1d_1; + interpolMatrices[0][1] = new IM_1d_2; + interpolMatrices[0][2] = NULL; + interpolMatrices[0][3] = NULL; + + interpolMatrices[1][0] = new IM_2d_1; + interpolMatrices[1][1] = new IM_2d_2; + interpolMatrices[1][2] = new IM_2d_3; + interpolMatrices[1][3] = new IM_2d_4; + + interpolMatrices[2][0] = new IM_3d_1; + interpolMatrices[2][1] = new IM_3d_2; + interpolMatrices[2][2] = new IM_3d_3; + interpolMatrices[2][3] = new IM_3d_4; + + restrictMatrices[0][0] = new InterpolRestrictMatrix(*(interpolMatrices[0][0]), true); + restrictMatrices[0][1] = new InterpolRestrictMatrix(*(interpolMatrices[0][1]), true); + restrictMatrices[0][2] = NULL; + restrictMatrices[0][3] = NULL; + + restrictMatrices[1][0] = new InterpolRestrictMatrix(*(interpolMatrices[1][0]), true); + restrictMatrices[1][1] = new InterpolRestrictMatrix(*(interpolMatrices[1][1]), true); + restrictMatrices[1][2] = new InterpolRestrictMatrix(*(interpolMatrices[1][2]), true); + restrictMatrices[1][3] = new InterpolRestrictMatrix(*(interpolMatrices[1][3]), true); + + restrictMatrices[2][0] = new InterpolRestrictMatrix(*(interpolMatrices[2][0]), true); + restrictMatrices[2][1] = new InterpolRestrictMatrix(*(interpolMatrices[2][1]), true); + restrictMatrices[2][2] = new InterpolRestrictMatrix(*(interpolMatrices[2][2]), true); + restrictMatrices[2][3] = new InterpolRestrictMatrix(*(interpolMatrices[2][3]), true); } template<typename MatrixType, typename VectorType> diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h index 19a158fd86e8b006233c6f8c9c08526e937c5ac8..3830691c3015fcb2a955f158e70bb88eba427d99 100644 --- a/AMDiS/src/OEMSolver.h +++ b/AMDiS/src/OEMSolver.h @@ -240,7 +240,7 @@ namespace AMDiS { /** \brief * Prints solve information. Uses \ref SOLVE_INFO. */ - int solve_info(const char *, int, double, double *); + int solve_info(const char *, int, double, double *); protected: /** \brief diff --git a/AMDiS/src/OEMSolver.hh b/AMDiS/src/OEMSolver.hh index a0a57c26c89091086cb1e756f271c56289f69613..aa2d25f0ae8fde13c303b551b577f892b606dd61 100644 --- a/AMDiS/src/OEMSolver.hh +++ b/AMDiS/src/OEMSolver.hh @@ -37,6 +37,9 @@ namespace AMDiS { OEMSolver<VectorType>::~OEMSolver() { FUNCNAME("OEMSolver::~OEMSolver"); + + if (vectorCreator) + DELETE vectorCreator; } diff --git a/AMDiS/src/Parameters.cc b/AMDiS/src/Parameters.cc index 4eeacde6a0351aa72987ce097cca2a9c2bdf76ef..e47f474b9795f44c660ce01764c94cfc028b6c76 100644 --- a/AMDiS/src/Parameters.cc +++ b/AMDiS/src/Parameters.cc @@ -456,7 +456,7 @@ namespace AMDiS { void Parameters::initIntern() { - if (NULL == singlett) { + if (singlett == NULL) { singlett = NEW Parameters; } } @@ -762,6 +762,13 @@ namespace AMDiS { fp.close(); } + void Parameters::clear() + { + if (singlett != NULL) { + DELETE singlett; + } + } + int Parameters::param::operator==(const param& aParam)const{ return key == aParam.key; } diff --git a/AMDiS/src/Parameters.h b/AMDiS/src/Parameters.h index 7e7c354f1af346533de2af10bf352929ebbe74ef..a564d3b857b5b6ab3dfd4593605a31b1ac30ce0d 100644 --- a/AMDiS/src/Parameters.h +++ b/AMDiS/src/Parameters.h @@ -208,12 +208,16 @@ namespace AMDiS { /** \brief * Returns specified info level */ - static int getMsgInfo() { return (singlett)?singlett->msgInfo:0; }; + static int getMsgInfo() { + return (singlett) ? singlett->msgInfo : 0; + } /** \brief * Returns specified wait value */ - static int getMsgWait() { return (singlett)?singlett->msgWait:0; }; + static int getMsgWait() { + return (singlett) ? singlett->msgWait : 0; + } /** \brief * Writes all up to now initialized parameters to file according to the @@ -228,11 +232,13 @@ namespace AMDiS { */ static bool initialized() { return (singlett != NULL); - }; + } static Parameters *getSingleton() { return singlett; - }; + } + + static void clear(); // ===== Serializable implementation ===== @@ -245,7 +251,7 @@ namespace AMDiS { for(i = 0; i < size; i++) { allParam[i].serialize(out); } - }; + } void deserialize(std::istream &in) { in.read(reinterpret_cast<char*>(¶mInfo), sizeof(int)); @@ -257,7 +263,7 @@ namespace AMDiS { for(i = 0; i < size; i++) { allParam[i].deserialize(in); } - }; + } private: static Parameters * singlett; @@ -275,9 +281,9 @@ namespace AMDiS { class param : public Serializable { public: - param() {}; + param() {} - virtual ~param() {}; + virtual ~param() {} param(const std::string& nkey, const std::string& nparameters, @@ -289,15 +295,17 @@ namespace AMDiS { filename(nfilename), funcName(nfuncName), lineNo(line) - {}; + {} - param(const std::string& k): key(k) {}; + param(const std::string& k) + : key(k) + {} int operator==(const class param& aParam) const; int operator!=(const class param& aParam) const { return !(aParam==*this); - }; + } // ===== Serializable implementation ===== @@ -307,7 +315,7 @@ namespace AMDiS { out << filename << std::endl; out << funcName << std::endl; out.write(reinterpret_cast<const char*>(&lineNo), sizeof(int)); - }; + } void deserialize(std::istream &in) { in >> key; in.get(); @@ -315,21 +323,24 @@ namespace AMDiS { in >> filename; in.get(); in >> funcName; in.get(); in.read(reinterpret_cast<char*>(&lineNo), sizeof(int)); - }; + } public: std::string key; std::string parameters; std::string filename; std::string funcName; - int lineNo; + int lineNo; }; static const char comment; std::string buffer; - inline bool isBlankChar(const char s) {return (s==' '||s=='\t'||s =='\f'||s=='\r'); }; + inline bool isBlankChar(const char s) { + return (s == ' ' || s == '\t' || s == '\f' || s == '\r'); + } + const char *getNextWord(std::string *s) const; void read(const std::string& filename,const std::string& maxKey=""); @@ -355,9 +366,9 @@ namespace AMDiS { : paramInfo(1), msgInfo(1), msgWait(1) - {}; + {} - virtual ~Parameters() {}; + virtual ~Parameters() {} static void initIntern(); diff --git a/AMDiS/src/ProblemInterpolScal.cc b/AMDiS/src/ProblemInterpolScal.cc index e419e3f6c5c9c1123a28b6067d7ceb1d6600ed3e..f0c8d9b104dc0ba438d20b43d54cb027f740beda 100644 --- a/AMDiS/src/ProblemInterpolScal.cc +++ b/AMDiS/src/ProblemInterpolScal.cc @@ -25,13 +25,14 @@ namespace AMDiS { } void ProblemInterpolScal::solve(AdaptInfo *adaptInfo) { - mesh_->dofCompress(); - solution_->interpol(interpolFct_); + mesh->dofCompress(); + solution->interpol(interpolFct_); } void ProblemInterpolScal::estimate(AdaptInfo *adaptInfo, double) { FUNCNAME("ProblemIterpolScal::estimate()"); - double errMax = 0.0, errSum = 0.0; + double errMax = 0.0; + double errSum = 0.0; int errorNorm = 0; int relErr = 0; @@ -45,10 +46,10 @@ namespace AMDiS { switch (errorNorm) { case 1: - errSum = Error<double>::H1Err(*grdInterpolFct_, *solution_, relErr, &errMax, true); + errSum = Error<double>::H1Err(*grdInterpolFct_, *solution, relErr, &errMax, true); break; case 2: - errSum = Error<double>::L2Err(*interpolFct_, *solution_, relErr, &errMax, true); + errSum = Error<double>::L2Err(*interpolFct_, *solution, relErr, &errMax, true); break; default: ERROR_EXIT("invalid error norm\n"); diff --git a/AMDiS/src/ProblemNonLin.cc b/AMDiS/src/ProblemNonLin.cc index 6a8309da4f3d020edf51a85060433bf06a0be918..171abe6857d6ce31dc931b9b89d7da9b8fe77915 100644 --- a/AMDiS/src/ProblemNonLin.cc +++ b/AMDiS/src/ProblemNonLin.cc @@ -28,7 +28,7 @@ namespace AMDiS { ((!initFlag.isSet(INIT_UPDATER))&& (adoptFlag.isSet(INIT_NONLIN_SOLVER))))) { - TEST_EXIT(updater_==NULL)("updater already created\n"); + TEST_EXIT(updater_ == NULL)("updater already created\n"); updater_ = dynamic_cast<ProblemNonLinScal*>(adoptProblem)->getUpdater(); } } @@ -63,7 +63,7 @@ namespace AMDiS { NonLinSolverCreator<DOFVector<double> > *nonLinSolverCreator = dynamic_cast<NonLinSolverCreator<DOFVector<double> >*>( CreatorMap<NonLinSolver<DOFVector<double> > >::getCreator(nonLinSolverType)); - nonLinSolverCreator->setLinearSolver(solver_); + nonLinSolverCreator->setLinearSolver(solver); nonLinSolverCreator->setName(name_ + "->nonlin solver"); nonLinSolverCreator->setNonLinUpdater(updater_); nonLinSolver_ = nonLinSolverCreator->create(); @@ -72,7 +72,7 @@ namespace AMDiS { void ProblemNonLinScal::solve(AdaptInfo *adaptInfo) { TEST_EXIT(nonLinSolver_)("no non-linear solver!\n"); - int iter = nonLinSolver_->solve(matVec_, solution_, rhs_, leftPrecon_, rightPrecon_); + int iter = nonLinSolver_->solve(matVec_, solution, rhs, leftPrecon, rightPrecon); adaptInfo->setSolverIterations(iter); } @@ -82,16 +82,16 @@ namespace AMDiS { MSG("%d DOFs for %s\n", feSpace_->getAdmin()->getUsedSize(), feSpace_->getName().c_str()); TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh_, -1, + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA); // for all elements ... - while(elInfo) { - if(solution_->getBoundaryManager()) - solution_->getBoundaryManager()->fillBoundaryConditions(elInfo, solution_); + while (elInfo) { + if (solution->getBoundaryManager()) + solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution); elInfo = stack.traverseNext(elInfo); } diff --git a/AMDiS/src/ProblemNonLin.h b/AMDiS/src/ProblemNonLin.h index 19e7538775582a36fa281fea3624c3b520615a44..045d6013233ff7125fab8aa85e144d2a9360b578 100644 --- a/AMDiS/src/ProblemNonLin.h +++ b/AMDiS/src/ProblemNonLin.h @@ -73,8 +73,8 @@ namespace AMDiS { * Used in \ref initialize(). */ virtual void createUpdater() { - updater_ = NEW NonLinUpdaterScal(systemMatrix_); - }; + updater_ = NEW NonLinUpdaterScal(systemMatrix); + } /** \brief * Used in \ref initialize(). @@ -100,8 +100,8 @@ namespace AMDiS { */ inline void setU0(AbstractFunction<double, WorldVector<double> > *u0Fct) { u0_ = u0Fct; - solution_->interpol(u0_); - }; + solution->interpol(u0_); + } /** \brief * Sets \ref nonLinSolver_. @@ -114,7 +114,7 @@ namespace AMDiS { * Sets \ref updater. */ inline void setUpdater(NonLinUpdaterScal *updater) { - updater_=updater; + updater_ = updater; } /** \brief @@ -169,7 +169,7 @@ namespace AMDiS { { u0_.resize(nComponents); u0_.set(NULL); - }; + } /** \brief * Sets \ref u0_ and interpolates it to \ref solution_. @@ -181,7 +181,7 @@ namespace AMDiS { TEST_EXIT(index < nComponents)("invalid index\n"); u0_[index] = u0Fct; solution_->getDOFVector(index)->interpol(u0Fct); - }; + } /** \brief * Destructor. @@ -200,7 +200,7 @@ namespace AMDiS { */ virtual void createUpdater() { updater_ = NEW NonLinUpdaterVec(systemMatrix_); - }; + } /** \brief * Used in \ref initialize(). @@ -237,14 +237,14 @@ namespace AMDiS { * Sets \ref nonLinSolver_. */ inline void setNonLinSolver(NonLinSolver<SystemVector> *nonLinSolver) { - nonLinSolver_=nonLinSolver; + nonLinSolver_ = nonLinSolver; } /** \brief * Sets \ref updater. */ inline void setUpdater(NonLinUpdaterVec *updater) { - updater_=updater; + updater_ = updater; } diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc index bed26d93c19c5e479c28e6ac8b156b2c470367c3..68cb2b3e68dcd8baf06f775e5e2c0ad2596ae6c3 100644 --- a/AMDiS/src/ProblemScal.cc +++ b/AMDiS/src/ProblemScal.cc @@ -29,53 +29,84 @@ namespace AMDiS { ProblemScal *ProblemScal::traversePtr_ = NULL; + ProblemScal::~ProblemScal() + { + FUNCNAME("ProblemScal::~ProblemScal()"); + + for (int i = 0; i < static_cast<int>(fileWriters.size()); i++) { + DELETE fileWriters[i]; + } + + if (systemMatrix) + DELETE systemMatrix; + if (rhs) + DELETE rhs; + if (solution) + DELETE solution; + if (estimator) + DELETE estimator; + if (marker) + DELETE marker; + if (solver) + DELETE solver; + if (leftPrecon) + DELETE leftPrecon; + if (rightPrecon) + DELETE rightPrecon; + if (mesh) + DELETE mesh; + + FiniteElemSpace::clear(); + Lagrange::clear(); + } + void ProblemScal::writeFiles(AdaptInfo *adaptInfo, bool force) { - for (int i = 0; i < static_cast<int>(fileWriters_.size()); i++) { - fileWriters_[i]->writeFiles(adaptInfo, force); + for (int i = 0; i < static_cast<int>(fileWriters.size()); i++) { + fileWriters[i]->writeFiles(adaptInfo, force); } } void ProblemScal::interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct) { - solution_->interpol(fct); + solution->interpol(fct); } void ProblemScal::addMatrixOperator(Operator *op, double *factor, double *estFactor) { - systemMatrix_->addOperator(op, factor, estFactor); + systemMatrix->addOperator(op, factor, estFactor); } void ProblemScal::addVectorOperator(Operator *op, double *factor, double *estFactor) { - rhs_->addOperator(op, factor, estFactor); + rhs->addOperator(op, factor, estFactor); } void ProblemScal::addDirichletBC(BoundaryType type, AbstractFunction<double, WorldVector<double> >* b) { DirichletBC *dirichlet = new DirichletBC(type, b, feSpace_); - if (systemMatrix_) - systemMatrix_->getBoundaryManager()->addBoundaryCondition(dirichlet); - if (rhs_) - rhs_->getBoundaryManager()->addBoundaryCondition(dirichlet); - if (solution_) - solution_->getBoundaryManager()->addBoundaryCondition(dirichlet); + if (systemMatrix) + systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet); + if (rhs) + rhs->getBoundaryManager()->addBoundaryCondition(dirichlet); + if (solution) + solution->getBoundaryManager()->addBoundaryCondition(dirichlet); } void ProblemScal::addDirichletBC(BoundaryType type, DOFVector<double> *vec) { DirichletBC *dirichlet = new DirichletBC(type, vec); - if (systemMatrix_) - systemMatrix_->getBoundaryManager()->addBoundaryCondition(dirichlet); - if (rhs_) - rhs_->getBoundaryManager()->addBoundaryCondition(dirichlet); - if (solution_) - solution_->getBoundaryManager()->addBoundaryCondition(dirichlet); + if (systemMatrix) + systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet); + if (rhs) + rhs->getBoundaryManager()->addBoundaryCondition(dirichlet); + if (solution) + solution->getBoundaryManager()->addBoundaryCondition(dirichlet); } void ProblemScal::addRobinBC(BoundaryType type, @@ -83,18 +114,18 @@ namespace AMDiS { AbstractFunction<double, WorldVector<double> > *r) { RobinBC *robin = new RobinBC(type, n, r, feSpace_); - if (rhs_) - rhs_->getBoundaryManager()->addBoundaryCondition(robin); - if (systemMatrix_) - systemMatrix_->getBoundaryManager()->addBoundaryCondition(robin); + if (rhs) + rhs->getBoundaryManager()->addBoundaryCondition(robin); + if (systemMatrix) + systemMatrix->getBoundaryManager()->addBoundaryCondition(robin); } void ProblemScal::addNeumannBC(BoundaryType type, AbstractFunction<double, WorldVector<double> > *n) { NeumannBC *neumann = new NeumannBC(type, n, feSpace_); - if(rhs_) - rhs_->getBoundaryManager()->addBoundaryCondition(neumann); + if(rhs) + rhs->getBoundaryManager()->addBoundaryCondition(neumann); } void ProblemScal::addRobinBC(BoundaryType type, @@ -102,20 +133,20 @@ namespace AMDiS { DOFVector<double> *r) { RobinBC *robin = new RobinBC(type, n, r, feSpace_); - if (rhs_) - rhs_->getBoundaryManager()->addBoundaryCondition(robin); - if (systemMatrix_) - systemMatrix_->getBoundaryManager()->addBoundaryCondition(robin); + if (rhs) + rhs->getBoundaryManager()->addBoundaryCondition(robin); + if (systemMatrix) + systemMatrix->getBoundaryManager()->addBoundaryCondition(robin); } void ProblemScal::addPeriodicBC(BoundaryType type) { PeriodicBC *periodic = new PeriodicBC(type, feSpace_); - if (systemMatrix_) - systemMatrix_->getBoundaryManager()->addBoundaryCondition(periodic); - if (rhs_) - rhs_->getBoundaryManager()->addBoundaryCondition(periodic); + if (systemMatrix) + systemMatrix->getBoundaryManager()->addBoundaryCondition(periodic); + if (rhs) + rhs->getBoundaryManager()->addBoundaryCondition(periodic); } void ProblemScal::createMesh() @@ -136,7 +167,7 @@ namespace AMDiS { TEST_EXIT(dim)("no problem dimension specified!\n"); // create the mesh - mesh_ = NEW Mesh(meshName, dim); + mesh = NEW Mesh(meshName, dim); switch (dim) { case 1: @@ -159,7 +190,7 @@ namespace AMDiS { void ProblemScal::setMeshFromProblemVec(ProblemVec* pv, int i) { - mesh_ = pv->getMesh(i); + mesh = pv->getMesh(i); coarseningManager_ = pv->getCoarseningManager(i); refinementManager_ = pv->getRefinementManager(i); } @@ -167,8 +198,8 @@ namespace AMDiS { Flag ProblemScal::markElements(AdaptInfo *adaptInfo) { - if (marker_) - return marker_->markMesh(adaptInfo, mesh_); + if (marker) + return marker->markMesh(adaptInfo, mesh); else WARNING("no marker\n"); @@ -178,14 +209,14 @@ namespace AMDiS { Flag ProblemScal::refineMesh(AdaptInfo *adaptInfo) { - return refinementManager_->refineMesh(mesh_); + return refinementManager_->refineMesh(mesh); } Flag ProblemScal::coarsenMesh(AdaptInfo *adaptInfo) { if (adaptInfo->isCoarseningAllowed(0)) { - return coarseningManager_->coarsenMesh(mesh_); + return coarseningManager_->coarsenMesh(mesh); } else { WARNING("coarsening not allowed\n"); return 0; @@ -195,7 +226,7 @@ namespace AMDiS { void ProblemScal::solve(AdaptInfo *adaptInfo, bool fixedMatrix) { FUNCNAME("Problem::solve()"); - if (!solver_) { + if (!solver) { WARNING("no solver\n"); return; } @@ -205,7 +236,7 @@ namespace AMDiS { #endif clock_t first = clock(); - int iter = solver_->solve(matVec_, solution_, rhs_, leftPrecon_, rightPrecon_); + int iter = solver->solve(matVec_, solution, rhs, leftPrecon, rightPrecon); #ifdef _OPENMP INFO(info_, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n", @@ -217,9 +248,9 @@ namespace AMDiS { #endif adaptInfo->setSolverIterations(iter); - adaptInfo->setMaxSolverIterations(solver_->getMaxIterations()); - adaptInfo->setSolverTolerance(solver_->getTolerance()); - adaptInfo->setSolverResidual(solver_->getResidual()); + adaptInfo->setMaxSolverIterations(solver->getMaxIterations()); + adaptInfo->setSolverTolerance(solver->getTolerance()); + adaptInfo->setSolverResidual(solver->getResidual()); } void ProblemScal::initialize(Flag initFlag, @@ -229,7 +260,7 @@ namespace AMDiS { FUNCNAME("Problem::initialize()"); // === create mesh === - if (mesh_) { + if (mesh) { WARNING("mesh already created\n"); } else { if (initFlag.isSet(CREATE_MESH) || @@ -242,14 +273,14 @@ namespace AMDiS { (adoptFlag.isSet(INIT_MESH) || adoptFlag.isSet(INIT_SYSTEM) || adoptFlag.isSet(INIT_FE_SPACE))) { - TEST_EXIT(!mesh_)("mesh already created\n"); - mesh_ = adoptProblem->getMesh(); + TEST_EXIT(!mesh)("mesh already created\n"); + mesh = adoptProblem->getMesh(); refinementManager_ = adoptProblem->refinementManager_; coarseningManager_ = adoptProblem->coarseningManager_; } } - if (!mesh_) + if (!mesh) WARNING("no mesh created\n"); // === create fespace === @@ -273,60 +304,60 @@ namespace AMDiS { createMatricesAndVectors(); } if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) { - TEST_EXIT(!solution_)("solution already created\n"); - TEST_EXIT(!rhs_)("rhs already created\n"); - TEST_EXIT(!systemMatrix_)("systemMatrix already created\n"); - solution_ = adoptProblem->getSolution(); - rhs_ = adoptProblem->getRHS(); - systemMatrix_ = adoptProblem->getSystemMatrix(); + TEST_EXIT(!solution)("solution already created\n"); + TEST_EXIT(!rhs)("rhs already created\n"); + TEST_EXIT(!systemMatrix)("systemMatrix already created\n"); + solution = adoptProblem->getSolution(); + rhs = adoptProblem->getRHS(); + systemMatrix = adoptProblem->getSystemMatrix(); } // === create solver === - if (solver_) { + if (solver) { WARNING("solver already created\n"); } else { if (initFlag.isSet(INIT_SOLVER)) { createSolver(); } if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) { - TEST_EXIT(!solver_)("solver already created\n"); - solver_ = adoptProblem->getSolver(); + TEST_EXIT(!solver)("solver already created\n"); + solver = adoptProblem->getSolver(); } } - if (!solver_) + if (!solver) WARNING("no solver created\n"); // === create estimator === - if (estimator_) { + if (estimator) { WARNING("estimator already created\n"); } else { if (initFlag.isSet(INIT_ESTIMATOR)) { createEstimator(); } if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) { - TEST_EXIT(!estimator_)("estimator already created\n"); - estimator_ = adoptProblem->getEstimator(); + TEST_EXIT(!estimator)("estimator already created\n"); + estimator = adoptProblem->getEstimator(); } } - if (!estimator_) + if (!estimator) WARNING("no estimator created\n"); // === create marker === - if (marker_) { + if (marker) { WARNING("marker already created\n"); } else { if (initFlag.isSet(INIT_MARKER)) { createMarker(); } if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) { - TEST_EXIT(!marker_)("marker already created\n"); - marker_ = adoptProblem->getMarker(); + TEST_EXIT(!marker)("marker already created\n"); + marker = adoptProblem->getMarker(); } } - if (!marker_) + if (!marker) WARNING("no marker created\n"); // === create file writer === @@ -363,27 +394,27 @@ namespace AMDiS { deserialize(in); in.close(); } else { - if (initFlag.isSet(INIT_MESH) && mesh_ && !mesh_->isInitialized()) { - mesh_->initialize(); + if (initFlag.isSet(INIT_MESH) && mesh && !mesh->isInitialized()) { + mesh->initialize(); } // === read value file and use it for the mesh values === std::string valueFilename(""); - GET_PARAMETER(0, mesh_->getName() + "->value file name", &valueFilename); + GET_PARAMETER(0, mesh->getName() + "->value file name", &valueFilename); if (valueFilename.length()) { ValueReader::readValue(valueFilename.c_str(), - mesh_, - solution_, - mesh_->getMacroFileInfo()); - mesh_->clearMacroFileInfo(); + mesh, + solution, + mesh->getMacroFileInfo()); + mesh->clearMacroFileInfo(); } // === do global refinements === if (initFlag.isSet(INIT_GLOBAL_REFINES)) { int globalRefinements = 0; - GET_PARAMETER(0, mesh_->getName() + "->global refinements", "%d", &globalRefinements); - refinementManager_->globalRefine(mesh_, globalRefinements); + GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinements); + refinementManager_->globalRefine(mesh, globalRefinements); } } } @@ -395,30 +426,30 @@ namespace AMDiS { int degree = 1; GET_PARAMETER(0, name_ + "->polynomial degree" ,"%d", °ree); feSpace_ = FiniteElemSpace::provideFESpace(NULL, - Lagrange::getLagrange(mesh_->getDim(), degree), - mesh_, + Lagrange::getLagrange(mesh->getDim(), degree), + mesh, name_ + "->feSpace"); // create dof admin for vertex dofs if neccessary - if (mesh_->getNumberOfDOFs(VERTEX) == 0) { - DimVec<int> ln_dof(mesh_->getDim(), DEFAULT_VALUE, 0); + if (mesh->getNumberOfDOFs(VERTEX) == 0) { + DimVec<int> ln_dof(mesh->getDim(), DEFAULT_VALUE, 0); ln_dof[VERTEX]= 1; - mesh_->createDOFAdmin("vertex dofs", ln_dof); + mesh->createDOFAdmin("vertex dofs", ln_dof); } } void ProblemScal::createMatricesAndVectors() { // === create vectors and system matrix === - systemMatrix_ = NEW DOFMatrix(feSpace_, feSpace_, "system matrix"); - rhs_ = NEW DOFVector<double>(feSpace_, "rhs"); - solution_ = NEW DOFVector<double>(feSpace_, "solution"); + systemMatrix = NEW DOFMatrix(feSpace_, feSpace_, "system matrix"); + rhs = NEW DOFVector<double>(feSpace_, "rhs"); + solution = NEW DOFVector<double>(feSpace_, "solution"); - solution_->setCoarsenOperation(COARSE_INTERPOL); - solution_->set(0.0); /* initialize u_h ! */ + solution->setCoarsenOperation(COARSE_INTERPOL); + solution->set(0.0); /* initialize u_h ! */ // === create matVec === - matVec_ = NEW StandardMatVec<DOFMatrix, DOFVector<double> >(systemMatrix_); + matVec_ = NEW StandardMatVec<DOFMatrix, DOFVector<double> >(systemMatrix); } void ProblemScal::createSolver() @@ -430,8 +461,8 @@ namespace AMDiS { dynamic_cast<OEMSolverCreator<DOFVector<double> >*>(CreatorMap<OEMSolver<DOFVector<double> > >::getCreator(solverType)); TEST_EXIT(solverCreator)("no solver type\n"); solverCreator->setName(name_ + "->solver"); - solver_ = solverCreator->create(); - solver_->initParameters(); + solver = solverCreator->create(); + solver->initParameters(); // === create preconditioners === std::string preconType("no"); @@ -448,8 +479,8 @@ namespace AMDiS { precon = preconCreator->create(); if (precon) { - dynamic_cast<PreconditionerScal*>(precon)->setMatrix(&systemMatrix_); - leftPrecon_ = precon; + dynamic_cast<PreconditionerScal*>(precon)->setMatrix(&systemMatrix); + leftPrecon = precon; } preconType.assign("no"); @@ -463,23 +494,19 @@ namespace AMDiS { precon = preconCreator->create(); if (precon) { - dynamic_cast<PreconditionerScal*>(precon)->setMatrix(&systemMatrix_); - rightPrecon_ = precon; + dynamic_cast<PreconditionerScal*>(precon)->setMatrix(&systemMatrix); + rightPrecon = precon; } // === create vector creator === - solver_->setVectorCreator(new DOFVector<double>::Creator(feSpace_)); + solver->setVectorCreator(new DOFVector<double>::Creator(feSpace_)); } void ProblemScal::createEstimator() { // create and set leaf data prototype - mesh_-> - setElementDataPrototype(NEW LeafDataEstimatable(NEW LeafDataCoarsenable)); - - // create estimator - //estimator = NEW ResidualEstimator(name + "->estimator"); + mesh->setElementDataPrototype(NEW LeafDataEstimatable(NEW LeafDataCoarsenable)); // === create estimator === std::string estimatorType("no"); @@ -490,28 +517,28 @@ namespace AMDiS { if (estimatorCreator) { estimatorCreator->setName(name_ + "->estimator"); if (estimatorType == "recovery") { - dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->setSolution(solution_); + dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->setSolution(solution); } - estimator_ = estimatorCreator->create(); + estimator = estimatorCreator->create(); // init estimator - estimator_->addSystem(systemMatrix_, solution_, rhs_); + estimator->addSystem(systemMatrix, solution, rhs); } } void ProblemScal::createMarker() { - marker_ = dynamic_cast<Marker*>(Marker::createMarker(name_ + "->marker", -1)); + marker = dynamic_cast<Marker*>(Marker::createMarker(name_ + "->marker", -1)); } void ProblemScal::createFileWriter() { - fileWriters_.push_back(NEW FileWriter(name_ + "->output", mesh_, solution_)); + fileWriters.push_back(NEW FileWriter(name_ + "->output", mesh, solution)); int writeSerialization = 0; GET_PARAMETER(0, name_ + "->output->write serialization", "%d", &writeSerialization); if (writeSerialization) { - fileWriters_.push_back(NEW Serializer<ProblemScal>(this)); + fileWriters.push_back(NEW Serializer<ProblemScal>(this)); } } @@ -519,22 +546,22 @@ namespace AMDiS { { FUNCNAME("Problem::estimate()"); - if (estimator_) { + if (estimator) { clock_t first = clock(); - estimator_->estimate(adaptInfo->getTimestep()); + estimator->estimate(adaptInfo->getTimestep()); INFO(info_,8)("estimation of the error needed %.5f seconds\n", TIME_USED(first,clock())); - adaptInfo->setEstSum(estimator_->getErrorSum(), 0); + adaptInfo->setEstSum(estimator->getErrorSum(), 0); adaptInfo-> - setTimeEstSum(estimator_->getTimeEst(), 0); + setTimeEstSum(estimator->getTimeEst(), 0); adaptInfo-> - setEstMax(estimator_->getErrorMax(), 0); + setEstMax(estimator->getErrorMax(), 0); adaptInfo-> - setTimeEstMax(estimator_->getTimeEstMax(), 0); + setTimeEstMax(estimator->getTimeEstMax(), 0); } else { WARNING("no estimator\n"); @@ -548,7 +575,7 @@ namespace AMDiS { clock_t first = clock(); - mesh_->dofCompress(); + mesh->dofCompress(); MSG("%d DOFs for %s\n", feSpace_->getAdmin()->getUsedSize(), @@ -556,28 +583,28 @@ namespace AMDiS { Flag assembleFlag = flag | - systemMatrix_->getAssembleFlag() | - rhs_->getAssembleFlag(); + systemMatrix->getAssembleFlag() | + rhs->getAssembleFlag(); if (useGetBound_) assembleFlag |= Mesh::FILL_BOUND; - systemMatrix_->clear(); - rhs_->set(0.0); + systemMatrix->clear(); + rhs->set(0.0); traversePtr_ = this; - mesh_->traverse(-1, assembleFlag, &buildAfterCoarsenFct); + mesh->traverse(-1, assembleFlag, &buildAfterCoarsenFct); // fill boundary conditions - if (systemMatrix_->getBoundaryManager()) - systemMatrix_->getBoundaryManager()->initMatrix(systemMatrix_); - if (rhs_->getBoundaryManager()) - rhs_->getBoundaryManager()->initVector(rhs_); - if (solution_->getBoundaryManager()) - solution_->getBoundaryManager()->initVector(solution_); + if (systemMatrix->getBoundaryManager()) + systemMatrix->getBoundaryManager()->initMatrix(systemMatrix); + if (rhs->getBoundaryManager()) + rhs->getBoundaryManager()->initVector(rhs); + if (solution->getBoundaryManager()) + solution->getBoundaryManager()->initVector(solution); TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh_, -1, + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS | @@ -587,23 +614,23 @@ namespace AMDiS { // for all elements ... while (elInfo) { - if (systemMatrix_->getBoundaryManager()) - systemMatrix_->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix_); - if (rhs_->getBoundaryManager()) - rhs_->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs_); - if (solution_->getBoundaryManager()) - solution_->getBoundaryManager()->fillBoundaryConditions(elInfo, solution_); + if (systemMatrix->getBoundaryManager()) + systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix); + if (rhs->getBoundaryManager()) + rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs); + if (solution->getBoundaryManager()) + solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution); elInfo = stack.traverseNext(elInfo); } - systemMatrix_->removeRowsWithDBC(systemMatrix_->getApplyDBCs()); + systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs()); - if (systemMatrix_->getBoundaryManager()) - systemMatrix_->getBoundaryManager()->exitMatrix(systemMatrix_); - if (rhs_->getBoundaryManager()) - rhs_->getBoundaryManager()->exitVector(rhs_); - if (solution_->getBoundaryManager()) - solution_->getBoundaryManager()->exitVector(solution_); + if (systemMatrix->getBoundaryManager()) + systemMatrix->getBoundaryManager()->exitMatrix(systemMatrix); + if (rhs->getBoundaryManager()) + rhs->getBoundaryManager()->exitVector(rhs); + if (solution->getBoundaryManager()) + solution->getBoundaryManager()->exitVector(solution); INFO(info_, 8)("buildAfterCoarsen needed %.5f seconds\n", TIME_USED(first,clock())); @@ -661,16 +688,16 @@ namespace AMDiS { { FUNCNAME("ProblemScal::serialize()"); - mesh_->serialize(out); - solution_->serialize(out); + mesh->serialize(out); + solution->serialize(out); } void ProblemScal::deserialize(std::istream &in) { FUNCNAME("ProblemScal::deserialize()"); - mesh_->deserialize(in); - solution_->deserialize(in); + mesh->deserialize(in); + solution->deserialize(in); } } diff --git a/AMDiS/src/ProblemScal.h b/AMDiS/src/ProblemScal.h index fb6e52cc13b1aec7fb8f6976d66b8d925fa28235..510d8e4f54f05fa6e43c704c3a5c01b30c2a3334 100644 --- a/AMDiS/src/ProblemScal.h +++ b/AMDiS/src/ProblemScal.h @@ -57,32 +57,29 @@ namespace AMDiS { ProblemScal(const char *name, ProblemIterationInterface *problemIteration = NULL) - : StandardProblemIteration(this) , + : StandardProblemIteration(this), name_(name), feSpace_(NULL), - mesh_(NULL), - marker_(NULL), - estimator_(NULL), - solver_(NULL), - solution_(NULL), - rhs_(NULL), - systemMatrix_(NULL), + mesh(NULL), + marker(NULL), + estimator(NULL), + solver(NULL), + solution(NULL), + rhs(NULL), + systemMatrix(NULL), matVec_(NULL), - leftPrecon_(NULL), - rightPrecon_(NULL), + leftPrecon(NULL), + rightPrecon(NULL), useGetBound_(true), refinementManager_(NULL), coarseningManager_(NULL), info_(10) - { - }; + {} /** \brief * Destructor */ - virtual ~ProblemScal() { - FUNCNAME("ProblemScal::~ProblemScal()"); - }; + virtual ~ProblemScal(); /** \brief * Initialisation of the problem. @@ -160,13 +157,13 @@ namespace AMDiS { * Implementation of ProblemStatBase::buildBeforeRefine(). * Does nothing here. */ - virtual void buildBeforeRefine(AdaptInfo *adaptInfo, Flag) {}; + virtual void buildBeforeRefine(AdaptInfo *adaptInfo, Flag) {} /** \brief * Implementation of ProblemStatBase::buildBeforeCoarsen(). * Does nothing here. */ - virtual void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag) {}; + virtual void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag) {} /** \brief * Implementation of ProblemStatBase::buildAfterCoarsen(). @@ -210,14 +207,14 @@ namespace AMDiS { interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct); /** \brief - * Adds an Operator to \ref systemMatrix_. + * Adds an Operator to \ref systemMatrix. */ void addMatrixOperator(Operator *op, double *factor = NULL, double *estFactor = NULL); /** \brief - * Adds an Operator to \ref rhs_. + * Adds an Operator to \ref rhs. */ void addVectorOperator(Operator *op, double *factor = NULL, @@ -265,99 +262,99 @@ namespace AMDiS { */ /** \brief - * Returns \ref solution_. + * Returns \ref solution. */ inline DOFVector<double>* getSolution() { - return solution_; - }; + return solution; + } /** \brief - * Returns \ref rhs_. + * Returns \ref rhs. */ inline DOFVector<double>* getRHS() { - return rhs_; - }; + return rhs; + } /** \brief - * Returns \ref systemMatrix_. + * Returns \ref systemMatrix. */ inline DOFMatrix *getSystemMatrix() { - return systemMatrix_; - }; + return systemMatrix; + } /** \brief - * Returns \ref mesh_ + * Returns \ref mesh */ inline Mesh* getMesh() { - return mesh_; - }; + return mesh; + } /** \brief * Returns \ref feSpace_. */ inline FiniteElemSpace* getFESpace() { return feSpace_; - }; + } /** \brief * Returns \ref estimator_. */ inline Estimator* getEstimator() { - return estimator_; - }; + return estimator; + } /** \brief * Returns \ref refinementManager_. */ inline RefinementManager* getRefinementManager() { return refinementManager_; - }; + } /** \brief - * Returns \ref solver_. + * Returns \ref solver. */ inline OEMSolver<DOFVector<double> >* getSolver() { - return solver_; - }; + return solver; + } /** \brief * Returns \ref marker_. */ inline Marker *getMarker() { - return marker_; - }; + return marker; + } /** \brief * Returns \ref name_. */ virtual inline const std::string& getName() { return name_; - }; + } /** \brief * Returns \ref useGetBound_. */ inline bool getBoundUsed() { return useGetBound_; - }; + } /** \brief - * Returns \ref leftPrecon_. + * Returns \ref leftPrecon. */ inline Preconditioner<DOFVector<double> > *getLeftPrecon() { - return leftPrecon_; - }; + return leftPrecon; + } /** \brief - * Returns \ref rightPrecon_. + * Returns \ref rightPrecon. */ inline Preconditioner<DOFVector<double> > *getRightPrecon() { - return rightPrecon_; - }; + return rightPrecon; + } inline CoarseningManager *getCoarseningManager() { return coarseningManager_; - }; + } /** \} */ @@ -368,7 +365,7 @@ namespace AMDiS { */ /** \brief - * Sets \ref mesh_ + * Sets \ref mesh */ void setMeshFromProblemVec(ProblemVec* pv, int i = 0); @@ -377,38 +374,40 @@ namespace AMDiS { */ inline void setFESpace(FiniteElemSpace* fe) { feSpace_ = fe; - }; + } /** \brief * Sets \ref estimator_. */ inline void setEstimator(Estimator* est) { - estimator_ = est; - }; + estimator = est; + } /** \brief - * Sets \ref solver_. + * Sets \ref solver. */ - inline void setSolver(OEMSolver<DOFVector<double> >* sol) { solver_ = sol; }; + inline void setSolver(OEMSolver<DOFVector<double> >* sol) { + solver = sol; + } /** \brief - * Sets \ref leftPrecon_. + * Sets \ref leftPrecon. */ inline void setLeftPrecon(Preconditioner<DOFVector<double> > *p) { - leftPrecon_ = p; - }; + leftPrecon = p; + } /** \brief - * Sets \ref rightPrecon_. + * Sets \ref rightPrecon. */ inline void setRightPrecon(Preconditioner<DOFVector<double> > *p) { - rightPrecon_ = p; - }; + rightPrecon = p; + } - inline void setMarker(Marker *marker) { - marker_ = marker; - }; + inline void setMarker(Marker *m) { + marker = m; + } /** \} */ @@ -427,8 +426,8 @@ namespace AMDiS { virtual void deserialize(std::istream &in); std::vector<FileWriterInterface*> getFileWriterList() { - return fileWriters_; - }; + return fileWriters; + } protected: /** \brief @@ -450,37 +449,37 @@ namespace AMDiS { /** \brief * Mesh of this problem. */ - Mesh *mesh_; + Mesh *mesh; /** \brief * Responsible for element marking. */ - Marker *marker_; + Marker *marker; /** \brief * Estimator of this problem. Used in \ref estimate(). */ - Estimator *estimator_; + Estimator *estimator; /** \brief * Linear solver of this problem. Used in \ref solve(). */ - OEMSolver<DOFVector<double> > *solver_; + OEMSolver<DOFVector<double> > *solver; /** \brief * DOFVector storing the calculated solution of the problem. */ - DOFVector<double> *solution_; + DOFVector<double> *solution; /** \brief * DOFVector for the right hand side */ - DOFVector<double> *rhs_; + DOFVector<double> *rhs; /** \brief * System matrix */ - DOFMatrix *systemMatrix_; + DOFMatrix *systemMatrix; /** \brief * Matrix-vector multiplication @@ -490,12 +489,12 @@ namespace AMDiS { /** \brief * Left preconditioner. Used in \ref solver. */ - Preconditioner<DOFVector<double> > *leftPrecon_; + Preconditioner<DOFVector<double> > *leftPrecon; /** \brief * Right preconditioner. Used in \ref solver. */ - Preconditioner<DOFVector<double> > *rightPrecon_; + Preconditioner<DOFVector<double> > *rightPrecon; /** \brief * Determines whether domain boundaries should be considered at assembling. @@ -505,7 +504,7 @@ namespace AMDiS { /** \brief * Writes the meshes and solution after the adaption loop. */ - std::vector<FileWriterInterface*> fileWriters_; + std::vector<FileWriterInterface*> fileWriters; /** \brief * All actions of mesh refinement are performed by refinementManager. diff --git a/AMDiS/src/RCNeighbourList.cc b/AMDiS/src/RCNeighbourList.cc index 24f1c17f2528ae0af620a336c3075e74f2925cc1..4ae838baf83d58cb720abd66a8c87bbb76e55734 100644 --- a/AMDiS/src/RCNeighbourList.cc +++ b/AMDiS/src/RCNeighbourList.cc @@ -16,15 +16,15 @@ namespace AMDiS { RCNeighbourList::RCNeighbourList(int maxEdgeNeigh) { rclist.resize(maxEdgeNeigh); - int i; - for(i = 0; i < maxEdgeNeigh; i++) + for (int i = 0; i < maxEdgeNeigh; i++) { rclist[i] = new RCListElement; + } } RCNeighbourList::~RCNeighbourList() { - int i, size = static_cast<int>(rclist.size()); - for(i = 0; i < size; i++) + for (int i = 0; i < static_cast<int>(rclist.size()); i++) { delete rclist[i]; + } } /****************************************************************************/ @@ -34,117 +34,105 @@ namespace AMDiS { bool RCNeighbourList::doCoarsePatch(int n_neigh) { - FUNCNAME("RCNeighbourList::doCoarsePatch"); - int i, j; - Element *lel; + FUNCNAME("RCNeighbourList::doCoarsePatch()"); - for (i = 0; i < n_neigh; i++) - { - lel =rclist[i]->el; + for (int i = 0; i < n_neigh; i++) { + Element *lel = rclist[i]->el; - if (lel->getMark() >= 0 || lel->isLeaf()) - { - /****************************************************************************/ - /* element must not be coarsend or element is a leaf element, reset the */ - /* the coarsening flag on all those elements that have to be coarsend with */ - /* this element */ - /****************************************************************************/ - lel->setMark(0); - for (j = 0; j < n_neigh; j++) - if (rclist[j]->flag) rclist[j]->el->setMark(0); - return(false); - } - else if (lel->getFirstChild()->getMark() >= 0 || - lel->getSecondChild()->getMark() >= 0) - { - /****************************************************************************/ - /* one of the element's children must not be coarsend; reset the coarsening*/ - /* flag on all those elements that have to be coarsend with this element */ - /****************************************************************************/ - lel->setMark(0); - for (j = 0; j < n_neigh; j++) - if (rclist[j]->flag) rclist[j]->el->setMark(0); - return(false); - } - else if (!lel->getFirstChild()->isLeaf() || - !lel->getSecondChild()->isLeaf()) - { - /****************************************************************************/ - /* one of the element's children is not a leaf element; */ - /* element may be coarsend after coarsening one of the children; try again */ - /****************************************************************************/ - coarseningManager->doMore = true; - return(false); - } - else - { - /****************************************************************************/ - /* either one element is a macro element or we can coarsen the patch */ - /****************************************************************************/ - if (rclist[i]->flag == 0) { + if (lel->getMark() >= 0 || lel->isLeaf()) { + /****************************************************************************/ + /* element must not be coarsend or element is a leaf element, reset the */ + /* the coarsening flag on all those elements that have to be coarsend with */ + /* this element */ + /****************************************************************************/ + lel->setMark(0); + for (int j = 0; j < n_neigh; j++) + if (rclist[j]->flag) + rclist[j]->el->setMark(0); - std::deque<MacroElement*>::const_iterator mel; + return false; + } else if (lel->getFirstChild()->getMark() >= 0 || + lel->getSecondChild()->getMark() >= 0) { + + /****************************************************************************/ + /* one of the element's children must not be coarsend; reset the coarsening*/ + /* flag on all those elements that have to be coarsend with this element */ + /****************************************************************************/ + lel->setMark(0); + for (int j = 0; j < n_neigh; j++) + if (rclist[j]->flag) + rclist[j]->el->setMark(0); + + return false; + } else if (!lel->getFirstChild()->isLeaf() || + !lel->getSecondChild()->isLeaf()) { + /****************************************************************************/ + /* one of the element's children is not a leaf element; */ + /* element may be coarsend after coarsening one of the children; try again */ + /****************************************************************************/ + coarseningManager->doMore = true; + + return false; + } else { + /****************************************************************************/ + /* either one element is a macro element or we can coarsen the patch */ + /****************************************************************************/ + if (rclist[i]->flag == 0) { + + std::deque<MacroElement*>::const_iterator mel; - Mesh* coarse_mesh = coarseningManager->getMesh(); - // set in Mesh::coarsen() - for (mel = coarse_mesh->firstMacroElement(); - mel!=coarse_mesh->endOfMacroElements(); - mel++) - if ((*mel)->getElement() == lel) break; - - TEST_EXIT(mel!=coarse_mesh->endOfMacroElements()) - ("incompatible coarsening patch found\n"); - } - } + Mesh* coarse_mesh = coarseningManager->getMesh(); + // set in Mesh::coarsen() + for (mel = coarse_mesh->firstMacroElement(); + mel!=coarse_mesh->endOfMacroElements(); + mel++) + if ((*mel)->getElement() == lel) + break; + + TEST_EXIT(mel!=coarse_mesh->endOfMacroElements()) + ("incompatible coarsening patch found\n"); + } } + } - return(true); + return true; } void RCNeighbourList::getNeighOnPatch(int n_neigh, int bound) { - FUNCNAME("RCNeighbourList::getNeighOnPatch"); - - int i, j, k, dir; - Element *el, *neigh; - - //Mesh *mesh = rclist[0]->el->getMesh(); - - for (i = 0; i < n_neigh; i++) - { - el = rclist[i]->el; - rclist[i]->no = i; - - for (dir = 0; dir < 2; dir++) - { - for (j = 0; j < n_neigh; j++) - { - if ((neigh = rclist[j]->el) == el) continue; - - for (k = 0; k < 2; k++) - { - if (neigh->getDOF(2+k) == el->getDOF(3-dir)) - { - rclist[i]->neigh[dir] = rclist[j]; - rclist[i]->oppVertex[dir] = 3-k; - break; - } - } - - if (k < 2) break; - } - - if (j >= n_neigh) - { - // TEST_EXIT(bound) - // ("neighbour of element %d in list not found\n", el->getIndex()); - rclist[i]->neigh[dir] = NULL; - rclist[i]->oppVertex[dir] = -1; - } + FUNCNAME("RCNeighbourList::getNeighOnPatch()"); + + int j, k, dir; + Element *el, *neigh; + + for (int i = 0; i < n_neigh; i++) { + el = rclist[i]->el; + rclist[i]->no = i; + + for (dir = 0; dir < 2; dir++) { + for (j = 0; j < n_neigh; j++) { + if ((neigh = rclist[j]->el) == el) + continue; + + for (k = 0; k < 2; k++) { + if (neigh->getDOF(2+k) == el->getDOF(3-dir)) { + rclist[i]->neigh[dir] = rclist[j]; + rclist[i]->oppVertex[dir] = 3-k; + break; + } } + + if (k < 2) + break; + } + + if (j >= n_neigh) { + rclist[i]->neigh[dir] = NULL; + rclist[i]->oppVertex[dir] = -1; + } } - return; + } } void RCNeighbourList::addDOFParent(int elIndex, DegreeOfFreedom* dof) // 3d @@ -348,17 +336,10 @@ namespace AMDiS { int *n_neigh, int *n_neigh_periodic) { - static RCNeighbourList *periodicList = NULL; - - if(periodicList) { - DELETE periodicList; - periodicList = NULL; - } - + RCNeighbourList *periodicList = NULL; *n_neigh_periodic = 0; - - int count = 0, n_neigh_old = *n_neigh; - + int count = 0; + int n_neigh_old = *n_neigh; bool secondPart = false; nextEdge[0] = NULL; @@ -367,13 +348,13 @@ namespace AMDiS { std::vector<RCListElement*>::iterator it = rclist.begin(); std::vector<RCListElement*>::iterator insertIt; - while(count < n_neigh_old) { + while (count < n_neigh_old) { DegreeOfFreedom *dof0 = const_cast<DegreeOfFreedom*>((*it)->el->getDOF(0)); DegreeOfFreedom *dof1 = const_cast<DegreeOfFreedom*>((*it)->el->getDOF(1)); - if(dof0 != edge[0] && dof0 != edge[1]) { + if (dof0 != edge[0] && dof0 != edge[1]) { secondPart = true; - if(!nextEdge[0]) { + if (!nextEdge[0]) { nextEdge[0] = dof0; nextEdge[1] = dof1; } @@ -382,7 +363,7 @@ namespace AMDiS { (*n_neigh)--; (*n_neigh_periodic)++; - if(!periodicList) { + if (!periodicList) { periodicList = NEW RCNeighbourList; insertIt = periodicList->rclist.end(); periodicList->coarseningManager = coarseningManager; @@ -394,7 +375,7 @@ namespace AMDiS { } } - if(secondPart) { + if (secondPart) { insertIt = periodicList->rclist.begin(); secondPart = false; } @@ -404,11 +385,10 @@ namespace AMDiS { it = rclist.erase(it); } - ++count; + count++; } - int i; - for(i = 0; i < *n_neigh_periodic; i++) { + for (int i = 0; i < *n_neigh_periodic; i++) { periodicList->rclist[i]->no = i; } diff --git a/AMDiS/src/RCNeighbourList.h b/AMDiS/src/RCNeighbourList.h index dd84cd6a0732f985ba090228520a9bcabc4e1e5f..a32dab3b0f8f263d931910b0f3634673ff1e3baf 100644 --- a/AMDiS/src/RCNeighbourList.h +++ b/AMDiS/src/RCNeighbourList.h @@ -63,7 +63,7 @@ namespace AMDiS { */ RCNeighbourList(int maxEdgeNeigh); - RCNeighbourList() {}; + RCNeighbourList() {} /** \brief * Destructor @@ -92,7 +92,7 @@ namespace AMDiS { */ inline int getNeighbourNr(int i, int j) const { return (rclist[i]->neigh[j])?rclist[i]->neigh[j]->no:-1; - }; + } /** \brief * If \ref rclist[i].neigh[j] is not a NULL pointer @@ -101,33 +101,36 @@ namespace AMDiS { */ inline Element* getNeighbourElement(int i, int j) const { return rclist[i]->neigh[j]? rclist[i]->neigh[j]->el : NULL; - }; + } /** \brief * Returns \ref rclist[i].el */ inline Element* getElement(int i) const { - if(static_cast<int>(rclist.size()) <= i) return NULL; + if (static_cast<int>(rclist.size()) <= i) + return NULL; return rclist[i]->el; - }; + } /** \brief * Sets \ref rclist[i].el to el and \ref rclist[i].flag to cp. */ - inline const Element* setElement(int i,const Element* el,bool cp=false) { - rclist[i]->el=const_cast<Element*>(el); - rclist[i]->flag=cp; + inline const Element* setElement(int i, const Element* el, bool cp = false) { + rclist[i]->el = const_cast<Element*>(el); + rclist[i]->flag = cp; return el; - }; + } /** \brief * Returns \ref rclist[i].elType */ - inline int getType(int i) const {return rclist[i]->elType;}; + inline int getType(int i) const { + return rclist[i]->elType; + } inline void setType(int i, int type) const { rclist[i]->elType = type; - }; + } /** \brief * If patch can be coarsend return true, else false and reset the element diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc index 4154998727cdf75a30a05a05f309fc1f098d15d5..0699e56235d12984c323dff7854a6b18e00deca7 100644 --- a/AMDiS/src/RefinementManager2d.cc +++ b/AMDiS/src/RefinementManager2d.cc @@ -26,8 +26,10 @@ namespace AMDiS { RCNeighbourList* refineList = NEW RCNeighbourList(2); - if (el_info->getElement()->getMark() <= 0) - return(el_info); /* element may not be refined */ + if (el_info->getElement()->getMark() <= 0) { + DELETE refineList; + return el_info; /* element may not be refined */ + } refineList->setElement(0, el_info->getElement()); n_neigh = 1; @@ -87,6 +89,8 @@ namespace AMDiS { newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound); + DELETE periodicList; + if (firstNewDOF == -1) { firstNewDOF = newDOF; } @@ -136,7 +140,7 @@ namespace AMDiS { DELETE refineList; - return(el_info); + return el_info; } @@ -181,8 +185,8 @@ namespace AMDiS { { FUNCNAME("RefinementManager2d::refinePatch()"); DegreeOfFreedom *dof[3] = {NULL, NULL, NULL}; - Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( refineList->getElement(0))); - Triangle *neigh = dynamic_cast<Triangle*>(const_cast<Element*>( refineList->getElement(1))); + Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(0))); + Triangle *neigh = dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(1))); /****************************************************************************/ /* there is one new vertex in the refinement edge */ diff --git a/AMDiS/src/SolutionDataStorage.h b/AMDiS/src/SolutionDataStorage.h index 3fb25ee5ae5cf9c1253c1ae729600136ed51accb..d54309c6974cf639e3d04209f62ef23444bec8a4 100644 --- a/AMDiS/src/SolutionDataStorage.h +++ b/AMDiS/src/SolutionDataStorage.h @@ -81,6 +81,10 @@ namespace AMDiS { return feSpaces[i]; } + bool isPoped() { + return poped; + } + protected: /** \brief * Deletes the fe space and all it content, i.e., also the dof admin, diff --git a/AMDiS/src/SolutionDataStorage.hh b/AMDiS/src/SolutionDataStorage.hh index 805b6de5a66a1d8c7a1c3de8333fbdc7ec75af1d..62fcb18eaefbb44e55b8deeef1d68317f60fe86c 100644 --- a/AMDiS/src/SolutionDataStorage.hh +++ b/AMDiS/src/SolutionDataStorage.hh @@ -23,6 +23,7 @@ namespace AMDiS { template<typename T> SolutionDataStorage<T>::~SolutionDataStorage() { + clear(); } template<typename T> @@ -47,7 +48,7 @@ namespace AMDiS { clear(); poped = false; } - + solutions.push_back(solution); timestamps.push_back(timestamp); @@ -132,13 +133,7 @@ namespace AMDiS { void SolutionDataStorage<T>::deleteFeSpace(FiniteElemSpace *feSpace) { if (feSpace) { - if (feSpace->getMesh()) { - DELETE feSpace->getMesh(); - } - if (feSpace->getAdmin()) { - DELETE feSpace->getAdmin(); - } - + DELETE feSpace->getMesh(); DELETE feSpace; } } diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h index 67fb780c9b68b89854b8aa0c82f86a353f30022e..0c549f6805cf733a28c8f0e4db37ef95b3b5af25 100644 --- a/AMDiS/src/Triangle.h +++ b/AMDiS/src/Triangle.h @@ -45,7 +45,9 @@ namespace AMDiS { */ Triangle(Mesh* aMesh) : Element(aMesh) - {}; + {} + + ~Triangle() {} /** \brief * implements Element::clone