diff --git a/CREDITS.txt b/CREDITS.txt
index 29ae579657e3785b3179e5a935b13df053efb63f..134c835f7fbf22e5baea415547119b20c88e57fc 100644
--- a/CREDITS.txt
+++ b/CREDITS.txt
@@ -32,9 +32,9 @@ Pierre-Alexandre Beaufort (HXT/reparam), Zhidong Han (LSDYNA export), Ismail
Badia (hierarchical basis functions), Jeremy Theler (X3D export), Thomas
Toulorge (high order mesh optimizer, new CGNS IO), Max Orok (binary PLY), Marek
Wojciechowski (PyPi packaging), Maxence Reberol (automatic transfinite, quad
-meshing tools), Michael Ermakov (Gambit export, Fortran API), Alex Krasner (X3D
-export). See comments in the sources for more information. If we forgot to list
-your contributions please send us an email!
+meshing tools), Michael Ermakov (Gambit export, Fortran API, transfinite3), Alex
+Krasner (X3D export). See comments in the sources for more information. If we
+forgot to list your contributions please send us an email!
Thanks to the following folks who have contributed by providing fresh ideas on
theoretical or programming topics, who have sent patches, requests for changes
diff --git a/src/geo/GFace.cpp b/src/geo/GFace.cpp
index a1205cdbf3e99897487180ec3dbad1ba20ef9271..2f28db174ec09c77460aa54086fbae9a0f99b948 100644
--- a/src/geo/GFace.cpp
+++ b/src/geo/GFace.cpp
@@ -237,6 +237,7 @@ void GFace::resetMeshAttributes()
meshAttributes.meshSizeFactor = 1.;
meshAttributes.algorithm = 0;
meshAttributes.meshSizeFromBoundary = -1;
+ meshAttributes.transfinite3 = false;
}
SBoundingBox3d GFace::bounds(bool fast)
diff --git a/src/geo/GFace.h b/src/geo/GFace.h
index 3266835a5e452bb335f8001f77340f8f0cc1c833..16fe0a5829aa5d3b616fddd66427cb0c48374d15 100644
--- a/src/geo/GFace.h
+++ b/src/geo/GFace.h
@@ -319,6 +319,8 @@ public:
int algorithm;
// do we force calculation of mesh size from boundary (if >= 0)
int meshSizeFromBoundary;
+ // do we use a specific transfinite algorithm for 3-sided surfaces
+ bool transfinite3;
} meshAttributes;
int getMeshingAlgo() const;
diff --git a/src/geo/GModelIO_OCC.cpp b/src/geo/GModelIO_OCC.cpp
index 67a339fcaa4900dae464ab16fe1ec20e804f8323..fab492d6dc88db0bbf7c8de3b440a35d39241efe 100644
--- a/src/geo/GModelIO_OCC.cpp
+++ b/src/geo/GModelIO_OCC.cpp
@@ -2054,6 +2054,11 @@ static bool makeTrimmedSurface(const Handle(Geom_Surface) &surf,
if(makeEdgeOnSurface(edge, surf, wire3D, edgeOnSurf))
w.Add(edgeOnSurf);
}
+ w.Build();
+ if(!w.IsDone()) {
+ Msg::Error("Could not create wire");
+ return false;
+ }
TopoDS_Wire wire = w.Wire();
wiresProj.push_back(wire);
}
diff --git a/src/geo/discreteFace.cpp b/src/geo/discreteFace.cpp
index 9b0f389ffabce4980dad6a5763001c344c758e24..c067ee676c2759568279e355ab56deb139a3bf0f 100644
--- a/src/geo/discreteFace.cpp
+++ b/src/geo/discreteFace.cpp
@@ -931,4 +931,5 @@ void discreteFace::resetMeshAttributes()
meshAttributes.reverseMesh = _s->ReverseMesh;
meshAttributes.algorithm = _s->MeshAlgorithm;
meshAttributes.meshSizeFromBoundary = _s->MeshSizeFromBoundary;
+ meshAttributes.transfinite3 = false;
}
diff --git a/src/geo/gmshFace.cpp b/src/geo/gmshFace.cpp
index 9024d0059abdb6a5658f5de00600f29864db25ec..5134df2c62587f19aff1910d74695c608ce46549 100644
--- a/src/geo/gmshFace.cpp
+++ b/src/geo/gmshFace.cpp
@@ -127,6 +127,7 @@ void gmshFace::resetMeshAttributes()
meshAttributes.reverseMesh = _s->ReverseMesh;
meshAttributes.algorithm = _s->MeshAlgorithm;
meshAttributes.meshSizeFromBoundary = _s->MeshSizeFromBoundary;
+ meshAttributes.transfinite3 = false;
}
Range<double> gmshFace::parBounds(int i) const { return Range<double>(0, 1); }
diff --git a/src/mesh/meshGFaceTransfinite.cpp b/src/mesh/meshGFaceTransfinite.cpp
index 9534703c3ba263a8190a76084da71bf07f2fa643..947aed8cf2b555452143e5c599a94ff9b390cf5e 100644
--- a/src/mesh/meshGFaceTransfinite.cpp
+++ b/src/mesh/meshGFaceTransfinite.cpp
@@ -2,6 +2,9 @@
//
// See the LICENSE.txt file in the Gmsh root directory for license information.
// Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
+//
+// Contributor(s):
+// Michael Ermakov (ermakov@ipmnet.ru)
#include <map>
#include <queue>
@@ -152,11 +155,11 @@ static void computeEdgeLoops(const GFace *gf,
double maxDistParam(const std::vector<double> &U, const std::vector<double> &V)
{
if(U.size() < 2 || (U.size() != V.size())) return 1e22;
- double d = sqrt(std::pow(U.back() - U.front(), 2) +
- std::pow(V.back() - V.front(), 2));
+ double d =
+ sqrt(std::pow(U.back() - U.front(), 2) + std::pow(V.back() - V.front(), 2));
for(std::size_t i = 1; i < U.size(); i++)
- d = std::max(d, sqrt(std::pow(U[i] - U[i - 1], 2) +
- std::pow(V[i] - V[i - 1], 2)));
+ d = std::max(
+ d, sqrt(std::pow(U[i] - U[i - 1], 2) + std::pow(V[i] - V[i - 1], 2)));
return d;
}
@@ -260,6 +263,8 @@ int MeshTransfiniteSurface(GFace *gf)
V = V2;
}
+ bool transfinite3 = false;
+
int N1 = N[0], N2 = N[1], N3 = N[2], N4 = N[3];
int L = N2 - N1, H = N3 - N2;
if(corners.size() == 4) {
@@ -273,14 +278,26 @@ int MeshTransfiniteSurface(GFace *gf)
}
else {
int Lb = m_vertices.size() - N3;
- if(Lb != L) {
- Msg::Error("Surface %d cannot be meshed using the transfinite algo "
- "(divisions %d != %d)",
- gf->tag(), L, Lb);
- return 0;
+#ifdef TFTria
+ if(Lb == L && H == L) {
+ transfinite3 = true;
+ Msg::Info("Using specific algorithm for 3-sided surface %d", gf->tag());
+ }
+ else {
+#endif
+ if(Lb != L) {
+ Msg::Error("Surface %d cannot be meshed using the transfinite algo "
+ "(divisions %d != %d)",
+ gf->tag(), L, Lb);
+ return 0;
+ }
+#ifdef TFTria
}
+#endif
}
+ gf->meshAttributes.transfinite3 = transfinite3;
+
/*
2L+H +------------+ L+H
| |
@@ -304,7 +321,8 @@ int MeshTransfiniteSurface(GFace *gf)
MVertex *v2 = m_vertices[i + 1];
double d1 = v1->distance(v2);
if(symmetric) {
- MVertex *v3 = m_vertices[(corners.size() == 3 && !i) ? 0 : (2 * L + H - i)];
+ MVertex *v3 =
+ m_vertices[(corners.size() == 3 && !i) ? 0 : (2 * L + H - i)];
MVertex *v4 = m_vertices[2 * L + H - i - 1];
double d2 = v3->distance(v4);
L_i += 0.5 * (d1 + d2);
@@ -329,7 +347,7 @@ int MeshTransfiniteSurface(GFace *gf)
double d2 = v3->distance(v4);
L_j += 0.5 * (d1 + d2);
}
- else{
+ else {
L_j += d1;
}
lengths_j.push_back(L_j);
@@ -397,6 +415,14 @@ int MeshTransfiniteSurface(GFace *gf)
}
}
else {
+ std::vector<double> u2, v2;
+ if(transfinite3) {
+ u2.reserve(H + 1);
+ for(int j = 0; j <= H; j++) u2.push_back(U[N2 + j]);
+ v2.reserve(H + 1);
+ for(int j = 0; j <= H; j++) v2.push_back(V[N2 + j]);
+ }
+
for(int i = 1; i < L; i++) {
double u = lengths_i[i] / L_i;
for(int j = 1; j < H; j++) {
@@ -406,8 +432,29 @@ int MeshTransfiniteSurface(GFace *gf)
int iP3 = ((N3 + N2) - i) % m_vertices.size();
double Up, Vp;
if(gf->geomType() != GEntity::RuledSurface) {
- Up = TRAN_TRI(U[iP1], U[iP2], U[iP3], UC1, UC2, UC3, u, v);
- Vp = TRAN_TRI(V[iP1], V[iP2], V[iP3], VC1, VC2, VC3, u, v);
+ if(!transfinite3) {
+ Up = TRAN_TRI(U[iP1], U[iP2], U[iP3], UC1, UC2, UC3, u, v);
+ Vp = TRAN_TRI(V[iP1], V[iP2], V[iP3], VC1, VC2, VC3, u, v);
+ }
+ else {
+ if(j >= i) {
+ tab[i][j] = tab[i][H];
+ continue;
+ }
+
+ double t = double(j) / double(i);
+ double v = t;
+ int k;
+ for(k = 1; k <= H; k++)
+ if(t < lengths_j[k] / L_j && t > lengths_j[k - 1] / L_j) break;
+
+ double a =
+ (t * L_j - lengths_j[k - 1]) / (lengths_j[k] - lengths_j[k - 1]);
+ double UP2 = u2[k - 1] + a * (u2[k] - u2[k - 1]);
+ double VP2 = v2[k - 1] + a * (v2[k] - v2[k - 1]);
+ Up = TRAN_TRI(U[iP1], UP2, U[iP3], UC1, UC2, UC3, u, v);
+ Vp = TRAN_TRI(V[iP1], VP2, V[iP3], VC1, VC2, VC3, u, v);
+ }
}
else {
// FIXME: to get nice meshes we would need to make the u,v
@@ -442,7 +489,7 @@ int MeshTransfiniteSurface(GFace *gf)
else if(gf->meshAttributes.transfiniteSmoothing > 0)
numSmooth = gf->meshAttributes.transfiniteSmoothing;
- if(corners.size() == 4 && numSmooth) {
+ if(corners.size() == 4 && numSmooth && !transfinite3) {
std::vector<std::vector<double> > u(L + 1), v(L + 1);
for(int i = 0; i <= L; i++) {
u[i].resize(H + 1);
@@ -543,31 +590,54 @@ int MeshTransfiniteSurface(GFace *gf)
}
}
else {
- for(int j = 0; j < H; j++) {
- MVertex *v1 = tab[0][0];
- MVertex *v2 = tab[1][j];
- MVertex *v3 = tab[1][j + 1];
- gf->triangles.push_back(new MTriangle(v1, v2, v3));
- }
- for(int i = 1; i < L; i++) {
+ if(!transfinite3) {
for(int j = 0; j < H; j++) {
+ MVertex *v1 = tab[0][0];
+ MVertex *v2 = tab[1][j];
+ MVertex *v3 = tab[1][j + 1];
+ gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ }
+ for(int i = 1; i < L; i++) {
+ for(int j = 0; j < H; j++) {
+ MVertex *v1 = tab[i][j];
+ MVertex *v2 = tab[i + 1][j];
+ MVertex *v3 = tab[i + 1][j + 1];
+ MVertex *v4 = tab[i][j + 1];
+ if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine)
+ gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
+ else if(gf->meshAttributes.transfiniteArrangement == 1 ||
+ (gf->meshAttributes.transfiniteArrangement == 2 &&
+ ((i % 2 == 0 && j % 2 == 1) ||
+ (i % 2 == 1 && j % 2 == 0))) ||
+ (gf->meshAttributes.transfiniteArrangement == -2 &&
+ ((i % 2 == 0 && j % 2 == 0) ||
+ (i % 2 == 1 && j % 2 == 1)))) {
+ gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ gf->triangles.push_back(new MTriangle(v3, v4, v1));
+ }
+ else {
+ gf->triangles.push_back(new MTriangle(v1, v2, v4));
+ gf->triangles.push_back(new MTriangle(v4, v2, v3));
+ }
+ }
+ }
+ }
+ else {
+ for(int i = 0; i < L; i++) {
+ int j = i;
MVertex *v1 = tab[i][j];
MVertex *v2 = tab[i + 1][j];
MVertex *v3 = tab[i + 1][j + 1];
- MVertex *v4 = tab[i][j + 1];
- if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine)
- gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
- else if(gf->meshAttributes.transfiniteArrangement == 1 ||
- (gf->meshAttributes.transfiniteArrangement == 2 &&
- ((i % 2 == 0 && j % 2 == 1) || (i % 2 == 1 && j % 2 == 0))) ||
- (gf->meshAttributes.transfiniteArrangement == -2 &&
- ((i % 2 == 0 && j % 2 == 0) || (i % 2 == 1 && j % 2 == 1)))) {
+ gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ }
+ for(int i = 1; i < L; i++) {
+ for(int j = 0; j < i; j++) {
+ MVertex *v1 = tab[i][j];
+ MVertex *v2 = tab[i + 1][j];
+ MVertex *v3 = tab[i + 1][j + 1];
+ MVertex *v4 = tab[i][j + 1];
+ gf->triangles.push_back(new MTriangle(v1, v3, v4));
gf->triangles.push_back(new MTriangle(v1, v2, v3));
- gf->triangles.push_back(new MTriangle(v3, v4, v1));
- }
- else {
- gf->triangles.push_back(new MTriangle(v1, v2, v4));
- gf->triangles.push_back(new MTriangle(v4, v2, v3));
}
}
}
diff --git a/src/mesh/meshGRegionTransfinite.cpp b/src/mesh/meshGRegionTransfinite.cpp
index 6a6be9343057f2e696a4fcc2f427db8be2199da3..dc6fe11f7285609ac1c6eb949029a56d572e4644 100644
--- a/src/mesh/meshGRegionTransfinite.cpp
+++ b/src/mesh/meshGRegionTransfinite.cpp
@@ -5,6 +5,7 @@
//
// Contributor(s):
// Trevor S. Strickler
+// Michael Ermakov (ermakov@ipmnet.ru)
//
#include <map>
@@ -68,6 +69,14 @@
tab[i + 1][j + 1][k + 1], tab[i][j + 1][k + 1], \
tab[i + 1][j][k + 1])
+#define CREATE_PRISM_3 \
+ new MPrism(tab[i][j][k], tab[i][j + 1][k], tab[i + 1][j + 1][k], \
+ tab[i][j][k + 1], tab[i][j + 1][k + 1], tab[i + 1][j + 1][k + 1])
+
+#define CREATE_PRISM_4 \
+ new MPrism(tab[i][j][k], tab[i + 1][j][k], tab[i + 1][j + 1][k], \
+ tab[i][j][k + 1], tab[i + 1][j][k + 1], tab[i + 1][j + 1][k + 1])
+
#define CREATE_SIM_1 \
new MTetrahedron(tab[i][j][k], tab[i + 1][j][k], tab[i][j + 1][k], \
tab[i][j][k + 1])
@@ -92,6 +101,30 @@
new MTetrahedron(tab[i + 1][j][k + 1], tab[i][j + 1][k + 1], \
tab[i + 1][j + 1][k + 1], tab[i + 1][j + 1][k])
+#define CREATE_SIM_7 \
+ new MTetrahedron(tab[i][j][k], tab[i][j + 1][k], tab[i][j][k + 1], \
+ tab[i + 1][j + 1][k])
+
+#define CREATE_SIM_8 \
+ new MTetrahedron(tab[i][j + 1][k], tab[i][j + 1][k + 1], tab[i][j][k + 1], \
+ tab[i + 1][j + 1][k])
+
+#define CREATE_SIM_9 \
+ new MTetrahedron(tab[i][j][k + 1], tab[i][j + 1][k + 1], \
+ tab[i + 1][j + 1][k + 1], tab[i + 1][j + 1][k])
+
+#define CREATE_SIM_10 \
+ new MTetrahedron(tab[i][j][k], tab[i + 1][j][k], tab[i + 1][j + 1][k], \
+ tab[i][j][k + 1])
+
+#define CREATE_SIM_11 \
+ new MTetrahedron(tab[i + 1][j][k], tab[i + 1][j + 1][k], tab[i][j][k + 1], \
+ tab[i + 1][j][k + 1])
+
+#define CREATE_SIM_12 \
+ new MTetrahedron(tab[i][j][k + 1], tab[i + 1][j][k + 1], \
+ tab[i + 1][j + 1][k], tab[i + 1][j + 1][k + 1])
+
static double transfiniteHex(double f1, double f2, double f3, double f4,
double f5, double f6, double c1, double c2,
double c3, double c4, double c5, double c6,
@@ -587,82 +620,28 @@ int MeshTransfiniteVolume(GRegion *gr)
}
}
else if(faces.size() == 5) {
- for(int j = 0; j < N_j - 1; j++) {
- for(int k = 0; k < N_k - 1; k++) {
-#if defined(HAVE_QUADTRI)
- if(gr->meshAttributes.QuadTri) {
- // create vertex array
- std::vector<MVertex *> verts;
- verts.resize(6);
- verts[0] = tab[0][j][k];
- verts[1] = tab[1][j][k];
- verts[2] = tab[1][j + 1][k];
- verts[3] = tab[0][j][k + 1];
- verts[4] = tab[1][j][k + 1];
- verts[5] = tab[1][j + 1][k + 1];
- if((!orientedFaces[0].recombined() && j == 0) ||
- (!orientedFaces[2].recombined() && j == N_j - 2) ||
- (!orientedFaces[4].recombined() && k == 0) ||
- (!orientedFaces[5].recombined() && k == N_k - 2)) {
- // make subdivided element
- meshTransfElemWithInternalVertex(gr, verts, &boundary_diags);
- }
- else
- gr->prisms.push_back(new MPrism(verts[0], verts[1], verts[2],
- verts[3], verts[4], verts[5]));
- // continue, skipping the rest which is for non-divided elements
- continue;
- }
-#endif
- if((orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
- orientedFaces[2].recombined() && orientedFaces[4].recombined() &&
- orientedFaces[5].recombined()) ||
- (orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
- orientedFaces[2].recombined() && !orientedFaces[4].recombined() &&
- !orientedFaces[5].recombined())) {
- gr->prisms.push_back(new MPrism(
- tab[0][j][k], tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1],
- tab[1][j][k + 1], tab[1][j + 1][k + 1]));
- }
- else if(!orientedFaces[0].recombined() &&
- !orientedFaces[1].recombined() &&
- !orientedFaces[2].recombined() &&
- !orientedFaces[4].recombined() &&
- !orientedFaces[5].recombined()) {
- gr->tetrahedra.push_back(new MTetrahedron(
- tab[0][j][k], tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1]));
- gr->tetrahedra.push_back(
- new MTetrahedron(tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1],
- tab[1][j][k + 1]));
- gr->tetrahedra.push_back(
- new MTetrahedron(tab[0][j][k + 1], tab[1][j + 1][k + 1],
- tab[1][j][k + 1], tab[1][j + 1][k]));
- }
- else {
- Msg::Error("Wrong surface recombination in transfinite volume %d",
- gr->tag());
- return 0;
- }
+ bool standard_algo = true;
+ for(int i = 0; i < 5; i++) {
+ if(orientedFaces[i].getSurface() &&
+ orientedFaces[i].getSurface()->meshAttributes.transfinite3) {
+ standard_algo = false;
+ break;
}
}
- for(int i = 1; i < N_i - 1; i++) {
+ if(standard_algo) {
for(int j = 0; j < N_j - 1; j++) {
for(int k = 0; k < N_k - 1; k++) {
#if defined(HAVE_QUADTRI)
if(gr->meshAttributes.QuadTri) {
- // create vertex array
std::vector<MVertex *> verts;
- verts.resize(8);
- verts[0] = tab[i][j][k];
- verts[1] = tab[i + 1][j][k];
- verts[2] = tab[i + 1][j + 1][k];
- verts[3] = tab[i][j + 1][k];
- verts[4] = tab[i][j][k + 1];
- verts[5] = tab[i + 1][j][k + 1];
- verts[6] = tab[i + 1][j + 1][k + 1];
- verts[7] = tab[i][j + 1][k + 1];
- if((!orientedFaces[1].recombined() && i == N_i - 2) ||
- (!orientedFaces[0].recombined() && j == 0) ||
+ verts.resize(6);
+ verts[0] = tab[0][j][k];
+ verts[1] = tab[1][j][k];
+ verts[2] = tab[1][j + 1][k];
+ verts[3] = tab[0][j][k + 1];
+ verts[4] = tab[1][j][k + 1];
+ verts[5] = tab[1][j + 1][k + 1];
+ if((!orientedFaces[0].recombined() && j == 0) ||
(!orientedFaces[2].recombined() && j == N_j - 2) ||
(!orientedFaces[4].recombined() && k == 0) ||
(!orientedFaces[5].recombined() && k == N_k - 2)) {
@@ -670,37 +649,125 @@ int MeshTransfiniteVolume(GRegion *gr)
meshTransfElemWithInternalVertex(gr, verts, &boundary_diags);
}
else
- gr->hexahedra.push_back(
- new MHexahedron(verts[0], verts[1], verts[2], verts[3],
- verts[4], verts[5], verts[6], verts[7]));
+ gr->prisms.push_back(new MPrism(verts[0], verts[1], verts[2],
+ verts[3], verts[4], verts[5]));
// continue, skipping the rest which is for non-divided elements
continue;
}
#endif
- if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
- orientedFaces[2].recombined() && orientedFaces[4].recombined() &&
- orientedFaces[5].recombined()) {
- gr->hexahedra.push_back(CREATE_HEX);
+ if((orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
+ orientedFaces[2].recombined() && orientedFaces[4].recombined() &&
+ orientedFaces[5].recombined()) ||
+ (orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
+ orientedFaces[2].recombined() && !orientedFaces[4].recombined() &&
+ !orientedFaces[5].recombined())) {
+ gr->prisms.push_back(new MPrism(
+ tab[0][j][k], tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1],
+ tab[1][j][k + 1], tab[1][j + 1][k + 1]));
}
- else if(orientedFaces[0].recombined() &&
- orientedFaces[1].recombined() &&
- orientedFaces[2].recombined() &&
+ else if(!orientedFaces[0].recombined() &&
+ !orientedFaces[1].recombined() &&
+ !orientedFaces[2].recombined() &&
!orientedFaces[4].recombined() &&
!orientedFaces[5].recombined()) {
- gr->prisms.push_back(CREATE_PRISM_1);
- gr->prisms.push_back(CREATE_PRISM_2);
+ gr->tetrahedra.push_back(new MTetrahedron(
+ tab[0][j][k], tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1]));
+ gr->tetrahedra.push_back(
+ new MTetrahedron(tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1],
+ tab[1][j][k + 1]));
+ gr->tetrahedra.push_back(
+ new MTetrahedron(tab[0][j][k + 1], tab[1][j + 1][k + 1],
+ tab[1][j][k + 1], tab[1][j + 1][k]));
+ }
+ else {
+ Msg::Error("Wrong surface recombination in transfinite volume %d",
+ gr->tag());
+ return 0;
+ }
+ }
+ }
+ for(int i = 1; i < N_i - 1; i++) {
+ for(int j = 0; j < N_j - 1; j++) {
+ for(int k = 0; k < N_k - 1; k++) {
+#if defined(HAVE_QUADTRI)
+ if(gr->meshAttributes.QuadTri) {
+ // create vertex array
+ std::vector<MVertex *> verts;
+ verts.resize(8);
+ verts[0] = tab[i][j][k];
+ verts[1] = tab[i + 1][j][k];
+ verts[2] = tab[i + 1][j + 1][k];
+ verts[3] = tab[i][j + 1][k];
+ verts[4] = tab[i][j][k + 1];
+ verts[5] = tab[i + 1][j][k + 1];
+ verts[6] = tab[i + 1][j + 1][k + 1];
+ verts[7] = tab[i][j + 1][k + 1];
+ if((!orientedFaces[1].recombined() && i == N_i - 2) ||
+ (!orientedFaces[0].recombined() && j == 0) ||
+ (!orientedFaces[2].recombined() && j == N_j - 2) ||
+ (!orientedFaces[4].recombined() && k == 0) ||
+ (!orientedFaces[5].recombined() && k == N_k - 2)) {
+ // make subdivided element
+ meshTransfElemWithInternalVertex(gr, verts, &boundary_diags);
+ }
+ else
+ gr->hexahedra.push_back(
+ new MHexahedron(verts[0], verts[1], verts[2], verts[3],
+ verts[4], verts[5], verts[6], verts[7]));
+ // continue, skipping the rest which is for non-divided elements
+ continue;
+ }
+#endif
+ if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
+ orientedFaces[2].recombined() && orientedFaces[4].recombined() &&
+ orientedFaces[5].recombined()) {
+ gr->hexahedra.push_back(CREATE_HEX);
+ }
+ else if(orientedFaces[0].recombined() &&
+ orientedFaces[1].recombined() &&
+ orientedFaces[2].recombined() &&
+ !orientedFaces[4].recombined() &&
+ !orientedFaces[5].recombined()) {
+ gr->prisms.push_back(CREATE_PRISM_1);
+ gr->prisms.push_back(CREATE_PRISM_2);
+ }
+ else if(!orientedFaces[0].recombined() &&
+ !orientedFaces[1].recombined() &&
+ !orientedFaces[2].recombined() &&
+ !orientedFaces[4].recombined() &&
+ !orientedFaces[5].recombined()) {
+ gr->tetrahedra.push_back(CREATE_SIM_1);
+ gr->tetrahedra.push_back(CREATE_SIM_2);
+ gr->tetrahedra.push_back(CREATE_SIM_3);
+ gr->tetrahedra.push_back(CREATE_SIM_4);
+ gr->tetrahedra.push_back(CREATE_SIM_5);
+ gr->tetrahedra.push_back(CREATE_SIM_6);
+ }
+ else {
+ Msg::Error("Wrong surface recombination in transfinite volume %d",
+ gr->tag());
+ return 0;
+ }
+ }
+ }
+ }
+ }
+ else { // if surfaces meshed with specific algo for 3-sided surfaces
+ for(int i = 0; i < N_i - 1; i++) {
+ int j = i;
+ for(int k = 0; k < N_k - 1; k++) {
+ if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
+ orientedFaces[2].recombined()) {
+ gr->prisms.push_back(CREATE_PRISM_4);
}
else if(!orientedFaces[0].recombined() &&
!orientedFaces[1].recombined() &&
!orientedFaces[2].recombined() &&
!orientedFaces[4].recombined() &&
!orientedFaces[5].recombined()) {
- gr->tetrahedra.push_back(CREATE_SIM_1);
- gr->tetrahedra.push_back(CREATE_SIM_2);
- gr->tetrahedra.push_back(CREATE_SIM_3);
- gr->tetrahedra.push_back(CREATE_SIM_4);
- gr->tetrahedra.push_back(CREATE_SIM_5);
- gr->tetrahedra.push_back(CREATE_SIM_6);
+ gr->tetrahedra.push_back(CREATE_SIM_10);
+ gr->tetrahedra.push_back(CREATE_SIM_11);
+ gr->tetrahedra.push_back(CREATE_SIM_12);
}
else {
Msg::Error("Wrong surface recombination in transfinite volume %d",
@@ -709,6 +776,34 @@ int MeshTransfiniteVolume(GRegion *gr)
}
}
}
+ for(int i = 1; i < N_i - 1; i++) {
+ for(int j = 0; j < i; j++) {
+ for(int k = 0; k < N_k - 1; k++) {
+ if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
+ orientedFaces[2].recombined()) {
+ gr->prisms.push_back(CREATE_PRISM_3);
+ gr->prisms.push_back(CREATE_PRISM_4);
+ }
+ else if(!orientedFaces[0].recombined() &&
+ !orientedFaces[1].recombined() &&
+ !orientedFaces[2].recombined() &&
+ !orientedFaces[4].recombined() &&
+ !orientedFaces[5].recombined()) {
+ gr->tetrahedra.push_back(CREATE_SIM_7);
+ gr->tetrahedra.push_back(CREATE_SIM_8);
+ gr->tetrahedra.push_back(CREATE_SIM_9);
+ gr->tetrahedra.push_back(CREATE_SIM_10);
+ gr->tetrahedra.push_back(CREATE_SIM_11);
+ gr->tetrahedra.push_back(CREATE_SIM_12);
+ }
+ else {
+ Msg::Error("Wrong surface recombination in transfinite volume %d",
+ gr->tag());
+ return 0;
+ }
+ }
+ }
+ }
}
}