From ab48be0c88df1b13a4ed67b94c6c7587d4994598 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski <thomas.witkowski@gmx.de> Date: Mon, 8 Sep 2008 13:52:57 +0000 Subject: [PATCH] * On the way to independent meshes ... --- AMDiS/bin/Makefile.am | 1 - AMDiS/bin/Makefile.in | 29 +++---- AMDiS/src/Assembler.cc | 2 +- AMDiS/src/CoarseningManager2d.cc | 109 +++++------------------- AMDiS/src/DOFAdmin.cc | 15 ++-- AMDiS/src/DOFAdmin.h | 8 +- AMDiS/src/DualTraverse.cc | 51 ++++++++---- AMDiS/src/DualTraverse.h | 41 ++++++--- AMDiS/src/ElInfo.cc | 8 +- AMDiS/src/ElInfo.h | 51 +++++++++--- AMDiS/src/ElInfo1d.cc | 49 +++++++---- AMDiS/src/ElInfo1d.h | 6 ++ AMDiS/src/ElInfo2d.cc | 139 ++++++++++++++++--------------- AMDiS/src/ElInfo2d.h | 15 ++++ AMDiS/src/ElInfo3d.cc | 73 ++++++++-------- AMDiS/src/ElInfo3d.h | 15 ++++ AMDiS/src/FixVec.h | 9 +- AMDiS/src/Flag.cc | 5 -- AMDiS/src/Mesh.cc | 11 ++- AMDiS/src/ProblemVec.cc | 1 + AMDiS/src/Traverse.cc | 17 +++- AMDiS/src/Traverse.h | 3 + 22 files changed, 363 insertions(+), 295 deletions(-) delete mode 100644 AMDiS/src/Flag.cc diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am index 91b0407e..acdfdb44 100644 --- a/AMDiS/bin/Makefile.am +++ b/AMDiS/bin/Makefile.am @@ -197,7 +197,6 @@ $(SOURCE_DIR)/ElInfo2d.cc \ $(SOURCE_DIR)/ElInfo3d.cc \ $(SOURCE_DIR)/FiniteElemSpace.cc \ $(SOURCE_DIR)/FixVec.cc \ -$(SOURCE_DIR)/Flag.cc \ $(SOURCE_DIR)/Global.cc \ $(SOURCE_DIR)/Lagrange.cc \ $(SOURCE_DIR)/Line.cc \ diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in index e1cf7751..a7441df7 100644 --- a/AMDiS/bin/Makefile.in +++ b/AMDiS/bin/Makefile.in @@ -227,11 +227,11 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ConditionalEstimator.h \ $(SOURCE_DIR)/Element.cc $(SOURCE_DIR)/ElInfo1d.cc \ $(SOURCE_DIR)/ElInfo2d.cc $(SOURCE_DIR)/ElInfo3d.cc \ $(SOURCE_DIR)/FiniteElemSpace.cc $(SOURCE_DIR)/FixVec.cc \ - $(SOURCE_DIR)/Flag.cc $(SOURCE_DIR)/Global.cc \ - $(SOURCE_DIR)/Lagrange.cc $(SOURCE_DIR)/Line.cc \ - $(SOURCE_DIR)/MacroElement.cc $(SOURCE_DIR)/MacroWriter.cc \ - $(SOURCE_DIR)/Parameters.cc $(SOURCE_DIR)/Parametric.cc \ - $(SOURCE_DIR)/Quadrature.cc $(SOURCE_DIR)/RCNeighbourList.cc \ + $(SOURCE_DIR)/Global.cc $(SOURCE_DIR)/Lagrange.cc \ + $(SOURCE_DIR)/Line.cc $(SOURCE_DIR)/MacroElement.cc \ + $(SOURCE_DIR)/MacroWriter.cc $(SOURCE_DIR)/Parameters.cc \ + $(SOURCE_DIR)/Parametric.cc $(SOURCE_DIR)/Quadrature.cc \ + $(SOURCE_DIR)/RCNeighbourList.cc \ $(SOURCE_DIR)/RefinementManager.cc \ $(SOURCE_DIR)/RefinementManager1d.cc \ $(SOURCE_DIR)/RefinementManager2d.cc \ @@ -303,11 +303,11 @@ am_libamdis_la_OBJECTS = $(am__objects_1) \ libamdis_la-Element.lo libamdis_la-ElInfo1d.lo \ libamdis_la-ElInfo2d.lo libamdis_la-ElInfo3d.lo \ libamdis_la-FiniteElemSpace.lo libamdis_la-FixVec.lo \ - libamdis_la-Flag.lo libamdis_la-Global.lo \ - libamdis_la-Lagrange.lo libamdis_la-Line.lo \ - libamdis_la-MacroElement.lo libamdis_la-MacroWriter.lo \ - libamdis_la-Parameters.lo libamdis_la-Parametric.lo \ - libamdis_la-Quadrature.lo libamdis_la-RCNeighbourList.lo \ + libamdis_la-Global.lo libamdis_la-Lagrange.lo \ + libamdis_la-Line.lo libamdis_la-MacroElement.lo \ + libamdis_la-MacroWriter.lo libamdis_la-Parameters.lo \ + libamdis_la-Parametric.lo libamdis_la-Quadrature.lo \ + libamdis_la-RCNeighbourList.lo \ libamdis_la-RefinementManager.lo \ libamdis_la-RefinementManager1d.lo \ libamdis_la-RefinementManager2d.lo \ @@ -640,7 +640,6 @@ $(SOURCE_DIR)/ElInfo2d.cc \ $(SOURCE_DIR)/ElInfo3d.cc \ $(SOURCE_DIR)/FiniteElemSpace.cc \ $(SOURCE_DIR)/FixVec.cc \ -$(SOURCE_DIR)/Flag.cc \ $(SOURCE_DIR)/Global.cc \ $(SOURCE_DIR)/Lagrange.cc \ $(SOURCE_DIR)/Line.cc \ @@ -804,7 +803,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-FiniteElemSpace.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-FirstOrderAssembler.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-FixVec.Plo@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Flag.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-GNUPlotWriter.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-GSSmoother.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Global.Plo@am__quote@ @@ -1487,13 +1485,6 @@ libamdis_la-FixVec.lo: $(SOURCE_DIR)/FixVec.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-FixVec.lo `test -f '$(SOURCE_DIR)/FixVec.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/FixVec.cc -libamdis_la-Flag.lo: $(SOURCE_DIR)/Flag.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-Flag.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-Flag.Tpo" -c -o libamdis_la-Flag.lo `test -f '$(SOURCE_DIR)/Flag.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/Flag.cc; \ -@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-Flag.Tpo" "$(DEPDIR)/libamdis_la-Flag.Plo"; else rm -f "$(DEPDIR)/libamdis_la-Flag.Tpo"; exit 1; fi -@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/Flag.cc' object='libamdis_la-Flag.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-Flag.lo `test -f '$(SOURCE_DIR)/Flag.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/Flag.cc - libamdis_la-Global.lo: $(SOURCE_DIR)/Global.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-Global.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-Global.Tpo" -c -o libamdis_la-Global.lo `test -f '$(SOURCE_DIR)/Global.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/Global.cc; \ @am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-Global.Tpo" "$(DEPDIR)/libamdis_la-Global.Plo"; else rm -f "$(DEPDIR)/libamdis_la-Global.Tpo"; exit 1; fi diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc index c64faa91..db81e6fd 100644 --- a/AMDiS/src/Assembler.cc +++ b/AMDiS/src/Assembler.cc @@ -131,7 +131,7 @@ namespace AMDiS { if (zeroOrderAssembler) { zeroOrderAssembler->calculateElementMatrix(rowElInfo, mat); } - + if (rememberElMat && userMat) { axpy(factor, *elementMatrix, *userMat); } diff --git a/AMDiS/src/CoarseningManager2d.cc b/AMDiS/src/CoarseningManager2d.cc index 4b9c2c89..a54e24c6 100644 --- a/AMDiS/src/CoarseningManager2d.cc +++ b/AMDiS/src/CoarseningManager2d.cc @@ -80,123 +80,58 @@ namespace AMDiS { DegreeOfFreedom *dof[3]; dof[0] = const_cast<int*>( el->getChild(0)->getDOF(2)); - if (mesh->getNumberOfDOFs(EDGE)) - { - dof[1] = const_cast<int*>( el->getChild(0)->getDOF(3)); - dof[2] = const_cast<int*>( el->getChild(1)->getDOF(4)); - } + if (mesh->getNumberOfDOFs(EDGE)) { + dof[1] = const_cast<int*>( el->getChild(0)->getDOF(3)); + dof[2] = const_cast<int*>( el->getChild(1)->getDOF(4)); + } int node = mesh->getNode(EDGE); - if (mesh->getNumberOfDOFs(EDGE)) - { - /****************************************************************************/ - /* get new dof on el at the midpoint of the coarsening edge */ - /****************************************************************************/ - if (!el->getDOF(node+2)) - { - el->setDOF(node+2, mesh->getDOF(EDGE)); - if (neigh) { - // // periodic boundary ? - // if(el->getDOF(0) != neigh->getDOF(0) && - // el->getDOF(1) != neigh->getDOF(0)) - // { - // neigh->setDOF(node+2, mesh->getDOF(EDGE)); - // } else { - neigh->setDOF(node+2, const_cast<int*>( el->getDOF(node+2))); - // } - } - } + if (mesh->getNumberOfDOFs(EDGE)) { + /****************************************************************************/ + /* get new dof on el at the midpoint of the coarsening edge */ + /****************************************************************************/ + if (!el->getDOF(node+2)) { + el->setDOF(node+2, mesh->getDOF(EDGE)); + if (neigh) { + neigh->setDOF(node+2, const_cast<int*>( el->getDOF(node+2))); + } } - - // // periodic boundary? - // DegreeOfFreedom *neighDOF[3] = {NULL, NULL, NULL}; - // RCNeighbourList *periodicCoarsenList = NULL; - // int n_neigh_periodic; - // if(neigh && - // (neigh->getDOF(0) != el->getDOF(1)) && - // (neigh->getDOF(0) != el->getDOF(0))) - // { - // neighDOF[0] = const_cast<int*>(neigh->getChild(0)->getDOF(2)); - // if (mesh->getNumberOfDOFs(EDGE)) { - // neighDOF[1] = const_cast<int*>(neigh->getChild(0)->getDOF(3)); - // neighDOF[2] = const_cast<int*>(neigh->getChild(1)->getDOF(4)); - // if (!neigh->getDOF(node+2)) { - // neigh->setDOF(node+2, mesh->getDOF(EDGE)); - // } - // } - - // DegreeOfFreedom *edge[2] = { - // const_cast<DegreeOfFreedom*>(el->getDOF(0)), - // const_cast<DegreeOfFreedom*>(el->getDOF(1)) - // }; - // DegreeOfFreedom *periodicEdge[2] = { - // const_cast<DegreeOfFreedom*>(neigh->getDOF(0)), - // const_cast<DegreeOfFreedom*>(neigh->getDOF(1)) - // }; - - // periodicCoarsenList = coarsenList->periodicSplit(edge, - // &n_neigh, - // periodicEdge, - // &n_neigh_periodic); - // } + } if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) { coarsenList->addDOFParents(n_neigh); - // if(periodicCoarsenList) { - // periodicCoarsenList->addDOFParents(n_neigh_periodic); - // } } /****************************************************************************/ /* restrict dof vectors to the parents on the patch */ /****************************************************************************/ - int iadmin; int nrAdmin = mesh->getNumberOfDOFAdmin(); - for(iadmin = 0; iadmin < nrAdmin; iadmin++) { + for (int iadmin = 0; iadmin < nrAdmin; iadmin++) { std::list<DOFIndexedBase*>::iterator it; DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin)); std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed(); - for(it = admin->beginDOFIndexed(); it != end; ++it) { + for (it = admin->beginDOFIndexed(); it != end; ++it) { (*it)->coarseRestrict(*coarsenList, n_neigh); - // if(periodicCoarsenList) { - // (*it)->coarseRestrict(*periodicCoarsenList, n_neigh_periodic); - // } } } - // if(periodicCoarsenList) DELETE periodicCoarsenList; - coarsenTriangle(el); - if (neigh) coarsenTriangle(neigh); + if (neigh) + coarsenTriangle(neigh); /****************************************************************************/ /* now, remove those dofs in the coarcening edge */ /****************************************************************************/ mesh->freeDOF(dof[0], VERTEX); - if (mesh->getNumberOfDOFs(EDGE)) - { - mesh->freeDOF(dof[1], EDGE); - mesh->freeDOF(dof[2], EDGE); - } + if (mesh->getNumberOfDOFs(EDGE)) { + mesh->freeDOF(dof[1], EDGE); + mesh->freeDOF(dof[2], EDGE); + } mesh->incrementNumberOfVertices(-1); mesh->incrementNumberOfEdges(-1); - // // periodic boundary? - // if(neigh && - // (neigh->getDOF(0) != el->getDOF(1)) && - // (neigh->getDOF(0) != el->getDOF(0))) - // { - // mesh->freeDOF(neighDOF[0], VERTEX); - // if (mesh->getNumberOfDOFs(EDGE)) { - // mesh->freeDOF(neighDOF[1], EDGE); - // mesh->freeDOF(neighDOF[2], EDGE); - // } - // mesh->incrementNumberOfVertices(-1); - // mesh->incrementNumberOfEdges(-1); - // } - return; } diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc index 9b53871b..fe9f303d 100755 --- a/AMDiS/src/DOFAdmin.cc +++ b/AMDiS/src/DOFAdmin.cc @@ -86,20 +86,21 @@ namespace AMDiS { std::list<DOFIndexedBase*>::iterator di; std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end(); - for(di = dofIndexedList.begin(); di != end; ++di) { + for (di = dofIndexedList.begin(); di != end; ++di) { (*di)->freeDOFContent(dof); } std::list<DOFContainer*>::iterator dc; std::list<DOFContainer*>::iterator dcend = dofContainerList.end(); - for(dc = dofContainerList.begin(); dc != dcend; ++dc) { + for (dc = dofContainerList.begin(); dc != dcend; ++dc) { (*dc)->freeDOFIndex(dof); } dofFree[dof] = true; - if (static_cast<int>(firstHole) > dof) firstHole = dof; + if (static_cast<int>(firstHole) > dof) + firstHole = dof; usedCount--; holeCount++; @@ -109,7 +110,7 @@ namespace AMDiS { int DOFAdmin::getDOFIndex() { - FUNCNAME("DOFAdmin::getDOFIndex"); + FUNCNAME("DOFAdmin::getDOFIndex()"); int dof = 0; // if there is a hole @@ -120,8 +121,8 @@ namespace AMDiS { dof = firstHole; // search new hole int dfsize = static_cast<int>(dofFree.size()); - int i = 0; - for (i = firstHole + 1; i < dfsize; i++) { + int i = firstHole + 1; + for (; i < dfsize; i++) { if (dofFree[i]) { break; } @@ -152,7 +153,7 @@ namespace AMDiS { void DOFAdmin::enlargeDOFLists(int minsize) { - FUNCNAME("DOFAdmin::enlargeDOFLists"); + FUNCNAME("DOFAdmin::enlargeDOFLists()"); int old = size; if (minsize > 0) { diff --git a/AMDiS/src/DOFAdmin.h b/AMDiS/src/DOFAdmin.h index 4e2d43b4..bb215a85 100644 --- a/AMDiS/src/DOFAdmin.h +++ b/AMDiS/src/DOFAdmin.h @@ -271,12 +271,16 @@ namespace AMDiS { /** \brief * Sets \ref name = n */ - inline void setName(const std::string& n) { name = n; }; + inline void setName(const std::string& n) { + name = n; + }; /** \brief * Sets \ref mesh = m */ - inline void setMesh(Mesh* m) { mesh = m; }; + inline void setMesh(Mesh* m) { + mesh = m; + }; /** \} */ diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc index 5a624099..ad1f1d43 100644 --- a/AMDiS/src/DualTraverse.cc +++ b/AMDiS/src/DualTraverse.cc @@ -69,7 +69,7 @@ namespace AMDiS { return false; } - rest_ = 1.0; + rest = 1.0; bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge); @@ -83,6 +83,8 @@ namespace AMDiS { return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge); } + fillSubElInfo(*elInfo1, *elInfo2, *elInfoSmall, *elInfoLarge); + return true; } @@ -92,12 +94,12 @@ namespace AMDiS { ElInfo **elInfoLarge) { // call standard traverse - if (inc1_) { + if (inc1) { do { *elInfo1 = stack1.traverseNext(*elInfo1); - } while(*elInfo1 != NULL && skipEl1(*elInfo1)); + } while (*elInfo1 != NULL && skipEl1(*elInfo1)); } - if (inc2_) { + if (inc2) { do { *elInfo2 = stack2.traverseNext(*elInfo2); } while (*elInfo2 != NULL && skipEl2(*elInfo2)); @@ -120,7 +122,9 @@ namespace AMDiS { (callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf())) { return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge); } - + + fillSubElInfo(*elInfo1, *elInfo2, *elInfoSmall, *elInfoLarge); + return true; } @@ -139,22 +143,39 @@ namespace AMDiS { *elInfo1 : *elInfo2; - int smallLevel = (*elInfoSmall)->getLevel(); - int largeLevel = (*elInfoLarge)->getLevel(); - // update rest - rest_ -= 1.0 / (1 << (smallLevel - largeLevel)); + rest -= 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel())); - if (rest_ < 1e-32) { + if (rest < 1e-32) { // large element finished -> increment both elements - rest_ = 1.0; - inc1_ = true; - inc2_ = true; + rest = 1.0; + inc1 = true; + inc2 = true; } else { // increment only small element - inc1_ = (*elInfo1 == *elInfoSmall) ? true : false; - inc2_ = (*elInfo2 == *elInfoSmall) ? true : false; + inc1 = (*elInfo1 == *elInfoSmall) ? true : false; + inc2 = (*elInfo2 == *elInfoSmall) ? true : false; } } + void DualTraverse::fillSubElInfo(ElInfo *elInfo1, + ElInfo *elInfo2, + ElInfo *elInfoSmall, + ElInfo *elInfoLarge) + { + if (!fillSubElemInfo) + return; + + VectorOfFixVecs<DimVec<double> > *subCoords = elInfoSmall->getSubElemCoords(); + if (!subCoords) { + subCoords = NEW VectorOfFixVecs<DimVec<double> >(1, 2, NO_INIT); + } + + if (elInfo1 == elInfoSmall) { + stack1.getCoordsInElem(elInfo2, subCoords); + } else { + stack2.getCoordsInElem(elInfo1, subCoords); + } + elInfoSmall->setSubElemCoords(subCoords); + } } diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h index ec321921..edd62d99 100644 --- a/AMDiS/src/DualTraverse.h +++ b/AMDiS/src/DualTraverse.h @@ -43,7 +43,11 @@ namespace AMDiS { public: MEMORY_MANAGED(DualTraverse); - virtual ~DualTraverse() {}; + DualTraverse() + : fillSubElemInfo(false) + {}; + + ~DualTraverse() {}; /** \brief * Start dual traversal @@ -67,23 +71,27 @@ namespace AMDiS { ElInfo **elInfoSmall, ElInfo **elInfoLarge); - virtual bool skipEl1(ElInfo *elInfo) { + bool skipEl1(ElInfo *elInfo) { return false; - }; + } - virtual bool skipEl2(ElInfo *elInfo) { + bool skipEl2(ElInfo *elInfo) { return false; - }; + } - virtual bool check(ElInfo **elInfo1, - ElInfo **elInfo2, - ElInfo **elInfoSmall, - ElInfo **elInfoLarge) + bool check(ElInfo **elInfo1, + ElInfo **elInfo2, + ElInfo **elInfoSmall, + ElInfo **elInfoLarge) { prepareNextStep(elInfo1, elInfo2, elInfoSmall, elInfoLarge); return true; - }; + } + inline void setFillSubElemInfo(bool b) { + fillSubElemInfo = b; + } + protected: /** \brief * Determines smaller and larger element, determines which element(s) has to @@ -94,6 +102,11 @@ namespace AMDiS { ElInfo **elInfoSmall, ElInfo **elInfoLarge); + void fillSubElInfo(ElInfo *elInfo1, + ElInfo *elInfo2, + ElInfo *elInfoSmall, + ElInfo *elInfoLarge); + protected: /** \brief * stack for mesh 1 @@ -109,17 +122,17 @@ namespace AMDiS { * used to determine whether all small elements belonging to the large * element are traversed. */ - double rest_; + double rest; /** \brief * true is element 1 should be incremented (set in prepareNextStep()) */ - bool inc1_; + bool inc1; /** \brief * true is element 2 should be incremented (set in prepareNextStep()) */ - bool inc2_; + bool inc2; /** \brief * for level traverse of mesh 1 @@ -140,6 +153,8 @@ namespace AMDiS { * for leaf element level traverse of mesh 2 */ bool callLeafElLevel2_; + + bool fillSubElemInfo; }; } diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc index aee887f8..8b46c346 100644 --- a/AMDiS/src/ElInfo.cc +++ b/AMDiS/src/ElInfo.cc @@ -22,6 +22,7 @@ namespace AMDiS { element_(NULL), parent_(NULL), macroElement_(NULL), + subElemCoords(NULL), coord_(mesh_->getDim(), NO_INIT), boundary_(mesh_->getDim(), DEFAULT_VALUE, INTERIOR), projection_(mesh_->getDim(), NO_INIT), @@ -29,7 +30,9 @@ namespace AMDiS { neighbour_(mesh_->getDim(), NO_INIT), neighbourCoord_(mesh_->getDim(), NO_INIT), oppVertex_(mesh_->getDim(), NO_INIT), - grdLambda_(mesh_->getDim(), NO_INIT) + grdLambda_(mesh_->getDim(), NO_INIT), + level_(0), + iChild(0) { projection_.set(NULL); @@ -42,6 +45,9 @@ namespace AMDiS { ElInfo::~ElInfo() { + if (subElemCoords) { + DELETE subElemCoords; + } } diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h index 5ac0f275..40e1e283 100644 --- a/AMDiS/src/ElInfo.h +++ b/AMDiS/src/ElInfo.h @@ -90,6 +90,8 @@ namespace AMDiS { parent_ = rhs.parent_; macroElement_ = rhs.macroElement_; fillFlag_ = rhs.fillFlag_; + level_ = rhs.level_; + iChild = rhs.iChild; coord_ = rhs.coord_; boundary_ = rhs.boundary_; oppCoord_ = rhs.oppCoord_; @@ -147,6 +149,12 @@ namespace AMDiS { return level_; }; + /** \brief + * Get ElInfo's \ref iChild + */ + inline int getIChild() const { + return iChild; + }; /** \brief * Get ElInfo's \ref coord_[i]. This is a WorldVector<double> filled with the @@ -204,7 +212,7 @@ namespace AMDiS { */ inline Element* getNeighbour(int i) const { return neighbour_[i]; - }; + } /** \brief * Get ElInfo's \ref neighbourCoord_[i] @@ -218,39 +226,43 @@ namespace AMDiS { */ inline unsigned char getOppVertex(int i) const { return oppVertex_[i]; - }; + } virtual int getSideOfNeighbour(int i) { return oppVertex_[i]; - }; + } /** \brief * Get ElInfo's \ref det_ */ inline double getDet() const { return det_; - }; + } /** \brief * Returns \ref grdLambda_ */ inline const DimVec<WorldVector<double> >& getGrdLambda() const { return grdLambda_; - }; + } /** \brief * Returns \ref projection_[i] */ inline Projection *getProjection(int i) const { return projection_[i]; - }; + } /** \brief * Returns \ref parametric_ */ inline bool getParametric() { return parametric_; - }; + } + + inline VectorOfFixVecs<DimVec<double> > *getSubElemCoords() { + return subElemCoords; + } /** \} */ @@ -321,28 +333,32 @@ namespace AMDiS { */ inline void setBoundary(int i, BoundaryType t) { boundary_[i] = newBound(boundary_[i], t); - }; + } /** \brief * Set \ref projection_[i] = p */ inline void setProjection(int i, Projection *p) { projection_[i] = p; - }; + } /** \brief * Set \ref det_ = d */ inline void setDet(double d) { det_ = d; - }; + } /** \brief * Set \ref parametric_ = param */ inline void setParametric(bool param) { parametric_ = param; - }; + } + + inline void setSubElemCoords(VectorOfFixVecs<DimVec<double> > *coords) { + subElemCoords = coords; + } /** \} */ @@ -393,7 +409,7 @@ namespace AMDiS { * pure virtual => must be overriden in sub-class. */ virtual const int worldToCoord(const WorldVector<double>& world, - DimVec<double>* lambda) const = 0; + DimVec<double>* lambda) const = 0; /** \brief * Fills this ElInfo with macro element information of mel. @@ -440,6 +456,10 @@ namespace AMDiS { return(0.0); }; + virtual void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const = 0; + + virtual void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords, + int iChild) const = 0; protected: @@ -477,6 +497,11 @@ namespace AMDiS { */ unsigned char level_; + /** \brief + * This ElInfo is the iChild-th child of the parent element. + */ + int iChild; + /** \brief * \ref coord_[i] is a WorldVector<double> storing the world coordinates of the * i-th vertex of element \ref element. @@ -542,6 +567,8 @@ namespace AMDiS { */ int dimOfWorld; + VectorOfFixVecs<DimVec<double> > *subElemCoords; + // ===== static public members ================================================ public: /** \brief diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc index b4737d41..9559b1a6 100644 --- a/AMDiS/src/ElInfo1d.cc +++ b/AMDiS/src/ElInfo1d.cc @@ -185,29 +185,30 @@ namespace AMDiS { } - void ElInfo1d::fillElInfo(int ichild, const ElInfo *elinfo_old) + void ElInfo1d::fillElInfo(int ichild, const ElInfo *elInfoOld) { FUNCNAME("ElInfo1d::fillElInfo()"); Element *nb; - Element *elem = elinfo_old->element_; + Element *elem = elInfoOld->element_; TEST_EXIT_DBG(elem->getChild(0))("no children?\n"); element_ = const_cast<Element*>(elem->getChild(ichild)); TEST_EXIT_DBG(element_)("missing child %d?\n", ichild); - macroElement_ = elinfo_old->macroElement_; - fillFlag_ = elinfo_old->fillFlag_; + macroElement_ = elInfoOld->macroElement_; + fillFlag_ = elInfoOld->fillFlag_; parent_ = elem; - level_ = elinfo_old->level_ + 1; + level_ = elInfoOld->level_ + 1; + iChild = ichild; int neighbours = mesh_->getGeo(NEIGH); if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) { - const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elinfo_old->coord_); + const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elInfoOld->coord_); coord_[ichild] = (*old_coord)[ichild]; if (elem->getNewCoord(-1)) { @@ -227,13 +228,13 @@ namespace AMDiS { if (i != ichild) { nb = const_cast<Element*>( elem->getChild(1-ichild)); if ( fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { - oppC = elinfo_old->coord_[i]; + oppC = elInfoOld->coord_[i]; } } else { - nb = const_cast<Element*>( elinfo_old->getNeighbour(i)); + nb = const_cast<Element*>( elInfoOld->getNeighbour(i)); if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { - oppC = elinfo_old->oppCoord_[i]; + oppC = elInfoOld->oppCoord_[i]; } } @@ -259,16 +260,36 @@ namespace AMDiS { } if (fillFlag_.isSet(Mesh::FILL_BOUND)) { - boundary_[ichild] = elinfo_old->getBoundary(ichild); - boundary_[1-ichild] = INTERIOR; + boundary_[ichild] = elInfoOld->getBoundary(ichild); + boundary_[1 - ichild] = INTERIOR; - if (elinfo_old->getProjection(0) && - elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) { - projection_[0] = elinfo_old->getProjection(0); + if (elInfoOld->getProjection(0) && + elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) { + projection_[0] = elInfoOld->getProjection(0); } } return; } + void ElInfo1d::getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const + { + (*coords)[0][0] = 0.0; + (*coords)[0][1] = 1.0; + (*coords)[1][0] = 1.0; + (*coords)[1][1] = 0.0; + } + + void ElInfo1d::getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords, + int iChild) const + { + if (iChild == 0) { + (*coords)[1][0] = ((*coords)[0][0] + (*coords)[1][0]) * 0.5; + (*coords)[1][1] = ((*coords)[0][1] + (*coords)[1][1]) * 0.5; + } else { + (*coords)[0][0] = ((*coords)[0][0] + (*coords)[1][0]) * 0.5; + (*coords)[0][1] = ((*coords)[0][1] + (*coords)[1][1]) * 0.5; + } + } + } diff --git a/AMDiS/src/ElInfo1d.h b/AMDiS/src/ElInfo1d.h index 9cb395c0..57b75cbe 100644 --- a/AMDiS/src/ElInfo1d.h +++ b/AMDiS/src/ElInfo1d.h @@ -78,6 +78,12 @@ namespace AMDiS { int getSideOfNeighbour(int i) { return (i + 1) % 2; }; + + void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const; + + void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords, + int iChild) const; + }; } diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc index 062bf6f2..a8ced02b 100644 --- a/AMDiS/src/ElInfo2d.cc +++ b/AMDiS/src/ElInfo2d.cc @@ -155,14 +155,14 @@ namespace AMDiS { /* fill ElInfo structure for one child of an element */ /****************************************************************************/ - void ElInfo2d::fillElInfo(int ichild, const ElInfo *elinfo_old) + void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld) { FUNCNAME("ElInfo::fillElInfo()"); - Element *elem = elinfo_old->element_; + Element *elem = elInfoOld->element_; Element *nb; - Flag fill_flag = elinfo_old->fillFlag_; + Flag fill_flag = elInfoOld->fillFlag_; TEST_EXIT_DBG(elem->getFirstChild())("no children?\n"); element_ = const_cast<Element*>((ichild == 0) ? @@ -170,10 +170,11 @@ namespace AMDiS { elem->getSecondChild()); TEST_EXIT_DBG(element_)("missing child %d?\n", ichild); - macroElement_ = elinfo_old->macroElement_; + macroElement_ = elInfoOld->macroElement_; fillFlag_ = fill_flag; parent_ = elem; - level_ = elinfo_old->level_ + 1; + level_ = elInfoOld->level_ + 1; + iChild = ichild; if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) || @@ -182,15 +183,15 @@ namespace AMDiS { if (elem->getNewCoord(-1)) { coord_[2] = *(elem->getNewCoord()); } else { - coord_[2].setMidpoint(elinfo_old->coord_[0], elinfo_old->coord_[1]); + coord_[2].setMidpoint(elInfoOld->coord_[0], elInfoOld->coord_[1]); } if (ichild == 0) { - coord_[0] = elinfo_old->coord_[2]; - coord_[1] = elinfo_old->coord_[0]; + coord_[0] = elInfoOld->coord_[2]; + coord_[1] = elInfoOld->coord_[0]; } else { - coord_[0] = elinfo_old->coord_[1]; - coord_[1] = elinfo_old->coord_[2]; + coord_[0] = elInfoOld->coord_[1]; + coord_[1] = elInfoOld->coord_[2]; } } @@ -201,12 +202,12 @@ namespace AMDiS { // Calculation of the neighbour 2, its oppCoords and the // cooresponding oppVertex. - neighbour_[2] = elinfo_old->neighbour_[1]; - oppVertex_[2] = elinfo_old->oppVertex_[1]; + neighbour_[2] = elInfoOld->neighbour_[1]; + oppVertex_[2] = elInfoOld->oppVertex_[1]; if (neighbour_[2] && fill_opp_coords) { - oppCoord_[2] = elinfo_old->oppCoord_[1]; - neighbourCoord_[2] = elinfo_old->neighbourCoord_[1]; + oppCoord_[2] = elInfoOld->oppCoord_[1]; + neighbourCoord_[2] = elInfoOld->neighbourCoord_[1]; } @@ -224,8 +225,8 @@ namespace AMDiS { if (elem->getSecondChild()->getNewCoord(-1)) { oppCoord_[1] = *(elem->getSecondChild()->getNewCoord()); } else { - oppCoord_[1].setMidpoint(elinfo_old->coord_[1], - elinfo_old->coord_[2]); + oppCoord_[1].setMidpoint(elInfoOld->coord_[1], + elInfoOld->coord_[2]); } neighbourCoord_[1][0] = coord_[0]; @@ -237,10 +238,10 @@ namespace AMDiS { oppVertex_[1] = 0; if (fill_opp_coords) { - oppCoord_[1] = elinfo_old->coord_[1]; + oppCoord_[1] = elInfoOld->coord_[1]; - neighbourCoord_[1][0] = elinfo_old->coord_[1]; - neighbourCoord_[1][1] = elinfo_old->coord_[2]; + neighbourCoord_[1][0] = elInfoOld->coord_[1]; + neighbourCoord_[1][1] = elInfoOld->coord_[2]; neighbourCoord_[1][2] = coord_[2]; } } @@ -249,9 +250,9 @@ namespace AMDiS { // Calculation of the neighbour 0, its oppCoords and the // cooresponding oppVertex. - nb = elinfo_old->neighbour_[2]; + nb = elInfoOld->neighbour_[2]; if (nb) { - TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n"); + TEST(elInfoOld->oppVertex_[2] == 2)("invalid neighbour\n"); TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n"); TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n"); @@ -264,13 +265,13 @@ namespace AMDiS { if (nb->getNewCoord(-1)) { oppCoord_[0] = *(nb->getNewCoord()); } else { - oppCoord_[0].setMidpoint(elinfo_old->neighbourCoord_[2][1], - elinfo_old->neighbourCoord_[2][2]); + oppCoord_[0].setMidpoint(elInfoOld->neighbourCoord_[2][1], + elInfoOld->neighbourCoord_[2][2]); } - neighbourCoord_[0][0].setMidpoint(elinfo_old->neighbourCoord_[2][0], - elinfo_old->neighbourCoord_[2][1]); - neighbourCoord_[0][1] = elinfo_old->neighbourCoord_[2][1]; + neighbourCoord_[0][0].setMidpoint(elInfoOld->neighbourCoord_[2][0], + elInfoOld->neighbourCoord_[2][1]); + neighbourCoord_[0][1] = elInfoOld->neighbourCoord_[2][1]; neighbourCoord_[0][2] = oppCoord_[0]; } @@ -279,12 +280,12 @@ namespace AMDiS { oppVertex_[0] = 1; if (fill_opp_coords) { - oppCoord_[0] = elinfo_old->oppCoord_[2]; + oppCoord_[0] = elInfoOld->oppCoord_[2]; - neighbourCoord_[0][0] = elinfo_old->neighbourCoord_[2][0]; - neighbourCoord_[0][1] = elinfo_old->neighbourCoord_[2][2]; - neighbourCoord_[0][2].setMidpoint(elinfo_old->neighbourCoord_[2][0], - elinfo_old->neighbourCoord_[2][1]); + neighbourCoord_[0][0] = elInfoOld->neighbourCoord_[2][0]; + neighbourCoord_[0][1] = elInfoOld->neighbourCoord_[2][2]; + neighbourCoord_[0][2].setMidpoint(elInfoOld->neighbourCoord_[2][0], + elInfoOld->neighbourCoord_[2][1]); } } } @@ -294,12 +295,12 @@ namespace AMDiS { // Calculation of the neighbour 2, its oppCoords and the // cooresponding oppVertex. - neighbour_[2] = elinfo_old->neighbour_[0]; - oppVertex_[2] = elinfo_old->oppVertex_[0]; + neighbour_[2] = elInfoOld->neighbour_[0]; + oppVertex_[2] = elInfoOld->oppVertex_[0]; if (neighbour_[2] && fill_opp_coords) { - oppCoord_[2] = elinfo_old->oppCoord_[0]; - neighbourCoord_[2] = elinfo_old->neighbourCoord_[0]; + oppCoord_[2] = elInfoOld->oppCoord_[0]; + neighbourCoord_[2] = elInfoOld->neighbourCoord_[0]; } @@ -314,8 +315,8 @@ namespace AMDiS { if (elem->getFirstChild()->getNewCoord(-1)) { oppCoord_[0] = *(elem->getFirstChild()->getNewCoord()); } else { - oppCoord_[0].setMidpoint(elinfo_old->coord_[0], - elinfo_old->coord_[2]); + oppCoord_[0].setMidpoint(elInfoOld->coord_[0], + elInfoOld->coord_[2]); } neighbourCoord_[0][0] = coord_[2]; @@ -327,10 +328,10 @@ namespace AMDiS { oppVertex_[0] = 1; if (fill_opp_coords) { - oppCoord_[0] = elinfo_old->coord_[0]; + oppCoord_[0] = elInfoOld->coord_[0]; - neighbourCoord_[0][0] = elinfo_old->coord_[2]; - neighbourCoord_[0][1] = elinfo_old->coord_[0]; + neighbourCoord_[0][0] = elInfoOld->coord_[2]; + neighbourCoord_[0][1] = elInfoOld->coord_[0]; neighbourCoord_[0][2] = coord_[2]; } } @@ -338,9 +339,9 @@ namespace AMDiS { // Calculation of the neighbour 1, its oppCoords and the // cooresponding oppVertex. - nb = elinfo_old->neighbour_[2]; + nb = elInfoOld->neighbour_[2]; if (nb) { - TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n"); + TEST(elInfoOld->oppVertex_[2] == 2)("invalid neighbour\n"); TEST((nb = nb->getFirstChild()))("missing child?\n"); if (nb->getFirstChild()) { @@ -350,13 +351,13 @@ namespace AMDiS { if (nb->getNewCoord(-1)) { oppCoord_[1] = *(nb->getNewCoord()); } else { - oppCoord_[1].setMidpoint(elinfo_old->neighbourCoord_[2][0], - elinfo_old->neighbourCoord_[2][2]); + oppCoord_[1].setMidpoint(elInfoOld->neighbourCoord_[2][0], + elInfoOld->neighbourCoord_[2][2]); } - neighbourCoord_[1][0] = elinfo_old->neighbourCoord_[2][0]; - neighbourCoord_[1][1].setMidpoint(elinfo_old->neighbourCoord_[2][0], - elinfo_old->neighbourCoord_[2][1]); + neighbourCoord_[1][0] = elInfoOld->neighbourCoord_[2][0]; + neighbourCoord_[1][1].setMidpoint(elInfoOld->neighbourCoord_[2][0], + elInfoOld->neighbourCoord_[2][1]); neighbourCoord_[1][2] = oppCoord_[1]; } nb = nb->getSecondChild(); @@ -365,12 +366,12 @@ namespace AMDiS { oppVertex_[1] = 0; if (fill_opp_coords) { - oppCoord_[1] = elinfo_old->oppCoord_[2]; + oppCoord_[1] = elInfoOld->oppCoord_[2]; - neighbourCoord_[1][0] = elinfo_old->neighbourCoord_[2][2]; - neighbourCoord_[1][1] = elinfo_old->neighbourCoord_[2][0]; - neighbourCoord_[1][2].setMidpoint(elinfo_old->neighbourCoord_[2][0], - elinfo_old->neighbourCoord_[2][1]); + neighbourCoord_[1][0] = elInfoOld->neighbourCoord_[2][2]; + neighbourCoord_[1][1] = elInfoOld->neighbourCoord_[2][0]; + neighbourCoord_[1][2].setMidpoint(elInfoOld->neighbourCoord_[2][0], + elInfoOld->neighbourCoord_[2][1]); } } } @@ -380,39 +381,39 @@ namespace AMDiS { if (fill_flag.isSet(Mesh::FILL_BOUND)) { - if (elinfo_old->getBoundary(2)) { - boundary_[5] = elinfo_old->getBoundary(2); + if (elInfoOld->getBoundary(2)) { + boundary_[5] = elInfoOld->getBoundary(2); } else { boundary_[5] = INTERIOR; } if (ichild == 0) { - boundary_[3] = elinfo_old->getBoundary(5); - boundary_[4] = elinfo_old->getBoundary(3); - boundary_[0] = elinfo_old->getBoundary(2); + boundary_[3] = elInfoOld->getBoundary(5); + boundary_[4] = elInfoOld->getBoundary(3); + boundary_[0] = elInfoOld->getBoundary(2); boundary_[1] = INTERIOR; - boundary_[2] = elinfo_old->getBoundary(1); + boundary_[2] = elInfoOld->getBoundary(1); } else { - boundary_[3] = elinfo_old->getBoundary(4); - boundary_[4] = elinfo_old->getBoundary(5); + boundary_[3] = elInfoOld->getBoundary(4); + boundary_[4] = elInfoOld->getBoundary(5); boundary_[0] = INTERIOR; - boundary_[1] = elinfo_old->getBoundary(2); - boundary_[2] = elinfo_old->getBoundary(0); + boundary_[1] = elInfoOld->getBoundary(2); + boundary_[2] = elInfoOld->getBoundary(0); } - if (elinfo_old->getProjection(0) && - elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) { + if (elInfoOld->getProjection(0) && + elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) { - projection_[0] = elinfo_old->getProjection(0); + projection_[0] = elInfoOld->getProjection(0); } else { // boundary projection if (ichild == 0) { - projection_[0] = elinfo_old->getProjection(2); + projection_[0] = elInfoOld->getProjection(2); projection_[1] = NULL; - projection_[2] = elinfo_old->getProjection(1); + projection_[2] = elInfoOld->getProjection(1); } else { projection_[0] = NULL; - projection_[1] = elinfo_old->getProjection(2); - projection_[2] = elinfo_old->getProjection(0); + projection_[1] = elInfoOld->getProjection(2); + projection_[2] = elInfoOld->getProjection(0); } } } diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h index 67b7291d..48ca680b 100644 --- a/AMDiS/src/ElInfo2d.h +++ b/AMDiS/src/ElInfo2d.h @@ -77,6 +77,21 @@ namespace AMDiS { */ double getElementNormal( WorldVector<double> &normal) const; + void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const + { + FUNCNAME("ElInfo2d::getRefSimplexCoords()"); + + ERROR_EXIT("NOT YET\n"); + } + + void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords, + int iChild) const + { + FUNCNAME("ElInfo2d::getSubElementCoords()"); + + ERROR_EXIT("NOT YET\n"); + } + protected: /** \brief * Temp vectors for function \ref calcGrdLambda. diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc index d92af234..2e41a93e 100644 --- a/AMDiS/src/ElInfo3d.cc +++ b/AMDiS/src/ElInfo3d.cc @@ -351,7 +351,7 @@ namespace AMDiS { - void ElInfo3d::fillElInfo(int ichild, const ElInfo *elinfo_old) + void ElInfo3d::fillElInfo(int ichild, const ElInfo *elInfoOld) { FUNCNAME("ElInfo3d::fillElInfo()"); @@ -360,20 +360,21 @@ namespace AMDiS { const int (*cvg)[4] = NULL; /* cvg = child_vertex[el_type] */ int *ce; /* ce = child_edge[el_type][ichild] */ Element *nb, *nbk; - Element *el_old = elinfo_old->element_; - Flag fillFlag__local = elinfo_old->fillFlag_; + Element *el_old = elInfoOld->element_; + Flag fillFlag__local = elInfoOld->fillFlag_; DegreeOfFreedom *dof; int ov = -1; - Mesh *mesh = elinfo_old->getMesh(); + Mesh *mesh = elInfoOld->getMesh(); TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n"); element_ = const_cast<Element*>( el_old->getChild(ichild)); - macroElement_ = elinfo_old->macroElement_; + macroElement_ = elInfoOld->macroElement_; fillFlag_ = fillFlag__local; parent_ = el_old; - level_ = elinfo_old->level_ + 1; - int el_type_local = (dynamic_cast<const ElInfo3d*>(elinfo_old))->getType(); + level_ = elInfoOld->level_ + 1; + iChild = ichild; + int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType(); el_type = (el_type_local + 1) % 3; TEST_EXIT_DBG(element_)("missing child %d?\n", ichild); @@ -389,14 +390,14 @@ namespace AMDiS { fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) { for (int i = 0; i < 3; i++) { for (int j = 0; j < dimOfWorld; j++) { - coord_[i][j] = elinfo_old->coord_[cv[i]][j]; + coord_[i][j] = elInfoOld->coord_[cv[i]][j]; } } if (el_old->getNewCoord()) { coord_[3] = *(el_old->getNewCoord()); } else { for (int j = 0; j < dimOfWorld; j++) { - coord_[3][j] = (elinfo_old->coord_[0][j] + elinfo_old->coord_[1][j]) / 2; + coord_[3][j] = (elInfoOld->coord_[0][j] + elInfoOld->coord_[1][j]) / 2; } } } @@ -405,7 +406,7 @@ namespace AMDiS { fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { FixVec<Element*, NEIGH> *neigh_local = &neighbour_; - const FixVec<Element*, NEIGH> *neigh_old = &elinfo_old->neighbour_; + const FixVec<Element*, NEIGH> *neigh_old = &elInfoOld->neighbour_; Flag fill_opp_coords; fill_opp_coords.setFlags(fillFlag__local & Mesh::FILL_OPP_COORDS); @@ -423,7 +424,7 @@ namespace AMDiS { } else { int k = cvg[ochild][1]; for (int j = 0; j < dimOfWorld; j++) { - oppCoord_[0][j] = (elinfo_old->coord_[ochild][j] + elinfo_old->coord_[k][j]) / 2; + oppCoord_[0][j] = (elInfoOld->coord_[ochild][j] + elInfoOld->coord_[k][j]) / 2; } } } @@ -432,7 +433,7 @@ namespace AMDiS { } else { if (fill_opp_coords.isAnySet()) { for (int j = 0; j < dimOfWorld; j++) { - oppCoord_[0][j] = elinfo_old->coord_[ochild][j]; + oppCoord_[0][j] = elInfoOld->coord_[ochild][j]; } } (*neigh_local)[0] = nb; @@ -455,7 +456,7 @@ namespace AMDiS { nbk = const_cast<Element*>( nb->getChild(k)); if (nbk->getDOF(0) == el_old->getDOF(ichild)) { /* opp. vertex */ - dof = const_cast<int*>(nb->getDOF(elinfo_old->oppVertex_[cv[i]])); + dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex_[cv[i]])); if (dof == nbk->getDOF(1)) { ov = 1; @@ -465,7 +466,7 @@ namespace AMDiS { oppCoord_[i] = *(nbk->getNewCoord()); } else { for (int j = 0; j < dimOfWorld; j++) { - oppCoord_[i][j] = (elinfo_old->oppCoord_[cv[i]][j] + elinfo_old->coord_[ichild][j]) / 2; + oppCoord_[i][j] = (elInfoOld->oppCoord_[cv[i]][j] + elInfoOld->coord_[ichild][j]) / 2; } } } @@ -483,7 +484,7 @@ namespace AMDiS { if (fill_opp_coords.isAnySet()) { for (int j = 0; j < dimOfWorld; j++) { - oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j]; + oppCoord_[i][j] = elInfoOld->oppCoord_[cv[i]][j]; } } (*neigh_local)[i] = nbk; @@ -500,7 +501,7 @@ namespace AMDiS { if (nbk->getDOF(0) == el_old->getDOF(ichild) || mesh->associated(nbk->getDOF(0, 0), el_old->getDOF(ichild, 0))) { /* opp. vertex */ - dof = const_cast<int*>(nb->getDOF(elinfo_old->oppVertex_[cv[i]])); + dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex_[cv[i]])); if (dof == nbk->getDOF(1) || mesh->associated(dof[0], nbk->getDOF(1, 0))) { @@ -511,7 +512,7 @@ namespace AMDiS { oppCoord_[i] = *(nbk->getNewCoord()); } else { for (int j = 0; j < dimOfWorld; j++) { - oppCoord_[i][j] = (elinfo_old->oppCoord_[cv[i]][j] + elinfo_old->coord_[ichild][j]) / 2; + oppCoord_[i][j] = (elInfoOld->oppCoord_[cv[i]][j] + elInfoOld->coord_[ichild][j]) / 2; } } } @@ -528,7 +529,7 @@ namespace AMDiS { if (fill_opp_coords.isAnySet()) { for (int j = 0; j < dimOfWorld; j++) { - oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j]; + oppCoord_[i][j] = elInfoOld->oppCoord_[cv[i]][j]; } } (*neigh_local)[i] = nbk; @@ -548,11 +549,11 @@ namespace AMDiS { /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/ if (((*neigh_local)[3] = (*neigh_old)[ochild])) { - oppVertex_[3] = elinfo_old->oppVertex_[ochild]; + oppVertex_[3] = elInfoOld->oppVertex_[ochild]; if (fill_opp_coords.isAnySet()) { for (int j = 0; j < dimOfWorld; j++) { - oppCoord_[3][j] = elinfo_old->oppCoord_[ochild][j]; + oppCoord_[3][j] = elInfoOld->oppCoord_[ochild][j]; } } } @@ -560,42 +561,42 @@ namespace AMDiS { if (fillFlag__local.isSet(Mesh::FILL_BOUND)) { for (int i = 0; i < 3; i++) { - boundary_[10 + i] = elinfo_old->getBoundary(10 + cv[i]); + boundary_[10 + i] = elInfoOld->getBoundary(10 + cv[i]); } - boundary_[13] = elinfo_old->getBoundary(4); + boundary_[13] = elInfoOld->getBoundary(4); boundary_[0] = INTERIOR; - boundary_[1] = elinfo_old->getBoundary(cv[1]); - boundary_[2] = elinfo_old->getBoundary(cv[2]); - boundary_[3] = elinfo_old->getBoundary(ochild); + boundary_[1] = elInfoOld->getBoundary(cv[1]); + boundary_[2] = elInfoOld->getBoundary(cv[2]); + boundary_[3] = elInfoOld->getBoundary(ochild); int geoFace = mesh_->getGeo(FACE); ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]); for (int iedge = 0; iedge < 4; iedge++) { - boundary_[geoFace + iedge] = elinfo_old->getBoundary(geoFace + ce[iedge]); + boundary_[geoFace + iedge] = elInfoOld->getBoundary(geoFace + ce[iedge]); } for (int iedge = 4; iedge < 6; iedge++) { int i = 5 - cv[iedge - 3]; /* old vertex opposite new edge */ - boundary_[geoFace + iedge] = elinfo_old->getBoundary(i); + boundary_[geoFace + iedge] = elInfoOld->getBoundary(i); } - if (elinfo_old->getProjection(0) && - elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) { + if (elInfoOld->getProjection(0) && + elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) { - projection_[0] = elinfo_old->getProjection(0); + projection_[0] = elInfoOld->getProjection(0); } else { // boundary projection projection_[0] = NULL; - projection_[1] = elinfo_old->getProjection(cv[1]); - projection_[2] = elinfo_old->getProjection(cv[2]); - projection_[3] = elinfo_old->getProjection(ochild); + projection_[1] = elInfoOld->getProjection(cv[1]); + projection_[2] = elInfoOld->getProjection(cv[2]); + projection_[3] = elInfoOld->getProjection(ochild); for (int iedge = 0; iedge < 4; iedge++) { - projection_[geoFace + iedge] = elinfo_old->getProjection(geoFace + ce[iedge]); + projection_[geoFace + iedge] = elInfoOld->getProjection(geoFace + ce[iedge]); } for (int iedge = 4; iedge < 6; iedge++) { - projection_[geoFace + iedge] = elinfo_old->getProjection(5 - cv[iedge - 3]); + projection_[geoFace + iedge] = elInfoOld->getProjection(5 - cv[iedge - 3]); } } } @@ -603,7 +604,7 @@ namespace AMDiS { if (fillFlag_.isSet(Mesh::FILL_ORIENTATION)) { orientation = - ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elinfo_old)))->orientation + ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elInfoOld)))->orientation * Tetrahedron::childOrientation[el_type_local][ichild]; } } diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h index 3c12ad81..1be858be 100644 --- a/AMDiS/src/ElInfo3d.h +++ b/AMDiS/src/ElInfo3d.h @@ -119,6 +119,21 @@ namespace AMDiS { orientation = o; }; + void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const + { + FUNCNAME("ElInfo3d::getRefSimplexCoords()"); + + ERROR_EXIT("NOT YET\n"); + } + + void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords, + int iChild) const + { + FUNCNAME("ElInfo3d::getSubElementCoords()"); + + ERROR_EXIT("NOT YET\n"); + } + protected: /** \brief diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h index ac7e312f..975f2e70 100644 --- a/AMDiS/src/FixVec.h +++ b/AMDiS/src/FixVec.h @@ -177,7 +177,7 @@ namespace AMDiS { TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n"); vec = GET_MEMORY(FixVecType*, size_); - for (FixVecType** i=&vec[0]; i < &vec[size_]; i++) { + for (FixVecType** i = &vec[0]; i < &vec[size_]; i++) { *i = NEW FixVecType(dim, NO_INIT); } } @@ -193,7 +193,7 @@ namespace AMDiS { TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n"); vec = GET_MEMORY(FixVecType*, size_); - for (FixVecType** i=&vec[0]; i < &vec[size_]; i++) { + for (FixVecType** i = &vec[0]; i < &vec[size_]; i++) { *i = NEW FixVecType(ini[i]); } } @@ -209,8 +209,7 @@ namespace AMDiS { TEST_EXIT_DBG(initType == DEFAULT_VALUE) ("wrong initType or wrong initializer\n"); vec = GET_MEMORY(FixVecType*, size_); - //vec = new FixVecType*[size]; - for(FixVecType** i=&vec[0]; i < &vec[size_]; i++) { + for (FixVecType** i = &vec[0]; i < &vec[size_]; i++) { *i = NEW FixVecType(ini); } } @@ -233,7 +232,7 @@ namespace AMDiS { */ virtual ~VectorOfFixVecs() { - for(FixVecType** i=&vec[0]; i<&vec[size_]; i++) { + for (FixVecType** i = &vec[0]; i<&vec[size_]; i++) { DELETE *i; } FREE_MEMORY(vec, FixVecType*, size_); diff --git a/AMDiS/src/Flag.cc b/AMDiS/src/Flag.cc deleted file mode 100644 index bd4afe36..00000000 --- a/AMDiS/src/Flag.cc +++ /dev/null @@ -1,5 +0,0 @@ -#include "Flag.h" - -namespace AMDiS { - -} diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index 475b13d2..9ade143d 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -43,9 +43,9 @@ namespace AMDiS { const Flag Mesh::FILL_ADD_ALL = 0X80L; - const Flag Mesh::FILL_ANY_1D= (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L); - const Flag Mesh::FILL_ANY_2D= (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L); - const Flag Mesh::FILL_ANY_3D= (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L); + const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L); + const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L); + const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L); //************************************************************************** // flags for Mesh traversal @@ -62,17 +62,16 @@ namespace AMDiS { // const Flag Mesh::USE_PARAMETRIC = 0X8000L ; // used in mg methods - // std::list<Mesh*> Mesh::meshes; DOFAdmin* Mesh::compressAdmin = NULL; Mesh* Mesh::traversePtr = NULL; int Mesh::iadmin = 0; std::vector<DegreeOfFreedom> Mesh::dof_used; - const int Mesh::MAX_DOF=100; + const int Mesh::MAX_DOF = 100; std::map<DegreeOfFreedom, DegreeOfFreedom*> Mesh::serializedDOFs; struct delmem { DegreeOfFreedom* ptr; - int len; + int len; }; diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index 4fc8beb9..aa3bb993 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/AMDiS/src/ProblemVec.cc @@ -994,6 +994,7 @@ namespace AMDiS { ElInfo *rowElInfo, *colElInfo; ElInfo *largeElInfo, *smallElInfo; + dualTraverse.setFillSubElemInfo(true); bool cont = dualTraverse.traverseFirst(rowMesh, colMesh, -1, -1, assembleFlag, assembleFlag, &rowElInfo, &colElInfo, diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc index e82a8c58..b54d11d5 100644 --- a/AMDiS/src/Traverse.cc +++ b/AMDiS/src/Traverse.cc @@ -47,7 +47,7 @@ namespace AMDiS { ElInfo* TraverseStack::traverseNext(ElInfo* elinfo_old) { FUNCNAME("TraverseStack::traverseNext()"); -\ + ElInfo *elinfo = NULL; ElInfo::traverseId = id; Parametric *parametric = traverse_mesh->getParametric(); @@ -112,7 +112,6 @@ namespace AMDiS { /* depending on the traverse_info->level variable */ /****************************************************************************/ - // int Traverse::recursive(ElInfo *elinfo) int Traverse::recursive(ElInfoStack *elInfoStack) { FUNCNAME("Traverse::recursive()"); @@ -1071,4 +1070,18 @@ namespace AMDiS { } } + void TraverseStack::getCoordsInElem(const ElInfo *upperElInfo, + VectorOfFixVecs<DimVec<double> > *coords) + { + int upperLevel = static_cast<int>(upperElInfo->getLevel()); + int lowerLevel = static_cast<int>(elinfo_stack[stack_used]->getLevel()); + int stackPos = stack_used; + upperElInfo->getRefSimplexCoords(coords); + + while (lowerLevel > upperLevel) { + upperElInfo->getSubElementCoords(coords, elinfo_stack[stackPos]->getIChild()); + stackPos--; + lowerLevel--; + } + } } diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h index ba18bf0f..32a2fc26 100644 --- a/AMDiS/src/Traverse.h +++ b/AMDiS/src/Traverse.h @@ -143,6 +143,9 @@ namespace AMDiS { */ void update(); + void getCoordsInElem(const ElInfo *upperElInfo, + VectorOfFixVecs<DimVec<double> > *coords); + /** \brief * */ -- GitLab