From 5b67c1819353d323474e0cfcaf032de9b00f65e6 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 15 Jun 2011 13:49:04 +0000
Subject: [PATCH] Lloyd's algorithm is accessible online (in norm 6) bug
removed in GFaceCompound.cpp (sorry for that, Tristan, my fault)
JF
---
Common/CommandLine.cpp | 5 +-
Geo/GFaceCompound.cpp | 3 +-
Geo/MVertex.cpp | 1 +
Mesh/Field.cpp | 41 ++++++-
Mesh/Generator.cpp | 16 ++-
Mesh/meshGFace.cpp | 16 +--
Mesh/meshGFace.h | 7 +-
Mesh/meshGFaceDelaunayInsertion.cpp | 1 +
Mesh/meshGFaceLloyd.cpp | 90 ++++++++------
Mesh/meshGRegion.cpp | 14 ++-
Mesh/meshGRegionMMG3D.cpp | 174 +++++++++++++++++++--------
benchmarks/step/capot.geo | 4 +-
contrib/mmg3d/build/sources/analar.c | 2 +-
13 files changed, 266 insertions(+), 108 deletions(-)
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 091f202422..4a35f7fdb9 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -279,7 +279,10 @@ void GetOptions(int argc, char *argv[])
}
else if(!strcmp(argv[i] + 1, "optimize_lloyd")) {
i++;
- CTX::instance()->mesh.optimizeLloyd = 1;
+ if(argv[i])
+ CTX::instance()->mesh.optimizeLloyd = atoi(argv[i++]);
+ else
+ Msg::Fatal("Missing number of lloyd iterations");
}
else if(!strcmp(argv[i] + 1, "nopopup")) {
CTX::instance()->noPopup = 1;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index ae9d869410..07ffc98953 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1469,7 +1469,8 @@ GPoint GFaceCompound::point(double par1, double par2) const
if (lt->gf->geomType() != GEntity::DiscreteSurface){
SPoint2 pParam = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
- return lt->gf->point(pParam);
+ GPoint pp = lt->gf->point(pParam);
+ return GPoint(pp.x(),pp.y(),pp.z(),this,par);
}
if(LINEARMESH){
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index bef8812b90..fd431b0e54 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -398,6 +398,7 @@ bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 ¶m,
else{
double uu, vv;
if(v->onWhat() == gf && v->getParameter(0, uu) && v->getParameter(1, vv)){
+ // printf("%d face %d pos %g %g\n",v->getNum(),gf->tag(),uu,vv);
param = SPoint2(uu, vv);
}
else {
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index b26a705e09..7f2a5304b0 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1724,7 +1724,6 @@ public:
if (pp.first.dim ==1){
GEdge *e = GModel::current()->getEdgeByTag(pp.first.ent);
-
// the tangent size at this point is the size of the
// 1D mesh at this point !
// hack : use curvature
@@ -1737,7 +1736,7 @@ public:
}
- if (dist < hwall_n){
+ if (dist < hwall_t){
L1 = lc_t;
L2 = lc_n;
L3 = lc_n;
@@ -1772,6 +1771,44 @@ public:
t1 = crossprod(t3,t2);
}
}
+ else {
+ GFace *gf = GModel::current()->getFaceByTag(pp.first.ent);
+
+ if (dist < hwall_t){
+ L1 = lc_n;
+ L2 = lc_t;
+ L3 = lc_t;
+ t1 = gf->normal(SPoint2(pp.first.u,pp.first.v));
+ t1.normalize();
+ if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
+ t2 = SVector3(1,0,0);
+ else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
+ t2 = SVector3(0,1,0);
+ else
+ t2 = SVector3(0,0,1);
+ t3 = crossprod(t1,t2);
+ t3.normalize();
+ t2 = crossprod(t3,t1);
+ // printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
+ }
+ else {
+ L1 = lc_t;
+ L2 = lc_n;
+ L3 = lc_t;
+ GPoint p = gf->point(SPoint2(pp.first.u,pp.first.v));
+ t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
+ if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
+ t1 = SVector3(1,0,0);
+ else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
+ t1 = SVector3(0,1,0);
+ else
+ t1 = SVector3(0,0,1);
+ t2.normalize();
+ t3 = crossprod(t1,t2);
+ t3.normalize();
+ t1 = crossprod(t3,t2);
+ }
+ }
}
else{
double delta = std::min(CTX::instance()->lc / 1e2, dist);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index abbc9828ee..5ff07f954a 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -21,6 +21,7 @@
#include "meshGEdge.h"
#include "meshGFace.h"
#include "meshGFaceOptimize.h"
+#include "meshGFaceLloyd.h"
#include "meshGFaceBDS.h"
#include "meshGRegion.h"
#include "BackgroundMesh.h"
@@ -465,11 +466,21 @@ static void Mesh2D(GModel *m)
}
if(!nbPending) break;
if(nIter++ > 10) break;
- }
+ }
}
+
// lloyd optimization
if (CTX::instance()->mesh.optimizeLloyd){
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ if((*it)->geomType()==GEntity::CompoundSurface){
+ smoothing smm(CTX::instance()->mesh.optimizeLloyd,6);
+ smm.optimize_face(*it);
+ int rec = (CTX::instance()->mesh.recombineAll || (*it)->meshAttributes.recombine);
+ if(rec) recombineIntoQuads(*it);
+ }
+ }
+ /*
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
GFace *gf = *it;
if(gf->geomType() == GEntity::DiscreteSurface) continue;
@@ -477,11 +488,10 @@ static void Mesh2D(GModel *m)
GFaceCompound *gfc = (GFaceCompound*) gf;
if(gfc->getNbSplit() != 0) continue;
}
- int rec = (CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine);
Msg::Info("Lloyd optimization for face %d", gf->tag());
gf->lloyd(25, rec);
- if(rec) recombineIntoQuads(gf);
}
+ */
}
// collapseSmallEdges(*m);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index d9a157434e..077de5a831 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -301,7 +301,8 @@ static bool algoDelaunay2D(GFace *gf)
if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY ||
CTX::instance()->mesh.algo2d == ALGO_2D_BAMG ||
CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL ||
- CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD)
+ CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD ||
+ CTX::instance()->mesh.algo2d == ALGO_2D_BAMG)
return true;
if(CTX::instance()->mesh.algo2d == ALGO_2D_AUTO && gf->geomType() == GEntity::Plane)
@@ -393,6 +394,7 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
// the domain, including embedded points and surfaces
static bool meshGenerator(GFace *gf, int RECUR_ITER,
bool repairSelfIntersecting1dMesh,
+ bool onlyInitialMesh,
bool debug = true)
{
@@ -406,7 +408,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
// replace edges by their compounds
// if necessary split compound and remesh parts
bool isMeshed = false;
- if(gf->geomType() == GEntity::CompoundSurface){
+ if(gf->geomType() == GEntity::CompoundSurface && !onlyInitialMesh){
isMeshed = checkMeshCompound((GFaceCompound*) gf, edges);
if (isMeshed) return true;
}
@@ -631,7 +633,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
delete m;
if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
return meshGenerator
- (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, debug);
+ (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh, debug);
return false;
}
@@ -825,7 +827,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
Msg::Debug("Starting to add internal points");
// start mesh generation
- if(!algoDelaunay2D(gf)){
+ if(!algoDelaunay2D(gf) && !onlyInitialMesh){
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true,
&recoverMapInv);
optimizeMeshBDS(gf, *m, 2);
@@ -883,7 +885,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
// the delaunay algo is based directly on internal gmsh structures
// BDS mesh is passed in order not to recompute local coordinates of
// vertices
- if(algoDelaunay2D(gf)){
+ if(algoDelaunay2D(gf) && !onlyInitialMesh){
if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
bowyerWatsonFrontal(gf);
else if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD)
@@ -1511,7 +1513,7 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
-int debugSurface = -1;
+int debugSurface = -100;
void meshGFace::operator() (GFace *gf)
{
@@ -1571,7 +1573,7 @@ void meshGFace::operator() (GFace *gf)
(!gf->periodic(0) && !gf->periodic(1))) &&
(noSeam(gf) || gf->getNativeType() == GEntity::GmshModel ||
gf->edgeLoops.empty())){
- meshGenerator(gf, 0, repairSelfIntersecting1dMesh,
+ meshGenerator(gf, 0, repairSelfIntersecting1dMesh, onlyInitialMesh,
debugSurface >= 0 || debugSurface == -100);
}
else {
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index 3f79ce359b..4619dbe766 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -18,9 +18,14 @@ class GFaceCompound;
class meshGFace {
const bool repairSelfIntersecting1dMesh;
int twoPassesMesh;
+ bool onlyInitialMesh;
public :
- meshGFace (bool r = true, int t = 0) : repairSelfIntersecting1dMesh(r), twoPassesMesh(t){}
+ meshGFace (bool r = true, int t = 0) : repairSelfIntersecting1dMesh(r), twoPassesMesh(t)
+ ,onlyInitialMesh(false)
+ {
+ }
void operator () (GFace *);
+ void setOnlyInitial() {onlyInitialMesh = true;}
};
// Destroy the mesh of the face
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 4824297c27..b6da957ef9 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1039,6 +1039,7 @@ void bowyerWatsonFrontal(GFace *gf)
ActiveTris.insert(*it);
else if ((*it)->getRadius() < LIMIT_)break;
}
+
// insert points
while (1){
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 8b7de441c2..18bbb11778 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -178,8 +178,9 @@ void smoothing::optimize_face(GFace* gf){
const double LC2D = sqrt((du.high()-du.low())*(du.high()-du.low()) +
(dv.high()-dv.low())*(dv.high()-dv.low()));
- //printf("Lloyd on face %d %d elements %d nodes LC %g\n", gf->tag(),
- // gf->getNumMeshElements(), (int)all.size(), LC2D);
+
+ printf("Lloyd on face %d %d elements %d nodes LC %g\n", gf->tag(),
+ gf->getNumMeshElements(), (int)all.size(), LC2D);
int i = 0;
for (std::set<MVertex*>::iterator it = all.begin(); it != all.end(); ++it){
@@ -197,8 +198,11 @@ void smoothing::optimize_face(GFace* gf){
triangulator.x(i) = p.x() + XX;
triangulator.y(i) = p.y() + YY;
triangulator.data(i++) = (*it);
+
+ // GPoint pp = gf->point(p);
+ // printf("point %d %g %g %g %g %g \n",(*it)->getNum(),(*it)->x(),(*it)->y(),(*it)->z(),p.x(),p.y());
}
-
+
triangulator.Voronoi();
triangulator.initialize();
int index,number,count,max;
@@ -209,13 +213,13 @@ void smoothing::optimize_face(GFace* gf){
for(int i=0;i<max;i++)
{
if(count>=number) break;
- index = (int)((triangulator.numPoints-1)*((double)rand()/(double)RAND_MAX));
- PointRecord& pt = triangulator.points[index];
- MVertex* v = (MVertex*)pt.data;
- if(v->onWhat()==gf && !triangulator.onHull(index)){
- flag = triangulator.remove_point(index);
- if(flag) count++;
- }
+ index = (int)((triangulator.numPoints-1)*((double)rand()/(double)RAND_MAX));
+ PointRecord& pt = triangulator.points[index];
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(index)){
+ flag = triangulator.remove_point(index);
+ if(flag) count++;
+ }
}
triangulator.remove_all();
@@ -224,11 +228,11 @@ void smoothing::optimize_face(GFace* gf){
delta = 0.01;
for(int i=0;i<triangulator.numPoints;i++){
PointRecord& pt = triangulator.points[i];
- MVertex* v = (MVertex*)pt.data;
- if(v->onWhat()==gf && !triangulator.onHull(i)){
- //triangulator.points[i].where.h = delta + (1.0-2.0*delta)*((double)rand()/(double)RAND_MAX);
- //triangulator.points[i].where.v = delta + (1.0-2.0*delta)*((double)rand()/(double)RAND_MAX);
- }
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ //triangulator.points[i].where.h = delta + (1.0-2.0*delta)*((double)rand()/(double)RAND_MAX);
+ //triangulator.points[i].where.v = delta + (1.0-2.0*delta)*((double)rand()/(double)RAND_MAX);
+ }
}
int index1 = -1;
@@ -241,18 +245,18 @@ void smoothing::optimize_face(GFace* gf){
int index8 = -1;
for(int i=0;i<triangulator.numPoints;i++){
PointRecord& pt = triangulator.points[i];
- MVertex* v = (MVertex*)pt.data;
- if(v->onWhat()==gf && !triangulator.onHull(i)){
- if(index1==-1) index1 = i;
- else if(index2==-1) index2 = i;
- else if(index3==-1) index3 = i;
- else if(index4==-1) index4 = i;
- else if(index5==-1) index5 = i;
- else if(index6==-1) index6 = i;
- else if(index7==-1) index7 = i;
- else if(index8==-1) index8 = i;
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ if(index1==-1) index1 = i;
+ else if(index2==-1) index2 = i;
+ else if(index3==-1) index3 = i;
+ else if(index4==-1) index4 = i;
+ else if(index5==-1) index5 = i;
+ else if(index6==-1) index6 = i;
+ else if(index7==-1) index7 = i;
+ else if(index8==-1) index8 = i;
- }
+ }
}
/*triangulator.points[index1].where.h = 0.01;
triangulator.points[index1].where.v = 0.01;
@@ -300,17 +304,17 @@ void smoothing::optimize_face(GFace* gf){
num_interior = 0;
for(int i=0;i<triangulator.numPoints;i++){
if(obj.interior(triangulator,gf,i)){
- num_interior++;
- }
+ num_interior++;
+ }
}
index = 0;
for(int i=0;i<triangulator.numPoints;i++){
- if(obj.interior(triangulator,gf,i)){
- initial_conditions[index] = triangulator.points[i].where.h;
- initial_conditions[index+num_interior] = triangulator.points[i].where.v;
- index++;
- }
+ if(obj.interior(triangulator,gf,i)){
+ initial_conditions[index] = triangulator.points[i].where.h;
+ initial_conditions[index+num_interior] = triangulator.points[i].where.v;
+ index++;
+ }
}
x.setcontent(2*num_interior, &initial_conditions[0]);
@@ -342,11 +346,11 @@ void smoothing::optimize_face(GFace* gf){
index = 0;
for(i=0;i<triangulator.numPoints;i++){
- if(obj.interior(triangulator,gf,i)){
- triangulator.points[i].where.h = x[index];
- triangulator.points[i].where.v = x[index + num_interior];
- index++;
- }
+ if(obj.interior(triangulator,gf,i)){
+ triangulator.points[i].where.h = x[index];
+ triangulator.points[i].where.v = x[index + num_interior];
+ index++;
+ }
}
triangulator.Voronoi();
@@ -368,9 +372,13 @@ void smoothing::optimize_face(GFace* gf){
// get the ith vertex
PointRecord &pt = triangulator.points[i];
MVertex *v = (MVertex*)pt.data;
+ if (v->onWhat() == gf && triangulator.onHull(i)){
+ // printf("strange !!\n");
+ }
if (v->onWhat() == gf && !triangulator.onHull(i)){
GPoint gp = gf->point (pt.where.h,pt.where.v);
- MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
+ // printf("%g %g vs %g %g\n",pt.where.h,pt.where.v,gp.u(),gp.v());
+ MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,pt.where.h,pt.where.v);
mesh_vertices.push_back(v);
}
}
@@ -378,13 +386,17 @@ void smoothing::optimize_face(GFace* gf){
// destroy the mesh
deMeshGFace killer;
killer(gf);
+ backgroundMesh::unset();
// put all additional vertices in the list of
// vertices
gf->additionalVertices = mesh_vertices;
// remesh the face with all those vertices in
Msg::Info("Lloyd remeshing of face %d ", gf->tag());
+
+
meshGFace mesher;
+ mesher.setOnlyInitial();
mesher(gf);
// assign those vertices to the face internal vertices
gf->mesh_vertices.insert(gf->mesh_vertices.begin(),
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 4fb32b0633..0452550af2 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -503,8 +503,18 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
char opts[128];
buildTetgenStructure(gr, in, numberedV);
//if (Msg::GetVerbosity() == 20) sprintf(opts, "peVvS0");
- sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q':
- (Msg::GetVerbosity() > 6) ? 'V': '\0');
+ if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL ||
+ CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_HEX ||
+ CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D ||
+ CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD ||
+ CTX::instance()->mesh.algo2d == ALGO_2D_BAMG){
+ sprintf(opts, "pY", (Msg::GetVerbosity() < 3) ? 'Q':
+ (Msg::GetVerbosity() > 6) ? 'V': '\0');
+ }
+ else {
+ sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q':
+ (Msg::GetVerbosity() > 6) ? 'V': '\0');
+ }
try{
tetrahedralize(opts, &in, &out);
}
diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp
index 9df01152a8..938fc8f38d 100644
--- a/Mesh/meshGRegionMMG3D.cpp
+++ b/Mesh/meshGRegionMMG3D.cpp
@@ -16,25 +16,27 @@ extern "C" {
}
void MMG2gmsh (GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gmsh){
- printf("%d vertices in MMG\n",mmg->np);
+ std::map<int,MVertex*> kToMVertex;
for (int k=1;k<= mmg->np ; k++){
MMG_pPoint ppt = &mmg->point[k];
if (ppt->tag & M_UNUSED) continue;
- if (mmg2gmsh.find(k) == mmg2gmsh.end()){
+ std::map<int,MVertex*>::iterator it = mmg2gmsh.find(ppt->ref);
+ if (it == mmg2gmsh.end()){
MVertex *v = new MVertex(ppt->c[0],ppt->c[1],ppt->c[2],gr);
- mmg2gmsh[k] = v;
gr->mesh_vertices.push_back(v);
+ kToMVertex[k] = v;
}
+ else kToMVertex[k] = it->second;
}
for (int k=1; k<=mmg->ne; k++) {
MMG_pTetra ptetra = &mmg->tetra[k];
if ( ptetra->v[0] ){
- MVertex *v1 = mmg2gmsh[ptetra->v[0]];
- MVertex *v2 = mmg2gmsh[ptetra->v[1]];
- MVertex *v3 = mmg2gmsh[ptetra->v[2]];
- MVertex *v4 = mmg2gmsh[ptetra->v[3]];
+ MVertex *v1 = kToMVertex[ptetra->v[0]];
+ MVertex *v2 = kToMVertex[ptetra->v[1]];
+ MVertex *v3 = kToMVertex[ptetra->v[2]];
+ MVertex *v4 = kToMVertex[ptetra->v[3]];
if (!v1 || !v2 || !v3 || !v4){
Msg::Error("Element %d Unknown Vertex in MMG2gmsh %d(%p) %d(%p) %d(%p) %d(%p)",k,ptetra->v[0],v1,ptetra->v[1],v2,ptetra->v[2],v3,ptetra->v[3],v4);
}
@@ -44,6 +46,7 @@ void MMG2gmsh (GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gmsh){
}
void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*> &mmg2gmsh){
+
mmg->ne = gr->tetrahedra.size();
std::set<MVertex*> allVertices;
for (int i=0;i< gr->tetrahedra.size() ; i++){
@@ -76,26 +79,25 @@ void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*>
sol->metold = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
std::map<MVertex*,std::pair<double,int> > LCS;
- if (CTX::instance()->mesh.lcExtendFromBoundary){
- for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
- for (int i=0;i<(*it)->triangles.size();i++){
- MTriangle *t = (*it)->triangles[i];
- double L = t->maxEdge();
- for (int k=0;k<3;k++){
- MVertex *v = t->getVertex(k);
- std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
- if (itv != LCS.end()){
- itv->second.first += L;
- itv->second.second ++;
- }
- else {
- LCS[v] = std::make_pair(L,1);
- }
+ for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
+ for (int i=0;i<(*it)->triangles.size();i++){
+ MTriangle *t = (*it)->triangles[i];
+ double L = t->maxEdge();
+ for (int k=0;k<3;k++){
+ MVertex *v = t->getVertex(k);
+ std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
+ if (itv != LCS.end()){
+ itv->second.first += L;
+ itv->second.second ++;
+ }
+ else {
+ LCS[v] = std::make_pair(L,1);
}
}
}
}
+ printf("%d vertices %d on faces\n",allVertices.size(), LCS.size());
int k=1;
int count = 1;//sol->offset;
@@ -103,12 +105,10 @@ void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*>
for (std::set<MVertex*>::iterator it = allVertices.begin() ; it != allVertices.end() ; ++it){
MMG_pPoint ppt = &mmg->point[k];
- mmg2gmsh[k] = *it;
-
ppt->c[0] = (*it)->x();
ppt->c[1] = (*it)->y();
ppt->c[2] = (*it)->z();
- ppt->ref = gr->tag();
+ ppt->ref = (*it)->getNum();
gmsh2mmg_num[(*it)->getNum()] = k;
MVertex *v = *it;
@@ -120,26 +120,27 @@ void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*>
v->getParameter(0,U);
v->getParameter(1,V);
}
- double lc = BGM_MeshSize(v->onWhat(), U,V,v->x(), v->y(), v->z());
- //SMetric3 m = BGM_MeshMetric(v->onWhat(), U,V,v->x(), v->y(), v->z());
+ //double lc = BGM_MeshSize(v->onWhat(), U,V,v->x(), v->y(), v->z());
+ SMetric3 m = BGM_MeshMetric(v->onWhat(), U,V,v->x(), v->y(), v->z());
- if (CTX::instance()->mesh.lcExtendFromBoundary){
- std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
- if (itv != LCS.end()){
+ std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
+ if (itv != LCS.end()){
+ mmg2gmsh[(*it)->getNum()] = *it;
+ if (CTX::instance()->mesh.lcExtendFromBoundary){
double LL = itv->second.first/itv->second.second;
- // SMetric3 l4(1./(LL*LL));
- // SMetric3 MM = intersection (l4, m);
- // m = MM;
- lc = std::min(LL,lc);
+ SMetric3 l4(1./(LL*LL));
+ SMetric3 MM = intersection (l4, m);
+ m = MM;
+ //lc = std::min(LL,lc);
}
}
- sol->met[count++] = 1/(lc*lc);
- sol->met[count++] = 0;
- sol->met[count++] = 0;
- sol->met[count++] = 1/(lc*lc);
- sol->met[count++] = 0;
- sol->met[count++] = 1/(lc*lc);
+ sol->met[count++] = m(0,0);
+ sol->met[count++] = m(1,0);
+ sol->met[count++] = m(2,0);
+ sol->met[count++] = m(1,1);
+ sol->met[count++] = m(2,1);
+ sol->met[count++] = m(2,2);
// printf("%g %g %g %g %g %g\n",m(0,0),m(0,1),m(0,2),m(1,1),m(1,2),m(2,2));
// for (int i=0; i<sol->offset; i++) {
@@ -173,23 +174,98 @@ void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*>
}
+void updateSizes (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol){
+ std::list<GFace*> f = gr->faces();
+ /*
+ std::map<MVertex*,std::pair<double,int> > LCS;
+ if (CTX::instance()->mesh.lcExtendFromBoundary){
+ for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
+ for (int i=0;i<(*it)->triangles.size();i++){
+ MTriangle *t = (*it)->triangles[i];
+ double L = t->maxEdge();
+ for (int k=0;k<3;k++){
+ MVertex *v = t->getVertex(k);
+ std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
+ if (itv != LCS.end()){
+ itv->second.first += L;
+ itv->second.second ++;
+ }
+ else {
+ LCS[v] = std::make_pair(L,1);
+ }
+ }
+ }
+ }
+ }
+ */
+
+ int count = 1;
+ for (int k=1 ; k<=mmg->np; k++){
+ MMG_pPoint ppt = &mmg->point[k];
+ if (ppt->tag & M_UNUSED) continue;
+
+ SMetric3 m = BGM_MeshMetric(gr, 0,0,ppt->c[0],ppt->c[1],ppt->c[2]);
+ /*
+ if (CTX::instance()->mesh.lcExtendFromBoundary){
+ std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
+ if (itv != LCS.end()){
+ double LL = itv->second.first/itv->second.second;
+ SMetric3 l4(1./(LL*LL));
+ SMetric3 MM = intersection (l4, m);
+ m = MM;
+ }
+ }
+ */
+ sol->met[count++] = m(0,0);
+ sol->met[count++] = m(1,0);
+ sol->met[count++] = m(2,0);
+ sol->met[count++] = m(1,1);
+ sol->met[count++] = m(2,1);
+ sol->met[count++] = m(2,2);
+ }
+ free(sol->metold);
+ sol->metold = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
+}
+
+
+void FREEMMG (MMG_pMesh mmgMesh, MMG_pSol mmgSol){
+ free(mmgMesh->point);
+ // free(mmgMesh->disp->alpha);
+ // free(mmgMesh->disp->mv);
+ free(mmgMesh->disp);
+ free(mmgMesh->tria);
+ free(mmgMesh->tetra);
+ free(mmgMesh);
+ if ( mmgSol->npfixe ){
+ free(mmgSol->met);
+ free(mmgSol->metold);
+ }
+ free(mmgSol);
+}
+
void refineMeshMMG(GRegion *gr){
MMG_pMesh mmg = (MMG_pMesh)calloc(1,sizeof(MMG_Mesh));
MMG_pSol sol = (MMG_pSol)calloc(1,sizeof(MMG_Sol));
std::map<int,MVertex*> mmg2gmsh;
gmsh2MMG (gr, mmg, sol,mmg2gmsh);
- int opt[9] = {1,0,64,0,0,(Msg::GetVerbosity() > 20) ? 3 : -1,0,0,0};
- mmg3d::MMG_mmg3dlib(opt,mmg,sol);
+
+ for (int ITER=0;ITER<2;ITER++){
+ int opt[9] = {1,0,64,0,0,(Msg::GetVerbosity() > 3) ? 335 : 3333,0,0,0};
+ Msg::Debug("-------- GMSH LAUNCHES MMG3D ---------------");
+ mmg3d::MMG_mmg3dlib(opt,mmg,sol);
+ Msg::Debug("-------- MG3D TERMINATED -------------------");
+ updateSizes (gr,mmg, sol);
+ }
+ MMG_saveMesh(mmg,"test.mesh");
+ gr->deleteVertexArrays();
for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i];
gr->tetrahedra.clear();
- for (int i=0;i<gr->mesh_vertices.size();++i)delete gr->mesh_vertices[i];
- gr->mesh_vertices.clear();
-
-
- MMG2gmsh (gr, mmg, mmg2gmsh);
- char test[] = "test.mesh";
- MMG_saveMesh(mmg, test);
+ // for (int i=0;i<gr->mesh_vertices.size();++i)delete gr->mesh_vertices[i];
+ //gr->mesh_vertices.clear();
+
+ MMG2gmsh (gr, mmg, mmg2gmsh);
+ FREEMMG (mmg,sol);
}
#else
diff --git a/benchmarks/step/capot.geo b/benchmarks/step/capot.geo
index d1e2eda5c1..ba6e70c436 100644
--- a/benchmarks/step/capot.geo
+++ b/benchmarks/step/capot.geo
@@ -1,9 +1,9 @@
Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
Mesh.RemeshAlgorithm=0; //(0) nosplit (1) automatic (2) split metis
-Mesh.CharacteristicLengthFactor=0.2;
+//Mesh.CharacteristicLengthFactor=0.2;
//Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
-//Mesh.RecombinationAlgorithm = 1;
+Mesh.RecombinationAlgorithm = 1;
Merge "capot.brep";
diff --git a/contrib/mmg3d/build/sources/analar.c b/contrib/mmg3d/build/sources/analar.c
index ee4aa1e5bb..4f56e8e3f6 100644
--- a/contrib/mmg3d/build/sources/analar.c
+++ b/contrib/mmg3d/build/sources/analar.c
@@ -171,7 +171,7 @@ MMG_bouffe = 0;
for (k=1; k<=nedep; k++) {
pt = &mesh->tetra[k];
if ( !pt->v[0] ) continue;
- else if ( pt->flag != base-1 ) continue;
+ // else if ( pt->flag != base-1 ) continue;
if ( pt->qual < declic ) continue;
pt->flag = base-2;
--
GitLab