diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index e05057fb1c4567173f900f403f16c21ad944b885..729ac6f88e853f91f1cd6687dc53799481049b71 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,3 +1,5 @@
+3.0.3 (June 27, 2017): new element quality measures; Block->Box; minor fixes.
+
3.0.2 (May 13, 2017): improved handling of meshing constraints and entity
numbering after boolean operations; improved handling of fast coarseness
transitions in MeshAdapt; new TIKZ export; small bug fixes.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9ebe2d1c9d9d5edcc66202f41e43071e32002d70..cfa3b2c2d87acc4bb6e6755c25adc13f09c86d73 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -93,7 +93,7 @@ opt(ZIPPER "Enable Zip file compression/decompression" OFF)
set(GMSH_MAJOR_VERSION 3)
set(GMSH_MINOR_VERSION 0)
-set(GMSH_PATCH_VERSION 2)
+set(GMSH_PATCH_VERSION 3)
set(GMSH_EXTRA_VERSION "" CACHE STRING "Gmsh extra version string")
set(GMSH_VERSION "${GMSH_MAJOR_VERSION}.${GMSH_MINOR_VERSION}")
diff --git a/Common/Context.h b/Common/Context.h
index 10675fc468bed2b046b7514f55906a0ab747206f..762b7380e327b7b00e29f6d67dea6047411c0dd0 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -279,7 +279,7 @@ class CTX {
int jpegQuality, jpegSmoothing, geoLabels, geoOnlyPhysicals;
int text, texAsEquation;
int gifDither, gifSort, gifInterlace, gifTransparent;
- int posElementary, posElement, posGamma, posEta, posRho, posDisto;
+ int posElementary, posElement, posGamma, posEta, posSICN, posSIGE, posDisto;
int compositeWindows, deleteTmpFiles, background;
int width, height;
double parameter, parameterFirst, parameterLast, parameterSteps;
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index f744e81f1ec5484322bf62cc01199dc48c4b51c5..7ec744426dacf7c2091fc665bf8c63b64f8c7d16 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -377,10 +377,10 @@ void CreateOutputFile(const std::string &fileName, int format,
case FORMAT_POS:
GModel::current()->writePOS
(name, CTX::instance()->print.posElementary,
- CTX::instance()->print.posElement, CTX::instance()->print.posGamma,
- CTX::instance()->print.posEta, CTX::instance()->print.posRho,
- CTX::instance()->print.posDisto, CTX::instance()->mesh.saveAll,
- CTX::instance()->mesh.scalingFactor);
+ CTX::instance()->print.posElement,
+ CTX::instance()->print.posSICN, CTX::instance()->print.posSIGE,
+ CTX::instance()->print.posGamma, CTX::instance()->print.posDisto,
+ CTX::instance()->mesh.saveAll, CTX::instance()->mesh.scalingFactor);
break;
case FORMAT_GEO:
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 430bee5ec9142a331bce6d6a6233cd27d71df1ca..7dbcb5cf43f9cd03ff9808538adc88debf199417 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1193,8 +1193,9 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "QualitySup" , opt_mesh_quality_sup , 0.0 ,
"Only display elements whose quality measure is smaller than QualitySup" },
{ F|O, "QualityType" , opt_mesh_quality_type , 2. ,
- "Type of quality measure (0=gamma~vol/sum_face/max_edge, "
- "1=eta~vol^(2/3)/sum_edge^2, 2=rho~min_edge/max_edge)" },
+ "Type of quality measure (0=SICN~signed inverse condition number, "
+ "1=SIGE~signed inverse gradient error, 2=gamma~vol/sum_face/max_edge, "
+ "3=Disto~minJ/maxJ"},
{ F|O, "RadiusInf" , opt_mesh_radius_inf , 0.0 ,
"Only display elements whose longest edge is greater than RadiusInf" },
@@ -1760,9 +1761,12 @@ StringXNumber PrintOptions_Number[] = {
{ F|O, "PostEta" , opt_print_pos_eta , 0. ,
"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, "PostSICN" , opt_print_pos_SICN , 0. ,
+ "Save SICN (signed inverse condition number) quality measure in mesh "
+ "statistics exported as post-processing views" },
+ { F|O, "PostSICN" , opt_print_pos_SIGE , 0. ,
+ "Save SIGE (signed inverse gradient error) 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" },
diff --git a/Common/GmshMessage.cpp b/Common/GmshMessage.cpp
index 24cbce9dc590e99ec85b3886c9e72d041052e5e0..0709a85dbb5597d112dd66ec8adcc9818860a1c0 100644
--- a/Common/GmshMessage.cpp
+++ b/Common/GmshMessage.cpp
@@ -773,10 +773,9 @@ void Msg::ProgressMeter(int n, int N, bool log, const char *fmt, ...)
(n > N - 1) ? "" : str2);
fflush(stdout);
}
-
-#if defined(_OPENMP)
-#pragma omp barrier
-#endif
+ //#if defined(_OPENMP)
+ //#pragma omp barrier
+ //#endif
}
}
diff --git a/Common/Options.cpp b/Common/Options.cpp
index da06739e14cd8711f6931759f3c73277bab8778d..6e0885b9a97e7c23521e2b42460d713a14e2ced9 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -9182,11 +9182,18 @@ double opt_print_pos_eta(OPT_ARGS_NUM)
return CTX::instance()->print.posEta;
}
-double opt_print_pos_rho(OPT_ARGS_NUM)
+double opt_print_pos_SICN(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
- CTX::instance()->print.posRho = (int)val;
- return CTX::instance()->print.posRho;
+ CTX::instance()->print.posSICN = (int)val;
+ return CTX::instance()->print.posSICN;
+}
+
+double opt_print_pos_SIGE(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX::instance()->print.posSIGE = (int)val;
+ return CTX::instance()->print.posSIGE;
}
double opt_print_pos_disto(OPT_ARGS_NUM)
diff --git a/Common/Options.h b/Common/Options.h
index bdf53ffa1f505f8ea39618e40da62806ca4b66ef..1acf4fd406fdfa04087474ead536ce216374b1d6 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -719,7 +719,8 @@ double opt_print_pos_elementary(OPT_ARGS_NUM);
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_SICN(OPT_ARGS_NUM);
+double opt_print_pos_SIGE(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);
diff --git a/Fltk/fileDialogs.cpp b/Fltk/fileDialogs.cpp
index b0e7ce1e1c36e328819816fbb2a319dede208f61..7250892adaed22d526919038710c9592c969923d 100644
--- a/Fltk/fileDialogs.cpp
+++ b/Fltk/fileDialogs.cpp
@@ -912,7 +912,7 @@ int meshStatFileDialog(const char *name)
{
struct _meshStatFileDialog{
Fl_Window *window;
- Fl_Check_Button *b[7];
+ Fl_Check_Button *b[8];
Fl_Button *ok, *cancel;
};
static _meshStatFileDialog *dialog = NULL;
@@ -921,7 +921,7 @@ int meshStatFileDialog(const char *name)
if(!dialog){
dialog = new _meshStatFileDialog;
- int h = 3 * WB + 8 * BH, w = 2 * BBB + 3 * WB, y = WB;
+ int h = 3 * WB + 9 * BH, w = 2 * BBB + 3 * WB, y = WB;
dialog->window = new Fl_Double_Window(w, h, "POS Options");
dialog->window->box(GMSH_WINDOW_BOX);
dialog->window->set_modal();
@@ -932,12 +932,14 @@ int meshStatFileDialog(const char *name)
dialog->b[2] = new Fl_Check_Button
(WB, y, 2 * BBB + WB, BH, "Print element numbers"); y += BH;
dialog->b[3] = new Fl_Check_Button
- (WB, y, 2 * BBB + WB, BH, "Print Gamma quality measure"); y += BH;
+ (WB, y, 2 * BBB + WB, BH, "Print SICN quality measure"); y += BH;
dialog->b[4] = new Fl_Check_Button
- (WB, y, 2 * BBB + WB, BH, "Print Eta quality measure"); y += BH;
+ (WB, y, 2 * BBB + WB, BH, "Print SIGE quality measure"); y += BH;
dialog->b[5] = new Fl_Check_Button
- (WB, y, 2 * BBB + WB, BH, "Print Rho quality measure"); y += BH;
+ (WB, y, 2 * BBB + WB, BH, "Print Gamma quality measure"); y += BH;
dialog->b[6] = new Fl_Check_Button
+ (WB, y, 2 * BBB + WB, BH, "Print Eta quality measure"); y += BH;
+ dialog->b[7] = new Fl_Check_Button
(WB, y, 2 * BBB + WB, BH, "Print Disto quality measure"); y += BH;
for(int i = 0; i < 6; i++)
dialog->b[i]->type(FL_TOGGLE_BUTTON);
@@ -950,9 +952,10 @@ int meshStatFileDialog(const char *name)
dialog->b[0]->value(CTX::instance()->mesh.saveAll ? 1 : 0);
dialog->b[1]->value(CTX::instance()->print.posElementary ? 1 : 0);
dialog->b[2]->value(CTX::instance()->print.posElement ? 1 : 0);
- dialog->b[3]->value(CTX::instance()->print.posGamma ? 1 : 0);
- dialog->b[4]->value(CTX::instance()->print.posEta ? 1 : 0);
- dialog->b[5]->value(CTX::instance()->print.posRho ? 1 : 0);
+ dialog->b[3]->value(CTX::instance()->print.posSICN ? 1 : 0);
+ dialog->b[4]->value(CTX::instance()->print.posSIGE ? 1 : 0);
+ dialog->b[5]->value(CTX::instance()->print.posGamma ? 1 : 0);
+ dialog->b[6]->value(CTX::instance()->print.posEta ? 1 : 0);
dialog->window->show();
while(dialog->window->shown()){
@@ -964,10 +967,13 @@ int meshStatFileDialog(const char *name)
opt_mesh_save_all(0, GMSH_SET | GMSH_GUI, dialog->b[0]->value() ? 1 : 0);
opt_print_pos_elementary(0, GMSH_SET | GMSH_GUI, dialog->b[1]->value() ? 1 : 0);
opt_print_pos_element(0, GMSH_SET | GMSH_GUI, dialog->b[2]->value() ? 1 : 0);
- 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);
+ opt_print_pos_SICN(0, GMSH_SET | GMSH_GUI,
+ dialog->b[3]->value() ? 1 : 0);
+ opt_print_pos_SIGE(0, GMSH_SET | GMSH_GUI,
+ dialog->b[4]->value() ? 1 : 0);
+ opt_print_pos_gamma(0, GMSH_SET | GMSH_GUI, dialog->b[5]->value() ? 1 : 0);
+ opt_print_pos_eta(0, GMSH_SET | GMSH_GUI, dialog->b[6]->value() ? 1 : 0);
+ opt_print_pos_disto(0, GMSH_SET | GMSH_GUI, dialog->b[7]->value() ? 1 : 0);
CreateOutputFile(name, FORMAT_POS);
dialog->window->hide();
return 1;
diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index cbd150631e3c0902108d43bba3996fe34660187c..f75ecfbacd3d28d0a62d90077a2d63ffb24ec9b7 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -2086,8 +2086,7 @@ static std::vector<std::string> getInfoStrings(MElement *ele)
std::ostringstream sstream;
sstream.precision(12);
sstream << " Quality: "
- << "gamma = " << ele->gammaShapeMeasure() << " "
- << "rho = " << ele->rhoShapeMeasure();
+ << "gamma = " << ele->gammaShapeMeasure();
info.push_back(sstream.str());
}
{
@@ -2098,6 +2097,14 @@ static std::vector<std::string> getInfoStrings(MElement *ele)
sstream << " SICN range: " << sICNMin << " " << sICNMax;
info.push_back(sstream.str());
}
+ {
+ std::ostringstream sstream;
+ sstream.precision(12);
+ double sIGEMin, sIGEMax;
+ ele->signedInvGradErrorRange(sIGEMin, sIGEMax);
+ sstream << " SIGE range: " << sIGEMin << " " << sIGEMax;
+ info.push_back(sstream.str());
+ }
{
std::ostringstream sstream;
sstream.precision(12);
diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index 5f30aa793af0bef8549ecc26b5243e33bacfb214..e2358c070dbbf896caa71cdb692e4f854238cacb 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -2518,8 +2518,8 @@ optionWindow::optionWindow(int deltaFontSize)
static Fl_Menu_Item menu_quality_type[] = {
{"SICN", 0, 0, 0},
+ {"SIGE", 0, 0, 0},
{"Gamma", 0, 0, 0},
- {"Rho", 0, 0, 0},
{"Disto", 0, 0, 0},
{0}
};
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 7f5e88acd600976d58dc231997b23040a1749989..29364f45086ee459462720522910263ff5fd2bb7 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -18,7 +18,9 @@
#include "OS.h"
#include "Field.h"
-enum QM_HISTO {QMH_SICN_XY, QMH_SICN_3D, QMH_GAMMA_XY, QMH_GAMMA_3D, QMH_RHO_XY, QMH_RHO_3D};
+enum QM_HISTO {QMH_SICN_XY, QMH_SICN_3D,
+ QMH_GAMMA_XY, QMH_GAMMA_3D,
+ QMH_SIGE_XY, QMH_SIGE_3D};
void statistics_cb(Fl_Widget *w, void *data)
{
@@ -50,12 +52,12 @@ static void statistics_histogram_cb(Fl_Widget *w, void *data)
}
new PView("Gamma", "# Elements", x, y);
}
- else if (qmh == QMH_RHO_XY) {
+ else if (qmh == QMH_SIGE_XY) {
for(int i = 0; i < 100; i++){
- x.push_back((double)i / 99);
- y.push_back(FlGui::instance()->stats->quality[2][i]);
+ x.push_back((double)(2*i-99) / 99);
+ y.push_back(FlGui::instance()->stats->quality[3][i]);
}
- new PView("Rho", "# Elements", x, y);
+ new PView("SIGE", "# Elements", x, y);
}
else {
std::vector<GEntity*> entities_;
@@ -69,13 +71,13 @@ static void statistics_histogram_cb(Fl_Widget *w, void *data)
d[e->getNum()].push_back(e->minSICNShapeMeasure());
else if (qmh == QMH_GAMMA_3D)
d[e->getNum()].push_back(e->gammaShapeMeasure());
- else if (qmh == QMH_RHO_3D)
- d[e->getNum()].push_back(e->rhoShapeMeasure());
+ else if (qmh == QMH_SIGE_3D)
+ d[e->getNum()].push_back(e->minSIGEShapeMeasure());
}
}
std::string name = (qmh == QMH_SICN_3D) ? "SICN" :
(qmh == QMH_GAMMA_3D) ? "Gamma" :
- (qmh == QMH_RHO_3D) ? "Rho" : "";
+ (qmh == QMH_SIGE_3D) ? "SIGE" : "";
new PView(name, "ElementData", GModel::current(), d);
}
@@ -128,8 +130,8 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
value[num]->tooltip("~ signed inverse condition number"); num++;
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 + 16 * BH, IW, BH, "Rho");
- value[num]->tooltip("~ min_edge_length / max_edge_length"); num++;
+ value[num] = new Fl_Output(2 * WB, 2 * WB + 16 * BH, IW, BH, "SIGE");
+ value[num]->tooltip("~ signed inverse error on gradient FE solution"); num++;
for(int i = 0; i < 3; i++){
int ww = 3 * FL_NORMAL_SIZE;
@@ -141,7 +143,7 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
(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;
+ qmh3 = QMH_GAMMA_3D, qmh4 = QMH_SIGE_XY, qmh5 = QMH_SIGE_3D;
butt[0]->callback(statistics_histogram_cb, (void*) &qmh0);
butt[1]->callback(statistics_histogram_cb, (void*) &qmh1);
butt[2]->callback(statistics_histogram_cb, (void*) &qmh2);
diff --git a/Fltk/viewButton.cpp b/Fltk/viewButton.cpp
index 36110cfca75235e354075de2bd579c212afcee0f..94af2572a2b25b3604a3e7a4f4894bc37b998881 100644
--- a/Fltk/viewButton.cpp
+++ b/Fltk/viewButton.cpp
@@ -32,7 +32,14 @@ static void view_toggle_cb(Fl_Widget *w, void *data)
int num = (intptr_t)data;
viewButton *but = FlGui::instance()->onelab->getViewButton(num);
if(but){
- opt_view_visible(num, GMSH_SET, but->value());
+ if(Fl::event_state(FL_SHIFT)){
+ for(int i = 0; i < PView::list.size(); i++){
+ if(i != num) opt_view_visible(i, GMSH_SET | GMSH_GUI, 0);
+ else opt_view_visible(i, GMSH_SET | GMSH_GUI, 1);
+ }
+ }
+ else
+ opt_view_visible(num, GMSH_SET, but->value());
drawContext::global()->draw();
}
}
diff --git a/Geo/GModel.h b/Geo/GModel.h
index b0f99a3b724d2937e9873201fb8016663bf4a8c4..e292adf43f030b3fc9c63415c21ad805693b6038 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -633,8 +633,9 @@ class GModel {
// mesh statistics (saved as a Gmsh post-processing view)
int writePOS(const std::string &name, bool printElementary,
- bool printElementNumber, bool printSICN, bool printGamma, bool printRho,
- bool printDisto, bool saveAll=false, double scalingFactor=1.0);
+ bool printElementNumber, bool printSICN, bool printSIGE,
+ bool printGamma, 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_POS.cpp b/Geo/GModelIO_POS.cpp
index 10b4fb39a3ad2fd0edc50d9a64f6d6a26ce82016..dfed5982be90e97f84e8a3c9354b09ce597bd951 100644
--- a/Geo/GModelIO_POS.cpp
+++ b/Geo/GModelIO_POS.cpp
@@ -9,8 +9,8 @@
#include "MElement.h"
int GModel::writePOS(const std::string &name, bool printElementary,
- bool printElementNumber, bool printSICN, bool printGamma,
- bool printRho, bool printDisto,
+ bool printElementNumber, bool printSICN, bool printSIGE,
+ bool printGamma, bool printDisto,
bool saveAll, double scalingFactor)
{
FILE *fp = Fopen(name.c_str(), "w");
@@ -36,32 +36,29 @@ int GModel::writePOS(const std::string &name, bool printElementary,
}
*/
- bool f[6] = {printElementary, printElementNumber, printSICN, printGamma, printRho,
- printDisto};
-
bool first = true;
std::string names;
- if(f[0]){
+ if(printElementary){
if(first) first = false; else names += ",";
names += "\"Elementary Entity\"";
}
- if(f[1]){
+ if(printElementNumber){
if(first) first = false; else names += ",";
names += "\"Element Number\"";
}
- if(f[2]){
+ if(printSICN){
if(first) first = false; else names += ",";
names += "\"SICN\"";
}
- if(f[3]){
+ if(printSIGE){
if(first) first = false; else names += ",";
- names += "\"Gamma\"";
+ names += "\"SIGE\"";
}
- if(f[4]){
+ if(printGamma){
if(first) first = false; else names += ",";
- names += "\"Rho\"";
+ names += "\"Gamma\"";
}
- if(f[5]){
+ if(printDisto){
if(first) first = false; else names += ",";
names += "\"Disto\"";
}
@@ -79,7 +76,8 @@ int GModel::writePOS(const std::string &name, bool printElementary,
if(saveAll || entities[i]->physicals.size())
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
entities[i]->getMeshElement(j)->writePOS
- (fp, f[0], f[1], f[2], f[3], f[4], f[5], scalingFactor, entities[i]->tag());
+ (fp, printElementary, printElementNumber, printSICN, printSIGE,
+ printGamma, printDisto, scalingFactor, entities[i]->tag());
fprintf(fp, "};\n");
fclose(fp);
diff --git a/Geo/GModelVertexArrays.cpp b/Geo/GModelVertexArrays.cpp
index 21bde1de96bd979bbe5a9348fde0d21f525d0f4a..02567d774ba73d1a3bd333482b7ec07ee3e0ec22 100644
--- a/Geo/GModelVertexArrays.cpp
+++ b/Geo/GModelVertexArrays.cpp
@@ -106,9 +106,9 @@ bool isElementVisible(MElement *ele)
if(CTX::instance()->mesh.qualityType == 3)
q = ele->distoShapeMeasure();
else if(CTX::instance()->mesh.qualityType == 2)
- q = ele->rhoShapeMeasure();
- else if(CTX::instance()->mesh.qualityType == 1)
q = ele->gammaShapeMeasure();
+ else if(CTX::instance()->mesh.qualityType == 1)
+ q = ele->minSIGEShapeMeasure();
else
q = ele->minSICNShapeMeasure();
if(q < CTX::instance()->mesh.qualityInf ||
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index e60ec2cd62ae8160da0ab252ed7d68402ccfaae9..8b6a0454c79a34c6edb285cb5ac1492fcae095fa 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2429,6 +2429,7 @@ static void RemoveDegenerateSurfaces()
List_Delete (ll);
List_Delete (ll2);
}
+ printf("Deleting Surface %d\n",s->Num);
DeleteSurface(s->Num);
}
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e76d4b04a7c5584cd355b831cb45ca8e69a01887..d71b0ff71dac7da7b3b55debdf944470e000bef8 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -116,16 +116,6 @@ double MElement::maxEdge()
return m;
}
-double MElement::rhoShapeMeasure()
-{
- double min = minEdge();
- double max = maxEdge();
- if(max)
- return min / max;
- else
- return 0.;
-}
-
double MElement::maxDistToStraight() const
{
const nodalBasis *lagBasis = getFunctionSpace();
@@ -152,7 +142,7 @@ double MElement::maxDistToStraight() const
double MElement::minIsotropyMeasure(bool knownValid, bool reversedOK)
{
#if defined(HAVE_MESH)
- return jacobianBasedQuality::minIsotropyMeasure(this, knownValid, reversedOK);
+ return jacobianBasedQuality::minICNMeasure(this, knownValid, reversedOK);
#else
return 0.;
#endif
@@ -161,7 +151,7 @@ double MElement::minIsotropyMeasure(bool knownValid, bool reversedOK)
double MElement::minScaledJacobian(bool knownValid, bool reversedOK)
{
#if defined(HAVE_MESH)
- return jacobianBasedQuality::minScaledJacobian(this, knownValid, reversedOK);
+ return jacobianBasedQuality::minIGEMeasure(this, knownValid, reversedOK);
#else
return 0.;
#endif
@@ -175,7 +165,7 @@ double MElement::specialQuality()
if (minJ <= 0.) return minJ;
// if (minJ < 0 && maxJ >= 0) return minJ/maxJ; // accept -inf as an answer
// if (minJ < 0 && maxJ < 0) return -std::numeric_limits<double>::infinity();
- return jacobianBasedQuality::minIsotropyMeasure(this, true);
+ return jacobianBasedQuality::minICNMeasure(this, true);
#else
return 0;
#endif
@@ -189,7 +179,7 @@ double MElement::specialQuality2()
if (minJ <= 0.) return minJ;
// if (minJ < 0 && maxJ >= 0) return minJ/maxJ; // accept -inf as an answer
// if (minJ < 0 && maxJ < 0) return -std::numeric_limits<double>::infinity();
- return jacobianBasedQuality::minScaledJacobian(this, true);
+ return jacobianBasedQuality::minIGEMeasure(this, true);
#else
return 0;
#endif
@@ -334,6 +324,13 @@ void MElement::signedInvCondNumRange(double &iCNMin, double &iCNMax, GEntity *ge
#endif
}
+void MElement::signedInvGradErrorRange(double &minSIGE, double &maxSIGE)
+{
+ jacobianBasedQuality::sampleIGEMeasure(this, getPolynomialOrder(),
+ minSIGE, maxSIGE);
+ return;
+}
+
void MElement::getNode(int num, double &u, double &v, double &w) const
{
// only for MElements that don't have a lookup table for this
@@ -1135,7 +1132,7 @@ void MElement::writeMSH2(FILE *fp, double version, bool binary, int num,
}
void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
- bool printSICN, bool printGamma, bool printRho,
+ bool printSICN, bool printSIGE, bool printGamma,
bool printDisto, double scalingFactor, int elementary)
{
const char *str = getStringForPOS();
@@ -1169,19 +1166,19 @@ void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
fprintf(fp, "%g", sICNMin);
}
}
- if(printGamma){
- double gamma = gammaShapeMeasure();
+ if(printSIGE){
+ double sIGEMin = minSIGEShapeMeasure();
for(int i = 0; i < n; i++){
if(first) first = false; else fprintf(fp, ",");
- fprintf(fp, "%g", gamma);
- //fprintf(fp, "%d", getVertex(i)->getNum());
+ fprintf(fp, "%g", sIGEMin);
}
}
- if(printRho){
- double rho = rhoShapeMeasure();
+ if(printGamma){
+ double gamma = gammaShapeMeasure();
for(int i = 0; i < n; i++){
if(first) first = false; else fprintf(fp, ",");
- fprintf(fp, "%g", rho);
+ fprintf(fp, "%g", gamma);
+ //fprintf(fp, "%d", getVertex(i)->getNum());
}
}
if(printDisto){
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7fa7470a99a63229ff20e7c51b50bf684b5f1cab..20b7f37ffee9b91ce44b6cb17b833d9015387537 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -199,21 +199,26 @@ class MElement
// get the quality measures
double skewness();
- virtual double rhoShapeMeasure();
virtual double gammaShapeMeasure(){ return 0.; }
virtual double etaShapeMeasure(){ return 0.; }
- double distoShapeMeasure()
- {
- double jmin, jmax;
- scaledJacRange(jmin, jmax);
- return jmin;
- }
double minSICNShapeMeasure()
{
double sICNMin, sICNMax;
signedInvCondNumRange(sICNMin, sICNMax);
return sICNMin;
}
+ double minSIGEShapeMeasure()
+ {
+ double minSIGE, maxSIGE;
+ signedInvGradErrorRange(minSIGE, maxSIGE);
+ return minSIGE;
+ }
+ double distoShapeMeasure()
+ {
+ double jmin, jmax;
+ scaledJacRange(jmin, jmax);
+ return jmin;
+ }
double minIsotropyMeasure(bool knownValid = false, bool reversedOk = false);
double minScaledJacobian(bool knownValid = false, bool reversedOk = false);
double specialQuality();
@@ -222,6 +227,7 @@ class MElement
virtual void scaledJacRange(double &jmin, double &jmax, GEntity *ge = 0) const;
virtual void idealJacRange(double &jmin, double &jmax, GEntity *ge = 0);
virtual void signedInvCondNumRange(double &iCNMin, double &iCNMax, GEntity *ge = 0);
+ virtual void signedInvGradErrorRange(double &minSIGE, double &maxSIGE);
// get the radius of the inscribed circle/sphere if it exists,
// otherwise get the minimum radius of all the circles/spheres
@@ -369,7 +375,7 @@ class MElement
int parentNum=0, int dom1Num = 0, int dom2Num = 0,
std::vector<short> *ghosts=0);
virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber,
- bool printSICN, bool printGamma, bool printRho,
+ bool printSICN, bool printSIGE, bool printGamma,
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);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 0b4ddd66a987bba024c8255dc9eba68c2167587a..ac0c7aa834457fb1a1c1389a64329de11b8ea1dd 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -123,7 +123,7 @@ template<class T>
static void GetQualityMeasure(std::vector<T*> &ele,
double &gamma, double &gammaMin, double &gammaMax,
double &minSICN, double &minSICNMin, double &minSICNMax,
- double &rho, double &rhoMin, double &rhoMax,
+ double &minSIGE, double &minSIGEMin, double &minSIGEMax,
double quality[3][100])
{
for(unsigned int i = 0; i < ele.size(); i++){
@@ -135,14 +135,14 @@ static void GetQualityMeasure(std::vector<T*> &ele,
minSICN += s;
minSICNMin = std::min(minSICNMin, s);
minSICNMax = std::max(minSICNMax, s);
- double r = ele[i]->rhoShapeMeasure();
- rho += r;
- rhoMin = std::min(rhoMin, r);
- rhoMax = std::max(rhoMax, r);
+ double e = ele[i]->minSIGEShapeMeasure();
+ minSIGE += e;
+ minSIGEMin = std::min(minSIGEMin, e);
+ minSIGEMax = std::max(minSIGEMax, e);
for(int j = 0; j < 100; j++){
if(s > (2*j-100) / 100. && s <= (2*j-98) / 100.) quality[0][j]++;
if(g > j / 100. && g <= (j + 1) / 100.) quality[1][j]++;
- if(r > j / 100. && r <= (j + 1) / 100.) quality[2][j]++;
+ if(e > (2*j-100) / 100. && e <= (2*j-98) / 100.) quality[2][j]++;
}
}
}
@@ -190,39 +190,45 @@ void GetStatistics(double stat[50], double quality[3][100])
stat[16] = CTX::instance()->meshTimer[2];
if(quality){
- for(int i = 0; i < 3; i++)
+ for(int i = 0; i < 4; i++)
for(int j = 0; j < 100; j++)
quality[i][j] = 0.;
double minSICN = 0., minSICNMin = 1., minSICNMax = -1.;
+ double minSIGE = 0., minSIGEMin = 1., minSIGEMax = -1.;
double gamma = 0., gammaMin = 1., gammaMax = 0.;
- double rho = 0., rhoMin = 1., rhoMax = 0.;
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,
- minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
+ minSICN, minSICNMin, minSICNMax,
+ minSIGE, minSIGEMin, minSIGEMax, quality);
GetQualityMeasure((*it)->hexahedra, gamma, gammaMin, gammaMax,
- minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
+ minSICN, minSICNMin, minSICNMax,
+ minSIGE, minSIGEMin, minSIGEMax, quality);
GetQualityMeasure((*it)->prisms, gamma, gammaMin, gammaMax,
- minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
+ minSICN, minSICNMin, minSICNMax,
+ minSIGE, minSIGEMin, minSIGEMax, quality);
GetQualityMeasure((*it)->pyramids, gamma, gammaMin, gammaMax,
- minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
+ minSICN, minSICNMin, minSICNMax,
+ minSIGE, minSIGEMin, minSIGEMax, quality);
}
}
else{ // 2D elements
N = stat[7] + stat[8];
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
GetQualityMeasure((*it)->quadrangles, gamma, gammaMin, gammaMax,
- minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
+ minSICN, minSICNMin, minSICNMax,
+ minSIGE, minSIGEMin, minSIGEMax, quality);
GetQualityMeasure((*it)->triangles, gamma, gammaMin, gammaMax,
- minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
+ minSICN, minSICNMin, minSICNMax,
+ minSIGE, minSIGEMin, minSIGEMax, quality);
}
}
if(N){
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;
+ stat[24] = minSIGE / N; stat[25] = minSIGEMin; stat[26] = minSIGEMax;
}
}
diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index 4681034471aafc50aa4a7074f90229720d64473d..f47a7eb60aa4fe85cb0775351a6f587782139d65 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -1124,7 +1124,8 @@ void delaunayTrgl (const unsigned int numThreads,
unsigned int Npts,
std::vector<Vert*> assignTo[],
tetContainer &allocator,
- edgeContainer *embeddedEdges)
+ edgeContainer *embeddedEdges,
+ bool filter)
{
#if defined(_VERBOSE)
double totSearchGlob=0;
diff --git a/Mesh/delaunay3d_private.h b/Mesh/delaunay3d_private.h
index c67f9480d0fc5983d6bd6866f4beaf47cd71c689..0c7871bb8f34db8c4fffdd4e1936d2160c5e80ee 100644
--- a/Mesh/delaunay3d_private.h
+++ b/Mesh/delaunay3d_private.h
@@ -405,7 +405,8 @@ void delaunayTrgl(const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
std::vector<Vert*> assignTo[],
- tetContainer &allocator, edgeContainer *embedded = 0);
+ tetContainer &allocator, edgeContainer *embedded = 0,
+ bool filter = false);
void edgeSwapPass (int numThreads, tetContainer &allocator, edgeContainer &embeddedEdges);
void vertexRelocationPass (int numThreads, std::vector<Vert*> &v);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index d7da61613f6dab1ccc7cd0874234f8c8be7ce7ce..41d2b7afa3845de83ce7ac55f1eeb51b8e3fbf89 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2415,7 +2415,7 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
-int debugSurface = -100; //-100;
+int debugSurface = -1; //-100;
void meshGFace::operator() (GFace *gf, bool print)
{
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 3616c8105757256abc1ae9192195069eb5e334fe..8ea7dfe613a992e8e5da51a59b3fb34b98498e68 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -433,7 +433,7 @@ void printTets (const char *fn, std::list<MTet4*> &cavity, bool force = false )
MTet4 *tet = *ittet;
if (force || !tet->isDeleted()){
MTetrahedron *t = tet->tet();
- t->writePOS (f, false,false,false,true,false,false);
+ t->writePOS (f, false,false,false,false,true,false);
}
ittet++;
}
@@ -855,12 +855,16 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
std::set<MEdge, Less_Edge> allEmbeddedEdges;
createAllEmbeddedEdges (gr, allEmbeddedEdges);
- // daaaaaaamn slow !!!
- // connectTets(allTets.begin(),allTets.end(),allEmbeddedFaces.empty() ? NULL : &allEmbeddedFaces);
+ if (allEmbeddedFaces.empty())
{
std::vector<faceXtet> conn;
connectTets_vector2(allTets, conn);
}
+ else
+ {
+ // daaaaaaamn slow !!!
+ connectTets(allTets.begin(),allTets.end(),&allEmbeddedFaces);
+ }
double t1 = Cpu();
std::vector<MTet4*> illegals;
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
index b4711b3f5a8bf0e36f241a9ea9b6daad17bc891d..4f5330e5bf5258f51aed56726486b3ee17c63597 100644
--- a/Mesh/meshGRegionRelocateVertex.cpp
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -27,7 +27,6 @@ static double objective_function (double xi, MVertex *ver,
else
// minQual = std::min((lt[i]->specialQuality()), minQual);
minQual = std::min(fabs(lt[i]->minSICNShapeMeasure())*.2, minQual);
-// minQual = std::min(lt[i]->minAnisotropyMeasure(), minQual);
}
ver->x() = x;
ver->y() = y;
diff --git a/Mesh/periodical.cpp b/Mesh/periodical.cpp
index 42e0c92972a70e0c81c86e7e0dcc84ebeb87a489..21833c682cbcf132e462f2bcd2c0c1133aa6a6a7 100644
--- a/Mesh/periodical.cpp
+++ b/Mesh/periodical.cpp
@@ -215,9 +215,9 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
*pointer = cell;
pointers.push_back(pointer);
generators.push_back(SPoint3(x,y,z));
- printf("%d %d (%f,%f,%f)\n",loopA.pid()+1,number+1,x,y,z);
- table.insert(std::pair<int,int>(loopA.pid(),number));
- number++;
+ // printf("%d %d (%f,%f,%f)\n",loopA.pid()+1,number+1,x,y,z);
+ table.insert(std::pair<int,int>(loopA.pid(),number));
+ number++;
}while(loopA.inc());
}
else{
@@ -229,9 +229,9 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
*pointer = cell;
pointers.push_back(pointer);
generators.push_back(SPoint3(x,y,z));
- printf("%d %d (%f,%f,%f)\n",loopB.pid()+1,number+1,x,y,z);
- table.insert(std::pair<int,int>(loopB.pid(),number));
- number++;
+ // printf("%d %d (%f,%f,%f)\n",loopB.pid()+1,number+1,x,y,z);
+ table.insert(std::pair<int,int>(loopB.pid(),number));
+ number++;
}while(loopB.inc());
}
@@ -257,7 +257,7 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
}
}
- printf("\nSquared root of smallest face area : %.9f\n\n",sqrt(min_area));
+ // printf("\nSquared root of smallest face area : %.9f\n\n",sqrt(min_area));
std::ofstream file("MicrostructurePolycrystal3D.pos");
file << "View \"test\" {\n";
@@ -407,6 +407,8 @@ int voroMetal3D::get_counter(){
}
void voroMetal3D::print_geo_point(int index,double x,double y,double z,std::ofstream& file){
+ file.precision(17);
+
file << "Point(" << index << ")={"
<< x << "," << y << "," << z
<< ",c};\n";
@@ -707,11 +709,12 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
file << "};\n";
- printf("\nNumber of exterior face periodicities : %d\n",2*count);
- printf("Total number of exterior faces : %zu\n\n",faces.size());
+ // printf("\nNumber of exterior face periodicities : %d\n",2*count);
+ // printf("Total number of exterior faces : %zu\n\n",faces.size());
std::ofstream file3;
file3.open("MicrostructurePolycrystal3D.geo",std::ios::out | std::ios::app);
+ file3.precision(17);
std::ofstream file4("MicrostructurePolycrystal3D2.pos");
file4 << "View \"test\" {\n";
@@ -733,6 +736,24 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
gf1 = pairs[i].first;
gf2 = pairs[i].second;
+ std::list<GVertex*> gv1 = gf1->vertices();
+ std::list<GVertex*> gv2 = gf2->vertices();
+
+ std::list<GVertex*>::iterator it1 = gv1.begin();
+ std::list<GVertex*>::iterator it2 = gv2.begin();
+
+ SPoint3 cg1 (0,0,0);
+ SPoint3 cg2 (0,0,0);
+
+ for (; it1 != gv1.end(); it1++,it2++){
+ cg1 += SPoint3((*it1)->x(),(*it1)->y(),(*it1)->z());
+ cg2 += SPoint3((*it2)->x(),(*it2)->y(),(*it2)->z());
+ }
+
+ SVector3 dx = (cg2-cg1) * (1./gv1.size());
+
+ // printf("%g %g %g\n",dx.x(),dx.y(),dx.z());
+
edges1 = gf1->edges();
edges2 = gf2->edges();
@@ -837,7 +858,7 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
if(indices1.size()!=indices2.size()){
printf("Error\n\n");
}
-
+ /*
file3 << "Periodic Surface " << gf1->tag() << " {";
for(j=0;j<indices1.size();j++){
@@ -852,6 +873,12 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
}
file3 << "};\n";
+ */
+
+
+ file3 << "Periodic Surface {" << gf1->tag() << " }={ " << gf2->tag() << " } Translate { " << -dx.x() << "," << -dx.y() << "," << -dx.z() << "};\n";
+
+
}
file4 << "};\n";
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index 81c82855b25d85d7248832da8b75f9528dd32ce8..97eb1da16e508001473a02a6e1d30d1e91510fe8 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -15,21 +15,162 @@
#include <iomanip>
#include "pointsGenerators.h"
#include "OS.h"
-#define TMBEG(n, what) static int I##n = 0, IND##n = _CoeffDataAnisotropy::mytimes.size(); \
-if (++I##n == 1) { \
- Msg::Info("Time%d = "#what, IND##n); \
- _CoeffDataAnisotropy::mytimes.push_back(0); \
- _CoeffDataAnisotropy::mynumbers.push_back(0); \
-} \
-_CoeffDataAnisotropy::mynumbers[IND##n] += 1; \
-double TM##n = Cpu();
-#define TMEND(n) _CoeffDataAnisotropy::mytimes[IND##n] += Cpu() - TM##n;
-namespace jacobianBasedQuality {
+static const double cTri = 2/std::sqrt(3);
+static const double cTet = std::sqrt(2);
+static const double cPyr = std::sqrt(2);
+
+static inline void computeCoeffLengthVectors_(const fullMatrix<double> &mat,
+ fullMatrix<double> &coeff,
+ int type, int numCoeff = -1)
+{
+ int sz1 = numCoeff > -1 ? numCoeff : mat.size1();
-double _CoeffDataScaledJac::cTri = 2/std::sqrt(3);
-double _CoeffDataScaledJac::cTet = std::sqrt(2);
-double _CoeffDataScaledJac::cPyr = std::sqrt(2);
+ switch (type) {
+ case TYPE_QUA: coeff.resize(sz1, 2); break;
+ case TYPE_TRI: coeff.resize(sz1, 3); break;
+ case TYPE_HEX: coeff.resize(sz1, 3); break;
+ case TYPE_PRI: coeff.resize(sz1, 4); break;
+ case TYPE_TET: coeff.resize(sz1, 6); break;
+ case TYPE_PYR: coeff.resize(sz1, 6); break;
+ default:
+ Msg::Error("Unkown type for IGE computation");
+ coeff.resize(0, 0);
+ return;
+ }
+
+ if (type != TYPE_PYR) {
+ for (int i = 0; i < sz1; i++) {
+ coeff(i, 0) = std::sqrt(pow_int(mat(i, 0), 2) +
+ pow_int(mat(i, 1), 2) +
+ pow_int(mat(i, 2), 2) );
+ coeff(i, 1) = std::sqrt(pow_int(mat(i, 3), 2) +
+ pow_int(mat(i, 4), 2) +
+ pow_int(mat(i, 5), 2) );
+ }
+ if (mat.size2() > 6) { // if 3D
+ for (int i = 0; i < sz1; i++) {
+ coeff(i, 2) = std::sqrt(pow_int(mat(i, 6), 2) +
+ pow_int(mat(i, 7), 2) +
+ pow_int(mat(i, 8), 2) );
+ }
+ }
+ else if (type == TYPE_TRI) {
+ for (int i = 0; i < sz1; i++) {
+ coeff(i, 2) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) +
+ pow_int(mat(i, 4) - mat(i, 1), 2) +
+ pow_int(mat(i, 5) - mat(i, 2), 2) );
+ }
+ }
+ switch (type) {
+ case TYPE_TET:
+ case TYPE_PRI:
+ for (int i = 0; i < sz1; i++) {
+ coeff(i, 3) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) +
+ pow_int(mat(i, 4) - mat(i, 1), 2) +
+ pow_int(mat(i, 5) - mat(i, 2), 2) );
+ }
+ }
+ if (type == TYPE_TET) {
+ for (int i = 0; i < sz1; i++) {
+ coeff(i, 4) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0), 2) +
+ pow_int(mat(i, 7) - mat(i, 1), 2) +
+ pow_int(mat(i, 8) - mat(i, 2), 2) );
+ coeff(i, 5) = std::sqrt(pow_int(mat(i, 6) - mat(i, 3), 2) +
+ pow_int(mat(i, 7) - mat(i, 4), 2) +
+ pow_int(mat(i, 8) - mat(i, 5), 2) );
+ }
+ }
+ }
+ else {
+ for (int i = 0; i < sz1; i++) {
+ coeff(i, 0) = std::sqrt(pow_int(2*mat(i, 0), 2) +
+ pow_int(2*mat(i, 1), 2) +
+ pow_int(2*mat(i, 2), 2) );
+ coeff(i, 1) = std::sqrt(pow_int(2*mat(i, 3), 2) +
+ pow_int(2*mat(i, 4), 2) +
+ pow_int(2*mat(i, 5), 2) );
+ coeff(i, 2) = std::sqrt(pow_int(mat(i, 6) + mat(i, 0) + mat(i, 3), 2) +
+ pow_int(mat(i, 7) + mat(i, 1) + mat(i, 4), 2) +
+ pow_int(mat(i, 8) + mat(i, 2) + mat(i, 5), 2) );
+ coeff(i, 3) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0) + mat(i, 3), 2) +
+ pow_int(mat(i, 7) - mat(i, 1) + mat(i, 4), 2) +
+ pow_int(mat(i, 8) - mat(i, 2) + mat(i, 5), 2) );
+ coeff(i, 4) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0) - mat(i, 3), 2) +
+ pow_int(mat(i, 7) - mat(i, 1) - mat(i, 4), 2) +
+ pow_int(mat(i, 8) - mat(i, 2) - mat(i, 5), 2) );
+ coeff(i, 5) = std::sqrt(pow_int(mat(i, 6) + mat(i, 0) - mat(i, 3), 2) +
+ pow_int(mat(i, 7) + mat(i, 1) - mat(i, 4), 2) +
+ pow_int(mat(i, 8) + mat(i, 2) - mat(i, 5), 2) );
+ }
+ }
+}
+
+static inline void computeIGE_(const fullVector<double> &det,
+ const fullMatrix<double> &v,
+ fullVector<double> &ige,
+ int type)
+{
+ int sz = std::min(det.size(), v.size1());
+ ige.resize(sz);
+
+ switch(type) {
+ case TYPE_QUA:
+ for (int i = 0; i < sz; i++) {
+ ige(i) = det(i)/v(i, 0)/v(i, 1);
+ }
+ break;
+ case TYPE_HEX:
+ for (int i = 0; i < sz; i++) {
+ ige(i) = det(i)/v(i, 0)/v(i, 1)/v(i, 2);
+ }
+ break;
+ case TYPE_TRI:
+ for (int i = 0; i < sz; i++) {
+ ige(i) = cTri * det(i) * (1/v(i,0)/v(i,1) +
+ 1/v(i,0)/v(i,2) +
+ 1/v(i,1)/v(i,2)) / 3;
+ }
+ break;
+ case TYPE_PRI:
+ for (int i = 0; i < sz; i++) {
+ ige(i) = cTri * det(i) * (1/v(i,0)/v(i,1)/v(i,2) +
+ 1/v(i,0)/v(i,3)/v(i,2) +
+ 1/v(i,1)/v(i,3)/v(i,2)) / 3;
+ }
+ break;
+ case TYPE_TET:
+ for (int i = 0; i < sz; i++) {
+ ige(i) = cTet * det(i) * (1/v(i,0)/v(i,5)/v(i,1) +
+ 1/v(i,0)/v(i,5)/v(i,2) +
+ 1/v(i,0)/v(i,5)/v(i,3) +
+ 1/v(i,0)/v(i,5)/v(i,4) +
+ 1/v(i,1)/v(i,4)/v(i,0) +
+ 1/v(i,1)/v(i,4)/v(i,2) +
+ 1/v(i,1)/v(i,4)/v(i,3) +
+ 1/v(i,1)/v(i,4)/v(i,5) +
+ 1/v(i,2)/v(i,3)/v(i,0) +
+ 1/v(i,2)/v(i,3)/v(i,1) +
+ 1/v(i,2)/v(i,3)/v(i,4) +
+ 1/v(i,2)/v(i,3)/v(i,5))/ 12;
+ }
+ break;
+ case TYPE_PYR:
+ for (int i = 0; i < sz; i++) {
+ ige(i) = cPyr * det(i) * (1/v(i,0)/v(i,1)/v(i,2) +
+ 1/v(i,0)/v(i,1)/v(i,3) +
+ 1/v(i,0)/v(i,1)/v(i,4) +
+ 1/v(i,0)/v(i,1)/v(i,5) +
+ 1/v(i,2)/v(i,3)/v(i,4) +
+ 1/v(i,2)/v(i,3)/v(i,5) +
+ 1/v(i,4)/v(i,5)/v(i,2) +
+ 1/v(i,4)/v(i,5)/v(i,3) ) / 8;
+ }
+ break;
+ }
+}
+
+namespace jacobianBasedQuality {
void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
const fullMatrix<double> *normals)
@@ -65,10 +206,12 @@ void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
}
}
-double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk)
+double minIGEMeasure(MElement *el, bool knownValid, bool reversedOk)
{
bool isReversed = false;
if (!knownValid) {
+ // Computation of the measure should never
+ // be performed to invalid elements (for which the measure is 0).
double jmin, jmax;
minMaxJacobianDeterminant(el, jmin, jmax);
if (jmax < 0) {
@@ -77,8 +220,6 @@ double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk)
}
else if (jmin <= 0) return 0;
}
- // Computation of the minimum (or maximum) of the scaled Jacobian should never
- // be performed to invalid elements (for which the measure is 0).
fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
el->getNodesCoord(nodesXYZ);
@@ -113,7 +254,7 @@ double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk)
jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse);
break;
default:
- Msg::Error("Scaled Jacobian not implemented for type of element %d", el->getType());
+ Msg::Error("IGE measure not implemented for type of element %d", el->getType());
return -1;
}
gradBasis = BasisFactory::getGradientBasis(jacMatSpace);
@@ -142,7 +283,7 @@ double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk)
std::vector<_CoeffData*> domains;
domains.push_back(
- new _CoeffDataScaledJac(coeffDetBez, coeffMatBez, jacBasis->getBezier(),
+ new _CoeffDataIGE(coeffDetBez, coeffMatBez, jacBasis->getBezier(),
gradBasis->getBezier(), 0, el->getType())
);
@@ -164,12 +305,14 @@ double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk)
return min;
}
-double minIsotropyMeasure(MElement *el,
- bool knownValid,
- bool reversedOk)
+double minICNMeasure(MElement *el,
+ bool knownValid,
+ bool reversedOk)
{
bool isReversed = false;
if (!knownValid) {
+ // Computation of the measure should never
+ // be performed to invalid elements (for which the measure is 0).
double jmin, jmax;
minMaxJacobianDeterminant(el, jmin, jmax);
if (jmax < 0) {
@@ -178,8 +321,6 @@ double minIsotropyMeasure(MElement *el,
}
else if (jmin <= 0) return 0;
}
- // Computation of the minimum (or maximum) of the scaled Jacobian should never
- // be performed to invalid elements (for which the measure is 0).
fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
el->getNodesCoord(nodesXYZ);
@@ -214,7 +355,7 @@ double minIsotropyMeasure(MElement *el,
jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse);
break;
default:
- Msg::Error("Isotropy not implemented for type of element %d", el->getType());
+ Msg::Error("ICN not implemented for type of element %d", el->getType());
return -1;
}
gradBasis = BasisFactory::getGradientBasis(jacMatSpace);
@@ -243,7 +384,7 @@ double minIsotropyMeasure(MElement *el,
std::vector<_CoeffData*> domains;
domains.push_back(
- new _CoeffDataIsotropy(coeffDetBez, coeffMatBez, jacBasis->getBezier(),
+ new _CoeffDataICN(coeffDetBez, coeffMatBez, jacBasis->getBezier(),
gradBasis->getBezier(), 0)
);
@@ -263,25 +404,59 @@ double minIsotropyMeasure(MElement *el,
return min;
}
-//double minSampledAnisotropyMeasure(MElement *el, int deg, bool writeInFile)//fordebug
-//{
-// fullMatrix<double> samplingPoints;
-// bool serendip = false;
-// gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
-//
-// double values[3];
-// double min = std::numeric_limits<double>::infinity();
-// for (int i = 0; i < samplingPoints.size1(); ++i) {
-// min = std::min(min, el->getEigenvaluesMetric(samplingPoints(i, 0),
-// samplingPoints(i, 1),
-// samplingPoints(i, 2),
-// values));
-// }
-//
-// return min;
-//}
+void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
+{
+ fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
+ el->getNodesCoord(nodesXYZ);
+
+ const bool serendipFalse = false;
+ FuncSpaceData jacMatSpace, jacDetSpace;
+
+ const int type = el->getType();
+ switch(type) {
+ case TYPE_TRI:
+ case TYPE_TET:
+ case TYPE_QUA:
+ case TYPE_HEX:
+ case TYPE_PRI:
+ jacMatSpace = FuncSpaceData(el, deg, &serendipFalse);
+ jacDetSpace = FuncSpaceData(el, deg, &serendipFalse);
+ break;
+ case TYPE_PYR:
+ jacMatSpace = FuncSpaceData(el, true, deg-1, 1, &serendipFalse);
+ jacDetSpace = FuncSpaceData(el, true, deg-1, 1, &serendipFalse);
+ break;
+ default:
+ Msg::Error("ICN not implemented for type of element %d", el->getType());
+ return;
+ }
+
+ const GradientBasis *gradBasis = BasisFactory::getGradientBasis(jacMatSpace);
+ const JacobianBasis *jacBasis = BasisFactory::getJacobianBasis(jacDetSpace);
+
+ fullVector<double> coeffDeterminant(jacBasis->getNumJacNodes());
+ jacBasis->getSignedJacobian(nodesXYZ, coeffDeterminant);
+
+ fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9);
+ gradBasis->getAllGradientsFromNodes(nodesXYZ, coeffMatLag);
+
+ fullMatrix<double> v;
+ computeCoeffLengthVectors_(coeffMatLag, v, type);
-double minSampledIsotropyMeasure(MElement *el, int deg, bool writeInFile)//fordebug
+ fullVector<double> ige;
+ computeIGE_(coeffDeterminant, v, ige, type);
+
+ min = std::numeric_limits<double>::infinity();
+ max = -min;
+ for (int i = 0; i < ige.size(); ++i) {
+ min = std::min(min, ige(i));
+ max = std::max(max, ige(i));
+ }
+
+ return;
+}
+
+double minSampledICNMeasure(MElement *el, int deg)//fordebug
{
fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
el->getNodesCoord(nodesXYZ);
@@ -310,7 +485,7 @@ double minSampledIsotropyMeasure(MElement *el, int deg, bool writeInFile)//forde
// jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse);
break;
default:
- Msg::Error("Isotropy not implemented for type of element %d", el->getType());
+ Msg::Error("ICN not implemented for type of element %d", el->getType());
return -1;
}
gradBasis = BasisFactory::getGradientBasis(jacMatSpace);
@@ -340,7 +515,7 @@ double minSampledIsotropyMeasure(MElement *el, int deg, bool writeInFile)//forde
return min;
}
-double minSampledScaledJacobian(MElement *el, int deg, bool writeInFile)//fordebug
+double minSampledIGEMeasure(MElement *el, int deg)//fordebug
{
fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
el->getNodesCoord(nodesXYZ);
@@ -369,7 +544,7 @@ double minSampledScaledJacobian(MElement *el, int deg, bool writeInFile)//fordeb
// jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse);
break;
default:
- Msg::Error("Isotropy not implemented for type of element %d", el->getType());
+ Msg::Error("ICN not implemented for type of element %d", el->getType());
return -1;
}
gradBasis = BasisFactory::getGradientBasis(jacMatSpace);
@@ -469,8 +644,8 @@ void _CoeffDataJac::getSubCoeff(std::vector<_CoeffData*> &v) const
}
}
-// Scaled Jacobian (quality of quads and hexes)
-_CoeffDataScaledJac::_CoeffDataScaledJac(fullVector<double> &det,
+// IGE measure (Inverse Gradient Error)
+_CoeffDataIGE::_CoeffDataIGE(fullVector<double> &det,
fullMatrix<double> &mat,
const bezierBasis *bfsDet,
const bezierBasis *bfsMat,
@@ -480,7 +655,7 @@ _CoeffDataScaledJac::_CoeffDataScaledJac(fullVector<double> &det,
_bfsDet(bfsDet), _bfsMat(bfsMat), _type(type)
{
if (!det.getOwnData() || !mat.getOwnData()) {
- Msg::Fatal("Cannot create an instance of _CoeffDataScaledJac from a "
+ Msg::Fatal("Cannot create an instance of _CoeffDataIGE from a "
"fullVector or a fullMatrix that does not own its data.");
}
// _coeffsJacDet and _coeffsJacMat reuse data, this avoid to allocate new
@@ -499,13 +674,13 @@ _CoeffDataScaledJac::_CoeffDataScaledJac(fullVector<double> &det,
// computation of _maxB not implemented for now
}
-bool _CoeffDataScaledJac::boundsOk(double minL, double maxL) const
+bool _CoeffDataIGE::boundsOk(double minL, double maxL) const
{
static double tol = 1e-3;
return minL - _minB < tol;
}
-void _CoeffDataScaledJac::getSubCoeff(std::vector<_CoeffData*> &v) const
+void _CoeffDataIGE::getSubCoeff(std::vector<_CoeffData*> &v) const
{
v.clear();
v.reserve(_bfsDet->getNumDivision());
@@ -522,74 +697,30 @@ void _CoeffDataScaledJac::getSubCoeff(std::vector<_CoeffData*> &v) const
fullMatrix<double> coeffM(szM1, szM2);
coeffD.copy(subCoeffD, i * szD, szD, 0);
coeffM.copy(subCoeffM, i * szM1, szM1, 0, szM2, 0, 0);
- _CoeffDataScaledJac *newData;
- newData = new _CoeffDataScaledJac(coeffD, coeffM, _bfsDet, _bfsMat,
+ _CoeffDataIGE *newData;
+ newData = new _CoeffDataIGE(coeffD, coeffM, _bfsDet, _bfsMat,
_depth+1, _type);
v.push_back(newData);
}
}
-void _CoeffDataScaledJac::_computeAtCorner(double &min, double &max) const
+void _CoeffDataIGE::_computeAtCorner(double &min, double &max) const
{
- min = std::numeric_limits<double>::infinity();
- max = -min;
-
fullMatrix<double> v;
- _getCoeffLengthVectors(v, true);
+ computeCoeffLengthVectors_(_coeffsJacMat, v, _type, _bfsDet->getNumLagCoeff());
- for (int i = 0; i < _bfsDet->getNumLagCoeff(); i++) {
- double sJ;
- switch(_type) {
- case TYPE_QUA:
- sJ = _coeffsJacDet(i)/v(i, 0)/v(i,1);
- break;
- case TYPE_HEX:
- sJ = _coeffsJacDet(i)/v(i,0)/v(i,1)/v(i,2);
- break;
- case TYPE_TRI:
- sJ = cTri * _coeffsJacDet(i) * (1/v(i,0)/v(i,1) +
- 1/v(i,0)/v(i,2) +
- 1/v(i,1)/v(i,2)) / 3;
- break;
- case TYPE_TET:
- sJ = cTet * _coeffsJacDet(i) * (1/v(i,0)/v(i,5)/v(i,1) +
- 1/v(i,0)/v(i,5)/v(i,2) +
- 1/v(i,0)/v(i,5)/v(i,3) +
- 1/v(i,0)/v(i,5)/v(i,4) +
- 1/v(i,1)/v(i,4)/v(i,0) +
- 1/v(i,1)/v(i,4)/v(i,2) +
- 1/v(i,1)/v(i,4)/v(i,3) +
- 1/v(i,1)/v(i,4)/v(i,5) +
- 1/v(i,2)/v(i,3)/v(i,0) +
- 1/v(i,2)/v(i,3)/v(i,1) +
- 1/v(i,2)/v(i,3)/v(i,4) +
- 1/v(i,2)/v(i,3)/v(i,5)) / 12;
- break;
- case TYPE_PRI:
- sJ = cTri * _coeffsJacDet(i) * (1/v(i,0)/v(i,1)/v(i,2) +
- 1/v(i,0)/v(i,3)/v(i,2) +
- 1/v(i,1)/v(i,3)/v(i,2)) / 3;
- break;
- case TYPE_PYR:
- sJ = cPyr * _coeffsJacDet(i) * (1/v(i,0)/v(i,1)/v(i,2) +
- 1/v(i,0)/v(i,1)/v(i,3) +
- 1/v(i,0)/v(i,1)/v(i,4) +
- 1/v(i,0)/v(i,1)/v(i,5) +
- 1/v(i,2)/v(i,3)/v(i,4) +
- 1/v(i,2)/v(i,3)/v(i,5) +
- 1/v(i,4)/v(i,5)/v(i,2) +
- 1/v(i,4)/v(i,5)/v(i,3) ) / 8;
- break;
- default:
- Msg::Error("Unkown type for scaled jac computation");
- return;
- }
- min = std::min(min, sJ);
- max = std::max(max, sJ);
+ fullVector<double> ige;
+ computeIGE_(_coeffsJacDet, v, ige, _type);
+
+ min = std::numeric_limits<double>::infinity();
+ max = -min;
+ for (int i = 0; i < ige.size(); ++i) {
+ min = std::min(min, ige(i));
+ max = std::max(max, ige(i));
}
}
-double _CoeffDataScaledJac::_computeLowerBound() const
+double _CoeffDataIGE::_computeLowerBound() const
{
// Speedup: If one coeff _coeffsJacDet is negative, without bounding
// J^2/(a^2+b^2), we would get with certainty a negative lower bound.
@@ -601,7 +732,7 @@ double _CoeffDataScaledJac::_computeLowerBound() const
}
fullMatrix<double> v;
- _getCoeffLengthVectors(v, false);
+ computeCoeffLengthVectors_(_coeffsJacMat, v, _type);
fullVector<double> prox[6];
for (int i = 0; i < v.size2(); ++i) {
@@ -706,12 +837,12 @@ double _CoeffDataScaledJac::_computeLowerBound() const
}
default:
- Msg::Info("Unknown type for scaled Jacobian (%d)", _type);
+ Msg::Info("Unknown type for IGE (%d)", _type);
return -1;
}
}
-void _CoeffDataScaledJac::_getCoeffLengthVectors(fullMatrix<double> &coeff,
+void _CoeffDataIGE::_getCoeffLengthVectors(fullMatrix<double> &coeff,
bool corners) const
{
int sz1 = corners ? _bfsDet->getNumLagCoeff() : _coeffsJacMat.size1();
@@ -724,7 +855,7 @@ void _CoeffDataScaledJac::_getCoeffLengthVectors(fullMatrix<double> &coeff,
case TYPE_TET: coeff.resize(sz1, 6); break;
case TYPE_PYR: coeff.resize(sz1, 6); break;
default:
- Msg::Error("Unkown type for scaled jac computation");
+ Msg::Error("Unkown type for IGE computation");
coeff.resize(0, 0);
return;
}
@@ -798,8 +929,8 @@ void _CoeffDataScaledJac::_getCoeffLengthVectors(fullMatrix<double> &coeff,
}
}
-// Isotropy measure
-_CoeffDataIsotropy::_CoeffDataIsotropy(fullVector<double> &det,
+// ICN measure (Inverse Condition Number)
+_CoeffDataICN::_CoeffDataICN(fullVector<double> &det,
fullMatrix<double> &mat,
const bezierBasis *bfsDet,
const bezierBasis *bfsMat,
@@ -809,7 +940,7 @@ _CoeffDataIsotropy::_CoeffDataIsotropy(fullVector<double> &det,
_bfsDet(bfsDet), _bfsMat(bfsMat)
{
if (!det.getOwnData() || !mat.getOwnData()) {
- Msg::Fatal("Cannot create an instance of _CoeffDataScaledJac from a "
+ Msg::Fatal("Cannot create an instance of _CoeffDataIGE from a "
"fullVector or a fullMatrix that does not own its data.");
}
// _coeffsJacDet and _coeffsMetric reuse data, this avoid to allocate new
@@ -832,13 +963,13 @@ _CoeffDataIsotropy::_CoeffDataIsotropy(fullVector<double> &det,
// _maxB not used for now
}
-bool _CoeffDataIsotropy::boundsOk(double minL, double maxL) const
+bool _CoeffDataICN::boundsOk(double minL, double maxL) const
{
static double tol = 1e-3;
return minL < tol*1e-3 || (_minB > minL-tol && _minB > minL*(1-100*tol));
}
-void _CoeffDataIsotropy::getSubCoeff(std::vector<_CoeffData*> &v) const
+void _CoeffDataICN::getSubCoeff(std::vector<_CoeffData*> &v) const
{
v.clear();
v.reserve(_bfsMat->getNumDivision());
@@ -855,13 +986,13 @@ void _CoeffDataIsotropy::getSubCoeff(std::vector<_CoeffData*> &v) const
fullMatrix<double> coeffM(szM1, szM2);
coeffD.copy(subCoeffD, i * szD, szD, 0);
coeffM.copy(subCoeffM, i * szM1, szM1, 0, szM2, 0, 0);
- _CoeffDataIsotropy *newData
- = new _CoeffDataIsotropy(coeffD, coeffM, _bfsDet, _bfsMat, _depth+1);
+ _CoeffDataICN *newData
+ = new _CoeffDataICN(coeffD, coeffM, _bfsDet, _bfsMat, _depth+1);
v.push_back(newData);
}
}
-void _CoeffDataIsotropy::_computeAtCorner(double &min, double &max) const
+void _CoeffDataICN::_computeAtCorner(double &min, double &max) const
{
min = std::numeric_limits<double>::infinity();
max = -min;
@@ -881,7 +1012,7 @@ void _CoeffDataIsotropy::_computeAtCorner(double &min, double &max) const
}
}
-double _CoeffDataIsotropy::_computeLowerBound() const
+double _CoeffDataICN::_computeLowerBound() const
{
// Speedup: If one coeff _coeffsJacDet is negative, we would get
// a negative lower bound. For now, returning 0.
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index 2ed23400530756962ad0bc52b2c3ff5d2875acc7..6c3b62992a1df10b80eed5c812d02207a1c11103 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -17,25 +17,22 @@ namespace jacobianBasedQuality {
void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
const fullMatrix<double> *normals = NULL);
-double minScaledJacobian(MElement *el,
- bool knownValid = false,
- bool reversedOk = false);
-double minIsotropyMeasure(MElement *el,
- bool knownValid = false,
- bool reversedOk = false);
-//double minSampledAnisotropyMeasure(MElement *el, int order,//fordebug
-// bool writeInFile = false);
-double minSampledIsotropyMeasure(MElement *el, int order,//fordebug
- bool writeInFile = false);
-double minSampledScaledJacobian(MElement *el, int order,//fordebug
- bool writeInFile = false);
+double minIGEMeasure(MElement *el,
+ bool knownValid = false,
+ bool reversedOk = false);
+double minICNMeasure(MElement *el,
+ bool knownValid = false,
+ bool reversedOk = false);
+void sampleIGEMeasure(MElement *el, int order, double &min, double &max);
+double minSampledICNMeasure(MElement *el, int order);//fordebug
+double minSampledIGEMeasure(MElement *el, int order);//fordebug
class _CoeffData
{
protected:
double _minL, _maxL; //Extremum of Jac at corners
double _minB, _maxB; //Extremum of bezier coefficients
- int _depth;
+ const int _depth;
public:
_CoeffData(int depth) : _minL(0), _maxL(0), _minB(0), _maxB(0),
@@ -75,24 +72,21 @@ public:
int getNumMeasure() const {return 1;}//fordebug
};
-class _CoeffDataScaledJac: public _CoeffData
+class _CoeffDataIGE: public _CoeffData
{
private:
const fullVector<double> _coeffsJacDet;
const fullMatrix<double> _coeffsJacMat;
const bezierBasis *_bfsDet, *_bfsMat;
- int _type;
- static double cTri;
- static double cTet;
- static double cPyr;
+ const int _type;
public:
- _CoeffDataScaledJac(fullVector<double> &det,
+ _CoeffDataIGE(fullVector<double> &det,
fullMatrix<double> &mat,
const bezierBasis *bfsDet,
const bezierBasis *bfsMat,
int depth, int type);
- ~_CoeffDataScaledJac() {}
+ ~_CoeffDataIGE() {}
bool boundsOk(double minL, double maxL) const;
void getSubCoeff(std::vector<_CoeffData*>&) const;
@@ -102,11 +96,9 @@ private:
void _computeAtCorner(double &min, double &max) const;
double _computeLowerBound() const;
void _getCoeffLengthVectors(fullMatrix<double> &, bool corners = false) const;
- void _getCoeffScaledJacobian(const fullMatrix<double> &coeffLengthVectors,
- fullMatrix<double> &coeffScaledJacobian) const;
};
-class _CoeffDataIsotropy: public _CoeffData
+class _CoeffDataICN: public _CoeffData
{
private:
const fullVector<double> _coeffsJacDet;
@@ -114,12 +106,12 @@ private:
const bezierBasis *_bfsDet, *_bfsMat;
public:
- _CoeffDataIsotropy(fullVector<double> &det,
+ _CoeffDataICN(fullVector<double> &det,
fullMatrix<double> &metric,
const bezierBasis *bfsDet,
const bezierBasis *bfsMet,
int depth);
- ~_CoeffDataIsotropy() {}
+ ~_CoeffDataICN() {}
bool boundsOk(double minL, double maxL) const;
void getSubCoeff(std::vector<_CoeffData*>&) const;
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 983fbf4695d3714b58685e0fed39829ab81f8ff3..12e5f2c7462b74d784f42036e84b390b5ab35c67 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -1477,7 +1477,7 @@ double Recombinator::min_scaled_jacobian(Hex &hex) {
/*MHexahedron *h1 = new MHexahedron(a, b, c, d, e, f, g, h);
MHexahedron *h2 = new MHexahedron(e, f, g, h, a, b, c, d);
double min1 = jacobianBasedQuality::minScaledJacobian(h1);
- double min2 = jacobianBasedQuality::minScaledJacobian(h2);
+ double min2 = jacobianBasedQuality::minIGEMeasure(h2);
for(i=0;i<8;i++){
file << jacobians[i] << " ";
}
@@ -2980,7 +2980,7 @@ double Supplementary::min_scaled_jacobian(Prism prism) {
/*MPrism *p1 = new MPrism(a, b, c, d, e, f);
MPrism *p2 = new MPrism(d, e, f, a, b, c);
double min1 = jacobianBasedQuality::minScaledJacobian(p1);
- double min2 = jacobianBasedQuality::minScaledJacobian(p2);
+ double min2 = jacobianBasedQuality::minIGEMeasure(p2);
for(i=0;i<6;i++){
file << jacobians[i] << " ";
}
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index 9025a5029bffbaf44f91a7643337c9f6b596f601..c1420bfe260aa1c3fde3e378751d55a5810c4073 100644
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -40,7 +40,7 @@ class Hex {
void compute_quality()
{
MHexahedron elt(vertices_);
- quality = jacobianBasedQuality::minScaledJacobian(&elt, false, true);
+ quality = jacobianBasedQuality::minIGEMeasure(&elt, false, true);
}
void initialize()
{
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index c85cafc660ac259413b928f439d81849f564dc1e..3a49ddb3afd1954d9d541c94c68006cdbd2d75d2 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -24,8 +24,8 @@ class bezierBasis;
StringXNumber CurvedMeshOptions_Number[] = {
{GMSH_FULLRC, "Jacobian determinant", NULL, 1},
- {GMSH_FULLRC, "Scaled Jacobian", NULL, 0},
- {GMSH_FULLRC, "Isotropy", NULL, 1},
+ {GMSH_FULLRC, "IGE measure", NULL, 1},
+ {GMSH_FULLRC, "ICN measure", NULL, 1},
{GMSH_FULLRC, "Hidding threshold", NULL, 9},
{GMSH_FULLRC, "Draw PView", NULL, 0},
{GMSH_FULLRC, "Recompute", NULL, 0},
@@ -54,36 +54,41 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
{
return "Plugin(AnalyseCurvedMesh) analyse all elements of a given dimension. "
"According to what is asked, it computes the minimum of the Jacobian "
- "determinant (J), of the scaled Jacobian and/or of the isotropy measure. "
- "Statistics are printed and if asked a Pview is created for each measure. "
+ "determinant (J), the IGE quality measure (Inverse Gradient Error) and/or "
+ "the ICN quality measure (Inverse Condition Number). "
+ "Statistics are printed and, if asked, a Pview is created for each measure. "
"The plugin hides elements for which the measure mu > 'Hidding threshold', "
- "where mu is the isotropy measure if asked otherwise the scaled Jacobian if "
+ "where mu is the ICN measure if asked otherwise the IGE measure if "
"asked otherwise the Jacobian determinant.\n"
"\n"
- "J is faster to compute but gives informations only on validity while the "
- "other measure gives also informations on quality.\n"
- "Warning: the scaled Jacobian is experimental for triangles, tetrahedra, "
- "prisms and pyramids. Computation may take a lot of time for those "
- "elements!\n"
+ "J is faster to compute but gives information only on validity while the "
+ "other measure gives also information on quality.\n"
+ "The IGE measure is related to the error on the gradient of the finite "
+ "element solution. It is the scaled Jacobian for quads and hexes and a new "
+ "measure for triangles and tetrahedra.\n"
+ "The ICN measure is related to the condition number of the stiffness matrix.\n"
+ "(See article \"Efficient computation of the minimum of shape quality "
+ "measures on curvilinear finite elements\" for details.)\n"
"\n"
"Parameters:\n"
"\n"
"- Jacobian determinant = {0, 1}\n"
- "- Scaled Jacobian = {0, 1}\n"
- "- Isotropy = {0, 1}\n"
"\n"
- "- Hidding threshold = [0, 1]: Does nothing if Isotropy == 0 and Scaled "
- "Jacobian == 0. Otherwise, hides all element for which min(mu) is "
- "strictly greater than the threshold, where mu is the isotropy if Isotropy "
- "== 1, otherwise it is the Scaled Jacobian. If threshold == 1, no effect, "
- "if == 0 hide all elements except invalid.\n"
+ "- IGE measure = {0, 1}\n"
"\n"
- "- Draw PView = {0, 1}: Creates a PView of min(J)/max(J), min(scaled Jac) "
- "and/or min(isotropy) according to what is asked. If 'Recompute' = 1, a "
- "new PView is redrawed.\n"
+ "- ICN measure = {0, 1}\n"
"\n"
- "- Recompute = {0,1}: If the mesh has changed, set to 1 to recompute the "
- "bounds.\n"
+ "- Hidding threshold = [0, 1]: Hides all element for which min(mu) is "
+ "strictly greater than the threshold, where mu is the ICN if ICN measure == 1, "
+ "otherwise mu is the IGE it IGE measure == 1, "
+ "otherwise mu is the Jacobian determinant.\n"
+ "If threshold == 0, hides all elements except invalid.\n"
+ "\n"
+ "- Draw PView = {0, 1}: Creates a PView of min(J)/max(J), min(IGE) "
+ "and/or min(ICN) according to what is asked. If 'Recompute' = 1, "
+ "new PViews are created.\n"
+ "\n"
+ "- Recompute = {0,1}: Should be 1 if the mesh has changed.\n"
"\n"
"- Dimension = {-1, 1, 2, 3, 4}: If == -1, analyse element of the greater "
"dimension. If == 4, analyse 2D and 3D elements.";
@@ -92,33 +97,33 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
{
_m = GModel::current();
- bool computeJac = static_cast<int>(CurvedMeshOptions_Number[0].def);
- bool computeScaledJac = static_cast<int>(CurvedMeshOptions_Number[1].def);
- bool computeIsotropy = static_cast<int>(CurvedMeshOptions_Number[2].def);
- _threshold = static_cast<double>(CurvedMeshOptions_Number[3].def);
- bool drawPView = static_cast<int>(CurvedMeshOptions_Number[4].def);
- bool recompute = static_cast<bool>(CurvedMeshOptions_Number[5].def);
- int askedDim = static_cast<int>(CurvedMeshOptions_Number[6].def);
+ bool computeJac = static_cast<int>(CurvedMeshOptions_Number[0].def);
+ bool computeIGE = static_cast<int>(CurvedMeshOptions_Number[1].def);
+ bool computeICN = static_cast<int>(CurvedMeshOptions_Number[2].def);
+ _threshold = static_cast<double>(CurvedMeshOptions_Number[3].def);
+ bool drawPView = static_cast<int>(CurvedMeshOptions_Number[4].def);
+ bool recompute = static_cast<bool>(CurvedMeshOptions_Number[5].def);
+ int askedDim = static_cast<int>(CurvedMeshOptions_Number[6].def);
if (askedDim < 0 || askedDim > 4) askedDim = _m->getDim();
if (recompute) {
_data.clear();
if (askedDim < 4) {
- _computedJ[askedDim-1] = false;
- _computedS[askedDim-1] = false;
- _computedI[askedDim-1] = false;
- _PViewJ[askedDim-1] = false;
- _PViewS[askedDim-1] = false;
- _PViewI[askedDim-1] = false;
+ _computedJac[askedDim-1] = false;
+ _computedIGE[askedDim-1] = false;
+ _computedICN[askedDim-1] = false;
+ _PViewJac[askedDim-1] = false;
+ _PViewIGE[askedDim-1] = false;
+ _PViewICN[askedDim-1] = false;
}
else {
- _computedJ[1] = _computedJ[2] = false;
- _computedS[1] = _computedS[2] = false;
- _computedI[1] = _computedI[2] = false;
- _PViewJ[1] = _PViewJ[2] = false;
- _PViewS[1] = _PViewS[2] = false;
- _PViewI[1] = _PViewI[2] = false;
+ _computedJac[1] = _computedJac[2] = false;
+ _computedIGE[1] = _computedIGE[2] = false;
+ _computedICN[1] = _computedICN[2] = false;
+ _PViewJac[1] = _PViewJac[2] = false;
+ _PViewIGE[1] = _PViewIGE[2] = false;
+ _PViewICN[1] = _PViewICN[2] = false;
}
}
@@ -128,7 +133,7 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
bool printStatI = false;
for (int dim = 1; dim <= 3 && dim <= _m->getDim(); ++dim) {
if ((askedDim == 4 && dim > 1) || dim == askedDim) {
- if (!_computedJ[dim-1]) {
+ if (!_computedJac[dim-1]) {
double time = Cpu();
Msg::StatusBar(true, "Computing Jacobian for %dD elements...", dim);
_computeMinMaxJandValidity(dim);
@@ -136,34 +141,34 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
dim, Cpu()-time);
printStatJ = true;
}
- if (computeScaledJac && !_computedS[dim-1]) {
+ if (computeIGE && !_computedIGE[dim-1]) {
double time = Cpu();
- Msg::StatusBar(true, "Computing scaled Jacobian for %dD elements...", dim);
- _computeMinScaledJac(dim);
- Msg::StatusBar(true, "Done computing scaled Jacobian for %dD elements (%g s)",
+ Msg::StatusBar(true, "Computing IGE for %dD elements...", dim);
+ _computeMinIGE(dim);
+ Msg::StatusBar(true, "Done computing IGE for %dD elements (%g s)",
dim, Cpu()-time);
printStatS = true;
}
- if (computeIsotropy && !_computedI[dim-1]) {
+ if (computeICN && !_computedICN[dim-1]) {
double time = Cpu();
- Msg::StatusBar(true, "Computing isotropy for %dD elements...", dim);
- _computeMinIsotropy(dim);
- Msg::StatusBar(true, "Done computing isotropy for %dD elements (%g s)",
+ Msg::StatusBar(true, "Computing ICN for %dD elements...", dim);
+ _computeMinICN(dim);
+ Msg::StatusBar(true, "Done computing ICN for %dD elements (%g s)",
dim, Cpu()-time);
printStatI = true;
}
}
}
if (printStatJ) _printStatJacobian();
- if (printStatS) _printStatScaledJac();
- if (printStatI) _printStatIsotropy();
+ if (printStatS) _printStatIGE();
+ if (printStatI) _printStatICN();
// Create PView
if (drawPView)
for (int dim = 1; dim <= 3; ++dim) {
if ((askedDim == 4 && dim > 1) || dim == askedDim) {
- if (!_PViewJ[dim-1] && computeJac) {
- _PViewJ[dim-1] = true;
+ if (!_PViewJac[dim-1] && computeJac) {
+ _PViewJac[dim-1] = true;
std::map<int, std::vector<double> > dataPV;
for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element();
@@ -176,8 +181,8 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
new PView(name.str().c_str(), "ElementData", _m, dataPV);
}
}
- if (!_PViewS[dim-1] && computeScaledJac) {
- _PViewS[dim-1] = true;
+ if (!_PViewIGE[dim-1] && computeIGE) {
+ _PViewIGE[dim-1] = true;
std::map<int, std::vector<double> > dataPV;
for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element();
@@ -186,12 +191,12 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
}
if (dataPV.size()) {
std::stringstream name;
- name << "Scaled Jacobian " << dim << "D";
+ name << "IGE " << dim << "D";
new PView(name.str().c_str(), "ElementData", _m, dataPV);
}
}
- if (!_PViewI[dim-1] && computeIsotropy) {
- _PViewI[dim-1] = true;
+ if (!_PViewICN[dim-1] && computeICN) {
+ _PViewICN[dim-1] = true;
std::map<int, std::vector<double> > dataPV;
for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element();
@@ -200,7 +205,7 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
}
if (dataPV.size()) {
std::stringstream name;
- name << "Isotropy " << dim << "D";
+ name << "ICN " << dim << "D";
new PView(name.str().c_str(), "ElementData", _m, dataPV);
}
}
@@ -210,8 +215,8 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
// Hide elements
#if defined(HAVE_OPENGL)
_hideWithThreshold(askedDim,
- computeIsotropy ? 2 :
- computeScaledJac ? 1 :
+ computeICN ? 2 :
+ computeIGE ? 1 :
computeJac ? 0 : -1);
CTX::instance()->mesh.changed = (ENT_ALL);
drawContext::global()->draw();
@@ -222,7 +227,7 @@ PView* GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
{
- if (_computedJ[dim-1]) return;
+ if (_computedJac[dim-1]) return;
std::set<GEntity*, GEntityLessThan> entities;
switch (dim) {
@@ -337,12 +342,12 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
}
delete normals;
}
- _computedJ[dim-1] = true;
+ _computedJac[dim-1] = true;
}
-void GMSH_AnalyseCurvedMeshPlugin::_computeMinScaledJac(int dim)
+void GMSH_AnalyseCurvedMeshPlugin::_computeMinIGE(int dim)
{
- if (_computedS[dim-1]) return;
+ if (_computedIGE[dim-1]) return;
MsgProgressStatus progress(_data.size());
@@ -353,17 +358,17 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinScaledJac(int dim)
_data[i].setMinS(0);
}
else {
- _data[i].setMinS(jacobianBasedQuality::minScaledJacobian(el, true));
+ _data[i].setMinS(jacobianBasedQuality::minIGEMeasure(el, true));
}
progress.next();
}
- _computedS[dim-1] = true;
+ _computedIGE[dim-1] = true;
}
-void GMSH_AnalyseCurvedMeshPlugin::_computeMinIsotropy(int dim)
+void GMSH_AnalyseCurvedMeshPlugin::_computeMinICN(int dim)
{
- if (_computedI[dim-1]) return;
+ if (_computedICN[dim-1]) return;
MsgProgressStatus progress(_data.size());
@@ -374,12 +379,12 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinIsotropy(int dim)
_data[i].setMinI(0);
}
else {
- _data[i].setMinI(jacobianBasedQuality::minIsotropyMeasure(el, true));
+ _data[i].setMinI(jacobianBasedQuality::minICNMeasure(el, true));
}
progress.next();
}
- _computedI[dim-1] = true;
+ _computedICN[dim-1] = true;
}
int GMSH_AnalyseCurvedMeshPlugin::_hideWithThreshold(int askedDim, int whichMeasure)
@@ -451,17 +456,17 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian()
avgratJ /= count;
avgratJc /= countc;
- Msg::Info("Min Jac. det. (minJ): %6.3g, %6.3g, %6.3g (= min, avg, max)",
+ Msg::Info("minJ = %8.3g, %8.3g, %8.3g (min, avg, max)",
infminJ, avgminJ, supminJ);
- Msg::Info("Ratio minJ/maxJ : %6.3f, %6.3f, %6.3f (= worst, avg, best)",
- infratJ, avgratJ, supratJ);
if (countc && countc < count)
- Msg::Info(" (avg = %.3f"
- " on the %d non-constant elements)",
+ Msg::Info("minJ/maxJ = %8.3f (avg on the %d "
+ "non-constant elements)",
avgratJc, countc);
+ Msg::Info("minJ/maxJ = %8.3f, %8.3f, %8.3f (worst, avg, best)",
+ infratJ, avgratJ, supratJ);
}
-void GMSH_AnalyseCurvedMeshPlugin::_printStatScaledJac()
+void GMSH_AnalyseCurvedMeshPlugin::_printStatIGE()
{
if (_data.empty()) {
Msg::Info("No stat to print.");
@@ -477,11 +482,11 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatScaledJac()
}
avgminS /= _data.size();
- Msg::Info("Scaled Jacobian : %6.3f, %6.3f, %6.3f (= worst, avg, best)",
- infminS, avgminS, supminS);
+ Msg::Info("IGE = %8.3f, %8.3f, %8.3f (worst, avg, best)",
+ infminS, avgminS, supminS);
}
-void GMSH_AnalyseCurvedMeshPlugin::_printStatIsotropy()
+void GMSH_AnalyseCurvedMeshPlugin::_printStatICN()
{
if (_data.empty()) {
Msg::Info("No stat to print.");
@@ -497,8 +502,8 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatIsotropy()
}
avgminI /= _data.size();
- Msg::Info("Isotropy : %6.3f, %6.3f, %6.3f (= worst, avg, best)",
- infminI, avgminI, supminI);
+ Msg::Info("ICN = %8.3f, %8.3f, %8.3f (worst, avg, best)",
+ infminI, avgminI, supminI);
}
#endif
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index 17aa641503509de110cc41d4c77580666cfcd786..d1aac032a5306c5db3438f49739f82779a8e130d 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -19,21 +19,21 @@ class data_elementMinMax
{
private:
MElement *_el;
- double _minJ, _maxJ, _minS, _minI;
+ double _minJac, _maxJac, _minIGE, _minICN;
public:
data_elementMinMax(MElement *e,
double minJ = 2,
double maxJ = 0,
- double minS = -1,
- double minI = -1)
- : _el(e), _minJ(minJ), _maxJ(maxJ), _minS(minS), _minI(minI) {}
- void setMinS(double r) { _minS = r; }
- void setMinI(double r) { _minI = r; }
+ double minIGE = -1,
+ double minICN = -1)
+ : _el(e), _minJac(minJ), _maxJac(maxJ), _minIGE(minIGE), _minICN(minICN) {}
+ void setMinS(double r) { _minIGE = r; }
+ void setMinI(double r) { _minICN = r; }
MElement* element() { return _el; }
- double minJ() { return _minJ; }
- double maxJ() { return _maxJ; }
- double minS() { return _minS; }
- double minI() { return _minI; }
+ double minJ() { return _minJac; }
+ double maxJ() { return _maxJac; }
+ double minS() { return _minIGE; }
+ double minI() { return _minICN; }
};
class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
@@ -43,8 +43,8 @@ private :
double _threshold;
// for 1d, 2d, 3d
- bool _computedJ[3], _computedS[3], _computedI[3];
- bool _PViewJ[3], _PViewS[3], _PViewI[3];
+ bool _computedJac[3], _computedIGE[3], _computedICN[3];
+ bool _PViewJac[3], _PViewIGE[3], _PViewICN[3];
std::vector<data_elementMinMax> _data;
@@ -54,18 +54,18 @@ public :
_m = NULL;
_threshold = -1;
for (int i = 0; i < 3; ++i) {
- _computedJ[i] = false;
- _computedS[i] = false;
- _computedI[i] = false;
- _PViewJ[i] = false;
- _PViewS[i] = false;
- _PViewI[i] = false;
+ _computedJac[i] = false;
+ _computedIGE[i] = false;
+ _computedICN[i] = false;
+ _PViewJac[i] = false;
+ _PViewIGE[i] = false;
+ _PViewICN [i] = false;
}
}
std::string getName() const { return "AnalyseCurvedMesh"; }
std::string getShortHelp() const
{
- return "Compute bounds on Jacobian and metric quality.";
+ return "Compute validity and quality of curved elements.";
}
std::string getHelp() const;
std::string getAuthor() const { return "Amaury Johnen"; }
@@ -75,12 +75,12 @@ public :
private :
void _computeMinMaxJandValidity(int dim);
- void _computeMinScaledJac(int dim);
- void _computeMinIsotropy(int dim);
+ void _computeMinIGE(int dim);
+ void _computeMinICN(int dim);
int _hideWithThreshold(int askedDim, int whichMeasure);
void _printStatJacobian();
- void _printStatScaledJac();
- void _printStatIsotropy();
+ void _printStatIGE();
+ void _printStatICN();
};
#endif
diff --git a/contrib/mobile/Android/app/src/main/AndroidManifest.xml b/contrib/mobile/Android/app/src/main/AndroidManifest.xml
index 463f90e4aad52c8ea88a63397068cec9dceefcca..5ed362e26a2101a858f3ca92505f340f2802e61f 100644
--- a/contrib/mobile/Android/app/src/main/AndroidManifest.xml
+++ b/contrib/mobile/Android/app/src/main/AndroidManifest.xml
@@ -1,7 +1,7 @@
<manifest xmlns:android="http://schemas.android.com/apk/res/android"
package="org.geuz.onelab"
android:versionName="2.0.0"
- android:versionCode="32"
+ android:versionCode="36"
android:installLocation="auto" >
<uses-sdk android:minSdkVersion="14"
diff --git a/contrib/mobile/androidUtils.cpp b/contrib/mobile/androidUtils.cpp
index 6620401bd53168c927de4a9f39ed26f91899947c..53bdb54bc22a987115660d39b2aa8d7f0407dfc1 100644
--- a/contrib/mobile/androidUtils.cpp
+++ b/contrib/mobile/androidUtils.cpp
@@ -48,10 +48,10 @@ public:
!gCallbackObject || (gJavaVM->AttachCurrentThread(&env, NULL)) < 0) return;
jstring jstr = env->NewStringUTF(message.c_str());
jclass jClass = env->FindClass("org/geuz/onelab/Gmsh");
- if(jClass == 0)
+ if(jClass == 0 || env->ExceptionCheck())
return;
jmethodID mid = env->GetMethodID(jClass, "ShowPopup", "(Ljava/lang/String;)V");
- if (mid == 0)
+ if (mid == 0 || env->ExceptionCheck())
return;
env->CallVoidMethod(gCallbackObject, mid, jstr);
env->DeleteLocalRef(jstr);
@@ -64,10 +64,10 @@ public:
!gCallbackObject || (gJavaVM->AttachCurrentThread(&env, NULL)) < 0) return;
jstring jstr = env->NewStringUTF(message.c_str());
jclass jClass = env->FindClass("org/geuz/onelab/Gmsh");
- if(jClass == 0)
+ if(jClass == 0 || env->ExceptionCheck())
return;
jmethodID mid = env->GetMethodID(jClass, "SetLoading", "(Ljava/lang/String;)V");
- if (mid == 0)
+ if (mid == 0 || env->ExceptionCheck())
return;
env->CallVoidMethod(gCallbackObject, mid, jstr);
env->DeleteLocalRef(jstr);
@@ -90,10 +90,10 @@ void requestRender()
if(gJavaVM->GetEnv(reinterpret_cast<void**>(&env), JNI_VERSION_1_6) != JNI_OK ||
!gCallbackObject || (gJavaVM->AttachCurrentThread(&env, NULL)) < 0) return;
jclass jClass = env->FindClass("org/geuz/onelab/Gmsh");
- if(jClass == 0)
+ if(jClass == 0 || env->ExceptionCheck())
return;
jmethodID mid = env->GetMethodID(jClass, "RequestRender", "()V");
- if (mid == 0)
+ if (mid == 0 || env->ExceptionCheck())
return;
env->CallVoidMethod(gCallbackObject, mid);
env->DeleteLocalRef(jClass);
@@ -106,19 +106,27 @@ void getBitmapFromString(const char *text, int textsize, unsigned char **map,
if(gJavaVM->GetEnv(reinterpret_cast<void**>(&env), JNI_VERSION_1_6) != JNI_OK ||
!gCallbackObject || (gJavaVM->AttachCurrentThread(&env, NULL)) < 0) return;
jclass jClass = env->FindClass("org/geuz/onelab/StringTexture");
- if(jClass == 0)
+ if(jClass == 0 || env->ExceptionCheck())
return;
jstring jtext = env->NewStringUTF(text);
jmethodID mid = env->GetStaticMethodID(jClass, "getHeightFromString",
"(Ljava/lang/String;I)I");
+ if(mid == 0 || env->ExceptionCheck())
+ return;
*height = env->CallStaticIntMethod(jClass, mid, jtext, textsize);
mid = env->GetStaticMethodID(jClass, "getWidthFromString", "(Ljava/lang/String;I)I");
+ if(mid == 0 || env->ExceptionCheck())
+ return;
*width = env->CallStaticIntMethod(jClass, mid, jtext, textsize);
if(realWidth != NULL){
mid = env->GetStaticMethodID(jClass, "getRealWidthFromString", "(Ljava/lang/String;I)I");
+ if(mid == 0 || env->ExceptionCheck())
+ return;
*realWidth = env->CallStaticIntMethod(jClass, mid, jtext, textsize);
}
mid = env->GetStaticMethodID(jClass, "getBytesFromString", "(Ljava/lang/String;I)[B");
+ if(mid == 0 || env->ExceptionCheck())
+ return;
jobject jbuffer = env->CallStaticObjectMethod(jClass, mid, jtext, textsize);
jbyteArray *jarray = reinterpret_cast<jbyteArray*>(&jbuffer);
int ms = (*height) * (*width);
diff --git a/doc/gmsh.html b/doc/gmsh.html
index d3a24e0137b5d6a3a820016ff537ec4eb1f4e297..1daa362300f77e4e1fc9afcd1cf7ebfdceec02d3 100644
--- a/doc/gmsh.html
+++ b/doc/gmsh.html
@@ -99,13 +99,13 @@ Public License (GPL)</a>:
<ul>
<li>
<p class="highlight">
- Current stable release (version 3.0.2, May 13 2017):
- <a href="bin/Windows/gmsh-3.0.2-Windows64.zip">Windows</a>
- (<a href="bin/Windows/gmsh-3.0.2-Windows32.zip">32 bit</a>),
- <a href="bin/Linux/gmsh-3.0.2-Linux64.tgz">Linux</a>,
- <a href="bin/MacOSX/gmsh-3.0.2-MacOSX.dmg">MacOS</a>
+ Current stable release (version 3.0.3, June 23 2017):
+ <a href="bin/Windows/gmsh-3.0.3-Windows64.zip">Windows</a>
+ (<a href="bin/Windows/gmsh-3.0.3-Windows32.zip">32 bit</a>),
+ <a href="bin/Linux/gmsh-3.0.3-Linux64.tgz">Linux</a>,
+ <a href="bin/MacOSX/gmsh-3.0.3-MacOSX.dmg">MacOS</a>
and
- <a href="src/gmsh-3.0.2-source.tgz">source code</a>
+ <a href="src/gmsh-3.0.3-source.tgz">source code</a>
</p>
<p>
<em>A <a href="doc/texinfo/gmsh.html#Tutorial"><strong>tutorial</strong></a>
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index d9902b9dee73dbc5ed2f9eb32b47e3cf5ea1acbb..07ce0734c8c89567e65e787c36e73823ac4d62d1 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -5483,13 +5483,15 @@ coefficients A,B,C,D of the equation A*x+B*y+C*y+D=0, which can be
adjusted interactively by dragging the mouse in the input
fields.
-@item What is the signification of Rho, Eta and Gamma in Tools->Statistics?
+@item What is the signification of SICN, Gamma and SIGE in Tools->Statistics?
They measure the quality of the tetrahedra in a mesh:
+SICN ~ signed inverse condition number
+
Gamma ~ inscribed_radius / circumscribed_radius
-Eta ~ volume^(2/3) / sum_edge_length^2
-Rho ~ min_edge_length / max_edge_length
+
+SIGE ~ signed inverse error on the gradient of FE solution
For the exact definitions, see Geo/MElement.cpp. The graphs plot the
the number of elements vs. the quality measure.
diff --git a/doc/texinfo/opt_geometry.texi b/doc/texinfo/opt_geometry.texi
index 4b9c983a425d66a68c0cd650d9e959d3daa1e807..d571f0bdb67797d5bd0364dcfb1bee3241ae0f5d 100644
--- a/doc/texinfo/opt_geometry.texi
+++ b/doc/texinfo/opt_geometry.texi
@@ -251,7 +251,7 @@ Saved in: @code{General.OptionsFileName}
@item Geometry.SurfaceType
Surface display type (0=cross, 1=wireframe, 2=solid)@*
-Default value: @code{2}@*
+Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}
@item Geometry.Tangents
diff --git a/doc/texinfo/opt_mesh.texi b/doc/texinfo/opt_mesh.texi
index 6e43a4da4b5659bac09688708090ea2c02af2838..a7965733d333f2d9b5783de8e82b6586c95bd830 100644
--- a/doc/texinfo/opt_mesh.texi
+++ b/doc/texinfo/opt_mesh.texi
@@ -500,7 +500,7 @@ Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}
@item Mesh.QualityType
-Type of quality measure (0=gamma~vol/sum_face/max_edge, 1=eta~vol^(2/3)/sum_edge^2, 2=rho~min_edge/max_edge)@*
+Type of quality measure (0=SICN~signed inverse condition number, 1=SIGE~signed inverse gradient error, 2=gamma~vol/sum_face/max_edge, 3=Disto~minJ/maxJ@*
Default value: @code{2}@*
Saved in: @code{General.OptionsFileName}
diff --git a/doc/texinfo/opt_plugin.texi b/doc/texinfo/opt_plugin.texi
index 9ad6a01c619369d16b9b3879a4fd360ef5607088..05f989ce6710ea7f634a5ed1a228eac0715adcda 100644
--- a/doc/texinfo/opt_plugin.texi
+++ b/doc/texinfo/opt_plugin.texi
@@ -5,31 +5,36 @@
@ftable @code
@item Plugin(AnalyseCurvedMesh)
-Plugin(AnalyseCurvedMesh) analyse all elements of a given dimension. According to what is asked, it computes the minimum of the Jacobian determinant (J), of the scaled Jacobian and/or of the isotropy measure. Statistics are printed and if asked a Pview is created for each measure. The plugin hides elements for which the measure mu > 'Hidding threshold', where mu is the isotropy measure if asked otherwise the scaled Jacobian if asked otherwise the Jacobian determinant.@*
+Plugin(AnalyseCurvedMesh) analyse all elements of a given dimension. According to what is asked, it computes the minimum of the Jacobian determinant (J), the IGE quality measure (Inverse Gradient Error) and/or the ICN quality measure (Inverse Condition Number). Statistics are printed and, if asked, a Pview is created for each measure. The plugin hides elements for which the measure mu > 'Hidding threshold', where mu is the ICN measure if asked otherwise the IGE measure if asked otherwise the Jacobian determinant.@*
@*
-J is faster to compute but gives informations only on validity while the other measure gives also informations on quality.@*
-Warning: the scaled Jacobian is experimental for triangles, tetrahedra, prisms and pyramids. Computation may take a lot of time for those elements!@*
+J is faster to compute but gives information only on validity while the other measure gives also information on quality.@*
+The IGE measure is related to the error on the gradient of the finite element solution. It is the scaled Jacobian for quads and hexes and a new measure for triangles and tetrahedra.@*
+The ICN measure is related to the condition number of the stiffness matrix.@*
+(See article "Efficient computation of the minimum of shape quality measures on curvilinear finite elements" for details.)@*
@*
Parameters:@*
@*
- Jacobian determinant = @{0, 1@}@*
-- Scaled Jacobian = @{0, 1@}@*
-- Isotropy = @{0, 1@}@*
@*
-- Hidding threshold = [0, 1]: Does nothing if Isotropy == 0 and Scaled Jacobian == 0. Otherwise, hides all element for which min(mu) is strictly greater than the threshold, where mu is the isotropy if Isotropy == 1, otherwise it is the Scaled Jacobian. If threshold == 1, no effect, if == 0 hide all elements except invalid.@*
+- IGE measure = @{0, 1@}@*
@*
-- Draw PView = @{0, 1@}: Creates a PView of min(J)/max(J), min(scaled Jac) and/or min(isotropy) according to what is asked. If 'Recompute' = 1, a new PView is redrawed.@*
+- ICN measure = @{0, 1@}@*
@*
-- Recompute = @{0,1@}: If the mesh has changed, set to 1 to recompute the bounds.@*
+- Hidding threshold = [0, 1]: Hides all element for which min(mu) is strictly greater than the threshold, where mu is the ICN if ICN measure == 1, otherwise mu is the IGE it IGE measure == 1, otherwise mu is the Jacobian determinant.@*
+If threshold == 0, hides all elements except invalid.@*
+@*
+- Draw PView = @{0, 1@}: Creates a PView of min(J)/max(J), min(IGE) and/or min(ICN) according to what is asked. If 'Recompute' = 1, new PViews are created.@*
+@*
+- Recompute = @{0,1@}: Should be 1 if the mesh has changed.@*
@*
- Dimension = @{-1, 1, 2, 3, 4@}: If == -1, analyse element of the greater dimension. If == 4, analyse 2D and 3D elements.
Numeric options:
@table @code
@item Jacobian determinant
Default value: @code{1}
-@item Scaled Jacobian
-Default value: @code{0}
-@item Isotropy
+@item IGE measure
+Default value: @code{1}
+@item ICN measure
Default value: @code{1}
@item Hidding threshold
Default value: @code{9}
diff --git a/doc/texinfo/opt_print.texi b/doc/texinfo/opt_print.texi
index b9e51319f6415f16e0c43a19eae57c07281b93ed..88ba43f8221ac98f3ea5d616a501d3d9bf8f36f1 100644
--- a/doc/texinfo/opt_print.texi
+++ b/doc/texinfo/opt_print.texi
@@ -164,8 +164,13 @@ Save Eta quality measure in mesh statistics exported as post-processing views@*
Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}
-@item Print.PostRho
-Save Rho quality measure in mesh statistics exported as post-processing views@*
+@item Print.PostSICN
+Save SICN (signed inverse condition number) quality measure in mesh statistics exported as post-processing views@*
+Default value: @code{0}@*
+Saved in: @code{General.OptionsFileName}
+
+@item Print.PostSICN
+Save SIGE (signed inverse gradient error) quality measure in mesh statistics exported as post-processing views@*
Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}
diff --git a/utils/misc/gource.sh b/utils/misc/gource.sh
index 83967f90b66a0e88f214e1f0dbe4b8407dc9ad31..6c9d541c9ea5e73ffb739ac0185040dd7c9a74e9 100755
--- a/utils/misc/gource.sh
+++ b/utils/misc/gource.sh
@@ -3,3 +3,6 @@
### -f --start-date '2006-01-01'
gource -f --time-scale 2 -b 000000 --seconds-per-day 0.0001 --hide filenames --file-filter projects --file-filter data --file-filter branches --file-filter benchmarks_private --file-filter benchmarks_kst --highlight-user geuzaine -1280x720 -o - | ffmpeg -y -r 60 -f image2pipe -vcodec ppm -i - -vcodec libx264 -preset ultrafast -pix_fmt yuv420p -crf 1 -threads 0 -bf 0 gource.mp4
+
+# resample to make it 4 x faster
+ffmpeg -i gource.mp4 -filter:v "setpts=0.25*PTS" gource_x4speed.mp4