Select Git revision
meshRefine.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
meshRefine.cpp 29.37 KiB
// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
//
// 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);
class MVertexLessThanParam{
public:
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(int i=el->getNumPrimaryVertices(); i<el->getNumVertices(); i++)
setBLData(el->getVertex(i));
return true;
}
static void Subdivide(GEdge *ge)
{
std::vector<MLine*> lines2;
for(unsigned int 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(unsigned int 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)
{
if(!splitIntoQuads && !splitIntoHexas){
std::vector<MTriangle*> triangles2;
for(unsigned int 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(unsigned int 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(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
if(t->getNumVertices() == 6){
SPoint2 pt;
SPoint3 ptx; t->pnt(0.5, 0.5, 0, ptx);
bool reparamOK = true;
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 (reparamOK){
GPoint gp = gf->point(pt);
newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, pt[0], pt[1]);
}
else {
newv = new MVertex(ptx.x(), ptx.y(), ptx.z(), gf);
}
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(unsigned int 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(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
MTetrahedron *t = gr->tetrahedra[i];
// FIXME: we should choose the template to maximize the quality
if(t->getNumVertices() == 10){
tetrahedra2.push_back
(new MTetrahedron(t->getVertex(0), t->getVertex(4), t->getVertex(7), t->getVertex(6)));
tetrahedra2.push_back
(new MTetrahedron(t->getVertex(1), t->getVertex(4), t->getVertex(5), t->getVertex(9)));
tetrahedra2.push_back
(new MTetrahedron(t->getVertex(2), t->getVertex(5), t->getVertex(6), t->getVertex(8)));
tetrahedra2.push_back
(new MTetrahedron(t->getVertex(3), t->getVertex(7), t->getVertex(9), t->getVertex(8)));
tetrahedra2.push_back
(new MTetrahedron(t->getVertex(5), t->getVertex(8), t->getVertex(7), t->getVertex(9)));
tetrahedra2.push_back
(new MTetrahedron(t->getVertex(5), t->getVertex(7), t->getVertex(4), t->getVertex(9))); tetrahedra2.push_back
(new MTetrahedron(t->getVertex(7), t->getVertex(8), t->getVertex(5), t->getVertex(6)));
tetrahedra2.push_back
(new MTetrahedron(t->getVertex(4), t->getVertex(7), t->getVertex(5), t->getVertex(6)));
setBLData(t);
}
delete t;
}
gr->tetrahedra = tetrahedra2;
}
// Split hexes into other hexes.
std::vector<MHexahedron*> hexahedra2;
for(unsigned int 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(unsigned int 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(unsigned int 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 DIVISION OF A PYRAMID INTO 88 HEXES !!!!!
// --------------------------------------------------
std::vector<MHexahedron*> dwarfs88;
for(unsigned int 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();
// --------------------------------------------------
// ---- THANKS TO TRISTAN CARRIER BAUDOUIN ----------
// --------------------------------------------------
}
gr->hexahedra = hexahedra2;
std::vector<MPrism*> prisms2;
for(unsigned int 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(unsigned int 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(unsigned int 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)
{
// splitIntoQuads = true;
// splitIntoHexas = true;
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);
//if(splitIntoQuads) recombineIntoQuads(*it, true, true);
}
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);
}
///------ 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)
{
unsigned int i;
int index1,index2,index3,index4;
int index5,index6,index7,index8;
SPoint3 point;
std::vector<MVertex*> v;
dwarfs88.resize(88);
v.resize(105);
for (int i=0;i<105;i++)v[i] = 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);
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(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]);
}
}
for(i=0;i<88;i++){
index1 = schneiders_connect(0,i);
index2 = schneiders_connect(1,i);
index3 = schneiders_connect(2,i);
index4 = schneiders_connect(3,i);
index5 = schneiders_connect(4,i);
index6 = schneiders_connect(5,i);
index7 = schneiders_connect(6,i);
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]));
}
}