diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index bc1350dcb4ef4cde59b1fcc20ba2e7964eac6087..f882267ababce80954ccd15f9caa41cc0569072e 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -162,6 +162,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MV
// SPoint2 pp = gf->parFromPoint(SPoint3 (v[j]->x(),v[j]->y(),v[j]->z()));
SPoint2 pp (0,0);
MFaceVertex *v1b = new MFaceVertex (v[j]->x(),v[j]->y(),v[j]->z(),gf,pp.x(),pp.y() );
+ gf->mesh_vertices.push_back(v1b);
numberedV[out.trifacelist[i * 3 + j] -1] = v1b;
delete v[j];
v[j] = v1b;
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 105e83a3ee2d94d3b80150fab6b0ac4e75924f71..57f367f3d6f46ba4ea1663a707f03377618b1033 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -3,6 +3,7 @@
#include "Numeric.h"
#include "Message.h"
#include <set>
+#include <map>
#include <algorithm>
int MTet4::inCircumSphere ( const double *p ) const
@@ -104,11 +105,10 @@ void recurFindCavity ( std::list<faceXtet> & shell,
}
}
-int III;
-
bool insertVertex (MVertex *v ,
MTet4 *t ,
- std::set<MTet4*,compareTet4Ptr> &allTets)
+ std::set<MTet4*,compareTet4Ptr> &allTets,
+ std::vector<double> & vSizes)
{
std::list<faceXtet> shell;
std::list<MTet4*> cavity;
@@ -132,7 +132,6 @@ bool insertVertex (MVertex *v ,
double oldVolume = 0;
// char name2[245];
- // sprintf(name2,"testv%d.pos",III);
// FILE *ff2 = fopen (name2,"w");
// fprintf(ff2,"View\"test\"{\n");
@@ -198,7 +197,7 @@ bool insertVertex (MVertex *v ,
// it->v[2]->y(),
// it->v[2]->z());
- MTet4 *t4 = new MTet4 ( t );
+ MTet4 *t4 = new MTet4 ( t , vSizes);
newTets[k++]=t4;
// all new tets are pushed front in order to
// ba able to destroy them if the cavity is not
@@ -242,45 +241,58 @@ bool insertVertex (MVertex *v ,
return false;
}
}
+
+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++)
+ {
+ MVertex *vi = t->getVertex(i);
+ MVertex *vj = t->getVertex(j);
+ double dx = vi->x()-vj->x();
+ double dy = vi->y()-vj->y();
+ double dz = vi->z()-vj->z();
+ double l = sqrt(dx*dx+dy*dy+dz*dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ if (iti==vSizes.end() || iti->second > l)vSizes[vi] = l;
+ if (itj==vSizes.end() || itj->second > l)vSizes[vj] = l;
+ }
+ }
+}
+
+
void insertVerticesInRegion (GRegion *gr)
{
std::set<MTet4*,compareTet4Ptr> allTets;
+ std::map<MVertex*,double> vSizesMap;
+ std::vector<double> vSizes;
- for (int i=0;i<gr->tetrahedra.size();i++)
+ for (int i=0;i<gr->tetrahedra.size();i++)setLcs ( gr->tetrahedra[i] , vSizesMap);
+
+ int NUM=0;
+ for (std::map<MVertex*,double>::iterator it = vSizesMap.begin();it!=vSizesMap.end();++it)
{
- allTets.insert ( new MTet4 ( gr->tetrahedra[i] ) );
+ it->first->setNum(NUM++);
+ vSizes.push_back(it->second);
}
+
+ for (int i=0;i<gr->tetrahedra.size();i++)
+ allTets.insert ( new MTet4 ( gr->tetrahedra[i] ,vSizes ) );
+
gr->tetrahedra.clear();
connectTets ( allTets.begin(), allTets.end() );
-// for (std::set<MTet4*,compareTet4Ptr>::iterator it = allTets.begin();it!=allTets.end();++it)
-// {
-// printf("tet %d %d %d %d neigh %p %p %p %p\n",
-// (*it)->tet()->getVertex(0)->getNum(),
-// (*it)->tet()->getVertex(1)->getNum(),
-// (*it)->tet()->getVertex(2)->getNum(),
-// (*it)->tet()->getVertex(3)->getNum(),
-// (*it)->getNeigh(0),(*it)->getNeigh(1),(*it)->getNeigh(2),(*it)->getNeigh(3));
-// }
-
-
Msg(INFO,"All %d tets were connected",allTets.size());
// here the classification should be done
- III = 0;
+ int ITER = 0;
while (1)
{
- // connectTets ( allTets.begin(), allTets.end() );
- // for (std::set<MTet4*,compareTet4Ptr>::iterator tit = allTets.begin();tit != allTets.end();++tit)
- // if(!(*tit)->assertNeigh())throw;
- // double vv = 0;
- // for (std::set<MTet4*,compareTet4Ptr>::iterator tit = allTets.begin();tit != allTets.end();++tit)
- // if (!(*tit)->isDeleted())vv+=fabs((*tit)->getVolume());
-
-
MTet4 *worst = *allTets.begin();
if (worst->isDeleted())
@@ -292,8 +304,9 @@ void insertVerticesInRegion (GRegion *gr)
}
else
{
- if(III%1000 ==0)Msg(INFO,"%d points created -- Worst tet radius is %g",gr->mesh_vertices.size(),worst->getRadius());
- if (worst->getRadius() < .015) break;
+ if(ITER++%5000 ==0)
+ Msg(INFO,"%d points created -- Worst tet radius is %g",vSizes.size(),worst->getRadius());
+ if (worst->getRadius() < 1) break;
double center[3];
worst->tet()->circumcenter(center);
double uvw[3];
@@ -301,19 +314,25 @@ void insertVerticesInRegion (GRegion *gr)
if (inside)
{
MVertex *v = new MVertex (center[0],center[1],center[2]);
+ v->setNum(NUM++);
+ double lc =
+ (1-uvw[0]-uvw[1]-uvw[2])*vSizes[worst->tet()->getVertex(0)->getNum()] +
+ (uvw[0])*vSizes[worst->tet()->getVertex(1)->getNum()] +
+ (uvw[1])*vSizes[worst->tet()->getVertex(2)->getNum()] +
+ (uvw[2])*vSizes[worst->tet()->getVertex(3)->getNum()];
+ vSizes.push_back(lc);
// compute mesh spacing there
- III++;
- if (!insertVertex ( v , worst, allTets))
+ if (!insertVertex ( v , worst, allTets,vSizes))
{
- // this one does not work
- // Msg(INFO,"Point is %g %g %g %d",uvw[0], uvw[1], uvw[2],worst->inCircumSphere ( v ));
allTets.erase(allTets.begin());
worst->forceRadius(0);
allTets.insert(worst);
delete v;
}
- else gr->mesh_vertices.push_back(v);
- // if(III == 15) break;
+ else
+ {
+ gr->mesh_vertices.push_back(v);
+ }
}
else
{
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index fcdbac13a2f68fee3343ca5e1a565f4997b6de5b..759bbacf97ad040567426eb581bc3d7b419fe043 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -17,7 +17,7 @@ class MTet4
void forceRadius (double r){circum_radius=r;}
double getRadius ()const {return circum_radius;}
- MTet4 ( MTetrahedron * t) : deleted(false), base (t)
+ MTet4 ( MTetrahedron * t, std::vector<double> & sizes) : deleted(false), base (t)
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
double center[3];
@@ -26,6 +26,12 @@ class MTet4
const double dy = base->getVertex(0)->y() - center[1];
const double dz = base->getVertex(0)->z() - center[2];
circum_radius = sqrt ( dx*dx + dy*dy + dz*dz);
+
+ double lc = 0.25*(sizes [base->getVertex(0)->getNum()]+
+ sizes [base->getVertex(1)->getNum()]+
+ sizes [base->getVertex(2)->getNum()]+
+ sizes [base->getVertex(3)->getNum()]);
+ circum_radius /= lc;
}
inline MTetrahedron * tet() const {return base;}
diff --git a/benchmarks/3d/Cube-04.geo b/benchmarks/3d/Cube-04.geo
index 357bd131cff1769057d24fcd49fd61327d15543a..47cc1819d9dbeb57935b629781f42467202c1208 100644
--- a/benchmarks/3d/Cube-04.geo
+++ b/benchmarks/3d/Cube-04.geo
@@ -7,7 +7,7 @@ Extrude Point {1, {1,0.0,0} };
Extrude Line {1, {0.0,0.0,1} };
aa[] = Extrude Surface {5, {0,1,0} };;
-Point(100) = {0.3,0.3,0.3,.02};
+Point(100) = {0.3,0.3,0.3,.1};
Extrude Point {100, {.4,0.0,0} };
Extrude Line {28, {0,0.4,0} };
Coherence;
diff --git a/benchmarks/3d/Torus.geo b/benchmarks/3d/Torus.geo
index 771b08b3aeca17b904eccd8054e31a098f6a3bd5..f4b965cb73c4e911210b5a384507fddfcf84f1db 100644
--- a/benchmarks/3d/Torus.geo
+++ b/benchmarks/3d/Torus.geo
@@ -1,4 +1,4 @@
-lc = .2;
+lc = .1;
Point(1) = {2.0,0.0,0.0,lc};
Point(2) = {2.0,1,0.0,lc};
Point(3) = {1,0,0.0,lc};