Select Git revision
linearSystemPETSc.h
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Jonathan Lambrechts authoredJonathan Lambrechts authored
meshRefine.cpp 33.01 KiB
// 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 "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. Contributed by Jose
// Paulo Moitinho de Almeida, 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();
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();
}
}
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]);
}
}