diff --git a/AMDiS/src/AdaptInfo.h b/AMDiS/src/AdaptInfo.h index 66bc9339a9476f3cb49f355940694f8af1078b5e..eb095322a6f8b07ca610f75b554f81092986faea 100644 --- a/AMDiS/src/AdaptInfo.h +++ b/AMDiS/src/AdaptInfo.h @@ -389,24 +389,12 @@ namespace AMDiS { /// Sets \ref est_sum. inline void setEstSum(double e, int index) { -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS - double send_est_sum = e; - double est_sum = 0.0; - MPI::COMM_WORLD.Allreduce(&send_est_sum, &est_sum, 1, MPI_DOUBLE, MPI_SUM); - e = est_sum; -#endif scalContents[index]->est_sum = e; } /// Sets \ref est_max. inline void setEstMax(double e, int index) { -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS - double send_est_max = e; - double est_max = 0.0; - MPI::COMM_WORLD.Allreduce(&send_est_max, &est_max, 1, MPI_DOUBLE, MPI_MAX); - e = est_max; -#endif scalContents[index]->est_max = e; } diff --git a/AMDiS/src/BoundaryObject.cc b/AMDiS/src/BoundaryObject.cc index 7890608b3ba2d048817fb7c4ca16baa1e2b313f1..7bc1615ad349c78d842fb9c96b527e25dcdddfa1 100644 --- a/AMDiS/src/BoundaryObject.cc +++ b/AMDiS/src/BoundaryObject.cc @@ -128,6 +128,19 @@ namespace AMDiS { } + bool BoundaryObject::operator<(const BoundaryObject& other) const + { + if (elIndex == other.elIndex) { + if (subObj == other.subObj) + return ithObj < other.ithObj; + + return subObj < other.subObj; + } + + return (elIndex < other.elIndex); + } + + bool AtomicBoundary::operator==(const AtomicBoundary& other) const { return (rankObj == other.rankObj && diff --git a/AMDiS/src/BoundaryObject.h b/AMDiS/src/BoundaryObject.h index d34cda375b7f12a5db7eeda672963d5b61914865..ec273251d1dfdf47685900332d343a87fbe237b4 100644 --- a/AMDiS/src/BoundaryObject.h +++ b/AMDiS/src/BoundaryObject.h @@ -55,6 +55,8 @@ namespace AMDiS { bool operator!=(const BoundaryObject& other) const; + bool operator<(const BoundaryObject& other) const; + /// The macro element to which the boundary element corresponds to. Element* el; @@ -64,11 +66,9 @@ namespace AMDiS { /// Element type index, only used in 3d. int elType; - /** \brief - * Defines the geometrical object at the boundary. It must be "a part" of the - * macro element \ref el, i.e., either 1 (a vertex), 2 (an edge) or 3 - * (a face). - */ + /// Defines the geometrical object at the boundary. It must be "a part" of + /// the macro element \ref el, i.e., either 1 (a vertex), 2 (an edge) or 3 + /// (a face). GeoIndex subObj; /** \brief diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h index d9967656ae67292047fcc0f0daa0c1e510e86ad5..f09ebc85942dcdbc7531a9ba4b7731e365643a0b 100644 --- a/AMDiS/src/ElInfo.h +++ b/AMDiS/src/ElInfo.h @@ -472,41 +472,29 @@ namespace AMDiS { /// \ref element is child of element parent Element *parent; - /** \brief - * \ref element is an element of the binary tree located at MacroElement - * macroElement - */ + /// \ref element is an element of the binary tree located at MacroElement + /// macroElement MacroElement *macroElement; - /** \brief - * Indicates wich elements will be called and wich information should be - * present while mesh traversal. - */ + /// Indicates wich elements will be called and wich information should be + /// present while mesh traversal. Flag fillFlag; - /** \brief - * Level of the element. The level is zero for macro elements and the level - * of the children is (level of the parent + 1). level_ is filled always by - * the traversal routines. - */ + /// Level of the element. The level is zero for macro elements and the level + /// of the children is (level of the parent + 1). level_ is filled always by + /// the traversal routines. unsigned char level; - /** \brief - * Elements type index. This is used only for 3d, where type can be either 0, 1 or - * 2. In all other cases type is not used and the variable is set to 0. - * In 3d, it is filled automatically by the traversal routines. - */ + /// Elements type index. This is used only for 3d, where type can be either 0, 1 or + /// 2. In all other cases type is not used and the variable is set to 0. + /// In 3d, it is filled automatically by the traversal routines. int elType; - /** \brief - * This ElInfo is the iChild-th child of the parent element. - */ + /// 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. - */ + /// \ref coord[i] is a WorldVector<double> storing the world coordinates of the + /// i-th vertex of element \ref element. FixVec<WorldVector<double>, VERTEX> coord; /** \brief @@ -518,34 +506,24 @@ namespace AMDiS { */ FixVec<BoundaryType, BOUNDARY> boundary; - /** \brief - * Vector storing pointers to projections for each face, edge, vertex. - */ + /// Vector storing pointers to projections for each face, edge, vertex. FixVec<Projection*, PROJECTION> projection; - /** \brief - * oppCoord[i] coordinates of the i-th neighbour vertex opposite the - * common edge/face. - */ + /// oppCoord[i] coordinates of the i-th neighbour vertex opposite the + /// common edge/face. FixVec<WorldVector<double>, NEIGH> oppCoord; - /** \brief - * neighbour[i] pointer to the element at the edge/face with local index i. - * It is a pointer to NULL for boundary edges/faces. - */ + /// neighbour[i] pointer to the element at the edge/face with local index i. + /// It is a pointer to NULL for boundary edges/faces. FixVec<Element*, NEIGH> neighbour; - /** \brief - * neighbourCoord[i][j] are the coordinate of the j-th vertex of the i-th - * neighbour element with the common edge/face. - */ + /// neighbourCoord[i][j] are the coordinate of the j-th vertex of the i-th + /// neighbour element with the common edge/face. FixVec<FixVec<WorldVector<double>, VERTEX>, NEIGH> neighbourCoord; - /** \brief - * oppVertex[i] is undefined if neighbour[i] is a pointer to NULL. - * Otherwise it is the local index of the neighbour's vertex opposite the - * common edge/face. - */ + /// oppVertex[i] is undefined if neighbour[i] is a pointer to NULL. + /// Otherwise it is the local index of the neighbour's vertex opposite the + /// common edge/face. FixVec<int, NEIGH> oppVertex; /// Elements determinant. @@ -569,18 +547,14 @@ namespace AMDiS { static std::vector<std::map<std::pair<int, unsigned long>, mtl::dense2D<double> > > subElemGradMatrices; - /** \brief - * child_vertex[el_type][child][i] = father's local vertex index of new - * vertex i. 4 stands for the newly generated vertex . - */ + /// child_vertex[el_type][child][i] = father's local vertex index of new + /// vertex i. 4 stands for the newly generated vertex . static const int childVertex[3][2][4]; - /** \brief - * child_edge[el_type][child][i] = father's local edge index of new edge i. - * new edge 2 is half of old edge 0, new edges 4,5 are really new edges, and - * value is different: child_edge[][][4,5] = index of same edge in other - * child. - */ + /// child_edge[el_type][child][i] = father's local edge index of new edge i. + /// new edge 2 is half of old edge 0, new edges 4,5 are really new edges, and + /// value is different: child_edge[][][4,5] = index of same edge in other + /// child. static const int childEdge[3][2][6]; friend class ElInfo1d; diff --git a/AMDiS/src/Marker.cc b/AMDiS/src/Marker.cc index d00b8de0c96b7f5edb5ece9ce256f50635d78a16..dae5790294b142fe00424106e74bbdc4ac734190 100644 --- a/AMDiS/src/Marker.cc +++ b/AMDiS/src/Marker.cc @@ -140,7 +140,7 @@ namespace AMDiS { markRLimit = ESThetaP * epsP / nLeaves; markCLimit = ESThetaCP * epsP / nLeaves; - + INFO(info, 2)("start mark_limits: %.3le %.3le nt = %d\n", markRLimit, markCLimit, nLeaves); } diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc index 6ae98684fc7982c45c8aa93d717aa3b3d88a35ee..a52db0219ef0ee70b0966108a4ddb57e53fae45d 100644 --- a/AMDiS/src/ProblemStat.cc +++ b/AMDiS/src/ProblemStat.cc @@ -1726,11 +1726,12 @@ namespace AMDiS { map<int, double> vec; TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(this->getMesh(comp), -1, - Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); + ElInfo *elInfo = + stack.traverseFirst(this->getMesh(comp), -1, Mesh::CALL_LEAF_EL); while (elInfo) { - vec[elInfo->getElement()->getIndex()] = elInfo->getElement()->getEstimation(comp); + vec[elInfo->getElement()->getIndex()] = + elInfo->getElement()->getEstimation(comp); elInfo = stack.traverseNext(elInfo); } diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h index e2a254ca39143a4fd5539ec5f1f8ffb39b218feb..8424a213ff6cb357e82837a01d7ce1bef6293f2a 100644 --- a/AMDiS/src/Quadrature.h +++ b/AMDiS/src/Quadrature.h @@ -157,19 +157,15 @@ namespace AMDiS { - /** \brief - * Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th - * quadrature point - */ + /// Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th + /// quadrature point inline double getLambda(int a, int b) const { return (lambda ? (*lambda)[a][b] : 0.0); } - /** \brief - * Returns \ref lambda[a] which is a DimVec<double> containing the - * coordiantes of the a-th quadrature point - */ + /// Returns \ref lambda[a] which is a DimVec<double> containing the + /// coordiantes of the a-th quadrature point inline const DimVec<double>& getLambda(int a) const { return (*lambda)[a]; @@ -206,10 +202,8 @@ namespace AMDiS { double *w; protected: - /** \brief - * Initialisation of all static Quadrature objects which will be returned - * by \ref provideQuadrature() - */ + /// Initialisation of all static Quadrature objects which will be returned + /// by \ref provideQuadrature() static void initStaticQuadratures(); /** \name static quadratures, used weights, and barycentric coords @@ -288,21 +282,15 @@ namespace AMDiS { - /** \brief - * Pre-compute the values of all basis functions at all quadrature nodes; - */ + /// Pre-compute the values of all basis functions at all quadrature nodes; const Flag INIT_PHI=1; - /** \brief - * Pre-compute the gradients (with respect to the barycentric coordinates) of - * all basis functions at all quadrature nodes - */ + /// Pre-compute the gradients (with respect to the barycentric coordinates) of + /// all basis functions at all quadrature nodes const Flag INIT_GRD_PHI=2; - /** \brief - * pre-compute all 2nd derivatives (with respect to the barycentric - * coordinates) of all basis functions at all quadrature nodes; - */ + /// pre-compute all 2nd derivatives (with respect to the barycentric + /// coordinates) of all basis functions at all quadrature nodes; const Flag INIT_D2_PHI=4; @@ -325,10 +313,8 @@ namespace AMDiS { class FastQuadrature { protected: - /** \brief - * Constructs a FastQuadrature for the given Quadrature, BasisFunction, and - * flag. - */ + /// Constructs a FastQuadrature for the given Quadrature, BasisFunction, and + /// flag. FastQuadrature(BasisFunction* basFcts, Quadrature* quad, Flag flag) : init_flag(flag), phi(0, 0), @@ -481,10 +467,8 @@ namespace AMDiS { } protected: - /** \brief - * Specifies which information should be pre-computed. Can be \ref INIT_PHI, - * \ref INIT_GRD_PHI, or \ref INIT_D2_PHI - */ + /// Specifies which information should be pre-computed. Can be \ref INIT_PHI, + /// \ref INIT_GRD_PHI, or \ref INIT_D2_PHI Flag init_flag; /** \brief @@ -513,16 +497,12 @@ namespace AMDiS { */ MatrixOfFixVecs<DimMat<double> > *D2Phi; - /** \brief - * List of all used FastQuadratures - */ + /// List of all used FastQuadratures static list<FastQuadrature*> fastQuadList; - /** \brief - * Maximal number of quadrature points for all yet initialised FastQuadrature - * objects. This value may change after a new initialisation of a - * FastQuadrature - */ + /// Maximal number of quadrature points for all yet initialised FastQuadrature + /// objects. This value may change after a new initialisation of a + /// FastQuadrature static int max_points; /// This FastQuadrature stores values for Quadrature quadrature diff --git a/AMDiS/src/est/ResidualEstimator.cc b/AMDiS/src/est/ResidualEstimator.cc index c96da9e627047320ed0835653bb8ced8de5eccfc..56deb13a0568174301ec84fac6f9d5fbbcef6326 100644 --- a/AMDiS/src/est/ResidualEstimator.cc +++ b/AMDiS/src/est/ResidualEstimator.cc @@ -173,7 +173,6 @@ namespace AMDiS { InteriorBoundary &intBoundary = MeshDistributor::globalMeshDistributor->getIntBoundary(); - RankToBoundMap& otherBound = intBoundary.getOther(); ElInfo *elInfo = mesh->createNewElInfo(); elInfo->setFillFlag(Mesh::FILL_COORDS); @@ -181,56 +180,66 @@ namespace AMDiS { StdMpi<vector<double> > stdMpiDet(MeshDistributor::globalMeshDistributor->getMpiComm()); StdMpi<vector<vector<WorldVector<double> > > > stdMpiGrdUh(MeshDistributor::globalMeshDistributor->getMpiComm()); - for (RankToBoundMap::iterator it = otherBound.begin(); - it != otherBound.end(); ++it) { + RankToBoundMap allBounds = intBoundary.getOther(); + allBounds.insert(intBoundary.getOwn().begin(), intBoundary.getOwn().end()); + + for (RankToBoundMap::iterator it = allBounds.begin(); + it != allBounds.end(); ++it) { + + vector<BoundaryObject> subBound; + for (unsigned int i = 0; i < it->second.size(); i++) { BoundaryObject &bObj = it->second[i].rankObj; if (bObj.subObj == VERTEX) continue; - vector<BoundaryObject> subBound; bObj.el->getSubBoundary(bObj, subBound); + } - WorldVector<int> faceIndices; - - for (unsigned int j = 0; j < subBound.size(); j++) { - Element *el = subBound[j].el; - int oppV = subBound[j].ithObj; - - elInfo->setElement(el); - el->sortFaceIndices(oppV, faceIndices); - - elInfo->getCoord(oppV) = coords[el->getDof(oppV, 0)]; - for (int k = 0; k < dim; k++) - elInfo->getCoord(faceIndices[k]) = coords[el->getDof(k, 0)]; - - double detNeigh = abs(elInfo->calcGrdLambda(*lambdaNeigh)); - stdMpiDet.getSendData(it->first).push_back(detNeigh); - - for (int system = 0; system < nSystems; system++) { - if (matrix[system] == NULL || secondOrderTerms[system] == false) - continue; - - uh[system]->getLocalVector(el, uhNeigh[system]); - - for (int iq = 0; iq < nPointsSurface; iq++) { + if (subBound.size() == 0) + continue; - (*lambda)[oppV] = 0.0; - for (int k = 0; k < dim; k++) - (*lambda)[faceIndices[k]] = surfaceQuad->getLambda(iq, k); + WorldVector<int> faceIndices; - basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]); - } + for (unsigned int i = 0; i < subBound.size(); i++) { + Element *el = subBound[i].el; + int oppV = subBound[i].ithObj; + + elInfo->setElement(el); + el->sortFaceIndices(oppV, faceIndices); + for (int k = 0; k <= dim; k++) + elInfo->getCoord(k) = coords[el->getDof(k, 0)]; + + double detNeigh = abs(elInfo->calcGrdLambda(*lambdaNeigh)); + stdMpiDet.getSendData(it->first).push_back(detNeigh); - stdMpiGrdUh.getSendData(it->first).push_back(grdUhNeigh); + for (int system = 0; system < nSystems; system++) { + if (matrix[system] == NULL || secondOrderTerms[system] == false) + continue; + + uh[system]->getLocalVector(el, uhNeigh[system]); + + for (int iq = 0; iq < nPointsSurface; iq++) { + + (*lambda)[oppV] = 0.0; + for (int k = 0; k < dim; k++) + (*lambda)[faceIndices[k]] = surfaceQuad->getLambda(iq, k); + + basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]); } + + stdMpiGrdUh.getSendData(it->first).push_back(grdUhNeigh); } } + + stdMpiDet.recv(it->first); + stdMpiGrdUh.recv(it->first); } stdMpiDet.updateSendDataSize(); stdMpiGrdUh.updateSendDataSize(); +#if 0 RankToBoundMap& ownBound = intBoundary.getOwn(); for (RankToBoundMap::iterator it = ownBound.begin(); @@ -245,13 +254,14 @@ namespace AMDiS { } } } +#endif stdMpiDet.startCommunication(); stdMpiGrdUh.startCommunication(); - for (RankToBoundMap::iterator it = ownBound.begin(); - it != ownBound.end(); ++it) { + for (RankToBoundMap::iterator it = allBounds.begin(); + it != allBounds.end(); ++it) { vector<BoundaryObject> subBound; for (unsigned int i = 0; i < it->second.size(); i++) { @@ -270,6 +280,11 @@ namespace AMDiS { TEST_EXIT_DBG(subBound.size() == stdMpiGrdUh.getRecvData(it->first).size()) ("Should not happen!\n"); + + for (unsigned int i = 0; i < subBound.size(); i++) { + elBoundDet[subBound[i]] = stdMpiDet.getRecvData(it->first)[i]; + elBoundGrdUhNeigh[subBound[i]] = stdMpiGrdUh.getRecvData(it->first)[i]; + } } MSG("INIT PARALLEL FINISCHED!\n"); @@ -482,20 +497,41 @@ namespace AMDiS { const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda(); double det = elInfo->getDet(); double h2 = h2_from_det(det, dim); + BoundaryObject blub; for (int face = 0; face < nNeighbours; face++) { Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face)); + bool parallelMode = false; + +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS + if (neigh == NULL) { + BoundaryObject testObj(el, elInfo->getType(), EDGE, face); + + if (elBoundDet.count(testObj)) { + parallelMode = true; + blub = testObj; + } else { + continue; + } + } + + if (!parallelMode && !(neigh && neigh->getMark())) + continue; +#else if (!(neigh && neigh->getMark())) continue; +#endif - int oppV = elInfo->getOppVertex(face); - + int oppV = elInfo->getOppVertex(face); el->sortFaceIndices(face, faceIndEl); - neigh->sortFaceIndices(oppV, faceIndNeigh); - neighInfo->setElement(const_cast<Element*>(neigh)); - neighInfo->setFillFlag(Mesh::FILL_COORDS); - neighInfo->getCoord(oppV) = elInfo->getOppCoord(face); + + if (!parallelMode) { + neigh->sortFaceIndices(oppV, faceIndNeigh); + neighInfo->setElement(const_cast<Element*>(neigh)); + neighInfo->setFillFlag(Mesh::FILL_COORDS); + neighInfo->getCoord(oppV) = elInfo->getOppCoord(face); + } // periodic leaf data ? ElementData *ldp = el->getElementData()->getElementData(PERIODIC); @@ -526,7 +562,7 @@ namespace AMDiS { } } // if (ldp) - if (!periodicCoords) { + if (!periodicCoords && !parallelMode) { for (int i = 0; i < dim; i++) { int i1 = faceIndEl[i]; int i2 = faceIndNeigh[i]; @@ -534,11 +570,19 @@ namespace AMDiS { } } + double detNeigh = 0.0; Parametric *parametric = mesh->getParametric(); - if (parametric) - neighInfo = parametric->addParametricInfo(neighInfo); + if (!parallelMode) { + if (parametric) + neighInfo = parametric->addParametricInfo(neighInfo); - double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh)); + detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh)); + } else { + TEST_EXIT_DBG(!parametric)("No yet implemented!\n"); + + detNeigh = elBoundDet[blub]; + } + for (int iq = 0; iq < nPointsSurface; iq++) jump[iq].set(0.0); @@ -548,7 +592,8 @@ namespace AMDiS { continue; uh[system]->getLocalVector(el, uhEl[system]); - uh[system]->getLocalVector(neigh, uhNeigh[system]); + if (!parallelMode) + uh[system]->getLocalVector(neigh, uhNeigh[system]); for (int iq = 0; iq < nPointsSurface; iq++) { (*lambda)[face] = 0.0; @@ -556,12 +601,18 @@ namespace AMDiS { (*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i); basFcts[system]->evalGrdUh(*lambda, grdLambda, uhEl[system], grdUhEl[iq]); + + if (!parallelMode) { + (*lambda)[oppV] = 0.0; + for (int i = 0; i < dim; i++) + (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i); + + basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]); + } else { + grdUhNeigh[iq] = elBoundGrdUhNeigh[blub][iq]; + } - (*lambda)[oppV] = 0.0; - for (int i = 0; i < dim; i++) - (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i); - - basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]); + grdUhEl[iq] -= grdUhNeigh[iq]; } // for iq @@ -599,11 +650,13 @@ namespace AMDiS { val *= C1 * h2_from_det(d, dim) * d; else val *= C1 * d; - - if (parametric) - neighInfo = parametric->removeParametricInfo(neighInfo); - - neigh->setEstimation(neigh->getEstimation(row) + val, row); + + if (!parallelMode) { + if (parametric) + neighInfo = parametric->removeParametricInfo(neighInfo); + neigh->setEstimation(neigh->getEstimation(row) + val, row); + } + result += val; } // for face diff --git a/AMDiS/src/est/ResidualEstimator.h b/AMDiS/src/est/ResidualEstimator.h index ec46b71a4b8e7e9ed9aef4a569b87ba95d70bf9e..bc7846ec06f45b0ffb7612efee8067f835a8b52b 100644 --- a/AMDiS/src/est/ResidualEstimator.h +++ b/AMDiS/src/est/ResidualEstimator.h @@ -188,6 +188,12 @@ namespace AMDiS { /// are used to ommit computations of the jump residual that is defined /// only on second order terms. std::vector<bool> secondOrderTerms; + +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS + map<BoundaryObject, double> elBoundDet; + + map<BoundaryObject, std::vector<WorldVector<double> > > elBoundGrdUhNeigh; +#endif }; }