From d906fc98b49718d333133bd2c11b7aa17d5baf84 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski Date: Mon, 23 Nov 2009 16:48:11 +0000 Subject: [PATCH] Code refactoring, so added some no bugs. --- AMDiS/src/CoarseningManager.cc | 83 ++-- AMDiS/src/CoarseningManager.h | 19 - AMDiS/src/Error.h | 21 - AMDiS/src/Error.hh | 340 +++++++-------- AMDiS/src/MacroReader.cc | 110 +++-- AMDiS/src/MacroReader.h | 4 +- AMDiS/src/Mesh.cc | 49 +-- AMDiS/src/Mesh.h | 22 - AMDiS/src/ProblemVec.cc | 3 + AMDiS/src/ProblemVecDbg.cc | 42 +- AMDiS/src/ProblemVecDbg.h | 1 + AMDiS/src/RefinementManager.cc | 29 +- AMDiS/src/RefinementManager.h | 14 - AMDiS/src/RefinementManager1d.cc | 44 +- AMDiS/src/RefinementManager1d.h | 4 +- AMDiS/src/RefinementManager2d.cc | 14 +- AMDiS/src/RefinementManager2d.h | 2 +- AMDiS/src/RefinementManager3d.cc | 683 ++++++++++++++----------------- AMDiS/src/TecPlotWriter.hh | 36 +- AMDiS/src/Traverse.cc | 206 +--------- AMDiS/src/Traverse.h | 3 - AMDiS/src/TraverseParallel.cc | 23 +- AMDiS/src/TraverseParallel.h | 10 +- 23 files changed, 708 insertions(+), 1054 deletions(-) diff --git a/AMDiS/src/CoarseningManager.cc b/AMDiS/src/CoarseningManager.cc index ed5dd54e..51f34507 100644 --- a/AMDiS/src/CoarseningManager.cc +++ b/AMDiS/src/CoarseningManager.cc @@ -10,57 +10,44 @@ namespace AMDiS { - CoarseningManager* CoarseningManager::traversePtr = NULL; - - /****************************************************************************/ - /* sets the mark on all elements that have to be coarsend */ - /****************************************************************************/ - - int CoarseningManager::coarsenMarkFunction(ElInfo *el_info) - { - el_info->getElement()->setMark(traversePtr->globalMark); - return 0; - } - /****************************************************************************/ /* tries to coarsen every element of mesh at least mark times */ /****************************************************************************/ Flag CoarseningManager::globalCoarsen(Mesh *aMesh, int mark) { - if (mark >= 0) return(0); - globalMark = mark; - traversePtr = this; - aMesh->traverse(-1, Mesh::CALL_LEAF_EL, coarsenMarkFunction); - return(coarsenMesh(aMesh)); - } - - int CoarseningManager::spreadCoarsenMarkFunction(ElInfo* el_info) - { - Element *el = el_info->getElement(); - signed char mark; - - if (el->getChild(0)) { - /****************************************************************************/ - /* interior node of the tree */ - /****************************************************************************/ - mark = max(el->getChild(0)->getMark(), el->getChild(1)->getMark()); - el->setMark(std::min(mark + 1, 0)); - } else { - /****************************************************************************/ - /* leaf node of the tree */ - /****************************************************************************/ - if (el->getMark() < 0) - el->setMark(el->getMark() - 1); + if (mark >= 0) + return 0; + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(aMesh, -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + elInfo->getElement()->setMark(mark); + elInfo = stack.traverseNext(elInfo); } - - return 0; + + return coarsenMesh(aMesh); } void CoarseningManager::spreadCoarsenMark() { - traversePtr = this; - mesh->traverse(-1, Mesh::CALL_EVERY_EL_POSTORDER, spreadCoarsenMarkFunction); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_POSTORDER); + while (elInfo) { + Element *el = elInfo->getElement(); + + if (el->getChild(0)) { + // interior node of the tree + signed char mark = max(el->getChild(0)->getMark(), el->getChild(1)->getMark()); + el->setMark(std::min(mark + 1, 0)); + } else { + // leaf node of the tree + if (el->getMark() < 0) + el->setMark(el->getMark() - 1); + } + + elInfo = stack.traverseNext(elInfo); + } } @@ -69,17 +56,15 @@ namespace AMDiS { /* resets the element marks */ /****************************************************************************/ - int CoarseningManager::cleanUpAfterCoarsenFunction(ElInfo *el_info) - { - Element *el = el_info->getElement(); - el->setMark(max(el->getMark(), 0)); - return 0; - } - void CoarseningManager::cleanUpAfterCoarsen() { - traversePtr = this; - mesh->traverse(-1, Mesh::CALL_LEAF_EL, cleanUpAfterCoarsenFunction); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + Element *el = elInfo->getElement(); + el->setMark(max(el->getMark(), 0)); + elInfo = stack.traverseNext(elInfo); + } } /****************************************************************************/ diff --git a/AMDiS/src/CoarseningManager.h b/AMDiS/src/CoarseningManager.h index 94ab64e1..f1103401 100644 --- a/AMDiS/src/CoarseningManager.h +++ b/AMDiS/src/CoarseningManager.h @@ -42,7 +42,6 @@ namespace AMDiS { CoarseningManager() : mesh(NULL), stack(NULL), - globalMark(0), doMore(0) {} @@ -92,15 +91,6 @@ namespace AMDiS { */ void spreadCoarsenMark(); - /// Used while traversal in spreadCoarsenMark - static int spreadCoarsenMarkFunction(ElInfo* el_info); - - /// Used while traversal in cleanUpAfterCoarsen - static int cleanUpAfterCoarsenFunction(ElInfo* el_info); - - /// Sets the mark on all elements that have to be coarsend - static int coarsenMarkFunction(ElInfo *el_info); - /// Resets the element marks void cleanUpAfterCoarsen(); @@ -111,18 +101,9 @@ namespace AMDiS { /// Used for non recursive mesh traversal. TraverseStack *stack; - /// Used by globalCoarsen to remember the given mark value - int globalMark; - /// Spezifies whether the coarsening operation is still in progress bool doMore; - /** \brief - * Used while mesh traversal to have a pointer to this CoarseningManager - * from a static method - */ - static CoarseningManager* traversePtr; - /// Spezifies how many DOFVectors should restricted while coarsening int callCoarseRestrict; diff --git a/AMDiS/src/Error.h b/AMDiS/src/Error.h index e92a0656..358662e7 100644 --- a/AMDiS/src/Error.h +++ b/AMDiS/src/Error.h @@ -68,13 +68,6 @@ namespace AMDiS { double* max, bool writeLeafData = false, int comp = 0); - - // methods for traversal - static int maxErrAtQpFct(ElInfo* elinfo); - static int H1ErrFct(ElInfo* elinfo); - static int L2ErrFct(ElInfo* elinfo); - static int relFct(ElInfo* elinfo); - public: static T errUFct(const DimVec& lambda); @@ -114,13 +107,6 @@ namespace AMDiS { static const AbstractFunction, WorldVector >* pGrdU; static const BasisFunction* basFct; static const DOFVector* errUh; - static double maxErr; - static double l2Err2; - static double l2Norm2; - static int relative; - static double relNorm2; - static double h1Err2; - static double h1Norm2; static bool writeInLeafData; static int component; }; @@ -131,13 +117,6 @@ namespace AMDiS { template const AbstractFunction, WorldVector >* Error::pGrdU = NULL; template const BasisFunction* Error::basFct = NULL; template const DOFVector* Error::errUh = NULL; - template double Error::maxErr = 0.0; - template double Error::l2Err2 = 0.0; - template double Error::l2Norm2 = 0.0; - template int Error::relative = 0; - template double Error::relNorm2 = 0.0; - template double Error::h1Err2 = 0.0; - template double Error::h1Norm2 = 0.0; template typename Error::AbstrFctErrU Error::errU; template typename Error::AbstrFctGrdErrU Error::grdErrU; template bool Error::writeInLeafData = false; diff --git a/AMDiS/src/Error.hh b/AMDiS/src/Error.hh index 3f0a2b38..d3d78896 100644 --- a/AMDiS/src/Error.hh +++ b/AMDiS/src/Error.hh @@ -1,6 +1,7 @@ #include "Mesh.h" #include "Parametric.h" #include "Quadrature.h" +#include "Traverse.h" namespace AMDiS { @@ -20,150 +21,80 @@ namespace AMDiS { return ((*pGrdU)(x)); } - - template - int Error::maxErrAtQpFct(ElInfo* el_info) - { - double err = 0.0; - const double *u_vec, *uh_vec; - - elinfo = el_info; - u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL); - uh_vec = errUh->getVecAtQPs(el_info, - NULL, - quadFast, - NULL); - - int numPoints = quadFast->getNumPoints(); - - for (int i = 0; i < numPoints; i++) { - err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i]; - maxErr = max(maxErr, err); - } - - return; - } - template double Error::maxErrAtQp(const AbstractFunction >& u, const DOFVector& uh, const Quadrature* q) { FUNCNAME("Error::maxErrAtQp()"); - const FiniteElemSpace *fe_space; + const FiniteElemSpace *fe_space; if (!(pU = &u)) { ERROR("no function u specified; doing nothing\n"); return(-1.0); } - - if (!(errUh = &uh) || !(fe_space = uh->getFESpace())) - { - ERROR("no discrete function or no fe_space for it; doing nothing\n"); - return(-1.0); - } - - if (!(basFct = fe_space->getBasisFcts())) - { - ERROR("no basis functions at discrete solution ; doing nothing\n"); - return(-1.0); - } - - int dim = fe_space->getMesh()->getDim(); + if (!(errUh = &uh) || !(fe_space = uh->getFESpace())) { + ERROR("no discrete function or no fe_space for it; doing nothing\n"); + return(-1.0); + } + if (!(basFct = fe_space->getBasisFcts())) { + ERROR("no basis functions at discrete solution ; doing nothing\n"); + return(-1.0); + } if (!q) - q = Quadrature::provideQuadrature(dim, - 2*fe_space->getBasisFcts()->getDegree() - - 2); - quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI); + q = Quadrature::provideQuadrature(fe_space->getMesh()->getDim(), + 2 * fe_space->getBasisFcts()->getDegree() - 2); - maxErr = 0.0; - - fe_space->getMesh()->traverse(-1, - Mesh::FILL_COORDS | - Mesh::CALL_LEAF_EL, - maxErrAtQpFct); - - return(maxErr); - } - - template - int Error::H1ErrFct(ElInfo* el_info) - { - int i, j; - double err, err_2, h1_err_el, norm_el, norm2, det, exact; - const WorldVector *grdu_vec, *grduh_vec; - - elinfo = el_info; - grdu_vec = quadFast->getQuadrature()->grdFAtQp(grdErrU, NULL); - det = el_info->getDet(); - grduh_vec = errUh->getGrdAtQPs(elinfo, NULL, quadFast, NULL); - - int numPoints = quadFast->getNumPoints(); - int dow = Global::getGeo(WORLD); - - for (h1_err_el = i = 0; i < numPoints; i++) - { - for (err_2 = j = 0; j < dow; j++) - { - err = grdu_vec[i][j] - grduh_vec[i][j]; - err_2 += sqr(err); - } - h1_err_el += quadFast->getWeight(i)*err_2; + quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI); + double maxErr = 0.0; + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1, + Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL); + while (elInfo) { + double err = 0.0; + const double *u_vec, *uh_vec; + + u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL); + uh_vec = errUh->getVecAtQPs(elInfo, NULL, quadFast, NULL); + + int nPoints = quadFast->getNumPoints(); + for (int i = 0; i < nPoints; i++) { + err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i]; + maxErr = max(maxErr, err); } - exact = det*h1_err_el; - h1Err2 += exact; - maxErr = max(maxErr, exact); - - if (writeInLeafData) - el_info->getElement()->setEstimation(exact, component); - - if (relative) { - for (norm_el = i = 0; i < numPoints; i++) { - for (norm2 = j = 0; j < dow; j++) - norm2 += sqr(grdu_vec[i][j]); - norm_el += quadFast->getWeight(i)*norm2; - } - h1Norm2 += det*norm_el; + elInfo = stack.traverseNext(elInfo); } - return 0; + return maxErr; } template - double Error::H1Err( - const AbstractFunction, WorldVector >& grdU, - const DOFVector& uh, - int relErr, - double* max, + double Error::H1Err(const AbstractFunction, WorldVector > &grdU, + const DOFVector &uh, + int relErr, + double* max, bool writeLeafData, int comp) { - FUNCNAME("Error::H1Err"); + FUNCNAME("Error::H1Err()"); + const FiniteElemSpace *fe_space; - writeInLeafData = writeLeafData; - component = comp; - Quadrature *q = NULL; - pGrdU = &grdU; - errUh = &uh; - if(!(fe_space = uh.getFESpace())) - { - ERROR("no fe_space for uh; doing nothing\n"); - return(0.0); - } - - if (!(basFct = fe_space->getBasisFcts())) - { - ERROR("no basis functions at discrete solution ; doing nothing\n"); - return(0.0); - } + if (!(fe_space = uh.getFESpace())) { + ERROR("no fe_space for uh; doing nothing\n"); + return(0.0); + } + if (!(basFct = fe_space->getBasisFcts())) { + ERROR("no basis functions at discrete solution ; doing nothing\n"); + return(0.0); + } int dim = fe_space->getMesh()->getDim(); int deg = grdU.getDegree(); @@ -174,69 +105,74 @@ namespace AMDiS { *q, INIT_GRD_PHI); - relative = relErr; - - maxErr = h1Err2 = h1Norm2 = 0.0; - + double relative = relErr; + double maxErr = 0.0, h1Err2 = 0.0, h1Norm2 = 0.0; + + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1, + Mesh::FILL_COORDS | + Mesh::CALL_LEAF_EL | + Mesh::FILL_DET | + Mesh::FILL_GRD_LAMBDA); + while (elInfo) { + int i, j; + double err, err_2, h1_err_el, norm_el, norm2, det, exact; + const WorldVector *grdu_vec, *grduh_vec; + + grdu_vec = quadFast->getQuadrature()->grdFAtQp(grdErrU, NULL); + det = elInfo->getDet(); + grduh_vec = errUh->getGrdAtQPs(elinfo, NULL, quadFast, NULL); + + int nPoints = quadFast->getNumPoints(); + int dow = Global::getGeo(WORLD); + + for (h1_err_el = i = 0; i < nPoints; i++) { + for (err_2 = j = 0; j < dow; j++) { + err = grdu_vec[i][j] - grduh_vec[i][j]; + err_2 += sqr(err); + } + h1_err_el += quadFast->getWeight(i)*err_2; + } + + exact = det*h1_err_el; + h1Err2 += exact; + maxErr = std::max(maxErr, exact); + + if (writeInLeafData) + elInfo->getElement()->setEstimation(exact, component); + + if (relative) { + for (norm_el = i = 0; i < nPoints; i++) { + for (norm2 = j = 0; j < dow; j++) + norm2 += sqr(grdu_vec[i][j]); + norm_el += quadFast->getWeight(i)*norm2; + } + h1Norm2 += det*norm_el; + } - fe_space->getMesh()->traverse(-1, - Mesh::FILL_COORDS | - Mesh::CALL_LEAF_EL | - Mesh::FILL_DET | - Mesh::FILL_GRD_LAMBDA, - H1ErrFct); + elInfo = stack.traverseNext(elInfo); + } if (relative) { - relNorm2 = h1Norm2 + 1.e-15; - - fe_space->getMesh()->traverse(-1, Mesh::CALL_LEAF_EL, relFct); + double relNorm2 = h1Norm2 + 1.e-15; + + elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + double exact = elInfo->getElement()->getEstimation(component) / relNorm2; + if (writeInLeafData) + elinfo->getElement()->setEstimation(exact, component); + elInfo = stack.traverseNext(elInfo); + } h1Err2 /= relNorm2; maxErr /= relNorm2; } - if (max) *max = maxErr; - - return(sqrt(h1Err2)); - } - - template - int Error::L2ErrFct(ElInfo* el_info) - { - int i; - double err, det, l2_err_el, norm_el, exact; - const double *u_vec, *uh_vec; - elinfo = el_info; - - u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL); - uh_vec = errUh->getVecAtQPs(el_info, - NULL, - quadFast, - NULL); - det = el_info->getDet(); - - int numPoints = quadFast->getNumPoints(); - - for (l2_err_el = i = 0; i < numPoints; i++) { - err = u_vec[i] - uh_vec[i]; - l2_err_el += quadFast->getWeight(i)*sqr(err); - } - - exact = det*l2_err_el; - l2Err2 += exact; - maxErr = max(maxErr, exact); - - if (writeInLeafData) { - el_info->getElement()->setEstimation(exact, component); - } - - if (relative) { - for (norm_el = i = 0; i < numPoints; i++) - norm_el += quadFast->getWeight(i)*sqr(u_vec[i]); - l2Norm2 += det*norm_el; - } + if (max) + *max = maxErr; - return 0; + return sqrt(h1Err2); } template @@ -278,20 +214,58 @@ namespace AMDiS { q = Quadrature::provideQuadrature(dim, degree); quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI); - relative = relErr; - - maxErr = l2Err2 = l2Norm2 = 0.0; + double relative = relErr; + double maxErr = 0.0, l2Err2 = 0.0, l2Norm2 = 0.0; + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1, + Mesh::FILL_COORDS | + Mesh::CALL_LEAF_EL | + Mesh::FILL_DET | + Mesh::FILL_GRD_LAMBDA); + while (elInfo) { + int i; + double err, det, l2_err_el, norm_el, exact; + const double *u_vec, *uh_vec; + + u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL); + uh_vec = errUh->getVecAtQPs(elInfo, NULL, quadFast, NULL); + det = elInfo->getDet(); + + int nPoints = quadFast->getNumPoints(); + + for (l2_err_el = i = 0; i < nPoints; i++) { + err = u_vec[i] - uh_vec[i]; + l2_err_el += quadFast->getWeight(i) * sqr(err); + } + + exact = det * l2_err_el; + l2Err2 += exact; + maxErr = std::max(maxErr, exact); + + if (writeInLeafData) + elInfo->getElement()->setEstimation(exact, component); + + if (relative) { + for (norm_el = i = 0; i < nPoints; i++) + norm_el += quadFast->getWeight(i) * sqr(u_vec[i]); + l2Norm2 += det * norm_el; + } - fe_space->getMesh()->traverse(-1, - Mesh::FILL_COORDS | - Mesh::FILL_DET | - Mesh::FILL_GRD_LAMBDA | - Mesh::CALL_LEAF_EL, - L2ErrFct); + elInfo = stack.traverseNext(elInfo); + } if (relative) { - relNorm2 = l2Norm2 + 1.e-15; - fe_space->getMesh()->traverse(-1, Mesh::CALL_LEAF_EL, relFct); + double relNorm2 = l2Norm2 + 1.e-15; + + elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + double exact = elInfo->getElement()->getEstimation(component) / relNorm2; + if (writeInLeafData) + elInfo->getElement()->setEstimation(exact, component); + elInfo = stack.traverseNext(elInfo); + } + l2Err2 /= relNorm2; } @@ -301,14 +275,4 @@ namespace AMDiS { return (sqrt(l2Err2)); } - template - int Error::relFct(ElInfo* elinfo) - { - double exact = elinfo->getElement()->getEstimation(component); - exact /= relNorm2; - if (writeInLeafData) - elinfo->getElement()->setEstimation(exact, component); - return 0; - } - } diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc index 9e58ec7a..8aedf92e 100644 --- a/AMDiS/src/MacroReader.cc +++ b/AMDiS/src/MacroReader.cc @@ -731,12 +731,12 @@ namespace AMDiS { if (projector && projector->getType() == VOLUME_PROJECTION) { mel[i]->setProjection(0, projector); } else { // boundary projection - for(j = 0; j < mesh->getGeo(NEIGH); j++) { + for (j = 0; j < mesh->getGeo(NEIGH); j++) { projector = Projection::getProjection((*ind)[j]); - if(projector) { + if (projector) { mel[i]->setProjection(j, projector); - if(dim > 2) { - for(k = 0; k < numEdgesAtBoundary; k++) { + if (dim > 2) { + for (k = 0; k < numEdgesAtBoundary; k++) { int edgeNr = Global::getReferenceElement(dim)->getEdgeOfFace(j, k); mel[i]->setProjection(numFaces + edgeNr, projector); } @@ -1091,16 +1091,14 @@ namespace AMDiS { DimVec periodic(dim, DEFAULT_VALUE, false); - if(ed) { + if (ed) { std::list &periodicInfos = dynamic_cast(ed)->getInfoList(); std::list::iterator it; std::list::iterator end = periodicInfos.end(); - for(it = periodicInfos.begin(); it != end; ++it) { - if(it->type != 0) { + for (it = periodicInfos.begin(); it != end; ++it) + if (it->type != 0) periodic[it->elementSide] = true; - } - } } for (k = 0; k < mesh->getGeo(EDGE); k++) { @@ -1155,11 +1153,10 @@ namespace AMDiS { default: ERROR_EXIT("invalid dim\n"); } - if (3 == dim) { + if (3 == dim) mesh->setMaxEdgeNeigh(std::max(8, 2 * max_n_neigh)); - } else { + else mesh->setMaxEdgeNeigh(dim - 1); - } } /* @@ -1222,9 +1219,8 @@ namespace AMDiS { if (test[(*macrolfd)->getIndex()] == 1) { macrolfd++; } else { - for (int i = 0; i < mesh->getNumberOfMacros(); i++) { - zykl[i] = 0; - } + for (int i = 0; i < mesh->getNumberOfMacros(); i++) + zykl[i] = 0; macro = macrolfd; flg = 2; @@ -1254,11 +1250,10 @@ namespace AMDiS { } } while(flg == 2); - if (flg == 1) { + if (flg == 1) macrolfd++; - } else { - macrolfd=mesh->endOfMacroElements(); - } + else + macrolfd=mesh->endOfMacroElements(); } } @@ -1357,7 +1352,7 @@ namespace AMDiS { opp_v = mel->getOppVertex(next_el[edge_no][0]); - if(mel->getBoundary(next_el[edge_no][0])) { + if (mel->getBoundary(next_el[edge_no][0])) { lbound = newBound(mel->getBoundary(next_el[edge_no][0]), lbound); lproject = mel->getProjection(next_el[edge_no][0]); } @@ -1369,7 +1364,7 @@ namespace AMDiS { if (nei->getElement()->getDOF(k) == dof[1]) break; // check for periodic boundary - if(j == 4 || k == 4) { + if (j == 4 || k == 4) { nei = NULL; break; } @@ -1396,7 +1391,7 @@ namespace AMDiS { nei->element->setDOF(node+edge_no,edge_dof); if (next_el[edge_no][0] != opp_v) { - if(nei->getBoundary(next_el[edge_no][0])) { + if (nei->getBoundary(next_el[edge_no][0])) { lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound); Projection *neiProject = nei->getProjection(next_el[edge_no][0]); if (!lproject) { @@ -1433,13 +1428,13 @@ namespace AMDiS { nei = mel->getNeighbour(next_el[edge_no][1]); opp_v = mel->getOppVertex(next_el[edge_no][1]); - if(mel->getBoundary(next_el[edge_no][1])) { + if (mel->getBoundary(next_el[edge_no][1])) { lbound = newBound(mel->getBoundary(next_el[edge_no][1]), lbound); Projection *neiProject = mel->getProjection(next_el[edge_no][1]); - if(!lproject) + if (!lproject) lproject = neiProject; else { - if(neiProject && (lproject->getID() < neiProject->getID())) { + if (neiProject && (lproject->getID() < neiProject->getID())) { lproject = neiProject; } } @@ -1452,7 +1447,7 @@ namespace AMDiS { if (nei->getElement()->getDOF(k) == dof[1]) break; // check for periodic boundary - if(j == 4 || k == 4) { + if (j == 4 || k == 4) { return false; } @@ -1480,7 +1475,7 @@ namespace AMDiS { ("dof %d on element %d is already set, but not on element %d\n", node + edge_no, nei->getIndex(), mel_index); - // if(periodic) { + // if (periodic) { // nei->element->setDOF(node+edge_no, mesh->getDOF(EDGE)); // } else { nei->element->setDOF(node+edge_no,edge_dof); @@ -1489,13 +1484,13 @@ namespace AMDiS { if (next_el[edge_no][0] != opp_v) { - if(nei->getBoundary(next_el[edge_no][0])) { + if (nei->getBoundary(next_el[edge_no][0])) { lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound); Projection *neiProject = nei->getProjection(next_el[edge_no][0]); - if(!lproject) + if (!lproject) lproject = neiProject; else { - if(neiProject &&( lproject->getID() < neiProject->getID())) { + if (neiProject &&( lproject->getID() < neiProject->getID())) { lproject = neiProject; } } @@ -1505,13 +1500,13 @@ namespace AMDiS { nei = nei->getNeighbour(next_el[edge_no][0]); } else { - if(nei->getBoundary(next_el[edge_no][1])) { + if (nei->getBoundary(next_el[edge_no][1])) { lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound); Projection *neiProject = nei->getProjection(next_el[edge_no][1]); - if(!lproject) + if (!lproject) lproject = neiProject; else { - if(neiProject && (lproject->getID() < neiProject->getID())) { + if (neiProject && (lproject->getID() < neiProject->getID())) { lproject = neiProject; } } @@ -1841,13 +1836,17 @@ namespace AMDiS { MSG("checking mesh\n"); - fill_flag = Mesh::CALL_EVERY_EL_INORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND; + fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND; - Mesh::traversePtr = mesh; - error_detected += mesh->traverse(-1, fill_flag|Mesh::FILL_ADD_ALL, basicCheckFct); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, fill_flag | Mesh::FILL_ADD_ALL); + while (elInfo) { + basicCheckFct(elInfo, mesh); + elInfo = stack.traverseNext(elInfo); + } if (mesh->preserveCoarseDOFs) - fill_flag = Mesh::CALL_EVERY_EL_INORDER; + fill_flag = Mesh::CALL_EVERY_EL_PREORDER; else fill_flag = Mesh::CALL_LEAF_EL; @@ -1857,9 +1856,9 @@ namespace AMDiS { localAdmin = mesh->admin[iadmin]; if (localAdmin->getSize() > 0) { - if (static_cast(mesh->dof_used.size()) < localAdmin->getSize()) { + if (static_cast(mesh->dof_used.size()) < localAdmin->getSize()) mesh->dof_used.resize(localAdmin->getSize() + 1000); - } + for (i = 0; i < static_cast(mesh->dof_used.size()); i++) mesh->dof_used[i] = 0; @@ -1884,7 +1883,7 @@ namespace AMDiS { DOFIteratorBase freeIt(localAdmin, FREE_DOFS); for (freeIt.reset(); !freeIt.end(); ++freeIt) { nfree++; - if(mesh->dof_used[freeIt.getDOFIndex()]) { + if (mesh->dof_used[freeIt.getDOFIndex()]) { error_detected++; MSG("dof[%d] used??\n",freeIt.getDOFIndex()); } @@ -1905,30 +1904,31 @@ namespace AMDiS { MSG("checking done; %d error%s detected\n", error_detected, error_detected == 1 ? "" : "s"); - Mesh::traversePtr = mesh; - mesh->traverse(-1, Mesh::CALL_EVERY_EL_INORDER, basicNodeFct); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_INORDER); + while (elInfo) { + basicNodeFct(elInfo, mesh); + elInfo = stack.traverseNext(elInfo); + } WAIT_REALLY; } } - int MacroReader::basicCheckFct(ElInfo* elinfo) + int MacroReader::basicCheckFct(ElInfo* elinfo, Mesh *mesh) { FUNCNAME("MacroReader::basicCheckFct()"); int j, k, opp_v; Element *el = elinfo->getElement(); const Element *neig; - int error_detected=0; - - Mesh* mesh = Mesh::traversePtr; - + int error_detected = 0; int dim = mesh->getDim(); elinfo->testFlag(Mesh::FILL_NEIGH); for (int i = 0; i < mesh->getGeo(NEIGH); i++) { if ((neig = elinfo->getNeighbour(i))) { - if(elinfo->getBoundary(i) > 0) { // < 0 => periodic boundary + if (elinfo->getBoundary(i) > 0) { // < 0 => periodic boundary if (!error_detected) MSG("error detected!!!\n"); error_detected++; @@ -1937,16 +1937,16 @@ namespace AMDiS { } opp_v = elinfo->getOppVertex(i); - if(opp_v < 0 || opp_v >= mesh->getGeo(NEIGH)) { + if (opp_v < 0 || opp_v >= mesh->getGeo(NEIGH)) { if (!error_detected) MSG("error detected!!!\n"); error_detected++; MSG("opp_v = %d\n", opp_v); } - if(elinfo->getBoundary(i) > 0) { // < 0 => periodic boundary - if(dim == 1) { - if(el->getDOF(i) != neig->getDOF(opp_v)) { + if (elinfo->getBoundary(i) > 0) { // < 0 => periodic boundary + if (dim == 1) { + if (el->getDOF(i) != neig->getDOF(opp_v)) { if (!error_detected) MSG("error detected!!!\n"); error_detected++; @@ -2122,18 +2122,14 @@ namespace AMDiS { } } - int MacroReader::basicNodeFct(ElInfo* elinfo) + void MacroReader::basicNodeFct(ElInfo* elinfo, Mesh *mesh) { FUNCNAME("MacroReader::basicNodeFct()"); Element *el = elinfo->getElement(); - Mesh *mesh = Mesh::traversePtr; - MSG("el %4d: ", el->getIndex()); for (int i = 0; i < mesh->getGeo(VERTEX); i++) Msg::print("%4d%s", el->getDOF(i,0), i < mesh->getGeo(VERTEX)-1 ? ", " : "\n"); - - return 0; } } diff --git a/AMDiS/src/MacroReader.h b/AMDiS/src/MacroReader.h index 0b0b430c..413a28fc 100644 --- a/AMDiS/src/MacroReader.h +++ b/AMDiS/src/MacroReader.h @@ -84,11 +84,11 @@ namespace AMDiS { static void checkMesh(Mesh *mesh); - static int basicCheckFct(ElInfo* elInfo); + static int basicCheckFct(ElInfo* elInfo, Mesh *mesh); static void basicDOFCheckFct(ElInfo* elInfo, Mesh *mesh, int iadmin); - static int basicNodeFct(ElInfo* elInfo); + static void basicNodeFct(ElInfo* elInfo, Mesh *mesh); friend class MacroInfo; }; diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index d56f11cc..5760b53e 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -64,7 +64,6 @@ namespace AMDiS { // const Flag Mesh::USE_PARAMETRIC = 0X8000L ; // used in mg methods DOFAdmin* Mesh::compressAdmin = NULL; - Mesh* Mesh::traversePtr = NULL; std::vector Mesh::dof_used; const int Mesh::MAX_DOF = 100; std::map, DegreeOfFreedom*> Mesh::serializedDOFs; @@ -353,34 +352,6 @@ namespace AMDiS { nVertices = nRemainDofs; } - int Mesh::traverse(int level, Flag flag, int (*el_fct)(ElInfo*)) - { - FUNCNAME("Mesh::traverse()"); - - ElInfoStack elInfoStack(this); - ElInfo* elinfo = elInfoStack.getNextElement(); - Traverse tinfo(this, flag, level, el_fct); - - elinfo->setFillFlag(flag); - - if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL) || - flag.isSet(Mesh::CALL_EL_LEVEL) || - flag.isSet(Mesh::CALL_MG_LEVEL)) { - TEST(level >= 0)("invalid level: %d\n", level); - } - - int sum = 0; - for (std::deque::iterator mel = macroElements.begin(); - mel != macroElements.end(); mel++) { - elinfo->fillMacroInfo(*mel); - sum += tinfo.recursive(&elInfoStack); - } - - elInfoStack.getBackElement(); - - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - void Mesh::addDOFAdmin(DOFAdmin *localAdmin) { FUNCNAME("Mesh::addDOFAdmin()"); @@ -464,15 +435,25 @@ namespace AMDiS { compressAdmin->compress(newDOF); - if (preserveCoarseDOFs) { + if (preserveCoarseDOFs) fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING; - } else { + else fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING; - } - traverse(-1, fill_flag, newDOFFct1); - traverse(-1, fill_flag, newDOFFct2); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag); + while (elInfo) { + newDOFFct1(elInfo); + elInfo = stack.traverseNext(elInfo); + } + + elInfo = stack.traverseFirst(this, -1, fill_flag); + while (elInfo) { + newDOFFct2(elInfo); + elInfo = stack.traverseNext(elInfo); + } + newDOF.resize(0); } } diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index bb3fff07..3b7ad11d 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -397,21 +397,6 @@ namespace AMDiS { /// Adds a DOFAdmin to the mesh virtual void addDOFAdmin(DOFAdmin *admin); - /** \brief - * Traverses the mesh. The argument level specifies the element level if - * CALL_EL_LEVEL or CALL_LEAF_EL_LEVEL, or the multigrid level if - * CALL_MG_LEVEL is set. Otherwise this variable is ignored. By the argument - * fillFlag the elements to be traversed and data to be filled into ElInfo is - * selected, using bitwise or of one CALL_... flag and several FILL_... - * flags. The argument elFct is a pointer to a function which is called on - * every element selected by the CALL_... part of fillFlag. - * It is possible to use the recursive mesh traversal recursively, by calling - * traverse() from elFct. - */ - int traverse(int level, - const Flag fillFlag, - int (*elFct)(ElInfo*)); - /// Recalculates the number of leave elements. void updateNumberOfLeaves(); @@ -752,13 +737,6 @@ namespace AMDiS { /// Needed during DOF compression (\ref DOFAdmin::compress). static DOFAdmin *compressAdmin; - /** \brief - * Used for recursive mesh traversal. Static pointer to the mesh - * that should be traversed. This allows access to the mesh even - * from the static traverse routines - */ - static Mesh* traversePtr; - /// Used by check functions static std::vector dof_used; diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index 8846c7cc..9b7abc31 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/AMDiS/src/ProblemVec.cc @@ -25,6 +25,7 @@ #include "TraverseParallel.h" #include "VtkWriter.h" #include "ValueReader.h" +#include "ProblemVecDbg.h" namespace AMDiS { @@ -619,6 +620,8 @@ namespace AMDiS { { FUNCNAME("ProblemVec::buildAfterCoarsen()"); + printOpenmpTraverseInfo(this, true); + // buildAfterCoarsen_sebastianMode(adaptInfo, flag); clock_t first = clock(); diff --git a/AMDiS/src/ProblemVecDbg.cc b/AMDiS/src/ProblemVecDbg.cc index 4702df5c..ec83a4e0 100644 --- a/AMDiS/src/ProblemVecDbg.cc +++ b/AMDiS/src/ProblemVecDbg.cc @@ -401,5 +401,45 @@ namespace AMDiS { elInfo = stack.traverseNext(elInfo); } } - + + void printOpenmpTraverseInfo(ProblemVec *prob, bool wait) + { + int nEl = 0; + std::vector > nElRank(9); + for (int i = 2; i <= 8; i++) { + nElRank[i].resize(i); + for (int j = 0; j < i; j++) + nElRank[i][j] = 0; + } + + Mesh *mesh = prob->getMesh(0); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + nEl++; + int elIndex = elInfo->getElement()->getIndex(); + for (int i = 2; i <= 8; i++) + for (int j = 0; j < i; j++) + if (elIndex % i == j) + nElRank[i][j]++; + + elInfo = stack.traverseNext(elInfo); + } + + std::cout << "==========================\n"; + std::cout << "= OpenMP traverse info =\n"; + std::cout << "==========================\n\n"; + std::cout << "nElements = " << nEl << "\n"; + + for (int i = 2; i <= 8; i++) { + std::cout << i << " threads: "; + for (int j = 0; j < i; j++) { + std::cout << "T" << j << " = " << nElRank[i][j] << " "; + } + std::cout << "\n"; + } + + if (wait) + WAIT_REALLY; + } } diff --git a/AMDiS/src/ProblemVecDbg.h b/AMDiS/src/ProblemVecDbg.h index 665ca261..80758f0b 100644 --- a/AMDiS/src/ProblemVecDbg.h +++ b/AMDiS/src/ProblemVecDbg.h @@ -65,6 +65,7 @@ namespace AMDiS { void createCoordToDofMap(CoordToDof &dofMap); }; + void printOpenmpTraverseInfo(ProblemVec *prob, bool wait = false); } #endif diff --git a/AMDiS/src/RefinementManager.cc b/AMDiS/src/RefinementManager.cc index 02c072cd..acdab746 100644 --- a/AMDiS/src/RefinementManager.cc +++ b/AMDiS/src/RefinementManager.cc @@ -12,19 +12,10 @@ namespace AMDiS { - RefinementManager* RefinementManager::traversePtr = NULL; bool RefinementManager::doMoreRecursiveRefine = false; int RefinementManager::callRefineInterpol = 0; RCNeighbourList* RefinementManager3d::refList = NULL; - int RefinementManager::globalMark = 0; - int RefinementManager::globalRefineFct(ElInfo* elinfo) - { - FUNCNAME("RefineManager::globalRefineFct"); - - elinfo->getElement()->setMark(globalMark); - return 0; - } Flag RefinementManager::globalRefine(Mesh *aMesh, int mark) { @@ -33,15 +24,16 @@ namespace AMDiS { if (mark <= 0) return static_cast(0); - globalMark = mark; - aMesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_BOUND, - globalRefineFct); - return refineMesh(aMesh); + TraverseStack stack; + ElInfo *elInfo = + stack.traverseFirst(aMesh, -1, + Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_BOUND); + while (elInfo) { + elInfo->getElement()->setMark(mark); + elInfo = stack.traverseNext(elInfo); + } } - - - Flag RefinementManager::refineMesh(Mesh *aMesh) { FUNCNAME("RefinementManager::refineMesh()"); @@ -55,8 +47,9 @@ namespace AMDiS { while (doMoreRecursiveRefine) { doMoreRecursiveRefine = false; - elInfo = stack->traverseFirst(mesh, -1, - Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND); + elInfo = + stack->traverseFirst(mesh, -1, + Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND); while (elInfo) { if (elInfo->getElement()->getMark() > 0) { doMoreRecursiveRefine = diff --git a/AMDiS/src/RefinementManager.h b/AMDiS/src/RefinementManager.h index bc88662a..4015eba2 100644 --- a/AMDiS/src/RefinementManager.h +++ b/AMDiS/src/RefinementManager.h @@ -106,11 +106,6 @@ namespace AMDiS { return stack; } - protected: - /// Used by \ref globalRefine - static int globalRefineFct(ElInfo* elinfo); - - protected: /// The Mesh to be refined Mesh *mesh; @@ -126,15 +121,6 @@ namespace AMDiS { /// Used for non recursive traversal TraverseStack* stack; - - /** \brief - * Used while mesh traversal to have a pointer to this RefinementManager - * from a static method - */ - static RefinementManager* traversePtr; - - /// Used by globalRefine to remember the given mark value - static int globalMark; }; } diff --git a/AMDiS/src/RefinementManager1d.cc b/AMDiS/src/RefinementManager1d.cc index f41c1c7b..dacffbb1 100644 --- a/AMDiS/src/RefinementManager1d.cc +++ b/AMDiS/src/RefinementManager1d.cc @@ -13,18 +13,18 @@ namespace AMDiS { - int RefinementManager1d::recursiveRefineFunction(ElInfo* elInfo) + void RefinementManager1d::recursiveRefineFunction(ElInfo* elInfo) { - Line *el = dynamic_cast(const_cast(elInfo->getElement())), *child[2]; + Line *el = + dynamic_cast(const_cast(elInfo->getElement())), *child[2]; Mesh* mesh = elInfo->getMesh(); - if (elInfo->getProjection(0)) { - traversePtr->newCoord(true); - } + if (elInfo->getProjection(0)) + newCoord(true); if (el->getMark() <= 0) - return 0; + return; child[0] = dynamic_cast(mesh->createNewElement(el)); child[1] = dynamic_cast(mesh->createNewElement(el)); @@ -43,7 +43,8 @@ namespace AMDiS { el->setFirstChild(child[0]); el->setSecondChild(child[1]); - if (child[0]->getMark() > 0) traversePtr->doMoreRecursiveRefine = true; + if (child[0]->getMark() > 0) + doMoreRecursiveRefine = true; DegreeOfFreedom *newVertexDOFs = mesh->getDOF(VERTEX); child[0]->setDOF(1, newVertexDOFs); @@ -92,8 +93,6 @@ namespace AMDiS { mesh->freeDOF(const_cast( el->getDOF(mesh->getNode(CENTER))), CENTER); el->setDOF(mesh->getNode(CENTER), NULL); } - - return 0; } Flag RefinementManager1d::refineMesh(Mesh *aMesh) @@ -104,21 +103,24 @@ namespace AMDiS { doMoreRecursiveRefine = true; while (doMoreRecursiveRefine) { doMoreRecursiveRefine = false; - traversePtr = this; - mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS, - recursiveRefineFunction); + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS); + while (elInfo) { + recursiveRefineFunction(elInfo); + elInfo = stack.traverseNext(elInfo); + } } nElements = mesh->getNumberOfLeaves() - nElements; - if (newCoords) { + if (newCoords) setNewCoords(); // call of sub-class method - } return (nElements ? MESH_REFINED : Flag(0)); } - int RefinementManager1d::newCoordsFct(ElInfo *elInfo) + void RefinementManager1d::newCoordsFct(ElInfo *elInfo) { Element *el = elInfo->getElement(); int dow = Global::getGeo(WORLD); @@ -134,15 +136,17 @@ namespace AMDiS { projector->project(*new_coord); el->setNewCoord(new_coord); } - - return 0; } void RefinementManager1d::setNewCoords() { - Flag fillFlag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_BOUND | Mesh::FILL_COORDS; - - mesh->traverse(-1, fillFlag, newCoordsFct); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, + Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_BOUND | Mesh::FILL_COORDS); + while (elInfo) { + newCoordsFct(elInfo); + elInfo = stack.traverseNext(elInfo); + } } } diff --git a/AMDiS/src/RefinementManager1d.h b/AMDiS/src/RefinementManager1d.h index 9a9bae8e..f960b1e0 100644 --- a/AMDiS/src/RefinementManager1d.h +++ b/AMDiS/src/RefinementManager1d.h @@ -47,10 +47,10 @@ namespace AMDiS { protected: /// Used by refineMesh while mesh traversal - static int recursiveRefineFunction(ElInfo* el_info); + void recursiveRefineFunction(ElInfo* el_info); /// Used by \ref setNewCoords - static int newCoordsFct(ElInfo *el_info); + void newCoordsFct(ElInfo *el_info); }; } diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc index 3b9db073..196e4527 100644 --- a/AMDiS/src/RefinementManager2d.cc +++ b/AMDiS/src/RefinementManager2d.cc @@ -141,7 +141,7 @@ namespace AMDiS { } - int RefinementManager2d::newCoordsFct(ElInfo *el_info) + void RefinementManager2d::newCoordsFct(ElInfo *el_info) { Element *el = el_info->getElement(); int dow = Global::getGeo(WORLD); @@ -158,17 +158,19 @@ namespace AMDiS { (*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j]) * 0.5; projector->project(*new_coord); - el->setNewCoord(new_coord); } - - return 0; } void RefinementManager2d::setNewCoords() { - Flag fillFlag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_BOUND | Mesh::FILL_COORDS; - mesh->traverse(-1, fillFlag, newCoordsFct); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, + Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_BOUND | Mesh::FILL_COORDS); + while (elInfo) { + newCoordsFct(elInfo); + elInfo = stack.traverseNext(elInfo); + } } diff --git a/AMDiS/src/RefinementManager2d.h b/AMDiS/src/RefinementManager2d.h index 1b9bc242..3e535eb9 100644 --- a/AMDiS/src/RefinementManager2d.h +++ b/AMDiS/src/RefinementManager2d.h @@ -41,7 +41,7 @@ namespace AMDiS { protected: /// Used by \ref setNewCoords - static int newCoordsFct(ElInfo *el_info); + void newCoordsFct(ElInfo *elInfo); /// Implements RefinementManager::setNewCoords void setNewCoords(); diff --git a/AMDiS/src/RefinementManager3d.cc b/AMDiS/src/RefinementManager3d.cc index 7a0997d7..e2cfa322 100644 --- a/AMDiS/src/RefinementManager3d.cc +++ b/AMDiS/src/RefinementManager3d.cc @@ -48,13 +48,12 @@ namespace AMDiS { child[0]->setDOF(n_vertices-1, dof[0]); child[1]->setDOF(n_vertices-1, dof[0]); - for (i = 0; i < n_vertices-1; i++) - { - child[0]-> - setDOF(i, const_cast( el->getDOF(Tetrahedron::childVertex[el_type][0][i]))); - child[1]-> - setDOF(i, const_cast( el->getDOF(Tetrahedron::childVertex[el_type][1][i]))); - } + for (i = 0; i < n_vertices-1; i++) { + child[0]-> + setDOF(i, const_cast( el->getDOF(Tetrahedron::childVertex[el_type][0][i]))); + child[1]-> + setDOF(i, const_cast( el->getDOF(Tetrahedron::childVertex[el_type][1][i]))); + } /****************************************************************************/ /* there is one more leaf element and two more hierachical elements */ /****************************************************************************/ @@ -67,193 +66,169 @@ namespace AMDiS { /* information */ /****************************************************************************/ - if (mesh->getNumberOfDOFs(EDGE)) - { - node = mesh->getNode(EDGE); + if (mesh->getNumberOfDOFs(EDGE)) { + node = mesh->getNode(EDGE); - /****************************************************************************/ - /* set pointers to those dof's that are handed on from the parant */ - /****************************************************************************/ + /****************************************************************************/ + /* set pointers to those dof's that are handed on from the parant */ + /****************************************************************************/ + + child[0]-> + setDOF(node, + const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][0][0]))); + child[1]-> + setDOF(node, + const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][1][0]))); + child[0]-> + setDOF(node+1, + const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][0][1]))); + child[1]-> + setDOF(node+1, + const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][1][1]))); + child[0]-> + setDOF(node+3, + const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][0][3]))); + child[1]-> + setDOF(node+3, + const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][1][3]))); + + /****************************************************************************/ + /* adjust pointers to the dof's in the refinement edge */ + /****************************************************************************/ + + if (el->getDOF(0) == edge[0]) { + child[0]->setDOF(node+2, dof[1]); + child[1]->setDOF(node+2, dof[2]); + } else { + child[0]->setDOF(node+2, dof[2]); + child[1]->setDOF(node+2, dof[1]); + } + } + + if (mesh->getNumberOfDOFs(FACE)) { + node = mesh->getNode(FACE); + + /****************************************************************************/ + /* set pointers to those dof's that are handed on from the parant */ + /****************************************************************************/ + + child[0]->setDOF(node+3, const_cast( el->getDOF(node+1))); + child[1]->setDOF(node+3, const_cast( el->getDOF(node+0))); + + /****************************************************************************/ + /* get new dof for the common face of child0 and child1 */ + /****************************************************************************/ + + DegreeOfFreedom *newDOF = mesh->getDOF(FACE); + child[0]->setDOF(node, static_cast( newDOF)); + child[1]->setDOF(node, static_cast( newDOF)); + } + + if (mesh->getNumberOfDOFs(CENTER)) { + node = mesh->getNode(CENTER); + child[0]->setDOF(node, const_cast( mesh->getDOF(CENTER))); + child[1]->setDOF(node, const_cast( mesh->getDOF(CENTER))); + } - child[0]-> - setDOF(node, - const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][0][0]))); - child[1]-> - setDOF(node, - const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][1][0]))); - child[0]-> - setDOF(node+1, - const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][0][1]))); - child[1]-> - setDOF(node+1, - const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][1][1]))); - child[0]-> - setDOF(node+3, - const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][0][3]))); - child[1]-> - setDOF(node+3, - const_cast( el->getDOF(node+Tetrahedron::childEdge[el_type][1][3]))); + if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(FACE)) + fillPatchConnectivity(ref_list, index); + } - /****************************************************************************/ - /* adjust pointers to the dof's in the refinement edge */ - /****************************************************************************/ + void RefinementManager3d::fillPatchConnectivity(RCNeighbourList* ref_list, + int index) + { + FUNCNAME("RefinementManager3d::fillPatchConnectivity"); + Element *el = ref_list->getElement(index), *neigh; + int dir, el_type = ref_list->getType(index), n_type = 0; + int adjc, i, j, i_neigh, j_neigh; + int node0, node1, opp_v = 0; - if (el->getDOF(0) == edge[0]) - { - child[0]->setDOF(node+2, dof[1]); - child[1]->setDOF(node+2, dof[2]); - } - else - { - child[0]->setDOF(node+2, dof[2]); - child[1]->setDOF(node+2, dof[1]); - } + for (dir = 0; dir < 2; dir++) { + neigh = ref_list->getNeighbourElement(index, dir); + if (neigh) { + n_type = ref_list->getType( ref_list->getNeighbourNr(index, dir) ); + opp_v = ref_list->getOppVertex(index, dir); } - if (mesh->getNumberOfDOFs(FACE)) - { - node = mesh->getNode(FACE); + if (!neigh || neigh->isLeaf()) { /****************************************************************************/ - /* set pointers to those dof's that are handed on from the parant */ + /* get new dof's in the midedge of the face of el and for the two midpoints*/ + /* of the sub-faces. If face is an interior face those pointers have to be */ + /* adjusted by the neighbour element also (see below) */ /****************************************************************************/ - child[0]->setDOF(node+3, const_cast( el->getDOF(node+1))); - child[1]->setDOF(node+3, const_cast( el->getDOF(node+0))); - + if (mesh->getNumberOfDOFs(EDGE)) { + node0 = node1 = mesh->getNode(EDGE); + node0 += Tetrahedron::nChildEdge[el_type][0][dir]; + node1 += Tetrahedron::nChildEdge[el_type][1][dir]; + DegreeOfFreedom *newDOF = mesh->getDOF(EDGE); + (const_cast(el->getFirstChild()))->setDOF(node0, newDOF); + (const_cast(el->getSecondChild()))->setDOF(node1, newDOF); + } + if (mesh->getNumberOfDOFs(FACE)) { + node0 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir]; + (const_cast(el->getFirstChild()))->setDOF(node0, mesh->getDOF(FACE)); + node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir]; + (const_cast(el->getSecondChild()))->setDOF(node1, mesh->getDOF(FACE)); + } + } else { /* if (!neigh || !neigh->child[0]) */ /****************************************************************************/ - /* get new dof for the common face of child0 and child1 */ + /* interior face and neighbour has been refined, look for position at the */ + /* refinement edge */ /****************************************************************************/ + + if (el->getDOF(0) == neigh->getDOF(0)) { + /****************************************************************************/ + /* same position at refinement edge */ + /****************************************************************************/ + adjc = 0; + } else { + /****************************************************************************/ + /* different position at refinement edge */ + /****************************************************************************/ + adjc = 1; + } - DegreeOfFreedom *newDOF = mesh->getDOF(FACE); - child[0]->setDOF(node, static_cast( newDOF)); - child[1]->setDOF(node, static_cast( newDOF)); - } - - if (mesh->getNumberOfDOFs(CENTER)) - { - node = mesh->getNode(CENTER); - child[0]->setDOF(node, const_cast( mesh->getDOF(CENTER))); - child[1]->setDOF(node, const_cast( mesh->getDOF(CENTER))); - } + for (i = 0; i < 2; i++) { + j = Tetrahedron::adjacentChild[adjc][i]; - if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(FACE)) - fillPatchConnectivity(ref_list, index); + i_neigh = Tetrahedron::nChildFace[el_type][i][dir]; + j_neigh = Tetrahedron::nChildFace[n_type][j][opp_v-2]; - // MSG("%d -> %d %d\n", - // el->getIndex(), - // child[0]->getIndex(), - // child[1]->getIndex()); - } + /****************************************************************************/ + /* adjust dof pointer in the edge in the common face of el and neigh and */ + /* the dof pointer in the sub-face child_i-child_j (allocated by neigh!) */ + /****************************************************************************/ - void RefinementManager3d::fillPatchConnectivity(RCNeighbourList* ref_list, - int index) - { - FUNCNAME("RefinementManager3d::fillPatchConnectivity"); - Element *el = ref_list->getElement(index), *neigh; - int dir, el_type = ref_list->getType(index), n_type = 0; - int adjc, i, j, i_neigh, j_neigh; - int node0, node1, opp_v = 0; + if (mesh->getNumberOfDOFs(EDGE)) { + node0 = + mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir]; + node1 = + mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][opp_v-2]; - for (dir = 0; dir < 2; dir++) - { - neigh = ref_list->getNeighbourElement(index, dir); - if (neigh) - { - n_type = ref_list->getType( ref_list->getNeighbourNr(index, dir) ); - opp_v = ref_list->getOppVertex(index, dir); - } + TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1)) + ("no dof on neighbour %d at node %d\n", + neigh->getChild(j)->getIndex(), node1); - if (!neigh || neigh->isLeaf()) - { + (const_cast( el->getChild(i)))-> + setDOF(node0, const_cast( neigh->getChild(j)->getDOF(node1))); + } + if (mesh->getNumberOfDOFs(FACE)) { + node0 = mesh->getNode(FACE) + i_neigh; + node1 = mesh->getNode(FACE) + j_neigh; - /****************************************************************************/ - /* get new dof's in the midedge of the face of el and for the two midpoints*/ - /* of the sub-faces. If face is an interior face those pointers have to be */ - /* adjusted by the neighbour element also (see below) */ - /****************************************************************************/ + TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1)) + ("no dof on neighbour %d at node %d\n", + neigh->getChild(j)->getIndex(), node1); - if (mesh->getNumberOfDOFs(EDGE)) - { - node0 = node1 = mesh->getNode(EDGE); - node0 += Tetrahedron::nChildEdge[el_type][0][dir]; - node1 += Tetrahedron::nChildEdge[el_type][1][dir]; - DegreeOfFreedom *newDOF = mesh->getDOF(EDGE); - (const_cast( el->getFirstChild()))->setDOF(node0, newDOF); - (const_cast( el->getSecondChild()))->setDOF(node1, newDOF); - } - if (mesh->getNumberOfDOFs(FACE)) - { - node0 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir]; - (const_cast( el->getFirstChild()))->setDOF(node0, mesh->getDOF(FACE)); - node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir]; - (const_cast( el->getSecondChild()))->setDOF(node1, mesh->getDOF(FACE)); - } + (const_cast( el->getChild(i)))-> + setDOF(node0, const_cast( neigh->getChild(j)->getDOF(node1))); } - else /* if (!neigh || !neigh->child[0]) */ - { - /****************************************************************************/ - /* interior face and neighbour has been refined, look for position at the */ - /* refinement edge */ - /****************************************************************************/ - - if (el->getDOF(0) == neigh->getDOF(0)) - { - /****************************************************************************/ - /* same position at refinement edge */ - /****************************************************************************/ - adjc = 0; - } - else - { - /****************************************************************************/ - /* different position at refinement edge */ - /****************************************************************************/ - adjc = 1; - } - for (i = 0; i < 2; i++) - { - j = Tetrahedron::adjacentChild[adjc][i]; - - i_neigh = Tetrahedron::nChildFace[el_type][i][dir]; - j_neigh = Tetrahedron::nChildFace[n_type][j][opp_v-2]; - - /****************************************************************************/ - /* adjust dof pointer in the edge in the common face of el and neigh and */ - /* the dof pointer in the sub-face child_i-child_j (allocated by neigh!) */ - /****************************************************************************/ - - if (mesh->getNumberOfDOFs(EDGE)) - { - node0 = - mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir]; - node1 = - mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][opp_v-2]; - - TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1)) - ("no dof on neighbour %d at node %d\n", - neigh->getChild(j)->getIndex(), node1); - - (const_cast( el->getChild(i)))-> - setDOF(node0, const_cast( neigh->getChild(j)->getDOF(node1))); - } - if (mesh->getNumberOfDOFs(FACE)) - { - node0 = mesh->getNode(FACE) + i_neigh; - node1 = mesh->getNode(FACE) + j_neigh; - - TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1)) - ("no dof on neighbour %d at node %d\n", - neigh->getChild(j)->getIndex(), node1); - - (const_cast( el->getChild(i)))-> - setDOF(node0, const_cast( neigh->getChild(j)->getDOF(node1))); - } - - } /* for (i = 0; i < 2; i++) */ - } /* else of if (!neigh || !neigh->child[0]) */ - } /* for (dir = 0; dir < 2; dir++) */ + } /* for (i = 0; i < 2; i++) */ + } /* else of if (!neigh || !neigh->child[0]) */ + } /* for (dir = 0; dir < 2; dir++) */ } @@ -268,7 +243,7 @@ namespace AMDiS { Projection *projector = el_info->getProjection(0); - if(!projector || projector->getType() != VOLUME_PROJECTION) { + if (!projector || projector->getType() != VOLUME_PROJECTION) { projector = el_info->getProjection(4); } @@ -288,23 +263,24 @@ namespace AMDiS { refList->setElement(0, el); refList->setElType(0, dynamic_cast(el_info)->getType()); n_neigh = 1; + for (i = 0; i < 2; i++) edge[i] = const_cast( el_info->getElement()->getDOF(i)); - if (getRefinePatch(&elinfo, edge, 0, refList, &n_neigh)) - { - /****************************************************************************/ - /* domain's boundary was reached while looping around the refinement edge */ - /****************************************************************************/ - getRefinePatch(&elinfo, edge, 1, refList, &n_neigh); - } - for (i = 1; i < n_neigh; i++) /* start with 1, as list[0]=el */ - { - TEST(!refList->getElement(i)->isNewCoordSet()) - ("non-nil new_coord in el %d ref_list[%d] el %d (n_neigh=%d)\n", - el->getIndex(), i, refList->getElement(i)->getIndex(), n_neigh); - refList->getElement(i)->setNewCoord(el->getNewCoord()); - } + if (getRefinePatch(&elinfo, edge, 0, refList, &n_neigh)) { + /****************************************************************************/ + /* domain's boundary was reached while looping around the refinement edge */ + /****************************************************************************/ + getRefinePatch(&elinfo, edge, 1, refList, &n_neigh); + } + + for (i = 1; i < n_neigh; i++) { /* start with 1, as list[0]=el */ + TEST(!refList->getElement(i)->isNewCoordSet()) + ("non-nil new_coord in el %d ref_list[%d] el %d (n_neigh=%d)\n", + el->getIndex(), i, refList->getElement(i)->getIndex(), n_neigh); + + refList->getElement(i)->setNewCoord(el->getNewCoord()); + } } return 0; @@ -348,53 +324,48 @@ namespace AMDiS { dof[0] = mesh->getDOF(VERTEX); mesh->incrementNumberOfVertices(1); - if (mesh->getNumberOfDOFs(EDGE)) - { - dof[1] = mesh->getDOF(EDGE); - dof[2] = mesh->getDOF(EDGE); - } + if (mesh->getNumberOfDOFs(EDGE)) { + dof[1] = mesh->getDOF(EDGE); + dof[2] = mesh->getDOF(EDGE); + } - for (i = 0; i < n_neigh; i++) { + for (i = 0; i < n_neigh; i++) bisectTetrahedron(refineList, i, dof, edge); - } /****************************************************************************/ /* if there are functions to interpolate data to the finer grid, do so */ /****************************************************************************/ int iadmin; int nrAdmin = mesh->getNumberOfDOFAdmin(); - for(iadmin = 0; iadmin < nrAdmin; iadmin++) { + for (iadmin = 0; iadmin < nrAdmin; iadmin++) { std::list::iterator it; DOFAdmin* admin = const_cast(&mesh->getDOFAdmin(iadmin)); std::list::iterator end = admin->endDOFIndexed(); - for(it = admin->beginDOFIndexed(); it != end; it++) + for (it = admin->beginDOFIndexed(); it != end; it++) (*it)->refineInterpol(*refineList, n_neigh); } - if (!mesh->queryCoarseDOFs()) - { + if (!mesh->queryCoarseDOFs()) { + /****************************************************************************/ + /* if there should be no dof information on interior leaf elements remove */ + /* dofs from edges, faces and the centers of parents */ + /****************************************************************************/ + if (mesh->getNumberOfDOFs(EDGE)) { /****************************************************************************/ - /* if there should be no dof information on interior leaf elements remove */ - /* dofs from edges, faces and the centers of parents */ + /* remove dof of the midpoint of the common refinement edge */ /****************************************************************************/ - if (mesh->getNumberOfDOFs(EDGE)) - { - /****************************************************************************/ - /* remove dof of the midpoint of the common refinement edge */ - /****************************************************************************/ - - el = dynamic_cast(const_cast( refineList->getElement(0))); - mesh->freeDOF(const_cast( el->getDOF(mesh->getNode(EDGE))), EDGE); - } - if (mesh->getNumberOfDOFs(EDGE) || - mesh->getNumberOfDOFs(FACE) || - mesh->getNumberOfDOFs(CENTER)) - { - for (i = 0; i < n_neigh; i++) - refineList->removeDOFParent(i); - } + el = dynamic_cast(const_cast( refineList->getElement(0))); + mesh->freeDOF(const_cast( el->getDOF(mesh->getNode(EDGE))), EDGE); } + + if (mesh->getNumberOfDOFs(EDGE) || + mesh->getNumberOfDOFs(FACE) || + mesh->getNumberOfDOFs(CENTER)) { + for (i = 0; i < n_neigh; i++) + refineList->removeDOFParent(i); + } + } /****************************************************************************/ @@ -402,18 +373,15 @@ namespace AMDiS { /* is a boundary edge or not */ /****************************************************************************/ - if (bound) - { - mesh->incrementNumberOfEdges(n_neigh + 2); - mesh->incrementNumberOfFaces(2*n_neigh + 1); - newCoords = true; - } - else - { - mesh->incrementNumberOfEdges(n_neigh + 1); - mesh->incrementNumberOfFaces(2*n_neigh); - } - + if (bound) { + mesh->incrementNumberOfEdges(n_neigh + 2); + mesh->incrementNumberOfFaces(2*n_neigh + 1); + newCoords = true; + } else { + mesh->incrementNumberOfEdges(n_neigh + 1); + mesh->incrementNumberOfFaces(2*n_neigh); + } + return dof[0][0]; } @@ -450,13 +418,13 @@ namespace AMDiS { if (neigh->getDOF(k) == edge[1]) break; - if(j > 3 || k > 3) { + if (j > 3 || k > 3) { for (j = 0; j < vertices; j++) if (mesh->associated(neigh->getDOF(j, 0), edge[0][0])) break; for (k = 0; k < vertices; k++) if (mesh->associated(neigh->getDOF(k, 0), edge[1][0])) break; - if(j > 3 || k > 3) { + if (j > 3 || k > 3) { for (j = 0; j < vertices; j++) if (mesh->indirectlyAssociated(neigh->getDOF(j, 0), edge[0][0])) break; for (k = 0; k < vertices; k++) @@ -473,15 +441,15 @@ namespace AMDiS { // LeafDataPeriodic *pd = // dynamic_cast(el->getElementData(PERIODIC)); - // if(pd) { + // if (pd) { // std::list::iterator it; // std::list::iterator end = // pd->getInfoList().end(); - // for(it = pd->getInfoList().begin(); it != end; ++it) { - // if(it->periodicMode != 0) { + // for (it = pd->getInfoList().begin(); it != end; ++it) { + // if (it->periodicMode != 0) { // PeriodicBC *periodicCondition = mesh->getPeriodicBCMap()[it->type]; - // if(periodicCondition) { + // if (periodicCondition) { // VertexVector *associated = mesh->getPeriodicAssociations()[it->type]; // for (j = 0; j < vertices; j++) // if (neigh->getDOF(j, 0) == (*associated)[edge[0][0]]) break; @@ -506,117 +474,108 @@ namespace AMDiS { edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0), neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0)); - if ((edge_no = Tetrahedron::edgeOfDOFs[j][k])) - { - /****************************************************************************/ - /* neigh has not a compatible refinement edge */ - /****************************************************************************/ - neigh->setMark(max(neigh->getMark(), 1)); - neigh_info = refineFunction(neigh_info); - - /****************************************************************************/ - /* now, go to a child at the edge and determine the opposite vertex for */ - /* this child; continue the looping around the edge with this element */ - /****************************************************************************/ - + if ((edge_no = Tetrahedron::edgeOfDOFs[j][k])) { + /****************************************************************************/ + /* neigh has not a compatible refinement edge */ + /****************************************************************************/ + neigh->setMark(max(neigh->getMark(), 1)); + neigh_info = refineFunction(neigh_info); + + /****************************************************************************/ + /* now, go to a child at the edge and determine the opposite vertex for */ + /* this child; continue the looping around the edge with this element */ + /****************************************************************************/ + + neigh_info = stack->traverseNext(neigh_info); + neigh_el_type = dynamic_cast(neigh_info)->getType(); + + switch (edge_no) { + case 1: + opp_v = opp_v == 1 ? 3 : 2; + break; + case 2: + opp_v = opp_v == 2 ? 1 : 3; + break; + case 3: neigh_info = stack->traverseNext(neigh_info); neigh_el_type = dynamic_cast(neigh_info)->getType(); - - switch (edge_no) - { - case 1: - opp_v = opp_v == 1 ? 3 : 2; - break; - case 2: - opp_v = opp_v == 2 ? 1 : 3; - break; - case 3: - neigh_info = stack->traverseNext(neigh_info); - neigh_el_type = dynamic_cast(neigh_info)->getType(); - if (neigh_el_type != 1) - opp_v = opp_v == 0 ? 3 : 2; - else - opp_v = opp_v == 0 ? 3 : 1; - break; - case 4: - neigh_info = stack->traverseNext(neigh_info); - neigh_el_type = dynamic_cast(neigh_info)->getType(); - if (neigh_el_type != 1) - opp_v = opp_v == 0 ? 3 : 1; - else - opp_v = opp_v == 0 ? 3 : 2; - break; - case 5: - if (neigh_el_type != 1) - { - neigh_info = stack->traverseNext(neigh_info); - neigh_el_type = (dynamic_cast(neigh_info))->getType(); - } - opp_v = 3; - break; - } - neigh = dynamic_cast(const_cast( neigh_info->getElement())); - } - else - { - /****************************************************************************/ - /* neigh is compatible devisible; put neigh to the list of patch elements; */ - /* go to next neighbour */ - /****************************************************************************/ - TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh()) - ("too many neighbours %d in refine patch\n", *n_neigh); - - - refineList->setElement(*n_neigh, neigh); - refineList->setElType(*n_neigh, neigh_el_type); - - /****************************************************************************/ - /* we have to go back to the starting element via opp_v values */ - /* correct information is produce by get_neigh_on_patch() */ - /****************************************************************************/ - refineList->setOppVertex(*n_neigh, 0, opp_v); - - ++*n_neigh; - - if (opp_v != 3) - i = 3; + if (neigh_el_type != 1) + opp_v = opp_v == 0 ? 3 : 2; else - i = 2; - - if (neigh_info->getNeighbour(i)) - { - opp_v = neigh_info->getOppVertex(i); - neigh_info = stack->traverseNeighbour3d(neigh_info, i); - neigh_el_type = (dynamic_cast(neigh_info))->getType(); - neigh = dynamic_cast(const_cast( neigh_info->getElement())); - } + opp_v = opp_v == 0 ? 3 : 1; + break; + case 4: + neigh_info = stack->traverseNext(neigh_info); + neigh_el_type = dynamic_cast(neigh_info)->getType(); + if (neigh_el_type != 1) + opp_v = opp_v == 0 ? 3 : 1; else - break; + opp_v = opp_v == 0 ? 3 : 2; + break; + case 5: + if (neigh_el_type != 1) { + neigh_info = stack->traverseNext(neigh_info); + neigh_el_type = (dynamic_cast(neigh_info))->getType(); + } + opp_v = 3; + break; } - } + neigh = dynamic_cast(const_cast( neigh_info->getElement())); + } else { + /****************************************************************************/ + /* neigh is compatible devisible; put neigh to the list of patch elements; */ + /* go to next neighbour */ + /****************************************************************************/ + TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh()) + ("too many neighbours %d in refine patch\n", *n_neigh); + + + refineList->setElement(*n_neigh, neigh); + refineList->setElType(*n_neigh, neigh_el_type); + + /****************************************************************************/ + /* we have to go back to the starting element via opp_v values */ + /* correct information is produce by get_neigh_on_patch() */ + /****************************************************************************/ + refineList->setOppVertex(*n_neigh, 0, opp_v); + + ++*n_neigh; - if (neigh == el) - { - (*el_info) = neigh_info; - return(0); + if (opp_v != 3) + i = 3; + else + i = 2; + + if (neigh_info->getNeighbour(i)) { + opp_v = neigh_info->getOppVertex(i); + neigh_info = stack->traverseNeighbour3d(neigh_info, i); + neigh_el_type = (dynamic_cast(neigh_info))->getType(); + neigh = dynamic_cast(const_cast( neigh_info->getElement())); + } else + break; } - + } + + if (neigh == el) { + (*el_info) = neigh_info; + return(0); + } + /****************************************************************************/ /* the domain's boundary is reached; loop back to the starting el */ /****************************************************************************/ - + i = *n_neigh-1; opp_v = refineList->getOppVertex(i, 0); - do - { - TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v) && i > 0) - ("while looping back domains boundary was reached or i == 0\n"); - opp_v = refineList->getOppVertex(i--, 0); - neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v); - } while(neigh_info->getElement() != el); + do { + TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v) && i > 0) + ("while looping back domains boundary was reached or i == 0\n"); + opp_v = refineList->getOppVertex(i--, 0); + neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v); + } while (neigh_info->getElement() != el); (*el_info) = neigh_info; - - return(1); + + return 1; } @@ -642,29 +601,24 @@ namespace AMDiS { /* give the refinement edge the right orientation */ /****************************************************************************/ - if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1,0)) - { - edge[0] = const_cast( el_info->getElement()->getDOF(0)); - edge[1] = const_cast( el_info->getElement()->getDOF(1)); - } - else - { - edge[1] = const_cast( el_info->getElement()->getDOF(0)); - edge[0] = const_cast( el_info->getElement()->getDOF(1)); - } + if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1,0)) { + edge[0] = const_cast( el_info->getElement()->getDOF(0)); + edge[1] = const_cast( el_info->getElement()->getDOF(1)); + } else { + edge[1] = const_cast( el_info->getElement()->getDOF(0)); + edge[0] = const_cast( el_info->getElement()->getDOF(1)); + } /****************************************************************************/ /* get the refinement patch */ /****************************************************************************/ - // MSG("index %d\n", el_info->getElement()->getIndex()); - if (getRefinePatch(&el_info, edge, 0, ref_list, &n_neigh)) - { - /****************************************************************************/ - /* domain's boundary was reached while looping around the refinement edge */ - /****************************************************************************/ - getRefinePatch(&el_info, edge, 1, ref_list, &n_neigh); - bound = true; - } + if (getRefinePatch(&el_info, edge, 0, ref_list, &n_neigh)) { + /****************************************************************************/ + /* domain's boundary was reached while looping around the refinement edge */ + /****************************************************************************/ + getRefinePatch(&el_info, edge, 1, ref_list, &n_neigh); + bound = true; + } /****************************************************************************/ /* fill neighbour information inside the patch in the refinement list */ @@ -688,11 +642,6 @@ namespace AMDiS { std::map::iterator it; std::map::iterator end = mesh->getPeriodicAssociations().end(); - // for(int i = 0; i < n_neigh; i++) { - // MSG("%d ", ref_list->getElement(i)->getIndex()); - // } - // MSG("\n"); - while(edge[0] != NULL) { periodicList = ref_list->periodicSplit(edge, next_edge, @@ -701,17 +650,15 @@ namespace AMDiS { TEST_EXIT_DBG(periodicList)("periodicList = NULL\n"); - newDOF = - refinePatch(edge, periodicList, n_neigh_periodic, bound); + newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound); - if(firstNewDOF == -1) { + if (firstNewDOF == -1) firstNewDOF = newDOF; - } - if(lastNewDOF != -1) { - for(it = mesh->getPeriodicAssociations().begin(); it != end; ++it) { - if(it->second) { - if(((*(it->second))[edge[0][0]] == last_edge[0][0] && + if (lastNewDOF != -1) { + for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) { + if (it->second) { + if (((*(it->second))[edge[0][0]] == last_edge[0][0] && (*(it->second))[edge[1][0]] == last_edge[1][0]) || ((*(it->second))[edge[0][0]] == last_edge[1][0] && (*(it->second))[edge[1][0]] == last_edge[0][0])) @@ -731,10 +678,10 @@ namespace AMDiS { edge[1] = next_edge[1]; } - if(lastNewDOF != firstNewDOF) { - for(it = mesh->getPeriodicAssociations().begin(); it != end; ++it) { - if(it->second) { - if(((*(it->second))[first_edge[0][0]] == last_edge[0][0] && + if (lastNewDOF != firstNewDOF) { + for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) { + if (it->second) { + if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] && (*(it->second))[first_edge[1][0]] == last_edge[1][0]) || ((*(it->second))[first_edge[0][0]] == last_edge[1][0] && (*(it->second))[first_edge[1][0]] == last_edge[0][0])) diff --git a/AMDiS/src/TecPlotWriter.hh b/AMDiS/src/TecPlotWriter.hh index f507af49..bf833935 100644 --- a/AMDiS/src/TecPlotWriter.hh +++ b/AMDiS/src/TecPlotWriter.hh @@ -28,13 +28,13 @@ namespace AMDiS { int dim = elinfo->getMesh()->getDim(); // int dow = Global::getGeo(WORLD); - for(int i=0; i < dim+1; i++) { + for (int i=0; i < dim+1; i++) { typename DOFCoords::iterator coords = find(dofCoords[dof[i][n0]].begin(), dofCoords[dof[i][n0]].end(), elinfo->getCoord(i)); TEST_EXIT(coords != dofCoords[dof[i][n0]].end())("coords not found"); (*outFile) << coords->vertex_index << " "; - if(dim == 1 && i == 1) { + if (dim == 1 && i == 1) { (*outFile) << coords->vertex_index; } } @@ -49,6 +49,8 @@ namespace AMDiS { const char* plotTitle, bool additional) { + FUNCNAME("TecPlotWriter::writeValues()"); + values = val; TEST_EXIT(values)("no values\n"); @@ -57,7 +59,7 @@ namespace AMDiS { TEST_EXIT(values->getFESpace())("no fe-space\n"); TEST_EXIT(values->getFESpace()->getMesh())("no mesh\n"); - if(!additional) + if (!additional) outFile = new std::ofstream(filename); else outFile = new std::ofstream(filename, std::ios::app); @@ -77,42 +79,44 @@ namespace AMDiS { writeCoords = !additional; - if(dim == 1 && dow == 1) { // use ordered plot - if(!additional) { + if (dim == 1 && dow == 1) { // use ordered plot + if (!additional) (*outFile) << "TITLE = \"" << plotTitle << "\"" << std::endl; - } // write file header - if(!additional) { + if (!additional) writeVarName(dim); - } + (*outFile) << "ZONE T=\"" << values->getName() << "\", I=" << mesh->getNumberOfVertices(); (*outFile) << ", F=POINT"; - if(additional) { + if (additional) (*outFile) << ", D=(1)"; - } + (*outFile) << std::endl << std::endl; // write data - mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS, writeValuesFct); + ERROR_EXIT("Reimplement mesh traverse for TecPlot writer!\n"); + // mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS, writeValuesFct); } else { // dow !=1 => use FE-plot - for(int i=0; i < 150; i++) + for (int i=0; i < 150; i++) (*outFile) << " "; (*outFile) << std::endl; // write data nv = 1; - mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS, writeValuesFct); + ERROR_EXIT("Reimplement mesh traverse for TecPlot writer!\n"); + // mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS, writeValuesFct); (*outFile) << std::endl; - mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS, writeIndicesFct); + ERROR_EXIT("Reimplement mesh traverse for TecPlot writer!\n"); + // mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS, writeIndicesFct); // write file header outFile->seekp(0); - if(!additional) { + if (!additional) { (*outFile) << "TITLE = \"" << plotTitle << "\"" << std::endl; writeVarName(dim); } @@ -121,7 +125,7 @@ namespace AMDiS { (*outFile) << ", E=" << mesh->getNumberOfLeaves(); (*outFile) << ", F=FEPOINT"; (*outFile) << ((dim==3) ? ", ET=TETRAHEDRON" : ", ET=TRIANGLE"); - if(additional) { + if (additional) { (*outFile) << ", D=("; (*outFile) << ((dim==3) ? "1,2,3" : "1,2"); (*outFile) << ",FECONNECT)"; diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc index b62d032f..d349485d 100644 --- a/AMDiS/src/Traverse.cc +++ b/AMDiS/src/Traverse.cc @@ -78,214 +78,12 @@ namespace AMDiS { } if (elinfo) { - if (parametric) { + if (parametric) elinfo = parametric->addParametricInfo(elinfo); - } elinfo->fillDetGrdLambda(); } - return(elinfo); - } - - - - /****************************************************************************/ - /* common (static) traversal routines for 2d and 3d */ - /* to be #included in traverse_r.c */ - /****************************************************************************/ - - /* traverse hierarchical mesh, call el_fct() for each/some element */ - - /****************************************************************************/ - /* recursive_traverse: */ - /* ------------------- */ - /* recursively traverse the mesh hierarchy tree */ - /* call the routine traverse_info->el_fct(el_info) with */ - /* - all tree leaves */ - /* - all tree leaves at a special level */ - /* - all tree elements in pre-/in-/post-order */ - /* depending on the traverse_info->level variable */ - /****************************************************************************/ - - int Traverse::recursive(ElInfoStack *elInfoStack) - { - FUNCNAME("Traverse::recursive()"); - - ElInfo *elinfo = elInfoStack->getCurrentElement(); - - Element *el = elinfo->getElement(); - int mg_level, sum = 0; - Parametric *parametric = mesh->getParametric(); - - if (flag.isSet(Mesh::CALL_LEAF_EL)) { - if (el->getFirstChild()) { - ElInfo* elinfo_new = elInfoStack->getNextElement(); - elinfo_new->fillElInfo(0, elinfo); - sum += recursive(elInfoStack); - elinfo_new->fillElInfo(1, elinfo); - sum += recursive(elInfoStack); - elInfoStack->getBackElement(); - } else { - if (el_fct != NULL) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - } - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - - if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) { - if (el->getFirstChild()) { - if ((elinfo->getLevel() >= level)) { - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - ElInfo* elinfo_new = elInfoStack->getNextElement(); - elinfo_new->fillElInfo(0, elinfo); - sum += recursive(elInfoStack); - elinfo->fillElInfo(1, elinfo); - sum += recursive(elInfoStack); - elInfoStack->getBackElement(); - } else { - if ((elinfo->getLevel() == level) && (el_fct!=NULL)) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - } - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - - if (flag.isSet(Mesh::CALL_EL_LEVEL)) { - if (elinfo->getLevel() == level) { - if (NULL != el_fct) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - } else { - if (elinfo->getLevel() > level){ - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } else if (el->getFirstChild()) { - ElInfo* elinfo_new = elInfoStack->getNextElement(); - elinfo_new->fillElInfo(0, elinfo); - sum += recursive(elInfoStack); - elinfo_new->fillElInfo(1, elinfo); - sum += recursive(elInfoStack); - elInfoStack->getBackElement(); - } - } - - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - - if (flag.isSet(Mesh::CALL_MG_LEVEL)) { - - mg_level = (elinfo->getLevel() + mesh->getDim() - 1) / mesh->getDim(); - - if (mg_level > level) { - return 0; - } - - if (!(el->getFirstChild())) { - if (el_fct != NULL) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - - if ((mg_level == level) && ((elinfo->getLevel() % mesh->getDim()) == 0)) { - if (el_fct != NULL) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - - ElInfo* elinfo_new = elInfoStack->getNextElement(); - - elinfo_new->fillElInfo(0, elinfo); - sum += recursive(elInfoStack); - - elinfo_new->fillElInfo(1, elinfo); - sum += recursive(elInfoStack); - - elInfoStack->getBackElement(); - - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; - } - - if (flag.isSet(Mesh::CALL_EVERY_EL_PREORDER)) { - if (el_fct != NULL) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - } - - if (el->getFirstChild()) { - ElInfo* elinfo_new = elInfoStack->getNextElement(); - - elinfo_new->fillElInfo(0, elinfo); - sum += recursive(elInfoStack); - - if (flag.isSet(Mesh::CALL_EVERY_EL_INORDER)) - if (el_fct!=NULL) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - elinfo_new->fillElInfo(1, elinfo); - sum += recursive(elInfoStack); - - elInfoStack->getBackElement(); - } else { - if (flag.isSet(Mesh::CALL_EVERY_EL_INORDER)) - if (el_fct != NULL) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - } - - if (flag.isSet(Mesh::CALL_EVERY_EL_POSTORDER)) - if (el_fct != NULL) { - if (parametric) - elinfo = parametric->addParametricInfo(elinfo); - elinfo->fillDetGrdLambda(); - sum += el_fct(elinfo); - if (parametric) - elinfo = parametric->removeParametricInfo(elinfo); - } - - return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0; + return elinfo; } void TraverseStack::enlargeTraverseStack() diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h index 2ecb504b..b4ded5bd 100644 --- a/AMDiS/src/Traverse.h +++ b/AMDiS/src/Traverse.h @@ -264,9 +264,6 @@ namespace AMDiS { TEST_EXIT(m)("No traverse without mesh!\n"); } - /// Performs the recursive traversal - int recursive(ElInfoStack*); - private: Mesh *mesh; diff --git a/AMDiS/src/TraverseParallel.cc b/AMDiS/src/TraverseParallel.cc index e3a198b2..4eedae94 100644 --- a/AMDiS/src/TraverseParallel.cc +++ b/AMDiS/src/TraverseParallel.cc @@ -6,7 +6,9 @@ namespace AMDiS { - TraverseParallelStack::TraverseParallelStack(int n) + TraverseParallelStack::TraverseParallelStack(int n, int mode) + : nThreads(0), + parallelMode(mode) { if (n == 0) nThreads = omp_get_overall_max_threads(); @@ -14,13 +16,20 @@ namespace AMDiS { nThreads = n; stacks.resize(nThreads); - int i = 0; - for (std::vector::iterator it = stacks.begin(); - it != stacks.end(); ++it) { - (*it) = new TraverseStack(); - (*it)->setMyThreadId(i++); - (*it)->setMaxThreads(nThreads); + if (parallelMode == 0) { + int i = 0; + for (std::vector::iterator it = stacks.begin(); + it != stacks.end(); ++it) { + + (*it) = new TraverseStack(); + (*it)->setMyThreadId(i++); + (*it)->setMaxThreads(nThreads); + } + } else { + for (std::vector::iterator it = stacks.begin(); + it != stacks.end(); ++it) + (*it) = new TraverseStack(); } } diff --git a/AMDiS/src/TraverseParallel.h b/AMDiS/src/TraverseParallel.h index 67077386..97844789 100644 --- a/AMDiS/src/TraverseParallel.h +++ b/AMDiS/src/TraverseParallel.h @@ -36,7 +36,7 @@ namespace AMDiS { class TraverseParallelStack { public: - TraverseParallelStack(int nThreads = 0); + TraverseParallelStack(int nThreads = 0, int mode = 0); ~TraverseParallelStack(); @@ -44,13 +44,19 @@ namespace AMDiS { inline ElInfo* traverseNext(ElInfo* elInfoOld) { - return stacks[omp_get_thread_num()]->traverseNext(elInfoOld); + if (parallelMode == 0) { + return stacks[omp_get_thread_num()]->traverseNext(elInfoOld); + } else { + + } } private: /// Number of threads using the stack in parallel. int nThreads; + int parallelMode; + std::vector stacks; }; -- GitLab