#include #include #include #include "time.h" #include "AdaptStationary.h" #include "AdaptInstationary.h" #include "FiniteElemSpace.h" #include "ElementData.h" #include "MacroElement.h" #include "MacroReader.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" 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 ; // used in mg methods // 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 Mesh::serializedDOFs; struct delmem { DegreeOfFreedom* ptr; int len; }; Mesh::Mesh(const 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), 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() { Element::deletedDOFs.clear(); for (std::deque::const_iterator it = macroElements.begin(); it != macroElements.end(); ++it) { (*it)->getElement()->deleteElementDOFs(); delete *it; } Element::deletedDOFs.clear(); 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::vector& macros) { FUNCNAME("Mesh::removeMacroElement()"); TEST_EXIT(dim == 2)("Not yet implemented!\n"); // Map that stores for each dof pointer (which may have a list of dofs) // all macro element indices that own the dof. std::map > dofsOwner; // Determine all dof owner macro elements. for (std::deque::iterator macroIt = macroElements.begin(); macroIt != macroElements.end(); ++macroIt) { Element *el = (*macroIt)->getElement(); for (int i = 0; i < 3; i++) dofsOwner[el->getDOF(i)].insert(*macroIt); } // Remove all the given macro elements. for (std::vector::iterator macroIt = macros.begin(); macroIt != macros.end(); ++macroIt) { bool found = false; // Remove the macro element from mesh's list of all macro elements. for (std::deque::iterator it = macroElements.begin(); it != macroElements.end(); ++it) { if (*it == *macroIt) { macroElements.erase(it, it + 1); found = true; break; } } TEST_EXIT(found)("Cannot find MacroElement that should be removed!\n"); // 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--; } } nLeaves--; nElements--; // Remove this macro element from the dof owner list. for (std::map >::iterator dofsIt = dofsOwner.begin(); dofsIt != dofsOwner.end(); ++dofsIt) { std::set::iterator mIt = dofsIt->second.find(*macroIt); if (mIt != dofsIt->second.end()) { dofsIt->second.erase(mIt); } } // And remove the macro element from memory delete *macroIt; } int nRemainDofs = 0; // Check now all the dofs, that have no owner anymore and therefore have to // be removed. for (std::map >::iterator dofsIt = dofsOwner.begin(); dofsIt != dofsOwner.end(); ++dofsIt) { if (dofsIt->second.size() == 0) { freeDOF(const_cast(dofsIt->first), VERTEX); } else { nRemainDofs++; } } nVertices = nRemainDofs; } int Mesh::traverse(int level, Flag flag, int (*el_fct)(ElInfo*)) { FUNCNAME("Mesh::traverse()"); std::deque::iterator mel; ElInfoStack elInfoStack(this); ElInfo* elinfo = elInfoStack.getNextElement(); Traverse tinfo(this, flag, level, el_fct); int sum = 0; 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); } for (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()"); 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()"); Flag fill_flag; for (int iadmin = 0; iadmin < static_cast(admin.size()); iadmin++) { compressAdmin = admin[iadmin]; TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin); int size = compressAdmin->getSize(); if (size < 1) continue; if (compressAdmin->getUsedDOFs() < 1) continue; if (compressAdmin->getHoleCount() < 1) continue; newDOF.resize(size); compressAdmin->compress(newDOF); if (preserveCoarseDOFs) { fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING; } else { fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING; } traverse(-1, fill_flag, newDOFFct1); traverse(-1, fill_flag, newDOFFct2); 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=NULL\n"); if (nNodeEl <= 0) return; delete [] ptrs; } const DOFAdmin *Mesh::createDOFAdmin(const 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() { 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 = NULL; 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); int k; while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) { if (mel->getNeighbour(k)) { mel = mel->getNeighbour(k); 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); if (outside >= 0) ERROR("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); if (outside >= 0) ERROR("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); if (outside >= 0) { ERROR("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); if (outside >= 0) { ERROR("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); int c_outside2; ElInfo *c_el_info2 = createNewElInfo(); c_el_info2->fillElInfo(1-ichild, el_info); 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; } void Mesh::setDiameter(const WorldVector& w) { diam = w; } void Mesh::setDiameter(int i, double w) { diam[i] = w; } int Mesh::newDOFFct1(ElInfo* ei) { ei->getElement()->newDOFFct1(compressAdmin); return 0; } int Mesh::newDOFFct2(ElInfo* ei) { ei->getElement()->newDOFFct2(compressAdmin); return 0; } void Mesh::serialize(std::ostream &out) { serializedDOFs.clear(); // write name out << name << "\n"; // write dim out.write(reinterpret_cast(&dim), sizeof(int)); // write nVertices out.write(reinterpret_cast(&nVertices), sizeof(int)); // write nEdges out.write(reinterpret_cast(&nEdges), sizeof(int)); // write nLeaves out.write(reinterpret_cast(&nLeaves), sizeof(int)); // write nElements out.write(reinterpret_cast(&nElements), sizeof(int)); // write nFaces out.write(reinterpret_cast(&nFaces), sizeof(int)); // write maxEdgeNeigh out.write(reinterpret_cast(&maxEdgeNeigh), sizeof(int)); // write diam diam.serialize(out); // write preserveCoarseDOFs out.write(reinterpret_cast(&preserveCoarseDOFs), sizeof(bool)); // write nDOFEl out.write(reinterpret_cast(&nDOFEl), sizeof(int)); // write nDOF nDOF.serialize(out); // write nNodeEl out.write(reinterpret_cast(&nNodeEl), sizeof(int)); // write node node.serialize(out); // write admins int size = static_cast(admin.size()); out.write(reinterpret_cast(&size), sizeof(int)); for (int i = 0; i < size; i++) admin[i]->serialize(out); // write macroElements size = static_cast(macroElements.size()); out.write(reinterpret_cast(&size), sizeof(int)); for (int i = 0; i < size; i++) macroElements[i]->serialize(out); // write elementIndex out.write(reinterpret_cast(&elementIndex), sizeof(int)); // write initialized out.write(reinterpret_cast(&initialized), sizeof(bool)); serializedDOFs.clear(); } void Mesh::deserialize(std::istream &in) { serializedDOFs.clear(); // read name in >> name; in.get(); // read dim int oldVal = dim; in.read(reinterpret_cast(&dim), sizeof(int)); TEST_EXIT_DBG((oldVal == 0) || (dim == oldVal))("invalid dimension\n"); // read nVertices in.read(reinterpret_cast(&nVertices), sizeof(int)); // read nEdges in.read(reinterpret_cast(&nEdges), sizeof(int)); // read nLeaves in.read(reinterpret_cast(&nLeaves), sizeof(int)); // read nElements in.read(reinterpret_cast(&nElements), sizeof(int)); // read nFaces in.read(reinterpret_cast(&nFaces), sizeof(int)); // read maxEdgeNeigh in.read(reinterpret_cast(&maxEdgeNeigh), sizeof(int)); // diam diam.deserialize(in); // read preserveCoarseDOFs in.read(reinterpret_cast(&preserveCoarseDOFs), sizeof(bool)); // read nDOFEl oldVal = nDOFEl; in.read(reinterpret_cast(&nDOFEl), sizeof(int)); TEST_EXIT_DBG((oldVal == 0) || (nDOFEl == oldVal))("invalid nDOFEl\n"); // read nDOF nDOF.deserialize(in); // read nNodeEl oldVal = nNodeEl; in.read(reinterpret_cast(&nNodeEl), sizeof(int)); TEST_EXIT_DBG((oldVal == 0) || (nNodeEl == oldVal))("invalid nNodeEl\n"); // read node node.deserialize(in); // read admins int size; in.read(reinterpret_cast(&size), sizeof(int)); admin.resize(size, NULL); for (int i = 0; i < size; i++) { if (!admin[i]) admin[i] = new DOFAdmin(this); admin[i]->deserialize(in); } // read macroElements in.read(reinterpret_cast(&size), sizeof(int)); std::vector< std::vector > neighbourIndices(size); for (int i = 0; i < static_cast(macroElements.size()); i++) { if (macroElements[i]) delete macroElements[i]; } macroElements.resize(size); for (int i = 0; i < size; i++) { macroElements[i] = new MacroElement(dim); macroElements[i]->writeNeighboursTo(&(neighbourIndices[i])); macroElements[i]->deserialize(in); } // read elementIndex in.read(reinterpret_cast(&elementIndex), sizeof(int)); // read initialized in.read(reinterpret_cast(&initialized), sizeof(bool)); // set neighbour pointer in macro elements int neighs = getGeo(NEIGH); for (int i = 0; i < static_cast(macroElements.size()); i++) { for (int j = 0; j < neighs; j++) { int index = neighbourIndices[i][j]; if(index != -1) { 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(); } void Mesh::initialize() { std::string macroFilename(""); std::string valueFilename(""); std::string periodicFile(""); 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", &periodicFile); GET_PARAMETER(0, name + "->check", "%d", &check); GET_PARAMETER(0, name + "->preserve coarse dofs", "%d", &preserveCoarseDOFs); if (macroFilename.length()) { macroFileInfo = MacroReader::readMacro(macroFilename.c_str(), this, periodicFile == "" ? NULL : periodicFile.c_str(), check); // 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; } 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; } 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; } }