diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 8abf311b23573c64579c3dc08a77eaf891548a8e..bc598b00a2d705be0c2ac26b46c2b9a5b67ed8db 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.43 2008-02-22 21:09:01 geuzaine Exp $
+// $Id: meshGRegion.cpp,v 1.44 2008-03-06 14:19:01 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -96,27 +96,25 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
tetgenio::facet *f = &in.facetlist[I];
- tetgenio::init(f);
+ tetgenio::init(f);
f->numberofholes = 0;
f->numberofpolygons = 1;
f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
tetgenio::polygon *p = &f->polygonlist[0];
- tetgenio::init(p);
+ tetgenio::init(p);
p->numberofvertices = 3;
p->vertexlist = new int[p->numberofvertices];
p->vertexlist[0] = t->getVertex(0)->getNum();
p->vertexlist[1] = t->getVertex(1)->getNum();
p->vertexlist[2] = t->getVertex(2)->getNum();
- in.facetmarkerlist[I] = gf->tag();
+ in.facetmarkerlist[I] = gf->tag();
++I;
}
++it;
}
}
-void TransferTetgenMesh(GRegion *gr,
- tetgenio &in,
- tetgenio &out,
+void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
std::vector<MVertex*> &numberedV)
{
for(int i = numberedV.size(); i < out.numberofpoints; i++){
@@ -141,13 +139,13 @@ void TransferTetgenMesh(GRegion *gr,
gf->triangles.clear();
gf->deleteVertexArrays();
++it;
- }
+ }
- // TODO: re-create 1D mesh
+ // TODO: re-create 1D mesh
for(int i = 0; i < out.numberofedges; i++){
MVertex *v[2];
- v[0] = numberedV[out.edgelist[i * 2 + 0] -1];
- v[1] = numberedV[out.edgelist[i * 2 + 1] -1];
+ v[0] = numberedV[out.edgelist[i * 2 + 0] - 1];
+ v[1] = numberedV[out.edgelist[i * 2 + 1] - 1];
}
// re-create the triangular meshes FIXME: this can lead to hanging
@@ -212,38 +210,40 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
allFaces.push_back(*it);
gr->set(allFaces);
- // mesh with tetgen, possibly changing the mesh on boundaries
- tetgenio in, out;
- std::vector<MVertex*> numberedV;
- char opts[128];
- buildTetgenStructure(gr, in, numberedV);
- sprintf(opts, "pe%c", (CTX.verbosity < 3) ? 'Q': (CTX.verbosity > 6)? 'V': '\0');
- try{
- tetrahedralize(opts, &in, &out);
- }
- catch (int error){
- Msg (WARNING, "Self intersecting Surface Mesh, computing intersections "
- "(this could take a while)");
- sprintf(opts, "dV");
+ // mesh with tetgen, possibly changing the mesh on boundaries (leave
+ // this in block, so in/out are destroyed before we refine the mesh)
+ {
+ tetgenio in, out;
+ std::vector<MVertex*> numberedV;
+ char opts[128];
+ buildTetgenStructure(gr, in, numberedV);
+ sprintf(opts, "pe%c", (CTX.verbosity < 3) ? 'Q': (CTX.verbosity > 6)? 'V': '\0');
try{
- tetrahedralize(opts, &in, &out);
- Msg(INFO,"%d faces self-intersect",out.numberoftrifaces);
- for (int i=0;i<out.numberoftrifaces;i++)
- {
+ tetrahedralize(opts, &in, &out);
+ }
+ catch (int error){
+ Msg (WARNING, "Self intersecting Surface Mesh, computing intersections "
+ "(this could take a while)");
+ sprintf(opts, "dV");
+ try{
+ tetrahedralize(opts, &in, &out);
+ Msg(INFO,"%d faces self-intersect",out.numberoftrifaces);
+ for (int i = 0; i < out.numberoftrifaces; i++){
Msg(INFO,"face (%d %d %d) on model face %d",
numberedV[out.trifacelist[i * 3 + 0] - 1]->getNum(),
numberedV[out.trifacelist[i * 3 + 1] - 1]->getNum(),
numberedV[out.trifacelist[i * 3 + 2] - 1]->getNum(),
out.trifacemarkerlist[i]);
}
+ }
+ catch (int error2){
+ Msg(GERROR, "Surface Mesh is wrong, cannot do the 3D mesh");
+ }
+ gr->set(faces);
+ return;
}
- catch (int error2){
- Msg (GERROR, "Surface Mesh is wrong, cannot do the 3D mesh");
- }
- gr->set(faces);
- return;
+ TransferTetgenMesh(gr, in, out, numberedV);
}
- TransferTetgenMesh(gr, in, out, numberedV);
// sort triangles in all model faces in order to be able to search in vectors
std::list<GFace*>::iterator itf = allFaces.begin();
@@ -257,7 +257,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
gr->set(faces);
// now do insertion of points
- insertVerticesInRegion(gr);
+ insertVerticesInRegion(gr);
Msg(INFO, "Gmsh 3D Delaunay has generated %d tets", gr->tetrahedra.size());
#endif
}
@@ -306,13 +306,13 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator it = faces.begin();
while(it != faces.end()){
- GFace *gf = (*it);
+ GFace *gf = (*it);
for(unsigned int i = 0; i< gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
int tmp[3];
tmp[0] = t->getVertex(0)->getNum();
tmp[1] = t->getVertex(1)->getNum();
- tmp[2] = t->getVertex(2)->getNum();
+ tmp[2] = t->getVertex(2)->getNum();
Ng_AddSurfaceElement(ngmesh, NG_TRIG, tmp);
}
++it;
@@ -326,8 +326,8 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
int tmp[4];
tmp[0] = t->getVertex(0)->getNum();
tmp[1] = t->getVertex(1)->getNum();
- tmp[2] = t->getVertex(2)->getNum();
- tmp[3] = t->getVertex(3)->getNum();
+ tmp[2] = t->getVertex(2)->getNum();
+ tmp[3] = t->getVertex(3)->getNum();
Ng_AddVolumeElement(ngmesh, NG_TET, tmp);
}
}
@@ -416,10 +416,10 @@ int intersect_line_triangle(double X[3], double Y[3], double Z[3] ,
b[1] = P[1] - Y[0];
b[2] = P[2] - Z[0];
- if(!sys3x3_with_tol(mat, b, res, &det))
+ if(!sys3x3_with_tol(mat, b, res, &det))
return 0;
- // Msg(INFO,"going there %g %g %g",res[0],res[1],res[2]);
+ // Msg(INFO, "going there %g %g %g", res[0], res[1], res[2]);
if(res[0] >= eps_prec && res[0] <= 1.0 - eps_prec &&
res[1] >= eps_prec && res[1] <= 1.0 - eps_prec &&
@@ -454,7 +454,7 @@ void meshNormalsPointOutOfTheRegion(GRegion *gr)
setRand(rrr);
while(it != faces.end()){
- GFace *gf = (*it);
+ GFace *gf = (*it);
int nb_intersect = 0;
for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
@@ -464,7 +464,7 @@ void meshNormalsPointOutOfTheRegion(GRegion *gr)
double P[3] = {(X[0]+X[1]+X[2])/3., (Y[0]+Y[1]+Y[2])/3., (Z[0]+Z[1]+Z[2])/3.};
double v1[3] = {X[0]-X[1], Y[0]-Y[1], Z[0]-Z[1]};
double v2[3] = {X[2]-X[1], Y[2]-Y[1], Z[2]-Z[1]};
- double N[3] ;
+ double N[3];
prodve(v1, v2, N);
norme(v1);
norme(v2);
@@ -611,7 +611,6 @@ void optimizeMeshGRegionGmsh::operator() (GRegion *gr)
// import mesh into netgen, including volume tets
gmshOptimizeMesh (gr, QMTET_2);
-
}
bool buildFaceSearchStructure(GModel *model, fs_cont &search)
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 79409189958ce4212eecd84793d8bb308f2855d7..2c0e68dfcfcad2be0d01db3c8aa76e60d35f6c61 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.38 2008-03-04 08:51:14 geuzaine Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.39 2008-03-06 14:19:01 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -76,7 +76,7 @@ struct faceXtet{
};
template <class ITER>
-void connectTets ( ITER beg, ITER end)
+void connectTets(ITER beg, ITER end)
{
std::set<faceXtet> conn;
while (beg != end){
@@ -131,11 +131,9 @@ void recurFindCavity(std::list<faceXtet> & shell,
}
}
-bool insertVertex(MVertex *v,
- MTet4 *t,
- MTet4Factory &myFactory,
- std::set<MTet4*,compareTet4Ptr> &allTets,
- std::vector<double> & vSizes)
+bool insertVertex(MVertex *v, MTet4 *t, MTet4Factory &myFactory,
+ std::set<MTet4*, compareTet4Ptr> &allTets,
+ std::vector<double> &vSizes)
{
std::list<faceXtet> shell;
std::list<MTet4*> cavity;
@@ -175,7 +173,7 @@ bool insertVertex(MVertex *v,
++ittet;
}
// fprintf(ff2,"};\n");
-// fclose (ff2);
+// fclose(ff2);
// Msg(INFO,"cavity of size %d volume %g",cavity.size(),oldVolume);
// create new tetrahedron using faces that are
// on the border of the cavity
@@ -186,7 +184,7 @@ bool insertVertex(MVertex *v,
// char name[245];
// sprintf(name,"test%d.pos",III);
-// FILE *ff = fopen (name,"w");
+// FILE *ff = fopen(name,"w");
// fprintf(ff,"View\"test\"{\n");
MTet4** newTets = new MTet4*[shell.size()];;
@@ -195,7 +193,7 @@ bool insertVertex(MVertex *v,
std::list<faceXtet>::iterator it = shell.begin();
while (it != shell.end()){
- MTetrahedron *tr = new MTetrahedron (it->v[0], it->v[1], it->v[2], v);
+ MTetrahedron *tr = new MTetrahedron(it->v[0], it->v[1], it->v[2], v);
// Msg(INFO,"shell %d %d %d",it->v[0]->getNum(),it->v[1]->getNum(),it->v[2]->getNum());
// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
// it->v[0]->x(),
@@ -229,8 +227,8 @@ bool insertVertex(MVertex *v,
// Msg(INFO,"new cavity of vol %g (%d boundaries)",newVolume,shell.size());
// OK, the cavity is star shaped
if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume){
- connectTets ( new_cavity.begin(),new_cavity.end());
- allTets.insert(newTets,newTets+shell.size());
+ connectTets(new_cavity.begin(), new_cavity.end());
+ allTets.insert(newTets, newTets + shell.size());
// ittet = cavity.begin();
// ittete = cavity.end();
@@ -254,7 +252,7 @@ bool insertVertex(MVertex *v,
}
}
-static void setLcs(MTetrahedron *t, std::map<MVertex*,double> &vSizes)
+static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes)
{
for (int i = 0; i < 4; i++){
for (int j = i + 1; j < 4; j++){
@@ -319,12 +317,9 @@ GRegion *getRegionFromBoundingFaces(GModel *model,
return 0;
}
-void recur_classify(MTet4 *t ,
- std::list<MTet4*> &theRegion,
- std::set<GFace *> &faces_bound,
- GRegion *bidon ,
- GModel *model,
- const fs_cont &search)
+void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
+ std::set<GFace*> &faces_bound, GRegion *bidon,
+ GModel *model, const fs_cont &search)
{
if (!t) Msg (GERROR,"a tet is not connected by a boundary face");
if (t->onWhat()) return; // should never return here...
@@ -368,7 +363,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets;
- for(unsigned int i=0;i<gr->tetrahedra.size();i++){
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
allTets.push_back(new MTet4(gr->tetrahedra[i], qm));
}
gr->tetrahedra.clear();
@@ -377,8 +372,8 @@ void adaptMeshGRegion::operator () (GRegion *gr)
double t1 = Cpu();
std::vector<MTet4*> illegals;
- const int nbRanges=10;
- int quality_ranges [nbRanges];
+ const int nbRanges = 10;
+ int quality_ranges[nbRanges];
{
double totalVolumeb = 0.0;
double worst = 1.0;
@@ -545,7 +540,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
//template <class CONTAINER, class DATA>
-void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
+void gmshOptimizeMesh(GRegion *gr, const gmshQualityMeasure4Tet &qm)
{
typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets;
@@ -726,19 +721,20 @@ void insertVerticesInRegion (GRegion *gr)
//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);
-
+ std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
int NUM = 0;
- for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
- it != vSizesMap.end(); ++it){
- it->first->setNum(NUM++);
- vSizes.push_back(it->second);
+
+ { // leave this in a block so the map gets deallocated directly
+ std::map<MVertex*, double> vSizesMap;
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
+ setLcs(gr->tetrahedra[i], vSizesMap);
+ for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
+ it != vSizesMap.end(); ++it){
+ it->first->setNum(NUM++);
+ vSizes.push_back(it->second);
+ }
}
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
@@ -751,7 +747,7 @@ void insertVerticesInRegion (GRegion *gr)
// Msg (INFO,"reclassifying %d tets", allTets.size());
fs_cont search;
- buildFaceSearchStructure(gr->model(), search );
+ buildFaceSearchStructure(gr->model(), search);
for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end(); ++it){
if(!(*it)->onWhat()){
@@ -812,7 +808,7 @@ void insertVerticesInRegion (GRegion *gr)
double uvw[3];
bool inside = worst->tet()->invertmapping(center, uvw);
if(inside){
- MVertex *v = new MVertex(center[0],center[1],center[2], worst->onWhat());
+ MVertex *v = new MVertex(center[0], center[1], center[2], worst->onWhat());
v->setNum(NUM++);
double lc1 =
(1 - uvw[0] - uvw[1] - uvw[2]) * vSizes[worst->tet()->getVertex(0)->getNum()] +
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 694ecc188eec813cf99fa673c5986bd54d0dcba5..3377e514e59bec584e8465dc23dba6c9e2504593 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -34,16 +34,28 @@ class GModel;
class MTet4Factory;
-// sizeof(MTet4) = 36 Bytes and sizeof(MTetrahedron) = 28 Bytes, so
-// normally it should take 36+28 = 64 MByte per million tets. We also
-// store 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 -> Extra cost of 20 Bytes/Tet, i.e., 84 MB per
-// million tets. A MVertex has a cost of 44 Bytes and there are about
-// 200000 of them per million tet, i.e., a new cost of 9MB per million
-// tets.
-
-// Grand total should be 92 MB per million tet (I observe 160M MB!)
+// Memory usage for 1 million tets:
+//
+// * sizeof(MTet4) = 36 Bytes and sizeof(MTetrahedron) = 28 Bytes
+// -> 64 MB
+// * rb tree containing all pointers sorted with respect to tet
+// radius: each bucket of the tree contains 4 pointers (16 Bytes)
+// plus the data -> 20 MB
+// * sizeof(MVertex) = 44 Bytes and there are about 200000 verts per
+// million tet -> 9MB
+// * vector of char lengths per vertex -> 1.6Mb
+// * vectors in GEntities to store the element and vertex pointers
+// -> 5Mb
+//
+// Grand total should thus be about 100 MB.
+//
+// The observed mem usage with "demos/cube.geo -clscale 0.61" is
+// 157MB. Where do the extra 57 MB come from?
+//
+// * surface mesh + all other overheads (model, etc.) is 19Mb
+// * tetgen initial mesh is about 20Mb, but it is deleted before mesh
+// refinement.
+// * ?
class MTet4
{
@@ -92,8 +104,8 @@ class MTet4
}
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 bool isDeleted () const { return deleted; }
+ inline void forceRadius (double r){ circum_radius = r; }
inline double getRadius () const { return circum_radius; }
inline double getQuality () const { return circum_radius; }
inline void setQuality (const double &q){ circum_radius = q; }
@@ -111,7 +123,7 @@ class MTet4
{
return inCircumSphere(v->x(), v->y(), v->z());
}
- double getVolume() const { return base->getVolume(); }
+ inline double getVolume() const { return base->getVolume(); }
inline void setDeleted(bool d)
{
deleted = d;
@@ -142,7 +154,7 @@ class compareTet4Ptr
{
if (a->getRadius() > b->getRadius()) return true;
if (a->getRadius() < b->getRadius()) return false;
- return a<b;
+ return a < b;
}
};
diff --git a/benchmarks/3d/periodic.geo b/benchmarks/3d/periodic.geo
index c82d1eff26f9af355931c6f4d89b5003a2f8b3af..f0cbc96243871ae038fbcb7aaad701ff9f076444 100644
--- a/benchmarks/3d/periodic.geo
+++ b/benchmarks/3d/periodic.geo
@@ -46,7 +46,7 @@ Physical Line(18) = {2,3};
Physical Surface(19) = {15};
If(use_prisms)
- Extrude Surface {15, {0.,0.,2.*R}}{Layers{nb_layers,83,1}; Recombine; };
+ Extrude Surface {15, {0.,0.,2.*R}}{ Layers{nb_layers}; Recombine; };
EndIf
If(!use_prisms)