diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index a2eba0339393a341034cd50a7f4d553d02fd80aa..defc56a3edacc85b28fc298bb74d29b2007e240f 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.106 2007-09-26 20:51:57 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.107 2007-11-04 21:03:16 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -90,6 +90,7 @@ void Print_Usage(char *name){
Msg(DIRECT, " -clscale float Set characteristic length scaling factor");
Msg(DIRECT, " -clcurv Compute characteristic lengths from curvatures");
Msg(DIRECT, " -epslc1d Set the accuracy of the evaluation of the LCFIELD for 1D mesh");
+ Msg(DIRECT, " -swapangle Set the treshold angle (in degree) between two adjacent faces below which a swap is allowed");
Msg(DIRECT, " -rand float Set random perturbation factor");
Msg(DIRECT, " -bgm file Load background mesh from file");
Msg(DIRECT, " -constrain Constrain background mesh with characteristic lengths");
@@ -368,6 +369,21 @@ void Get_Options(int argc, char *argv[])
exit(1);
}
}
+ else if(!strcmp(argv[i] + 1, "swapangle")) {
+ i++;
+ if(argv[i] != NULL) {
+ CTX.mesh.allow_swap_edge_angle = atof(argv[i++]);
+ if(CTX.mesh.allow_swap_edge_angle <= 0.0) {
+ fprintf(stderr, ERROR_STR
+ "Treshold angle for edge swap must be > 0\n");
+ exit(1);
+ }
+ }
+ else {
+ fprintf(stderr, ERROR_STR "Missing number\n");
+ exit(1);
+ }
+ }
else if(!strcmp(argv[i] + 1, "clcurv")) {
CTX.mesh.lc_from_curvature = 1;
i++;
diff --git a/Common/Context.h b/Common/Context.h
index 858fc3e8e9ed06ea62f754ac3052e7cbbddd27a6..ca3d3d55cc83a277adce3ffb400382f8f3329296 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -195,6 +195,7 @@ public :
char *triangle_options;
int smooth_normals, reverse_all_normals;
double angle_smooth_normals;
+ double allow_swap_edge_angle;
} mesh;
// post processing options
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 0a9c4d0476b2aff17bc67823dc29bb288ba19561..53253eb6a25f4fc404be8d0c8523627fd1af954e 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -888,10 +888,10 @@ StringXNumber MeshOptions_Number[] = {
"3D mesh algorithm (1=Tetgen+Delaunay, 4=Netgen)" },
{ F|O, "AngleSmoothNormals" , opt_mesh_angle_smooth_normals , 30.0 ,
"Threshold angle below which normals are not smoothed" },
-
+ { F|O, "AllowSwapAngle" , opt_mesh_allow_swap_edge_angle , 10.0 ,
+ "Treshold angle (in degrees) between faces normals under which we allow an edge swap" },
{ F|O, "BdfFieldFormat" , opt_mesh_bdf_field_format , 1. ,
"Field format for Nastran BDF files (0=free, 1=small, 2=large)" },
-
{ F|O, "CharacteristicLengthFactor" , opt_mesh_lc_factor , 1.0 ,
"Factor applied to all characteristic lengths" },
{ F|O, "CharacteristicLengthFromCurvature" , opt_mesh_lc_from_curvature , 0. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 7a53f9232e51574f78a410e54e7df3f0836481dc..90742cbcfa208192403187cefdc68d41fe3d9f6d 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.365 2007-10-14 09:51:17 geuzaine Exp $
+// $Id: Options.cpp,v 1.366 2007-11-04 21:03:16 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -4865,6 +4865,13 @@ double opt_mesh_min_circ_points(OPT_ARGS_NUM)
return CTX.mesh.min_circ_points;
}
+double opt_mesh_allow_swap_edge_angle(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.mesh.allow_swap_edge_angle = val;
+ return CTX.mesh.allow_swap_edge_angle;
+}
+
double opt_mesh_min_curv_points(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 2449d31589e8292e8114541ef6897c64aaedf995..864b3aa168d074db3fb91d76678a2619b59c52ba 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -476,6 +476,7 @@ double opt_mesh_algo3d(OPT_ARGS_NUM);
double opt_mesh_lc_integration_precision(OPT_ARGS_NUM);
double opt_mesh_recombine_algo(OPT_ARGS_NUM);
double opt_mesh_min_circ_points(OPT_ARGS_NUM);
+double opt_mesh_allow_swap_edge_angle(OPT_ARGS_NUM);
double opt_mesh_min_curv_points(OPT_ARGS_NUM);
double opt_mesh_constrained_bgmesh(OPT_ARGS_NUM);
double opt_mesh_order(OPT_ARGS_NUM);
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 635b2e8abf4588d42e5303ca0fbdbfa60ea03494..7adce0a57ff53f95b4f2468f30465d7a5f561cf4 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.38 2007-10-16 20:00:06 geuzaine Exp $
+// $Id: GFace.cpp,v 1.39 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -39,6 +39,7 @@ extern Context_T CTX;
GFace::GFace(GModel *model, int tag) : GEntity(model, tag), r1(0), r2(0)
{
+ meshStatistics.status = GFace::PENDING;
resetMeshAttributes();
}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 6305156aebf818388241f033099aff2a9175afb7..9fbf063fe41f2f960f3eb208615c3f0ba21674b3 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -179,6 +179,15 @@ class GFace : public GEntity
// edge loops
} meshAttributes ;
+ typedef enum {PENDING,DONE,FAILED} meshGenerationStatus ;
+ struct {
+ meshGenerationStatus status;
+ double worst_element_shape, best_element_shape, average_element_shape;
+ double smallest_edge_length, longest_edge_length, efficiency_index;
+ int nbEdge, nbTriangle;
+ int nbGoodQuality,nbGoodLength;
+ } meshStatistics;
+
// a crude graphical representation using a "cross" defined by pairs
// of start/end points
std::vector<SPoint3> cross;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 442984294ded35f32e96313dd54f361e8607844c..eabea2713ffc605b1c1114456f7ef0282bc9ff7c 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.44 2007-10-03 19:40:41 geuzaine Exp $
+// $Id: MElement.cpp,v 1.45 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -26,6 +26,7 @@
#include "Numeric.h"
#include "Message.h"
#include "Context.h"
+#include "qualityMeasures.h"
extern Context_T CTX;
@@ -96,6 +97,11 @@ double MElement::rhoShapeMeasure()
return 0.;
}
+double MTriangle::gammaShapeMeasure()
+{
+ return qmTriangle(this,QMTRI_RHO);
+}
+
double MTetrahedron::gammaShapeMeasure()
{
double p0[3] = { _v[0]->x(), _v[0]->y(), _v[0]->z() };
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 79c9e77e0cfb1e32075a89bf098b4501abaf4a23..1b4630561e5475162585cfc193f9544f488b137f 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -326,6 +326,7 @@ class MTriangle : public MElement {
mat[1][1] = _v[2]->y() - _v[0]->y();
}
void circumcenterXY(double *res) const;
+ virtual double gammaShapeMeasure();
void circumcenterUV(GFace*,double *res);
static void circumcenterXYZ(double *p1, double *p2, double *p3, double *res,
double *uv = 0);
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 67e1e0072854befcb0c393886d51cd15a56ba52f..cb5c930fbd3983f51563fed7a91b83a27874e61c 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.24 2007-10-14 09:51:17 geuzaine Exp $
+// $Id: OCCEdge.cpp,v 1.25 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -58,6 +58,7 @@ void OCCEdge::setTrimmed (OCCFace *f)
}
}
+// Plugin Netgen Salome !
SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
{
const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
@@ -106,9 +107,11 @@ int OCCEdge::isSeam(GFace *face) const
{
const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
BRepAdaptor_Surface surface(*s);
- if(surface.IsUPeriodic() || surface.IsVPeriodic()){
+ // printf("asking if edge %d is a seam of face %d\n",tag(),face->tag());
+ // printf("periodic %d %d\n",surface.IsUPeriodic(),surface.IsVPeriodic());
+ // if(surface.IsUPeriodic() || surface.IsVPeriodic()){
return BRep_Tool::IsClosed(c, *s);
- }
+ // }
return 0;
}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index f41a834dd2a1abc93934b3ca1646d6763af441f4..ab92212aa10d4b8f36b28ccfa3d685c00fa955eb 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.83 2007-10-14 19:54:16 remacle Exp $
+// $Id: BDS.cpp,v 1.84 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -26,6 +26,7 @@
#include "BDS.h"
#include "Message.h"
#include "GFace.h"
+#include "qualityMeasures.h"
bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Face *t);
@@ -82,12 +83,7 @@ void normal_triangle(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3,
double c[3])
{
vector_triangle(p1, p2, p3, c);
- double l = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
- if(l == 0)
- return;
- c[0] /= l;
- c[1] /= l;
- c[2] /= l;
+ norme(c);
}
double surface_triangle(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3)
@@ -104,21 +100,6 @@ double surface_triangle_param(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3)
return (0.5 * c);
}
-double quality_triangle(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3)
-{
- double e12 = ((p1->X - p2->X) * (p1->X - p2->X) +
- (p1->Y - p2->Y) * (p1->Y - p2->Y) +
- (p1->Z - p2->Z) * (p1->Z - p2->Z));
- double e22 = ((p1->X - p3->X) * (p1->X - p3->X) +
- (p1->Y - p3->Y) * (p1->Y - p3->Y) +
- (p1->Z - p3->Z) * (p1->Z - p3->Z));
- double e32 = ((p3->X - p2->X) * (p3->X - p2->X) +
- (p3->Y - p2->Y) * (p3->Y - p2->Y) +
- (p3->Z - p2->Z) * (p3->Z - p2->Z));
- double a = surface_triangle(p1, p2, p3);
- return a / (e12 + e22 + e32);
-}
-
void BDS_Point::getTriangles(std::list < BDS_Face * >&t) const
{
t.clear();
@@ -530,7 +511,7 @@ BDS_Mesh ::~ BDS_Mesh ()
}
-void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
+bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
{
/*
p1
@@ -552,6 +533,13 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
+ e->oppositeof(op);
+
+// double l1 = sqrt((op[0]->X-op[1]->X) *(op[0]->X-op[1]->X) +
+// (op[0]->Y-op[1]->Y) *(op[0]->Y-op[1]->Y) +
+// (op[0]->Z-op[1]->Z) *(op[0]->Z-op[1]->Z) );
+// if (l1 < 0.5* mid->lc()) return false;
+
BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
@@ -568,7 +556,6 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
}
// we should project
- e->oppositeof(op);
BDS_GeomEntity *g1 = 0, *g2 = 0, *ge = e->g;
BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
@@ -642,6 +629,7 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
p2->config_modified = true;
op[0]->config_modified = true;
op[1]->config_modified = true;
+ return true;
}
@@ -1160,16 +1148,32 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf)
V /= (3.*sTot);
LC /= (3.*sTot);
+ GPoint gp = gf->point(U*scalingU,V*scalingV);
+
+ const double oldX = p->X;
+ const double oldY = p->Y;
+ const double oldZ = p->Z;
+
+ double oldWorst=1.;
+ double newWorst=1.;
+
it = ts.begin();
while(it != ite) {
BDS_Face *t = *it;
if (!test_move_point_parametric_triangle ( p, U, V, t))
- return false;
+ return false;
+// // p->X = gp.x();
+// // p->Y = gp.y();
+// // p->Z = gp.z();
+// // newWorst = std::min(newWorst,qmTriangle(*it,QMTRI_RHO));
+// // p->X = oldX;
+// // p->Y = oldY;
+// // p->Z = oldZ;
+// // oldWorst = std::min(oldWorst,qmTriangle(*it,QMTRI_RHO));
++it;
}
+ // if (oldWorst > newWorst)return false;
-
- GPoint gp = gf->point(U*scalingU,V*scalingV);
p->u = U;
p->v = V;
p->lc() = LC;
@@ -1190,10 +1194,10 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
if (!p->config_modified)return false;
if(p->g && p->g->classif_degree <= 1)
return false;
-
+
double U = 0;
double V = 0;
- double tot_length = 0;
+ double tot_length = 0;
double LC = 0;
std::list < BDS_Edge * >::iterator eit = p->edges.begin();
@@ -1208,8 +1212,6 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
++eit;
}
- //U /= (p->edges.size());
- //V /= (p->edges.size());
U /= tot_length;
V /= tot_length;
LC /= p->edges.size();
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index d3b7c743103dc5e2595acf797392b84fd12df4ce..6043487b1888c5eaa5859a6e947bbabd33a4aacd 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -45,7 +45,6 @@ void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
-double quality_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
class BDS_GeomEntity
@@ -456,7 +455,7 @@ public:
bool smooth_point_parametric(BDS_Point * p, GFace *gf);
bool smooth_point_centroid(BDS_Point * p, GFace *gf);
bool move_point(BDS_Point *p , double X, double Y, double Z);
- void split_edge(BDS_Edge *, BDS_Point *);
+ bool split_edge(BDS_Edge *, BDS_Point *);
void saturate_edge(BDS_Edge *, std::vector<BDS_Point *>&);
bool edge_constraint ( BDS_Point *p1, BDS_Point *p2 );
bool recombine_edge ( BDS_Edge *e );
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 85604ef54ba52c025984d4821f1797a899e5ac45..1ac4bb2b3b9d60e840d8c472bff278157fda4cbb 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.123 2007-09-10 04:47:04 geuzaine Exp $
+// $Id: Generator.cpp,v 1.124 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -177,6 +177,42 @@ void Mesh1D(GModel *m)
Msg(STATUS1, "Mesh");
}
+void PrintMesh2dStatistics (GModel *m)
+{
+ 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;
+ Msg(INFO,"2D Mesh Statistics :");
+ for (GModel::fiter it = m->firstFace() ; it!=m->lastFace(); ++it)
+ {
+ worst = std::min((*it)->meshStatistics.worst_element_shape,worst);
+ best = std::max((*it)->meshStatistics.best_element_shape,best);
+ avg += (*it)->meshStatistics.average_element_shape * (*it)->meshStatistics.nbTriangle ;
+
+ e_avg += (*it)->meshStatistics.efficiency_index ;//* (*it)->meshStatistics.nbEdge;
+
+ 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());
+
+ nTotT += (*it)->meshStatistics.nbTriangle ;
+ nTotE += (*it)->meshStatistics.nbEdge ;
+ nTotGoodLength += (*it)->meshStatistics.nbGoodLength ;
+ nTotGoodQuality+= (*it)->meshStatistics.nbGoodQuality ;
+ }
+
+ 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);
+
+
+}
+
void Mesh2D(GModel *m)
{
if(TooManyElements(m, 2)) return;
@@ -198,6 +234,8 @@ void Mesh2D(GModel *m)
// field generated from the surface mesh of the source surfaces
if(!Mesh2DWithBoundaryLayers(m))
std::for_each(m->firstFace(), m->lastFace(), meshGFace());
+
+ PrintMesh2dStatistics(m);
double t2 = Cpu();
CTX.mesh_timer[1] = t2 - t1;
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 7131843c4332c1cc2e879825bd1d5ce042c6963e..f802a394610c51e312fbffe38363d2e7dfe9f94d 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.184 2007-09-22 20:35:18 geuzaine Exp $
+# $Id: Makefile,v 1.185 2007-11-04 21:03:17 remacle Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -46,6 +46,7 @@ SRC = Generator.cpp \
meshGRegionCarveHole.cpp \
DivideAndConquer.cpp \
BackgroundMesh.cpp \
+ qualityMeasures.cpp \
BoundaryLayer.cpp \
BDS.cpp \
HighOrder.cpp
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 60cd78c686aea8c13b77ace5d0778964ef0e9188..500f16f1b384398e4d1313d3e08f9cc415a9fe7e 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.99 2007-10-16 20:00:07 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.100 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -34,6 +34,7 @@
#include "Message.h"
#include "Numeric.h"
#include "BDS.h"
+#include "qualityMeasures.h"
extern Context_T CTX;
@@ -263,35 +264,6 @@ double NewGetLc(BDS_Point *p1,BDS_Point *p2, GFace *f)
return 2*linearLength / (l1 + l2);
}
-bool edgeSwapTest(BDS_Edge *e,GFace *gf)
-{
- BDS_Point *op[2];
-
- if(!e->p1->config_modified && ! e->p2->config_modified) return false;
-
- if(e->numfaces() != 2) return false;
-
- e->oppositeof (op);
-
- double edgeLength1, edgeLength2;
-
- edgeLength1 = computeEdgeLinearLength(e,gf);
- edgeLength2 = computeEdgeLinearLength(op[0], op[1],gf);
- double lp1 = NewGetLc(e->p1);
- double lp2 = NewGetLc(e->p1);
- double lo1 = NewGetLc(op[0]);
- double lo2 = NewGetLc(op[1]);
-
-
- double el1 = 2*edgeLength1 / (lp1 + lp2);
- double el2 = 2*edgeLength2 / (lo1 + lo2);
-
- double q1 = fabs(1-el1);
- double q2 = fabs(1-el2);
-
- return q2 < 0.5*q1;
-}
-
void fourthPoint (double *p1, double *p2, double *p3, double *p4)
{
double c[3];
@@ -308,11 +280,49 @@ void fourthPoint (double *p1, double *p2, double *p3, double *p4)
p4[2] = c[2] + R * vz[2];
}
-bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
+
+bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
{
+ BDS_Face *f1 = e->faces(0);
+ BDS_Face *f2 = e->faces(1);
+ BDS_Point *n1[4];
+ BDS_Point *n2[4];
+ f1->getNodes(n1);
+ f2->getNodes(n2);
+ double norm1[3];
+ double norm2[3];
+ normal_triangle(n1[0],n1[1],n1[2],norm1);
+ normal_triangle(n2[0],n2[1],n2[2],norm2);
+ double cosa;prosca (norm1,norm2,&cosa);
+ return cosa > min_cos;
+}
+
+int edgeSwapTestQuality(BDS_Edge *e, double fact = 1.1)
+{
+ BDS_Point *op[2];
+
+ if(!e->p1->config_modified && ! e->p2->config_modified) return 0;
+
+ if(e->numfaces() != 2) return 0;
+
+ e->oppositeof (op);
- if (CTX.mesh.algo2d == ALGO_2D_MESHADAPT)
- return edgeSwapTest (e,gf);
+ if (! edgeSwapTestAngle(e, cos(CTX.mesh.allow_swap_edge_angle*M_PI/180.)) ) return -1;
+
+
+ double qa1 = qmTriangle(e->p1, e->p2, op[0],QMTRI_RHO);
+ double qa2 = qmTriangle(e->p1, e->p2, op[1],QMTRI_RHO);
+ double qb1 = qmTriangle(e->p1, op[0], op[1],QMTRI_RHO);
+ double qb2 = qmTriangle(e->p2, op[0], op[1],QMTRI_RHO);
+ double qa = std::min(qa1,qa2);
+ double qb = std::min(qb1,qb2);
+ if(qb > fact*qa) return 1;
+ if(qb < qa/fact) return -1;
+ return 0;
+}
+
+bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
+{
BDS_Point *op[2];
@@ -334,28 +344,6 @@ bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
return result > 0.;
}
-
-int edgeSwapTestQuality(BDS_Edge *e, double fact = 1.1)
-{
- BDS_Point *op[2];
-
- if(!e->p1->config_modified && ! e->p2->config_modified) return false;
-
- if(e->numfaces() != 2) return false;
-
- e->oppositeof (op);
-
- double qa1 = quality_triangle(e->p1, e->p2, op[0]);
- double qa2 = quality_triangle(e->p1, e->p2, op[1]);
- double qb1 = quality_triangle(e->p1, op[0], op[1]);
- double qb2 = quality_triangle(e->p2, op[0], op[1]);
- double qa = (qa1<qa2) ? qa1 : qa2;
- double qb = (qb1<qb2) ? qb1 : qb2;
- if(qb > fact*qa) return 1;
- if(qb < qa/fact) return -1;
- return 0;
-}
-
void OptimizeMesh(GFace *gf, BDS_Mesh &m, const int NIT)
{
// optimize
@@ -387,10 +375,12 @@ void OptimizeMesh(GFace *gf, BDS_Mesh &m, const int NIT)
if (NN2++ >= NN1)break;
if (!(*it)->deleted)
{
- int result = edgeSwapTestQuality(*it,5);
- if (result >= 0)
- if(edgeSwapTestDelaunay(*it,gf))
- m.swap_edge ( *it , BDS_SwapEdgeTestParametric());
+ const double qual = CTX.mesh.algo2d == ALGO_2D_MESHADAPT ? 1 : 5;
+ int result = edgeSwapTestQuality(*it,qual);
+ if (CTX.mesh.algo2d == ALGO_2D_MESHADAPT && result == 1)
+ m.swap_edge ( *it , BDS_SwapEdgeTestParametric());
+ else if ( result >= 0 && edgeSwapTestDelaunay(*it,gf))
+ m.swap_edge ( *it , BDS_SwapEdgeTestParametric());
}
++it;
}
@@ -412,13 +402,12 @@ void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
// result = 1 => oblige to swap because the quality is greatly improved
if (!(*it)->deleted)
{
- int result = edgeSwapTestQuality(*it,5);
- if (result >= 0)
- {
- if(edgeSwapTestDelaunay(*it,gf) || result > 0)
- if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
- nb_swap++;
- }
+ const double qual = CTX.mesh.algo2d == ALGO_2D_MESHADAPT ? 1 : 5;
+ int result = edgeSwapTestQuality(*it,qual);
+ if (CTX.mesh.algo2d == ALGO_2D_MESHADAPT && result == 1)
+ { if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))nb_swap++; }
+ else if ( result >= 0 && edgeSwapTestDelaunay(*it,gf))
+ { if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric())) nb_swap++; }
}
++it;
}
@@ -455,8 +444,8 @@ void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
mid->X,mid->Y,mid->Z);
//mid->lc() = 2./ ( 1./(*it)->p1->lc() + 1./(*it)->p2->lc() );
mid->lc() = 0.5 * ( (*it)->p1->lc() + (*it)->p2->lc() );
- m.split_edge ( *it, mid );
- nb_split++;
+ if(!m.split_edge ( *it, mid )) m.del_point(mid);
+ else nb_split++;
}
}
++it;
@@ -536,6 +525,77 @@ void collapseEdgePass ( GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb
}
+void computeMeshSizeFieldAccuracy (GFace *gf, BDS_Mesh &m, double &avg, double &max_e, double &min_e, int &nE, int &GS)
+{
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ avg=0;
+ min_e = 1.e22;
+ max_e = 0;
+ nE = 0;
+ GS = 0;
+ double oneoversqr2 = 1./sqrt(2);
+ double sqr2 = sqrt(2);
+ while (it!= m.edges.end())
+ {
+ if (!(*it)->deleted)
+ {
+ double lone = NewGetLc ( *it,gf);
+ if (lone > oneoversqr2 && lone < sqr2 ) GS++;
+
+ avg+=lone>1 ? (1./lone) - 1. : lone - 1.;
+ max_e = std::max(max_e,lone);
+ min_e = std::min(min_e,lone);
+ nE++;
+ }
+ ++it;
+ }
+}
+
+bool computeElementShapes (GFace *gf, BDS_Mesh &m, double &worst, double &avg, double &best, int &nT, int &nbGQ)
+{
+ std::list<BDS_Face*>::iterator it = m.triangles.begin();
+ worst = 1.e22;
+ avg = 0.0;
+ best = 0.0;
+ nT = 0;
+ nbGQ = 0;
+ while (it!= m.triangles.end())
+ {
+ if (!(*it)->deleted)
+ {
+ double q = qmTriangle(*it,QMTRI_RHO);
+ if (q>.9) nbGQ++;
+ avg+=q;
+ worst = std::min(worst,q);
+ best = std::max(best,q);
+ nT++;
+ }
+ ++it;
+ }
+ avg /= nT;
+}
+
+
+bool computeElementShapes (GFace *gf, double &worst, double &avg, double &best, int &nT, int &greaterThan )
+{
+ worst = 1.e22;
+ avg = 0.0;
+ best = 0.0;
+ nT = 0;
+ greaterThan = 0;
+ for (int i=0;i<gf->triangles.size();i++)
+ {
+ double q = qmTriangle(gf->triangles[i],QMTRI_RHO);
+ if (q>.9)greaterThan++;
+ avg+=q;
+ worst = std::min(worst,q);
+ best = std::max(best,q);
+ nT++;
+ }
+ avg /= nT;
+}
+
+
void smoothVertexPass ( GFace *gf, BDS_Mesh &m, int &nb_smooth)
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
@@ -563,11 +623,16 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
double L = 1.e22;
+ int ne = 0;
while(it!=ite){
double l = (*it)->length();
- if (l<L && (*it)->g && (*it)->g->classif_degree == 1)L=l;
+ if ((*it)->g && (*it)->g->classif_degree == 1){
+ L=std::min(L,l);
+ ne++;
+ }
++it;
}
+ if (!ne) L = 1.e22;
if(!CTX.mesh.constrained_bgmesh)
(*itp)->lc() = L;
(*itp)->lcBGM() = L;
@@ -581,6 +646,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
double t_spl=0, t_sw=0,t_col=0,t_sm=0;
const double MINE_ = 0.67, MAXE_=1.4;
+
+ double mesh_quality = 0;
while (1)
{
// we count the number of local mesh modifs.
@@ -622,7 +689,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
splitEdgePass ( gf, m, maxE, nb_split);
//saturateEdgePass ( gf, m, maxE, nb_split);
clock_t t2 = clock();
- // swapEdgePass ( gf, m, nb_swap);
+ swapEdgePass ( gf, m, nb_swap);
clock_t t3 = clock();
collapseEdgePass ( gf, m, minE , MAXNP, nb_collaps);
clock_t t4 = clock();
@@ -634,7 +701,6 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
clock_t t7 = clock();
// clean up the mesh
-
t_spl += (double)(t2-t1)/CLOCKS_PER_SEC;
t_sw += (double)(t3-t2)/CLOCKS_PER_SEC;
t_sw += (double)(t5-t4)/CLOCKS_PER_SEC;
@@ -642,11 +708,15 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
t_col += (double)(t4-t3)/CLOCKS_PER_SEC;
t_sm += (double)(t6-t5)/CLOCKS_PER_SEC;
+ double smallest, longest, old_mesh_quality = mesh_quality;int nE,NG;
+ computeMeshSizeFieldAccuracy (gf, m, mesh_quality, smallest, longest,nE,NG);
+ double worst,avg,best; int nT,GT; computeElementShapes (gf, m, worst, avg,best,nT,GT);
m.cleanup();
IT++;
- Msg(DEBUG1," iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : %6d splits, %6d swaps, %6d collapses, %6d moves",IT,minL,minE,maxL,maxE,nb_split,nb_swap,nb_collaps,nb_smooth);
+ 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;
}
double t_total = t_spl + t_sw + t_col + t_sm;
@@ -1139,9 +1209,9 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
if (!AlgoDelaunay2D ( gf ))
{
RefineMesh (gf,*m, CTX.mesh.refine_steps);
-// OptimizeMesh(gf, *m, 2);
-// RefineMesh (gf,*m, -CTX.mesh.refine_steps);
-// OptimizeMesh(gf, *m, 2);
+ OptimizeMesh(gf, *m, 2);
+ RefineMesh (gf,*m, -CTX.mesh.refine_steps);
+
if (gf->meshAttributes.recombine)
{
m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
@@ -1154,6 +1224,9 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
// outputScalarField(m->triangles, name,0);
// fill the small gmsh structures
+ computeMeshSizeFieldAccuracy (gf,*m, gf->meshStatistics.efficiency_index,gf->meshStatistics.longest_edge_length,gf->meshStatistics.smallest_edge_length,gf->meshStatistics.nbEdge,gf->meshStatistics.nbGoodLength);
+
+ gf->meshStatistics.status = GFace::DONE;
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
while (itp != m->points.end())
@@ -1216,6 +1289,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
delete [] U_;
delete [] V_;
+ computeElementShapes (gf, gf->meshStatistics.worst_element_shape,gf->meshStatistics.average_element_shape,gf->meshStatistics.best_element_shape,gf->meshStatistics.nbTriangle,gf->meshStatistics.nbGoodQuality);
+
return true;
}
@@ -1239,7 +1314,7 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
SBoundingBox3d &bbox,
BDS_Mesh *m,
std::map<BDS_Point*,MVertex*> &recover_map,
- int &count)
+ int &count, double tol)
{
// for each edge, we build a list of points that
@@ -1249,7 +1324,9 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
// printf("new edge loop\n");
std::map<GEntity*,std::vector<SPoint2> > meshes;
- std::map<GEntity*,std::vector<SPoint2> > meshes_seam;
+ std::map<GEntity*,std::vector<SPoint2> > meshes_seam;
+
+ result.clear();
GEdgeLoop::iter it = gel.begin();
@@ -1294,8 +1371,6 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
SPoint2 last_coord(0,0);
int counter = 0;
- double tol = 1.e-7;
-
while (unordered.size())
{
std::list<GEdgeSigned>::iterator it = unordered.begin();
@@ -1325,7 +1400,7 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
{
SPoint2 first_coord = mesh1d[0];
double d = dist2(last_coord,first_coord);
- // printf("d = %12.5E %d\n",d, coords.size());
+ // printf("d = %12.5E %d\n",d, coords.size());
if (d < tol)
{
coords.clear();
@@ -1336,10 +1411,10 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
}
SPoint2 first_coord_reversed = mesh1d_reversed[0];
double d_reversed = dist2(last_coord,first_coord_reversed);
- // printf("d_r = %12.5E\n",d_reversed);
+ // printf("d_r = %12.5E\n",d_reversed);
if (d_reversed < tol)
{
- // printf("d_r = %12.5E\n",d_reversed);
+ // printf("d_r = %12.5E\n",d_reversed);
coords.clear();
coords = mesh1d_reversed;
found = (GEdgeSigned(-1,ge));
@@ -1351,7 +1426,7 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
SPoint2 first_coord_seam = mesh1d_seam[0];
SPoint2 first_coord_seam_reversed = mesh1d_seam_reversed[0];
double d_seam = dist2(last_coord,first_coord_seam);
- // printf("d_seam = %12.5E\n",d_seam);
+ // printf("d_seam = %12.5E\n",d_seam);
if (d_seam < tol)
{
coords.clear();
@@ -1361,7 +1436,7 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
goto Finalize;
}
double d_seam_reversed = dist2(last_coord,first_coord_seam_reversed);
- // printf("d_seam_reversed = %12.5E\n",d_seam_reversed);
+ // printf("d_seam_reversed = %12.5E\n",d_seam_reversed);
if (d_seam_reversed < tol)
{
coords.clear();
@@ -1377,8 +1452,6 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
}
Finalize:
-
-
if (coords.size() == 0)return false;
std::vector<MVertex*> edgeLoop;
@@ -1455,27 +1528,26 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
{
- // if(gf->tag() != 6)return true;
-
+ // if (gf->tag() != 32) return true;
std::map<BDS_Point*,MVertex*> recover_map;
Range<double> rangeU = gf->parBounds(0);
Range<double> rangeV = gf->parBounds(1);
- const double du = rangeU.high() -rangeU.low();
- const double dv = rangeV.high() -rangeV.low();
+ double du = rangeU.high() -rangeU.low();
+ double dv = rangeV.high() -rangeV.low();
const double LC2D = sqrt ( du*du + dv*dv );
- // printf("LC2D %g (%g,%g), (%g,%g)\n",LC2D,rangeU.high(),rangeU.low(),rangeV.high(),rangeV.low());
+ // printf("LC2D %g du %g (%g %g) dv %g\n",LC2D,du,rangeU.high(),rangeU.low(),dv);
// Buid a BDS_Mesh structure that is convenient for doing the actual meshing procedure
BDS_Mesh *m = new BDS_Mesh;
- m->scalingU = fabs(du);
- m->scalingV = fabs(dv);
- SCALINGU = m->scalingU;
- SCALINGV = m->scalingV;
+ m->scalingU = fabs(du);
+ m->scalingV = fabs(dv);
+ SCALINGU = m->scalingU;
+ SCALINGV = m->scalingV;
std::vector< std::vector<BDS_Point* > > edgeLoops_BDS;
SBoundingBox3d bbox;
int nbPointsTotal = 0;
@@ -1483,7 +1555,14 @@ 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)==false)return false;
+ // if(buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-7*LC2D)==false)
+ // if(buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-6*LC2D)==false)
+ if(buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-7*LC2D)==false)
+ {
+ gf->meshStatistics.status = GFace::FAILED;
+ Msg(GERROR,"The 1D Mesh seems not to be forming a closed loop");
+ return false;
+ }
edgeLoops_BDS.push_back(edgeLoop_BDS);
}
}
@@ -1554,7 +1633,6 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
free (doc.delaunay);
-
// Recover the boundary edges
// and compute characteristic lenghts using mesh edge spacing
@@ -1593,6 +1671,7 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
if (!e)
{
Msg(GERROR,"impossible to recover the edge %d %d",edgeLoop_BDS[j]->iD,edgeLoop_BDS[(j+1)%edgeLoop_BDS.size()]->iD);
+ gf->meshStatistics.status = GFace::FAILED;
SCALINGU = SCALINGV = 1;
return false;
}
@@ -1693,10 +1772,13 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
{
m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
}
+ // compute mesh statistics
+ computeMeshSizeFieldAccuracy (gf,*m, gf->meshStatistics.efficiency_index,gf->meshStatistics.longest_edge_length,gf->meshStatistics.smallest_edge_length,gf->meshStatistics.nbEdge,gf->meshStatistics.nbGoodLength);
+ gf->meshStatistics.status = GFace::DONE;
}
- // fill the small gmsh structures
-
+
+ // fill the small gmsh structures
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
while (itp != m->points.end())
@@ -1760,6 +1842,7 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
delete m;
SCALINGU = SCALINGV = 1;
+ computeElementShapes (gf, gf->meshStatistics.worst_element_shape,gf->meshStatistics.average_element_shape,gf->meshStatistics.best_element_shape,gf->meshStatistics.nbTriangle,gf->meshStatistics.nbGoodQuality);
return true;
}
@@ -1775,6 +1858,10 @@ void deMeshGFace::operator() (GFace *gf)
for (unsigned int i=0;i<gf->quadrangles.size();i++) delete gf->quadrangles[i];
gf->quadrangles.clear();
gf->deleteVertexArrays();
+
+ gf->meshStatistics.status = GFace::PENDING;
+ gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
+
}
void meshGFace::operator() (GFace *gf)
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 1c5cad48be07250633a89b628a5a9bfe33c04a4d..5af64a19655bd5c0ba08414cfa52fc37c6d7d482 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.33 2007-09-12 20:14:35 geuzaine Exp $
+// $Id: meshGRegion.cpp,v 1.34 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -229,6 +229,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// now do insertion of points
insertVerticesInRegion(gr);
+ Msg(INFO, "Gmsh 3D Delaunay has generated %d tets", gr->tetrahedra.size());
#endif
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 7018ced6575d705915a8090e5791e083fb86acb3..ff243ebb1e6ec2fb0071debb858be9eba3877986 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.18 2007-09-04 13:47:02 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.19 2007-11-04 21:03:17 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -130,6 +130,7 @@ void recurFindCavity ( std::list<faceXtet> & shell,
bool insertVertex (MVertex *v ,
MTet4 *t ,
+ MTet4Factory &myFactory,
std::set<MTet4*,compareTet4Ptr> &allTets,
std::vector<double> & vSizes)
{
@@ -211,7 +212,7 @@ bool insertVertex (MVertex *v ,
// it->v[2]->y(),
// it->v[2]->z());
- MTet4 *t4 = new MTet4 ( tr , vSizes);
+ MTet4 *t4 = myFactory.Create ( tr , vSizes);
t4->setOnWhat(t->onWhat());
newTets[k++]=t4;
// all new tets are pushed front in order to
@@ -237,14 +238,21 @@ bool insertVertex (MVertex *v ,
{
connectTets ( new_cavity.begin(),new_cavity.end() );
allTets.insert(newTets,newTets+shell.size());
- //connectTets ( allTets.begin(),allTets.end() );
+
+// ittet = cavity.begin();
+// ittete = cavity.end();
+// while ( ittet != ittete )
+// {
+// myFactory.Free (*ittet);
+// ++ittet;
+// }
delete [] newTets;
return true;
}
// The cavity is NOT star shaped
else
{
- for (unsigned int i=0;i<shell.size();i++)delete newTets[i];
+ for (unsigned int i=0;i<shell.size();i++)myFactory.Free(newTets[i]);
delete [] newTets;
ittet = cavity.begin();
ittete = cavity.end();
@@ -410,10 +418,12 @@ void recur_classify ( MTet4 *t ,
void insertVerticesInRegion (GRegion *gr)
{
- std::set<MTet4*,compareTet4Ptr> allTets;
- std::set<MTet4*> voidTets;
+ printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",sizeof(MTet4),sizeof(MTetrahedron), sizeof(MVertex));
+
std::map<MVertex*,double> vSizesMap;
std::vector<double> vSizes;
+ MTet4Factory myFactory(1600000) ;
+ std::set<MTet4*,compareTet4Ptr> &allTets = myFactory.getAllTets();
for (unsigned int i=0;i<gr->tetrahedra.size();i++)setLcs ( gr->tetrahedra[i] , vSizesMap);
@@ -425,14 +435,14 @@ void insertVerticesInRegion (GRegion *gr)
}
for (unsigned int i=0;i<gr->tetrahedra.size();i++)
- allTets.insert ( new MTet4 ( gr->tetrahedra[i] ,vSizes ) );
+ allTets.insert(myFactory.Create ( gr->tetrahedra[i] ,vSizes ));
gr->tetrahedra.clear();
connectTets ( allTets.begin(), allTets.end() );
// classify the tets on the right region
// Msg (INFO,"reclassifying %d tets",allTets.size());
- for (std::set<MTet4*,compareTet4Ptr>::iterator it = allTets.begin();it!=allTets.end();++it)
+ for (MTet4Factory::iterator it = allTets.begin();it!=allTets.end();++it)
{
if (!(*it)->onWhat())
{
@@ -452,12 +462,12 @@ void insertVerticesInRegion (GRegion *gr)
}
- for (std::set<MTet4*,compareTet4Ptr>::iterator it2 = allTets.begin();it2!=allTets.end();++it2)
+ for (MTet4Factory::iterator it = allTets.begin();it!=allTets.end();++it)
{
- (*it2)->setNeigh(0,0);
- (*it2)->setNeigh(1,0);
- (*it2)->setNeigh(2,0);
- (*it2)->setNeigh(3,0);
+ (*it)->setNeigh(0,0);
+ (*it)->setNeigh(1,0);
+ (*it)->setNeigh(2,0);
+ (*it)->setNeigh(3,0);
}
connectTets ( allTets.begin(), allTets.end() );
Msg(DEBUG,"All %d tets were connected",allTets.size());
@@ -472,13 +482,12 @@ void insertVerticesInRegion (GRegion *gr)
Msg(GERROR, "No tetrahedra in region %d", gr->tag());
break;
}
-
+
MTet4 *worst = *allTets.begin();
-
+
if (worst->isDeleted())
{
- delete worst->tet();
- delete worst;
+ myFactory.Free(worst);
allTets.erase(allTets.begin());
// Msg(INFO,"Worst tet is deleted");
}
@@ -503,41 +512,48 @@ void insertVerticesInRegion (GRegion *gr)
double lc = std::min(lc1,BGM_MeshSize(gr,0,0,center[0],center[1],center[2]));
vSizes.push_back(lc);
// compute mesh spacing there
- if (!insertVertex ( v , worst, allTets,vSizes))
+ if (!insertVertex ( v , worst, myFactory,allTets,vSizes))
{
- allTets.erase(allTets.begin());
- worst->forceRadius(0);
- allTets.insert(worst);
+ myFactory.changeTetRadius(allTets.begin(),0.0);
delete v;
}
else
- {
- v->onWhat()->mesh_vertices.push_back(v);
- }
+ v->onWhat()->mesh_vertices.push_back(v);
}
else
+ myFactory.changeTetRadius(allTets.begin(),0.0);
+ }
+ // Normally, a tet mesh contains about 5 to 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
+ if(allTets.size() > 7 * vSizes.size())
+ {
+ int n1 = allTets.size();
+ for(std::set<MTet4*,compareTet4Ptr>::iterator itd = allTets.begin() ; itd !=allTets.end() ; ++itd)
{
- // Msg(INFO,"Point is outside");
- allTets.erase(allTets.begin());
- worst->forceRadius(0);
- allTets.insert(worst);
+ if ((*itd)->isDeleted())
+ {
+ myFactory.Free((*itd));
+ std::set<MTet4*,compareTet4Ptr>::iterator itdb = itd;
+ itd++;
+ allTets.erase(itdb);
+ }
}
+ Msg(INFO,"cleaning up the memory %d -> %d ",n1,allTets.size());
}
}
-
+
while (1)
{
if (allTets.begin() == allTets.end() ) break;
MTet4 *worst = *allTets.begin();
- if (worst->isDeleted())
- {
- delete worst->tet();
- }
- else
+ if (!worst->isDeleted())
{
worst->onWhat()->tetrahedra.push_back(worst->tet());
+ worst->tet() = 0;
}
- delete worst;
+ myFactory.Free(worst);
allTets.erase(allTets.begin());
}
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 11e972fa40878a915f9f95f3afd13f0fc30142a8..9309625ea2c4f683b443d0ed7d5899f4b21bf249 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -23,28 +23,46 @@
#include "MElement.h"
#include <list>
#include <set>
+#include <stack>
+//#define _GMSH_PRE_ALLOCATE_STRATEGY_ 1
class GRegion;
+class MTet4Factory;
+
class MTet4
{
+ friend class MTet4Factory;
+ // a total of 36 Bytes for a MTet4
+ // Normally it should take 36 MByte in excess per million of tets
+ // 36 MB for 10e6 Tets
+ // Each MTetrahedron has a size of 28 Bytes
+ // The total memory required for a tet is 64 Bytes
+ // i.e. 64 MB per 10e6 Tet
+ // Yet, we store also a rb tree containing all pointers
+ // sorted with respect to tet radius.
+ // each bucket of the tree contain 4 pointers, i.e. 16 Bytes plus
+ // the data. We have therefor an extra cost of 20 Bytes/Tet
+ // i.e. 84 MB/ million tets
+ // A MVertex has a cost of 40 Bytes and there are about 200000 of them
+ // per million tet: a new cost of 8MB/10e6 tet : total 92 MB/10e6 tet
+ // (I observe 160M MB !!!)
+
bool deleted;
double circum_radius;
MTetrahedron *base;
MTet4 *neigh[4];
GRegion *gr;
- public :
- inline GRegion * onWhat () const {return gr;}
- inline void setOnWhat (GRegion *g) {gr=g;}
-
- bool isDeleted () const {return deleted;}
- void forceRadius (double r){circum_radius=r;}
- inline double getRadius ()const {return circum_radius;}
-
+ ~MTet4 (){}
+ MTet4 () : deleted(false), base (0), gr(0),circum_radius(0.0)
+ {
+ neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
+ }
- MTet4 ( MTetrahedron * t, std::vector<double> & sizes) : deleted(false), base (t), gr(0)
+ void setup ( MTetrahedron * t, std::vector<double> & sizes)
{
+ base = t;
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
double center[3];
base->circumcenter(center);
@@ -58,9 +76,22 @@ class MTet4
sizes [base->getVertex(2)->getNum()]+
sizes [base->getVertex(3)->getNum()]);
circum_radius /= lc;
+ deleted = false;
}
+ public :
+
+ inline GRegion * onWhat () const {return gr;}
+ inline void setOnWhat (GRegion *g) {gr=g;}
+
+ bool isDeleted () const {return deleted;}
+ void forceRadius (double r){circum_radius=r;}
+ inline double getRadius ()const {return circum_radius;}
+
+
+
inline MTetrahedron * tet() const {return base;}
+ inline MTetrahedron * &tet() {return base;}
inline void setNeigh (int iN , MTet4 *n) {neigh[iN]=n;}
inline MTet4 *getNeigh (int iN ) const {return neigh[iN];}
int inCircumSphere ( const double *p ) const;
@@ -77,6 +108,7 @@ class MTet4
double getVolume () const { return base -> getVolume() ; };
inline void setDeleted (bool d)
{
+ // circum_radius = d ? fabs(circum_radius) : fabs(circum_radius)
deleted = d;
}
inline bool assertNeigh() const
@@ -110,4 +142,81 @@ class compareTet4Ptr
}
};
+class MTet4Factory
+{
+public:
+ typedef std::set<MTet4*,compareTet4Ptr> container;
+ typedef container::iterator iterator;
+private:
+ container allTets;
+#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
+ MTet4* allSlots;
+ int s_last, s_alloc;
+ std::stack<MTet4*> emptySlots;
+
+ inline MTet4 * getANewSlot()
+ {
+ if (s_last >= s_alloc)return 0;
+ MTet4 * t = &(allSlots[s_last]);
+ s_last++;
+ return t;
+ }
+ inline MTet4 * getAnEmptySlot()
+ {
+ if(!emptySlots.empty())
+ {
+ MTet4* t = emptySlots.top();
+ emptySlots.pop();
+ return t;
+ }
+ return getANewSlot();
+ };
+#endif
+
+
+ public :
+ MTet4Factory (int _size = 1000000)
+ {
+#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
+ s_last=0; s_alloc=_size;
+ allSlots = new MTet4 [s_alloc];
+#endif
+ }
+ ~MTet4Factory () {
+#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
+ delete [] allSlots;
+#endif
+ }
+ MTet4 * Create (MTetrahedron * t, std::vector<double> & sizes)
+ {
+#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
+ MTet4 * t4 = getAnEmptySlot();
+#else
+ MTet4 * t4 = new MTet4;
+#endif
+ t4->setup(t,sizes);
+ return t4;
+ }
+ void Free (MTet4* t)
+ {
+ if (t->tet()) delete t->tet();
+ t->tet() = 0;
+#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
+ emptySlots.push(t);
+ t->setDeleted(true);
+#else
+ delete t;
+#endif
+ }
+
+ void changeTetRadius ( iterator it , double r)
+ {
+ allTets.erase(it);
+ (*it)->forceRadius(r);
+ allTets.insert(*it);
+ }
+ container & getAllTets () {return allTets;}
+
+};
+
#endif
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a97e3d811d579486668d52a8594d3e6dfe033cea
--- /dev/null
+++ b/Mesh/qualityMeasures.cpp
@@ -0,0 +1,66 @@
+#include "qualityMeasures.h"
+#include "BDS.h"
+#include "MVertex.h"
+#include "MElement.h"
+#include "Numeric.h"
+
+
+double qmTriangle ( const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3, const gmshQualityMeasure4Triangle &cr)
+{
+ return qmTriangle (p1->X,p1->Y,p1->Z,p2->X,p2->Y,p2->Z,p3->X,p3->Y,p3->Z,cr);
+}
+
+double qmTriangle ( BDS_Face *t, const gmshQualityMeasure4Triangle &cr)
+{
+ BDS_Point *n[4];
+ t->getNodes(n);
+ return qmTriangle (n[0],n[1],n[2],cr);
+}
+
+double qmTriangle ( MTriangle*t, const gmshQualityMeasure4Triangle &cr)
+{
+ return qmTriangle (t->getVertex(0),t->getVertex(1),t->getVertex(2),cr);
+}
+
+double qmTriangle ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const gmshQualityMeasure4Triangle &cr)
+{
+ return qmTriangle (v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(),v3->x(),v3->y(),v3->z(),cr);
+}
+
+// Triangle abc
+// quality is between 0 and 1
+
+double qmTriangle ( const double &xa, const double &ya, const double &za,
+ const double &xb, const double &yb, const double &zb,
+ const double &xc, const double &yc, const double &zc,
+ const gmshQualityMeasure4Triangle &cr)
+{
+ double quality;
+ switch(cr)
+ {
+ case QMTRI_RHO:
+ {
+ // quality = rho / R = 2 * inscribed radius / circumradius
+
+ double a [3] = {xc-xb,yc-yb,zc-zb};
+ double b [3] = {xa-xc,ya-yc,za-zc};
+ double c [3] = {xb-xa,yb-ya,zb-za};
+
+ const double la = norme (a);
+ const double lb = norme (b);
+ const double lc = norme (c);
+
+ double pva [3]; prodve (b,c,pva); const double sina = norm3 (pva);
+ double pvb [3]; prodve (c,a,pvb); const double sinb = norm3 (pvb);
+ double pvc [3]; prodve (a,b,pvc); const double sinc = norm3 (pvc);
+
+ if (sina == 0.0 && sinb == 0.0 && sinc == 0.0) quality = 0.0;
+ else quality = 2 * (2*sina*sinb*sinc/(sina + sinb + sinc) );
+ }
+ break;
+ default:
+ throw;
+ }
+ return quality;
+}
+
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
new file mode 100644
index 0000000000000000000000000000000000000000..02429b4fb89da054386de2dde1ab5887fa3bd27e
--- /dev/null
+++ b/Mesh/qualityMeasures.h
@@ -0,0 +1,20 @@
+#ifndef _QUALITY_MEASURES_H_
+#define _QUALITY_MEASURES_H_
+
+class BDS_Point;
+class BDS_Face;
+class MVertex;
+class MTriangle;
+enum gmshQualityMeasure4Triangle {QMTRI_RHO};
+
+double qmTriangle ( MTriangle *f, const gmshQualityMeasure4Triangle &cr);
+double qmTriangle ( BDS_Face *f, const gmshQualityMeasure4Triangle &cr);
+double qmTriangle ( const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3, const gmshQualityMeasure4Triangle &cr);
+double qmTriangle ( const MVertex *v1, const MVertex *v2, const MVertex *v3, const gmshQualityMeasure4Triangle &cr);
+double qmTriangle ( const double *d1, const double *d2, const double *d3, const gmshQualityMeasure4Triangle &cr);
+double qmTriangle ( 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 gmshQualityMeasure4Triangle &cr);
+
+#endif