#include "Element.h" #include "DOFAdmin.h" #include "Mesh.h" #include "CoarseningManager.h" #include "FixVec.h" #include "ElementRegion_ED.h" namespace AMDiS { std::map Element::deletedDOFs; int Element::getRegion() const { if (!elementData) return -1; ElementRegion_ED* red = dynamic_cast(elementData->getElementData(ELEMENT_REGION)); if (red) return red->getRegion(); return -1; } void Element::setDOFPtrs() { FUNCNAME("Element::setDOFPtrs()"); TEST_EXIT_DBG(mesh)("no mesh!\n"); dof = mesh->createDOFPtrs(); } Element::Element(Mesh *aMesh) { mesh = aMesh; index = mesh ? mesh->getNextElementIndex() : -1; child[0] = NULL; child[1] = NULL; newCoord = NULL; elementData = NULL; mark = 0; if (mesh) { setDOFPtrs(); } else { mesh = NULL; } } // call destructor through Mesh::freeElement !!! Element::~Element() { if (child[0]) delete child[0]; if (child[1]) delete child[1]; if (newCoord) delete newCoord; if (elementData) { elementData->deleteDecorated(); delete elementData; } } bool Element::deleteElementData(int typeID) { FUNCNAME("Element::deleteElementData()"); if (elementData) { if (elementData->isOfType(typeID)) { ElementData *tmp = elementData; elementData = elementData->getDecorated(); delete tmp; tmp = NULL; return true; } else { return elementData->deleteDecorated(typeID); } } return false; } void Element::deleteElementDOFs() { int dim = mesh->getDim(); int j = 0; for (int pos = 0; pos <= dim; pos++) { GeoIndex position = INDEX_OF_DIM(pos, dim); int ndof = 0; for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) { ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position); } if (ndof > 0) { for (int i = 0; i < mesh->getGeo(position); i++) { if (dof[j]) { if (deletedDOFs.count(dof[j]) == 0) { deletedDOFs[dof[j]] = true; delete [] dof[j]; } } j++; } } } FREE_MEMORY(dof, DegreeOfFreedom*, mesh->getNumberOfNodes()); if (child[0]) child[0]->deleteElementDOFs(); if (child[1]) child[1]->deleteElementDOFs(); } Element* Element::cloneWithDOFs() { Element *el; if (isLine()) { el = new Line(NULL); } else if (isTriangle()) { el = new Triangle(NULL); } else { el = new Tetrahedron(NULL); } el->mesh = mesh; el->index = index; el->mark = mark; if (newCoord) { WorldVector *nc = new WorldVector(); *nc = *newCoord; el->newCoord = nc; } /* =========== And here we clone the DOFs =========== */ el->dof = new DegreeOfFreedom*[mesh->getNumberOfNodes()]; int dim = mesh->getDim(); int j = 0; for (int pos = 0; pos <= dim; pos++) { GeoIndex position = INDEX_OF_DIM(pos, dim); int ndof = 0; for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position); if (ndof > 0) { for (int i = 0; i < mesh->getGeo(position); i++) { if (dof[j] != NULL) { if (Mesh::serializedDOFs[dof[j][0]] == NULL) { el->dof[j] = new DegreeOfFreedom[ndof]; for (int k = 0; k < ndof; k++) el->dof[j][k] = dof[j][k]; Mesh::serializedDOFs[dof[j][0]] = el->dof[j]; } else { el->dof[j] = Mesh::serializedDOFs[dof[j][0]]; } } else { el->dof[j] = NULL; } j++; } } } /* =========== And clone the children ============= */ if (child[0]) { el->child[0] = child[0]->cloneWithDOFs(); } if (child[1]) { el->child[1] = child[1]->cloneWithDOFs(); } return el; } /****************************************************************************/ /* ATTENTION: */ /* new_dof_fct() destroys new_dof !!!!!!!!!! */ /* should be used only at the end of dof_compress()!!!!! */ /****************************************************************************/ /* CHANGE_DOFS_1 changes old dofs to NEGATIVE new dofs */ #define CHANGE_DOFS_1(el) \ ldof = el->dof[n0 + i] + nd0; \ for (j = 0; j < nd; j++) { \ if ((k = ldof[j]) >= 0) { \ /* do it only once! (dofs are visited more than once) */ \ ldof[j] = - admin->getMesh()->newDOF[k] - 1; \ } } /* CHANGE_DOFS_2 changes NEGATIVE new dofs to POSITIVE */ #define CHANGE_DOFS_2(el) \ ldof = el->dof[n0+i] + nd0; \ for (j = 0; j < nd; j++) { \ if ((k = ldof[j]) < 0) { \ /* do it only once! (dofs are visited more than once) */ \ ldof[j] = - k - 1; \ } } void Element::newDOFFct1(const DOFAdmin* admin) { int j, k, n0, nd, nd0; DegreeOfFreedom *ldof; int vertices = mesh->getGeo(VERTEX); int edges = mesh->getGeo(EDGE); int faces = mesh->getGeo(FACE); if ((nd = admin->getNumberOfDOFs(VERTEX))) { nd0 = admin->getNumberOfPreDOFs(VERTEX); n0 = admin->getMesh()->getNode(VERTEX); for (int i = 0; i < vertices; i++) { CHANGE_DOFS_1(this); } } if (mesh->getDim() > 1) { if ((nd = admin->getNumberOfDOFs(EDGE))) { nd0 = admin->getNumberOfPreDOFs(EDGE); n0 = admin->getMesh()->getNode(EDGE); for (int i = 0; i < edges; i++) { CHANGE_DOFS_1(this); } } } if (mesh->getDim() == 3) { if ((nd = admin->getNumberOfDOFs(FACE))) { nd0 = admin->getNumberOfPreDOFs(FACE); n0 = admin->getMesh()->getNode(FACE); for (int i = 0; i < faces; i++) { CHANGE_DOFS_1(this); } } } if ((nd = admin->getNumberOfDOFs(CENTER))) { nd0 = admin->getNumberOfPreDOFs(CENTER); n0 = admin->getMesh()->getNode(CENTER); int i = 0; /* only one center */ CHANGE_DOFS_1(this); } } void Element::newDOFFct2(const DOFAdmin* admin) { int i, j, k, n0, nd, nd0; DegreeOfFreedom *ldof; int vertices = mesh->getGeo(VERTEX); int edges = mesh->getGeo(EDGE); int faces = mesh->getGeo(FACE); if ((nd = admin->getNumberOfDOFs(VERTEX))) { nd0 = admin->getNumberOfPreDOFs(VERTEX); n0 = admin->getMesh()->getNode(VERTEX); for (i = 0; i < vertices; i++) { CHANGE_DOFS_2(this); } } if (mesh->getDim() > 1) { if ((nd = admin->getNumberOfDOFs(EDGE))) { nd0 = admin->getNumberOfPreDOFs(EDGE); n0 = admin->getMesh()->getNode(EDGE); for (i = 0; i < edges; i++) { CHANGE_DOFS_2(this); } } } if (mesh->getDim() == 3) { if ((nd = admin->getNumberOfDOFs(FACE))) { nd0 = admin->getNumberOfPreDOFs(FACE); n0 = admin->getMesh()->getNode(FACE); for (i = 0; i < faces; i++) { CHANGE_DOFS_2(this); } } } if ((nd = admin->getNumberOfDOFs(CENTER))) { nd0 = admin->getNumberOfPreDOFs(CENTER); n0 = admin->getMesh()->getNode(CENTER); i = 0; /* only one center */ CHANGE_DOFS_2(this); } } #undef CHANGE_DOFS_1 #undef CHANGE_DOFS_2 /****************************************************************************/ /* opp_vertex checks whether the face with vertices dof[0],..,dof[DIM-1] is */ /* part of mel's boundary. returns the opposite vertex if true, -1 else */ /****************************************************************************/ int Element::oppVertex(FixVec pdof) const { int nv = 0; int ov = 0; int vertices = mesh->getGeo(VERTEX); int dim = mesh->getDim(); for (int i = 0; i < vertices; i++) { if (nv < i - 1) return(-1); for (int j = 0; j < dim; j++) { if (dof[i] == pdof[j]) { /****************************************************************************/ /* i is a common vertex */ /****************************************************************************/ ov += i; nv++; break; } } } if (nv != mesh->getDim()) return(-1); /****************************************************************************/ /* the opposite vertex is 3(6) - (sum of indices of common vertices) in */ /* 2d(3d) */ /****************************************************************************/ switch(mesh->getDim()) { case 1: return ov; break; case 2: return 3 - ov; break; case 3: return 6 - ov; break; default: ERROR_EXIT("invalid dim\n"); return 0; } } void Element::eraseNewCoord() { if (newCoord != NULL) { delete newCoord; newCoord = NULL; } } void Element::serialize(std::ostream &out) { // write children if (child[0]) { out << child[0]->getTypeName() << "\n"; child[0]->serialize(out); child[1]->serialize(out); } else { out << "NULL\n"; } // write dofs int dim = mesh->getDim(); int nodes = mesh->getNumberOfNodes(); int j = 0; out.write(reinterpret_cast(&nodes), sizeof(int)); for (int pos = 0; pos <= dim; pos++) { GeoIndex position = INDEX_OF_DIM(pos, dim); int ndof = 0; for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) { ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position); } if (ndof > 0) { for (int i = 0; i < mesh->getGeo(position); i++) { if (dof[j] != NULL) { if (Mesh::serializedDOFs[dof[j][0]] == NULL) { Mesh::serializedDOFs[dof[j][0]] = dof[j]; out.write(reinterpret_cast(&ndof), sizeof(int)); out.write(reinterpret_cast(dof[j]), ndof * sizeof(DegreeOfFreedom)); } else { int minusOne = -1; out.write(reinterpret_cast(&minusOne), sizeof(int)); out.write(reinterpret_cast(&(dof[j][0])), sizeof(DegreeOfFreedom)); } } else { int zero = 0; out.write(reinterpret_cast(&zero), sizeof(int)); } j++; } } } // write index out.write(reinterpret_cast(&index), sizeof(int)); // write mark out.write(reinterpret_cast(&mark), sizeof(signed char)); // write newCoord if (newCoord) { out << "WorldVector\n"; newCoord->serialize(out); } else { out << "NULL\n"; } // write element data if (elementData) { out << elementData->getTypeName() << "\n"; elementData->serialize(out); } else { out << "NULL\n"; } } void Element::deserialize(std::istream &in) { std::string typeName; // read children in >> typeName; in.get(); if (typeName != "NULL") { if (typeName == "Line") { child[0] = new Line(NULL); child[1] = new Line(NULL); } if (typeName == "Triangle") { child[0] = new Triangle(NULL); child[1] = new Triangle(NULL); } if (typeName == "Tetrahedron") { child[0] = new Tetrahedron(NULL); child[1] = new Tetrahedron(NULL); } child[0]->deserialize(in); child[1]->deserialize(in); } else { child[0] = child[1] = NULL; } // read dofs int nodes; in.read(reinterpret_cast(&nodes), sizeof(int)); dof = new DegreeOfFreedom*[nodes]; for (int i = 0; i < nodes; i++) { int dofs; in.read(reinterpret_cast(&dofs), sizeof(int)); if (dofs) { if (dofs != -1) { dof[i] = new DegreeOfFreedom[dofs]; in.read(reinterpret_cast(dof[i]), dofs * sizeof(DegreeOfFreedom)); if (Mesh::serializedDOFs[dof[i][0]] != NULL) { DegreeOfFreedom *dofPtr = Mesh::serializedDOFs[dof[i][0]]; delete [] dof[i]; dof[i] = dofPtr; } else { Mesh::serializedDOFs[dof[i][0]] = dof[i]; } } else { DegreeOfFreedom index; in.read(reinterpret_cast(&index), sizeof(DegreeOfFreedom)); dof[i] = Mesh::serializedDOFs[index]; } } else { dof[i] = NULL; } } // read index in.read(reinterpret_cast(&index), sizeof(int)); // read mark in.read(reinterpret_cast(&mark), sizeof(signed char)); // read newCoord in >> typeName; in.get(); if (typeName != "NULL") { if (typeName == "WorldVector") { newCoord = new WorldVector; newCoord->deserialize(in); } else { ERROR_EXIT("unexpected type name\n"); } } else { newCoord = NULL; } // read element data in >> typeName; in.get(); if (typeName != "NULL") { elementData = CreatorMap::getCreator(typeName)->create(); if (elementData) { elementData->deserialize(in); } else { ERROR_EXIT("unexpected type name\n"); } } else { elementData = NULL; } } int Element::calcMemoryUsage() { int result = 0; result += sizeof(Element); result += mesh->getNumberOfNodes() * sizeof(DegreeOfFreedom*); if (child[0]) { result += child[0]->calcMemoryUsage() + child[1]->calcMemoryUsage(); } return result; } }