From 0628eba005681e4e6e89406f5f085dffe83038c9 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 10 Oct 2007 08:49:34 +0000
Subject: [PATCH] *** empty log message ***
---
Geo/MVertex.h | 7 +-
Mesh/BackgroundMesh.cpp | 20 +-
Mesh/BackgroundMesh.h | 3 +-
Mesh/meshGEdge.cpp | 8 +-
Mesh/meshGFace.cpp | 288 ++++++++++++++++------------
Mesh/meshGFaceDelaunayInsertion.cpp | 4 +-
benchmarks/2d/Square-Attr2.geo | 2 +-
7 files changed, 198 insertions(+), 134 deletions(-)
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 67311529f7..9406a296fb 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -117,15 +117,16 @@ class MVertex{
class MEdgeVertex : public MVertex{
protected:
- double _u;
+ double _u, _lc;
public :
- MEdgeVertex(double x, double y, double z, GEntity *ge, double u)
- : MVertex(x, y, z, ge), _u(u)
+ MEdgeVertex(double x, double y, double z, GEntity *ge, double u, double lc = -1.0)
+ : MVertex(x, y, z, ge), _u(u),_lc(lc)
{
}
virtual ~MEdgeVertex(){}
virtual bool getParameter(int i, double &par) const{ par = _u; return true; }
virtual bool setParameter(int i, double par){ _u = par; return true; }
+ double getLc () const {return _lc;}
};
class MFaceVertex : public MVertex{
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 826ccb2b16..321199353e 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.25 2007-09-10 04:47:03 geuzaine Exp $
+// $Id: BackgroundMesh.cpp,v 1.26 2007-10-10 08:49:34 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -183,7 +183,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
double l4 = lc_field.empty() ? MAX_LC : lc_field(X, Y, Z);
// use the field unconstrained by other characteristic lengths
- if(l4 < MAX_LC && !CTX.mesh.constrained_bgmesh)
+ if(!lc_field.empty() && !CTX.mesh.constrained_bgmesh)
return l4 * CTX.mesh.lc_factor;
if(CTX.mesh.lc_from_curvature && ge->dim() < 3)
@@ -195,5 +195,21 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
// printf("l1 = %12.5E l2 = %12.5E l4 = %12.5E\n",l1,l2,l4);
double lc = std::min(std::min(std::min(l1, l2), l3), l4);
+
return lc * CTX.mesh.lc_factor;
}
+
+// we extend the 1d mesh in surfaces if no background mesh exists
+// in this case, it is the only way to have something smooth
+// we do it also if CTX.mesh.constrained_bgmesh is true;
+bool Extend1dMeshIn2dSurfaces ()
+{
+ if ( lc_field.empty()) return true;
+ if ( CTX.mesh.constrained_bgmesh == true) return true;
+ return false;
+}
+
+bool Extend2dMeshIn3dVolumes ()
+{
+ return Extend1dMeshIn2dSurfaces ();
+}
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 2d761aceb8..ec7958c8cb 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -26,5 +26,6 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
bool BGMExists();
void BGMAddField(Field *field);
void BGMReset();
-
+bool Extend1dMeshIn2dSurfaces ();
+bool Extend2dMeshIn3dVolumes ();
#endif
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index b7887242c0..cccfc7e1f4 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.43 2007-09-24 08:14:29 geuzaine Exp $
+// $Id: meshGEdge.cpp,v 1.44 2007-10-10 08:49:34 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -365,10 +365,12 @@ void meshGEdge::operator() (GEdge *ge)
const double d = (double)NUMP * b;
if((fabs(P2.p) >= fabs(d)) && (fabs(P1.p) < fabs(d))) {
double dt = P2.t - P1.t;
+ double dlc = P2.lc - P1.lc;
double dp = P2.p - P1.p;
- double t = P1.t + dt / dp * (d - P1.p);
+ double t = P1.t + dt / dp * (d - P1.p);
+ double lc = P1.lc + dlc / dp * (d - P1.p);
GPoint V = ge->point(t);
- ge->mesh_vertices[NUMP - 1] = new MEdgeVertex(V.x(), V.y(), V.z(), ge, t);
+ ge->mesh_vertices[NUMP - 1] = new MEdgeVertex(V.x(), V.y(), V.z(), ge, t, lc);
NUMP++;
}
else {
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 7dde59d0fb..1b0f5a6c16 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.89 2007-09-24 08:14:29 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.90 2007-10-10 08:49:34 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -106,7 +106,7 @@ double F_LC_ANALY(double xx, double yy, double zz)
double NewGetLc(BDS_Point *p)
{
- return std::min(p->lc(),p->lcBGM());
+ return Extend1dMeshIn2dSurfaces () ? std::min(p->lc(),p->lcBGM()) : p->lcBGM();
}
inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
@@ -297,6 +297,102 @@ void OptimizeMesh(GFace *gf, BDS_Mesh &m, const int NIT)
}
}
+void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
+{
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1)
+ {
+ if (NN2++ >= NN1)break;
+ // result = -1 => forbid swap because too badly shaped elements
+ // result = 0 => whatever
+ // 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++;
+ ++it;
+ }
+ }
+}
+
+void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
+{
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1)
+ {
+ if (NN2++ >= NN1)break;
+ if (!(*it)->deleted)
+ {
+ double lone = NewGetLc ( *it,gf);
+ if ((*it)->numfaces() == 2 && (lone > MAXE_))
+ {
+ const double coord = 0.5;
+ BDS_Point *mid ;
+ mid = m.add_point(++m.MAXPOINTNUMBER,
+ coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
+ coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v,gf);
+ double l1;
+ // if (BGMExists())
+ mid->lcBGM() = BGM_MeshSize(gf,
+ (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
+ (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
+ 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++;
+ }
+ }
+ ++it;
+ }
+}
+
+
+void collapseEdgePass ( GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
+{
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1)
+ {
+ if (NN2++ >= NN1)break;
+ double lone = NewGetLc ( *it,gf);
+ // if (lone < 1.e-10 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-5) lone = 2;
+
+ if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < MINE_ )
+ {
+ bool res = false;
+ if ( (*it)->p1->iD > MAXNP )
+ res =m.collapse_edge_parametric ( *it, (*it)->p1);
+ else if ( (*it)->p2->iD > MAXNP )
+ res =m.collapse_edge_parametric ( *it, (*it)->p2);
+ if (res)
+ nb_collaps ++;
+ }
+ ++it;
+ }
+}
+
+
+void smoothVertexPass ( GFace *gf, BDS_Mesh &m, int &nb_smooth)
+{
+ std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
+ while (itp != m.points.end())
+ {
+
+ if(m.smooth_point_centroid(*itp,gf))
+ nb_smooth ++;
+ ++itp;
+ }
+}
+
void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
{
int IT =0;
@@ -305,9 +401,9 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
int MAXNP = m.MAXPOINTNUMBER;
- // computecharacteristic lengths using mesh edge spacing
+ // computecharacteristic lengths using 1D mesh edge spacing
// those lengths will propagate in the 2D mesh.
- if (NIT > 0)
+ if (0 && NIT > 0)
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
@@ -320,8 +416,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
if (l<L && (*it)->g && (*it)->g->classif_degree == 1)L=l;
++it;
}
- if(!CTX.mesh.constrained_bgmesh)
- (*itp)->lc() = L;
+ (*itp)->lc() = L;
(*itp)->lcBGM() = L;
++itp;
}
@@ -329,6 +424,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
double OLDminL=1.E22,OLDmaxL=0;
+ double t_spl=0, t_sw=0,t_col=0,t_sm=0;
+
const double MINE_ = 0.7, MAXE_=1.4;
while (1)
{
@@ -342,6 +439,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
// split long edges
double minL=1.E22,maxL=0;
+
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
@@ -364,116 +462,45 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT)))break;
- NN1 = m.edges.size();
- NN2 = 0;
- it = m.edges.begin();
- while (1)
- {
- if (NN2++ >= NN1)break;
- if (!(*it)->deleted)
- {
- double lone = NewGetLc ( *it,gf);
- if ((*it)->numfaces() == 2 && (lone > MAXE_))
- {
- const double coord = 0.5;
- BDS_Point *mid ;
- mid = m.add_point(++m.MAXPOINTNUMBER,
- coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
- coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v,gf);
- // if (BGMExists())
- mid->lcBGM() = BGM_MeshSize(gf,
- (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
- (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
- 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++;
- }
- }
- ++it;
- }
-
- // swap edges that provide a better configuration
-
- NN2 = 0;
- it = m.edges.begin();
- while (1)
- {
- if (NN2++ >= NN1)break;
- // result = -1 => forbid swap because too badly shaped elements
- // result = 0 => whatever
- // 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++;
- ++it;
- }
- }
-
- // collapse short edges
- NN1 = m.edges.size();
- NN2 = 0;
- it = m.edges.begin();
- while (1)
- {
- if (NN2++ >= NN1)break;
- double lone = NewGetLc ( *it,gf);
- // if (lone < 1.e-10 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-5) lone = 2;
+ double maxE = MAXE_;//std::max(MAXE_,maxL / 5);
+ double minE = MINE_;//std::min(MINE_,minL * 1.2);
+ clock_t t1 = clock();
+ splitEdgePass ( gf, m, maxE, nb_split);
+ clock_t t2 = clock();
+ swapEdgePass ( gf, m, nb_swap);
+ clock_t t3 = clock();
+ collapseEdgePass ( gf, m, minE , MAXNP, nb_collaps);
+ clock_t t4 = clock();
+ swapEdgePass ( gf, m, nb_swap);
+ clock_t t5 = clock();
+ smoothVertexPass ( gf, m, nb_smooth);
+ clock_t t6 = clock();
+ // clean up the mesh
- if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < MINE_ )
- {
- bool res = false;
- if ( (*it)->p1->iD > MAXNP )
- res =m.collapse_edge_parametric ( *it, (*it)->p1);
- else if ( (*it)->p2->iD > MAXNP )
- res =m.collapse_edge_parametric ( *it, (*it)->p2);
- if (res)
- nb_collaps ++;
- }
- ++it;
- }
- // swap edges that provide a better configuration
- NN1 = m.edges.size();
- NN2 = 0;
- it = m.edges.begin();
- while (1)
- {
- if (NN2++ >= NN1)break;
- 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++;
- ++it;
- }
- }
- // smooth resulting mesh
- {
- std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
- while (itp != m.points.end())
- {
-
- if(m.smooth_point_centroid(*itp,gf))
- nb_smooth ++;
- ++itp;
- }
- }
+ t_spl += (double)(t2-t1)/CLOCKS_PER_SEC;
+ t_sw += (double)(t3-t2)/CLOCKS_PER_SEC;
+ t_sw += (double)(t5-t4)/CLOCKS_PER_SEC;
+ t_col += (double)(t4-t3)/CLOCKS_PER_SEC;
+ t_sm += (double)(t6-t5)/CLOCKS_PER_SEC;
- // clean up the mesh
m.cleanup();
IT++;
- Msg(DEBUG1," iter %3d minL %8.3f maxL %8.3f : %6d splits, %6d swaps, %6d collapses, %6d moves",IT,minL,maxL,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",IT,minL,minE,maxL,maxE,nb_split,nb_swap,nb_collaps,nb_smooth);
if (nb_split==0 && nb_collaps == 0)break;
}
+ double t_total = t_spl + t_sw + t_col + t_sm;
+ Msg(DEBUG1," ---------------------------------------");
+ Msg(DEBUG1," CPU Report ");
+ Msg(DEBUG1," ---------------------------------------");
+ Msg(DEBUG1," CPU SWAP %12.5E sec (%3.0f %%)",t_sw,100 * t_sw/t_total);
+ Msg(DEBUG1," CPU SPLIT %12.5E sec (%3.0f %%) ",t_spl,100 * t_spl/t_total);
+ Msg(DEBUG1," CPU COLLAPS %12.5E sec (%3.0f %%) ",t_col,100 * t_col/t_total);
+ Msg(DEBUG1," CPU SMOOTH %12.5E sec (%3.0f %%) ",t_sm,100 * t_sm/t_total);
+ Msg(DEBUG1," ---------------------------------------");
+ Msg(DEBUG1," CPU TOTAL %12.5E sec ",t_total);
+ Msg(DEBUG1," ---------------------------------------");
}
@@ -760,7 +787,21 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
V = V_[num];
}
- m->add_point ( num, U,V, gf);
+ BDS_Point *pp = m->add_point ( num, U,V, gf);
+
+ GEntity *ge = here->onWhat();
+ if(ge->dim() == 1)
+ {
+ double t;
+ here->getParameter(0,t);
+ pp->lcBGM() = BGM_MeshSize(ge,t,-12,here->x(),here->y(),here->z());
+ }
+ else
+ {
+ MEdgeVertex *eve = (MEdgeVertex*) here;
+ pp->lcBGM() = eve->getLc();
+ }
+ pp->lc() = pp->lcBGM();
}
Msg(DEBUG1,"Meshing of the convex hull (%d points) done",all_vertices.size());
@@ -1189,19 +1230,22 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
double U,V;
SPoint2 param = coords [i];
U = param.x() / m->scalingU ;
- V = param.y() / m->scalingV;
- BDS_Point *pp;
- pp = m->add_point ( count, U,V,gf );
-// if(ge->dim() == 1)
-// {
-// double t;
-// here->getParameter(0,t);
-// pp->lc() = BGM_MeshSize(ge,t,-12,here->x(),here->y(),here->z());
-// }
-// else
-// {
-// pp->lc() = BGM_MeshSize(ge,param.x(),param.y(),here->x(),here->y(),here->z());
-// }
+ V = param.y() / m->scalingV;
+ BDS_Point *pp = m->add_point ( count, U,V,gf );
+ if(ge->dim() == 1)
+ {
+ double t;
+ here->getParameter(0,t);
+ pp->lcBGM() = BGM_MeshSize(ge,t,-12,here->x(),here->y(),here->z());
+ pp->lc() = pp->lcBGM();
+ }
+ else
+ {
+ MEdgeVertex *eve = (MEdgeVertex*) here;
+ // pp->lc() = BGM_MeshSize(ge,param.x(),param.y(),here->x(),here->y(),here->z());
+ pp->lc() = eve->getLc();
+ pp->lcBGM() = eve->getLc();
+ }
m->add_geom (ge->tag(), ge->dim());
BDS_GeomEntity *g = m->get_geom(ge->tag(),ge->dim());
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 3d99547357..99b1cadd6c 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.5 2007-09-10 13:37:21 remacle Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.6 2007-10-10 08:49:34 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -319,7 +319,7 @@ bool insertVertex (MVertex *v ,
vSizesBGM [t->getVertex(1)->getNum()]+
vSizesBGM [t->getVertex(2)->getNum()]);
- MTri3 *t4 = new MTri3 ( t , std::min(lc,lcBGM));
+ MTri3 *t4 = new MTri3 ( t , Extend1dMeshIn2dSurfaces() ? std::min(lc,lcBGM) : lcBGM );
// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
// Us [t->getVertex(0)->getNum()],
// Vs [t->getVertex(0)->getNum()],
diff --git a/benchmarks/2d/Square-Attr2.geo b/benchmarks/2d/Square-Attr2.geo
index a69af148a7..9ca3898285 100644
--- a/benchmarks/2d/Square-Attr2.geo
+++ b/benchmarks/2d/Square-Attr2.geo
@@ -36,4 +36,4 @@ Line(15) = {19,20};
Attractor Line{7,8,9,10,11,12,13,14,15} = {.04,lc/30,lc,1000,7} ;
// Argh:
//Attractor Line{7,8,9,10,11,12,13,14,15} = {.2,0.002,10} ;
-Mesh.Algorithm = 2;
+//Mesh.Algorithm = 2;
--
GitLab