From ec1f8f27c994f718c96b51afe9df8d8e053292d8 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 28 Nov 2007 14:18:10 +0000
Subject: [PATCH] 3D Optimization is now implemented !
---
Common/Context.h | 2 +-
Common/DefaultOptions.h | 2 +
Common/Options.cpp | 13 ++-
Common/Options.h | 1 +
Fltk/Callbacks.cpp | 17 ++-
Fltk/Callbacks.h | 1 +
Fltk/GUI.cpp | 17 +--
Mesh/Generator.cpp | 19 +++-
Mesh/Generator.h | 1 +
Mesh/meshGRegion.cpp | 19 +++-
Mesh/meshGRegion.h | 10 +-
Mesh/meshGRegionDelaunayInsertion.cpp | 144 +++++++++++---------------
Mesh/meshGRegionDelaunayInsertion.h | 18 +++-
Mesh/meshGRegionLocalMeshMod.cpp | 50 +++++----
Mesh/meshGRegionLocalMeshMod.h | 3 +-
15 files changed, 189 insertions(+), 128 deletions(-)
diff --git a/Common/Context.h b/Common/Context.h
index 66ca89f9d2..da82dddd23 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -178,7 +178,7 @@ public :
double label_frequency;
int point_type; // flat or 3D
double point_size, line_width;
- int optimize, refine_steps;
+ int optimize, optimizeNetgen, refine_steps;
int quality_type, label_type;
double quality_inf, quality_sup, radius_inf, radius_sup;
double scaling_factor, lc_factor, rand_factor, lc_integration_precision;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 8b10e3d11f..dccec3c36f 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -981,6 +981,8 @@ StringXNumber MeshOptions_Number[] = {
"Display size of normal vectors (in pixels)" },
{ F|O, "Optimize" , opt_mesh_optimize , 0. ,
+ "Optimize the mesh to improve the quality of tetrahedral elements" },
+ { F|O, "OptimizeNetgen" , opt_mesh_optimize_netgen , 0. ,
"Optimize the mesh using Netgen to improve the quality of tetrahedral elements" },
{ F|O, "Points" , opt_mesh_points , 0. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 7af769f1e0..f90161a6a1 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.368 2007-11-26 14:34:09 remacle Exp $
+// $Id: Options.cpp,v 1.369 2007-11-28 14:18:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -4209,6 +4209,17 @@ double opt_mesh_optimize(OPT_ARGS_NUM)
return CTX.mesh.optimize;
}
+double opt_mesh_optimize_netgen(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.mesh.optimize =(int) val;
+#if defined(HAVE_FLTK)
+ if(WID && (action & GMSH_GUI))
+ WID->mesh_butt[2]->value(CTX.mesh.optimizeNetgen);
+#endif
+ return CTX.mesh.optimize;
+}
+
double opt_mesh_refine_steps(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 864b3aa168..aac39194c4 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -426,6 +426,7 @@ double opt_geometry_snap1(OPT_ARGS_NUM);
double opt_geometry_snap2(OPT_ARGS_NUM);
double opt_mesh_label_frequency(OPT_ARGS_NUM);
double opt_mesh_optimize(OPT_ARGS_NUM);
+double opt_mesh_optimize_netgen(OPT_ARGS_NUM);
double opt_mesh_refine_steps(OPT_ARGS_NUM);
double opt_mesh_normals(OPT_ARGS_NUM);
double opt_mesh_tangents(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 2bc63a6bf1..3e64e74a13 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.552 2007-10-03 19:40:40 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.553 2007-11-28 14:18:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -1071,6 +1071,7 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
opt_mesh_reverse_all_normals(0, GMSH_SET, WID->mesh_butt[0]->value());
opt_mesh_lc_from_curvature(0, GMSH_SET, WID->mesh_butt[1]->value());
opt_mesh_optimize(0, GMSH_SET, WID->mesh_butt[2]->value());
+ opt_mesh_optimize_netgen(0, GMSH_SET, WID->mesh_butt[24]->value());
opt_mesh_order(0, GMSH_SET, WID->mesh_value[3]->value());
opt_mesh_smooth_internal_edges(0, GMSH_SET, WID->mesh_butt[3]->value());
opt_mesh_second_order_incomplete(0, GMSH_SET, WID->mesh_butt[4]->value());
@@ -3860,6 +3861,20 @@ void mesh_optimize_cb(CALLBACK_ARGS)
Msg(STATUS2N, " ");
}
+void mesh_optimize_netgen_cb(CALLBACK_ARGS)
+{
+ if(CTX.threads_lock) {
+ Msg(INFO, "I'm busy! Ask me that later...");
+ return;
+ }
+ CTX.threads_lock = 1;
+ OptimizeMeshNetgen(GModel::current());
+ CTX.threads_lock = 0;
+ CTX.mesh.changed = ENT_LINE | ENT_SURFACE | ENT_VOLUME;
+ Draw();
+ Msg(STATUS2N, " ");
+}
+
void mesh_define_length_cb(CALLBACK_ARGS)
{
diff --git a/Fltk/Callbacks.h b/Fltk/Callbacks.h
index 49b4b14fd7..72ff2afc10 100644
--- a/Fltk/Callbacks.h
+++ b/Fltk/Callbacks.h
@@ -282,6 +282,7 @@ void mesh_update_edges_cb(CALLBACK_ARGS);
void mesh_parameterize_cb(CALLBACK_ARGS);
void mesh_degree_cb(CALLBACK_ARGS);
void mesh_optimize_cb(CALLBACK_ARGS);
+void mesh_optimize_netgen_cb(CALLBACK_ARGS);
void mesh_classify_cb(CALLBACK_ARGS);
void mesh_define_length_cb (CALLBACK_ARGS);
void mesh_define_recombine_cb (CALLBACK_ARGS);
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index e2c30e0844..c4282222c5 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.644 2007-11-08 19:30:30 geuzaine Exp $
+// $Id: GUI.cpp,v 1.645 2007-11-28 14:18:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -288,8 +288,9 @@ Context_Item menu_mesh[] = {
{"3D", (Fl_Callback *)mesh_3d_cb} ,
{"First order", (Fl_Callback *)mesh_degree_cb, (void*)1 } ,
{"Second order", (Fl_Callback *)mesh_degree_cb, (void*)2 } ,
+ {"Optimize", (Fl_Callback *)mesh_optimize_cb} ,
#if defined(HAVE_NETGEN)
- {"Optimize quality", (Fl_Callback *)mesh_optimize_cb} ,
+ {"Optimize (netgen)", (Fl_Callback *)mesh_optimize_netgen_cb} ,
#endif
#if defined(HAVE_FOURIER_MODEL)
{"Reparameterize", (Fl_Callback *)mesh_parameterize_cb} ,
@@ -2457,16 +2458,20 @@ void GUI::create_option_window()
mesh_butt[2] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 3 * BH, BW, BH, "Optimize quality of tetrahedra");
mesh_butt[2]->type(FL_TOGGLE_BUTTON);
+ mesh_butt[2]->callback(mesh_options_ok_cb);
+
+ mesh_butt[24] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 4 * BH, BW, BH, "Optimize quality of tetrahedra with Netgen");
+ mesh_butt[24]->type(FL_TOGGLE_BUTTON);
#if !defined(HAVE_NETGEN)
- mesh_butt[2]->deactivate();
+ mesh_butt[24]->deactivate();
#endif
- mesh_butt[2]->callback(mesh_options_ok_cb);
+ mesh_butt[24]->callback(mesh_options_ok_cb);
- mesh_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 4 * BH, BW, BH, "Optimize high order mesh (only for 2D-plane)");
+ mesh_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 5 * BH, BW, BH, "Optimize high order mesh (only for 2D-plane)");
mesh_butt[3]->type(FL_TOGGLE_BUTTON);
mesh_butt[3]->callback(mesh_options_ok_cb);
- mesh_butt[21] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 5 * BH, BW, BH, "Impose C1 continuity (only for 2D-plane, order 2 and 3)");
+ mesh_butt[21] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 6 * BH, BW, BH, "Impose C1 continuity (only for 2D-plane, order 2 and 3)");
mesh_butt[21]->type(FL_TOGGLE_BUTTON);
mesh_butt[21]->callback(mesh_options_ok_cb);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index f8e6968612..a64bac22a7 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.126 2007-11-27 16:45:27 geuzaine Exp $
+// $Id: Generator.cpp,v 1.127 2007-11-28 14:18:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -307,12 +307,24 @@ void Mesh3D(GModel *m)
Msg(STATUS1, "Mesh");
}
+void OptimizeMeshNetgen(GModel *m)
+{
+ Msg(STATUS1, "Optimizing 3D with Netgen...");
+ double t1 = Cpu();
+
+ std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionNetgen());
+
+ double t2 = Cpu();
+ Msg(INFO, "Mesh 3D optimization with Netgen complete (%g s)", t2 - t1);
+ Msg(STATUS1, "Mesh");
+}
+
void OptimizeMesh(GModel *m)
{
Msg(STATUS1, "Optimizing 3D...");
double t1 = Cpu();
- std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegion());
+ std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionGmsh());
double t2 = Cpu();
Msg(INFO, "Mesh 3D optimization complete (%g s)", t2 - t1);
@@ -359,6 +371,9 @@ void GenerateMesh(int ask)
// Optimize quality
if(m->getMeshStatus() == 3 && CTX.mesh.optimize)
OptimizeMesh(m);
+ // Optimize quality with netgen
+ if(m->getMeshStatus() == 3 && CTX.mesh.optimizeNetgen)
+ OptimizeMeshNetgen(m);
// Create high order elements
if(m->getMeshStatus() && CTX.mesh.order > 1)
diff --git a/Mesh/Generator.h b/Mesh/Generator.h
index ad54f320e1..a39d30ff9e 100644
--- a/Mesh/Generator.h
+++ b/Mesh/Generator.h
@@ -25,5 +25,6 @@ class GModel;
void GetStatistics(double stat[50], double quality[3][100]=0);
void GenerateMesh(int dimension);
void OptimizeMesh(GModel *m);
+void OptimizeMeshNetgen(GModel *m);
#endif
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 552f8545e4..0635f2abdd 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.36 2007-11-12 10:41:52 colignon Exp $
+// $Id: meshGRegion.cpp,v 1.37 2007-11-28 14:18:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -565,7 +565,7 @@ void meshGRegion::operator() (GRegion *gr)
}
}
-void optimizeMeshGRegion::operator() (GRegion *gr)
+void optimizeMeshGRegionNetgen::operator() (GRegion *gr)
{
if(gr->geomType() == GEntity::DiscreteVolume) return;
@@ -594,3 +594,18 @@ void optimizeMeshGRegion::operator() (GRegion *gr)
Ng_Exit();
#endif
}
+
+void optimizeMeshGRegionGmsh::operator() (GRegion *gr)
+{
+ if(gr->geomType() == GEntity::DiscreteVolume) return;
+
+ // don't optimize extruded meshes
+ ExtrudeParams *ep = gr->meshAttributes.extrude;
+ if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY) return;
+
+ Msg(STATUS2, "Optimizing volume %d", gr->tag());
+ // import mesh into netgen, including volume tets
+
+ gmshOptimizeMesh (gr, QMTET_2);
+
+}
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index cc9081fb47..6744f4bf6f 100644
--- a/Mesh/meshGRegion.h
+++ b/Mesh/meshGRegion.h
@@ -38,8 +38,14 @@ class meshGRegionExtruded {
void operator () (GRegion *);
};
-// Optimize the mesh of the region
-class optimizeMeshGRegion {
+// Optimize the mesh of the region using gmsh's algo
+class optimizeMeshGRegionGmsh {
+ public :
+ void operator () (GRegion *);
+};
+
+// Optimize the mesh of the region using netgen's algo
+class optimizeMeshGRegionNetgen {
public :
void operator () (GRegion *);
};
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 3a13696d0e..16c8064b82 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.23 2007-11-28 11:47:36 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.24 2007-11-28 14:18:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -432,32 +432,53 @@ void recur_classify ( MTet4 *t ,
}
}
-template <class CONTAINER, class DATA>
-void gmshOptimizeMesh (CONTAINER &allTets, DATA &vSizes)
+//template <class CONTAINER, class DATA>
+void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
{
+ typedef std::list<MTet4 *> CONTAINER ;
+ CONTAINER allTets;
+ for (int i=0;i<gr->tetrahedra.size();i++){
+ allTets.push_back(new MTet4(gr->tetrahedra[i],qm));
+ }
+ gr->tetrahedra.clear();
+
+ connectTets(allTets.begin(),allTets.end());
+
+
+
+
double t1 = Cpu();
std::vector<MTet4*> illegals;
const int nbRanges=10;
int quality_ranges [nbRanges];
-
-
{
double totalVolumeb = 0.0;
double worst = 1.0;
double avg = 0;
int count=0;
- for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
{
if (!(*it)->isDeleted()){
- double vol;
- double qual = qmTet((*it)->tet(),QMTET_2,&vol);
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
worst = std::min(qual,worst);
avg+=qual;
count++;
totalVolumeb+=vol;
+ for (int i=0;i<nbRanges;i++){
+ double low = (double)i/(nbRanges);
+ double high = (double)(i+1)/(nbRanges);
+ if (qual >= low && qual < high)quality_ranges[i]++;
+ }
}
}
Msg(INFO,"Opti : START with %12.5E QBAD %12.5E QAVG %12.5E",totalVolumeb,worst,avg/count);
+ for (int i=0;i<nbRanges;i++){
+ double low = (double)i/(nbRanges);
+ double high = (double)(i+1)/(nbRanges);
+ Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
+ }
}
double qMin = 0.5;
@@ -466,35 +487,29 @@ void gmshOptimizeMesh (CONTAINER &allTets, DATA &vSizes)
int nbESwap=0, nbFSwap=0, nbReloc=0;
while (1){
- std::vector<MTet4*> newTets;
-
- // get a bad tet and try to swap each of its edges and faces
-
- for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
+ std::vector<MTet4*> newTets;
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
if (!(*it)->isDeleted()){
double vol;
- double qq = qmTet((*it)->tet(),QMTET_2,&vol);
+ double qq = (*it)->getQuality();
if (qq < qMin){
for (int i=0;i<4;i++){
- if (gmshFaceSwap(newTets,*it,i,QMTET_2)){nbFSwap++;break;}
+ if (gmshFaceSwap(newTets,*it,i,qm)){nbFSwap++;break;}
}
}
}
}
- connectTets(allTets.begin(),allTets.end());
-
illegals.clear();
for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
- for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
{
if (!(*it)->isDeleted()){
- double vol;
- double qq = qmTet((*it)->tet(),QMTET_2,&vol);
+ double qq = (*it)->getQuality();
if (qq < qMin)
for (int i=0;i<6;i++){
- if (gmshEdgeSwap(newTets,*it,i,QMTET_2)) {nbESwap++;break;}
+ if (gmshEdgeSwap(newTets,*it,i,qm)) {nbESwap++;break;}
}
if (!(*it)->isDeleted()){
if (qq < sliverLimit)illegals.push_back(*it);
@@ -504,29 +519,15 @@ void gmshOptimizeMesh (CONTAINER &allTets, DATA &vSizes)
if (qq >= low && qq < high)quality_ranges[i]++;
}
}
- // if (newTets.size())break;
}
}
- // if no new tet is created, leave
-
if (!newTets.size())break;
// add all the new tets in the container
for (int i=0;i<newTets.size();i++){
if (!newTets[i]->isDeleted()){
- // it was bugged here because the setup fct removes
- // the neighbors. Now it seems to work fine
- MTet4 *t0 = newTets[i]->getNeigh(0);
- MTet4 *t1 = newTets[i]->getNeigh(1);
- MTet4 *t2 = newTets[i]->getNeigh(2);
- MTet4 *t3 = newTets[i]->getNeigh(3);
- newTets[i]->setup(newTets[i]->tet(),vSizes);
- newTets[i]->setNeigh(0,t0);
- newTets[i]->setNeigh(1,t1);
- newTets[i]->setNeigh(2,t2);
- newTets[i]->setNeigh(3,t3);
- allTets.insert(newTets[i]);
+ allTets.push_back(newTets[i]);
}
else{
delete newTets[i]->tet();
@@ -534,56 +535,14 @@ void gmshOptimizeMesh (CONTAINER &allTets, DATA &vSizes)
}
}
-// int nbNeigh = 0;
-// int k =0;
-// MTet4 *test [10000][4];
-// for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
-// {
-// if (!(*it)->isDeleted()){
-// for (int i=0;i<4;i++){
-// MTet4 *neigh = (*it)->getNeigh(i);
-// test[k][i] = neigh;
-// if (neigh){
-// nbNeigh++;
-// }
-// }
-// k++;
-// }
-// }
-
-// printf("nbNeigh = %d --> ",nbNeigh);
-
- // connectTets(allTets.begin(),allTets.end());
-
-// nbNeigh = 0;
-// k = 0;
-// for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
-// {
-// if (!(*it)->isDeleted()){
-// for (int i=0;i<4;i++){
-// MTet4 *neigh = (*it)->getNeigh(i);
-// if (test[k][i] != neigh)
-// {
-// printf("connexion tet %d %p , item %d differs %p %p\n",k,(*it),i,neigh,test[k][i]);
-// }
-// if (neigh){
-// nbNeigh++;
-// }
-// }
-// k++;
-// }
-// }
-
-// printf("nbNeigh = %d\n",nbNeigh);
// relocate vertices
- for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
{
if (!(*it)->isDeleted()){
- double vol;
- double qq = qmTet((*it)->tet(),QMTET_2,&vol);
- if (qq < .5)
+ double qq = (*it)->getQuality();
+ if (qq < qMin)
for (int i=0;i<4;i++){
- if (gmshSmoothVertex(*it,i))nbReloc++;
+ if (gmshSmoothVertex(*it,i,qm))nbReloc++;
}
}
}
@@ -592,11 +551,11 @@ void gmshOptimizeMesh (CONTAINER &allTets, DATA &vSizes)
double worst = 1.0;
double avg = 0;
int count=0;
- for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
{
if (!(*it)->isDeleted()){
- double vol;
- double qual = qmTet((*it)->tet(),QMTET_2,&vol);
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
worst = std::min(qual,worst);
avg+=qual;
count++;
@@ -623,6 +582,19 @@ void gmshOptimizeMesh (CONTAINER &allTets, DATA &vSizes)
double high = (double)(i+1)/(nbRanges);
Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
}
+
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->isDeleted()){
+ gr->tetrahedra.push_back((*it)->tet());
+ delete *it;
+ }
+ else{
+ delete (*it)->tet();
+ delete *it;
+ }
+ }
+
}
@@ -760,7 +732,7 @@ void insertVerticesInRegion (GRegion *gr)
Msg(INFO,"cleaning up the memory %d -> %d ",n1,allTets.size());
}
}
- gmshOptimizeMesh (allTets,vSizes);
+
while (1)
{
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 846385298a..d10f0c91de 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -21,6 +21,7 @@
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "MElement.h"
+#include "qualityMeasures.h"
#include <list>
#include <set>
#include <stack>
@@ -62,10 +63,16 @@ class MTet4
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
}
- MTet4 (MTetrahedron * t) : deleted(false), gr(0),circum_radius(0.0),base(t)
+ MTet4 (MTetrahedron * t, double qual) : deleted(false), gr(0),circum_radius(qual),base(t)
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
}
+ MTet4 (MTetrahedron * t, const gmshQualityMeasure4Tet &qm) : deleted(false), gr(0),base(t)
+ {
+ neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
+ double vol;
+ circum_radius = qmTet(t,qm,&vol);
+ }
void setup ( MTetrahedron * t, std::vector<double> & sizes)
{
@@ -91,10 +98,9 @@ class MTet4
bool isDeleted () const {return deleted;}
void forceRadius (double r){circum_radius=r;}
- inline double getRadius ()const {return circum_radius;}
-
-
-
+ inline double getRadius ()const {return circum_radius;}
+ inline double getQuality ()const {return circum_radius;}
+ inline double setQuality (const double &q){circum_radius=q;}
inline MTetrahedron * tet() const {return base;}
inline MTetrahedron * &tet() {return base;}
inline void setNeigh (int iN , MTet4 *n) {neigh[iN]=n;}
@@ -225,4 +231,6 @@ private:
};
+void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm);
+
#endif
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index 448be11299..e0017a26b5 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -203,15 +203,14 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
bool closed = gmshBuildEdgeCavity ( tet,iLocalEdge,&v1,&v2,cavity,outside,ring);
-
if (!closed)return false;
double volumeRef = 0.0;
double tetQualityRef = 1;
for (int i=0;i<cavity.size();i++){
- double vol;
- tetQualityRef = std::min(qmTet (cavity[i]->tet() , cr , &vol) , tetQualityRef);
- volumeRef += fabs(vol);
+ double vol = fabs(cavity[i]->tet()->getVolume());
+ tetQualityRef = std::min(tetQualityRef,cavity[i]->getQuality());
+ volumeRef += vol;
}
// build swap patterns
@@ -288,8 +287,8 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
pv2,
pv1,
v2);
- MTet4 *t41 = new MTet4 ( tr1 );
- MTet4 *t42 = new MTet4 ( tr2 );
+ MTet4 *t41 = new MTet4 ( tr1 , tetQuality1[iT]);
+ MTet4 *t42 = new MTet4 ( tr2 , tetQuality2[iT]);
t41->setOnWhat(cavity[0]->onWhat());
t42->setOnWhat(cavity[0]->onWhat());
outside.push_back(t41);
@@ -346,10 +345,10 @@ bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
// printf("%p %p -- %p %p %p\n",v1,v2,f1,f2,f3);
- double vol1;
- double q1 = qmTet(t1->tet(),cr,&vol1);
- double vol2;
- double q2 = qmTet(t2->tet(),cr,&vol2);
+ double vol1 = fabs(t1->tet()->getVolume());
+ double q1 = t1->getQuality();
+ double vol2 = fabs(t2->tet()->getVolume());
+ double q2 = t2->getQuality();
double vol3;
double q3 = qmTet(f1,f2,v1,v2,cr,&vol3);
@@ -392,9 +391,9 @@ bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
MTetrahedron *tr1 = new MTetrahedron ( f1,f2,v1,v2);
MTetrahedron *tr2 = new MTetrahedron ( f2,f3,v1,v2);
MTetrahedron *tr3 = new MTetrahedron ( f3,f1,v1,v2);
- MTet4 *t41 = new MTet4 ( tr1 );
- MTet4 *t42 = new MTet4 ( tr2 );
- MTet4 *t43 = new MTet4 ( tr3 );
+ MTet4 *t41 = new MTet4 ( tr1 , q3);
+ MTet4 *t42 = new MTet4 ( tr2 , q4);
+ MTet4 *t43 = new MTet4 ( tr3 , q5);
t41->setOnWhat(t1->onWhat());
t42->setOnWhat(t1->onWhat());
t43->setOnWhat(t1->onWhat());
@@ -446,7 +445,8 @@ void gmshBuildVertexCavity_recur ( MTet4 *t,
bool gmshSmoothVertex ( MTet4 *t,
- int iVertex){
+ int iVertex,
+ const gmshQualityMeasure4Tet &cr){
if(t->isDeleted())throw;
if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3)return false;
@@ -455,14 +455,13 @@ bool gmshSmoothVertex ( MTet4 *t,
cavity.push_back(t);
gmshBuildVertexCavity_recur (t,t->tet()->getVertex(iVertex),cavity);
- // printf("cavity o size %d\n",cavity.size());
-
double xcg=0,ycg=0,zcg=0;
double vTot=0;
double worst = 1.0;
+
for (int i=0 ; i< cavity.size() ; i++){
- double volume;
- double q = qmTet(cavity[i]->tet(),QMTET_2,&volume);
+ double volume = fabs(cavity[i]->tet()->getVolume());
+ double q = cavity[i]->getQuality();
worst = std::min(worst,q);
xcg += 0.25*(cavity[i]->tet()->getVertex(0)->x()+
cavity[i]->tet()->getVertex(1)->x()+
@@ -492,11 +491,13 @@ bool gmshSmoothVertex ( MTet4 *t,
t->tet()->getVertex(iVertex)->y() = ycg;
t->tet()->getVertex(iVertex)->z() = zcg;
double worstAfter = 1.0;
+ double newQuals[2000];
+ if (cavity.size() > 2000) throw;
for (int i=0 ; i< cavity.size() ; i++){
double volume;
- double q = qmTet(cavity[i]->tet(),QMTET_2,&volume);
+ newQuals[i] = qmTet(cavity[i]->tet(),cr,&volume);
volumeAfter += volume;
- worstAfter =std::min(worstAfter,q);
+ worstAfter =std::min(worstAfter,newQuals[i]);
}
if (fabs(volumeAfter-vTot) > 1.e-10 * vTot || worstAfter < worst){
@@ -505,7 +506,14 @@ bool gmshSmoothVertex ( MTet4 *t,
t->tet()->getVertex(iVertex)->z() = z;
return false;
}
- return true;
+ else{
+ // restore new quality
+ for (int i=0 ; i< cavity.size() ; i++){
+ cavity[i]->setQuality(newQuals[i]);
+ }
+
+ return true;
+ }
}
diff --git a/Mesh/meshGRegionLocalMeshMod.h b/Mesh/meshGRegionLocalMeshMod.h
index cad25536b3..7a4ef0dfc1 100644
--- a/Mesh/meshGRegionLocalMeshMod.h
+++ b/Mesh/meshGRegionLocalMeshMod.h
@@ -32,7 +32,8 @@ bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
int iLocalFace,
const gmshQualityMeasure4Tet &cr);
bool gmshSmoothVertex ( MTet4 *t,
- int iLocalVertex);
+ int iLocalVertex,
+ const gmshQualityMeasure4Tet &cr);
#endif
--
GitLab