Hello Trying to understand the quality of the meshes obtained when using the "Refine by splitting" option I looked at the code in Mesh/meshRefine.cpp
After some digging I came to realise why the comment at line 198 is there, >> // FIXME: we should choose the template to maximize the quality at first thinking that the "right choice" would imply some geometry calculations. Then I found > Bey, Jürgen. "Simplicial grid refinement: on Freudenthal's algorithm and the > optimal number of congruence classes." Numerische Mathematik 85.1 (2000): > 1-29. and updated the template using the sequence in Algorithm RedRefinement3D therein. At first the result was not correct, but a simple switch of indices made it work. Starting with a mesh where all elements are equal we obtain refined meshes with elements of 3 shapes and for a non-regular mesh the properties of the elements stabilise. The code is attached, I hope it is useful. (I know that this could be done in a more efficient way with git. Sorry I will learn how to do that next week..). Regards ZP
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // issues on https://gitlab.onelab.info/gmsh/gmsh/issues. // // Contributor(s): // Brian Helenbrook // #include "HighOrder.h" #include "MLine.h" #include "MTriangle.h" #include "MQuadrangle.h" #include "MTetrahedron.h" #include "MHexahedron.h" #include "MPrism.h" #include "MPyramid.h" #include "GmshMessage.h" #include "OS.h" #include "Context.h" #include "meshGFaceOptimize.h" void subdivide_pyramid(MElement *element, GRegion *gr, faceContainer &faceVertices, std::vector<MHexahedron *> &dwarfs88); struct MVertexLessThanParam { bool operator()(const MVertex *v1, const MVertex *v2) const { double u1 = 0., u2 = 1.; v1->getParameter(0, u1); v2->getParameter(0, u2); return u1 < u2; } }; // Set BM data on vertex static void setBLData(MVertex *v) { switch(v->onWhat()->dim()) { case 1: { MEdgeVertex *ve = dynamic_cast<MEdgeVertex *>(v); if(ve) ve->bl_data = new MVertexBoundaryLayerData(); break; } case 2: { MFaceVertex *vf = dynamic_cast<MFaceVertex *>(v); if(vf) vf->bl_data = new MVertexBoundaryLayerData(); break; } } } // If all low-order nodes in are marked as BL, then mark high-order nodes as BL // (only works in 2D) static bool setBLData(MElement *el) { // Check whether all low-order nodes are marked as BL nodes (only works in 2D) for(int i = 0; i < el->getNumPrimaryVertices(); i++) { MVertex *v = el->getVertex(i); bool isBL = false; switch(v->onWhat()->dim()) { case 0: isBL = true; break; case 1: { MEdgeVertex *ve = dynamic_cast<MEdgeVertex *>(v); if(ve && ve->bl_data) isBL = true; break; } case 2: { MFaceVertex *vf = dynamic_cast<MFaceVertex *>(v); if(vf && vf->bl_data) isBL = true; break; } } if(!isBL) return false; } // Mark high-order nodes as BL nodes (only works in 2D) for(std::size_t i = el->getNumPrimaryVertices(); i < el->getNumVertices(); i++) setBLData(el->getVertex(i)); return true; } static void Subdivide(GEdge *ge) { std::vector<MLine *> lines2; for(std::size_t i = 0; i < ge->lines.size(); i++) { MLine *l = ge->lines[i]; if(l->getNumVertices() == 3) { lines2.push_back(new MLine(l->getVertex(0), l->getVertex(2))); lines2.push_back(new MLine(l->getVertex(2), l->getVertex(1))); setBLData(l); } delete l; } ge->lines = lines2; // 2nd order meshing destroyed the ordering of the vertices on the edge std::sort(ge->mesh_vertices.begin(), ge->mesh_vertices.end(), MVertexLessThanParam()); for(std::size_t i = 0; i < ge->mesh_vertices.size(); i++) ge->mesh_vertices[i]->setPolynomialOrder(1); ge->deleteVertexArrays(); } static void Subdivide(GFace *gf, bool splitIntoQuads, bool splitIntoHexas, faceContainer &faceVertices, bool linear) { if(!splitIntoQuads && !splitIntoHexas) { std::vector<MTriangle *> triangles2; for(std::size_t i = 0; i < gf->triangles.size(); i++) { MTriangle *t = gf->triangles[i]; if(t->getNumVertices() == 6) { triangles2.push_back( new MTriangle(t->getVertex(0), t->getVertex(3), t->getVertex(5))); triangles2.push_back( new MTriangle(t->getVertex(3), t->getVertex(4), t->getVertex(5))); triangles2.push_back( new MTriangle(t->getVertex(3), t->getVertex(1), t->getVertex(4))); triangles2.push_back( new MTriangle(t->getVertex(5), t->getVertex(4), t->getVertex(2))); setBLData(t); } delete t; } gf->triangles = triangles2; } std::vector<MQuadrangle *> quadrangles2; for(std::size_t i = 0; i < gf->quadrangles.size(); i++) { MQuadrangle *q = gf->quadrangles[i]; if(q->getNumVertices() == 9) { quadrangles2.push_back(new MQuadrangle(q->getVertex(0), q->getVertex(4), q->getVertex(8), q->getVertex(7))); quadrangles2.push_back(new MQuadrangle(q->getVertex(4), q->getVertex(1), q->getVertex(5), q->getVertex(8))); quadrangles2.push_back(new MQuadrangle(q->getVertex(8), q->getVertex(5), q->getVertex(2), q->getVertex(6))); quadrangles2.push_back(new MQuadrangle(q->getVertex(7), q->getVertex(8), q->getVertex(6), q->getVertex(3))); setBLData(q); } delete q; } if(splitIntoQuads || splitIntoHexas) { for(std::size_t i = 0; i < gf->triangles.size(); i++) { MTriangle *t = gf->triangles[i]; if(t->getNumVertices() == 6) { SPoint2 pt; SPoint3 ptx; t->pnt(1. / 3., 1. / 3., 0., ptx); bool reparamOK = true; if(!linear){ for(int k = 0; k < 6; k++) { SPoint2 temp; reparamOK &= reparamMeshVertexOnFace(t->getVertex(k), gf, temp); pt[0] += temp[0] / 6.; pt[1] += temp[1] / 6.; } } MVertex *newv; if(linear || !reparamOK) { newv = new MVertex(ptx.x(), ptx.y(), ptx.z(), gf); } else { GPoint gp = gf->point(pt); newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, pt[0], pt[1]); } gf->mesh_vertices.push_back(newv); if(splitIntoHexas) faceVertices[t->getFace(0)].push_back(newv); quadrangles2.push_back(new MQuadrangle(t->getVertex(0), t->getVertex(3), newv, t->getVertex(5))); quadrangles2.push_back(new MQuadrangle(t->getVertex(3), t->getVertex(1), t->getVertex(4), newv)); quadrangles2.push_back(new MQuadrangle( t->getVertex(5), newv, t->getVertex(4), t->getVertex(2))); if(setBLData(t)) setBLData(newv); delete t; } } gf->triangles.clear(); } gf->quadrangles = quadrangles2; for(std::size_t i = 0; i < gf->mesh_vertices.size(); i++) gf->mesh_vertices[i]->setPolynomialOrder(1); gf->deleteVertexArrays(); } static void Subdivide(GRegion *gr, bool splitIntoHexas, faceContainer &faceVertices) { if(!splitIntoHexas) { // Split tets into other tets std::vector<MTetrahedron *> tetrahedra2; for(std::size_t i = 0; i < gr->tetrahedra.size(); i++) { MTetrahedron *t = gr->tetrahedra[i]; // Use a template that maximizes the quality, which is a modification of // Algorithm RedRefinement3D in: // Bey, Jürgen. "Simplicial grid refinement: on Freudenthal's algorithm // and the optimal number of congruence classes." Numerische Mathematik 85.1 (2000): 1-29.) // JPM April 2019 if(t->getNumVertices() == 10) { tetrahedra2.push_back(new MTetrahedron( t->getVertex(0), t->getVertex(4), t->getVertex(6), t->getVertex(7))); tetrahedra2.push_back(new MTetrahedron( t->getVertex(4), t->getVertex(1), t->getVertex(5), t->getVertex(9))); tetrahedra2.push_back(new MTetrahedron( t->getVertex(6), t->getVertex(5), t->getVertex(2), t->getVertex(8))); tetrahedra2.push_back(new MTetrahedron( t->getVertex(7), t->getVertex(9), t->getVertex(8), t->getVertex(3))); tetrahedra2.push_back(new MTetrahedron( t->getVertex(4), t->getVertex(6), t->getVertex(7), t->getVertex(9))); tetrahedra2.push_back(new MTetrahedron( t->getVertex(4), t->getVertex(9), t->getVertex(5), t->getVertex(6))); tetrahedra2.push_back(new MTetrahedron( t->getVertex(6), t->getVertex(7), t->getVertex(9), t->getVertex(8))); tetrahedra2.push_back(new MTetrahedron( t->getVertex(6), t->getVertex(8), t->getVertex(9), t->getVertex(5))); setBLData(t); } delete t; } gr->tetrahedra = tetrahedra2; } // Split hexes into other hexes. std::vector<MHexahedron *> hexahedra2; for(std::size_t i = 0; i < gr->hexahedra.size(); i++) { MHexahedron *h = gr->hexahedra[i]; if(h->getNumVertices() == 27) { hexahedra2.push_back(new MHexahedron(h->getVertex(0), h->getVertex(8), h->getVertex(20), h->getVertex(9), h->getVertex(10), h->getVertex(21), h->getVertex(26), h->getVertex(22))); hexahedra2.push_back(new MHexahedron( h->getVertex(10), h->getVertex(21), h->getVertex(26), h->getVertex(22), h->getVertex(4), h->getVertex(16), h->getVertex(25), h->getVertex(17))); hexahedra2.push_back(new MHexahedron(h->getVertex(8), h->getVertex(1), h->getVertex(11), h->getVertex(20), h->getVertex(21), h->getVertex(12), h->getVertex(23), h->getVertex(26))); hexahedra2.push_back(new MHexahedron( h->getVertex(21), h->getVertex(12), h->getVertex(23), h->getVertex(26), h->getVertex(16), h->getVertex(5), h->getVertex(18), h->getVertex(25))); hexahedra2.push_back(new MHexahedron(h->getVertex(9), h->getVertex(20), h->getVertex(13), h->getVertex(3), h->getVertex(22), h->getVertex(26), h->getVertex(24), h->getVertex(15))); hexahedra2.push_back(new MHexahedron( h->getVertex(22), h->getVertex(26), h->getVertex(24), h->getVertex(15), h->getVertex(17), h->getVertex(25), h->getVertex(19), h->getVertex(7))); hexahedra2.push_back(new MHexahedron(h->getVertex(20), h->getVertex(11), h->getVertex(2), h->getVertex(13), h->getVertex(26), h->getVertex(23), h->getVertex(14), h->getVertex(24))); hexahedra2.push_back(new MHexahedron( h->getVertex(26), h->getVertex(23), h->getVertex(14), h->getVertex(24), h->getVertex(25), h->getVertex(18), h->getVertex(6), h->getVertex(19))); setBLData(h); } delete h; } // Split tets into other hexes. if(splitIntoHexas) { for(std::size_t i = 0; i < gr->tetrahedra.size(); i++) { MTetrahedron *t = gr->tetrahedra[i]; if(t->getNumVertices() == 10) { std::vector<MVertex *> newv; for(int j = 0; j < t->getNumFaces(); j++) { MFace face = t->getFace(j); faceContainer::iterator fIter = faceVertices.find(face); if(fIter != faceVertices.end()) { newv.push_back(fIter->second[0]); } else { SPoint3 pc = face.barycenter(); newv.push_back(new MVertex(pc.x(), pc.y(), pc.z(), gr)); faceVertices[face].push_back(newv.back()); gr->mesh_vertices.push_back(newv.back()); } } SPoint3 pc = t->barycenter(); newv.push_back(new MVertex(pc.x(), pc.y(), pc.z(), gr)); gr->mesh_vertices.push_back(newv.back()); hexahedra2.push_back(new MHexahedron( t->getVertex(0), t->getVertex(4), newv[0], t->getVertex(6), t->getVertex(7), newv[1], newv[4], newv[2])); hexahedra2.push_back( new MHexahedron(t->getVertex(4), t->getVertex(1), t->getVertex(5), newv[0], newv[1], t->getVertex(9), newv[3], newv[4])); hexahedra2.push_back(new MHexahedron( t->getVertex(6), newv[0], t->getVertex(5), t->getVertex(2), newv[2], newv[4], newv[3], t->getVertex(8))); hexahedra2.push_back(new MHexahedron( t->getVertex(3), t->getVertex(9), newv[1], t->getVertex(7), t->getVertex(8), newv[3], newv[4], newv[2])); if(setBLData(t)) { setBLData(newv[0]); setBLData(newv[1]); setBLData(newv[2]); setBLData(newv[3]); setBLData(newv[4]); } delete t; } } gr->tetrahedra.clear(); for(std::size_t i = 0; i < gr->prisms.size(); i++) { MPrism *p = gr->prisms[i]; if(p->getNumVertices() == 18) { std::vector<MVertex *> newv; for(int j = 0; j < 2; j++) { MFace face = p->getFace(j); faceContainer::iterator fIter = faceVertices.find(face); if(fIter != faceVertices.end()) { newv.push_back(fIter->second[0]); } else { SPoint3 pc = face.barycenter(); newv.push_back(new MVertex(pc.x(), pc.y(), pc.z(), gr)); faceVertices[face].push_back(newv.back()); gr->mesh_vertices.push_back(newv.back()); } } SPoint3 pc = p->barycenter(); newv.push_back(new MVertex(pc.x(), pc.y(), pc.z(), gr)); gr->mesh_vertices.push_back(newv.back()); hexahedra2.push_back(new MHexahedron( p->getVertex(0), p->getVertex(6), newv[0], p->getVertex(7), p->getVertex(8), p->getVertex(15), newv[2], p->getVertex(16))); hexahedra2.push_back(new MHexahedron( p->getVertex(1), p->getVertex(9), newv[0], p->getVertex(6), p->getVertex(10), p->getVertex(17), newv[2], p->getVertex(15))); hexahedra2.push_back(new MHexahedron( p->getVertex(2), p->getVertex(7), newv[0], p->getVertex(9), p->getVertex(11), p->getVertex(16), newv[2], p->getVertex(17))); hexahedra2.push_back(new MHexahedron( p->getVertex(8), p->getVertex(15), newv[2], p->getVertex(16), p->getVertex(3), p->getVertex(12), newv[1], p->getVertex(13))); hexahedra2.push_back(new MHexahedron( p->getVertex(10), p->getVertex(17), newv[2], p->getVertex(15), p->getVertex(4), p->getVertex(14), newv[1], p->getVertex(12))); hexahedra2.push_back(new MHexahedron( p->getVertex(11), p->getVertex(16), newv[2], p->getVertex(17), p->getVertex(5), p->getVertex(13), newv[1], p->getVertex(14))); if(setBLData(p)) { setBLData(newv[0]); setBLData(newv[1]); setBLData(newv[2]); } } } gr->prisms.clear(); // Yamakawa subdivision of a pyramid into 88 hexes (thanks to Tristan // Carrier Baudouin!) std::vector<MHexahedron *> dwarfs88; for(std::size_t i = 0; i < gr->pyramids.size(); i++) { MPyramid *p = gr->pyramids[i]; if(p->getNumVertices() == 14) { subdivide_pyramid(p, gr, faceVertices, dwarfs88); for(int j = 0; j < 88; j++) hexahedra2.push_back(dwarfs88[j]); } } gr->pyramids.clear(); } gr->hexahedra = hexahedra2; std::vector<MPrism *> prisms2; for(std::size_t i = 0; i < gr->prisms.size(); i++) { MPrism *p = gr->prisms[i]; if(p->getNumVertices() == 18) { prisms2.push_back(new MPrism(p->getVertex(0), p->getVertex(6), p->getVertex(7), p->getVertex(8), p->getVertex(15), p->getVertex(16))); prisms2.push_back(new MPrism(p->getVertex(8), p->getVertex(15), p->getVertex(16), p->getVertex(3), p->getVertex(12), p->getVertex(13))); prisms2.push_back(new MPrism(p->getVertex(6), p->getVertex(1), p->getVertex(9), p->getVertex(15), p->getVertex(10), p->getVertex(17))); prisms2.push_back(new MPrism(p->getVertex(15), p->getVertex(10), p->getVertex(17), p->getVertex(12), p->getVertex(4), p->getVertex(14))); prisms2.push_back(new MPrism(p->getVertex(7), p->getVertex(9), p->getVertex(2), p->getVertex(16), p->getVertex(17), p->getVertex(11))); prisms2.push_back(new MPrism(p->getVertex(16), p->getVertex(17), p->getVertex(11), p->getVertex(13), p->getVertex(14), p->getVertex(5))); prisms2.push_back(new MPrism(p->getVertex(9), p->getVertex(7), p->getVertex(6), p->getVertex(17), p->getVertex(16), p->getVertex(15))); prisms2.push_back(new MPrism(p->getVertex(17), p->getVertex(16), p->getVertex(15), p->getVertex(14), p->getVertex(13), p->getVertex(12))); setBLData(p); } delete p; } gr->prisms = prisms2; std::vector<MPyramid *> pyramids2; for(std::size_t i = 0; i < gr->pyramids.size(); i++) { if(splitIntoHexas) { Msg::Error("Full hexahedron subdivision is not implemented for pyramids"); return; } MPyramid *p = gr->pyramids[i]; if(p->getNumVertices() == 14) { // Base pyramids2.push_back(new MPyramid(p->getVertex(0), p->getVertex(5), p->getVertex(13), p->getVertex(6), p->getVertex(7))); pyramids2.push_back(new MPyramid(p->getVertex(5), p->getVertex(1), p->getVertex(8), p->getVertex(13), p->getVertex(9))); pyramids2.push_back(new MPyramid(p->getVertex(13), p->getVertex(8), p->getVertex(2), p->getVertex(10), p->getVertex(11))); pyramids2.push_back(new MPyramid(p->getVertex(6), p->getVertex(13), p->getVertex(10), p->getVertex(3), p->getVertex(12))); // Split remaining into tets // Top gr->tetrahedra.push_back((new MTetrahedron( p->getVertex(7), p->getVertex(9), p->getVertex(12), p->getVertex(4)))); gr->tetrahedra.push_back((new MTetrahedron( p->getVertex(9), p->getVertex(11), p->getVertex(12), p->getVertex(4)))); // Upside down one gr->tetrahedra.push_back( (new MTetrahedron(p->getVertex(9), p->getVertex(12), p->getVertex(11), p->getVertex(13)))); gr->tetrahedra.push_back((new MTetrahedron( p->getVertex(7), p->getVertex(12), p->getVertex(9), p->getVertex(13)))); // Four tets around bottom perimeter gr->tetrahedra.push_back((new MTetrahedron( p->getVertex(7), p->getVertex(9), p->getVertex(5), p->getVertex(13)))); gr->tetrahedra.push_back((new MTetrahedron( p->getVertex(9), p->getVertex(11), p->getVertex(8), p->getVertex(13)))); gr->tetrahedra.push_back( (new MTetrahedron(p->getVertex(12), p->getVertex(10), p->getVertex(11), p->getVertex(13)))); gr->tetrahedra.push_back((new MTetrahedron( p->getVertex(7), p->getVertex(6), p->getVertex(12), p->getVertex(13)))); setBLData(p); } delete p; } gr->pyramids = pyramids2; for(std::size_t i = 0; i < gr->mesh_vertices.size(); i++) gr->mesh_vertices[i]->setPolynomialOrder(1); gr->deleteVertexArrays(); } void RefineMesh(GModel *m, bool linear, bool splitIntoQuads, bool splitIntoHexas) { Msg::StatusBar(true, "Refining mesh..."); double t1 = Cpu(); // Create 2nd order mesh (using "2nd order complete" elements) to // generate vertex positions SetOrderN(m, 2, linear, false); // only used when splitting tets into hexes faceContainer faceVertices; // Subdivide the second order elements to create the refined linear // mesh for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) Subdivide(*it); for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) Subdivide(*it, splitIntoQuads, splitIntoHexas, faceVertices, linear); for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) Subdivide(*it, splitIntoHexas, faceVertices); // Check all 3D elements for negative volume and reverse if needed m->setAllVolumesPositive(); CTX::instance()->mesh.changed = ENT_ALL; double t2 = Cpu(); Msg::StatusBar(true, "Done refining mesh (%g s)", t2 - t1); } void BarycentricRefineMesh(GModel *m) { Msg::StatusBar(true, "Barycentrically refining mesh..."); double t1 = Cpu(); m->destroyMeshCaches(); // Only update triangles in 2D, only update tets in 3D if(m->getNumRegions() == 0) { for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) { GFace *gf = *it; std::size_t numt = gf->triangles.size(); if(!numt) continue; std::vector<MTriangle *> triangles2(3 * numt); for(std::size_t i = 0; i < numt; i++) { MTriangle *t = gf->triangles[i]; SPoint3 bary = t->barycenter(); // FIXME: create an MFaceVertex (with correct parametric coordinates)? MVertex *v = new MVertex(bary.x(), bary.y(), bary.z(), gf); triangles2[3 * i] = new MTriangle(t->getVertex(0), t->getVertex(1), v); triangles2[3 * i + 1] = new MTriangle(t->getVertex(1), t->getVertex(2), v); triangles2[3 * i + 2] = new MTriangle(t->getVertex(2), t->getVertex(0), v); delete t; gf->mesh_vertices.push_back(v); } gf->triangles = triangles2; gf->deleteVertexArrays(); } } else { for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) { GRegion *gr = *it; std::size_t numt = gr->tetrahedra.size(); if(!numt) continue; std::vector<MTetrahedron *> tetrahedra2(4 * numt); for(std::size_t i = 0; i < numt; i++) { MTetrahedron *t = gr->tetrahedra[i]; SPoint3 bary = t->barycenter(); // FIXME: create an MFaceVertex (with correct parametric coordinates)? MVertex *v = new MVertex(bary.x(), bary.y(), bary.z(), gr); tetrahedra2[4 * i] = new MTetrahedron(t->getVertex(0), t->getVertex(1), t->getVertex(2), v); tetrahedra2[4 * i + 1] = new MTetrahedron( t->getVertex(1), t->getVertex(2), t->getVertex(3), v); tetrahedra2[4 * i + 2] = new MTetrahedron( t->getVertex(2), t->getVertex(3), t->getVertex(0), v); tetrahedra2[4 * i + 3] = new MTetrahedron( t->getVertex(3), t->getVertex(0), t->getVertex(1), v); delete t; gr->mesh_vertices.push_back(v); } gr->tetrahedra = tetrahedra2; gr->deleteVertexArrays(); } } CTX::instance()->mesh.changed = ENT_ALL; double t2 = Cpu(); Msg::StatusBar(true, "Done barycentrically refining mesh (%g s)", t2 - t1); } // Tristan Carrier Baudouin's contribution on Full Hex Meshing static double schneiders_x(int i) { double coord[105] = { 0.500000, 0.666667, 0.500000, 1.000000, 1.000000, 1.000000, 0.289057, 0.324970, 0.276710, 0.337200, 0.364878, 0.325197, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -0.289057, -0.324970, -0.276710, -0.337200, -0.364878, -0.325197, -0.500000, -0.666667, -0.500000, -1.000000, -1.000000, -1.000000, 0.084599, 0.263953, 0.442960, 0.310954, 0.000000, 0.000000, 0.118244, 0.212082, 0.244049, 0.213940, 0.040495, 0.110306, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -0.118244, -0.212082, -0.244049, -0.213940, -0.040495, -0.110306, -0.084599, -0.263953, -0.442960, -0.310954, 0.650493, 0.454537, 0.000000, 0.000000, 0.320508, 0.264129, 0.063695, 0.092212, 0.000000, 0.000000, 0.000000, 0.000000, -0.320508, -0.264129, -0.063695, -0.092212, -0.650493, -0.454537, 0.619616, 0.000000, 0.277170, 0.124682, 0.000000, 0.000000, -0.277170, -0.124682, -0.619616, 0.128101, 0.000000, 0.176104, 0.084236, 0.000000, 0.000000, -0.176104, -0.084236, -0.128101, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000}; return coord[i]; } static double schneiders_y(int i) { double coord[105] = { 0.707107, 0.471405, 0.707107, 0.000000, 0.000000, 0.000000, 0.474601, 0.421483, 0.463847, 0.367090, 0.344843, 0.361229, 0.429392, 0.410339, 0.435001, 0.407674, 0.401208, 0.404164, 0.474601, 0.421483, 0.463847, 0.367090, 0.344843, 0.361229, 0.707107, 0.471405, 0.707107, 0.000000, 0.000000, 0.000000, 0.536392, 0.697790, 0.569124, 0.742881, 0.558045, 0.946690, 0.497378, 0.532513, 0.500133, 0.541479, 0.530503, 0.666252, 0.463274, 0.465454, 0.430972, 0.451135, 0.515274, 0.589713, 0.497378, 0.532513, 0.500133, 0.541479, 0.530503, 0.666252, 0.536392, 0.697790, 0.569124, 0.742881, 0.080018, 0.104977, 0.216381, 0.200920, 0.326416, 0.339933, 0.280915, 0.305725, 0.396502, 0.394423, 0.310617, 0.337499, 0.326416, 0.339933, 0.280915, 0.305725, 0.080018, 0.104977, 0.081443, 0.204690, 0.354750, 0.334964, 0.403611, 0.367496, 0.354750, 0.334964, 0.081443, 0.501199, 0.538575, 0.447454, 0.486224, 0.431723, 0.470065, 0.447454, 0.486224, 0.501199, 0.488701, 0.471405, 0.017336, 0.000000, 0.452197, 0.471405, 0.057682, 0.000000, 1.414213, 0.015731, 0.000000}; return coord[i]; } static double schneiders_z(int i) { double coord[105] = { 0.500000, 0.000000, -0.500000, 1.000000, 0.000000, -1.000000, 0.051666, -0.058015, -0.148140, 0.071987, -0.057896, -0.181788, -0.016801, -0.054195, -0.104114, -0.015276, -0.054392, -0.110605, 0.051666, -0.058015, -0.148140, 0.071987, -0.057896, -0.181788, 0.500000, 0.000000, -0.500000, 1.000000, 0.000000, -1.000000, -0.208673, -0.162945, 0.021476, 0.389516, -0.157646, 0.159885, -0.142645, -0.073557, -0.032793, 0.060339, -0.136482, 0.043449, -0.111103, -0.079664, -0.047879, -0.008734, -0.124554, 0.008560, -0.142645, -0.073557, -0.032793, 0.060339, -0.136482, 0.043449, -0.208673, -0.162945, 0.021476, 0.389516, -0.041899, -0.680880, -0.103504, -0.392255, -0.065989, -0.212535, -0.093142, -0.227139, -0.056201, -0.124443, -0.087185, -0.182164, -0.065989, -0.212535, -0.093142, -0.227139, -0.041899, -0.680880, 0.786284, 0.443271, 0.104202, 0.144731, -0.005330, 0.073926, 0.104202, 0.144731, 0.786284, -0.364254, -0.282882, -0.189794, -0.182143, -0.127036, -0.148665, -0.189794, -0.182143, -0.364254, 0.642222, 0.666667, 0.959658, 1.000000, -0.455079, -0.666667, -0.844452, -1.000000, 0.000000, -0.009020, 0.000000}; return coord[i]; } static int schneiders_connect(int i, int j) { int n0[88] = {0, 1, 6, 7, 12, 13, 18, 19, 41, 39, 37, 41, 36, 40, 47, 45, 43, 47, 42, 46, 53, 51, 49, 53, 48, 52, 35, 57, 55, 35, 54, 34, 34, 35, 65, 63, 64, 62, 69, 67, 68, 66, 73, 71, 72, 70, 61, 75, 60, 74, 60, 61, 79, 78, 81, 80, 83, 82, 77, 84, 77, 88, 87, 90, 89, 92, 91, 86, 93, 86, 57, 24, 95, 85, 2, 99, 84, 74, 97, 104, 104, 97, 26, 35, 35, 24, 35, 24}; int n1[88] = {1, 2, 7, 8, 13, 14, 19, 20, 39, 6, 38, 37, 37, 36, 45, 12, 44, 43, 43, 42, 51, 18, 50, 49, 49, 48, 57, 24, 56, 55, 55, 54, 40, 41, 63, 11, 62, 10, 67, 17, 66, 16, 71, 23, 70, 22, 75, 29, 74, 28, 64, 65, 78, 9, 80, 15, 82, 21, 84, 27, 79, 87, 8, 89, 14, 91, 20, 93, 26, 88, 35, 57, 94, 86, 85, 98, 77, 60, 96, 103, 28, 27, 99, 55, 102, 25, 31, 102}; int n2[88] = {4, 5, 10, 11, 16, 17, 22, 23, 38, 7, 7, 36, 8, 87, 44, 13, 13, 42, 14, 89, 50, 19, 19, 48, 20, 91, 56, 25, 25, 54, 26, 93, 46, 47, 62, 10, 78, 9, 66, 16, 80, 15, 70, 22, 82, 21, 74, 28, 84, 27, 68, 69, 39, 6, 45, 12, 51, 18, 57, 24, 81, 63, 11, 67, 17, 71, 23, 75, 29, 90, 33, 94, 33, 93, 98, 93, 76, 58, 76, 58, 74, 84, 2, 26, 2, 26, 2, 0}; int n3[88] = {3, 4, 9, 10, 15, 16, 21, 22, 37, 38, 8, 40, 87, 88, 43, 44, 14, 46, 89, 90, 49, 50, 20, 52, 91, 92, 55, 56, 26, 34, 93, 86, 52, 53, 64, 62, 79, 78, 68, 66, 81, 80, 72, 70, 83, 82, 60, 74, 77, 84, 72, 73, 41, 39, 47, 45, 53, 51, 35, 57, 83, 65, 63, 69, 67, 73, 71, 61, 75, 92, 94, 95, 0, 98, 99, 26, 96, 103, 3, 4, 103, 96, 102, 102, 31, 102, 102, 95}; int n4[88] = {6, 7, 12, 13, 18, 19, 24, 25, 35, 33, 31, 35, 30, 34, 41, 39, 37, 41, 36, 40, 47, 45, 43, 47, 42, 46, 53, 51, 49, 53, 48, 52, 86, 34, 61, 59, 60, 58, 65, 63, 64, 62, 69, 67, 68, 66, 73, 71, 72, 70, 77, 60, 77, 76, 79, 78, 81, 80, 83, 82, 35, 86, 85, 88, 87, 90, 89, 92, 91, 61, 84, 27, 97, 59, 5, 101, 74, 75, 104, 101, 101, 104, 93, 34, 34, 57, 33, 57}; int n5[88] = {7, 8, 13, 14, 19, 20, 25, 26, 33, 0, 32, 31, 31, 30, 39, 6, 38, 37, 37, 36, 45, 12, 44, 43, 43, 42, 51, 18, 50, 49, 49, 48, 88, 40, 59, 5, 58, 4, 63, 11, 62, 10, 67, 17, 66, 16, 71, 23, 70, 22, 79, 64, 76, 3, 78, 9, 80, 15, 82, 21, 41, 85, 2, 87, 8, 89, 14, 91, 20, 65, 77, 84, 96, 61, 59, 100, 60, 61, 103, 100, 29, 28, 98, 54, 86, 56, 32, 35}; int n6[88] = {10, 11, 16, 17, 22, 23, 28, 29, 32, 1, 1, 30, 2, 85, 38, 7, 7, 36, 8, 87, 44, 13, 13, 42, 14, 89, 50, 19, 19, 48, 20, 91, 90, 46, 58, 4, 76, 3, 62, 10, 78, 9, 66, 16, 80, 15, 70, 22, 82, 21, 81, 68, 33, 0, 39, 6, 45, 12, 51, 18, 47, 59, 5, 63, 11, 67, 17, 71, 23, 69, 76, 96, 76, 75, 100, 75, 58, 59, 58, 59, 75, 74, 85, 93, 85, 55, 1, 33}; int n7[88] = {9, 10, 15, 16, 21, 22, 27, 28, 31, 32, 2, 34, 85, 86, 37, 38, 8, 40, 87, 88, 43, 44, 14, 46, 89, 90, 49, 50, 20, 52, 91, 92, 92, 52, 60, 58, 77, 76, 64, 62, 79, 78, 68, 66, 81, 80, 72, 70, 83, 82, 83, 72, 35, 33, 41, 39, 47, 45, 53, 51, 53, 61, 59, 65, 63, 69, 67, 73, 71, 73, 96, 97, 3, 100, 101, 29, 103, 100, 4, 5, 100, 103, 86, 86, 30, 35, 0, 94}; if(i == 0) { return n0[j]; } else if(i == 1) { return n1[j]; } else if(i == 2) { return n2[j]; } else if(i == 3) { return n3[j]; } else if(i == 4) { return n4[j]; } else if(i == 5) { return n5[j]; } else if(i == 6) { return n6[j]; } else { return n7[j]; } } void subdivide_pyramid(MElement *element, GRegion *gr, faceContainer &faceVertices, std::vector<MHexahedron *> &dwarfs88) { std::vector<MVertex *> v(105, (MVertex *)NULL); v[29] = element->getVertex(0); v[27] = element->getVertex(1); v[3] = element->getVertex(2); v[5] = element->getVertex(3); v[102] = element->getVertex(4); v[28] = element->getVertex(5); v[97] = element->getVertex(8); v[4] = element->getVertex(10); v[101] = element->getVertex(6); v[26] = element->getVertex(7); v[24] = element->getVertex(9); v[0] = element->getVertex(11); v[2] = element->getVertex(12); // the one in the middle of rect face v[104] = element->getVertex(13); SPoint3 point; faceContainer::iterator fIter; fIter = faceVertices.find(MFace(v[29], v[27], v[102])); if(fIter != faceVertices.end()) v[25] = fIter->second[0]; else { element->pnt(0.0, -0.666667, 0.471405 / 1.414213, point); v[25] = new MVertex(point.x(), point.y(), point.z(), gr); gr->addMeshVertex(v[25]); faceVertices[MFace(v[29], v[27], v[102])].push_back(v[25]); } fIter = faceVertices.find(MFace(v[27], v[3], v[102])); if(fIter != faceVertices.end()) v[95] = fIter->second[0]; else { element->pnt(0.666667, 0.0, 0.471405 / 1.414213, point); v[95] = new MVertex(point.x(), point.y(), point.z(), gr); gr->addMeshVertex(v[95]); faceVertices[MFace(v[27], v[3], v[102])].push_back(v[95]); } fIter = faceVertices.find(MFace(v[3], v[5], v[102])); if(fIter != faceVertices.end()) v[1] = fIter->second[0]; else { element->pnt(0.0, 0.666667, 0.471405 / 1.414213, point); v[1] = new MVertex(point.x(), point.y(), point.z(), gr); gr->addMeshVertex(v[1]); faceVertices[MFace(v[3], v[5], v[102])].push_back(v[1]); } fIter = faceVertices.find(MFace(v[5], v[29], v[102])); if(fIter != faceVertices.end()) v[99] = fIter->second[0]; else { element->pnt(-0.666667, 0.0, 0.471405 / 1.414213, point); v[99] = new MVertex(point.x(), point.y(), point.z(), gr); gr->addMeshVertex(v[99]); faceVertices[MFace(v[5], v[29], v[102])].push_back(v[99]); } for(int i = 0; i < 105; i++) { if(!v[i]) { element->pnt(schneiders_z(i), schneiders_x(i), schneiders_y(i) / 1.414213, point); v[i] = new MVertex(point.x(), point.y(), point.z(), gr); gr->addMeshVertex(v[i]); } } dwarfs88.resize(88); for(int i = 0; i < 88; i++) { int const index1 = schneiders_connect(0, i); int const index2 = schneiders_connect(1, i); int const index3 = schneiders_connect(2, i); int const index4 = schneiders_connect(3, i); int const index5 = schneiders_connect(4, i); int const index6 = schneiders_connect(5, i); int const index7 = schneiders_connect(6, i); int const index8 = schneiders_connect(7, i); dwarfs88[i] = new MHexahedron(v[index1], v[index2], v[index3], v[index4], v[index5], v[index6], v[index7], v[index8]); } }
_______________________________________________ gmsh mailing list gmsh@onelab.info http://onelab.info/mailman/listinfo/gmsh