diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 49f9c3f2832e7fd9c726973e0629832dd55d7a11..fb991814faae2776ded3293fc16147d2e5be5de7 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1623,9 +1623,9 @@ void bowyerWatsonParallelograms(GFace *gf)
std::vector<MVertex*> packed;
std::vector<SMetric3> metrics;
- printf("creating the points\n");
- packingOfParallelograms(gf, packed);//, metrics);
- printf("points created\n");
+ // printf("creating the points\n");
+ packingOfParallelograms(gf, packed, metrics);
+ // printf("points created\n");
buildMeshGenerationDataStructures
(gf, AllTris, vSizes, vSizesBGM, vMetricsBGM,Us, Vs);
@@ -1668,7 +1668,7 @@ void bowyerWatsonParallelograms(GFace *gf)
}
if(1.0* AllTris.size() > 2.5 * vSizes.size()){
- int n1 = AllTris.size();
+ // int n1 = AllTris.size();
std::set<MTri3*,compareTri3Ptr>::iterator itd = AllTris.begin();
while(itd != AllTris.end()){
if((*itd)->isDeleted()){
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index be25549f3e2e3127ef9582f21101968ac67ff565..8a544bf960301ec99ec0cfe8d2cd75b386bc4471 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -671,14 +671,14 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D ||
CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD ||
CTX::instance()->mesh.algo2d == ALGO_2D_BAMG){
- sprintf(opts, "-q1.3pY%c", (Msg::GetVerbosity() < 3) ? 'Q':
+ sprintf(opts, "Ype%c", (Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
// removed -q because mesh sizes at vertices were wrong...
// sprintf(opts, "-q1.5pY%c", (Msg::GetVerbosity() < 3) ? 'Q':
// (Msg::GetVerbosity() > 6) ? 'V': '\0');
}
else {
- sprintf(opts, "-q1.5Ype%c", (Msg::GetVerbosity() < 3) ? 'Q':
+ sprintf(opts, "Ype%c", (Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
}
try{
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index cce8ac8dd02e4e58f976f6f4ee8a2ec95958567c..7ba8d7827f1d09899576ed1de29d07dbdc8f37f9 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -235,20 +235,17 @@ void nonrecurFindCavity(std::list<faceXtet> & shell,
}
-bool insertVertex(MVertex *v,
- MTet4 *t,
- MTet4Factory &myFactory,
- std::set<MTet4*,compareTet4Ptr> &allTets,
- std::vector<double> & vSizes,
- std::vector<double> & vSizesBGM,
- std::set<MTet4*,compareTet4Ptr> *activeTets = 0 )
+bool insertVertexB(std::list<faceXtet> &shell,
+ std::list<MTet4*> &cavity,
+ MVertex *v,
+ MTet4 *t,
+ MTet4Factory &myFactory,
+ std::set<MTet4*,compareTet4Ptr> &allTets,
+ std::vector<double> & vSizes,
+ std::vector<double> & vSizesBGM,
+ std::set<MTet4*,compareTet4Ptr> *activeTets = 0 )
{
- std::list<faceXtet> shell;
- std::list<MTet4*> cavity;
std::list<MTet4*> new_cavity;
-
- recurFindCavity(shell, cavity, v, t);
-
// check that volume is conserved
double newVolume = 0;
double oldVolume = 0;
@@ -341,6 +338,24 @@ bool insertVertex(MVertex *v,
}
}
+bool insertVertex(MVertex *v,
+ MTet4 *t,
+ MTet4Factory &myFactory,
+ std::set<MTet4*,compareTet4Ptr> &allTets,
+ std::vector<double> & vSizes,
+ std::vector<double> & vSizesBGM,
+ std::set<MTet4*,compareTet4Ptr> *activeTets = 0 )
+{
+ std::list<faceXtet> shell;
+ std::list<MTet4*> cavity;
+
+ recurFindCavity(shell, cavity, v, t);
+
+ return insertVertexB(shell,cavity,v,t,myFactory,allTets,vSizes,vSizesBGM,activeTets);
+
+}
+
+
static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes)
{
for (int i = 0; i < 4; i++){
@@ -1019,6 +1034,7 @@ void insertVerticesInRegion (GRegion *gr)
double center[3];
double uvw[3];
MTetrahedron *base = worst->tet();
+
double pa[3] = {base->getVertex(0)->x(),
base->getVertex(0)->y(),
base->getVertex(0)->z()};
@@ -1032,22 +1048,30 @@ void insertVerticesInRegion (GRegion *gr)
base->getVertex(3)->y(),
base->getVertex(3)->z()};
- // double UU,VV,WW;
tetcircumcenter(pa,pb,pc,pd, center,&uvw[0],&uvw[1],&uvw[2] );
- // worst->tet()->xyz2uvw(center, uvw);
- // printf("%12.5E %12.5E %12.5E -- %12.5E %12.5E %12.5E \n",
- // UU,VV,WW,uvw[0],uvw[1],uvw[2]);
- /*
- if (!worst->tet()->isInside(uvw[0], uvw[1], uvw[2]) &&
- worst->getRadius() > 20){
- uvw[0] = uvw[1] = uvw[2] = 0.25;
- center[0] = 0.25*(pa[0]+pb[0]+pc[0]+pd[0]);
- center[1] = 0.25*(pa[1]+pb[1]+pc[1]+pd[1]);
- center[2] = 0.25*(pa[2]+pb[2]+pc[2]+pd[2]);
- }
- */
- if(worst->tet()->isInside(uvw[0], uvw[1], uvw[2])){
+
+ //// A TEST !!!
+ std::list<faceXtet> shell;
+ std::list<MTet4*> cavity;
+ MVertex vv (center[0], center[1], center[2], worst->onWhat());
+ recurFindCavity(shell, cavity, &vv, worst);
+ bool FOUND = false;
+ for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
+ MTetrahedron *toto = (*itc)->tet();
+ // (*itc)->setDeleted(false);
+ toto->xyz2uvw(center,uvw);
+ if (toto->isInside(uvw[0], uvw[1], uvw[2])){
+ worst = (*itc);
+ FOUND = true;
+ break;
+ }
+ }
+ /// END TETS
+
+ // worst->tet()->xyz2uvw(center,uvw);
+
+ if(FOUND){
MVertex *v = new MVertex(center[0], center[1], center[2], worst->onWhat());
v->setIndex(NUM++);
double lc1 =
@@ -1060,7 +1084,7 @@ void insertVerticesInRegion (GRegion *gr)
vSizes.push_back(lc1);
vSizesBGM.push_back(lc);
// compute mesh spacing there
- if(!insertVertex(v, worst, myFactory, allTets, vSizes,vSizesBGM)){
+ if(!insertVertexB(shell,cavity,v, worst, myFactory, allTets, vSizes,vSizesBGM)){
COUNT_MISS_1++;
myFactory.changeTetRadius(allTets.begin(), 0.);
delete v;
@@ -1072,6 +1096,7 @@ void insertVerticesInRegion (GRegion *gr)
// printf("point outside %12.5E %12.5E %12.5E %12.5E %12.5E\n",VV,uvw[0], uvw[1], uvw[2],1-uvw[0]-uvw[1]-uvw[2]);
myFactory.changeTetRadius(allTets.begin(), 0.0);
COUNT_MISS_2++;
+ for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc) (*itc)->setDeleted(false);
}
}
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index c637e4626eded60c3823a311f8be95c378e49123..8c5d44206412d03a5b97a59ef068496ee6733c04 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -1,5 +1,7 @@
#include "GmshConfig.h"
#include "surfaceFiller.h"
+#include "Field.h"
+#include "GModel.h"
#include <queue>
#include <stack>
@@ -31,6 +33,7 @@ struct surfacePointWithExclusionRegion {
SPoint2 _center;
SPoint2 _p[4][NUMDIR];
SPoint2 _q[4];
+ SMetric3 _meshMetric;
/*
+ p3
p4 |
@@ -40,7 +43,8 @@ struct surfacePointWithExclusionRegion {
*/
- surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR]){
+ surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], SMetric3 & meshMetric){
+ _meshMetric = meshMetric;
_v = v;
_center = (p[0][0]+p[1][0]+p[2][0]+p[3][0])*.25;
for (int i=0;i<4;i++)_q[i] = _center + (p[i][0]+p[(i+1)%4][0]-_center*2)*FACTOR;
@@ -137,7 +141,8 @@ bool inExclusionZone (SPoint2 &p,
bool compute4neighbors (GFace *gf, // the surface
MVertex *v_center, // the wertex for which we wnt to generate 4 neighbors
bool goNonLinear, // do we compute the position in the real surface which is nonlinear
- SPoint2 newP[4][NUMDIR]) // look into other directions
+ SPoint2 newP[4][NUMDIR], // look into other directions
+ SMetric3 &metricField) // the mesh metric
{
// we assume that v is on surface gf
@@ -145,8 +150,19 @@ bool compute4neighbors (GFace *gf, // the surface
SPoint2 midpoint;
reparamMeshVertexOnFace(v_center, gf, midpoint);
- double size_1 = backgroundMesh::current()->operator()(midpoint[0],midpoint[1],0.0);
- double size_2 = size_1;
+ double L = backgroundMesh::current()->operator()(midpoint[0],midpoint[1],0.0);
+ metricField = SMetric3(1./(L*L));
+ FieldManager *fields = gf->model()->getFields();
+ if(fields->getBackgroundField() > 0){
+ Field *f = fields->get(fields->getBackgroundField());
+ if (!f->isotropic()){
+ (*f)(v_center->x(),v_center->y(),v_center->z(), metricField,gf);
+ }
+ else {
+ L = (*f)(v_center->x(),v_center->y(),v_center->z(), gf);
+ metricField = SMetric3(1./(L*L));
+ }
+ }
// get the unit normal at that point
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(midpoint[0],midpoint[1]));
@@ -166,6 +182,9 @@ bool compute4neighbors (GFace *gf, // the surface
SVector3 t2 = crossprod(t1,n);
t2.normalize();
+ double size_1 = sqrt(1. / dot(t1,metricField,t1));
+ double size_2 = sqrt(1. / dot(t2,metricField,t2));
+
// compute the first fundamental form i.e. the metric tensor at the point
// M_{ij} = s_i \cdot s_j
double M = dot(s1,s1);
@@ -203,14 +222,14 @@ bool compute4neighbors (GFace *gf, // the surface
// printf("%12.5E %12.5E %g %g %g %g %g %g %g %g %g %g %g\n",l1,l2,size_1,size_2,size_param_1,size_param_2,M,N,E,s1.x(),s1.y(),s2.x(),s2.y());
// this is the rectangle in the parameter plane.
- double r1 = 1.e-8*(double)rand() / RAND_MAX;
- double r2 = 1.e-8*(double)rand() / RAND_MAX;
- double r3 = 1.e-8*(double)rand() / RAND_MAX;
- double r4 = 1.e-8*(double)rand() / RAND_MAX;
- double r5 = 1.e-8*(double)rand() / RAND_MAX;
- double r6 = 1.e-8*(double)rand() / RAND_MAX;
- double r7 = 1.e-8* (double)rand() / RAND_MAX;
- double r8 = 1.e-8*(double)rand() / RAND_MAX;
+ double r1 = 0*1.e-8*(double)rand() / RAND_MAX;
+ double r2 = 0*1.e-8*(double)rand() / RAND_MAX;
+ double r3 = 0*1.e-8*(double)rand() / RAND_MAX;
+ double r4 = 0*1.e-8*(double)rand() / RAND_MAX;
+ double r5 = 0*1.e-8*(double)rand() / RAND_MAX;
+ double r6 = 0*1.e-8*(double)rand() / RAND_MAX;
+ double r7 = 0*1.e-8* (double)rand() / RAND_MAX;
+ double r8 = 0*1.e-8*(double)rand() / RAND_MAX;
double newPoint[4][2] = {{midpoint[0] - covar1[0] * size_param_1 +r1,
midpoint[1] - covar1[1] * size_param_1 +r2},
{midpoint[0] - covar2[0] * size_param_2 +r3,
@@ -248,7 +267,7 @@ bool compute4neighbors (GFace *gf, // the surface
// fills a surface with points in order to build a nice
// quad mesh ------------
-void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed){
+void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics){
#if defined(HAVE_RTREE)
// get all the boundary vertices
@@ -262,21 +281,17 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed){
}
// put boundary vertices in a fifo queue
-
- // let us start without r-tree for now
-
std::queue<surfacePointWithExclusionRegion*> fifo;
std::vector<surfacePointWithExclusionRegion*> vertices;
// put the RTREE
RTree<surfacePointWithExclusionRegion*,double,2,double> rtree;
-
- SVector3 t1;
+ SMetric3 metricField(1.0);
SPoint2 newp[4][NUMDIR];
std::set<MVertex*>::iterator it = bnd_vertices.begin() ;
for (; it != bnd_vertices.end() ; ++it){
- compute4neighbors (gf, *it, false, newp);
+ compute4neighbors (gf, *it, false, newp, metricField);
surfacePointWithExclusionRegion *sp =
- new surfacePointWithExclusionRegion (*it, newp);
+ new surfacePointWithExclusionRegion (*it, newp, metricField);
fifo.push(sp);
vertices.push_back(sp);
double _min[2],_max[2];
@@ -303,9 +318,9 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed){
GPoint gp = gf->point(parent->_p[i][dir]);
MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
// printf(" %g %g %g %g\n",parent._center.x(),parent._center.y(),gp.u(),gp.v());
- compute4neighbors (gf, v, false, newp);
+ compute4neighbors (gf, v, false, newp, metricField);
surfacePointWithExclusionRegion *sp =
- new surfacePointWithExclusionRegion (v, newp);
+ new surfacePointWithExclusionRegion (v, newp, metricField);
fifo.push(sp);
vertices.push_back(sp);
double _min[2],_max[2];
@@ -327,6 +342,7 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed){
for (int i=0;i<vertices.size();i++){
if(vertices[i]->_v->onWhat() == gf) {
packed.push_back(vertices[i]->_v);
+ metrics.push_back(vertices[i]->_meshMetric);
SPoint2 midpoint;
reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint);
// fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",vertices[i]._v->x(),vertices[i]._v->y(),vertices[i]._v->z());
diff --git a/Mesh/surfaceFiller.h b/Mesh/surfaceFiller.h
index c5b215be935fcb336f9e6123ae6ef0a1c4247599..9354629babaa3ff9e9b456e6118d4477b9653330 100644
--- a/Mesh/surfaceFiller.h
+++ b/Mesh/surfaceFiller.h
@@ -4,7 +4,8 @@
// bugs and problems to <gmsh@geuz.org>.
//
+#include "STensor3.h"
#include <vector>
class GFace;
class MVertex;
-void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed);
+void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );