diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 886e474f09210e722c879f5fb33b342785301dcb..0fca4b945c169faad9c93e51989c03be79851ab3 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -407,6 +407,10 @@ void GetOptions(int argc, char *argv[])
i++;
opt_mesh_partition_pyr_weight(0,GMSH_SET,atoi(argv[i]));
}
+ else if (!strcmp(argv[i],"trihedron")) {
+ i++;
+ opt_mesh_partition_trih_weight(0,GMSH_SET,atoi(argv[i]));
+ }
else if (!strcmp(argv[i],"hex")) {
i++;
opt_mesh_partition_hex_weight(0,GMSH_SET,atoi(argv[i]));
diff --git a/Common/Context.cpp b/Common/Context.cpp
index 67f4277662693f53de6ce49569066ec47a065c6d..2657781376285aca2282c75b1dad3b59a9ec004e 100644
--- a/Common/Context.cpp
+++ b/Common/Context.cpp
@@ -83,7 +83,7 @@ CTX::CTX() : gamepad(0)
mesh.qualityInf = mesh.qualitySup = mesh.qualityType = 0;
mesh.radiusInf = mesh.radiusSup = 0;
mesh.lines = mesh.triangles = mesh.tetrahedra = mesh.quadrangles = 0;
- mesh.prisms = mesh.pyramids = mesh.hexahedra = 0;
+ mesh.prisms = mesh.pyramids = mesh.hexahedra = mesh.trihedra = 0;
mesh.volumesEdges = mesh.volumesFaces = mesh.surfacesEdges = mesh.surfacesFaces = 0;
mesh.volumesFaces = mesh.surfacesEdges = mesh.surfacesFaces = 0;
mesh.hoOptimize = mesh.smoothNormals = mesh.reverseAllNormals = 0;
@@ -93,7 +93,7 @@ CTX::CTX() : gamepad(0)
mesh.ignorePartBound = 0;
mesh.saveTri = 0;
color.mesh.tangents = color.mesh.tetrahedron = color.mesh.triangle = 0;
- color.mesh.prism = color.mesh.pyramid = color.mesh.hexahedron = 0;
+ color.mesh.prism = color.mesh.pyramid = color.mesh.hexahedron = color.mesh.trihedron = 0;
color.mesh.tangents = color.mesh.line = color.mesh.quadrangle = 0;
for (int i = 0; i < 20; i++)
color.mesh.carousel[i] = 0;
diff --git a/Common/Context.h b/Common/Context.h
index 0913deddcc1c0c5cd7624ae3fc4fc12e38c6a362..c8dad77afbebd6b9776b8992cd2dae111d7ba924 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -20,7 +20,7 @@ class GamePad;
struct contextMeshOptions {
int draw, changed, light, lightTwoSide, lightLines, pointType;
- int points, lines, triangles, quadrangles, tetrahedra, hexahedra, prisms, pyramids;
+ int points, lines, triangles, quadrangles, tetrahedra, hexahedra, prisms, pyramids, trihedra;
int surfacesEdges, surfacesFaces, volumesEdges, volumesFaces, numSubEdges;
int pointsNum, linesNum, surfacesNum, volumesNum, qualityType, labelType;
int optimize, optimizeNetgen, optimizeLloyd, smoothCrossField, refineSteps;
@@ -293,7 +293,7 @@ class CTX {
} geom;
struct{
unsigned int vertex, vertexSup, line, triangle, quadrangle;
- unsigned int tetrahedron, hexahedron, prism, pyramid;
+ unsigned int tetrahedron, hexahedron, prism, pyramid, trihedron;
unsigned int carousel[20];
unsigned int tangents, normals;
} mesh;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 5a203ea8b323a6722544367a24100723d9a15ee9..10c36cba74e01445cba47b66e5813dabb7756f0c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1096,6 +1096,8 @@ StringXNumber MeshOptions_Number[] = {
"Weight of prismatic element (wedge) for METIS load balancing" },
{ F|O, "PartitionPyramidWeight" , opt_mesh_partition_pyr_weight , 1 ,
"Weight of pyramidal element for METIS load balancing" },
+ { F|O, "PartitionTrihedronWeight" , opt_mesh_partition_trih_weight , 0 ,
+ "Weight of trihedron element for METIS load balancing" },
{ F|O, "PartitionQuadWeight" , opt_mesh_partition_qua_weight , 1 ,
"Weight of quadrangle for METIS load balancing" },
{ F|O, "PartitionTetWeight" , opt_mesh_partition_tet_weight , 1 ,
@@ -1113,6 +1115,8 @@ StringXNumber MeshOptions_Number[] = {
"Number of prisms in the current mesh (read-only)" },
{ F, "NbPyramids" , opt_mesh_nb_pyramids , 0. ,
"Number of pyramids in the current mesh (read-only)" },
+ { F, "NbTrihedra" , opt_mesh_nb_trihedra , 0. ,
+ "Number of trihedra in the current mesh (read-only)" },
{ F, "NbQuadrangles" , opt_mesh_nb_quadrangles , 0. ,
"Number of quadrangles in the current mesh (read-only)" },
{ F, "NbTetrahedra" , opt_mesh_nb_tetrahedra , 0. ,
@@ -1144,6 +1148,8 @@ StringXNumber MeshOptions_Number[] = {
"Display mesh prisms?" },
{ F|O, "Pyramids" , opt_mesh_pyramids , 1. ,
"Display mesh pyramids?" },
+ { F|O, "Trihedra" , opt_mesh_trihedra , 1. ,
+ "Display mesh trihedra?" },
{ F|O, "Quadrangles" , opt_mesh_quadrangles , 1. ,
"Display mesh quadrangles?" },
@@ -1445,6 +1451,8 @@ StringXNumber ViewOptions_Number[] = {
"Display post-processing prisms?" },
{ F|O, "DrawPyramids" , opt_view_draw_pyramids , 1. ,
"Display post-processing pyramids?" },
+ { F|O, "DrawTrihedra" , opt_view_draw_trihedra , 1. ,
+ "Display post-processing trihedra?" },
{ F|O, "DrawQuadrangles" , opt_view_draw_quadrangles , 1. ,
"Display post-processing quadrangles?" },
{ F|O, "DrawScalars" , opt_view_draw_scalars , 1. ,
@@ -1798,6 +1806,7 @@ StringXColor GeometryOptions_Color[] = {
#define COLQ {130, 120, 225, 255}
#define COLP {232, 210, 23, 255}
#define COLY {217, 113, 38, 255}
+#define COLR {20, 255, 0, 255}
#define COL0 {255, 120, 0, 255}
#define COL2 {255, 160, 0, 255}
@@ -1848,6 +1857,9 @@ StringXColor MeshOptions_Color[] = {
{ F|O, "Pyramids" , opt_mesh_color_pyramid ,
COLY, COLY, COLW, COLY,
"Mesh pyramid color (if Mesh.ColorCarousel=0)" },
+ { F|O, "Trihedra" , opt_mesh_color_trihedron ,
+ COLR, COLR, COLW, COLR,
+ "Mesh trihedron color (if Mesh.ColorCarousel=0)" },
{ F|O, "Tangents" , opt_mesh_color_tangents ,
{255, 255, 0, 255}, {255, 255, 0, 255}, {0, 0, 0, 255}, {255, 255, 0, 255},
"Tangent mesh vector color" },
@@ -1916,6 +1928,7 @@ StringXColor ViewOptions_Color[] = {
{ F|O, "Hexahedra" , opt_view_color_hexahedra , ELECOL, "Hexahedron color" },
{ F|O, "Prisms" , opt_view_color_prisms , ELECOL, "Prism color" },
{ F|O, "Pyramids" , opt_view_color_pyramids , ELECOL, "Pyramid color" },
+ { F|O, "Trihedra" , opt_view_color_trihedra , ELECOL, "Trihedron color" },
{ F|O, "Tangents" , opt_view_color_tangents ,
{255, 255, 0, 255}, {255, 255, 0, 255}, {0, 0, 0, 255}, {255, 255, 0, 255},
"Tangent vector color" },
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index a36ee969052d2b394d6e1cae439d98d8a966ca93..d68ec33f65abc395e944a4cd0780eaa6e24b05c2 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -67,6 +67,7 @@
#define TYPE_POLYH 10
#define TYPE_XFEM 11
#define TYPE_MINI 12
+#define TYPE_TRIH 13
// Element types in .msh file format (numbers should not be changed)
#define MSH_LIN_2 1
@@ -215,8 +216,9 @@
#define MSH_TET_16 137
#define MSH_TRI_MINI 138
#define MSH_TET_MINI 139
+#define MSH_TRIH_4 140
-#define MSH_NUM_TYPE 139
+#define MSH_NUM_TYPE 140
// Geometric entities
#define ENT_NONE 0
diff --git a/Common/Options.cpp b/Common/Options.cpp
index e4ec17a867f81c1d007dcf44ebfd7406139e0148..5e0170cc31ce5eb55b7dcdcc31347259c7b3cbea 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5251,6 +5251,25 @@ double opt_mesh_pyramids(OPT_ARGS_NUM)
return CTX::instance()->mesh.pyramids;
}
+double opt_mesh_trihedra(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET) {
+ if(CTX::instance()->mesh.trihedra != val)
+ CTX::instance()->mesh.changed |= ENT_VOLUME;
+ CTX::instance()->mesh.trihedra = (int)val;
+ }
+#if defined(HAVE_FLTK)
+ if(FlGui::available() && (action & GMSH_GUI)){
+ if(CTX::instance()->mesh.trihedra)
+ ((Fl_Menu_Item*)FlGui::instance()->options->mesh.menu->menu())[6].set();
+ else
+ ((Fl_Menu_Item*)FlGui::instance()->options->mesh.menu->menu())[6].clear();
+ }
+#endif
+ return CTX::instance()->mesh.trihedra;
+}
+
+
double opt_mesh_surfaces_edges(OPT_ARGS_NUM)
{
if(action & GMSH_SET) {
@@ -5562,6 +5581,13 @@ double opt_mesh_partition_pyr_weight(OPT_ARGS_NUM)
return CTX::instance()->partitionOptions.pyrWeight;
}
+double opt_mesh_partition_trih_weight(OPT_ARGS_NUM)
+{
+ if (action & GMSH_SET)
+ CTX::instance()->partitionOptions.trihWeight = (int) val;
+ return CTX::instance()->partitionOptions.trihWeight;
+}
+
double opt_mesh_partition_qua_weight(OPT_ARGS_NUM)
{
if (action & GMSH_SET)
@@ -6129,11 +6155,18 @@ double opt_mesh_nb_pyramids(OPT_ARGS_NUM)
return s[12];
}
+double opt_mesh_nb_trihedra(OPT_ARGS_NUM)
+{
+ double s[50];
+ GetStatistics(s);
+ return s[13];
+}
+
double opt_mesh_cpu_time(OPT_ARGS_NUM)
{
double s[50];
GetStatistics(s);
- return s[13] + s[14] + s[15];
+ return s[14] + s[15] + s[16];
}
double opt_mesh_partition_partitioner(OPT_ARGS_NUM)
@@ -8134,6 +8167,28 @@ double opt_view_draw_pyramids(OPT_ARGS_NUM)
#endif
}
+double opt_view_draw_trihedra(OPT_ARGS_NUM)
+{
+#if defined(HAVE_POST)
+ GET_VIEWo(0.);
+ if(action & GMSH_SET) {
+ opt->drawTrihedra = (int)val;
+ if(view) view->setChanged(true);
+ }
+#if defined(HAVE_FLTK)
+ if(_gui_action_valid(action, num)){
+ if(opt->drawTrihedra)
+ ((Fl_Menu_Item*)FlGui::instance()->options->view.menu[1]->menu())[8].set();
+ else
+ ((Fl_Menu_Item*)FlGui::instance()->options->view.menu[1]->menu())[8].clear();
+ }
+#endif
+ return opt->drawTrihedra;
+#else
+ return 0.;
+#endif
+}
+
double opt_view_draw_scalars(OPT_ARGS_NUM)
{
#if defined(HAVE_POST)
@@ -9428,13 +9483,29 @@ unsigned int opt_mesh_color_pyramid(OPT_ARGS_COL)
return CTX::instance()->color.mesh.pyramid;
}
+unsigned int opt_mesh_color_trihedron(OPT_ARGS_COL)
+{
+ if(action & GMSH_SET) {
+ // vertex arrays need to be regenerated only when we color by
+ // element type
+ if(CTX::instance()->color.mesh.trihedron != val &&
+ CTX::instance()->mesh.colorCarousel == 0)
+ CTX::instance()->mesh.changed |= ENT_VOLUME;
+ CTX::instance()->color.mesh.trihedron = val;
+ }
+#if defined(HAVE_FLTK)
+ CCC(CTX::instance()->color.mesh.trihedron, FlGui::instance()->options->mesh.color[9]);
+#endif
+ return CTX::instance()->color.mesh.trihedron;
+}
+
unsigned int opt_mesh_color_tangents(OPT_ARGS_COL)
{
if(action & GMSH_SET) {
CTX::instance()->color.mesh.tangents = val;
}
#if defined(HAVE_FLTK)
- CCC(CTX::instance()->color.mesh.tangents, FlGui::instance()->options->mesh.color[9]);
+ CCC(CTX::instance()->color.mesh.tangents, FlGui::instance()->options->mesh.color[10]);
#endif
return CTX::instance()->color.mesh.tangents;
}
@@ -9445,7 +9516,7 @@ unsigned int opt_mesh_color_normals(OPT_ARGS_COL)
CTX::instance()->color.mesh.normals = val;
}
#if defined(HAVE_FLTK)
- CCC(CTX::instance()->color.mesh.normals, FlGui::instance()->options->mesh.color[10]);
+ CCC(CTX::instance()->color.mesh.normals, FlGui::instance()->options->mesh.color[11]);
#endif
return CTX::instance()->color.mesh.normals;
}
@@ -9461,7 +9532,7 @@ unsigned int opt_mesh_color_(int i, OPT_ARGS_COL)
CTX::instance()->color.mesh.carousel[i] = val;
}
#if defined(HAVE_FLTK)
- CCC(CTX::instance()->color.mesh.carousel[i], FlGui::instance()->options->mesh.color[11+i]);
+ CCC(CTX::instance()->color.mesh.carousel[i], FlGui::instance()->options->mesh.color[12+i]);
#endif
return CTX::instance()->color.mesh.carousel[i];
}
@@ -9639,6 +9710,25 @@ unsigned int opt_view_color_pyramids(OPT_ARGS_COL)
#endif
}
+unsigned int opt_view_color_trihedra(OPT_ARGS_COL)
+{
+#if defined(HAVE_POST)
+ GET_VIEWo(0);
+ if(action & GMSH_SET) {
+ opt->color.trihedron = val;
+ if(view) view->setChanged(true);
+ }
+#if defined(HAVE_FLTK)
+ if(_gui_action_valid(action, num)){
+ CCC(opt->color.trihedron, FlGui::instance()->options->view.color[8]);
+ }
+#endif
+ return opt->color.trihedron;
+#else
+ return 0;
+#endif
+}
+
unsigned int opt_view_color_tangents(OPT_ARGS_COL)
{
#if defined(HAVE_POST)
@@ -9649,7 +9739,7 @@ unsigned int opt_view_color_tangents(OPT_ARGS_COL)
}
#if defined(HAVE_FLTK)
if(_gui_action_valid(action, num)){
- CCC(opt->color.tangents, FlGui::instance()->options->view.color[8]);
+ CCC(opt->color.tangents, FlGui::instance()->options->view.color[9]);
}
#endif
return opt->color.tangents;
@@ -9668,7 +9758,7 @@ unsigned int opt_view_color_normals(OPT_ARGS_COL)
}
#if defined(HAVE_FLTK)
if(_gui_action_valid(action, num)){
- CCC(opt->color.normals, FlGui::instance()->options->view.color[9]);
+ CCC(opt->color.normals, FlGui::instance()->options->view.color[10]);
}
#endif
return opt->color.normals;
diff --git a/Common/Options.h b/Common/Options.h
index 1da21151e4e0935df324ce360264f850f4195214..0f3f076fad9436b46c723eb566feaa11683fc7e3 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -423,6 +423,7 @@ double opt_mesh_tetrahedra(OPT_ARGS_NUM);
double opt_mesh_hexahedra(OPT_ARGS_NUM);
double opt_mesh_prisms(OPT_ARGS_NUM);
double opt_mesh_pyramids(OPT_ARGS_NUM);
+double opt_mesh_trihedra(OPT_ARGS_NUM);
double opt_mesh_surfaces_edges(OPT_ARGS_NUM);
double opt_mesh_surfaces_faces(OPT_ARGS_NUM);
double opt_mesh_volumes_edges(OPT_ARGS_NUM);
@@ -449,6 +450,7 @@ double opt_mesh_msh_file_partitioned(OPT_ARGS_NUM);
double opt_mesh_partition_hex_weight(OPT_ARGS_NUM);
double opt_mesh_partition_pri_weight(OPT_ARGS_NUM);
double opt_mesh_partition_pyr_weight(OPT_ARGS_NUM);
+double opt_mesh_partition_trih_weight(OPT_ARGS_NUM);
double opt_mesh_partition_qua_weight(OPT_ARGS_NUM);
double opt_mesh_partition_tet_weight(OPT_ARGS_NUM);
double opt_mesh_partition_tri_weight(OPT_ARGS_NUM);
@@ -503,6 +505,7 @@ double opt_mesh_nb_tetrahedra(OPT_ARGS_NUM);
double opt_mesh_nb_hexahedra(OPT_ARGS_NUM);
double opt_mesh_nb_prisms(OPT_ARGS_NUM);
double opt_mesh_nb_pyramids(OPT_ARGS_NUM);
+double opt_mesh_nb_trihedra(OPT_ARGS_NUM);
double opt_mesh_cpu_time(OPT_ARGS_NUM);
double opt_mesh_partition_partitioner(OPT_ARGS_NUM);
double opt_mesh_partition_num(OPT_ARGS_NUM);
@@ -652,6 +655,7 @@ double opt_view_draw_tetrahedra(OPT_ARGS_NUM);
double opt_view_draw_hexahedra(OPT_ARGS_NUM);
double opt_view_draw_prisms(OPT_ARGS_NUM);
double opt_view_draw_pyramids(OPT_ARGS_NUM);
+double opt_view_draw_trihedra(OPT_ARGS_NUM);
double opt_view_draw_scalars(OPT_ARGS_NUM);
double opt_view_draw_vectors(OPT_ARGS_NUM);
double opt_view_draw_tensors(OPT_ARGS_NUM);
@@ -747,6 +751,7 @@ unsigned int opt_mesh_color_tetrahedra(OPT_ARGS_COL);
unsigned int opt_mesh_color_hexahedra(OPT_ARGS_COL);
unsigned int opt_mesh_color_prisms(OPT_ARGS_COL);
unsigned int opt_mesh_color_pyramid(OPT_ARGS_COL);
+unsigned int opt_mesh_color_trihedron(OPT_ARGS_COL);
unsigned int opt_mesh_color_tangents(OPT_ARGS_COL);
unsigned int opt_mesh_color_normals(OPT_ARGS_COL);
unsigned int opt_mesh_color_0(OPT_ARGS_COL);
@@ -777,6 +782,7 @@ unsigned int opt_view_color_tetrahedra(OPT_ARGS_COL);
unsigned int opt_view_color_hexahedra(OPT_ARGS_COL);
unsigned int opt_view_color_prisms(OPT_ARGS_COL);
unsigned int opt_view_color_pyramids(OPT_ARGS_COL);
+unsigned int opt_view_color_trihedra(OPT_ARGS_COL);
unsigned int opt_view_color_tangents(OPT_ARGS_COL);
unsigned int opt_view_color_normals(OPT_ARGS_COL);
unsigned int opt_view_color_text2d(OPT_ARGS_COL);
diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index dbf7ae25cddd3e687206ebf67d3b2568b2610d00..77df50ea95b7c64514da31968e18b78895bea766 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -496,6 +496,7 @@ static void mesh_options_ok_cb(Fl_Widget *w, void *data)
opt_mesh_hexahedra(0, GMSH_SET, o->mesh.menu->menu()[3].value() ? 1 : 0);
opt_mesh_prisms(0, GMSH_SET, o->mesh.menu->menu()[4].value() ? 1 : 0);
opt_mesh_pyramids(0, GMSH_SET, o->mesh.menu->menu()[5].value() ? 1 : 0);
+ opt_mesh_trihedra(0, GMSH_SET, o->mesh.menu->menu()[6].value() ? 1 : 0);
opt_mesh_surfaces_edges(0, GMSH_SET, o->mesh.butt[8]->value());
opt_mesh_surfaces_faces(0, GMSH_SET, o->mesh.butt[9]->value());
opt_mesh_volumes_edges(0, GMSH_SET, o->mesh.butt[10]->value());
@@ -2503,6 +2504,7 @@ optionWindow::optionWindow(int deltaFontSize)
{"Hexahedra", 0, 0, 0, FL_MENU_TOGGLE},
{"Prisms", 0, 0, 0, FL_MENU_TOGGLE},
{"Pyramids", 0, 0, 0, FL_MENU_TOGGLE},
+ {"Trihedra", 0, 0, 0, FL_MENU_TOGGLE},
{0}
};
@@ -3687,6 +3689,7 @@ void optionWindow::updateViewGroup(int index)
opt_view_draw_hexahedra(index, GMSH_GUI, 0);
opt_view_draw_prisms(index, GMSH_GUI, 0);
opt_view_draw_pyramids(index, GMSH_GUI, 0);
+ opt_view_draw_trihedra(index, GMSH_GUI, 0);
opt_view_draw_scalars(index, GMSH_GUI, 0);
opt_view_draw_vectors(index, GMSH_GUI, 0);
opt_view_draw_tensors(index, GMSH_GUI, 0);
@@ -3792,6 +3795,7 @@ void optionWindow::updateViewGroup(int index)
opt_view_color_hexahedra(index, GMSH_GUI, 0);
opt_view_color_prisms(index, GMSH_GUI, 0);
opt_view_color_pyramids(index, GMSH_GUI, 0);
+ opt_view_color_trihedra(index, GMSH_GUI, 0);
opt_view_color_tangents(index, GMSH_GUI, 0);
opt_view_color_normals(index, GMSH_GUI, 0);
opt_view_color_text2d(index, GMSH_GUI, 0);
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 42d552c343496f6874c4ae006a5c341200d8873f..439a30ab3c357f46cd0f51983815de2316d1e138 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -89,7 +89,7 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
int num = 0;
int width = 26 * FL_NORMAL_SIZE;
- int height = 5 * WB + 17 * BH;
+ int height = 5 * WB + 18 * BH;
win = new paletteWindow
(width, height, CTX::instance()->nonModalWindows ? true : false, "Statistics");
@@ -118,26 +118,27 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
value[num++] = new Fl_Output(2 * WB, 2 * WB + 7 * BH, IW, BH, "Hexahedra");
value[num++] = new Fl_Output(2 * WB, 2 * WB + 8 * BH, IW, BH, "Prisms");
value[num++] = new Fl_Output(2 * WB, 2 * WB + 9 * BH, IW, BH, "Pyramids");
+ value[num++] = new Fl_Output(2 * WB, 2 * WB + 10 * BH, IW, BH, "Trihedra");
- value[num++] = new Fl_Output(2 * WB, 2 * WB + 10 * BH, IW, BH, "Time for 1D mesh");
- value[num++] = new Fl_Output(2 * WB, 2 * WB + 11 * BH, IW, BH, "Time for 2D mesh");
- value[num++] = new Fl_Output(2 * WB, 2 * WB + 12 * BH, IW, BH, "Time for 3D mesh");
+ value[num++] = new Fl_Output(2 * WB, 2 * WB + 11 * BH, IW, BH, "Time for 1D mesh");
+ value[num++] = new Fl_Output(2 * WB, 2 * WB + 12 * BH, IW, BH, "Time for 2D mesh");
+ value[num++] = new Fl_Output(2 * WB, 2 * WB + 13 * BH, IW, BH, "Time for 3D mesh");
- value[num] = new Fl_Output(2 * WB, 2 * WB + 13 * BH, IW, BH, "SICN");
+ value[num] = new Fl_Output(2 * WB, 2 * WB + 14 * BH, IW, BH, "SICN");
value[num]->tooltip("~ signed inverse condition number"); num++;
- value[num] = new Fl_Output(2 * WB, 2 * WB + 14 * BH, IW, BH, "Gamma");
+ value[num] = new Fl_Output(2 * WB, 2 * WB + 15 * BH, IW, BH, "Gamma");
value[num]->tooltip("~ inscribed_radius / circumscribed_radius (simplices)"); num++;
- value[num] = new Fl_Output(2 * WB, 2 * WB + 15 * BH, IW, BH, "Rho");
+ value[num] = new Fl_Output(2 * WB, 2 * WB + 16 * BH, IW, BH, "Rho");
value[num]->tooltip("~ min_edge_length / max_edge_length"); num++;
for(int i = 0; i < 3; i++){
int ww = 3 * FL_NORMAL_SIZE;
new Fl_Box
- (FL_NO_BOX, width - 3 * ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "Plot");
+ (FL_NO_BOX, width - 3 * ww - 2 * WB, 2 * WB + (14 + i) * BH, ww, BH, "Plot");
butt[2 * i] = new Fl_Button
- (width - 2 * ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "X-Y");
+ (width - 2 * ww - 2 * WB, 2 * WB + (14 + i) * BH, ww, BH, "X-Y");
butt[2 * i + 1] = new Fl_Button
- (width - ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "3D");
+ (width - ww - 2 * WB, 2 * WB + (14 + i) * BH, ww, BH, "3D");
}
static const QM_HISTO qmh0 = QMH_SICN_XY, qmh1 = QMH_SICN_3D, qmh2 = QMH_GAMMA_XY,
qmh3 = QMH_GAMMA_3D, qmh4 = QMH_RHO_XY, qmh5 = QMH_RHO_3D;
@@ -162,7 +163,8 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
value[num++] = new Fl_Output(2 * WB, 2 * WB + 7 * BH, IW, BH, "Hexahedra");
value[num++] = new Fl_Output(2 * WB, 2 * WB + 8 * BH, IW, BH, "Prisms");
value[num++] = new Fl_Output(2 * WB, 2 * WB + 9 * BH, IW, BH, "Pyramids");
- value[num++] = new Fl_Output(2 * WB, 2 * WB + 10 * BH, IW, BH, "Strings");
+ value[num++] = new Fl_Output(2 * WB, 2 * WB + 10 * BH, IW, BH, "Trihedra");
+ value[num++] = new Fl_Output(2 * WB, 2 * WB + 11 * BH, IW, BH, "Strings");
group[2]->end();
}
o->end();
@@ -216,10 +218,11 @@ void statisticsWindow::compute(bool elementQuality)
sprintf(label[num], "%g", s[10]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[11]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[12]); value[num]->value(label[num]); num++;
-
sprintf(label[num], "%g", s[13]); value[num]->value(label[num]); num++;
+
sprintf(label[num], "%g", s[14]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[15]); value[num]->value(label[num]); num++;
+ sprintf(label[num], "%g", s[16]); value[num]->value(label[num]); num++;
if(!elementQuality){
for(int i = 0; i < 6; i += 2) butt[i]->deactivate();
@@ -235,19 +238,18 @@ void statisticsWindow::compute(bool elementQuality)
}
else{
for(int i = 0; i < 6; i += 2) butt[i]->activate();
- sprintf(label[num], "%.4g (%.4g->%.4g)", s[17], s[18], s[19]);
+ sprintf(label[num], "%.4g (%.4g->%.4g)", s[18], s[19], s[20]);
value[num]->activate();
value[num]->value(label[num]); num++;
- sprintf(label[num], "%.4g (%.4g->%.4g)", s[20], s[21], s[22]);
+ sprintf(label[num], "%.4g (%.4g->%.4g)", s[21], s[22], s[23]);
value[num]->activate();
value[num]->value(label[num]); num++;
- sprintf(label[num], "%.4g (%.4g->%.4g)", s[23], s[24], s[25]);
+ sprintf(label[num], "%.4g (%.4g->%.4g)", s[24], s[25], s[26]);
value[num]->activate();
value[num]->value(label[num]); num++;
}
// post
- sprintf(label[num], "%g", s[26]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[27]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[28]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[29]); value[num]->value(label[num]); num++;
@@ -257,6 +259,8 @@ void statisticsWindow::compute(bool elementQuality)
sprintf(label[num], "%g", s[33]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[34]); value[num]->value(label[num]); num++;
sprintf(label[num], "%g", s[35]); value[num]->value(label[num]); num++;
+ sprintf(label[num], "%g", s[36]); value[num]->value(label[num]); num++;
+ sprintf(label[num], "%g", s[37]); value[num]->value(label[num]); num++;
static char mem[256];
long m = GetMemoryUsage();
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index e91f0bef8da156e8b2a2187a0db76e0cd4dbd04e..b1d3a1d9d2d281e3058bc801e92700aed7b49124 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -43,7 +43,7 @@ set(SRC
MFace.cpp
MElement.cpp MElementOctree.cpp
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
- MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp MSubElement.cpp
+ MHexahedron.cpp MPrism.cpp MPyramid.cpp MTrihedron.cpp MElementCut.cpp MSubElement.cpp
MZone.cpp MZoneBoundary.cpp
Cell.cpp CellComplex.cpp ChainComplex.cpp Homology.cpp Chain.cpp
Curvature.cpp
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 36309fa2cb8bebf5dc041ac5a57683f8890fb378..b11f31a513ac45dcf7a9e6c9471152de8bebbb2b 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -21,6 +21,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MElementCut.h"
#include "MElementOctree.h"
#include "discreteRegion.h"
@@ -587,6 +588,128 @@ bool GModel::setAllVolumesPositive()
return ok;
}
+static void addToMap(std::multimap< MFace , MElement *, Less_Face> &faceToElement, std::map< MElement *, std::vector < std::pair <MElement *, bool> > > &elToNeighbors, const MFace &face, MElement *el){
+ std::map< MFace , MElement *, Less_Face>::iterator fit = faceToElement.find(face);
+ if (fit == faceToElement.end()){
+ faceToElement.insert(std::pair< MFace , MElement *>(face, el));
+ } else { //We found the neighbor face outFace
+ faceToElement.insert(std::pair< MFace , MElement *>(face, el));
+ if (faceToElement.count(face) > 2)
+ Msg::Fatal("Topological fault: Face sharing two other faces. Element %i. Number of nodes %i. Count of faces: %i Three first nodes %i %i %i",el->getNum(),face.getNumVertices(),faceToElement.count(face),face.getVertex(0)->getNum(),face.getVertex(1)->getNum(),face.getVertex(2)->getNum());
+ MFace outFace = fit->first;
+ std::vector<std::pair<MElement *, bool> > &neigh = elToNeighbors[el];
+ for (size_t iN = 0; iN < neigh.size(); iN++)
+ if (neigh[iN].first == fit->second)
+ return;
+ int i0 = -1;
+ while (face.getVertex(0) != outFace.getVertex(++i0));
+ bool sameOrientation = face.getVertex(1) == outFace.getVertex((i0+1)%face.getNumVertices());
+ // if (sameOrientation) printf("Non-matching orientation between el %i and el %i for face IN %i %i %i out %i %i %i\n",el->getNum(),fit->second->getNum(),face.getVertex(0)->getNum(),face.getVertex(1)->getNum(),face.getVertex(2)->getNum(),outFace.getVertex(0)->getNum(),outFace.getVertex(1)->getNum(),outFace.getVertex(2)->getNum());
+ neigh.push_back(std::make_pair(fit->second, !sameOrientation));
+ elToNeighbors[fit->second].push_back(std::make_pair(el, !sameOrientation));
+ }
+}
+
+//static void checkConformity(std::multimap< MFace , MElement *, Less_Face> &faceToElement, std::map< MElement *, std::vector < std::pair <MElement *, bool> > > &elToNeighbors, const MFace &face, MElement *el){
+// //Check conformity
+// if (face.getNumVertices() == 4){
+// MVertex *v0 = face.getVertex(0);
+// MVertex *v1 = face.getVertex(1);
+// MVertex *v2 = face.getVertex(2);
+// MVertex *v3 = face.getVertex(3);
+// MFace f[5] = {MFace(v0, v1, v2), MFace(v0, v2, v3), MFace(v0, v1, v3), MFace(v1, v2, v3)};
+// if (faceToElement.find(f[0])!= faceToElement.end() || faceToElement.find(f[1])!= faceToElement.end() || faceToElement.find(f[2])!= faceToElement.end() || faceToElement.find(f[3])!= faceToElement.end() ){
+// //Check if the element is a trihedron used to recover conformity:
+// //A single element contains the two triangles, the quad and no other face
+// std::multiset<size_t> facesNeighborsId;
+// int nFaces=0;
+// for (int iFace = 0; iFace<4; iFace++){
+// int nbCount = faceToElement.count(f[iFace]);
+// if (nbCount != 0){
+// nFaces++;
+// if (nFaces > 3) Msg::Fatal("Topological fault: A quad face shares more than 3 triangles (element %i)",el->getNum());
+// std::pair <std::multimap< MFace , MElement *, Less_Face>::iterator, std::multimap< MFace , MElement *, Less_Face>::iterator> faceNeighbors;
+// faceNeighbors = faceToElement.equal_range(f[iFace]);
+// for (std::multimap< MFace , MElement *, Less_Face>::iterator itEl=faceNeighbors.first; itEl!=faceNeighbors.second; ++itEl)
+// facesNeighborsId.insert(itEl->second->getNum());
+// }
+// }
+// //In some cases, there can be 3 triangles matching a quad face and result in a valid trihedron (1 or 4 triangles not allowed)
+// if ((nFaces != 2 || nFaces != 3) && faceToElement.count(face)!= 2) Msg::Fatal("Nonconforming element which is not a trihedron (element %i nFaces %i nodes %i %i %i %i)",el->getNum(),nFaces,v0->getNum(),v1->getNum(),v2->getNum(),v3->getNum());
+// size_t maxOccurence=0;
+// for (std::multiset<size_t>::iterator iElId = facesNeighborsId.begin(); iElId != facesNeighborsId.end(); ++iElId){
+// maxOccurence = std::max(facesNeighborsId.count(*iElId), maxOccurence);
+// }
+// if (maxOccurence != 3) //A trihedron should have 3 and only 3 neighbors
+// Msg::Fatal("Nonconforming face for element: %i vertices %i %i %i %i. Number of neighbors: %i",el->getNum(), v0->getNum(), v1->getNum(), v2->getNum(), v3->getNum(),maxOccurence);
+// }
+// }
+//}
+
+static void checkConformity(std::multimap< MFace , MElement *, Less_Face> &faceToElement, std::map< MElement *, std::vector < std::pair <MElement *, bool> > > &elToNeighbors, MFace face, MElement *el){
+ int connectivity = faceToElement.count(face);
+ if (ElementType::ParentTypeFromTag(el->getType()) == TYPE_TRIH){
+ //Each face of a trihedron should exist twice (no face on the boundary)
+ if (connectivity != 2) Msg::Fatal("Non conforming trihedron %i (nb connections for a face %i)",el->getNum(), faceToElement.count(face));
+ }else{
+ //A face can exist twice (inside) or once (boundary)
+ if (connectivity != 2){
+ for (int iV = 0; iV < face.getNumVertices(); iV++)
+ if (face.getVertex(iV)->onWhat()->dim() == 3 || connectivity != 1){
+ for (int jV = 0; jV < face.getNumVertices(); jV++)
+ Msg::Info("Vertex %i",face.getVertex(jV)->getNum());
+ Msg::Fatal("Non conforming element %i (nb connections for a face %i vertex %i)",el->getNum(), connectivity, face.getVertex(iV)->getNum());
+ }
+ }
+ }
+}
+
+void GModel::setAllVolumesPositiveTopology()
+{
+ Msg::Info("Orienting volumes according to topology");
+ std::map< MElement *, std::vector < std::pair < MElement *, bool> > > elToNeighbors;
+ std::multimap< MFace , MElement *, Less_Face> faceToElement;
+
+ MElement *el;
+ for(riter it = regions.begin(); it != regions.end(); ++it){
+ for (unsigned int iEl = 0; iEl < (*it)->getNumMeshElements(); ++iEl) {
+ el = (*it)->getMeshElement(iEl);
+ for (int iFace = 0; iFace < el->getNumFaces(); iFace++){
+ addToMap(faceToElement, elToNeighbors, el->getFace(iFace), el);
+ }
+ }
+ }
+ for(riter it = regions.begin(); it != regions.end(); ++it){
+ for (unsigned int iEl = 0; iEl < (*it)->getNumMeshElements(); ++iEl) {
+ el = (*it)->getMeshElement(iEl);
+ for (int iFace = 0; iFace < el->getNumFaces(); iFace++){
+ checkConformity(faceToElement, elToNeighbors, el->getFace(iFace), el);
+ }
+ }
+ }
+ std::vector< std::pair <MElement *, bool > > queue;
+ std::set<MElement*> queued;
+ if ( (*regions.begin())->tetrahedra.size() == 0)
+ Msg::Fatal("setAllVolumePositiveTopology needs at least one tetrahedron to start");
+ el = (*regions.begin())->tetrahedra[0];
+ queue.push_back(std::make_pair(el, true));
+ for (size_t i = 0; i < queue.size(); i++){
+ el = queue[i].first;
+ if (!queue[i].second){
+ el->reverse();
+ // Msg::Info("Reverted element %i of type %i", el->getNum(), el->getType());
+ }
+ const std::vector < std::pair <MElement *, bool> > &neigh = elToNeighbors[el];
+ for (size_t iN = 0; iN < neigh.size(); iN++)
+ if (queued.count(neigh[iN].first) == 0){
+ queue.push_back(std::make_pair(neigh[iN].first, neigh[iN].second == queue[i].second));
+ // if (!(neigh[iN].second == queue[i].second))
+ // Msg::Info("Queuing element %i (%i) from el %i (%i)", neigh[iN].first->getNum(), neigh[iN].second, el->getNum(),queue[i].second);
+ queued.insert(neigh[iN].first);
+ }
+ }
+}
+
int GModel::adaptMesh(std::vector<int> technique,
std::vector<simpleFunction<double>* > f,
std::vector<std::vector<double> > parameters,
@@ -758,7 +881,7 @@ int GModel::getMeshStatus(bool countDiscrete)
(*it)->meshAttributes.method != MESH_NONE)) &&
((*it)->tetrahedra.size() ||(*it)->hexahedra.size() ||
(*it)->prisms.size() || (*it)->pyramids.size() ||
- (*it)->polyhedra.size())) return 3;
+ (*it)->polyhedra.size() || (*it)->trihedra.size())) return 3;
for(fiter it = firstFace(); it != lastFace(); ++it)
if((countDiscrete || ((*it)->geomType() != GEntity::DiscreteSurface &&
(*it)->meshAttributes.method != MESH_NONE)) &&
@@ -803,12 +926,12 @@ int GModel::getNumMeshParentElements()
return n;
}
-int GModel::getNumMeshElements(unsigned c[5])
+int GModel::getNumMeshElements(unsigned c[6])
{
- c[0] = 0; c[1] = 0; c[2] = 0; c[3] = 0; c[4] = 0;
+ c[0] = 0; c[1] = 0; c[2] = 0; c[3] = 0; c[4] = 0; c[5] = 0;
for(riter it = firstRegion(); it != lastRegion(); ++it)
(*it)->getNumMeshElements(c);
- if(c[0] + c[1] + c[2] + c[3] + c[4]) return 3;
+ if(c[0] + c[1] + c[2] + c[3] + c[4] + c[5]) return 3;
for(fiter it = firstFace(); it != lastFace(); ++it)
(*it)->getNumMeshElements(c);
if(c[0] + c[1] + c[2]) return 2;
@@ -970,6 +1093,7 @@ void GModel::removeInvisibleElements()
removeInvisible((*it)->hexahedra, all);
removeInvisible((*it)->prisms, all);
removeInvisible((*it)->pyramids, all);
+ removeInvisible((*it)->trihedra, all);
removeInvisible((*it)->polyhedra, all);
(*it)->deleteVertexArrays();
}
@@ -1148,7 +1272,7 @@ void GModel::_storeElementsInEntities(std::map< int, std::vector<MElement* > >&
else _addElements(f->polygons, it->second);
}
break;
- case TYPE_TET: case TYPE_HEX: case TYPE_PYR: case TYPE_PRI: case TYPE_POLYH:
+ case TYPE_TET: case TYPE_HEX: case TYPE_PYR: case TYPE_TRIH: case TYPE_PRI: case TYPE_POLYH:
{
GRegion *r = getRegionByTag(it->first);
if(!r){
@@ -1159,6 +1283,7 @@ void GModel::_storeElementsInEntities(std::map< int, std::vector<MElement* > >&
else if(type == TYPE_HEX) _addElements(r->hexahedra, it->second);
else if(type == TYPE_PRI) _addElements(r->prisms, it->second);
else if(type == TYPE_PYR) _addElements(r->pyramids, it->second);
+ else if(type == TYPE_TRIH) _addElements(r->trihedra, it->second);
else _addElements(r->polyhedra, it->second);
}
break;
@@ -1190,6 +1315,7 @@ void GModel::_associateEntityWithMeshVertices()
_associateEntityWithElementVertices(*it, (*it)->hexahedra);
_associateEntityWithElementVertices(*it, (*it)->prisms);
_associateEntityWithElementVertices(*it, (*it)->pyramids);
+ _associateEntityWithElementVertices(*it, (*it)->trihedra);
_associateEntityWithElementVertices(*it, (*it)->polyhedra);
}
for(fiter it = firstFace(); it != lastFace(); ++it){
@@ -1254,6 +1380,7 @@ void GModel::pruneMeshVertexAssociations()
_associateEntityWithElementVertices(*it, (*it)->hexahedra, true);
_associateEntityWithElementVertices(*it, (*it)->prisms, true);
_associateEntityWithElementVertices(*it, (*it)->pyramids, true);
+ _associateEntityWithElementVertices(*it, (*it)->trihedra, true);
_associateEntityWithElementVertices(*it, (*it)->polyhedra, true);
}
for(fiter it = _chainFaces.begin(); it != _chainFaces.end(); ++it){
@@ -1355,7 +1482,7 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
}
GModel *gm = new GModel();
- std::map<int, std::vector<MElement*> > elements[10];
+ std::map<int, std::vector<MElement*> > elements[11];
std::map<int, std::map<int, std::string> > physicals[4];
std::vector<MVertex*> vertexVector;
@@ -1428,8 +1555,9 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
case TYPE_HEX : elements[5][elementary[i]].push_back(e); break;
case TYPE_PRI : elements[6][elementary[i]].push_back(e); break;
case TYPE_PYR : elements[7][elementary[i]].push_back(e); break;
- case TYPE_POLYG: elements[8][elementary[i]].push_back(e); break;
- case TYPE_POLYH: elements[9][elementary[i]].push_back(e); break;
+ case TYPE_TRIH : elements[8][elementary[i]].push_back(e); break;
+ case TYPE_POLYG: elements[9][elementary[i]].push_back(e); break;
+ case TYPE_POLYH: elements[10][elementary[i]].push_back(e); break;
default : Msg::Error("Wrong type of element"); delete gm; return 0;
}
int dim = e->getDim();
@@ -1578,7 +1706,7 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
return 0;
}
- std::map<int, std::vector<MElement*> > elements[10];
+ std::map<int, std::vector<MElement*> > elements[11];
for(unsigned int i = 0; i < entities.size(); i++){
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(j);
@@ -1601,8 +1729,9 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
case TYPE_HEX: elements[5][entities[i]->tag()].push_back(e2); break;
case TYPE_PRI: elements[6][entities[i]->tag()].push_back(e2); break;
case TYPE_PYR: elements[7][entities[i]->tag()].push_back(e2); break;
- case TYPE_POLYG: elements[8][entities[i]->tag()].push_back(e2); break;
- case TYPE_POLYH: elements[9][entities[i]->tag()].push_back(e2); break;
+ case TYPE_TRIH: elements[8][entities[i]->tag()].push_back(e2); break;
+ case TYPE_POLYG: elements[9][entities[i]->tag()].push_back(e2); break;
+ case TYPE_POLYH: elements[10][entities[i]->tag()].push_back(e2); break;
}
}
else
@@ -1824,6 +1953,7 @@ void GModel::makeDiscreteRegionsSimplyConnected()
case TYPE_HEX: r->hexahedra.push_back((MHexahedron*)e2); break;
case TYPE_PRI: r->prisms.push_back((MPrism*)e2); break;
case TYPE_PYR: r->pyramids.push_back((MPyramid*)e2); break;
+ case TYPE_TRIH: r->trihedra.push_back((MTrihedron*)e2); break;
}
}
r->mesh_vertices.insert
@@ -2299,6 +2429,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
std::vector<MHexahedron*> newHexahedra;
std::vector<MPrism*> newPrisms;
std::vector<MPyramid*> newPyramids;
+ std::vector<MTrihedron*> newTrihedra;
for (unsigned int i = 0; i < gr->getNumMeshElements(); ++i){
MElement *e = gr->getMeshElement(i);
std::vector<MVertex *> v;
@@ -2317,6 +2448,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
case TYPE_HEX: newHexahedra.push_back((MHexahedron*)e2); break;
case TYPE_PRI: newPrisms.push_back((MPrism*)e2); break;
case TYPE_PYR: newPyramids.push_back((MPyramid*)e2); break;
+ case TYPE_TRIH: newTrihedra.push_back((MTrihedron*)e2); break;
}
}
gr->deleteVertexArrays();
@@ -2324,16 +2456,18 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
for(unsigned int i = 0; i < gr->hexahedra.size(); i++) delete gr->hexahedra[i];
for(unsigned int i = 0; i < gr->prisms.size(); i++) delete gr->prisms[i];
for(unsigned int i = 0; i < gr->pyramids.size(); i++) delete gr->pyramids[i];
+ for(unsigned int i = 0; i < gr->trihedra.size(); i++) delete gr->trihedra[i];
gr->tetrahedra = newTetrahedra;
gr->hexahedra = newHexahedra;
gr->prisms = newPrisms;
gr->pyramids = newPyramids;
+ gr->trihedra = newTrihedra;
}
Msg::Debug("Done creating topology for edges");
}
-void makeSimplyConnected(std::map<int, std::vector<MElement*> > elements[10])
+void makeSimplyConnected(std::map<int, std::vector<MElement*> > elements[11])
{
//only for tetras and triangles
Msg::Info("Make simply connected regions and surfaces.");
@@ -2502,7 +2636,7 @@ GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
else
CTX::instance()->mesh.saveTri = 0;
- std::map<int, std::vector<MElement*> > elements[10];
+ std::map<int, std::vector<MElement*> > elements[11];
std::map<int, std::map<int, std::string> > physicals[4];
std::map<int, MVertex*> vertexMap;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 9996a42c4149be4115fb99b2e8f46cf65a297c41..11411531bbfee66df6c74329fd7142ca4901d82c 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -362,7 +362,7 @@ class GModel
// get the number of each type of element in the mesh at the largest
// dimension and return the dimension
- int getNumMeshElements(unsigned c[5]);
+ int getNumMeshElements(unsigned c[6]);
// access a mesh element by coordinates (using an octree search)
MElement *getMeshElementByCoord(SPoint3 &p, int dim=-1, bool strict=true);
@@ -467,6 +467,7 @@ class GModel
// Ensure that the Jacobian of all volume elements is positive
bool setAllVolumesPositive();
+ void setAllVolumesPositiveTopology();
// make the mesh a high order mesh at order N
// linear is 1 if the high order points are not placed on the geometry of the model
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index 0ee171ed0ed307882a1d477b6a4ab5c449c500e8..6562f53a3a88089f774b9207d41e24ac801102b3 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/Geo/GModelIO_MSH.cpp
@@ -17,6 +17,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "StringUtils.h"
#include "discreteVertex.h"
#include "discreteEdge.h"
@@ -170,7 +171,7 @@ int GModel::readMSH(const std::string &name)
double version = 0.;
bool binary = false, swap = false, postpro = false;
int minVertex = 0;
- std::map<int, std::vector<MElement*> > elements[10];
+ std::map<int, std::vector<MElement*> > elements[11];
while(1) {
@@ -406,8 +407,9 @@ int GModel::readMSH(const std::string &name)
case TYPE_HEX: elements[5][entity].push_back(element); break;
case TYPE_PRI: elements[6][entity].push_back(element); break;
case TYPE_PYR: elements[7][entity].push_back(element); break;
- case TYPE_POLYG: elements[8][entity].push_back(element); break;
- case TYPE_POLYH: elements[9][entity].push_back(element); break;
+ case TYPE_TRIH: elements[8][entity].push_back(element); break;
+ case TYPE_POLYG: elements[9][entity].push_back(element); break;
+ case TYPE_POLYH: elements[10][entity].push_back(element); break;
}
if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
@@ -665,6 +667,9 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, this, *it, (*it)->pyramids, saveAll, saveSinglePartition,
binary);
+ for(riter it = firstRegion(); it != lastRegion(); ++it)
+ writeElementsMSH(fp, this, *it, (*it)->trihedra, saveAll, saveSinglePartition,
+ binary);
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, this, *it, (*it)->triangles, saveAll, saveSinglePartition,
binary);
diff --git a/Geo/GModelIO_MSH2.cpp b/Geo/GModelIO_MSH2.cpp
index a99fd2f459fa305b03a2c39abba6fdb867cf5c91..db97baf5987f4c80aca3195f16738ba1a35e41e1 100644
--- a/Geo/GModelIO_MSH2.cpp
+++ b/Geo/GModelIO_MSH2.cpp
@@ -847,6 +847,7 @@ static int getNumElementsMSH(GModel *m, bool saveAll, int saveSinglePartition)
if((*it)->polyhedra[i]->ownsParent())
n += (saveAll ? 1 : (*it)->physicals.size());
}
+ n -= (*it)->trihedra.size();
}
return n;
}
diff --git a/Geo/GModelVertexArrays.cpp b/Geo/GModelVertexArrays.cpp
index d3a6c0992482f700320e8188cdb092e3e851b76e..f45a2f7f1a999fc1a053a604942090cbbdfe31ec 100644
--- a/Geo/GModelVertexArrays.cpp
+++ b/Geo/GModelVertexArrays.cpp
@@ -14,6 +14,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MElementCut.h"
#include "Context.h"
#include "VertexArray.h"
@@ -56,6 +57,7 @@ static unsigned int getColorByElement(MElement *ele)
case TYPE_HEX: return CTX::instance()->color.mesh.hexahedron;
case TYPE_PRI: return CTX::instance()->color.mesh.prism;
case TYPE_PYR: return CTX::instance()->color.mesh.pyramid;
+ case TYPE_TRIH: return CTX::instance()->color.mesh.trihedron;
default: return CTX::instance()->color.mesh.vertex;
}
}
@@ -353,7 +355,7 @@ class initMeshGRegion {
for(unsigned int i = 0; i < r->polyhedra.size(); i++)
numLP += 2 * r->polyhedra[i]->getNumEdges();
num += (12 * r->tetrahedra.size() + 24 * r->hexahedra.size() +
- 18 * r->prisms.size() + 16 * r->pyramids.size() + numLP) / 4;
+ 18 * r->prisms.size() + 16 * r->pyramids.size() + 10 * r->trihedra.size() + numLP) / 4;
num = _estimateIfClipped(num);
if(CTX::instance()->mesh.explode != 1.) num *= 4;
if(_curved) num *= 2;
@@ -368,7 +370,7 @@ class initMeshGRegion {
for(unsigned int i = 0; i < r->polyhedra.size(); i++)
numFP += r->polyhedra[i]->getNumFaces();
num += (4 * r->tetrahedra.size() + 12 * r->hexahedra.size() +
- 8 * r->prisms.size() + 6 * r->pyramids.size() + numFP) / 2;
+ 8 * r->prisms.size() + 6 * r->pyramids.size() + 4 * r->trihedra.size() + numFP) / 2;
num = _estimateIfClipped(num);
if(CTX::instance()->mesh.explode != 1.) num *= 2;
if(_curved) num *= 4;
@@ -384,7 +386,8 @@ class initMeshGRegion {
(CTX::instance()->mesh.tetrahedra && areAllElementsVisible(r->tetrahedra) &&
CTX::instance()->mesh.hexahedra && areAllElementsVisible(r->hexahedra) &&
CTX::instance()->mesh.prisms && areAllElementsVisible(r->prisms) &&
- CTX::instance()->mesh.pyramids && areAllElementsVisible(r->pyramids));
+ CTX::instance()->mesh.pyramids && areAllElementsVisible(r->pyramids) &&
+ CTX::instance()->mesh.trihedra && areAllElementsVisible(r->trihedra));
bool edg = CTX::instance()->mesh.volumesEdges;
bool fac = CTX::instance()->mesh.volumesFaces;
@@ -392,13 +395,15 @@ class initMeshGRegion {
_curved = (areSomeElementsCurved(r->tetrahedra) ||
areSomeElementsCurved(r->hexahedra) ||
areSomeElementsCurved(r->prisms) ||
- areSomeElementsCurved(r->pyramids));
+ areSomeElementsCurved(r->pyramids) ||
+ areSomeElementsCurved(r->trihedra));
r->va_lines = new VertexArray(2, _estimateNumLines(r));
r->va_triangles = new VertexArray(3, _estimateNumTriangles(r));
if(CTX::instance()->mesh.tetrahedra) addElementsInArrays(r, r->tetrahedra, edg, fac);
if(CTX::instance()->mesh.hexahedra) addElementsInArrays(r, r->hexahedra, edg, fac);
if(CTX::instance()->mesh.prisms) addElementsInArrays(r, r->prisms, edg, fac);
if(CTX::instance()->mesh.pyramids) addElementsInArrays(r, r->pyramids, edg, fac);
+ if(CTX::instance()->mesh.trihedra) addElementsInArrays(r, r->trihedra, edg, fac);
addElementsInArrays(r, r->polyhedra, edg, fac);
r->va_lines->finalize();
r->va_triangles->finalize();
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 08de35b82016fc2aceaa6c6cbf9d20fd4c59048a..778a04dfcd94df1fcac94f0bf4013a361bf65d99 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -11,6 +11,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MElementCut.h"
#include "GmshMessage.h"
#include "VertexArray.h"
@@ -45,8 +46,11 @@ void GRegion::deleteMesh()
prisms.clear();
for(unsigned int i = 0; i < pyramids.size(); i++) delete pyramids[i];
pyramids.clear();
+ for(unsigned int i = 0; i < trihedra.size(); i++) delete trihedra[i];
+ trihedra.clear();
for(unsigned int i = 0; i < polyhedra.size(); i++) delete polyhedra[i];
polyhedra.clear();
+
deleteVertexArrays();
model()->destroyMeshCaches();
}
@@ -54,7 +58,7 @@ void GRegion::deleteMesh()
unsigned int GRegion::getNumMeshElements()
{
return tetrahedra.size() + hexahedra.size() + prisms.size() + pyramids.size() +
- polyhedra.size();
+ trihedra.size() + polyhedra.size();
}
unsigned int GRegion::getNumMeshParentElements()
@@ -72,7 +76,8 @@ void GRegion::getNumMeshElements(unsigned *const c) const
c[1] += hexahedra.size();
c[2] += prisms.size();
c[3] += pyramids.size();
- c[4] += polyhedra.size();
+ c[4] += trihedra.size();
+ c[5] += polyhedra.size();
}
MElement *const *GRegion::getStartElementType(int type) const
@@ -91,6 +96,9 @@ MElement *const *GRegion::getStartElementType(int type) const
if(pyramids.empty()) return 0; // msvc would throw an exception
return reinterpret_cast<MElement *const *>(&pyramids[0]);
case 4:
+ if(trihedra.empty()) return 0;
+ return reinterpret_cast<MElement *const *>(&trihedra[0]);
+ case 5:
if(polyhedra.empty()) return 0;
return reinterpret_cast<MElement *const *>(&polyhedra[0]);
}
@@ -109,9 +117,14 @@ MElement *GRegion::getMeshElement(unsigned int index)
pyramids.size())
return pyramids[index - tetrahedra.size() - hexahedra.size() - prisms.size()];
else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
- pyramids.size() + polyhedra.size())
+ pyramids.size() + trihedra.size())
+ return trihedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() -
+ pyramids.size()];
+ else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
+ pyramids.size() + trihedra.size() + polyhedra.size())
return polyhedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() -
- pyramids.size()];
+ pyramids.size() - trihedra.size()];
+
return 0;
}
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index d05cdf6e845737cf5277528ae5895e90b32d93a0..2812a8a8f7c8fb253295c24095c65e1401a9a944 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -12,6 +12,7 @@
#include <stdio.h>
#include "GEntity.h"
#include "boundaryLayersData.h"
+#include "GmshDefines.h"
class MElement;
class MTetrahedron;
@@ -19,6 +20,7 @@ class MHexahedron;
class MPrism;
class MPyramid;
class MPolyhedron;
+class MTrihedron;
class ExtrudeParams;
class GRegionCompound;
class BoundaryLayerColumns;
@@ -90,7 +92,7 @@ class GRegion : public GEntity {
virtual void writeGEO(FILE *fp);
// number of types of elements
- int getNumElementTypes() const { return 5; }
+ int getNumElementTypes() const { return 6; }
// get total/by-type number of elements in the mesh
unsigned int getNumMeshElements();
@@ -131,13 +133,41 @@ class GRegion : public GEntity {
std::vector<MHexahedron*> hexahedra;
std::vector<MPrism*> prisms;
std::vector<MPyramid*> pyramids;
+ std::vector<MTrihedron*> trihedra;
std::vector<MPolyhedron*> polyhedra;
void addTetrahedron(MTetrahedron *t){ tetrahedra.push_back(t); }
+
void addHexahedron(MHexahedron *h){ hexahedra.push_back(h); }
void addPrism(MPrism *p){ prisms.push_back(p); }
void addPyramid(MPyramid *p){ pyramids.push_back(p); }
void addPolyhedron(MPolyhedron *p){ polyhedra.push_back(p); }
+ void addTrihedron(MTrihedron *t){ trihedra.push_back(t); }
+ void addElement(int type, MElement *e){
+ switch (type){
+ case TYPE_TET:
+ addTetrahedron((MTetrahedron*) e);
+ break;
+ case TYPE_HEX:
+ addHexahedron((MHexahedron*) e);
+ break;
+ case TYPE_PRI:
+ addPrism((MPrism*) e);
+ break;
+ case TYPE_PYR:
+ addPyramid((MPyramid*) e);
+ break;
+ case TYPE_TRIH:
+ addTrihedron((MTrihedron*) e);
+ break;
+ case TYPE_POLYH:
+ addPolyhedron((MPolyhedron*) e);
+ break;
+ default:
+ Msg::Fatal("Trying to add unsupported element");
+ }
+ }
+
// get the boundary layer columns
BoundaryLayerColumns *getColumns () {return &_columns;}
};
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index a678a2be8f98d42b05d5b190fadcb884a556ace6..61e6e3d2a5f68a11b209a4f4bcf67e8cf6c0abed 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/Geo/GeomMeshMatcher.cpp
@@ -22,6 +22,7 @@
#include "MVertex.h"
#include "MLine.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MPrism.h"
#include "MPoint.h"
#include "MHexahedron.h"
@@ -722,6 +723,7 @@ void copy_elements (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_me
copy_elements<MHexahedron> (dest->hexahedra ,orig->hexahedra ,_mesh_to_geom);
copy_elements<MPrism> (dest->prisms ,orig->prisms ,_mesh_to_geom);
copy_elements<MPyramid> (dest->pyramids ,orig->pyramids ,_mesh_to_geom);
+ copy_elements<MTrihedron> (dest->trihedra ,orig->trihedra ,_mesh_to_geom);
}
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 0424d058c4fadc1b9d3399ce5a6646b3a79c3a91..533d163e7e1fa182c7c355f954270c14f702eac7 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -17,6 +17,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MElementCut.h"
#include "MSubElement.h"
#include "GEntity.h"
@@ -116,6 +117,9 @@ MElement* MElement::createElement(int tag, const std::vector<MVertex*> &vertices
else
return new MPrismN(vertices, order, num, part);
+ case TYPE_TRIH:
+ return new MTrihedron(vertices, num, part);
+
case TYPE_HEX:
if (order == 1)
return new MHexahedron(vertices, num, part);
@@ -1547,6 +1551,7 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_PYR_53 : if(name) *name = "Pyramid 53"; return 5 + 8*6;
case MSH_PYR_61 : if(name) *name = "Pyramid 61"; return 5 + 8*7;
case MSH_PYR_69 : if(name) *name = "Pyramid 69"; return 5 + 8*8;
+ case MSH_TRIH_4 : if(name) *name = "Trihedron 4"; return 4;
case MSH_POLYH_ : if(name) *name = "Polyhedron"; return 0;
case MSH_PNT_SUB : if(name) *name = "Point Xfem"; return 1;
case MSH_LIN_SUB : if(name) *name = "Line Xfem"; return 2;
@@ -1757,6 +1762,7 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_PYR_204: return new MPyramidN(v, 7, num, part);
case MSH_PYR_285: return new MPyramidN(v, 8, num, part);
case MSH_PYR_385: return new MPyramidN(v, 9, num, part);
+ case MSH_TRIH_4: return new MTrihedron(v, num, part);
default: return 0;
}
}
diff --git a/Geo/MTrihedron.cpp b/Geo/MTrihedron.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d9b2193145fbf1c815421da008665d6f4cf60b09
--- /dev/null
+++ b/Geo/MTrihedron.cpp
@@ -0,0 +1,3 @@
+#include "GmshConfig.h"
+#include "MTrihedron.h"
+
diff --git a/Geo/MTrihedron.h b/Geo/MTrihedron.h
new file mode 100644
index 0000000000000000000000000000000000000000..94a3f9c6d985e60dceedcfb38347dc433c86e32f
--- /dev/null
+++ b/Geo/MTrihedron.h
@@ -0,0 +1,162 @@
+#ifndef _MTRIHEDRON_H_
+#define _MTRIHEDRON_H_
+
+#include "MElement.h"
+
+/*
+ * MTrihedron
+ * A MTrihedron is a plane element composed of
+ * a quadrangle and two triangles.
+ * It serves as an interface between two non-conforming
+ * elements
+ *
+ * v
+ * ^
+ * |
+ * 3-----------2
+ * |'\ | |
+ * | '\ | |
+ * | +---- | --> u
+ * | '\ |
+ * | '\ |
+ * 0-----------1
+ *
+ */
+
+class MTrihedron : public MElement {
+ protected:
+ MVertex *_v[4];
+ void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v[0] = _v[edges_trihedron(num, 0)];
+ v[1] = _v[edges_trihedron(num, 1)];
+ }
+ void _getFaceVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ if(num > 0) {
+ v[0] = _v[faces_trihedron(num, 0)];
+ v[1] = _v[faces_trihedron(num, 1)];
+ v[2] = _v[faces_trihedron(num, 2)];
+ }
+ else {
+ v[0] = _v[0];
+ v[1] = _v[1];
+ v[2] = _v[2];
+ v[3] = _v[3];
+ }
+ }
+
+ public:
+ MTrihedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0)
+ : MElement(num, part)
+ {
+ _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3;
+ }
+ MTrihedron(const std::vector<MVertex*> &v, int num=0, int part=0)
+ : MElement(num, part)
+ {
+ for(int i = 0; i < 4; i++) _v[i] = v[i];
+ }
+ ~MTrihedron(){}
+ virtual int getDim() const { return 3; } //Can have a volume...
+ virtual int getNumVertices() const { return 4; }
+ virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual const MVertex *getVertex(int num) const { return _v[num]; }
+ virtual void setVertex(int num, MVertex *v){ _v[num] = v; }
+ virtual int getNumEdges()const{ return 5; }
+ virtual MEdge getEdge(int num) const
+ {
+ return MEdge(_v[edges_trihedron(num, 0)], _v[edges_trihedron(num, 1)]);
+ }
+ virtual int getNumEdgesRep(bool curved){ return 5; }
+ virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
+ {
+ MEdge e(getEdge(num));
+ _getEdgeRep(e.getVertex(0), e.getVertex(1), x, y, z, n, 0);
+ }
+ virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v.resize(2);
+ _getEdgeVertices(num, v);
+ }
+ virtual int getNumFaces(){ return 3; }
+ virtual MFace getFace(int num)
+ {
+ if(num > 0)
+ return MFace(_v[faces_trihedron(num, 0)],
+ _v[faces_trihedron(num, 1)],
+ _v[faces_trihedron(num, 2)]);
+ else
+ return MFace(_v[0], _v[1], _v[2], _v[3]);
+ }
+
+ virtual int getNumFacesRep(bool curved){ return 2; }
+ virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
+ {
+ static const int f[2][3] = {
+ {0, 1, 3}, {1, 2, 3}
+ };
+ _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
+ x, y, z, n);
+ }
+ virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v.resize((num == 0) ? 4 : 3);
+ _getFaceVertices(num, v);
+ }
+ virtual int getType() const { return TYPE_TRIH; }
+ virtual int getTypeForMSH() const { return MSH_TRIH_4; }
+
+ virtual void reverse()
+ {
+ MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
+ }
+ virtual int getVolumeSign(){return 0;};
+ virtual double getVolume(){return 0;};
+ virtual bool setVolumePositive(){return false;};
+ virtual void getNode(int num, double &u, double &v, double &w) const
+ {
+ w = 0;
+ switch(num) {
+ case 0 : u = -1.; v = -1.; break;
+ case 1 : u = 1.; v = -1.; break;
+ case 2 : u = 1.; v = 1.; break;
+ case 3 : u = -1.; v = 1.; break;
+ default: u = 0.; v = 0.; break;
+ }
+ }
+ virtual SPoint3 barycenterUVW() const
+ {
+ return SPoint3(0., 0., 0.);
+ }
+ virtual bool isInside(double u, double v, double w) const
+ {
+ double tol = _isInsideTolerance;
+ if(u < -(1. + tol) || v < -(1. + tol) || u > (1. + tol) || v > (1. + tol) ||
+ fabs(w) > tol)
+ return false;
+ return true;
+ }
+ static int edges_trihedron(const int edge, const int vert)
+ {
+ static const int e[5][2] = {
+ {0, 1},
+ {1, 2},
+ {2, 3},
+ {3, 0},
+ {1, 3},
+ };
+ return e[edge][vert];
+ }
+ static int faces_trihedron(const int face, const int vert)
+ {
+ static const int f[3][4] = {
+ {0, 1, 2, 3},
+ {0, 3, 1, -1},
+ {1, 3, 2, -1},
+ };
+ return f[face][vert];
+ }
+};
+
+#endif
diff --git a/Graphics/drawMesh.cpp b/Graphics/drawMesh.cpp
index 13410eeb23f58be3cd1e0de831cbd4588d6d1257..3c91e8517259a58966dba7746321e79a0ef7a381 100644
--- a/Graphics/drawMesh.cpp
+++ b/Graphics/drawMesh.cpp
@@ -17,6 +17,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MElementCut.h"
#include "Context.h"
#include "OS.h"
@@ -585,6 +586,10 @@ class drawMeshGRegion {
drawElementLabels(_ctx, r, r->pyramids, CTX::instance()->mesh.volumesFaces ||
CTX::instance()->mesh.surfacesFaces,
CTX::instance()->color.mesh.line);
+ if(CTX::instance()->mesh.trihedra)
+ drawElementLabels(_ctx, r, r->trihedra, CTX::instance()->mesh.volumesFaces ||
+ CTX::instance()->mesh.surfacesFaces,
+ CTX::instance()->color.mesh.line);
drawElementLabels(_ctx, r, r->polyhedra, CTX::instance()->mesh.volumesFaces ||
CTX::instance()->mesh.surfacesFaces,
CTX::instance()->color.mesh.line);
@@ -598,6 +603,7 @@ class drawMeshGRegion {
if(CTX::instance()->mesh.hexahedra) drawVerticesPerElement(_ctx, r, r->hexahedra);
if(CTX::instance()->mesh.prisms) drawVerticesPerElement(_ctx, r, r->prisms);
if(CTX::instance()->mesh.pyramids) drawVerticesPerElement(_ctx, r, r->pyramids);
+ if(CTX::instance()->mesh.trihedra) drawVerticesPerElement(_ctx, r, r->trihedra);
drawVerticesPerElement(_ctx, r, r->polyhedra);
}
}
@@ -607,6 +613,7 @@ class drawMeshGRegion {
if(CTX::instance()->mesh.hexahedra) drawBarycentricDual(r->hexahedra);
if(CTX::instance()->mesh.prisms) drawBarycentricDual(r->prisms);
if(CTX::instance()->mesh.pyramids) drawBarycentricDual(r->pyramids);
+ if(CTX::instance()->mesh.trihedra) drawBarycentricDual(r->trihedra);
drawBarycentricDual(r->polyhedra);
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 779b4885ef9f71449497886f2674a778d8b2d4dd..82b2008ab7e7583f51648368cd5d05b2798f01ed 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -111,11 +111,12 @@ void GetStatistics(double stat[50], double quality[3][100])
stat[10] += (*it)->hexahedra.size();
stat[11] += (*it)->prisms.size();
stat[12] += (*it)->pyramids.size();
+ stat[13] += (*it)->trihedra.size();
}
- stat[13] = CTX::instance()->meshTimer[0];
- stat[14] = CTX::instance()->meshTimer[1];
- stat[15] = CTX::instance()->meshTimer[2];
+ stat[14] = CTX::instance()->meshTimer[0];
+ stat[15] = CTX::instance()->meshTimer[1];
+ stat[16] = CTX::instance()->meshTimer[2];
if(quality){
for(int i = 0; i < 3; i++)
@@ -125,7 +126,7 @@ void GetStatistics(double stat[50], double quality[3][100])
double gamma = 0., gammaMin = 1., gammaMax = 0.;
double rho = 0., rhoMin = 1., rhoMax = 0.;
- double N = stat[9] + stat[10] + stat[11] + stat[12];
+ double N = stat[9] + stat[10] + stat[11] + stat[12] + stat[13];
if(N){ // if we have 3D elements
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
GetQualityMeasure((*it)->tetrahedra, gamma, gammaMin, gammaMax,
@@ -148,25 +149,25 @@ void GetStatistics(double stat[50], double quality[3][100])
}
}
if(N){
- stat[17] = minSICN / N; stat[18] = minSICNMin; stat[19] = minSICNMax;
- stat[20] = gamma / N; stat[21] = gammaMin; stat[22] = gammaMax;
- stat[23] = rho / N; stat[24] = rhoMin; stat[25] = rhoMax;
+ stat[18] = minSICN / N; stat[19] = minSICNMin; stat[20] = minSICNMax;
+ stat[21] = gamma / N; stat[22] = gammaMin; stat[23] = gammaMax;
+ stat[25] = rho / N; stat[25] = rhoMin; stat[26] = rhoMax;
}
}
#if defined(HAVE_POST)
- stat[26] = PView::list.size();
+ stat[27] = PView::list.size();
for(unsigned int i = 0; i < PView::list.size(); i++) {
PViewData *data = PView::list[i]->getData(true);
- stat[27] += data->getNumPoints();
- stat[28] += data->getNumLines();
- stat[29] += data->getNumTriangles();
- stat[30] += data->getNumQuadrangles();
- stat[31] += data->getNumTetrahedra();
- stat[32] += data->getNumHexahedra();
- stat[33] += data->getNumPrisms();
- stat[34] += data->getNumPyramids();
- stat[35] += data->getNumStrings2D() + data->getNumStrings3D();
+ stat[28] += data->getNumPoints();
+ stat[29] += data->getNumLines();
+ stat[30] += data->getNumTriangles();
+ stat[31] += data->getNumQuadrangles();
+ stat[32] += data->getNumTetrahedra();
+ stat[33] += data->getNumHexahedra();
+ stat[34] += data->getNumPrisms();
+ stat[35] += data->getNumPyramids();
+ stat[36] += data->getNumStrings2D() + data->getNumStrings3D();
}
#endif
}
@@ -534,9 +535,9 @@ static void Mesh3D(GModel *m)
Recombinator rec;
rec.execute(gr);
Supplementary sup;
- sup.execute(gr);
+ // sup.execute(gr);
PostOp post;
- post.execute(gr,0);
+ post.execute(gr,0, true); //0: no pyramid, 1: single-step, 2: two-steps (conforming), true: fill non-conformities with trihedra
nb_elements_recombination += post.get_nb_elements();
nb_hexa_recombination += post.get_nb_hexahedra();
vol_element_recombination += post.get_vol_elements();
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 887b30eb1bd74a90940c5f44601da9bceb41a743..bdb33df6fb0d4fba956e9d5a4b44016e31f443f4 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -19,6 +19,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MElementCut.h"
#include "MPoint.h"
#include "partitionVertex.h"
@@ -219,6 +220,8 @@ int RenumberMeshElements(std::vector<MElement*> &elements, meshPartitionOptions
gr->prisms.push_back((MPrism*)(*it));
else if ((*it)->getType() == TYPE_PYR)
gr->pyramids.push_back((MPyramid*)(*it));
+ else if ((*it)->getType() == TYPE_TRIH)
+ gr->trihedra.push_back((MTrihedron*)(*it));
}
tmp_model->add(gr);
RenumberMesh(tmp_model, options, elements);
@@ -556,7 +559,8 @@ int MakeGraph(GModel *const model, Graph &graph, meshPartitionOptions &options,
ElemTypeHexa = 1,
ElemTypePrism = 2,
ElemTypePyramid = 3,
- ElemTypePolyh = 4
+ ElemTypeTrih = 4,
+ ElemTypePolyh = 5
};
enum {
ElemTypeTri = 0,
@@ -574,7 +578,7 @@ int MakeGraph(GModel *const model, Graph &graph, meshPartitionOptions &options,
//--Get the dimension of the mesh and count the numbers of elements
- unsigned numElem[5];
+ unsigned numElem[6];
const int meshDim = model->getNumMeshElements(numElem);
if(meshDim < 2) {
Msg::Error("No mesh elements were found");
@@ -613,11 +617,11 @@ int MakeGraph(GModel *const model, Graph &graph, meshPartitionOptions &options,
const int numGrVert =
numElem[ElemTypeTetra] + numElem[ElemTypeHexa] +
numElem[ElemTypePrism] + numElem[ElemTypePyramid] +
- numElem[ElemTypePolyh];
+ + numElem[ElemTypeTrih] + numElem[ElemTypePolyh];
const int maxGrEdge =
(numElem[ElemTypeTetra]*4 + numElem[ElemTypeHexa]*6 +
numElem[ElemTypePrism]*5 + numElem[ElemTypePyramid]*5 +
- numElem[ElemTypePolyh]*5)/2;
+ numElem[ElemTypeTrih]*3 + numElem[ElemTypePolyh]*5)/2;
graph.allocate(numGrVert, maxGrEdge);
// Make the graph
MakeGraphDIM<3>(model->firstRegion(), model->lastRegion(),
@@ -1204,6 +1208,7 @@ int CreatePartitionBoundaries(GModel *model, bool createGhostCells, bool createA
fillit_(faceToElement, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
fillit_(faceToElement, (*it)->prisms.begin(), (*it)->prisms.end());
fillit_(faceToElement, (*it)->pyramids.begin(), (*it)->pyramids.end());
+ fillit_(faceToElement, (*it)->trihedra.begin(), (*it)->trihedra.end());
fillit_(faceToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
}
std::multimap<MFace, MElement*, Less_Face>::iterator it = faceToElement.begin();
@@ -1235,6 +1240,7 @@ int CreatePartitionBoundaries(GModel *model, bool createGhostCells, bool createA
fillit_(edgeToElement, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
fillit_(edgeToElement, (*it)->prisms.begin(), (*it)->prisms.end());
fillit_(edgeToElement, (*it)->pyramids.begin(), (*it)->pyramids.end());
+ fillit_(edgeToElement, (*it)->trihedra.begin(), (*it)->trihedra.end());
fillit_(edgeToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
}
}
@@ -1268,6 +1274,7 @@ int CreatePartitionBoundaries(GModel *model, bool createGhostCells, bool createA
fillit_(vertexToElement, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
fillit_(vertexToElement, (*it)->prisms.begin(), (*it)->prisms.end());
fillit_(vertexToElement, (*it)->pyramids.begin(), (*it)->pyramids.end());
+ fillit_(vertexToElement, (*it)->trihedra.begin(), (*it)->trihedra.end());
fillit_(vertexToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
}
}
diff --git a/Mesh/meshPartitionOptions.h b/Mesh/meshPartitionOptions.h
index a3c96fa8771103c41e37fa2a72dc87fbf4910224..52c05574d1f10f00b1060bb589f144b6a6074616 100644
--- a/Mesh/meshPartitionOptions.h
+++ b/Mesh/meshPartitionOptions.h
@@ -78,6 +78,7 @@ class meshPartitionOptions
int tetWeight;
int priWeight;
int pyrWeight;
+ int trihWeight;
int hexWeight;
// NODAL WEIGHT
@@ -130,6 +131,7 @@ class meshPartitionOptions
tetWeight = 1;
priWeight = 1;
pyrWeight = 1;
+ trihWeight = 0;
hexWeight = 1;
}
};
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index fcef9661f769f368e5a2ed38092632a8e76fa6e9..7a22b59c944109b04d7ca3556a5b08c3ac3188e6 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -19,10 +19,15 @@
#include "MQuadrangle.h"
#include "MTriangle.h"
#include "MPyramid.h"
+#include "MTrihedron.h"
#include "MPrism.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
+#include "CondNumBasis.h"
+
+
+
void export_gregion_mesh(GRegion *gr, string filename)
{
// FIXME: use MElement::writeMSH
@@ -2297,6 +2302,111 @@ bool Recombinator::valid(Hex &hex,const std::set<MElement*>& parts){
}
}
+//A face with 4 nodes on the boundary should be close to planar
+//A face with 3 nodes on the boundary should not exist
+bool validFace(MVertex *a, MVertex *b, MVertex *c, MVertex *d){
+ //A face should not have exactly 3 nodes on the same surface
+ //or the face tag is undefined
+ //A better option would be to topologically check that the fase has not 4 nodes
+ //on the surface while the face is not no the surface (face2elems map)
+ MVertex *v[4] = {a, b, c, d};
+ std::map<GFace*, int> nbTagOccurences;
+ int maxOccurences = 0;
+ for (int iV = 0; iV < 4; iV++){
+ if (v[iV]->onWhat()->dim() < 3){
+ std::list<GFace*> faces;
+ if (v[iV]->onWhat()->dim() == 2)
+ faces.push_back((GFace*) v[iV]->onWhat());
+ else
+ faces = v[iV]->onWhat()->faces();
+ // printf("nbFaces %i\n",faces.size());
+ for (std::list<GFace*>::iterator it=faces.begin(); it != faces.end(); it++){
+ GFace* idFace = (*it);
+ // printf("Id is %p %i\n",idFace,idFace->getNativeInt());
+ if (nbTagOccurences.find(idFace)==nbTagOccurences.end())
+ nbTagOccurences[idFace] = 1;
+ else
+ nbTagOccurences[idFace]++;
+ maxOccurences = std::max(nbTagOccurences[idFace], maxOccurences);
+ }
+ }
+ }
+ if (maxOccurences == 3) return false;
+
+
+
+ //Geometry: A face with 4 nodes on the boundary should be close enough to planar
+ int nbBndNodes = 0;
+ if (a->onWhat()->dim() < 3) nbBndNodes++;
+ if (b->onWhat()->dim() < 3) nbBndNodes++;
+ if (c->onWhat()->dim() < 3) nbBndNodes++;
+ if (d->onWhat()->dim() < 3) nbBndNodes++;
+ if (nbBndNodes == 4){
+ SVector3 vec1 = SVector3(b->x()-a->x(),b->y()-a->y(),b->z()-a->z()).unit();
+ SVector3 vec2 = SVector3(c->x()-a->x(),c->y()-a->y(),c->z()-a->z()).unit();
+ SVector3 vec3 = SVector3(d->x()-a->x(),d->y()-a->y(),d->z()-a->z()).unit();
+
+ SVector3 crossVec1Vec2 = crossprod(vec1, vec2);
+ double angle = fabs(acos(dot(crossVec1Vec2, vec3))*180/M_PI);
+ double maxAngle = 15;
+ if (fabs(angle-90) > maxAngle) return false;
+ }
+ //A face with 3 nodes on the boundary should not exist
+ //Constraint should be relaxed
+ if (nbBndNodes == 3){
+ return false;
+ }
+ return true;
+}
+
+bool validFaces(Hex &hex){
+ bool v1, v2, v3, v4, v5, v6;
+ MVertex *a,*b,*c,*d;
+ MVertex *e,*f,*g,*h;
+
+ a = hex.get_a();
+ b = hex.get_b();
+ c = hex.get_c();
+ d = hex.get_d();
+ e = hex.get_e();
+ f = hex.get_f();
+ g = hex.get_g();
+ h = hex.get_h();
+
+ v1 = validFace(a,b,c,d); //SHOULD CHECK GEOMETRY
+ v2 = validFace(e,f,g,h);
+ v3 = validFace(a,b,f,e);
+ v4 = validFace(b,c,g,f);
+ v5 = validFace(d,c,g,h);
+ v6 = validFace(d,a,e,h);
+
+ return v1 && v2 && v3 && v4 && v5 && v6;
+}
+
+bool validCondNum(Hex &hex){
+ return true;
+ CondNumBasis basis(MSH_HEX_8);
+ fullMatrix<double> nodesXYZ(8,3);
+ for (int iV = 0; iV < 8; iV++){
+ nodesXYZ(iV,0) = hex.getVertex(iV)->x();
+ nodesXYZ(iV,1) = hex.getVertex(iV)->y();
+ nodesXYZ(iV,2) = hex.getVertex(iV)->z();
+ }
+ fullMatrix<double> normals(8,3);
+ fullVector<double> invCondNum(basis.getNumCondNumNodes());
+ basis.getSignedInvCondNum(nodesXYZ, normals, invCondNum);
+ for (int i=0; i<invCondNum.size(); i++){
+ if (invCondNum(i) * invCondNum(0) <= 0){
+ // nodesXYZ.print("HEX nodes");
+ // invCondNum.print("HEX invCondNum");
+ // Msg::Error("Wrong Hex");
+ return false;
+ }
+ }
+ return true;
+}
+
+
// renvoie true si le "MQuadrangle::etaShapeMeasure" des 6 faces est plus grand que 0.000001
bool Recombinator::valid(Hex &hex){
double k;
@@ -2324,10 +2434,10 @@ bool Recombinator::valid(Hex &hex){
eta6 = eta(d,c,g,h);
if(eta1>k && eta2>k && eta3>k && eta4>k && eta5>k && eta6>k){
- return 1;
+ return (validFaces(hex) && validCondNum(hex));
}
else{
- return 0;
+ return false;
}
}
@@ -3982,6 +4092,65 @@ bool Supplementary::valid(Prism prism,const std::set<MElement*>& parts){
}
}
+
+bool validFaces(Prism &prism){
+ bool v1, v2, v3;
+ MVertex *a,*b,*c;
+ MVertex *d, *e,*f;
+
+ a = prism.get_a();
+ b = prism.get_b();
+ c = prism.get_c();
+ d = prism.get_d();
+ e = prism.get_e();
+ f = prism.get_f();
+
+ v1 = validFace(a,d,f,c); //SHOULD CHECK GEOMETRY
+ v2 = validFace(a,d,e,b);
+ v3 = validFace(b,c,e,f);
+
+ return v1 && v2 && v3;
+}
+
+
+bool validCondNum(Prism &prism){
+ return true;
+ CondNumBasis basis(MSH_PRI_6);
+ fullMatrix<double> nodesXYZ(6,3);
+ nodesXYZ(0,0) = prism.get_a()->x();
+ nodesXYZ(0,1) = prism.get_a()->y();
+ nodesXYZ(0,2) = prism.get_a()->z();
+ nodesXYZ(1,0) = prism.get_b()->x();
+ nodesXYZ(1,1) = prism.get_b()->y();
+ nodesXYZ(1,2) = prism.get_b()->z();
+ nodesXYZ(2,0) = prism.get_c()->x();
+ nodesXYZ(2,1) = prism.get_c()->y();
+ nodesXYZ(2,2) = prism.get_c()->z();
+ nodesXYZ(3,0) = prism.get_d()->x();
+ nodesXYZ(3,1) = prism.get_d()->y();
+ nodesXYZ(3,2) = prism.get_d()->z();
+ nodesXYZ(4,0) = prism.get_e()->x();
+ nodesXYZ(4,1) = prism.get_e()->y();
+ nodesXYZ(4,2) = prism.get_e()->z();
+ nodesXYZ(5,0) = prism.get_f()->x();
+ nodesXYZ(5,1) = prism.get_f()->y();
+ nodesXYZ(5,2) = prism.get_f()->z();
+
+ fullMatrix<double> normals(6,3);
+ fullVector<double> invCondNum(basis.getNumCondNumNodes());
+ basis.getSignedInvCondNum(nodesXYZ, normals, invCondNum);
+ for (int i=0; i<invCondNum.size(); i++){
+ if (invCondNum(i) * invCondNum(0) <= 0){
+ // nodesXYZ.print("HEX nodes");
+ // invCondNum.print("HEX invCondNum");
+ // Msg::Error("Wrong PRISM");
+ return false;
+ }
+ }
+ return true;
+}
+
+
bool Supplementary::valid(Prism prism){
double k;
double eta1,eta2,eta3;
@@ -4002,7 +4171,7 @@ bool Supplementary::valid(Prism prism){
eta3 = eta(b,c,f,e);
if(eta1>k && eta2>k && eta3>k){
- return 1;
+ return (validFaces(prism) && validCondNum(prism));
}
else{
return 0;
@@ -4719,7 +4888,7 @@ PostOp::PostOp(){}
PostOp::~PostOp(){}
-void PostOp::execute(bool flag){
+void PostOp::execute(int level, bool insertTrihedra){
GRegion* gr;
GModel* model = GModel::current();
GModel::riter it;
@@ -4728,24 +4897,26 @@ void PostOp::execute(bool flag){
{
gr = *it;
if(gr->getNumMeshElements()>0){
- execute(gr,flag);
+ execute(gr,level, insertTrihedra);
}
}
}
-void PostOp::execute(GRegion* gr,bool flag){
+void PostOp::execute(GRegion* gr,int level,bool insertTrihedra){
printf("................PYRAMIDS................\n");
estimate1 = 0;
estimate2 = 0;
iterations = 0;
build_tuples(gr);
- init_markings(gr);
- build_vertex_to_tetrahedra(gr);
- pyramids1(gr);
- rearrange(gr);
+ if (level >= 1){
+ init_markings(gr);
+ build_vertex_to_tetrahedra(gr);
+ pyramids1(gr);
+ rearrange(gr);
+ }
- if(flag){
+ if(level >= 2){
init_markings(gr);
build_vertex_to_tetrahedra(gr);
build_vertex_to_pyramids(gr);
@@ -4753,6 +4924,15 @@ void PostOp::execute(GRegion* gr,bool flag){
rearrange(gr);
}
+ if (insertTrihedra){
+ init_markings(gr);
+ build_vertex_to_tetrahedra(gr);
+ build_vertex_to_pyramids(gr);
+ trihedra(gr);
+ rearrange(gr);
+ }
+
+
statistics(gr);
modify_surfaces(gr);
@@ -4807,11 +4987,11 @@ void PostOp::pyramids1(GRegion* gr){
g = element->getVertex(6);
h = element->getVertex(7);
- pyramids1(a,b,c,d,gr);
+ pyramids1(b,a,d,c,gr);
pyramids1(e,f,g,h,gr);
pyramids1(a,b,f,e,gr);
pyramids1(b,c,g,f,gr);
- pyramids1(d,c,g,h,gr);
+ pyramids1(c,d,h,g,gr);
pyramids1(d,a,e,h,gr);
}
@@ -4844,6 +5024,7 @@ void PostOp::pyramids1(GRegion* gr){
}
}
+
void PostOp::pyramids2(GRegion* gr){
unsigned int i;
MVertex *a,*b,*c,*d;
@@ -4870,7 +5051,7 @@ void PostOp::pyramids2(GRegion* gr){
for(i=0;i<hexahedra.size();i++){
element = hexahedra[i];
-
+
a = element->getVertex(0);
b = element->getVertex(1);
c = element->getVertex(2);
@@ -4933,36 +5114,53 @@ void PostOp::pyramids2(GRegion* gr){
void PostOp::pyramids1(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
MVertex* vertex;
std::set<MElement*> bin;
- std::set<MElement*> bin1;
- std::set<MElement*> bin2;
+ std::set<MElement*> bin1a;
+ std::set<MElement*> bin2a;
+ std::set<MElement*> bin1b;
+ std::set<MElement*> bin2b;
std::set<MElement*>::iterator it;
std::map<MElement*,bool>::iterator it1;
std::map<MElement*,bool>::iterator it2;
- bin1.clear();
- bin2.clear();
- find_tetrahedra(a,c,bin1);
- find_tetrahedra(b,d,bin2);
+ bin1a.clear();
+ bin1b.clear();
+ bin2a.clear();
+ bin2b.clear();
+ find_tetrahedra(a,c,d,bin1a);
+ find_tetrahedra(a,b,c,bin1b);
+ find_tetrahedra(b,d,a,bin2a);
+ find_tetrahedra(b,c,d,bin2b);
+
+// bin1.clear();
+// bin2.clear();
+// find_tetrahedra(a,c,bin1);
+// find_tetrahedra(b,d,bin2);
bin.clear();
- for(it=bin1.begin();it!=bin1.end();it++){
+ for(it=bin1a.begin();it!=bin1a.end();it++){
bin.insert(*it);
}
- for(it=bin2.begin();it!=bin2.end();it++){
+ for(it=bin1b.begin();it!=bin1b.end();it++){
+ bin.insert(*it);
+ }
+ for(it=bin2a.begin();it!=bin2a.end();it++){
+ bin.insert(*it);
+ }
+ for(it=bin2b.begin();it!=bin2b.end();it++){
bin.insert(*it);
}
- if(bin.size()==2){
- it = bin.begin();
- it1 = markings.find(*it);
+ if(bin.size()==2){ //2 tetrahedra on face
+ it = bin.begin();
+ it1 = markings.find(*it); //1st tetrahedra
it++;
- it2 = markings.find(*it);
+ it2 = markings.find(*it); //2nd tetrahedra
if(it1->second==0 && it2->second==0){
vertex = find(a,b,c,d,*it);
-
- if(vertex!=0){
- gr->addPyramid(new MPyramid(a,b,c,d,vertex));
+ MVertex* vertex1 = find(a,b,c,d,*bin.begin());
+ if(vertex!=0 && vertex == vertex1){//Seb added this 2nd condition. Is it necessary??
+ gr->addPyramid(new MPyramid(a,d,c,b,vertex));
it1->second = 1;
it2->second = 1;
}
@@ -4970,6 +5168,93 @@ void PostOp::pyramids1(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
}
}
+
+void PostOp::trihedra(GRegion* gr){
+ unsigned int i;
+ MVertex *a,*b,*c,*d;
+ MVertex *e,*f,*g,*h;
+ MElement* element;
+ std::vector<MElement*> hexahedra;
+ std::vector<MElement*> prisms;
+ std::vector<MTetrahedron*> opt1;
+ std::vector<MPyramid*> opt2;
+ std::map<MElement*,bool>::iterator it;
+
+ hexahedra.clear();
+ prisms.clear();
+
+ for(i=0;i<gr->getNumMeshElements();i++){
+ element = gr->getMeshElement(i);
+ if(eight(element)){
+ hexahedra.push_back(element);
+ }
+ else if(six(element)){
+ prisms.push_back(element);
+ }
+ }
+
+ for(i=0;i<hexahedra.size();i++){
+ element = hexahedra[i];
+
+ a = element->getVertex(0);
+ b = element->getVertex(1);
+ c = element->getVertex(2);
+ d = element->getVertex(3);
+ e = element->getVertex(4);
+ f = element->getVertex(5);
+ g = element->getVertex(6);
+ h = element->getVertex(7);
+
+ trihedra(a,b,c,d,gr);
+ trihedra(e,f,g,h,gr);
+ trihedra(a,b,f,e,gr);
+ trihedra(b,c,g,f,gr);
+ trihedra(d,c,g,h,gr);
+ trihedra(d,a,e,h,gr);
+ }
+
+ for(i=0;i<prisms.size();i++){
+ element = prisms[i];
+
+ a = element->getVertex(0);
+ b = element->getVertex(1);
+ c = element->getVertex(2);
+ d = element->getVertex(3);
+ e = element->getVertex(4);
+ f = element->getVertex(5);
+
+ trihedra(a,d,f,c,gr);
+ trihedra(a,b,e,d,gr);
+ trihedra(b,c,f,e,gr);
+ }
+
+}
+
+void PostOp::trihedra(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
+ std::set<MElement*> diag1a;
+ std::set<MElement*> diag1b;
+ std::set<MElement*> diag2a;
+ std::set<MElement*> diag2b;
+
+ find_tetrahedra(a,b,c,diag1a);
+ find_tetrahedra(a,c,d,diag1b);
+ find_tetrahedra(b,c,d,diag2a);
+ find_tetrahedra(a,b,d,diag2b);
+ find_pyramids_from_tri(a,b,c,diag1a);
+ find_pyramids_from_tri(a,c,d,diag1b);
+ find_pyramids_from_tri(b,c,d,diag2a);
+ find_pyramids_from_tri(a,b,d,diag2b);
+ if(diag1a.size()==1 && diag1b.size()==1){
+ MTrihedron *trih = new MTrihedron(b, c, d, a);
+ // if (diag1a.size() != 1 || diag1b.size() != 1) Msg::Error("Quad face neighbor with %i+%i triangular faces (other diagonal: %i+%i) Trihedron: %i",diag1a.size(),diag1b.size(),diag2a.size(),diag2b.size(),trih->getNum());
+ gr->addTrihedron(trih);
+ }else if(diag2a.size()==1 && diag2b.size()==1){
+ MTrihedron *trih = new MTrihedron(a, b, c, d);
+ // if (diag2a.size() != 1 || diag2b.size() != 1) Msg::Error("Quad face neighbor with %i+%i triangular faces (other diagonal: %i+%i) Trihedron: %i",diag2a.size(),diag2b.size(),diag1a.size(),diag1b.size(),trih->getNum());
+ gr->addTrihedron(trih);
+ }
+}
+
void PostOp::pyramids2(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
bool flag;
double x,y,z;
@@ -5055,6 +5340,7 @@ void PostOp::pyramids2(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
gr->addMeshVertex(mid);
temp2 = new MPyramid(a,b,c,d,mid);
+
gr->addPyramid(temp2);
markings.insert(std::pair<MElement*,bool>(temp2,false));
build_vertex_to_pyramids(temp2);
@@ -5173,6 +5459,7 @@ void PostOp::statistics(GRegion* gr){
nbr6 = 0;
nbr5 = 0;
nbr4 = 0;
+ nbr4Trih = 0;
vol = 0.0;
vol8 = 0.0;
vol6 = 0.0;
@@ -5202,6 +5489,10 @@ void PostOp::statistics(GRegion* gr){
vol4 = vol4 + element->getVolume();
}
+ if(fourTrih(element)){
+ nbr4Trih = nbr4Trih + 1;
+ }
+
nbr = nbr + 1;
if(!five(element)){
@@ -5218,6 +5509,7 @@ void PostOp::statistics(GRegion* gr){
printf(" percentage of prisms : %.2f\n",nbr6*100.0/nbr);
printf(" percentage of pyramids : %.2f\n",nbr5*100.0/nbr);
printf(" percentage of tetrahedra : %.2f\n",nbr4*100.0/nbr);
+ printf(" percentage of trihedra : %.2f\n",nbr4Trih*100.0/nbr);
printf("Volume :\n");
printf(" percentage of hexahedra : %.2f\n",vol8*100.0/vol);
printf(" percentage of prisms : %.2f\n",vol6*100.0/vol);
@@ -5412,7 +5704,12 @@ void PostOp::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
}
bool PostOp::four(MElement* element){
- if(element->getNumVertices()==4) return 1;
+ if(element->getNumVertices()==4 && element->getNumFaces()==4) return 1;
+ else return 0;
+}
+
+bool PostOp::fourTrih(MElement* element){
+ if(element->getNumVertices()==4 && element->getNumFaces()==3) return 1;
else return 0;
}
@@ -5610,6 +5907,72 @@ void PostOp::find_tetrahedra(MVertex* v1,MVertex* v2,std::set<MElement*>& final)
}
}
+void PostOp::find_tetrahedra(MVertex* v1,MVertex* v2,MVertex* v3,std::set<MElement*>& final){
+ std::map<MVertex*,std::set<MElement*> >::iterator it1;
+ std::map<MVertex*,std::set<MElement*> >::iterator it2;
+ std::map<MVertex*,std::set<MElement*> >::iterator it3;
+
+ it1 = vertex_to_tetrahedra.find(v1);
+ it2 = vertex_to_tetrahedra.find(v2);
+ it3 = vertex_to_tetrahedra.find(v3);
+
+ std::set<MElement*> buf;
+
+ if(it1!=vertex_to_tetrahedra.end() && it2!=vertex_to_tetrahedra.end() && it3!=vertex_to_tetrahedra.end()){
+ intersection(it1->second,it2->second,buf);
+ intersection(buf,it3->second,final);
+
+ }
+}
+
+void PostOp::find_pyramids_from_tri(MVertex* v1,MVertex* v2,MVertex* v3,std::set<MElement*>& final){
+ bool flag;
+ std::set<MElement*>::iterator it;
+ std::map<MVertex*,std::set<MElement*> >::iterator it1;
+ std::map<MVertex*,std::set<MElement*> >::iterator it2;
+ std::map<MVertex*,std::set<MElement*> >::iterator it3;
+ std::set<MElement*> temp1, temp2;
+
+ it1 = vertex_to_pyramids.find(v1);
+ it2 = vertex_to_pyramids.find(v2);
+ it3 = vertex_to_pyramids.find(v3);
+
+ temp1.clear();
+ temp2.clear();
+
+ if(it1!=vertex_to_pyramids.end() && it2!=vertex_to_pyramids.end() && it3!=vertex_to_pyramids.end()){
+ intersection(it1->second,it2->second,temp1);
+ intersection(temp1,it3->second,temp2);
+ }
+ //Now temp2 contains pyramids containing the 3 vertices
+ //We check that one vertex is on the summit and the others are neighbors in the base
+ for(it=temp2.begin();it!=temp2.end();it++){
+ flag = false;
+ if (v1 == (*it)->getVertex(4)){
+ flag = (equal(v2,v3,(*it)->getVertex(0),(*it)->getVertex(1)) ||
+ equal(v2,v3,(*it)->getVertex(1),(*it)->getVertex(2)) ||
+ equal(v2,v3,(*it)->getVertex(2),(*it)->getVertex(3)) ||
+ equal(v2,v3,(*it)->getVertex(3),(*it)->getVertex(0)));
+ }
+ if (v2 == (*it)->getVertex(4)){
+ flag = (equal(v1,v3,(*it)->getVertex(0),(*it)->getVertex(1)) ||
+ equal(v1,v3,(*it)->getVertex(1),(*it)->getVertex(2)) ||
+ equal(v1,v3,(*it)->getVertex(2),(*it)->getVertex(3)) ||
+ equal(v1,v3,(*it)->getVertex(3),(*it)->getVertex(0)));
+ }
+ if (v3 == (*it)->getVertex(4)){
+ flag = (equal(v1,v2,(*it)->getVertex(0),(*it)->getVertex(1)) ||
+ equal(v1,v2,(*it)->getVertex(1),(*it)->getVertex(2)) ||
+ equal(v1,v2,(*it)->getVertex(2),(*it)->getVertex(3)) ||
+ equal(v1,v2,(*it)->getVertex(3),(*it)->getVertex(0)));
+ }
+ if(flag){
+ final.insert(*it);
+ }
+ }
+}
+
+
void PostOp::find_pyramids(MVertex* v1,MVertex* v2,std::set<MElement*>& final){
bool flag1,flag2,flag3,flag4;
bool flag5,flag6,flag7,flag8;
@@ -5642,6 +6005,39 @@ void PostOp::find_pyramids(MVertex* v1,MVertex* v2,std::set<MElement*>& final){
}
}
+//void PostOp::find_pyramids(MVertex* v1,MVertex* v2,std::set<MElement*>& final){
+// bool flag1,flag2,flag3,flag4;
+// bool flag5,flag6,flag7,flag8;
+// std::set<MElement*>::iterator it;
+// std::map<MVertex*,std::set<MElement*> >::iterator it1;
+// std::map<MVertex*,std::set<MElement*> >::iterator it2;
+// std::set<MElement*> temp;
+//
+// it1 = vertex_to_pyramids.find(v1);
+// it2 = vertex_to_pyramids.find(v2);
+//
+// temp.clear();
+//
+// if(it1!=vertex_to_pyramids.end() && it2!=vertex_to_pyramids.end()){
+// intersection(it1->second,it2->second,temp);
+// }
+//
+// for(it=temp.begin();it!=temp.end();it++){
+// flag1 = equal(v1,v2,(*it)->getVertex(0),(*it)->getVertex(1));
+// flag2 = equal(v1,v2,(*it)->getVertex(1),(*it)->getVertex(2));
+// flag3 = equal(v1,v2,(*it)->getVertex(2),(*it)->getVertex(3));
+// flag4 = equal(v1,v2,(*it)->getVertex(3),(*it)->getVertex(0));
+// flag5 = equal(v1,v2,(*it)->getVertex(0),(*it)->getVertex(4));
+// flag6 = equal(v1,v2,(*it)->getVertex(1),(*it)->getVertex(4));
+// flag7 = equal(v1,v2,(*it)->getVertex(2),(*it)->getVertex(4));
+// flag8 = equal(v1,v2,(*it)->getVertex(3),(*it)->getVertex(4));
+// if(flag1 || flag2 || flag3 || flag4 || flag5 || flag6 || flag7 || flag8){
+// final.insert(*it);
+// }
+// }
+//}
+
+
void PostOp::intersection(const std::set<MElement*>& bin1,const std::set<MElement*>& bin2,std::set<MElement*>& final){
std::set_intersection(bin1.begin(),bin1.end(),bin2.begin(),bin2.end(),std::inserter(final,final.end()));
}
@@ -6221,7 +6617,7 @@ Recombinator_Graph::~Recombinator_Graph(){
void Recombinator_Graph::createBlossomInfo(){
- throw;
+ throw;
GRegion* gr;
@@ -6606,7 +7002,8 @@ void Recombinator_Graph::merge_clique(GRegion* gr, cliques_losses_graph<Hex*> &c
// inserting the hex
quality = quality + current_hex->get_quality();
- gr->addHexahedron(new MHexahedron(current_hex->get_a(), current_hex->get_b(),current_hex->get_c(),current_hex->get_d(),current_hex->get_e(),current_hex->get_f(),current_hex->get_g(),current_hex->get_h()));
+ MHexahedron *h=new MHexahedron(current_hex->get_a(), current_hex->get_b(),current_hex->get_c(),current_hex->get_d(),current_hex->get_e(),current_hex->get_f(),current_hex->get_g(),current_hex->get_h());
+ gr->addHexahedron(h);
count++;
// hex_to_export.insert(current_hex);
// cout << " inserting " << current_hex << " made of ";
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index d3e005603dec6ee34bd040abbf76f2a33e7d3d23..ec7b9b07ae33ba572cd5209c131eea315785a291 100644
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -609,7 +609,7 @@ public:
class PostOp{
private:
- int nbr,nbr8,nbr6,nbr5,nbr4;
+ int nbr,nbr8,nbr6,nbr5,nbr4,nbr4Trih;
double vol,vol8,vol6,vol5,vol4;
int estimate1;
int estimate2;
@@ -623,8 +623,9 @@ public:
PostOp();
~PostOp();
- void execute(bool);
- void execute(GRegion*,bool);
+ void execute(int, bool);
+ //0: no pyramid, 1: single-step, 2: two-steps (conforming), true: fill non-conformities with trihedra
+ void execute(GRegion*,int level, bool addTrihedra);
inline int get_nb_hexahedra()const{return nbr8;};
inline double get_vol_hexahedra()const{return vol8;};
@@ -634,8 +635,10 @@ public:
void init_markings(GRegion*);
void pyramids1(GRegion*);
void pyramids2(GRegion*);
+ void trihedra(GRegion*);
void pyramids1(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
void pyramids2(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
+ void trihedra(MVertex*,MVertex*,MVertex*,MVertex*,GRegion*);
void rearrange(GRegion*);
void statistics(GRegion*);
void build_tuples(GRegion*);
@@ -643,6 +646,7 @@ public:
void modify_surfaces(MVertex*,MVertex*,MVertex*,MVertex*);
bool four(MElement*);
+ bool fourTrih(MElement*);
bool five(MElement*);
bool six(MElement*);
bool eight(MElement*);
@@ -656,6 +660,8 @@ public:
MVertex* find(MVertex*,MVertex*,MVertex*,MVertex*,MElement*);
void find_tetrahedra(MVertex*,MVertex*,std::set<MElement*>&);
+ void find_tetrahedra(MVertex*,MVertex*,MVertex*,std::set<MElement*>&);
+ void find_pyramids_from_tri(MVertex*,MVertex*,MVertex*,std::set<MElement*>&);
void find_pyramids(MVertex*,MVertex*,std::set<MElement*>&);
void intersection(const std::set<MElement*>&,const std::set<MElement*>&,std::set<MElement*>&);
diff --git a/Numeric/CondNumBasis.cpp b/Numeric/CondNumBasis.cpp
index cf46578ae835c1f96609b4cb66e8db53e113ece2..06a62aeb856f88946968faff6475f9839f359c37 100644
--- a/Numeric/CondNumBasis.cpp
+++ b/Numeric/CondNumBasis.cpp
@@ -315,10 +315,16 @@ CondNumBasis::CondNumBasis(int tag, int cnOrder) :
_tag(tag), _dim(ElementType::DimensionFromTag(tag)),
_condNumOrder(cnOrder >= 0 ? cnOrder : condNumOrder(tag))
{
+ if ( ElementType::ParentTypeFromTag(tag) == TYPE_TRIH){
+ _nCondNumNodes = 1;
+ _nMapNodes = 4;
+ _nPrimMapNodes = 4;
+ return;
+ }
const bool serendip = false;
- FuncSpaceData data = ( ElementType::ParentTypeFromTag(tag) == TYPE_PYR )? FuncSpaceData(true, tag, true, _condNumOrder+2, _condNumOrder, &serendip) : FuncSpaceData(true, tag, _condNumOrder, &serendip);
- _gradBasis = BasisFactory::getGradientBasis(data);
+ FuncSpaceData data = ( ElementType::ParentTypeFromTag(tag) == TYPE_PYR )? FuncSpaceData(true, tag, true, _condNumOrder+1, _condNumOrder, &serendip) : FuncSpaceData(true, tag, _condNumOrder, &serendip);
+
fullMatrix<double> lagPoints; // Sampling points
gmshGeneratePoints(data, lagPoints);
_nCondNumNodes = lagPoints.size1();
@@ -381,6 +387,7 @@ int CondNumBasis::condNumOrder(int parentType, int order)
case TYPE_PRI : return order;
case TYPE_HEX : return order;
case TYPE_PYR : return (order == 1) ? 0 : order;
+ case TYPE_TRIH : return 0;
default :
Msg::Error("Unknown element type %d, return order 0", parentType);
return 0;
@@ -398,7 +405,7 @@ inline void CondNumBasis::getInvCondNumGeneral(int nCondNumNodes,
const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals,
fullVector<double> &condNum) const
-{
+{
switch (_dim) {
case 0 : {
@@ -426,6 +433,10 @@ inline void CondNumBasis::getInvCondNumGeneral(int nCondNumNodes,
}
case 3 : {
+ if (ElementType::ParentTypeFromTag(_tag) == TYPE_TRIH){
+ for (int i = 0; i < nCondNumNodes; i++) condNum(i) = 1.;
+ break;
+ }
fullMatrix<double> dxyzdX(nCondNumNodes, 3), dxyzdY(nCondNumNodes, 3), dxyzdZ(nCondNumNodes, 3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
@@ -519,6 +530,17 @@ inline void CondNumBasis::getInvCondNumAndGradientsGeneral(int nCondNumNodes,
}
case 3 : {
+ if (ElementType::ParentTypeFromTag(_tag) == TYPE_TRIH){
+ for (int i = 0; i < nCondNumNodes; i++) {
+ for (int j = 0; j < _nMapNodes; j++) {
+ IDI(i,j) = 0.;
+ IDI(i,j+1*_nMapNodes) = 0.;
+ IDI(i,j+2*_nMapNodes) = 0.;
+ }
+ IDI(i,3*_nMapNodes) = 1.;
+ }
+ break;
+ }
fullMatrix<double> dxyzdX(nCondNumNodes,3), dxyzdY(nCondNumNodes,3), dxyzdZ(nCondNumNodes,3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
diff --git a/Numeric/ElementType.cpp b/Numeric/ElementType.cpp
index 4cf8a204774617dbe5ecda6fb1d46e6d32773b11..6da8698049deeb417b4548e7e337247d21d19c86 100644
--- a/Numeric/ElementType.cpp
+++ b/Numeric/ElementType.cpp
@@ -93,6 +93,9 @@ int ElementType::ParentTypeFromTag(int tag)
case(MSH_TRI_MINI):
case(MSH_TET_MINI):
return TYPE_MINI;
+ case(MSH_TRIH_4):
+ return TYPE_TRIH;
+
default:
Msg::Error("Unknown element tag %i for parent type, returning -1.", tag);
return -1;
@@ -230,6 +233,7 @@ int ElementType::OrderFromTag(int tag)
case MSH_PYR_53 : return 7;
case MSH_PYR_61 : return 8;
case MSH_PYR_69 : return 9;
+ case MSH_TRIH_4 : return 1;
case MSH_POLYH_ : return 1;
default :
Msg::Warning("Unknown element tag %d, assuming order 1.",tag);
@@ -320,6 +324,8 @@ int ElementType::DimensionFromTag(int tag)
case(MSH_HEX_68): case(MSH_HEX_80):
case(MSH_HEX_92): case(MSH_HEX_104):
+ case(MSH_TRIH_4):
+
case(MSH_POLYH_):
return 3;
@@ -393,6 +399,8 @@ int ElementType::SerendipityFromTag(int tag)
case MSH_PYR_140 : case MSH_PYR_204 :
case MSH_PYR_285 : case MSH_PYR_385 :
+ case MSH_TRIH_4 :
+
return 0; // Not Serendipity
@@ -549,6 +557,8 @@ int ElementType::getTag(int parentTag, int order, bool serendip)
default : Msg::Error("pyramid order %i unknown", order); return 0;
}
break;
+ case TYPE_TRIH :
+ return MSH_TRIH_4;
default : Msg::Warning("unknown element type %i, returning 0", parentTag); return 0;
}
}
diff --git a/Post/PViewData.h b/Post/PViewData.h
index fff10746d1df20541006dd7d088419c52029bdc2..b6112e30d5241b02169d293b3f52117e21e015ce 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -113,6 +113,7 @@ class PViewData {
virtual int getNumHexahedra(int step=-1){ return 0; }
virtual int getNumPrisms(int step=-1){ return 0; }
virtual int getNumPyramids(int step=-1){ return 0; }
+ virtual int getNumTrihedra(int step=-1){ return 0; }
virtual int getNumPolyhedra(int step=-1){ return 0; }
// return the number of geometrical entities in the view
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 3b00fc528b30660bd27cacd173390fec3129183e..5fccf53ebf8eedbe3395fe75c401b06d2aa410f3 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -449,6 +449,16 @@ int PViewDataGModel::getNumPyramids(int step)
return n;
}
+int PViewDataGModel::getNumTrihedra(int step)
+{
+ if(_steps.empty()) return 0;
+ GModel *m = _steps[0]->getModel(); // to generalize
+ int n = 0;
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
+ n += (*it)->trihedra.size();
+ return n;
+}
+
int PViewDataGModel::getNumPolyhedra(int step)
{
if(_steps.empty()) return 0;
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 3e4e8dfab55eb6ab1892ca6456681849f294c68c..6ccd0fc8ebb69edfe47d898cdaff59d0c425a97c 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -202,6 +202,7 @@ class PViewDataGModel : public PViewData {
int getNumHexahedra(int step=-1);
int getNumPrisms(int step=-1);
int getNumPyramids(int step=-1);
+ int getNumTrihedra(int step=-1);
int getNumPolyhedra(int step=-1);
int getNumEntities(int step=-1);
int getNumElements(int step=-1, int ent=-1);
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index fe39970349785d9e39e518f3203732c8177614a6..e971ebeb390d87c99c663959ae9d8903be6f20e0 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -20,12 +20,13 @@ PViewDataList::PViewDataList(bool isAdapted)
NbSG(0), NbVG(0), NbTG(0),
NbSS(0), NbVS(0), NbTS(0), NbSH(0), NbVH(0), NbTH(0),
NbSI(0), NbVI(0), NbTI(0), NbSY(0), NbVY(0), NbTY(0),
- NbSD(0), NbVD(0), NbTD(0), NbT2(0), NbT3(0),
+ NbSR(0), NbVR(0), NbTR(0), NbSD(0), NbVD(0), NbTD(0),
+ NbT2(0), NbT3(0),
_lastElement(-1), _lastDimension(-1), _lastNumNodes(-1),
_lastNumComponents(-1), _lastNumValues(-1), _lastNumEdges(-1),
_lastType(-1), _lastXYZ(0), _lastVal(0), _isAdapted(isAdapted)
{
- for(int i = 0; i < 30; i++) _index[i] = 0;
+ for(int i = 0; i < 33; i++) _index[i] = 0;
polyTotNumNodes[0] = 0.; polyTotNumNodes[1] = 0.;
polyAgNumNodes[0].push_back(0.); polyAgNumNodes[1].push_back(0.);
}
@@ -74,6 +75,8 @@ bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolatio
_stat(TI, 9, NbTI, 6, TYPE_PRI);
_stat(SY, 1, NbSY, 5, TYPE_PYR); _stat(VY, 3, NbVY, 5, TYPE_PYR);
_stat(TY, 9, NbTY, 5, TYPE_PYR);
+ _stat(SY, 1, NbSR, 4, TYPE_TRIH); _stat(VY, 3, NbVR, 4, TYPE_TRIH);
+ _stat(TY, 9, NbTR, 4, TYPE_TRIH);
_stat(SG, 1, NbSG, 3, TYPE_POLYG); _stat(VG, 3, NbVG, 3, TYPE_POLYG);
_stat(TG, 9, NbTG, 3, TYPE_POLYG);
_stat(SD, 1, NbSD, 4, TYPE_POLYH); _stat(VD, 3, NbVD, 4, TYPE_POLYH);
@@ -87,11 +90,11 @@ bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolatio
}
// compute starting element indices
- int nb[30] = {NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT,
+ int nb[33] = {NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT,
NbSQ, NbVQ, NbTQ, NbSS, NbVS, NbTS, NbSH, NbVH, NbTH,
- NbSI, NbVI, NbTI, NbSY, NbVY, NbTY, NbSG, NbVG, NbTG,
- NbSD, NbVD, NbTD};
- for(int i = 0; i < 30; i++){
+ NbSI, NbVI, NbTI, NbSY, NbVY, NbTY, NbSR, NbVR, NbTR,
+ NbSG, NbVG, NbTG, NbSD, NbVD, NbTD};
+ for(int i = 0; i < 33; i++){
_index[i] = 0;
for(int j = 0; j <= i; j++)
_index[i] += nb[j];
@@ -104,17 +107,17 @@ bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolatio
int PViewDataList::getNumScalars(int step)
{
- return NbSP + NbSL + NbST + NbSQ + NbSS + NbSH + NbSI + NbSY + NbSG + NbSD;
+ return NbSP + NbSL + NbST + NbSQ + NbSS + NbSH + NbSI + NbSY + NbSR + NbSG + NbSD;
}
int PViewDataList::getNumVectors(int step)
{
- return NbVP + NbVL + NbVT + NbVQ + NbVS + NbVH + NbVI + NbVY + NbVG + NbVD;
+ return NbVP + NbVL + NbVT + NbVQ + NbVS + NbVH + NbVI + NbVY + NbVR + NbVG + NbVD;
}
int PViewDataList::getNumTensors(int step)
{
- return NbTP + NbTL + NbTT + NbTQ + NbTS + NbTH + NbTI + NbTY + NbTG + NbTD;
+ return NbTP + NbTL + NbTT + NbTQ + NbTS + NbTH + NbTI + NbTY + NbTR + NbTG + NbTD;
}
int PViewDataList::getNumElements(int step, int ent)
@@ -334,17 +337,22 @@ void PViewDataList::_setLast(int ele)
else if(ele < _index[22]) _setLast(ele - _index[21], 3, 5, 3, 8, TYPE_PYR, VY, NbVY);
else _setLast(ele - _index[22], 3, 5, 9, 8, TYPE_PYR, TY, NbTY);
}
- else if(ele < _index[26]){ // polygons
- int nN = polyNumNodes[0][ele - _index[23]];
- if(ele < _index[24]) _setLast(ele - _index[23], 2, nN, 1, nN, TYPE_POLYG, SG, NbSG);
- else if(ele < _index[25]) _setLast(ele - _index[24], 2, nN, 3, nN, TYPE_POLYG, VG, NbVG);
- else _setLast(ele - _index[25], 2, nN, 9, nN, TYPE_POLYG, TG, NbTG);
+ else if(ele < _index[26]){ // trihedra
+ if(ele < _index[24]) _setLast(ele - _index[23], 3, 4, 1, 5, TYPE_TRIH, SR, NbSR);
+ else if(ele < _index[25]) _setLast(ele - _index[24], 3, 4, 3, 5, TYPE_TRIH, VR, NbVR);
+ else _setLast(ele - _index[25], 3, 4, 9, 5, TYPE_TRIH, TR, NbTR);
}
- else if(ele < _index[29]){ // polyhedra
- int nN = polyNumNodes[1][ele - _index[26]];
- if(ele < _index[27]) _setLast(ele - _index[26], 3, nN, 1, nN*1.5, TYPE_POLYH, SD, NbSD);
- else if(ele < _index[28]) _setLast(ele - _index[27], 3, nN, 3, nN*1.5, TYPE_POLYH, VD, NbVD);
- else _setLast(ele - _index[28], 3, nN, 9, nN*1.5, TYPE_POLYH, TD, NbTD);
+ else if(ele < _index[29]){ // polygons
+ int nN = polyNumNodes[0][ele - _index[26]];
+ if(ele < _index[27]) _setLast(ele - _index[26], 2, nN, 1, nN, TYPE_POLYG, SG, NbSG);
+ else if(ele < _index[28]) _setLast(ele - _index[27], 2, nN, 3, nN, TYPE_POLYG, VG, NbVG);
+ else _setLast(ele - _index[28], 2, nN, 9, nN, TYPE_POLYG, TG, NbTG);
+ }
+ else if(ele < _index[32]){ // polyhedra
+ int nN = polyNumNodes[1][ele - _index[29]];
+ if(ele < _index[30]) _setLast(ele - _index[29], 3, nN, 1, nN*1.5, TYPE_POLYH, SD, NbSD);
+ else if(ele < _index[32]) _setLast(ele - _index[30], 3, nN, 3, nN*1.5, TYPE_POLYH, VD, NbVD);
+ else _setLast(ele - _index[31], 3, nN, 9, nN*1.5, TYPE_POLYH, TD, NbTD);
}
}
@@ -576,12 +584,12 @@ void PViewDataList::smooth()
std::vector<double> *list = 0;
int *nbe = 0, nbc, nbn;
- for(int i = 0; i < 24; i++){
+ for(int i = 0; i < 27; i++){
_getRawData(i, &list, &nbe, &nbc, &nbn);
if(nbn > 1)
generateConnectivities(*list, *nbe, NbTimeStep, nbn, nbc, data);
}
- for(int i = 0; i < 24; i++){
+ for(int i = 0; i < 27; i++){
_getRawData(i, &list, &nbe, &nbc, &nbn);
if(nbn > 1)
smoothList(*list, *nbe, NbTimeStep, nbn, nbc, data);
@@ -603,6 +611,7 @@ double PViewDataList::getMemoryInMb()
b += (SH.size() + VH.size() + TH.size()) * sizeof(double);
b += (SI.size() + VI.size() + TI.size()) * sizeof(double);
b += (SY.size() + VY.size() + TY.size()) * sizeof(double);
+ b += (SR.size() + VR.size() + TR.size()) * sizeof(double);
b += (SD.size() + VD.size() + TD.size()) * sizeof(double);
b += (T2D.size() + T3D.size()) * sizeof(double);
return b / 1024. / 1024.;
@@ -656,7 +665,8 @@ bool PViewDataList::combineSpace(nameData &nd)
dVecMerge(l->SI, SI); NbSI += l->NbSI; dVecMerge(l->VI, VI); NbVI += l->NbVI;
dVecMerge(l->TI, TI); NbTI += l->NbTI; dVecMerge(l->SY, SY); NbSY += l->NbSY;
dVecMerge(l->VY, VY); NbVY += l->NbVY; dVecMerge(l->TY, TY); NbTY += l->NbTY;
-
+ dVecMerge(l->VR, VR); NbVR += l->NbVR; dVecMerge(l->TR, TR); NbTR += l->NbTR;
+
// merge strings
for(unsigned int i = 0; i < l->T2D.size(); i += 4){
T2D.push_back(l->T2D[i]);
@@ -725,7 +735,7 @@ bool PViewDataList::combineTime(nameData &nd)
std::vector<double> *list = 0, *list2 = 0;
// use the first data set as the reference
- for(int i = 0; i < 24; i++){
+ for(int i = 0; i < 27; i++){
_getRawData(i, &list, &nbe, &nbc, &nbn);
data[0]->_getRawData(i, &list2, &nbe2, &nbc2, &nbn2);
*nbe = *nbe2;
@@ -739,7 +749,7 @@ bool PViewDataList::combineTime(nameData &nd)
_interpolation[it->first].push_back(new fullMatrix<double>(*it->second[i]));
// merge values for all element types
- for(int i = 0; i < 24; i++){
+ for(int i = 0; i < 27; i++){
_getRawData(i, &list, &nbe, &nbc, &nbn);
for(int j = 0; j < *nbe; j++){
for(unsigned int k = 0; k < data.size(); k++){
@@ -853,7 +863,7 @@ int PViewDataList::_getRawData(int idxtype, std::vector<double> **l, int **ne,
{
int type = 0;
// No constant nn for polygons!
- if(idxtype > 23 && idxtype < 30)
+ if(idxtype > 26 && idxtype < 33)
Msg::Warning("No constant number of nodes for polygons and polyhedra");
switch(idxtype){
case 0 : *l = &SP; *ne = &NbSP; *nc = 1; *nn = 1; type = TYPE_PNT; break;
@@ -880,12 +890,15 @@ int PViewDataList::_getRawData(int idxtype, std::vector<double> **l, int **ne,
case 21: *l = &SY; *ne = &NbSY; *nc = 1; *nn = 5; type = TYPE_PYR; break;
case 22: *l = &VY; *ne = &NbVY; *nc = 3; *nn = 5; type = TYPE_PYR; break;
case 23: *l = &TY; *ne = &NbTY; *nc = 9; *nn = 5; type = TYPE_PYR; break;
- case 24: *l = &SG; *ne = &NbSG; *nc = 1; *nn = 3; type = TYPE_POLYG; break;
- case 25: *l = &VG; *ne = &NbVG; *nc = 3; *nn = 3; type = TYPE_POLYG; break;
- case 26: *l = &TG; *ne = &NbTG; *nc = 9; *nn = 3; type = TYPE_POLYG; break;
- case 27: *l = &SD; *ne = &NbSD; *nc = 1; *nn = 4; type = TYPE_POLYH; break;
- case 28: *l = &VD; *ne = &NbVD; *nc = 3; *nn = 4; type = TYPE_POLYH; break;
- case 29: *l = &TD; *ne = &NbTD; *nc = 9; *nn = 4; type = TYPE_POLYH; break;
+ case 24: *l = &SR; *ne = &NbSR; *nc = 1; *nn = 4; type = TYPE_TRIH; break;
+ case 25: *l = &VR; *ne = &NbVR; *nc = 3; *nn = 4; type = TYPE_TRIH; break;
+ case 26: *l = &TR; *ne = &NbTR; *nc = 9; *nn = 4; type = TYPE_TRIH; break;
+ case 27: *l = &SG; *ne = &NbSG; *nc = 1; *nn = 3; type = TYPE_POLYG; break;
+ case 28: *l = &VG; *ne = &NbVG; *nc = 3; *nn = 3; type = TYPE_POLYG; break;
+ case 29: *l = &TG; *ne = &NbTG; *nc = 9; *nn = 3; type = TYPE_POLYG; break;
+ case 30: *l = &SD; *ne = &NbSD; *nc = 1; *nn = 4; type = TYPE_POLYH; break;
+ case 31: *l = &VD; *ne = &NbVD; *nc = 3; *nn = 4; type = TYPE_POLYH; break;
+ case 32: *l = &TD; *ne = &NbTD; *nc = 9; *nn = 4; type = TYPE_POLYH; break;
default: Msg::Error("Wrong type in PViewDataList"); break;
}
@@ -963,6 +976,11 @@ std::vector<double> *PViewDataList::incrementList(int numComp, int type, int num
else if(numComp == 3){ NbVY++; return &VY; }
else if(numComp == 9){ NbTY++; return &TY; }
break;
+ case TYPE_TRIH:
+ if (numComp == 1){ NbSR++; return &SR; }
+ else if(numComp == 3){ NbVR++; return &VR; }
+ else if(numComp == 9){ NbTR++; return &TR; }
+ break;
case TYPE_POLYG:
polyNumNodes[0].push_back(numNodes);
nb = (polyAgNumNodes[0].size()) ? polyAgNumNodes[0].back() : 0;
diff --git a/Post/PViewDataList.h b/Post/PViewDataList.h
index 97c1f74881c71d44b0f7615c0d191388331e71ba..59d340f3165f3561e0440a3a0ec371d1db885ee5 100644
--- a/Post/PViewDataList.h
+++ b/Post/PViewDataList.h
@@ -41,6 +41,8 @@ class PViewDataList : public PViewData {
std::vector<double> SI, VI, TI; // prisms
int NbSY, NbVY, NbTY;
std::vector<double> SY, VY, TY; // pyramids
+ int NbSR, NbVR, NbTR;
+ std::vector<double> SR, VR, TR; // trihedra
int NbSD, NbVD, NbTD;
std::vector<double> SD, VD, TD; // polyhedra
int NbT2, NbT3;
@@ -50,7 +52,7 @@ class PViewDataList : public PViewData {
std::vector<int> polyAgNumNodes[2];
int polyTotNumNodes[2];
private:
- int _index[30];
+ int _index[33];
int _lastElement, _lastDimension;
int _lastNumNodes, _lastNumComponents, _lastNumValues, _lastNumEdges, _lastType;
double *_lastXYZ, *_lastVal;
@@ -90,6 +92,7 @@ class PViewDataList : public PViewData {
int getNumHexahedra(int step=-1){ return NbSH + NbVH + NbTH; }
int getNumPrisms(int step=-1){ return NbSI + NbVI + NbTI; }
int getNumPyramids(int step=-1){ return NbSY + NbVY + NbTY; }
+ int getNumTrihedra(int step=-1){ return NbSR + NbVR + NbTR; }
int getNumPolyhedra(int step=-1){ return NbSD + NbVD + NbTD; }
int getNumEntities(int step=-1){ return 1; }
int getNumElements(int step=-1, int ent=-1);
diff --git a/Post/PViewOptions.cpp b/Post/PViewOptions.cpp
index 6e602a96255e60d778a9ca45eae55a11727d582b..241df2ee59eb26e68948d2038d2aa1a32e27fcba 100644
--- a/Post/PViewOptions.cpp
+++ b/Post/PViewOptions.cpp
@@ -137,6 +137,7 @@ bool PViewOptions::skipElement(int type)
case TYPE_HEX: return !drawHexahedra;
case TYPE_PRI: return !drawPrisms;
case TYPE_PYR: return !drawPyramids;
+ case TYPE_TRIH: return !drawTrihedra;
case TYPE_POLYH: return false;
default: return true;
}
diff --git a/Post/PViewOptions.h b/Post/PViewOptions.h
index a5ebaa5224415e578fc513d914815c5543c1f653..99c8c25aa6612aa78249f365449c977f67956d5e 100644
--- a/Post/PViewOptions.h
+++ b/Post/PViewOptions.h
@@ -82,7 +82,7 @@ class PViewOptions {
double currentTime;
int drawStrings;
int drawPoints, drawLines, drawTriangles, drawQuadrangles, drawPolygons;
- int drawTetrahedra, drawHexahedra, drawPrisms, drawPyramids, drawPolyhedra;
+ int drawTetrahedra, drawHexahedra, drawPrisms, drawPyramids, drawTrihedra, drawPolyhedra;
int drawScalars, drawVectors, drawTensors;
int boundary, pointType, lineType, drawSkinOnly;
double pointSize, lineWidth;
@@ -103,7 +103,7 @@ class PViewOptions {
int closed;
struct{
unsigned int point, line, triangle, quadrangle;
- unsigned int tetrahedron, hexahedron, prism, pyramid;
+ unsigned int tetrahedron, hexahedron, prism, pyramid, trihedron;
unsigned int tangents, normals;
unsigned int text2d, text3d, axes, background2d;
} color;
diff --git a/Post/PViewVertexArrays.cpp b/Post/PViewVertexArrays.cpp
index 8f8970075c68e51df25afff66731fe6acae001ff..e3df2253e539908baa88562a7bd4a5061cc3138a 100644
--- a/Post/PViewVertexArrays.cpp
+++ b/Post/PViewVertexArrays.cpp
@@ -857,6 +857,20 @@ static void addScalarPyramid(PView *p, double **xyz,
addScalarTetrahedron(p, xyz, val, pre, is[i][0], is[i][1], is[i][2], is[i][3]);
}
+static void addOutlineTrihedron(PView *p, double **xyz,
+ unsigned int color, bool pre)
+{
+ addOutlineQuadrangle(p, xyz, color, pre, 0, 1, 2, 3);
+}
+
+static void addScalarTrihedron(PView *p, double **xyz,
+ double **val, bool pre, int i0=0,
+ int i1=1, int i2=2, int i3=3, bool unique=false)
+{
+ addScalarQuadrangle(p, xyz, val, pre, i0, i1, i2, i3, unique);
+}
+
+
static void addOutlinePolyhedron(PView *p, double **xyz,
unsigned int color, bool pre, int numNodes)
{
@@ -915,6 +929,7 @@ static void addOutlineElement(PView *p, int type, double **xyz,
case TYPE_HEX: addOutlineHexahedron(p, xyz, opt->color.hexahedron, pre); break;
case TYPE_PRI: addOutlinePrism(p, xyz, opt->color.prism, pre); break;
case TYPE_PYR: addOutlinePyramid(p, xyz, opt->color.pyramid, pre); break;
+ case TYPE_TRIH: addOutlineTrihedron(p, xyz, opt->color.pyramid, pre); break;
case TYPE_POLYH: addOutlinePolyhedron(p, xyz, opt->color.pyramid, pre, numNodes); break;
}
}
@@ -932,6 +947,7 @@ static void addScalarElement(PView *p, int type, double **xyz,
case TYPE_HEX: addScalarHexahedron(p, xyz, val, pre); break;
case TYPE_PRI: addScalarPrism(p, xyz, val, pre); break;
case TYPE_PYR: addScalarPyramid(p, xyz, val, pre); break;
+ case TYPE_TRIH: addScalarTrihedron(p, xyz, val, pre); break;
case TYPE_POLYH: addScalarPolyhedron(p, xyz, val, pre, numNodes); break;
}
}
@@ -1342,6 +1358,7 @@ class initPView {
int tets = data->getNumTetrahedra(opt->timeStep);
int prisms = data->getNumPrisms(opt->timeStep);
int pyrs = data->getNumPyramids(opt->timeStep);
+ int trihs = data->getNumTrihedra(opt->timeStep);
int hexas = data->getNumHexahedra(opt->timeStep);
int polyhs = data->getNumPolyhedra(opt->timeStep);
int heuristic = 0;
@@ -1349,10 +1366,10 @@ class initPView {
heuristic = (tets + prisms + pyrs + hexas + polyhs) / 10;
else if(opt->intervalsType == PViewOptions::Continuous)
heuristic = (tris + 2 * quads + 3 * polygs + 6 * tets +
- 8 * prisms + 6 * pyrs + 12 * hexas + 10 * polyhs);
+ 8 * prisms + 6 * pyrs + 2 * trihs + 12 * hexas + 10 * polyhs);
else if(opt->intervalsType == PViewOptions::Discrete)
heuristic = (tris + 2 * quads + 3 * polygs + 6 * tets +
- 8 * prisms + 6 * pyrs + 12 * hexas + 10 * polyhs) * 2;
+ 8 * prisms + 6 * pyrs + 2 * trihs + 12 * hexas + 10 * polyhs) * 2;
heuristic = _estimateIfClipped(p, heuristic);
return heuristic + 10000;
}
diff --git a/benchmarks/3d/Submarine/Submarine.geo b/benchmarks/3d/Submarine/Submarine.geo
index ec10af692e6dd02ccd235a9f077390f73873829d..e592c80a88f7f5ac337495486b786a565c9ffe3c 100644
--- a/benchmarks/3d/Submarine/Submarine.geo
+++ b/benchmarks/3d/Submarine/Submarine.geo
@@ -6,32 +6,37 @@ Mesh.Algorithm3D=9;
Mesh.Smoothing=0;
Mesh.Recombine3DAll=1;
Mesh.SaveParametric = 1;
-MeshAlgorithm Surface {12} 6;
+MeshAlgorithm Surface {12} 5; //Not in final mesh
MeshAlgorithm Surface {10} 5;
MeshAlgorithm Surface {11} 5;
//Mesh.Algorithm=5;
//Mesh.Algorithm3D=5;
-Mesh.Optimize = 1;
+Mesh.Optimize = 0;
//Mesh.CharacteristicLengthMax = 25;
-Characteristic Length {14, 13, 15, 16, 10, 9, 11, 12, 6, 5, 19, 8, 4, 3, 7, 2, 1} = 100;
-Characteristic Length {21} = 200;
-Characteristic Length {17, 18} = 100;
-Characteristic Length {20} = 200;
+LcMin = 39;
+LcMax = 78;
+//LcMin = 80;
+//LcMax = 160;
+Characteristic Length {14, 13, 15, 16, 10, 9, 11, 12, 6, 5, 19, 8, 4, 3, 7, 2, 1} = LcMax;
+Characteristic Length {21} = LcMin;
+Characteristic Length {17, 18} = LcMin;
+Characteristic Length {20} = LcMax;
Compound Line(34) = {22, 24, 25};
Compound Surface(13) = {10};
Compound Surface(14) = {11};
+Compound Surface(15) = {12};
submarineLine[] = {15, 21, 17, 26, 11, 20, 13, 27, 29, 8, 19, 10, 34};
Physical Line(1) = submarineLine[];
submarineSurf[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 13};
Physical Surface(2) = submarineSurf[];
-Physical Surface(3) = {12,14};
+Physical Surface(3) = {14,15};
Physical Volume(4) = {1};
Physical Line(5) = {30, 31};
diff --git a/benchmarks/3d/Submarine/optimize.py b/benchmarks/3d/Submarine/optimize.py
index 3093216b562489ee8b59da5356090e195673cb01..929a5c10b645a43d44a8f63de296ec22d9a17d1e 100644
--- a/benchmarks/3d/Submarine/optimize.py
+++ b/benchmarks/3d/Submarine/optimize.py
@@ -1,26 +1,38 @@
from gmshpy import *
import sys
import pickle
-
-name = "Submarine2"
+import preprofix
+name = "Submarine"
g = GModel()
g.load(name + ".geo")
-#g.mesh(3)
-#g.save(name + ".msh")
-g.load(name + ".msh")
+g.mesh(3)
+g.writeMSH(name + ".msh", 3.0)
+#g.load(name + ".msh")
+g.setAllVolumesPositiveTopology()
+
+g.writeMSH(name + "_topology.msh", 3.0)
+preprofix.fixTetras(g)
+preprofix.fixPyramids(g)
+preprofix.fixPrisms(g)
+preprofix.fixHexs(g)
+g.writeMSH(name + "_FIX.msh", 3.0)
+#exit()
+#preprofix.fixPyramids(g)
+
+#g.writeMSH(name + "Fix.msh",3.0)
OH = MeshQualOptParameters()
OH.onlyVisible = False
OH.dim = 3
-OH.fixBndNodes = True # Fix boundary nodes or not
-OH.strategy = 1 # 0 = Connected blobs, 1 = Adaptive one-by-one (recommended in 3D)
+OH.fixBndNodes = False # Fix boundary nodes or not
+OH.strategy = 0 # 0 = Connected blobs, 1 = Adaptive one-by-one (recommended in 3D)
OH.excludeHex = False
OH.excludePrism = False
-OH.nbLayers = 2 # Nb. of layers around invalid element to construct blob
-OH.distanceFactor = 2 # Distance factor to construct blob
+OH.nbLayers = 3 # Nb. of layers around invalid element to construct blob
+OH.distanceFactor = 3 # Distance factor to construct blob
OH.maxPatchAdapt = 5 # Number of blob adaption iteration
OH.maxLayersAdaptFact = 3 # Factor to multiply nb. of layers when adapting
OH.distanceAdaptFact = 3 # Factor to multiply distance factor when adapting
@@ -36,7 +48,7 @@ OH.nCurses = 1
OH.logFileName = "log"
OH.maxOptIter = 20 # Nb of optimixation iterations
-OH.maxBarrierUpdates = 30 # Nb. of optimization passes
+OH.maxBarrierUpdates = 15 # Nb. of optimization passes
#print("minTargetIdealJac = %g" % OH.minTargetIdealJac)
print("minTargetInvCondNum = %g" % OH.minTargetInvCondNum)
@@ -55,6 +67,11 @@ MeshQualityOptimizer(g, OH)
print("RESULT: minInvCondNum = %g, maxInvCondNum = %g, CPU = %g" %
(OH.minInvCondNum, OH.maxInvCondNum, OH.CPU))
-g.save(name + "_opt.msh")
+g.writeMSH(name + "_opt.msh", 3.0)
+pOpt = meshPartitionOptions()
+pOpt.setNumOfPartitions(5)
+PartitionMesh(g, pOpt)
+g.writeMSH(name + "_opt_part.msh", 3.0)
+
#g.writeMSH(name + "_OK.msh",3.0,False,False,False,1.0,0,0)
diff --git a/benchmarks/3d/Submarine/preprofix.py b/benchmarks/3d/Submarine/preprofix.py
index 51264e1b7c12a9b598187ae50217a2c3cc48ba90..51207d5ba5c327cad01e48f74ff7479f617bafd8 100644
--- a/benchmarks/3d/Submarine/preprofix.py
+++ b/benchmarks/3d/Submarine/preprofix.py
@@ -1,48 +1,60 @@
import gmshpy
-import sys
-m = gmshpy.GModel()
-if (len(sys.argv) < 3) :
- print("please specify input and output files")
- exit()
-m.load(sys.argv[1])
+import math
-# generate a list of boundary elements
-
-def isBoundaryFace(f) :
- vent = [f.getVertex(i).onWhat() for i in range(f.getNumVertices())]
- return all([e.dim() < 3 for e in vent])
def areCoplanarFaces(f0, f1, tol) :
n0 = f0.normal()
n1 = f1.normal()
- return abs(gmshpy.dot(n0, n1)) > tol
+ return (abs(math.acos(gmshpy.dot(n0, n1))*180/math.pi) < tol)
+
+def nbVerticesOnSurface(el) :
+ vent = [el.getVertex(i).onWhat().dim() for i in range(el.getNumVertices())]
+ return len([e for e in vent if e < 3])
-def isBadPyramid(p) :
- fs = [p.getFace(i) for i in range(p.getNumFaces())]
- bndfs = [f for f in fs if (isBoundaryFace(f) and f.getNumVertices() == 3)]
+def isBadElement(el) :
+ fs = [el.getFace(i) for i in range(el.getNumFaces())]
bad = False
- for i in range(len(bndfs)) :
+ for i in range(len(fs)) :
for j in range(i):
- bad |= areCoplanarFaces(bndfs[0], bndfs[1], 0.9)
- return bad
-
-
-for r in m.bindingsGetRegions() :
- badp = [p for p in r.pyramids if isBadPyramid(p)]
- for p in badp :
- cog = p.barycenter()
- nv = gmshpy.MVertex(cog.x(), cog.y(), cog.z(), r);
- nv.thisown = False
- r.addMeshVertex(nv)
- for i in range(4) :
- f = p.getFace(i)
- vs = [f.getVertex(i) for i in range(3)] + [nv]
- t = gmshpy.MTetrahedron(*vs)
- t.thisown = False
+ bad |= areCoplanarFaces(fs[i], fs[j], 10) #10 degrees min angle allowed
+ return bad
+
+def fixTetras(m):
+ for r in m.bindingsGetRegions() :
+ goodt = [t for t in r.tetrahedra if (not ((nbVerticesOnSurface(t) == 4) and isBadElement(t)))]
+ oldSize = r.tetrahedra.size()
+ r.tetrahedra.clear()
+ for t in goodt :
r.addTetrahedron(t)
- p.setVertex(4, nv)
-
-m.save(sys.argv[2])
+ print("-- Removed %i tetrahedra over %i"%(oldSize-r.tetrahedra.size(), oldSize))
+
+
+def fixPyramids(m):
+ for r in m.bindingsGetRegions() :
+ goodp = [p for p in r.pyramids if (not ((nbVerticesOnSurface(p) == 5) and isBadElement(p)))]
+ oldSize = r.pyramids.size()
+ r.pyramids.clear()
+ for p in goodp :
+ r.addPyramid(p)
+ print("-- Removed %i pyramids over %i"%(oldSize-r.pyramids.size(), oldSize))
+
+def fixPrisms(m):
+ for r in m.bindingsGetRegions() :
+ goodp = [p for p in r.prisms if (not ((nbVerticesOnSurface(p) == 6) and isBadElement(p)))]
+ oldSize = r.prisms.size()
+ r.prisms.clear()
+ for p in goodp :
+ r.addPrism(p)
+ print("-- Removed %i prisms over %i"%(oldSize-r.prisms.size(), oldSize))
+
+def fixHexs(m):
+ for r in m.bindingsGetRegions() :
+ goodh = [h for h in r.hexahedra if (not ((nbVerticesOnSurface(h) == 8) and isBadElement(p)))]
+ oldSize = r.hexahedra.size()
+ r.hexahedra.clear()
+ for h in goodh :
+ r.addHexahedron(h)
+ print("-- Removed %i hexahedra over %i"%(oldSize-r.hexahedra.size(), oldSize))
diff --git a/demos/title.script b/demos/title.script
index 7e49d3964ab60266cefbd6418a10125216a09a2f..c2ae942435f9aef0dcb03fc91d72baf042e6d299 100644
--- a/demos/title.script
+++ b/demos/title.script
@@ -29,7 +29,7 @@ Plugin(Annotate).Run ;
// if we have a mesh, print some statistics
nbelm = Mesh.NbTriangles + Mesh.NbQuadrangles + Mesh.NbTetrahedra +
- Mesh.NbHexahedra + Mesh.NbPrisms + Mesh.NbPyramids;
+ Mesh.NbHexahedra + Mesh.NbPrisms + Mesh.NbPyramids + Mesh.NbTrihedra;
nbnod = Mesh.NbNodes;
If(nbelm && nbnod)
y += 20;
diff --git a/wrappers/gmshpy/gmshGeo.i b/wrappers/gmshpy/gmshGeo.i
index 500a3d00992bd9408e603f604761f718b7daded7..d195c204bf5a99ef5a1ff26ff2c77405bf632713 100644
--- a/wrappers/gmshpy/gmshGeo.i
+++ b/wrappers/gmshpy/gmshGeo.i
@@ -31,6 +31,7 @@
#include "MTetrahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+ #include "MTrihedron.h"
#include "MHexahedron.h"
#include "MQuadrangle.h"
#include "MLine.h"
@@ -58,6 +59,9 @@ namespace std {
%template(MTetrahedronVector) vector< MTetrahedron *,std::allocator< MTetrahedron * > >;
%template(MQuadrangleVector) vector< MQuadrangle *,std::allocator< MQuadrangle * > >;
%template(MPyramidVector) vector< MPyramid *,std::allocator< MPyramid * > >;
+ %template(MTrihedronVector) vector< MTrihedron *,std::allocator< MTrihedron * > >;
+ %template(MHexahedronVector) vector< MHexahedron *,std::allocator< MHexahedron * > >;
+ %template(MPrismVector) vector< MPrism *,std::allocator< MPrism * > >;
%template(GEdgeVectorVector) vector< std::vector< GEdge *,std::allocator< GEdge * > >,std::allocator< std::vector< GEdge *,std::allocator< GEdge * > > > >;
%template(GFaceVectorVector) vector< std::vector< GFace *,std::allocator< GFace * > >,std::allocator< std::vector< GFace *,std::allocator< GFace * > > > >;
%template(GFaceList) list<GFace*, std::allocator<GFace*> >;
@@ -125,6 +129,7 @@ namespace std {
%include "MTetrahedron.h"
%include "MPrism.h"
%include "MPyramid.h"
+%include "MTrihedron.h"
%include "MHexahedron.h"
%include "MQuadrangle.h"
%include "MLine.h"
@@ -195,6 +200,11 @@ namespace std {
$self->xyz2uvw(xyz, &uvw[0]);
return uvw;
}
+ std::vector<MVertex*> getVertices() {
+ std::vector <MVertex*> vertices;
+ $self->getVertices(vertices);
+ return vertices;
+ }
}
%extend GEdge {
diff --git a/wrappers/gmshpy/gmshNumeric.i b/wrappers/gmshpy/gmshNumeric.i
index 3f3df01d41a0d83362a07cd4dba7186973e781eb..e58e8f7795193f8689736a9cf67fbb9278eda1ff 100644
--- a/wrappers/gmshpy/gmshNumeric.i
+++ b/wrappers/gmshpy/gmshNumeric.i
@@ -15,6 +15,7 @@
#include "fullMatrix.h"
#include "nodalBasis.h"
#include "BasisFactory.h"
+ #include "CondNumBasis.h"
#include "simpleFunction.h"
#include "simpleFunctionPython.h"
#include "polynomialBasis.h"
@@ -37,6 +38,7 @@
%include "polynomialBasis.h"
%include "pyramidalBasis.h"
%include "BasisFactory.h"
+%include "CondNumBasis.h"
%extend nodalBasis {
fullMatrix<double> F(const fullMatrix<double> &xi) {
fullMatrix<double> psi;