#include #include #include #include "time.h" #include "AdaptStationary.h" #include "AdaptInstationary.h" #include "FiniteElemSpace.h" #include "ElementData.h" #include "ElementDofIterator.h" #include "MacroElement.h" #include "MacroReader.h" #include "MacroInfo.h" #include "Mesh.h" #include "Traverse.h" #include "Parameters.h" #include "FixVec.h" #include "DOFVector.h" #include "CoarseningManager.h" #include "DOFIterator.h" #include "VertexVector.h" #include "MacroWriter.h" #include "PeriodicMap.h" #include "Projection.h" #include "ElInfoStack.h" #include "Serializer.h" #include "Lagrange.h" namespace AMDiS { #define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC) //************************************************************************** // flags, which information should be present in the elInfo structure //************************************************************************** const Flag Mesh::FILL_NOTHING = 0X00L; const Flag Mesh::FILL_COORDS = 0X01L; const Flag Mesh::FILL_BOUND = 0X02L; const Flag Mesh::FILL_NEIGH = 0X04L; const Flag Mesh::FILL_OPP_COORDS = 0X08L; const Flag Mesh::FILL_ORIENTATION= 0X10L; const Flag Mesh::FILL_DET = 0X20L; const Flag Mesh::FILL_GRD_LAMBDA = 0X40L; const Flag Mesh::FILL_ADD_ALL = 0X80L; const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L); const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L); const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L); //************************************************************************** // flags for Mesh traversal //************************************************************************** const Flag Mesh::CALL_EVERY_EL_PREORDER = 0X0100L; const Flag Mesh::CALL_EVERY_EL_INORDER = 0X0200L; const Flag Mesh::CALL_EVERY_EL_POSTORDER = 0X0400L; const Flag Mesh::CALL_LEAF_EL = 0X0800L; const Flag Mesh::CALL_LEAF_EL_LEVEL = 0X1000L; const Flag Mesh::CALL_EL_LEVEL = 0X2000L; const Flag Mesh::CALL_MG_LEVEL = 0X4000L; const Flag Mesh::CALL_REVERSE_MODE = 0X8000L; std::vector Mesh::dof_used; const int Mesh::MAX_DOF = 100; std::map, DegreeOfFreedom*> Mesh::serializedDOFs; struct delmem { DegreeOfFreedom* ptr; int len; }; Mesh::Mesh(std::string aName, int dimension) : name(aName), dim(dimension), nVertices(0), nEdges(0), nLeaves(0), nElements(0), parametric(NULL), preserveCoarseDOFs(false), nDOFEl(0), nDOF(dimension, DEFAULT_VALUE, 0), nNodeEl(0), node(dimension, DEFAULT_VALUE, 0), elementPrototype(NULL), elementDataPrototype(NULL), elementIndex(-1), initialized(false), macroFileInfo(NULL), changeIndex(0), final_lambda(dimension, DEFAULT_VALUE, 0.0) { FUNCNAME("Mesh::Mesh()"); // set default element prototype switch(dim) { case 1: elementPrototype = new Line(this); break; case 2: elementPrototype = new Triangle(this); break; case 3: elementPrototype = new Tetrahedron(this); break; default: ERROR_EXIT("invalid dimension\n"); } elementPrototype->setIndex(-1); elementIndex = 0; } Mesh::~Mesh() { deleteMeshStructure(); if (macroFileInfo != NULL) delete macroFileInfo; if (elementPrototype) delete elementPrototype; if (elementDataPrototype) delete elementDataPrototype; for (int i = 0; i < static_cast(admin.size()); i++) delete admin[i]; } Mesh& Mesh::operator=(const Mesh& m) { FUNCNAME("Mesh::operator=()"); if (this == &m) return *this; TEST_EXIT(dim == m.dim)("operator= works only on meshes with equal dim!\n"); name = m.name; nVertices = m.nVertices; nEdges = m.nEdges; nLeaves = m.nLeaves; nElements = m.nElements; nFaces = m.nFaces; maxEdgeNeigh = m.maxEdgeNeigh; diam = m.diam; parametric = NULL; preserveCoarseDOFs = m.preserveCoarseDOFs; nDOFEl = m.nDOFEl; nDOF = m.nDOF; nNodeEl = m.nNodeEl; node = m.node; newDOF = m.newDOF; elementIndex = m.elementIndex; initialized = m.initialized; final_lambda = m.final_lambda; /* ====================== Create new DOFAdmins ================== */ admin.resize(m.admin.size()); for (int i = 0; i < static_cast(admin.size()); i++) { admin[i] = new DOFAdmin(this); *(admin[i]) = *(m.admin[i]); admin[i]->setMesh(this); } /* ====================== Copy macro elements =================== */ // mapIndex[i] is the index of the MacroElement element in the vector // macroElements, for which holds: element->getIndex() = i std::map mapIndex; // We use this map for coping the DOFs of the Elements within the // MacroElements objects. Mesh::serializedDOFs.clear(); int insertCounter = 0; macroElements.clear(); // Go through all MacroElements of mesh m, and create for every a new // MacroElement in this mesh. for (std::deque::const_iterator it = m.macroElements.begin(); it != m.macroElements.end(); ++it, insertCounter++) { // Create new MacroElement. MacroElement *el = new MacroElement(dim); // Use copy operator to copy all the data to the new MacroElement. *el = **it; // Make a copy of the Element data, together with all DOFs el->setElement((*it)->getElement()->cloneWithDOFs()); // Insert the new MacroElement in the vector of all MacroElements. macroElements.push_back(el); // Update the index map. mapIndex.insert(std::pair(el->getIndex(), insertCounter)); } // Now we have to go through all the new MacroElements, and update the neighbour // connections. insertCounter = 0; for (std::deque::const_iterator it = m.macroElements.begin(); it != m.macroElements.end(); ++it, insertCounter++) { // Go through all neighbours. for (int i = 0; i < dim; i++) { // 1. Get index of the old MacroElement for its i-th neighbour. // 2. Because the index in the new MacroElement is the same, search // for the vector index the corresponding element is stored in. // 3. Get this element from macroElements, and set it as the i-th // neighbour for the current element. macroElements[insertCounter]-> setNeighbour(i, macroElements[mapIndex[(*it)->getNeighbour(i)->getIndex()]]); } } // Cleanup Mesh::serializedDOFs.clear(); /* ================== Things will be done when required ============ */ TEST_EXIT(elementDataPrototype == NULL)("TODO\n"); TEST_EXIT(m.parametric == NULL)("TODO\n"); TEST_EXIT(periodicAssociations.size() == 0)("TODO\n"); return *this; } void Mesh::updateNumberOfLeaves() { nLeaves = 0; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL); while (elInfo) { nLeaves++; elInfo = stack.traverseNext(elInfo); } } void Mesh::addMacroElement(MacroElement* me) { macroElements.push_back(me); me->setIndex(macroElements.size()); } void Mesh::removeMacroElements(std::set& macros, const FiniteElemSpace *feSpace) { FUNCNAME("Mesh::removeMacroElement()"); typedef std::map > DofElMap; typedef std::map DofPosMap; TEST_EXIT(admin.size() == 1)("Not yet implemented for multiple admins!\n"); TEST_EXIT(admin[0])("There is something wrong!\n"); ElementDofIterator elDofIter(feSpace); // Map that stores for each dof pointer (which may have a list of dofs) // all macro element indices that own this dof. DofElMap dofsOwner; DofPosMap dofsPosIndex; // === Determine to all dofs the macro elements where the dof is part of. === for (std::deque::iterator macroIt = macroElements.begin(); macroIt != macroElements.end(); ++macroIt) { elDofIter.reset((*macroIt)->getElement()); do { dofsOwner[elDofIter.getDofPtr()].insert(*macroIt); dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex(); } while (elDofIter.nextStrict()); } // === Remove macro elements from mesh macro element list. === // Removing arbitrary elements from an std::deque is very slow. Therefore, we // create a new deque with all macro elements that should not be deleted. The // macro element deque is than replaced by the new created one. std::deque newMacroElements; for (std::deque::iterator elIter = macroElements.begin(); elIter != macroElements.end(); ++elIter) { // If the current mesh macro element should not be deleted, i.e., it is not a // member of the list of macro elements to be deleted, is is inserted to the new // macro element list. if (macros.find(*elIter) == macros.end()) newMacroElements.push_back(*elIter); } // And replace the macro element list with the new one. macroElements.clear(); macroElements = newMacroElements; // === For all macro elements to be deleted, delete them also to be neighbours === // === of some other macro elements. === for (std::set::iterator macroIt = macros.begin(); macroIt != macros.end(); ++macroIt) { // Go through all neighbours of the macro element and remove this macro element // to be neighbour of some other macro element. for (int i = 0; i <= dim; i++) { if ((*macroIt)->getNeighbour(i)) { for (int j = 0; j <= dim; j++) if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt) (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL); } else { // There is no neighbour at this edge, so we have to decrease the number // of edges in the mesh. nEdges--; } } // Decrease also the number of elements of the mesh. nLeaves--; nElements--; } // === Check now all the dofs, that have no owner anymore and therefore have === // === to be removed. === int nRemainDofs = 0; for (DofElMap::iterator dofsIt = dofsOwner.begin(); dofsIt != dofsOwner.end(); ++dofsIt) { bool deleteDof = true; for (std::set::iterator elIter = dofsIt->second.begin(); elIter != dofsIt->second.end(); ++elIter) { std::set::iterator mIt = macros.find(*elIter); if (mIt == macros.end()) { deleteDof = false; break; } } if (deleteDof) freeDof(const_cast(dofsIt->first), dofsPosIndex[dofsIt->first]); else nRemainDofs++; } // === Finally, remove the macro elements from memory. === for (std::set::iterator macroIt = macros.begin(); macroIt != macros.end(); ++macroIt) delete *macroIt; nVertices = nRemainDofs; } void Mesh::addDOFAdmin(DOFAdmin *localAdmin) { FUNCNAME("Mesh::addDOFAdmin()"); localAdmin->setMesh(this); std::vector::iterator dai = std::find(admin.begin(), admin.end(), localAdmin); TEST_EXIT(dai == admin.end()) ("admin %s is already associated to mesh %s\n", localAdmin->getName().c_str(), this->getName().c_str()); // if this will be required, see the untested code in revision < 224 // TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n"); admin.push_back(localAdmin); nDOFEl = 0; localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]); nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX); nDOFEl += getGeo(VERTEX) * nDOF[VERTEX]; if (dim > 1) { localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]); nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE); nDOFEl += getGeo(EDGE) * nDOF[EDGE]; } localAdmin->setNumberOfPreDOFs(CENTER,nDOF[CENTER]); nDOF[CENTER] += localAdmin->getNumberOfDOFs(CENTER); nDOFEl += nDOF[CENTER]; TEST_EXIT_DBG(nDOF[VERTEX] > 0)("no vertex dofs\n"); node[VERTEX] = 0; nNodeEl = getGeo(VERTEX); if (dim > 1) { node[EDGE] = nNodeEl; if (nDOF[EDGE] > 0) nNodeEl += getGeo(EDGE); } if (dim == 3) { localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]); nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE); nDOFEl += getGeo(FACE) * nDOF[FACE]; node[FACE] = nNodeEl; if (nDOF[FACE] > 0) nNodeEl += getGeo(FACE); } node[CENTER] = nNodeEl; if (nDOF[CENTER] > 0) nNodeEl += 1; } void Mesh::dofCompress() { FUNCNAME("Mesh::dofCompress()"); for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) { DOFAdmin* compressAdmin = admin[iadmin]; TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin); int size = compressAdmin->getSize(); if (size < 1 || compressAdmin->getUsedDOFs() < 1 || compressAdmin->getHoleCount() < 1) continue; newDOF.resize(size); compressAdmin->compress(newDOF); Flag fill_flag = (preserveCoarseDOFs ? Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING : Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag); while (elInfo) { elInfo->getElement()->newDOFFct1(compressAdmin); elInfo = stack.traverseNext(elInfo); } elInfo = stack.traverseFirst(this, -1, fill_flag); while (elInfo) { elInfo->getElement()->newDOFFct2(compressAdmin); elInfo = stack.traverseNext(elInfo); } newDOF.resize(0); } } DegreeOfFreedom *Mesh::getDof(GeoIndex position) { FUNCNAME("Mesh::getDof()"); TEST_EXIT_DBG(position >= CENTER && position <= FACE) ("unknown position %d\n", position); int ndof = getNumberOfDOFs(position); if (ndof <= 0) return NULL; DegreeOfFreedom *dof = new DegreeOfFreedom[ndof]; for (int i = 0; i < getNumberOfDOFAdmin(); i++) { const DOFAdmin *localAdmin = &getDOFAdmin(i); TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i); int n = localAdmin->getNumberOfDOFs(position); int n0 = localAdmin->getNumberOfPreDOFs(position); TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof); for (int j = 0; j < n; j++) dof[n0 + j] = const_cast(localAdmin)->getDOFIndex(); } return dof; } DegreeOfFreedom **Mesh::createDofPtrs() { FUNCNAME("Mesh::createDofPtrs()"); if (nNodeEl <= 0) return NULL; DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl]; for (int i = 0; i < nNodeEl; i++) ptrs[i] = NULL; return ptrs; } void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs) { FUNCNAME("Mesh::freeDOFPtrs()"); TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n"); if (nNodeEl <= 0) return; delete [] ptrs; } const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec lnDOF) { FUNCNAME("Mesh::createDOFAdmin()"); DOFAdmin *localAdmin = new DOFAdmin(this, lname); for (int i = 0; i < dim + 1; i++) localAdmin->setNumberOfDOFs(i, lnDOF[i]); addDOFAdmin(localAdmin); return localAdmin; } const DOFAdmin* Mesh::getVertexAdmin() const { const DOFAdmin *localAdmin = NULL; for (int i = 0; i < static_cast(admin.size()); i++) { if (admin[i]->getNumberOfDOFs(VERTEX)) { if (!localAdmin) localAdmin = admin[i]; else if (admin[i]->getSize() < localAdmin->getSize()) localAdmin = admin[i]; } } return localAdmin; } void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position) { FUNCNAME("Mesh::freeDof()"); TEST_EXIT_DBG(position >= CENTER && position <= FACE) ("unknown position %d\n", position); int ndof = nDOF[position]; if (ndof) { if (!dof) { MSG("dof = NULL, but ndof = %d\n", ndof); return; } } else { if (dof) MSG("dof != NULL, but ndof = 0\n"); return; } TEST_EXIT_DBG(ndof <= MAX_DOF) ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF); for (int i = 0; i < static_cast(admin.size()); i++) { DOFAdmin *localAdmin = admin[i]; int n = localAdmin->getNumberOfDOFs(position); int n0 = localAdmin->getNumberOfPreDOFs(position); TEST_EXIT_DBG(n + n0 <= ndof) ("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof); for (int j = 0; j < n; j++) localAdmin->freeDofIndex(dof[n0 + j]); } delete [] dof; } void Mesh::freeElement(Element* el) { freeDOFPtrs(const_cast(el->getDOF())); delete el; } Element* Mesh::createNewElement(Element *parent) { FUNCNAME("Mesh::createNewElement()"); TEST_EXIT_DBG(elementPrototype)("no element prototype\n"); Element *el = parent ? parent->clone() : elementPrototype->clone(); if (!parent && elementDataPrototype) el->setElementData(elementDataPrototype->clone()); else el->setElementData(NULL); // must be done in ElementData::refineElementData() return el; } ElInfo* Mesh::createNewElInfo() { FUNCNAME("Mesh::createNewElInfo()"); switch(dim) { case 1: return new ElInfo1d(this); break; case 2: return new ElInfo2d(this); break; case 3: return new ElInfo3d(this); break; default: ERROR_EXIT("invalid dim\n"); return NULL; } } bool Mesh::findElInfoAtPoint(const WorldVector& xy, ElInfo *el_info, DimVec& bary, const MacroElement *start_mel, const WorldVector *xy0, double *sp) { static const MacroElement *mel = NULL; DimVec lambda(dim, NO_INIT); ElInfo *mel_info = createNewElInfo(); if (start_mel != NULL) mel = start_mel; else if (mel == NULL || mel->getElement()->getMesh() != this) mel = *(macroElements.begin()); mel_info->setFillFlag(Mesh::FILL_COORDS); g_xy = &xy; g_xy0 = xy0; g_sp = sp; mel_info->fillMacroInfo(mel); // We have the care about not to visite a macro element twice. In this case the // function would end up in an infinite loop. If a macro element is visited a // second time, what can happen with periodic boundary conditions, the point is // not within the mesh! std::set macrosVisited; macrosVisited.insert(mel->getIndex()); int k; while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) { if (mel->getNeighbour(k)) { mel = mel->getNeighbour(k); if (macrosVisited.count(mel->getIndex())) return false; macrosVisited.insert(mel->getIndex()); mel_info->fillMacroInfo(mel); continue; } break; } /* now, descend in tree to find leaf element at point */ bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info); for (int i = 0; i <= dim; i++) bary[i] = final_lambda[i]; delete mel_info; return inside; } bool Mesh::findElementAtPoint(const WorldVector& xy, Element **elp, DimVec& bary, const MacroElement *start_mel, const WorldVector *xy0, double *sp) { ElInfo *el_info = createNewElInfo(); int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp); *elp = el_info->getElement(); delete el_info; return val; } bool Mesh::findElementAtPointRecursive(ElInfo *el_info, const DimVec& lambda, int outside, ElInfo* final_el_info) { FUNCNAME("Mesh::findElementAtPointRecursive()"); Element *el = el_info->getElement(); DimVec c_lambda(dim, NO_INIT); int inside; int ichild, c_outside; if (el->isLeaf()) { *final_el_info = *el_info; if (outside < 0) { for (int i = 0; i <= dim; i++) final_lambda[i] = lambda[i]; return true; } else { /* outside */ if (g_xy0) { /* find boundary point of [xy0, xy] */ el_info->worldToCoord(*(g_xy0), &c_lambda); double s = lambda[outside] / (lambda[outside] - c_lambda[outside]); for (int i = 0; i <= dim; i++) final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i]; if (g_sp) *(g_sp) = s; if (dim == 3) MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s); return false; /* ??? */ } else { return false; } } } ElInfo *c_el_info = createNewElInfo(); if (dim == 1) { if (lambda[0] >= lambda[1]) { c_el_info->fillElInfo(0, el_info); if (outside >= 0) { outside = el_info->worldToCoord(*(g_xy), &c_lambda); TEST_EXIT(outside == 0)("point outside domain\n"); } else { c_lambda[0] = lambda[0] - lambda[1]; c_lambda[1] = 2.0 * lambda[1]; } } else { c_el_info->fillElInfo(1, el_info); if (outside >= 0) { outside = el_info->worldToCoord(*(g_xy), &c_lambda); TEST_EXIT(outside == 0)("point outside domain\n"); } else { c_lambda[1] = lambda[1] - lambda[0]; c_lambda[0] = 2.0 * lambda[0]; } } } /* DIM == 1 */ if (dim == 2) { if (lambda[0] >= lambda[1]) { c_el_info->fillElInfo(0, el_info); if (el->isNewCoordSet()) { outside = c_el_info->worldToCoord(*(g_xy), &c_lambda); TEST_EXIT(outside == 0)("outside curved boundary child 0\n"); } else { c_lambda[0] = lambda[2]; c_lambda[1] = lambda[0] - lambda[1]; c_lambda[2] = 2.0 * lambda[1]; } } else { c_el_info->fillElInfo(1, el_info); if (el->isNewCoordSet()) { outside = c_el_info->worldToCoord(*(g_xy), &c_lambda); TEST_EXIT(outside == 0)("outside curved boundary child 1\n"); } else { c_lambda[0] = lambda[1] - lambda[0]; c_lambda[1] = lambda[2]; c_lambda[2] = 2.0 * lambda[0]; } } } /* DIM == 2 */ if (dim == 3) { if (el->isNewCoordSet()) { if (lambda[0] >= lambda[1]) ichild = 0; else ichild = 1; c_el_info->fillElInfo(ichild, el_info); c_outside = c_el_info->worldToCoord(*(g_xy), &c_lambda); if (c_outside>=0) { /* test is other child is better... */ DimVec c_lambda2(dim, NO_INIT); ElInfo *c_el_info2 = createNewElInfo(); c_el_info2->fillElInfo(1-ichild, el_info); int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2); MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n", ichild, c_outside, c_lambda[0], c_lambda[1], c_lambda[2], c_lambda[3]); MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n", 1 - ichild, c_outside2, c_lambda2[0], c_lambda2[1], c_lambda2[2], c_lambda2[3]); if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) { for (int i = 0; i <= dim; i++) c_lambda[i] = c_lambda2[i]; c_outside = c_outside2; *c_el_info = *c_el_info2; ichild = 1 - ichild; } delete c_el_info2; } outside = c_outside; } else { /* no new_coord */ if (lambda[0] >= lambda[1]) { c_el_info->fillElInfo(0, el_info); c_lambda[0] = lambda[0] - lambda[1]; c_lambda[1] = lambda[Tetrahedron::childVertex[(dynamic_cast(el_info))-> getType()][0][1]]; c_lambda[2] = lambda[Tetrahedron::childVertex[(dynamic_cast(el_info))-> getType()][0][2]]; c_lambda[3] = 2.0 * lambda[1]; } else { c_el_info->fillElInfo(1, el_info); c_lambda[0] = lambda[1] - lambda[0]; c_lambda[1] = lambda[Tetrahedron::childVertex[(dynamic_cast(el_info))-> getType()][1][1]]; c_lambda[2] = lambda[Tetrahedron::childVertex[(dynamic_cast(el_info))-> getType()][1][2]]; c_lambda[3] = 2.0 * lambda[0]; } } } /* DIM == 3 */ inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info); delete c_el_info; return inside; } bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, const FiniteElemSpace* feSpace, WorldVector& coords) { DimVec* baryCoords; bool found = false; TraverseStack stack; std::vector dofVec(feSpace->getBasisFcts()->getNumber()); ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while (elInfo) { feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec); for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) { if (dofVec[i] == dof) { baryCoords = feSpace->getBasisFcts()->getCoords(i); elInfo->coordToWorld(*baryCoords, coords); found = true; break; } } if (found) break; elInfo = stack.traverseNext(elInfo); } return found; } void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace, DOFVector >& coords) { const BasisFunction* basFcts = feSpace->getBasisFcts(); int nBasFcts = basFcts->getNumber(); std::vector dofVec(nBasFcts); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while (elInfo) { basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec); for (int i = 0; i < nBasFcts; i++) { DimVec *baryCoords = basFcts->getCoords(i); elInfo->coordToWorld(*baryCoords, coords[dofVec[i]]); } elInfo = stack.traverseNext(elInfo); } } void Mesh::setDiameter(const WorldVector& w) { diam = w; } void Mesh::setDiameter(int i, double w) { diam[i] = w; } void Mesh::serialize(std::ostream &out) { serializedDOFs.clear(); out << name << "\n"; SerUtil::serialize(out, dim); SerUtil::serialize(out, nVertices); SerUtil::serialize(out, nEdges); SerUtil::serialize(out, nLeaves); SerUtil::serialize(out, nElements); SerUtil::serialize(out, nFaces); SerUtil::serialize(out, maxEdgeNeigh); diam.serialize(out); SerUtil::serialize(out, preserveCoarseDOFs); SerUtil::serialize(out, nDOFEl); nDOF.serialize(out); SerUtil::serialize(out, nNodeEl); node.serialize(out); // === Write admins. === int size = static_cast(admin.size()); SerUtil::serialize(out, size); for (int i = 0; i < size; i++) admin[i]->serialize(out); // === Write macroElements. === size = static_cast(macroElements.size()); SerUtil::serialize(out, size); for (int i = 0; i < size; i++) macroElements[i]->serialize(out); SerUtil::serialize(out, elementIndex); SerUtil::serialize(out, initialized); // === Write periodic associations. === int mapSize = periodicAssociations.size(); SerUtil::serialize(out, mapSize); for (std::map::iterator it = periodicAssociations.begin(); it != periodicAssociations.end(); ++it) { BoundaryType b = it->first; // Check which DOFAdmin is used for the current VertexVector we want to serialize. int ithAdmin = -1; for (int i = 0; i < static_cast(admin.size()); i++) { if (admin[i] == it->second->getAdmin()) { ithAdmin = i; break; } } TEST_EXIT(ithAdmin >= 0) ("No DOFAdmin found for serialization of periodic associations!\n"); SerUtil::serialize(out, b); SerUtil::serialize(out, ithAdmin); it->second->serialize(out); } serializedDOFs.clear(); } void Mesh::deserialize(std::istream &in) { FUNCNAME("Mesh::deserialize()"); serializedDOFs.clear(); in >> name; in.get(); int oldVal = dim; SerUtil::deserialize(in, dim); TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n"); SerUtil::deserialize(in, nVertices); SerUtil::deserialize(in, nEdges); SerUtil::deserialize(in, nLeaves); SerUtil::deserialize(in, nElements); SerUtil::deserialize(in, nFaces); SerUtil::deserialize(in, maxEdgeNeigh); diam.deserialize(in); SerUtil::deserialize(in, preserveCoarseDOFs); oldVal = nDOFEl; SerUtil::deserialize(in, nDOFEl); TEST_EXIT_DBG(oldVal == 0 || nDOFEl == oldVal)("Invalid nDOFEl!\n"); nDOF.deserialize(in); oldVal = nNodeEl; SerUtil::deserialize(in, nNodeEl); TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n"); node.deserialize(in); // === Read admins. === int size; SerUtil::deserialize(in, size); admin.resize(size, NULL); for (int i = 0; i < size; i++) { if (!admin[i]) admin[i] = new DOFAdmin(this); admin[i]->deserialize(in); } SerUtil::deserialize(in, size); std::vector< std::vector > neighbourIndices(size); deleteMeshStructure(); // All macro elements are stored in the list \ref macroElements with continous // index from 0 to n - 1. But the macro element index may be different and may // contain holes (for example when macro elements were removed because of domain // decomposition based parallelization. Therefore we create a temporary map // from macro element indices to the continous index of \ref macroElements. This // will be used later to find the correct neighbours of the macro elements. std::map elIndexVecIndex; macroElements.resize(size); for (int i = 0; i < size; i++) { macroElements[i] = new MacroElement(dim); macroElements[i]->writeNeighboursTo(&(neighbourIndices[i])); macroElements[i]->deserialize(in); elIndexVecIndex[macroElements[i]->getIndex()] = i; } SerUtil::deserialize(in, elementIndex); SerUtil::deserialize(in, initialized); // set neighbour pointer in macro elements int nNeighbour = getGeo(NEIGH); for (int i = 0; i < static_cast(macroElements.size()); i++) { for (int j = 0; j < nNeighbour; j++) { int index = neighbourIndices[i][j]; if (index != -1) { TEST_EXIT_DBG(elIndexVecIndex.count(index) == 1) ("Cannot find correct index from neighbouring macro element!\n"); index = elIndexVecIndex[index]; macroElements[i]->setNeighbour(j, macroElements[index]); } else { macroElements[i]->setNeighbour(j, NULL); } } } // set mesh pointer in elements TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER); while (elInfo) { elInfo->getElement()->setMesh(this); elInfo = stack.traverseNext(elInfo); } serializedDOFs.clear(); /// === Read periodic assoications. === int mapSize = 0; SerUtil::deserialize(in, mapSize); for (int i = 0; i < mapSize; i++) { BoundaryType b = 0; int ithAdmin = 0; SerUtil::deserialize(in, b); SerUtil::deserialize(in, ithAdmin); VertexVector *tmpvec = new VertexVector(admin[ithAdmin], ""); tmpvec->deserialize(in); periodicAssociations[b] = tmpvec; } } void Mesh::initialize() { FUNCNAME("Mesh::initialize()"); std::string macroFilename(""); std::string valueFilename(""); std::string periodicFilename(""); int check = 1; GET_PARAMETER(0, name + "->macro file name", ¯oFilename); GET_PARAMETER(0, name + "->value file name", &valueFilename); GET_PARAMETER(0, name + "->periodic file", &periodicFilename); GET_PARAMETER(0, name + "->check", "%d", &check); GET_PARAMETER(0, name + "->preserve coarse dofs", "%d", &preserveCoarseDOFs); if (macroFilename.length()) { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS if (checkParallelMacroFile(macroFilename, periodicFilename, check)) { std::string newPeriodicFilename = ""; if (periodicFilename != "") newPeriodicFilename = periodicFilename + ".tmp"; macroFileInfo = MacroReader::readMacro(macroFilename + ".tmp", this, newPeriodicFilename, check); } else { macroFileInfo = MacroReader::readMacro(macroFilename, this, periodicFilename, check); } #else // In sequentiel computations just read the macro file to the mesh. macroFileInfo = MacroReader::readMacro(macroFilename, this, periodicFilename, check); #endif // If there is no value file which should be written, we can delete // the information of the macro file. if (!valueFilename.length()) clearMacroFileInfo(); } initialized = true; } #ifdef HAVE_PARALLEL_DOMAIN_AMDIS bool Mesh::checkParallelMacroFile(std::string macroFilename, std::string periodicFilename, int check) { FUNCNAME("Mesh::checkParallelMacroFile()"); // In parallel computations, first the mesh must be checked if it is refined // in an apropriate way. TEST_EXIT(admin.size() == 1)("Not yet implemented!\n"); // === Create a temporary mesh and load the macro file to it. === Mesh testMesh(name, dim); DOFAdmin *localAdmin = new DOFAdmin(&testMesh, admin[0]->getName()); localAdmin->setNumberOfDOFs(admin[0]->getNumberOfDOFs()); testMesh.addDOFAdmin(localAdmin); MacroInfo *testMacroInfo = MacroReader::readMacro(macroFilename, &testMesh, periodicFilename, check); // === Check the mesh structure. === int nMacroElements = 0; int elType = -1; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(&testMesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { if (elType == -1) { elType = elInfo->getType(); } else { TEST_EXIT(elType == elInfo->getType()) ("All elements in mesh must have the same element type!\n"); } nMacroElements++; elInfo = stack.traverseNext(elInfo); } // === Calculate the number of global refinements such that the new mesh === // === would fulfill all requirements. === // There should be at least 10 macro Elements per processor, therefore: // nMacroElements * 2^gr >= nProcs * 10 // => gr = log_2(nProcs * 10 / nMacroElements) double tmp = 10.0 * MPI::COMM_WORLD.Get_size() / nMacroElements; int nrRefine = ceil(log(tmp) / log(2)); if (nrRefine < 0) nrRefine = 0; if (dim == 3) { int newElType = (elType + nrRefine) % 3; switch (newElType) { case 1: if (nrRefine > 0) nrRefine--; else nrRefine = 2; break; case 2: nrRefine++; break; } TEST_EXIT((elType + nrRefine) % 3 == 0)("This should not happen!\n"); } // === If we do not need to refine the mesh, return back. === if (nrRefine == 0) return false; // === If macro weights are explicitly given, we cannot change the mesh. === std::string macroWeightsFilename = ""; GET_PARAMETER(0, name + "->macro weights", ¯oWeightsFilename); if (macroWeightsFilename != "") { ERROR_EXIT("Should not happen!\n"); } // === Rank 0 creates a new mesh file. === if (MPI::COMM_WORLD.Get_rank() == 0) { RefinementManager *refManager; if (dim == 2) refManager = new RefinementManager2d(); else refManager = new RefinementManager3d(); refManager->globalRefine(&testMesh, nrRefine); delete refManager; Lagrange* basFcts = Lagrange::getLagrange(dim, 1); FiniteElemSpace *feSpace = FiniteElemSpace::provideFeSpace(localAdmin, basFcts, &testMesh, "tmp"); DataCollector dc(feSpace); MacroWriter::writeMacro(&dc, (macroFilename + ".tmp").c_str()); if (periodicFilename != "") MacroWriter::writePeriodicFile(&dc, (periodicFilename + ".tmp").c_str()); } // === All ranks must wait until rank 0 has created the new macro mesh file. === MPI::COMM_WORLD.Barrier(); // === We have refined the mesh, so reduce the number of global refinements. === int globalRefinements = 0; GET_PARAMETER(0, name + "->global refinements", "%d", &globalRefinements); if (globalRefinements < nrRefine) globalRefinements = 0; else globalRefinements -= nrRefine; std::stringstream oss; oss << globalRefinements; ADD_PARAMETER(0, name + "->global refinements", oss.str().c_str()); // === Print a note to the screen that another mesh file will be used. === MSG("The macro mesh file \"%s\" was refined %d times and stored to file \"%s\".\n", macroFilename.c_str(), nrRefine, (macroFilename + ".tmp").c_str()); return true; } #endif bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) { std::map::iterator it; std::map::iterator end = periodicAssociations.end(); for (it = periodicAssociations.begin(); it != end; ++it) if ((*(it->second))[dof1] == dof2) return true; return false; } bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) { std::vector associatedToDOF1; std::map::iterator it; std::map::iterator end = periodicAssociations.end(); DegreeOfFreedom dof, assDOF; associatedToDOF1.push_back(dof1); for (it = periodicAssociations.begin(); it != end; ++it) { int size = static_cast(associatedToDOF1.size()); for (int i = 0; i < size; i++) { dof = associatedToDOF1[i]; assDOF = (*(it->second))[dof]; if (assDOF == dof2) { return true; } else { if (assDOF != dof) associatedToDOF1.push_back(assDOF); } } } return false; } void Mesh::clearMacroFileInfo() { macroFileInfo->clear(); delete macroFileInfo; macroFileInfo = NULL; } void Mesh::computeMatrixFillin(const FiniteElemSpace *feSpace, std::vector &nnzInRow, int &overall, int &average) { std::map dofCounter; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL); ElementDofIterator elDofIter(feSpace); while (elInfo) { elDofIter.reset(elInfo->getElement()); do { DegreeOfFreedom dof = elDofIter.getDof(); if (dofCounter.count(dof) == 0) { dofCounter[dof] = 1; } else { dofCounter[dof]++; } } while (elDofIter.next()); elInfo = stack.traverseNext(elInfo); } overall = 0; for (std::map::iterator it = dofCounter.begin(); it != dofCounter.end(); ++it) overall += it->second * 15; } int Mesh::calcMemoryUsage() { int result = sizeof(Mesh); result += nDOFEl; for (int i = 0; i < static_cast(admin.size()); i++) { result += admin[i]->calcMemoryUsage(); result += admin[i]->getUsedSize() * sizeof(DegreeOfFreedom); } for (int i = 0; i < static_cast(macroElements.size()); i++) result += macroElements[i]->calcMemoryUsage(); return result; } void Mesh::deleteMeshStructure() { Element::deletedDOFs.clear(); for (std::deque::const_iterator it = macroElements.begin(); it != macroElements.end(); ++it) { (*it)->getElement()->deleteElementDOFs(); delete *it; } Element::deletedDOFs.clear(); } }