diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 46313931e14be03b38ba1712687f3a81c1383013..c55eacd260275a58e16de5969204c3dc9ebaa3f7 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -537,6 +537,10 @@ void Get_Options(int argc, char *argv[])
CTX.mesh.dual = 1;
i++;
}
+ else if(!strcmp(argv[i] + 1, "voronoi")) {
+ CTX.mesh.voronoi = 1;
+ i++;
+ }
else if(!strcmp(argv[i] + 1, "noview")) {
opt_view_visible(0, GMSH_SET, 0);
i++;
diff --git a/Common/Context.h b/Common/Context.h
index 81a26ba0b78d8211da215f21f936df3772f65e50..3e4161a615c5dfdec903055efee20f8743d32a8d 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -160,7 +160,7 @@ class Context_T {
double quality_inf, quality_sup, radius_inf, radius_sup;
double scaling_factor, lc_factor, rand_factor, lc_integration_precision, lc_min, lc_max;
int lc_from_points, lc_from_curvature, lc_extend_from_boundary;
- int dual, draw_skin_only;
+ int dual, voronoi, draw_skin_only;
int light, light_two_side, light_lines;
int format, nb_smoothing, algo2d, algo3d, algo_recombine;
int order, second_order_linear, second_order_incomplete;
@@ -202,7 +202,7 @@ class Context_T {
int jpeg_quality, jpeg_smoothing;
int gif_dither, gif_sort, gif_interlace, gif_transparent;
int geo_labels;
- int pos_elementary, pos_element, pos_gamma, pos_eta, pos_rho;
+ int pos_elementary, pos_element, pos_gamma, pos_eta, pos_rho,pos_disto;
int text, tex_as_equation;
} print;
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index 043e0253a9264c45c5c9d4ecd864c6f683114281..f65c4fa0e37657b8cd9d689ec9ff98b8938a0df4 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -185,7 +185,7 @@ void CreateOutputFile(const char *filename, int format)
case FORMAT_POS:
GModel::current()->writePOS(name, CTX.print.pos_elementary, CTX.print.pos_element,
CTX.print.pos_gamma, CTX.print.pos_eta, CTX.print.pos_rho,
- CTX.mesh.save_all, CTX.mesh.scaling_factor);
+ CTX.print.pos_disto,CTX.mesh.save_all, CTX.mesh.scaling_factor);
break;
case FORMAT_GEO:
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 90401d760ffb0cab5adf70574d3f981042cb8bad..9f854006e7203b3ea4757f6365be6a818cdd7b62 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1102,6 +1102,9 @@ StringXNumber MeshOptions_Number[] = {
"Display faces of volume mesh?" },
{ F|O, "VolumeNumbers" , opt_mesh_volumes_num , 0. ,
"Display volume mesh element numbers?" },
+ { F|O, "Voronoi" , opt_mesh_voronoi , 0. ,
+ "Display the voronoi diagram" },
+
{ F|O, "ZoneDefinition" , opt_mesh_zone_definition , 2. ,
"Method for defining a zone (0=single zone, 1=by partition, 2=by physical)" },
@@ -1474,6 +1477,8 @@ StringXNumber PrintOptions_Number[] = {
"Save Eta quality measure in mesh statistics exported as post-processing views" },
{ F|O, "PostRho" , opt_print_pos_rho , 0. ,
"Save Rho quality measure in mesh statistics exported as post-processing views" },
+ { F|O, "PostDisto" , opt_print_pos_disto , 0. ,
+ "Save Disto quality measure in mesh statistics exported as post-processing views" },
{ F|O, "TexAsEquation" , opt_print_tex_as_equation , 0. ,
"Print all TeX strings as equations" },
diff --git a/Common/GmshMatrix.h b/Common/GmshMatrix.h
index 49c0ec3e1c15d09f1e7838f4c055dbafc5caac60..f5c03b178b0b1c5246680b619c6439eabe163a3c 100644
--- a/Common/GmshMatrix.h
+++ b/Common/GmshMatrix.h
@@ -293,6 +293,11 @@ class GSL_Matrix
{
gsl_blas_dgemv(CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
}
+ inline void blas_dgemm(const GSL_Matrix & x, const GSL_Matrix& b,
+ const double c_a = 1.0, const double c_b = 1.0)
+ {
+ gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, c_a, x.data, b.data, c_b, data);
+ }
inline gsl_matrix_view touchSubmatrix(int i0, int ni, int j0, int nj)
{
return gsl_matrix_submatrix(data, i0, j0, ni, nj);
diff --git a/Common/Options.cpp b/Common/Options.cpp
index bf984c53fcec5675be3a7466e9943b7f89c4a6b2..f70550c8c6906640470b0a6c83785c19c9af0eb5 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -4546,7 +4546,7 @@ double opt_mesh_quality_type(OPT_ARGS_NUM)
if(CTX.mesh.quality_type != val)
CTX.mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
CTX.mesh.quality_type = (int)val;
- if(CTX.mesh.quality_type < 0 || CTX.mesh.quality_type > 2)
+ if(CTX.mesh.quality_type < 0 || CTX.mesh.quality_type > 3)
CTX.mesh.quality_type = 0;
}
#if defined(HAVE_FLTK)
@@ -5198,6 +5198,14 @@ double opt_mesh_dual(OPT_ARGS_NUM)
return CTX.mesh.dual;
}
+double opt_mesh_voronoi(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET) {
+ CTX.mesh.voronoi= (int)val;
+ }
+ return CTX.mesh.voronoi;
+}
+
double opt_mesh_draw_skin_only(OPT_ARGS_NUM)
{
if(action & GMSH_SET) {
@@ -7878,6 +7886,13 @@ double opt_print_pos_rho(OPT_ARGS_NUM)
return CTX.print.pos_rho;
}
+double opt_print_pos_disto(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.print.pos_disto = (int)val;
+ return CTX.print.pos_disto;
+}
+
double opt_print_gif_dither(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 8116a499641c51d0b56df04cc0a9c53a6f1810c7..d46301269a01616dfe2702230bf402b02d319643 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -486,6 +486,7 @@ double opt_mesh_second_order_experimental(OPT_ARGS_NUM);
double opt_mesh_second_order_linear(OPT_ARGS_NUM);
double opt_mesh_second_order_incomplete(OPT_ARGS_NUM);
double opt_mesh_dual(OPT_ARGS_NUM);
+double opt_mesh_voronoi(OPT_ARGS_NUM);
double opt_mesh_draw_skin_only(OPT_ARGS_NUM);
double opt_mesh_save_all(OPT_ARGS_NUM);
double opt_mesh_save_groups_of_nodes(OPT_ARGS_NUM);
@@ -675,6 +676,7 @@ double opt_print_pos_element(OPT_ARGS_NUM);
double opt_print_pos_gamma(OPT_ARGS_NUM);
double opt_print_pos_eta(OPT_ARGS_NUM);
double opt_print_pos_rho(OPT_ARGS_NUM);
+double opt_print_pos_disto(OPT_ARGS_NUM);
double opt_print_gif_dither(OPT_ARGS_NUM);
double opt_print_gif_sort(OPT_ARGS_NUM);
double opt_print_gif_interlace(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 48bb2c9b03a91c1bad89cc02f591dbaea2a1553d..e44a72d74298bd31bda100bdc5fea88eb706ee1c 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1844,8 +1844,10 @@ void statistics_histogram_cb(CALLBACK_ARGS)
type = 0;
else if(!strcmp(name, "Eta"))
type = 1;
- else
+ else if(!strcmp(name, "Rho"))
type = 2;
+ else
+ type = 3;
std::vector<double> x, y;
for(int i = 0; i < 100; i++) y.push_back(WID->quality[type][i]);
new PView(name, "# Elements", x, y);
@@ -3681,6 +3683,7 @@ void mesh_inspect_cb(CALLBACK_ARGS)
Msg::Direct(" Rho: %g", elements[0]->rhoShapeMeasure());
Msg::Direct(" Gamma: %g", elements[0]->gammaShapeMeasure());
Msg::Direct(" Eta: %g", elements[0]->etaShapeMeasure());
+ Msg::Direct(" Disto: %g", elements[0]->distoShapeMeasure());
CTX.mesh.changed = ENT_ALL;
Draw();
WID->create_message_window();
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 86083c9a7be7f02db4a91ce3ae727db875e6abc8..6caa0418e3cfb07b572833c7e1a491d48b02d7e8 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -2617,6 +2617,7 @@ void GUI::create_option_window()
{"Gamma", 0, 0, 0},
{"Eta", 0, 0, 0},
{"Rho", 0, 0, 0},
+ {"Disto", 0, 0, 0},
{0}
};
mesh_choice[6] = new Fl_Choice(L + 2 * WB + IW / 2, 2 * WB + 8 * BH, IW/2, BH, "Quality range");
@@ -4008,7 +4009,7 @@ void GUI::create_statistics_window()
}
int width = 26 * fontsize;
- int height = 5 * WB + 17 * BH;
+ int height = 5 * WB + 18 * BH;
stat_window = new Dialog_Window(width, height, CTX.non_modal_windows, "Statistics");
stat_window->box(GMSH_WINDOW_BOX);
@@ -4045,6 +4046,8 @@ void GUI::create_statistics_window()
stat_value[num]->tooltip("~ volume^(2/3) / sum_edge_length^2"); num++;
stat_value[num] = new Fl_Output(2 * WB, 2 * WB + 15 * BH, IW, BH, "Rho");
stat_value[num]->tooltip("~ min_edge_length / max_edge_length"); num++;
+ stat_value[num] = new Fl_Output(2 * WB, 2 * WB + 16 * BH, IW, BH, "Disto");
+ stat_value[num]->tooltip("~ min (J0/J, J/J0)"); num++;
stat_butt[0] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 13 * BH, BB, BH, "Graph");
stat_butt[0]->callback(statistics_histogram_cb, (void *)"Gamma");
@@ -4052,6 +4055,8 @@ void GUI::create_statistics_window()
stat_butt[1]->callback(statistics_histogram_cb, (void *)"Eta");
stat_butt[2] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 15 * BH, BB, BH, "Graph");
stat_butt[2]->callback(statistics_histogram_cb, (void *)"Rho");
+ stat_butt[3] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 16 * BH, BB, BH, "Graph");
+ stat_butt[3]->callback(statistics_histogram_cb, (void *)"Disto");
g[1]->end();
}
@@ -4124,7 +4129,10 @@ void GUI::set_statistics(bool compute_quality)
sprintf(label[num], "%g", s[15]); stat_value[num]->value(label[num]); num++;
if(!compute_quality){
- for(int i = 0; i < 3; i++) stat_butt[i]->deactivate();
+ for(int i = 0; i < 4; i++) stat_butt[i]->deactivate();
+ sprintf(label[num], "Press Update");
+ stat_value[num]->deactivate();
+ stat_value[num]->value(label[num]); num++;
sprintf(label[num], "Press Update");
stat_value[num]->deactivate();
stat_value[num]->value(label[num]); num++;
@@ -4136,7 +4144,7 @@ void GUI::set_statistics(bool compute_quality)
stat_value[num]->value(label[num]); num++;
}
else{
- for(int i = 0; i < 3; i++) stat_butt[i]->activate();
+ for(int i = 0; i < 4; i++) stat_butt[i]->activate();
sprintf(label[num], "%.4g (%.4g->%.4g)", s[17], s[18], s[19]);
stat_value[num]->activate();
stat_value[num]->value(label[num]); num++;
@@ -4146,6 +4154,9 @@ void GUI::set_statistics(bool compute_quality)
sprintf(label[num], "%.4g (%.4g->%.4g)", s[23], s[24], s[25]);
stat_value[num]->activate();
stat_value[num]->value(label[num]); num++;
+ sprintf(label[num], "%.4g (%.4g->%.4g)", s[46], s[47], s[48]);
+ stat_value[num]->activate();
+ stat_value[num]->value(label[num]); num++;
}
// post
diff --git a/Fltk/GUI.h b/Fltk/GUI.h
index 48a5336919ca0ed754a9a8325773843237d5f6c7..7f44a933798b6b66a603007a277632bcf92341db 100644
--- a/Fltk/GUI.h
+++ b/Fltk/GUI.h
@@ -229,8 +229,8 @@ public:
// statistics window
Fl_Window *stat_window;
Fl_Output *stat_value[50];
- Fl_Button *stat_butt[3];
- double quality[3][100];
+ Fl_Button *stat_butt[4];
+ double quality[4][100];
// message window
Fl_Window *msg_window;
diff --git a/Fltk/GUI_Extras.cpp b/Fltk/GUI_Extras.cpp
index 896c0625f056b896b599dae9a45611a1d440e161..3c034ca49405503e4b2c94d5ecc0e6b8cefc3ed2 100644
--- a/Fltk/GUI_Extras.cpp
+++ b/Fltk/GUI_Extras.cpp
@@ -763,7 +763,7 @@ int pos_dialog(const char *name)
{
struct _pos_dialog{
Fl_Window *window;
- Fl_Check_Button *b[6];
+ Fl_Check_Button *b[7];
Fl_Button *ok, *cancel;
};
static _pos_dialog *dialog = NULL;
@@ -774,7 +774,7 @@ int pos_dialog(const char *name)
if(!dialog){
dialog = new _pos_dialog;
- int h = 3 * WB + 7 * BH, w = 2 * BB + 3 * WB, y = WB;
+ int h = 3 * WB + 8 * BH, w = 2 * BB + 3 * WB, y = WB;
// not a "Dialog_Window" since it is modal
dialog->window = new Fl_Double_Window(w, h, "POS Options");
dialog->window->box(GMSH_WINDOW_BOX);
@@ -784,7 +784,8 @@ int pos_dialog(const char *name)
dialog->b[3] = new Fl_Check_Button(WB, y, 2 * BB + WB, BH, "Print Gamma quality measure"); y += BH;
dialog->b[4] = new Fl_Check_Button(WB, y, 2 * BB + WB, BH, "Print Eta quality measure"); y += BH;
dialog->b[5] = new Fl_Check_Button(WB, y, 2 * BB + WB, BH, "Print Rho quality measure"); y += BH;
- for(int i = 0; i < 5; i++)
+ dialog->b[6] = new Fl_Check_Button(WB, y, 2 * BB + WB, BH, "Print Disto quality measure"); y += BH;
+ for(int i = 0; i < 6; i++)
dialog->b[i]->type(FL_TOGGLE_BUTTON);
dialog->ok = new Fl_Return_Button(WB, y + WB, BB, BH, "OK");
dialog->cancel = new Fl_Button(2 * WB + BB, y + WB, BB, BH, "Cancel");
@@ -813,6 +814,7 @@ int pos_dialog(const char *name)
opt_print_pos_gamma(0, GMSH_SET | GMSH_GUI, dialog->b[3]->value() ? 1 : 0);
opt_print_pos_eta(0, GMSH_SET | GMSH_GUI, dialog->b[4]->value() ? 1 : 0);
opt_print_pos_rho(0, GMSH_SET | GMSH_GUI, dialog->b[5]->value() ? 1 : 0);
+ opt_print_pos_disto(0, GMSH_SET | GMSH_GUI, dialog->b[6]->value() ? 1 : 0);
CreateOutputFile(name, FORMAT_POS);
dialog->window->hide();
return 1;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index a0bfcff95b4141141952a8a03aa61cb1799de16c..a98a560b8a030e57e698a18c4844da6cea1713a1 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -270,7 +270,7 @@ class GModel
// Mesh statistics (as Gmsh post-processing views)
int writePOS(const std::string &name, bool printElementary,
bool printElementNumber, bool printGamma, bool printEta, bool printRho,
- bool saveAll=false, double scalingFactor=1.0);
+ bool printDisto,bool saveAll=false, double scalingFactor=1.0);
// Stereo lithography format
int readSTL(const std::string &name, double tolerance=1.e-3);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 6483527b4a41d5d9a4448fec1678003b20ce8bbb..7499a99a2729a850e16e0c59375ed5e76f254e60 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -613,7 +613,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
int GModel::writePOS(const std::string &name, bool printElementary,
bool printElementNumber, bool printGamma, bool printEta,
- bool printRho, bool saveAll, double scalingFactor)
+ bool printRho, bool printDisto, bool saveAll, double scalingFactor)
{
FILE *fp = fopen(name.c_str(), "w");
if(!fp){
@@ -621,7 +621,7 @@ int GModel::writePOS(const std::string &name, bool printElementary,
return 0;
}
- bool f[5] = {printElementary, printElementNumber, printGamma, printEta, printRho};
+ bool f[6] = {printElementary, printElementNumber, printGamma, printEta, printRho,printDisto};
bool first = true;
std::string names;
@@ -645,6 +645,10 @@ int GModel::writePOS(const std::string &name, bool printElementary,
if(first) first = false; else names += ",";
names += "\"Rho\"";
}
+ if(f[5]){
+ if(first) first = false; else names += ",";
+ names += "\"Disto\"";
+ }
if(names.empty()) return 0;
@@ -656,7 +660,7 @@ int GModel::writePOS(const std::string &name, bool printElementary,
for(eiter it = firstEdge(); it != lastEdge(); ++it) {
if(saveAll || (*it)->physicals.size()){
for(unsigned int i = 0; i < (*it)->lines.size(); i++)
- (*it)->lines[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4],
+ (*it)->lines[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
scalingFactor, (*it)->tag());
}
}
@@ -664,10 +668,10 @@ int GModel::writePOS(const std::string &name, bool printElementary,
for(fiter it = firstFace(); it != lastFace(); ++it) {
if(saveAll || (*it)->physicals.size()){
for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
- (*it)->triangles[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4],
+ (*it)->triangles[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
scalingFactor, (*it)->tag());
for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
- (*it)->quadrangles[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4],
+ (*it)->quadrangles[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
scalingFactor, (*it)->tag());
}
}
@@ -675,16 +679,16 @@ int GModel::writePOS(const std::string &name, bool printElementary,
for(riter it = firstRegion(); it != lastRegion(); ++it) {
if(saveAll || (*it)->physicals.size()){
for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
- (*it)->tetrahedra[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4],
+ (*it)->tetrahedra[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
scalingFactor, (*it)->tag());
for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
- (*it)->hexahedra[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4],
+ (*it)->hexahedra[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
scalingFactor, (*it)->tag());
for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
- (*it)->prisms[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4],
+ (*it)->prisms[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
scalingFactor, (*it)->tag());
for(unsigned int i = 0; i < (*it)->pyramids.size(); i++)
- (*it)->pyramids[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4],
+ (*it)->pyramids[i]->writePOS(fp, f[0], f[1], f[2], f[3], f[4], f[5],
scalingFactor, (*it)->tag());
}
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index bc476d6eb4e7ad9fa04cfeb17c444789f8a1a44f..20a7f3f10e039e754dd044808495a942c6e89cb2 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -18,6 +18,8 @@
# include "Message.h"
# include "Context.h"
# include "qualityMeasures.h"
+# include "meshGFaceDelaunayInsertion.h"
+# include "meshGRegionDelaunayInsertion.h"
#endif
#define SQU(a) ((a)*(a))
@@ -96,6 +98,51 @@ void MElement::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
Msg::Error("No integration points defined for this type of element");
}
+SPoint3 MTriangle::circumcenter()
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return 0.;
+#else
+ double p1[3] = {_v[0]->x(),_v[0]->y(),_v[0]->z()};
+ double p2[3] = {_v[1]->x(),_v[1]->y(),_v[1]->z()};
+ double p3[3] = {_v[2]->x(),_v[2]->y(),_v[2]->z()};
+ double res[3];
+ circumCenterXYZ(p1,p2,p3,res);
+ return SPoint3(res[0],res[1],res[2]);
+#endif
+}
+
+SPoint3 MTetrahedron::circumcenter()
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return 0.;
+#else
+ MTet4 t(this,0);
+ double res[3];
+ t.circumcenter(res);
+ return SPoint3(res[0],res[1],res[2]);
+#endif
+}
+
+
+double MTriangle::distoShapeMeasure()
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return 1.;
+#else
+ return qmDistorsionOfMapping(this);
+#endif
+}
+
+double MTetrahedron::distoShapeMeasure()
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return 1.;
+#else
+ return qmDistorsionOfMapping(this);
+#endif
+}
+
double MTriangle::gammaShapeMeasure()
{
#if defined(HAVE_GMSH_EMBEDDED)
@@ -420,7 +467,7 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
bool printGamma, bool printEta, bool printRho,
- double scalingFactor, int elementary)
+ bool printDisto, double scalingFactor, int elementary)
{
const char *str = getStringForPOS();
if(!str) return;
@@ -467,6 +514,13 @@ void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
fprintf(fp, "%g", rho);
}
}
+ if(printDisto){
+ double disto = distoShapeMeasure();
+ for(int i = 0; i < n; i++){
+ if(first) first = false; else fprintf(fp, ",");
+ fprintf(fp, "%g", disto);
+ }
+ }
fprintf(fp, "};\n");
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 2012099606a3725b6b3fb24fb23ba7ddcaa60ef9..c976e651db2d9c28b2e6923c0343c65483171b7f 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -138,6 +138,7 @@ class MElement
virtual double rhoShapeMeasure();
virtual double gammaShapeMeasure(){ return 0.; }
virtual double etaShapeMeasure(){ return 0.; }
+ virtual double distoShapeMeasure() {return 1.0;}
// compute the barycenter
virtual SPoint3 barycenter();
@@ -193,7 +194,7 @@ class MElement
int num=0, int elementary=1, int physical=1);
virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber,
bool printGamma, bool printEta, bool printRho,
- double scalingFactor=1.0, int elementary=1);
+ bool printDisto,double scalingFactor=1.0, int elementary=1);
virtual void writeSTL(FILE *fp, bool binary=false, double scalingFactor=1.0);
virtual void writeVRML(FILE *fp);
virtual void writeUNV(FILE *fp, int num=0, int elementary=1, int physical=1);
@@ -459,6 +460,7 @@ class MTriangle : public MElement {
~MTriangle(){}
virtual int getDim(){ return 2; }
virtual double gammaShapeMeasure();
+ virtual double distoShapeMeasure();
virtual int getNumVertices(){ return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
@@ -550,6 +552,7 @@ class MTriangle : public MElement {
pnt(1, 0, u, v, w, p);
}
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual SPoint3 circumcenter();
private:
int edges_tri(const int edge, const int vert) const
{
@@ -1185,6 +1188,7 @@ class MTetrahedron : public MElement {
virtual double getVolume();
virtual int getVolumeSign(){ return (getVolume() >= 0) ? 1 : -1; }
virtual double gammaShapeMeasure();
+ virtual double distoShapeMeasure();
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
virtual void getShapeFunction(int num, double u, double v, double w, double &s)
@@ -1236,6 +1240,7 @@ class MTetrahedron : public MElement {
pnt(1, 0, u, v, w, p);
}
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual SPoint3 circumcenter();
private:
int edges_tetra(const int edge, const int vert) const
{
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 62a8268debbc8be09e66fc5ffca3e23dfb642888..e34dbf68bdcb5104d997f371069fc50a25bb0580 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -99,7 +99,9 @@ static bool isElementVisible(MElement *ele)
if(!ele->getVisibility()) return false;
if(CTX.mesh.quality_sup) {
double q;
- if(CTX.mesh.quality_type == 2)
+ if(CTX.mesh.quality_type == 3)
+ q = ele->distoShapeMeasure();
+ else if(CTX.mesh.quality_type == 2)
q = ele->rhoShapeMeasure();
else if(CTX.mesh.quality_type == 1)
q = ele->etaShapeMeasure();
@@ -363,6 +365,59 @@ static void drawBarycentricDual(std::vector<T*> &elements)
gl2psDisable(GL2PS_LINE_STIPPLE);
}
+template<class T>
+static void drawVoronoiDual(std::vector<T*> &elements)
+{
+ glColor4ubv((GLubyte *) & CTX.color.fg);
+ glEnable(GL_LINE_STIPPLE);
+ glLineStipple(1, 0x0F0F);
+ gl2psEnable(GL2PS_LINE_STIPPLE);
+ glBegin(GL_LINES);
+ for(unsigned int i = 0; i < elements.size(); i++){
+ T *ele = elements[i];
+ if(!isElementVisible(ele)) continue;
+ SPoint3 pc = ele->circumcenter();
+ if(ele->getDim() == 2){
+ for(int j = 0; j < ele->getNumEdges(); j++){
+ MEdge e = ele->getEdge(j);
+ SVector3 p2p1 ( e.getVertex(1)->x() - e.getVertex(0)->x(),
+ e.getVertex(1)->y() - e.getVertex(0)->y(),
+ e.getVertex(1)->z() - e.getVertex(0)->z());
+ SVector3 pcp1 ( pc.x() - e.getVertex(0)->x(),
+ pc.y() - e.getVertex(0)->y(),
+ pc.z() - e.getVertex(0)->z());
+
+ double alpha = dot(pcp1,p2p1) / dot(p2p1,p2p1);
+
+ SPoint3 p ((1.-alpha)*e.getVertex(0)->x() + alpha*e.getVertex(1)->x(),
+ (1.-alpha)*e.getVertex(0)->y() + alpha*e.getVertex(1)->y(),
+ (1.-alpha)*e.getVertex(0)->z() + alpha*e.getVertex(1)->z());
+ glVertex3d(pc.x(), pc.y(), pc.z());
+ glVertex3d(p.x(), p.y(), p.z());
+ }
+ }
+ else if(ele->getDim() == 3){
+ for(int j = 0; j < ele->getNumFaces(); j++){
+ MFace f = ele->getFace(j);
+ SPoint3 p = f.barycenter();
+ glVertex3d(pc.x(), pc.y(), pc.z());
+ glVertex3d(p.x(), p.y(), p.z());
+ for(int k = 0; k < f.getNumVertices(); k++){
+ MEdge e(f.getVertex(k), (k == f.getNumVertices() - 1) ?
+ f.getVertex(0) : f.getVertex(k + 1));
+ SPoint3 pe = e.barycenter();
+ glVertex3d(p.x(), p.y(), p.z());
+ glVertex3d(pe.x(), pe.y(), pe.z());
+ }
+ }
+ }
+ }
+ glEnd();
+ glDisable(GL_LINE_STIPPLE);
+ gl2psDisable(GL2PS_LINE_STIPPLE);
+}
+
+
// Routine to fill the smooth normal container
template<class T>
@@ -696,6 +751,10 @@ class drawMeshGFace {
if(CTX.mesh.triangles) drawBarycentricDual(f->triangles);
if(CTX.mesh.quadrangles) drawBarycentricDual(f->quadrangles);
}
+ else if(CTX.mesh.voronoi) {
+ if(CTX.mesh.triangles) drawVoronoiDual(f->triangles);
+ }
+
if(CTX.render_mode == GMSH_SELECT) {
glPopName();
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 7384b21347a2c4e23f9319a14f517f3c2d1532f7..71d3e7a76d24c4893c55f593fc1edc65b894a263 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -30,7 +30,8 @@ static void GetQualityMeasure(std::vector<T*> &ele,
double &gamma, double &gammaMin, double &gammaMax,
double &eta, double &etaMin, double &etaMax,
double &rho, double &rhoMin, double &rhoMax,
- double quality[3][100])
+ double &disto, double &distoMin, double &distoMax,
+ double quality[4][100])
{
for(unsigned int i = 0; i < ele.size(); i++){
double g = ele[i]->gammaShapeMeasure();
@@ -45,15 +46,20 @@ static void GetQualityMeasure(std::vector<T*> &ele,
rho += r;
rhoMin = std::min(rhoMin, r);
rhoMax = std::max(rhoMax, r);
+ double d = ele[i]->distoShapeMeasure();
+ disto += d;
+ distoMin = std::min(distoMin, d);
+ distoMax = std::max(distoMax, d);
for(int j = 0; j < 100; j++){
if(g > j / 100. && g <= (j + 1) / 100.) quality[0][j]++;
if(e > j / 100. && e <= (j + 1) / 100.) quality[1][j]++;
if(r > j / 100. && r <= (j + 1) / 100.) quality[2][j]++;
+ if(d > j / 100. && d <= (j + 1) / 100.) quality[3][j]++;
}
}
}
-void GetStatistics(double stat[50], double quality[3][100])
+void GetStatistics(double stat[50], double quality[4][100])
{
for(int i = 0; i < 50; i++) stat[i] = 0.;
@@ -99,15 +105,32 @@ void GetStatistics(double stat[50], double quality[3][100])
double gamma=0., gammaMin=1., gammaMax=0.;
double eta=0., etaMin=1., etaMax=0.;
double rho=0., rhoMin=1., rhoMax=0.;
- for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
- GetQualityMeasure((*it)->tetrahedra, gamma, gammaMin, gammaMax,
- eta, etaMin, etaMax, rho, rhoMin, rhoMax, quality);
- GetQualityMeasure((*it)->hexahedra, gamma, gammaMin, gammaMax,
- eta, etaMin, etaMax, rho, rhoMin, rhoMax, quality);
- GetQualityMeasure((*it)->prisms, gamma, gammaMin, gammaMax,
- eta, etaMin, etaMax, rho, rhoMin, rhoMax, quality);
- GetQualityMeasure((*it)->pyramids, gamma, gammaMin, gammaMax,
- eta, etaMin, etaMax, rho, rhoMin, rhoMax, quality);
+ double disto=0., distoMin=1., distoMax=0.;
+ if (m->firstRegion() == m->lastRegion()){
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ GetQualityMeasure((*it)->quadrangles, gamma, gammaMin, gammaMax,
+ eta, etaMin, etaMax, rho, rhoMin, rhoMax,
+ disto, distoMin,distoMax,quality);
+ GetQualityMeasure((*it)->triangles, gamma, gammaMin, gammaMax,
+ eta, etaMin, etaMax, rho, rhoMin, rhoMax,
+ disto, distoMin,distoMax,quality);
+ }
+ }
+ else{
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
+ GetQualityMeasure((*it)->tetrahedra, gamma, gammaMin, gammaMax,
+ eta, etaMin, etaMax, rho, rhoMin, rhoMax,
+ disto, distoMin,distoMax,quality);
+ GetQualityMeasure((*it)->hexahedra, gamma, gammaMin, gammaMax,
+ eta, etaMin, etaMax, rho, rhoMin, rhoMax,
+ disto, distoMin,distoMax,quality);
+ GetQualityMeasure((*it)->prisms, gamma, gammaMin, gammaMax,
+ eta, etaMin, etaMax, rho, rhoMin, rhoMax,
+ disto, distoMin,distoMax,quality);
+ GetQualityMeasure((*it)->pyramids, gamma, gammaMin, gammaMax,
+ eta, etaMin, etaMax, rho, rhoMin, rhoMax,
+ disto, distoMin,distoMax,quality);
+ }
}
double N = stat[9] + stat[10] + stat[11] + stat[12];
stat[17] = N ? gamma / N : 0.;
@@ -119,6 +142,9 @@ void GetStatistics(double stat[50], double quality[3][100])
stat[23] = N ? rho / N : 0;
stat[24] = rhoMin;
stat[25] = rhoMax;
+ stat[46] = N ? disto / N : 0;
+ stat[47] = distoMin;
+ stat[48] = distoMax;
}
#if !defined(HAVE_NO_POST)
diff --git a/Mesh/Generator.h b/Mesh/Generator.h
index 8c6959fec3934b0daa13b84601dc1417e01fb940..bef74f51ad754824e39c1a7331a262ed2c754eb8 100644
--- a/Mesh/Generator.h
+++ b/Mesh/Generator.h
@@ -8,7 +8,7 @@
class GModel;
-void GetStatistics(double stat[50], double quality[3][100]=0);
+void GetStatistics(double stat[50], double quality[4][100]=0);
void AdaptMesh(GModel *m);
void GenerateMesh(GModel *m, int dimension);
void OptimizeMesh(GModel *m);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index abc6b5b3ea27d36e71a56c6dfabfee0bd78c43f1..985c21489a52f60a080439b595b5667499b507cf 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -709,7 +709,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
// vertices
if(AlgoDelaunay2D(gf)){
gmshBowyerWatson(gf);
- laplaceSmoothing(gf);
+ for (int i=0;i<CTX.mesh.nb_smoothing;i++)laplaceSmoothing(gf);
}
else if (debug){
char name[256];
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 1193db51e6248fef4c0de96600ec6b9f4989f0ce..f3b45dfc8e8845d24a271b747ee24e4985ebdc20 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -171,3 +171,56 @@ double qmTet(const double &x1, const double &y1, const double &z1,
return 0.;
}
}
+
+
+static double mesh_functional_distorsion(MTriangle *t, double u, double v)
+{
+ // compute uncurved element jacobian d_u x and d_v x
+ double mat[2][3];
+ t->jac(1, 0, 0, 0, 0, mat);
+ double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
+ double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
+ double normal1[3];
+ prodve(v1, v2, normal1);
+ double nn = sqrt(normal1[0]*normal1[0] +
+ normal1[1]*normal1[1] +
+ normal1[2]*normal1[2]);
+
+ // compute uncurved element jacobian d_u x and d_v x
+ t->jac(u, v, 0, mat);
+ double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
+ double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
+ double normal[3];
+ prodve(v1b, v2b, normal);
+
+ double sign;
+ prosca(normal1, normal, &sign);
+ double det = norm3(normal) * (sign > 0 ? 1. : -1.) / nn;
+
+ // compute distorsion
+ double dist = std::min(1. / det, det);
+ return dist;
+}
+
+
+double qmDistorsionOfMapping (MTriangle *e)
+{
+ if (e->getPolynomialOrder() == 1)return 1.0;
+ IntPt *pts;
+ int npts;
+ e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
+ double dmin;
+ for (int i=0;i<npts;i++){
+ const double u = pts[i].pt[0];
+ const double v = pts[i].pt[1];
+ const double w = pts[i].pt[2];
+ const double di = mesh_functional_distorsion (e,u,v);
+ dmin = (i==0)? di : std::min(dmin,di);
+ }
+ return dmin;
+}
+
+double qmDistorsionOfMapping (MTetrahedron *e)
+{
+ return 1.0;
+}
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index 06ba20f05cba540cf20464e184eaf91f53b57491..198c4188c6f43a38ba3d24205620a51a1a059ded 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -14,6 +14,10 @@ class MTetrahedron;
enum gmshQualityMeasure4Triangle {QMTRI_RHO, QMTRI_COND};
enum gmshQualityMeasure4Tet{QMTET_1, QMTET_2, QMTET_3, QMTET_ONE, QMTET_COND};
+
+double qmDistorsionOfMapping (MTriangle *e);
+double qmDistorsionOfMapping (MTetrahedron *e);
+
double qmTriangle(MTriangle *f, const gmshQualityMeasure4Triangle &cr);
double qmTriangle(BDS_Face *f, const gmshQualityMeasure4Triangle &cr);
double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,