diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index defc56a3edacc85b28fc298bb74d29b2007e240f..0d0a618a6607c52f2abae9da1a663038969bed92 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.107 2007-11-04 21:03:16 remacle Exp $
+// $Id: CommandLine.cpp,v 1.108 2007-11-26 14:34:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -410,6 +410,30 @@ void Get_Options(int argc, char *argv[])
i++;
opt_mesh_c1(0, GMSH_SET, 1);
}
+ else if(!strcmp(argv[i] + 1, "statreport")) {
+ i++;
+ CTX.create_append_statreport = 1;
+ if(argv[i] != NULL) {
+ strcpy(CTX.statreport,argv[i]);
+ i++;
+ }
+ else {
+ fprintf(stderr, ERROR_STR "Missing argument\n");
+ exit(1);
+ }
+ }
+ else if(!strcmp(argv[i] + 1, "append_statreport")) {
+ i++;
+ CTX.create_append_statreport = 2;
+ if(argv[i] != NULL) {
+ strcpy(CTX.statreport,argv[i]);
+ i++;
+ }
+ else {
+ fprintf(stderr, ERROR_STR "Missing argument\n");
+ exit(1);
+ }
+ }
else if(!strcmp(argv[i] + 1, "optimize_hom")) {
i++;
opt_mesh_smooth_internal_edges(0, GMSH_SET, 1);
diff --git a/Common/Context.h b/Common/Context.h
index ca3d3d55cc83a277adce3ffb400382f8f3329296..66ca89f9d25e037be04feb7366ab55831051c74b 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -46,6 +46,12 @@ public :
char *error_filename, error_filename_fullpath[256];
// the name of the error file
+ // gmsh enables to output some mesh stats in a test suite file
+ // mesh algorithms are very sensitive to small parameter modifications
+ // and a test suite is indeed useful ;-)
+ char statreport[256];
+ int create_append_statreport; // do nothing 0 create 1 append 2
+
int session_save, options_save; // save session/option file on exit
int confirm_overwrite; // confirm overwrite when file->save as
char *display; // forced display host:0.0 under X11
diff --git a/Common/Options.cpp b/Common/Options.cpp
index ff4aa64cf9fafa2f8d0311910b8d5b6ac1b79671..7af769f1e0c65766227d8b6fdcaabdebeae3d928 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.367 2007-11-08 19:30:30 geuzaine Exp $
+// $Id: Options.cpp,v 1.368 2007-11-26 14:34:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -106,6 +106,9 @@ void Init_Options(int num)
else
strcpy(CTX.home_dir, "");
+ // By defaults, no stat report
+ CTX.create_append_statreport = 0;
+
int len = strlen(CTX.home_dir);
if(len && CTX.home_dir[len-1] != '/')
strcat(CTX.home_dir, "/");
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index eabea2713ffc605b1c1114456f7ef0282bc9ff7c..682897fc2a0b83eae82960416ba58ea4a20ae912 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.45 2007-11-04 21:03:17 remacle Exp $
+// $Id: MElement.cpp,v 1.46 2007-11-26 14:34:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -104,29 +104,14 @@ double MTriangle::gammaShapeMeasure()
double MTetrahedron::gammaShapeMeasure()
{
- double p0[3] = { _v[0]->x(), _v[0]->y(), _v[0]->z() };
- double p1[3] = { _v[1]->x(), _v[1]->y(), _v[1]->z() };
- double p2[3] = { _v[2]->x(), _v[2]->y(), _v[2]->z() };
- double p3[3] = { _v[3]->x(), _v[3]->y(), _v[3]->z() };
- double s1 = fabs(triangle_area(p0, p1, p2));
- double s2 = fabs(triangle_area(p0, p2, p3));
- double s3 = fabs(triangle_area(p0, p1, p3));
- double s4 = fabs(triangle_area(p1, p2, p3));
- double rhoin = 3. * fabs(getVolume()) / (s1 + s2 + s3 + s4);
- return 12. * rhoin / (sqrt(6.) * maxEdge());
+ double vol;
+ return qmTet(this,QMTET_2,&vol);
}
double MTetrahedron::etaShapeMeasure()
{
- double lij2 = 0.;
- for(int i = 0; i <= 3; i++) {
- for(int j = i + 1; j <= 3; j++) {
- double lij = _v[i]->distance(_v[j]);
- lij2 += lij * lij;
- }
- }
- double v = fabs(getVolume());
- return 12. * pow(0.9 * v * v, 1./3.) / lij2;
+ double vol;
+ return qmTet(this,QMTET_3,&vol);
}
SPoint3 MElement::barycenter()
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index e528a2d8820a98806103994a53a273d728256e92..77779abe43ff25424b5d256e37465268498092e3 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.26 2007-11-11 19:53:57 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.27 2007-11-26 14:34:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -193,10 +193,12 @@ int OCCEdge::minimumMeshSegments() const
int np;
if(geomType() == Line)
np= GEdge::minimumMeshSegments();
- else if(geomType() == Circle)
- np= CTX.mesh.min_circ_points - 1;
else
np=CTX.mesh.min_curv_points - 1;
+
+ // if the edge is closed, ensure that at least 3 points are generated
+ if (getBeginVertex() == getEndVertex()) np=std::max(2,np);
+
return np;
}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index df32a0832302f142684957e8c919839994bfed81..bb5ad7cf6c1f6862cecfaeeb2c94bd140e30b3c9 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.85 2007-11-11 19:53:57 remacle Exp $
+// $Id: BDS.cpp,v 1.86 2007-11-26 14:34:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -201,7 +201,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
Msg(DEBUG," edge %d %d has to be recovered",num1,num2);
int ix = 0;
- int ixMax = 100;
+ int ixMax = 300;
while(1)
{
std::vector<BDS_Edge *> intersected;
@@ -220,9 +220,9 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
{
std::set<EdgeToRecover>::iterator itr1 = e2r->find(EdgeToRecover(e->p1->iD,e->p2->iD,0));
std::set<EdgeToRecover>::iterator itr2 = e2r->find(EdgeToRecover(num1,num2,0));
- Msg(WARNING," edge %d %d on model edge %d cannot be recovered because it intersects %d %d on model edge %d",
- num1,num2,itr2->ge->tag(),
- e->p1->iD,e->p2->iD,itr1->ge->tag());
+ // Msg(WARNING," edge %d %d on model edge %d cannot be recovered because it intersects %d %d on model edge %d",
+ // num1,num2,itr2->ge->tag(),
+ // e->p1->iD,e->p2->iD,itr1->ge->tag());
// now throw a class that contains the diagnostic
not_recovered->insert (EdgeToRecover( num1 , num2, itr2->ge));
@@ -234,13 +234,21 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
++it;
}
- if (!intersected.size() || ix > ixMax)
+// if (ix > 300){
+// Msg(WARNING," edge %d %d cannot be recovered after %d iterations, trying again",
+// num1,num2,ix);
+// ix = 0;
+// }
+
+ if (!intersected.size() || ix > 1000)
{
BDS_Edge *eee = find_edge (num1, num2);
if (!eee)
{
outputScalarField(triangles, "debugp.pos",1);
outputScalarField(triangles, "debugr.pos",0);
+ Msg(GERROR," edge %d %d cannot be recovered at all, look at debugp.pos and debugr.pos",
+ num1,num2);
return 0;
}
return eee;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 20d9ca97a3baf36c44f4b87d6939a40099c03b81..08a2dcb922e227b1f37cd9f9e543dc1bfb55725e 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.29 2007-11-11 19:53:57 remacle Exp $
+// $Id: BackgroundMesh.cpp,v 1.30 2007-11-26 14:34:09 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -182,19 +182,22 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
double l3 = CTX.lc;
double l4 = lc_field.empty() ? MAX_LC : lc_field(X, Y, Z);
+
double lc;
if(l4 < MAX_LC && !CTX.mesh.constrained_bgmesh){
// use the fields unconstrained by other characteristic lengths
lc = l4 * CTX.mesh.lc_factor;
}
else{
- if(CTX.mesh.lc_from_curvature && ge->dim() < 3)
- l1 = LC_MVertex_CURV(ge, U, V);
if(ge->dim() < 2)
l2 = LC_MVertex_PNTS(ge, U, V);
lc = std::min(std::min(std::min(l1, l2), l3), l4) * CTX.mesh.lc_factor;
}
-
+
+
+ if(CTX.mesh.lc_from_curvature && ge->dim() < 3)
+ lc = std::min (lc,LC_MVertex_CURV(ge, U, V));
+
if(lc <= 0.){
Msg(GERROR, "Incorrect char. length lc = %g: using default instead", lc);
return l3 * CTX.mesh.lc_factor;
@@ -208,7 +211,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
// we do it also if CTX.mesh.constrained_bgmesh is true;
bool Extend1dMeshIn2dSurfaces()
{
- if(lc_field.empty()) return true;
+ if(lc_field.empty() && !CTX.mesh.lc_from_curvature) return true;
if(CTX.mesh.constrained_bgmesh) return true;
return false;
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 1ac4bb2b3b9d60e840d8c472bff278157fda4cbb..e16d722148b7a4d5f49ad30f07049fd6d2ce1f61 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.124 2007-11-04 21:03:17 remacle Exp $
+// $Id: Generator.cpp,v 1.125 2007-11-26 14:34:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -179,9 +179,16 @@ void Mesh1D(GModel *m)
void PrintMesh2dStatistics (GModel *m)
{
+ FILE *statreport = 0;
+ if (CTX.create_append_statreport == 1)
+ statreport = fopen (CTX.statreport,"w");
+ else if (CTX.create_append_statreport == 2)
+ statreport = fopen (CTX.statreport,"a");
+ else return;
+
double worst=1,best=0,avg=0;
double e_long=0,e_short=1.e22,e_avg=0;
- int nTotT=0,nTotE=0,nTotGoodLength=0,nTotGoodQuality=0;
+ int nTotT=0,nTotE=0,nTotGoodLength=0,nTotGoodQuality=0,nUnmeshed=0,numFaces=0;
Msg(INFO,"2D Mesh Statistics :");
for (GModel::fiter it = m->firstFace() ; it!=m->lastFace(); ++it)
{
@@ -194,22 +201,25 @@ void PrintMesh2dStatistics (GModel *m)
e_long = std::max((*it)->meshStatistics.longest_edge_length,e_long);
e_short = std::min((*it)->meshStatistics.smallest_edge_length,e_short);
- if ((*it)->meshStatistics.worst_element_shape < 0.05)
- Msg(INFO,"Badly Shaped Element (rho = %12.5E) on GFace %d (%d triangles)",(*it)->meshStatistics.worst_element_shape,(*it)->tag(),(*it)->meshStatistics.nbTriangle);
- if ((*it)->meshStatistics.status == GFace::FAILED)
- Msg(INFO,"GMSH was unable to mesh GFace %d ",(*it)->tag());
+ if ((*it)->meshStatistics.status == GFace::FAILED || (*it)->meshStatistics.status == GFace::PENDING)nUnmeshed++;
nTotT += (*it)->meshStatistics.nbTriangle ;
nTotE += (*it)->meshStatistics.nbEdge ;
nTotGoodLength += (*it)->meshStatistics.nbGoodLength ;
nTotGoodQuality+= (*it)->meshStatistics.nbGoodQuality ;
+ numFaces ++;
}
- Msg(INFO,"Element Quality (%d triangles) : avg %8.7f best %8.7f worst %8.7f greaterthan90 %d (%12.5E\\%)",
- nTotT,avg/(double)nTotT,best,worst,nTotGoodQuality,(double)nTotGoodQuality/nTotT);
- Msg(INFO,"Size Field Accuracy (%d edges) : %8.7f (should be as close as 1 as possible) in good brackets %d (%12.5E\%)",
- nTotE,exp(e_avg/(double)nTotE),nTotGoodLength,(double)nTotGoodLength/nTotE);
+ if (CTX.create_append_statreport == 1){
+ fprintf(statreport,"2D stats\tname\t\t#faces\t\t#fail\t\t#t\t\tQavg\t\tQbest\t\tQworst\t\t#Q>90\t\t#Q>90/#t\t#e\t\ttau\t\t#Egood\t\t#Egood/#e\tCPU\n");
+ }
+ fprintf(statreport,"\t%16s\t%d\t\t%d\t\t",CTX.base_filename,numFaces, nUnmeshed);
+ fprintf(statreport,"%d\t\t%8.7f\t%8.7f\t%8.7f\t%d\t\t%8.7f\t",
+ nTotT,avg/(double)nTotT,best,worst,nTotGoodQuality,(double)nTotGoodQuality/nTotT);
+ fprintf(statreport,"%d\t\t%8.7f\t%d\t\t%8.7f\t%8.1f\n",
+ nTotE,exp(e_avg/(double)nTotE),nTotGoodLength,(double)nTotGoodLength/nTotE,CTX.mesh_timer[1]);
+ fclose (statreport);
}
@@ -232,15 +242,29 @@ void Mesh2D(GModel *m)
// boundary layers are special: their generation (including vertices
// and curve meshes) is global as it depends on a smooth normal
// field generated from the surface mesh of the source surfaces
- if(!Mesh2DWithBoundaryLayers(m))
- std::for_each(m->firstFace(), m->lastFace(), meshGFace());
-
- PrintMesh2dStatistics(m);
+ if(!Mesh2DWithBoundaryLayers(m)){
+ std::for_each(m->firstFace(), m->lastFace(), meshGFace());
+ int nIter = 0;
+ while(1)
+ {
+ meshGFace mesher;
+ int nbPending = 0;
+ for (GModel::fiter it = m->firstFace() ; it!=m->lastFace(); ++it)
+ {
+ if ((*it)->meshStatistics.status == GFace::PENDING){mesher(*it);nbPending++;}
+ }
+ if (!nbPending)break;
+ if (nIter > 10)break;
+ }
+ }
+
double t2 = Cpu();
CTX.mesh_timer[1] = t2 - t1;
Msg(INFO, "Mesh 2D complete (%g s)", CTX.mesh_timer[1]);
Msg(STATUS1, "Mesh");
+
+ PrintMesh2dStatistics(m);
}
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 6ee88cba6dd77e4baf97b339cc2841c71b166bce..075f8a8a0a01a4e6e5396c705d2ddf5e624f2aac 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.187 2007-11-12 10:49:16 geuzaine Exp $
+# $Id: Makefile,v 1.188 2007-11-26 14:34:10 remacle Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -44,6 +44,7 @@ SRC = Generator.cpp \
meshGRegionTransfinite.cpp \
meshGRegionExtruded.cpp \
meshGRegionCarveHole.cpp \
+ meshGRegionLocalMeshMod.cpp\
DivideAndConquer.cpp \
BackgroundMesh.cpp \
qualityMeasures.cpp \
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e0589a529504013d42f434ba2a46476bb062d413..28a1cb1d01b3dcf7d0a2d865a4e4e0cb1be2edbd 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.101 2007-11-11 19:53:57 remacle Exp $
+// $Id: meshGFace.cpp,v 1.102 2007-11-26 14:34:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -337,7 +337,7 @@ bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
void OptimizeMesh(GFace *gf, BDS_Mesh &m, const int NIT)
{
// optimize
- // if (0)
+ if (0)
{
for(int i = 0 ; i < NIT ; i++){
{
@@ -602,30 +602,41 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
int MAXNP = m.MAXPOINTNUMBER;
- if (NIT > 0)
+ // if (NIT > 0)
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
{
- std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
- std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
- double L=0;
- int ne = 0;
- while(it!=ite){
- double l = (*it)->length();
- if ((*it)->g && (*it)->g->classif_degree == 1){
- L=ne?std::max(L,l):l;
- // L=ne?std::min(L,l):l;
- // L+=l;
- ne++;
+ if (NIT > 0)
+ {
+ std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
+ std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
+ double L=0;
+ int ne = 0;
+ while(it!=ite){
+ double l = (*it)->length();
+ if ((*it)->g && (*it)->g->classif_degree == 1){
+ L=ne?std::max(L,l):l;
+ // L=ne?std::min(L,l):l;
+ // L+=l;
+ ne++;
+ }
+ ++it;
+ }
+ if (!ne) L = 1.e22;
+ // else L/=ne;
+ if(!CTX.mesh.constrained_bgmesh)
+ (*itp)->lc() = L;
+ (*itp)->lcBGM() = L;
}
- ++it;
- }
- if (!ne) L = 1.e22;
- // else L/=ne;
- if(!CTX.mesh.constrained_bgmesh)
- (*itp)->lc() = L;
- (*itp)->lcBGM() = L;
+// else if (CTX.mesh.lc_from_curvature)
+// {
+// double Crv = gf->curvature(SPoint2((*itp)->u, (*itp)->v));
+// double lc = Crv > 0 ? 2*M_PI / Crv / CTX.mesh.min_circ_points : 1.e22;
+// lc *= CTX.mesh.lc_factor;
+// printf("%d = %d %12.5e %12.5E %12.5E\n",gf->tag(),CTX.mesh.min_circ_points,Crv,lc,(*itp)->lcBGM());
+// (*itp)->lc() = std::min(lc,(*itp)->lc());
+// }
++itp;
}
}
@@ -708,7 +719,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
Msg(DEBUG1," iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : %6d splits, %6d swaps, %6d collapses, %6d moves, accuracy %7.5f, shapes (worst %7.6f, avg %7.6f, best %7.6f) ",IT,minL,minE,maxL,maxE,nb_split,nb_swap,nb_collaps,nb_smooth, exp(mesh_quality/nE), worst,avg/nT,best);
if (nb_split==0 && nb_collaps == 0)break;
- if (mesh_quality < old_mesh_quality && worst > 0.1 && maxL < 1.45) break;
+ // if (mesh_quality < old_mesh_quality && worst > 0.1 && maxL < 1.45) break;
}
double t_total = t_spl + t_sw + t_col + t_sm;
@@ -1185,14 +1196,6 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
m->cleanup();
- if (debug)
- {
- char name[245];
- sprintf(name,"surface%d-recovered-real.pos",gf->tag());
- outputScalarField(m->triangles, name,0);
- sprintf(name,"surface%d-recovered-param.pos",gf->tag());
- outputScalarField(m->triangles, name,1);
- }
{
@@ -1213,12 +1216,24 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
}
}
+
+
m->cleanup();
m->del_point(m->find_point(-1));
m->del_point(m->find_point(-2));
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
+
+ if (debug)
+ {
+ char name[245];
+ sprintf(name,"surface%d-recovered-real.pos",gf->tag());
+ outputScalarField(m->triangles, name,0);
+ sprintf(name,"surface%d-recovered-param.pos",gf->tag());
+ outputScalarField(m->triangles, name,1);
+ }
+
// start mesh generation
if (!AlgoDelaunay2D ( gf ))
{
@@ -1597,8 +1612,8 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
for (std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin() ; it != gf->edgeLoops.end() ; it++)
{
std::vector<BDS_Point* > edgeLoop_BDS;
- // if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-7*LC2D))
- if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-5))
+ if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-7*LC2D))
+ if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-5*LC2D))
if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-3*LC2D))
{
gf->meshStatistics.status = GFace::FAILED;
@@ -1949,15 +1964,20 @@ void meshGFace::operator() (GFace *gf)
computeEdgeLoops(gf, points, indices);
// temp fix until we create MEdgeLoops in gmshFace
- Msg(DEBUG1, "Generating the mesh");
- if(noseam (gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
- gmsh2DMeshGenerator(gf,0, false);
- }
- else{
- if(!gmsh2DMeshGeneratorPeriodic(gf,false))
- Msg(GERROR, "Impossible to mesh face %d", gf->tag());
- }
-
+ // if (gf->tag() == 2)
+ {
+ Msg(DEBUG1, "Generating the mesh");
+ if(noseam (gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
+ gmsh2DMeshGenerator(gf,0, false);
+ }
+ else{
+ if(!gmsh2DMeshGeneratorPeriodic(gf,false))
+ Msg(GERROR, "Impossible to mesh face %d", gf->tag());
+ }
+ }
+// else
+// gf->meshStatistics.status = GFace::DONE;
+
Msg(DEBUG1, "type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index de2f45f678c06d9aa2100fa13375041ae71d20c3..83e3eab036246451b7c7911176464713e735c88e 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.20 2007-11-11 19:53:57 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.21 2007-11-26 14:34:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -19,7 +19,9 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
+#include "OS.h"
#include "BackgroundMesh.h"
+#include "edgeSwap.h"
#include "meshGRegionDelaunayInsertion.h"
#include "GModel.h"
#include "GRegion.h"
@@ -92,6 +94,8 @@ void connectTets ( ITER beg, ITER end)
}
}
}
+void connectTets ( std::list<MTet4*> & l) {connectTets(l.begin(),l.end());}
+void connectTets ( std::vector<MTet4*> &l ) {connectTets(l.begin(),l.end());}
void recurFindCavity ( std::list<faceXtet> & shell,
@@ -428,6 +432,98 @@ void recur_classify ( MTet4 *t ,
}
}
+template <class CONTAINER, class DATA>
+void gmshOptimizeMesh (CONTAINER &allTets, DATA &vSizes)
+{
+ double t1 = Cpu();
+ {
+ 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)
+ {
+ if (!(*it)->isDeleted()){
+ double vol;
+ double qual = qmTet((*it)->tet(),QMTET_2,&vol);
+ worst = std::min(qual,worst);
+ avg+=qual;
+ count++;
+ totalVolumeb+=vol;
+ }
+ }
+ Msg(INFO,"Optimization Vol = %12.5E QBAD %12.5E QAVG %12.5E",totalVolumeb,worst,avg/count);
+ }
+
+
+ while (1){
+ std::vector<MTet4*> newTets;
+
+ // get a bad tet and try to swap each of its edges
+
+ for (typename CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->isDeleted()){
+ double vol;
+ double qq = qmTet((*it)->tet(),QMTET_2,&vol);
+ // gmshSmoothVertex(*it,IED%4);
+ if (qq < .4)
+ for (int i=0;i<6;i++){
+ gmshEdgeSwap(newTets,*it,i,QMTET_2);
+ if ((*it)->isDeleted())i=10;
+ }
+ }
+ }
+
+ // relocate vertices
+// for (typename 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)
+// for (int i=0;i<4;i++){
+// gmshSmoothVertex(*it,i);
+// }
+// }
+// }
+
+
+ // 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()){
+ newTets[i]->setup(newTets[i]->tet(),vSizes);
+ allTets.insert(newTets[i]);
+ }
+ else{
+ delete newTets[i];
+ }
+ }
+
+ 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)
+ {
+ if (!(*it)->isDeleted()){
+ double vol;
+ double qual = qmTet((*it)->tet(),QMTET_2,&vol);
+ worst = std::min(qual,worst);
+ avg+=qual;
+ count++;
+ totalVolumeb+=vol;
+ }
+ }
+ double t2 = Cpu();
+ Msg(INFO,"Optimization Vol = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",totalVolumeb,worst,avg/count,t2-t1);
+ }
+}
+
+
void insertVerticesInRegion (GRegion *gr)
{
@@ -542,7 +638,7 @@ void insertVerticesInRegion (GRegion *gr)
else
myFactory.changeTetRadius(allTets.begin(),0.0);
}
- // Normally, a tet mesh contains about 5 to 6 times more tets than
+ // Normally, a tet mesh contains about 6 times more tets than
// vertices
// This allows to clean up the set of tets when lots of deleted ones
// are present in the mesh
@@ -562,6 +658,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 9309625ea2c4f683b443d0ed7d5899f4b21bf249..846385298af98f44a406f0e31396091de5df3899 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -54,11 +54,18 @@ class MTet4
MTet4 *neigh[4];
GRegion *gr;
+
+ public :
+
~MTet4 (){}
MTet4 () : deleted(false), base (0), gr(0),circum_radius(0.0)
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
}
+ MTet4 (MTetrahedron * t) : deleted(false), gr(0),circum_radius(0.0),base(t)
+ {
+ neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
+ }
void setup ( MTetrahedron * t, std::vector<double> & sizes)
{
@@ -79,8 +86,6 @@ class MTet4
deleted = false;
}
- public :
-
inline GRegion * onWhat () const {return gr;}
inline void setOnWhat (GRegion *g) {gr=g;}
@@ -129,6 +134,7 @@ class MTet4
};
void connectTets ( std::list<MTet4*> & );
+void connectTets ( std::vector<MTet4*> & );
void insertVerticesInRegion (GRegion *gr) ;
class compareTet4Ptr
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..87a4757c84d51eaf98547ec4723877e74bc3efc5
--- /dev/null
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -0,0 +1,418 @@
+#include "meshGRegionLocalMeshMod.h"
+
+static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
+static int efaces[6][2] = {{0,2},{0,1},{1,2},{0,3},{2,3},{1,3}};
+static int enofaces[6][2] = {{1,3},{2,3},{0,3},{1,2},{0,1},{0,2}};
+static int facesXedges[4][3] = {{0,1,3},{1,2,5},{0,2,4},{3,4,5}};
+static int faces[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
+
+// as input, we give a tet and an edge, as return, we get
+// all tets that share this edge and all vertices that are
+// forming the outer ring of the cavity
+// we return true if the cavity is closed and false if it is open
+bool gmshBuildEdgeCavity ( MTet4 *t,
+ int iLocalEdge,
+ MVertex **v1,MVertex **v2,
+ std::vector<MTet4*> &cavity,
+ std::vector<MTet4*> &outside,
+ std::vector<MVertex*> &ring ){
+ cavity.clear();
+ ring.clear();
+ outside.clear();
+
+ // FILE *f = fopen ("pts.pos","w");
+ // fprintf(f,"View \"\"{\n");
+
+ *v1 = t->tet()->getVertex(edges[iLocalEdge][0]);
+ *v2 = t->tet()->getVertex(edges[iLocalEdge][1]);
+
+ // printf("trying to swap %p %p (%p,%p,%p,%p)\n",(*v1),(*v2),t->tet()->getVertex(0),t->tet()->getVertex(1),t->tet()->getVertex(2),t->tet()->getVertex(3));
+
+ MVertex *lastinring = t->tet()->getVertex(edges[5-iLocalEdge][0]);
+
+ // fprintf(f,"SP(%g,%g,%g){%d};\n",(*v1)->x(),(*v1)->y(),(*v1)->z(),-1);
+ // fprintf(f,"SP(%g,%g,%g){%d};\n",(*v2)->x(),(*v2)->y(),(*v2)->z(),-2);
+
+ // fprintf(f,"SP(%g,%g,%g){%d};\n",lastinring->x(),lastinring->y(),lastinring->z(),ring.size());
+
+ // printf("the ring starts with %p \n",lastinring);
+
+ ring.push_back(lastinring);
+ cavity.push_back(t);
+
+ while (1){
+ MVertex *ov1 = t->tet()->getVertex(edges[5-iLocalEdge][0]);
+ MVertex *ov2 = t->tet()->getVertex(edges[5-iLocalEdge][1]);
+ int K = ov1 == lastinring ? 1 : 0;
+ lastinring = ov1 == lastinring ? ov2 : ov1;
+ // look in the 2 faces sharing this edge the one that has vertex
+ // ov2 i.e. edges[5-iLocalEdge][1]
+
+ int iFace;
+// int iFaceOK;
+// for (int i=0;i<4;i++)
+// {
+// MVertex *f1 = t->tet()->getVertex(faces[i][0]);
+// MVertex *f2 = t->tet()->getVertex(faces[i][1]);
+// MVertex *f3 = t->tet()->getVertex(faces[i][2]);
+// if ((f1 == *v1 && f2 == *v2 && f3 == lastinring) ||
+// (f1 == *v1 && f3 == *v2 && f2 == lastinring) ||
+// (f2 == *v1 && f1 == *v2 && f3 == lastinring) ||
+// (f2 == *v1 && f3 == *v2 && f1 == lastinring) ||
+// (f3 == *v1 && f2 == *v2 && f1 == lastinring) ||
+// (f3 == *v1 && f1 == *v2 && f2 == lastinring) )iFaceOK = i;
+// }
+
+ int iFace1 = efaces[iLocalEdge][0];
+ int iFace2 = efaces[iLocalEdge][1];
+ if (faces[iFace1][0] == edges[5-iLocalEdge][K] ||
+ faces[iFace1][1] == edges[5-iLocalEdge][K] ||
+ faces[iFace1][2] == edges[5-iLocalEdge][K] ) iFace = iFace1;
+ else if (faces[iFace2][0] == edges[5-iLocalEdge][K] ||
+ faces[iFace2][1] == edges[5-iLocalEdge][K] ||
+ faces[iFace2][2] == edges[5-iLocalEdge][K] ) iFace = iFace2;
+ else {printf("error of connexion\n");throw;}
+ // else iFace = iFace2;
+
+ // printf("iFaceOK %d iFace %d\n",iFaceOK,iFace);
+
+ if (!t->getNeigh(iFace))return false;
+ t=t->getNeigh(iFace);
+ if (t->isDeleted())throw;
+ // printf("next tet (%p,%p,%p,%p)\n",t->tet()->getVertex(0),t->tet()->getVertex(1),t->tet()->getVertex(2),t->tet()->getVertex(3));
+ // int iNoFace1 = enofaces[iLocalEdge][0];
+ // int iNoFace2 = enofaces[iLocalEdge][1];
+
+ if (t==*(cavity.begin())) break;
+ // fprintf(f,"SP(%g,%g,%g){%d};\n",lastinring->x(),lastinring->y(),lastinring->z(),ring.size());
+ ring.push_back(lastinring);
+ // printf("the ring continues with %p \n",lastinring);
+
+ cavity.push_back(t);
+ iLocalEdge = -1;
+ for (int i=0;i<6;i++)
+ {
+ MVertex *a = t->tet()->getVertex(edges[i][0]);
+ MVertex *b = t->tet()->getVertex(edges[i][1]);
+ if ((a == *v1 && b == *v2) || (a == *v2 && b == *v1)){
+ iLocalEdge = i;
+ break;
+ }
+ }
+ if (iLocalEdge == -1){
+ printf("ERROR : loc = %d\n",iLocalEdge);
+ throw;
+ }
+ }
+ // fprintf(f,"};\n");
+ // fclose(f);
+ // getchar();
+
+ for (int i=0;i<cavity.size();i++){
+ for (int j=0;j<4;j++){
+ MTet4 * neigh = cavity[i]->getNeigh(j);
+ if (neigh)
+ {
+ bool found = false;
+ for (int k=0;k<outside.size();k++)if(outside[k]==neigh)found=true;
+ if(!found)outside.push_back(neigh);
+ }
+ }
+ }
+
+ return true;
+}
+
+typedef struct {
+ int nbr_triangles ; /* number of different triangles */
+ int (*triangles)[3] ; /* triangles array */
+ int nbr_trianguls ; /* number of different triangulations */
+ int nbr_triangles_2 ; /* number of triangles / triangulation */
+ int (*trianguls)[5] ; /* retriangulations array */
+} SwapPattern ;
+
+
+void BuildSwapPattern3(SwapPattern *sc)
+{
+ static int trgl[][3] = { 0,1,2 } ;
+ static int trgul[][5] = { 0,-1,-1,-1,-1 } ;
+
+ sc->nbr_triangles = 1 ;
+ sc->nbr_triangles_2 = 1 ;
+ sc->nbr_trianguls = 1 ;
+ sc->triangles = trgl ;
+ sc->trianguls = trgul ;
+}
+
+void BuildSwapPattern4(SwapPattern *sc)
+{
+ static int trgl[][3] =
+ { 0,1,2, 0,2,3, 0,1,3, 1,2,3 } ;
+ static int trgul[][5] =
+ { 0,1,-1,-1,-1, 2,3,-1,-1,-1 } ;
+
+ sc->nbr_triangles = 4 ;
+ sc->nbr_triangles_2 = 2 ;
+ sc->nbr_trianguls = 2 ;
+ sc->triangles = trgl ;
+ sc->trianguls = trgul ;
+}
+
+
+void BuildSwapPattern5(SwapPattern *sc)
+{
+ static int trgl[][3] =
+ { 0,1,2, 0,2,3, 0,3,4, 0,1,4, 1,3,4, 1,2,3, 2,3,4, 0,2,4, 0,1,3, 1,2,4 } ;
+ static int trgul[][5] =
+ { 0,1,2,-1,-1, 3,4,5,-1,-1, 0,6,7,-1,-1, 2,5,8,-1,-1, 3,6,9,-1,-1 } ;
+
+ sc->nbr_triangles = 10 ;
+ sc->nbr_triangles_2 = 3 ;
+ sc->nbr_trianguls = 5 ;
+ sc->triangles = trgl ;
+ sc->trianguls = trgul ;
+}
+
+void BuildSwapPattern6(SwapPattern *sc)
+{
+ static int trgl[][3] =
+ { 0,1,2, 0,2,3, 0,3,4, 0,4,5, 0,2,5, 2,4,5, 2,3,4, 0,3,5,
+ 3,4,5, 0,2,4, 2,3,5, 1,2,3, 0,1,3, 0,1,5, 1,4,5, 1,3,4,
+ 0,1,4, 1,3,5, 1,2,4, 1,2,5 } ;
+ static int trgul[][5] =
+ { 0,1,2,3,-1, 0,4,5,6,-1, 0,1,7,8,-1, 0,3,6,9,-1, 0,4,8,10,-1,
+ 2,3,11,12,-1, 11,13,14,15,-1, 7,8,11,12,-1, 3,11,15,16,-1,
+ 8,11,13,17,-1, 6,13,14,18,-1, 3,6,16,18,-1, 5,6,13,19,-1,
+ 8,10,13,19,-1 } ;
+
+ sc->nbr_triangles = 20 ;
+ sc->nbr_triangles_2 = 4 ;
+ sc->nbr_trianguls = 14 ;
+ sc->triangles = trgl ;
+ sc->trianguls = trgul ;
+}
+
+void BuildSwapPattern7(SwapPattern *sc)
+{
+ static int trgl[][3] =
+ { 0,1,2, 0,2,3, 0,3,4, 0,4,5, 0,5,6, 0,3,6, 3,5,6, 3,4,5, 0,4,6,
+ 4,5,6, 0,3,5, 3,4,6, 0,2,4, 2,3,4, 0,2,6, 2,5,6, 2,4,5, 0,2,5,
+ 2,4,6, 2,3,5, 2,3,6, 0,1,3, 1,2,3, 0,1,4, 1,3,4, 0,1,6, 1,5,6,
+ 1,4,5, 0,1,5, 1,4,6, 1,3,5, 1,3,6, 1,2,4, 1,2,5, 1,2,6 };
+ static int trgul[][5] =
+ { 0,1,2,3,4, 0,1,5,6,7, 0,1,2,8,9, 0,1,4,7,10, 0,1,5,9,11, 0,3,4,12,13,
+ 0,13,14,15,16, 0,8,9,12,13, 0,4,13,16,17, 0,9,13,14,18, 0,7,14,15,19,
+ 0,4,7,17,19, 0,6,7,14,20, 0,9,11,14,20, 2,3,4,21,22, 5,6,7,21,22,
+ 2,8,9,21,22, 4,7,10,21,22, 5,9,11,21,22, 3,4,22,23,24, 22,24,25,26,27,
+ 8,9,22,23,24, 4,22,24,27,28, 9,22,24,25,29, 7,22,25,26,30, 4,7,22,28,30,
+ 6,7,22,25,31, 9,11,22,25,31, 3,4,13,23,32, 13,25,26,27,32, 8,9,13,23,32,
+ 4,13,27,28,32, 9,13,25,29,32, 13,16,25,26,33, 4,13,16,28,33,
+ 13,15,16,25,34, 9,13,18,25,34, 7,19,25,26,33, 4,7,19,28,33,
+ 7,15,19,25,34, 6,7,20,25,34, 9,11,20,25,34 } ;
+
+ sc->nbr_triangles = 35 ;
+ sc->nbr_triangles_2 = 5 ;
+ sc->nbr_trianguls = 42 ;
+ sc->triangles = trgl ;
+ sc->trianguls = trgul ;
+}
+
+
+void gmshEdgeSwap (std::vector<MTet4 *> &newTets,
+ MTet4 *tet,
+ int iLocalEdge,
+ const gmshQualityMeasure4Tet &cr){
+ std::vector<MTet4*> cavity;
+ std::vector<MTet4*> outside;
+ std::vector<MVertex*> ring;
+ MVertex *v1,*v2;
+ // printf("swapping an edge\n");
+
+ bool closed = gmshBuildEdgeCavity ( tet,iLocalEdge,&v1,&v2,cavity,outside,ring);
+
+
+ if (!closed)return;
+
+ 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);
+ // printf("%p %p -> %p %p %p %p\n",v1,v2,cavity[i]->tet()->getVertex(0),cavity[i]->tet()->getVertex(1),cavity[i]->tet()->getVertex(2),cavity[i]->tet()->getVertex(3));
+ volumeRef += fabs(vol);
+ }
+
+ // build swap patterns
+
+ SwapPattern sp;
+ switch (ring.size()){
+ case 3 : BuildSwapPattern3 (&sp); break;
+ case 4 : BuildSwapPattern4 (&sp); break;
+ case 5 : BuildSwapPattern5 (&sp); break;
+ case 6 : BuildSwapPattern6 (&sp); break;
+ case 7 : BuildSwapPattern7 (&sp); break;
+ default : return;
+ }
+
+ // printf("the cavity contains %d tets the ring is of size %d volume %12.5E qual %12.5E\n",cavity.size(),ring.size(),volumeRef,tetQualityRef);
+
+ // compute qualities of all tets that appear in the patterns
+ double tetQuality1[100],tetQuality2[100];
+ double volume1[100],volume2[100];
+ for (int i=0;i<sp.nbr_triangles;i++){
+ int p1 = sp.triangles[i][0];
+ int p2 = sp.triangles[i][1];
+ int p3 = sp.triangles[i][2];
+ tetQuality1[i] = qmTet ( ring[p1], ring[p2], ring[p3], v1 , cr , &(volume1[i]));
+ tetQuality2[i] = qmTet ( ring[p1], ring[p2], ring[p3], v2 , cr , &(volume2[i]));
+ // printf("points %d %d %d vol %g %g qual %g %g\n",p1,p2,p3,volume1[i],volume2[i],tetQuality1[i],tetQuality2[i]);
+ }
+
+
+
+ // look for the best triangulation, i.e. the one
+ // that maximize the minimum element quality
+ double minQuality[100];
+ // for all triangulations
+ for (int i=0;i<sp.nbr_trianguls;i++){
+ // for all triangles in a triangulation
+ minQuality[i] = 1;
+ double volume=0;
+ for (int j=0;j<sp.nbr_triangles_2;j++){
+ int iT = sp.trianguls[i][j];
+ minQuality[i] = std::min (minQuality[i],tetQuality1[iT]);
+ minQuality[i] = std::min (minQuality[i],tetQuality2[iT]);
+ volume += (volume1[iT] + volume2[iT] );
+ }
+ // printf("config %3d qual %12.5E volume %12.5E\n",i,minQuality[i],volume);
+ if (fabs(volume-volumeRef) > 1.e-12 * (volume+volumeRef))minQuality[i] = 0;
+ }
+
+ int iBest;
+ double best = 0.0;
+ for (int i=0;i<sp.nbr_trianguls;i++){
+ if(minQuality[i] > best)
+ {
+ best = minQuality[i];
+ iBest = i;
+ }
+ }
+
+ // if there exist no swap that enhance the quality
+ if (best <= tetQualityRef) return;
+
+// printf("swap is performed\n");
+
+// return;
+
+ // we have the best configuration, so we swap
+ // printf("outside size = %d\n",outside.size());
+
+ double volumeAssert = 0;
+ for (int j=0;j<sp.nbr_triangles_2;j++){
+ int iT = sp.trianguls[iBest][j];
+ int p1 = sp.triangles[iT][0];
+ int p2 = sp.triangles[iT][1];
+ int p3 = sp.triangles[iT][2];
+ MVertex *pv1 = ring[p1];
+ MVertex *pv2 = ring[p2];
+ MVertex *pv3 = ring[p3];
+ MTetrahedron *tr1 = new MTetrahedron ( pv1,
+ pv2,
+ pv3,
+ v1);
+ MTetrahedron *tr2 = new MTetrahedron ( pv3,
+ pv2,
+ pv1,
+ v2);
+ volumeAssert += (fabs(tr1->getVolume()) +fabs(tr2->getVolume()) );
+
+ MTet4 *t41 = new MTet4 ( tr1 );
+ MTet4 *t42 = new MTet4 ( tr2 );
+ t41->setOnWhat(cavity[0]->onWhat());
+ t42->setOnWhat(cavity[0]->onWhat());
+ outside.push_back(t41);
+ outside.push_back(t42);
+ newTets.push_back(t41);
+ newTets.push_back(t42);
+ }
+ if (fabs(volumeRef-volumeAssert) > 1.e-10 * volumeRef) printf("volumes = %12.5E %12.5E\n",volumeRef,volumeAssert);
+
+ for (int i=0;i<cavity.size();i++){
+ cavity[i]->setDeleted(true);
+ }
+
+ connectTets ( outside );
+
+}
+
+
+void gmshBuildVertexCavity_recur ( MTet4 *t,
+ MVertex *v,
+ std::vector<MTet4*> &cavity){
+ for (int i=0; i<4; i++){
+ MTet4 *neigh = t->getNeigh(i);
+ if (neigh){
+ bool foundv = false;
+ for (int j=0;j<4;j++){if (t->tet()->getVertex(j) == v){foundv = true; break;}}
+ if (foundv)
+ {
+ bool found = false;
+ for (int j=0;j<cavity.size();j++){if (cavity[j] == neigh){found = true; break;}}
+ if (!found){
+ cavity.push_back(neigh);
+ gmshBuildVertexCavity_recur ( neigh,v,cavity);
+ }
+ }
+ }
+ }
+}
+void gmshSmoothVertex ( MTet4 *t,
+ int iVertex){
+ std::vector<MTet4*> cavity;
+ 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);
+ worst = std::min(worst,q);
+ xcg += 0.25*(cavity[i]->tet()->getVertex(0)->x()+cavity[i]->tet()->getVertex(1)->x()+cavity[i]->tet()->getVertex(2)->x()+cavity[i]->tet()->getVertex(3)->x())*volume;
+ ycg += 0.25*(cavity[i]->tet()->getVertex(0)->y()+cavity[i]->tet()->getVertex(1)->y()+cavity[i]->tet()->getVertex(2)->y()+cavity[i]->tet()->getVertex(3)->y())*volume;
+ zcg += 0.25*(cavity[i]->tet()->getVertex(0)->z()+cavity[i]->tet()->getVertex(1)->z()+cavity[i]->tet()->getVertex(2)->z()+cavity[i]->tet()->getVertex(3)->z())*volume;
+ vTot += volume;
+
+ }
+ xcg /= (vTot);
+ ycg /= (vTot);
+ zcg /= (vTot);
+ double volumeAfter = 0.0;
+
+ double x = t->tet()->getVertex(iVertex)->x();
+ double y = t->tet()->getVertex(iVertex)->y();
+ double z = t->tet()->getVertex(iVertex)->z();
+
+ t->tet()->getVertex(iVertex)->x() = xcg;
+ t->tet()->getVertex(iVertex)->y() = ycg;
+ t->tet()->getVertex(iVertex)->z() = zcg;
+ double worstAfter = 1.0;
+ for (int i=0 ; i< cavity.size() ; i++){
+ double volume;
+ double q = qmTet(cavity[i]->tet(),QMTET_2,&volume);
+ volumeAfter += volume;
+ worstAfter =std::min(worstAfter,q);
+ }
+
+ if (fabs(volumeAfter-vTot) > 1.e-10 * vTot || worstAfter < worst){
+ t->tet()->getVertex(iVertex)->x() = x;
+ t->tet()->getVertex(iVertex)->y() = y;
+ t->tet()->getVertex(iVertex)->z() = z;
+ }
+}
+
+
diff --git a/Mesh/meshGRegionLocalMeshMod.h b/Mesh/meshGRegionLocalMeshMod.h
new file mode 100644
index 0000000000000000000000000000000000000000..e0dca76adacf917988801ed9160a4f677a657ee3
--- /dev/null
+++ b/Mesh/meshGRegionLocalMeshMod.h
@@ -0,0 +1,39 @@
+#ifndef _MESHGREGIONLOCALMESHMOD_H_
+#define _MESHGREGIONLOCALMESHMOD_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "meshGRegionDelaunayInsertion.h"
+#include "qualityMeasures.h"
+
+void gmshEdgeSwap (std::vector<MTet4 *> &newTets,
+ MTet4 *tet,
+ int iLocalEdge,
+ const gmshQualityMeasure4Tet &cr);
+void gmshFaceSwap (std::vector<MTet4 *> &newTets,
+ MTet4 *tet,
+ int iLocalFace,
+ const gmshQualityMeasure4Tet &cr);
+void gmshSmoothVertex ( MTet4 *t,
+ int iLocalVertex);
+
+#endif
+
+
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index ba640923dfd99f285f8ed166011dcce091c24407..10530237031ba4c6f540e7577c1f68fda0ae25ba 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -66,3 +66,77 @@ double qmTriangle ( const double &xa, const double &ya, const double &z
return quality;
}
+double qmTet ( MTetrahedron*t, const gmshQualityMeasure4Tet &cr, double *volume)
+{
+ return qmTet (t->getVertex(0),t->getVertex(1),t->getVertex(2),t->getVertex(3),cr,volume);
+}
+
+double qmTet ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const MVertex *v4, const gmshQualityMeasure4Tet &cr, double *volume)
+{
+ return qmTet (v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(),v3->x(),v3->y(),v3->z(),v4->x(),v4->y(),v4->z(),cr,volume);
+}
+
+double qmTet ( const double &x1, const double &y1, const double &z1,
+ const double &x2, const double &y2, const double &z2,
+ const double &x3, const double &y3, const double &z3,
+ const double &x4, const double &y4, const double &z4,
+ const gmshQualityMeasure4Tet &cr, double *volume ){
+ double quality;
+ switch(cr)
+ {
+ case QMTET_3:
+ {
+ double mat[3][3];
+ mat[0][0] = x2 - x1;
+ mat[0][1] = x3 - x1;
+ mat[0][2] = x4 - x1;
+ mat[1][0] = y2 - y1;
+ mat[1][1] = y3 - y1;
+ mat[1][2] = y4 - y1;
+ mat[2][0] = z2 - z1;
+ mat[2][1] = z3 - z1;
+ mat[2][2] = z4 - z1;
+ *volume = fabs(det3x3(mat)) / 6.;
+ double l = ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
+ l += ((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1));
+ l +=((x4-x1)*(x4-x1)+(y4-y1)*(y4-y1)+(z4-z1)*(z4-z1));
+ l +=((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2));
+ l +=((x4-x2)*(x4-x2)+(y4-y2)*(y4-y2)+(z4-z2)*(z4-z2));
+ l +=((x3-x4)*(x3-x4)+(y3-y4)*(y3-y4)+(z3-z4)*(z3-z4));
+ return 12. * pow(3*fabs(*volume),2./3.)/l;
+ }
+ case QMTET_2:
+ {
+ double mat[3][3];
+ mat[0][0] = x2 - x1;
+ mat[0][1] = x3 - x1;
+ mat[0][2] = x4 - x1;
+ mat[1][0] = y2 - y1;
+ mat[1][1] = y3 - y1;
+ mat[1][2] = y4 - y1;
+ mat[2][0] = z2 - z1;
+ mat[2][1] = z3 - z1;
+ mat[2][2] = z4 - z1;
+ *volume = fabs(det3x3(mat)) / 6.;
+ double p0[3] = { x1,y1,z1};
+ double p1[3] = { x2,y2,z2};
+ double p2[3] = { x3,y3,z3};
+ double p3[3] = { x4,y4,z4};
+ double s1 = fabs(triangle_area(p0, p1, p2));
+ double s2 = fabs(triangle_area(p0, p2, p3));
+ double s3 = fabs(triangle_area(p0, p1, p3));
+ double s4 = fabs(triangle_area(p1, p2, p3));
+ double rhoin = 3. * fabs(*volume) / (s1 + s2 + s3 + s4);
+ double l = sqrt ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
+ l = std::max(l,sqrt ((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1)));
+ l = std::max(l,sqrt ((x4-x1)*(x4-x1)+(y4-y1)*(y4-y1)+(z4-z1)*(z4-z1)));
+ l = std::max(l,sqrt ((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2)));
+ l = std::max(l,sqrt ((x4-x2)*(x4-x2)+(y4-y2)*(y4-y2)+(z4-z2)*(z4-z2)));
+ l = std::max(l,sqrt ((x3-x4)*(x3-x4)+(y3-y4)*(y3-y4)+(z3-z4)*(z3-z4)));
+ return 2. * sqrt(6.) * rhoin / l;
+ }
+ break;
+ default:
+ throw;
+ }
+}
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index 02429b4fb89da054386de2dde1ab5887fa3bd27e..6b4a974e03365aeccb159aecba3267bc5785bdf4 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -5,7 +5,9 @@ class BDS_Point;
class BDS_Face;
class MVertex;
class MTriangle;
+class MTetrahedron;
enum gmshQualityMeasure4Triangle {QMTRI_RHO};
+enum gmshQualityMeasure4Tet {QMTET_1,QMTET_2,QMTET_3};
double qmTriangle ( MTriangle *f, const gmshQualityMeasure4Triangle &cr);
double qmTriangle ( BDS_Face *f, const gmshQualityMeasure4Triangle &cr);
@@ -16,5 +18,12 @@ double qmTriangle ( const double &x1, const double &y1, const double &z
const double &x2, const double &y2, const double &z2,
const double &x3, const double &y3, const double &z3,
const gmshQualityMeasure4Triangle &cr);
+double qmTet ( MTetrahedron *t, const gmshQualityMeasure4Tet &cr, double *volume = 0);
+double qmTet ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const MVertex *v4, const gmshQualityMeasure4Tet &cr, double *volume = 0);
+double qmTet ( const double &x1, const double &y1, const double &z1,
+ const double &x2, const double &y2, const double &z2,
+ const double &x3, const double &y3, const double &z3,
+ const double &x4, const double &y4, const double &z4,
+ const gmshQualityMeasure4Tet &cr, double *volume = 0);
#endif
diff --git a/Plugin/CutMap.cpp b/Plugin/CutMap.cpp
index 46a9338cad29c439d437264bd8b9ffc802ed4cc3..ba8e8cffa02f3508ff9d4f87cb90db7abb9c0d73 100644
--- a/Plugin/CutMap.cpp
+++ b/Plugin/CutMap.cpp
@@ -1,4 +1,4 @@
-// $Id: CutMap.cpp,v 1.55 2007-09-24 08:14:29 geuzaine Exp $
+// $Id: CutMap.cpp,v 1.56 2007-11-26 14:34:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -25,7 +25,7 @@
extern Context_T CTX;
StringXNumber CutMapOptions_Number[] = {
- {GMSH_FULLRC, "A", GMSH_CutMapPlugin::callbackA, 1.},
+ {GMSH_FULLRC, "A", GMSH_CutMapPlugin::callbackA, 0.},
{GMSH_FULLRC, "dTimeStep", NULL, -1.},
{GMSH_FULLRC, "dView", NULL, -1.},
{GMSH_FULLRC, "ExtractVolume", GMSH_CutMapPlugin::callbackVol, 0.},
diff --git a/Plugin/Levelset.cpp b/Plugin/Levelset.cpp
index c61617fa443f886690fc2876a5a700f9cc384765..76ef9d10acae547fb3548a495b4fa6cb4743bb30 100644
--- a/Plugin/Levelset.cpp
+++ b/Plugin/Levelset.cpp
@@ -1,4 +1,4 @@
-// $Id: Levelset.cpp,v 1.37 2007-09-24 08:14:30 geuzaine Exp $
+// $Id: Levelset.cpp,v 1.38 2007-11-26 14:34:10 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -781,7 +781,7 @@ static bool recur_sign_change(adapt_triangle *t, double val,
double v1 = plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
double v2 = plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
double v3 = plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
- if(v1 * v2 > 0 && v1 * v3 > 0)
+ if(v1 * v2 > 0 && v1 * v3 > 0 && v2 * v3 > 0)
t->visible = false;
else
t->visible = true;
diff --git a/Post/AdaptiveViews.cpp b/Post/AdaptiveViews.cpp
index 1ecc325f6f8480ae0331cdd199efc0e61148986c..c78bd253ed1fe2d469f8b33e6f6cfcc03cbdb561 100644
--- a/Post/AdaptiveViews.cpp
+++ b/Post/AdaptiveViews.cpp
@@ -756,10 +756,10 @@ int Adaptive_Post_View::zoomElement(int ielem,
itt = ELEM::all_elems.begin();
- for(; itt != itte; itt++) {
- if((*itt)->visible && !(*itt)->e[0] && level != levelmax)
- return 0;
- }
+ for(; itt != itte; itt++) {
+ if((*itt)->visible && !(*itt)->e[0] && level != levelmax)
+ return 0;
+ }
itt = ELEM::all_elems.begin();
adapt_point **p;
@@ -930,8 +930,7 @@ void Adaptive_Post_View::setAdaptiveResolutionLevel_TEMPL(int level, int levelma
kk++;
}
- for(int i = 0; i < nbelm; ++i) {
- if(!done[i])
+ for(int i = 0; i < nbelm; ++i) {// if(!done[i])
done[i] = zoomElement<ELEM>(i, level, levelmax, plug, *myList, counter);
}
}