Skip to content
Snippets Groups Projects
Commit 0628eba0 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent cb7d16a9
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/MVertex.h
+ 4
− 3
View file @ 0628eba0
......@@ -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{
......
Mesh/BackgroundMesh.cpp
+ 18
− 2
View file @ 0628eba0
// $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 ();
}
Mesh/BackgroundMesh.h
+ 2
− 1
View file @ 0628eba0
......@@ -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
Mesh/meshGEdge.cpp
+ 5
− 3
View file @ 0628eba0
// $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 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 {
......
Mesh/meshGFace.cpp
+ 166
− 122
View file @ 0628eba0
// $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,76 +297,34 @@ void OptimizeMesh(GFace *gf, BDS_Mesh &m, const int NIT)
}
}
void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
{
int IT =0;
// printf("lc (1,1) = %g\n",Attractor::lc(1,1,0));
int MAXNP = m.MAXPOINTNUMBER;
// computecharacteristic lengths using mesh edge spacing
// those lengths will propagate in the 2D mesh.
if (NIT > 0)
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
{
std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
double L = 1.e22;
while(it!=ite){
double l = (*it)->length();
if (l<L && (*it)->g && (*it)->g->classif_degree == 1)L=l;
++it;
}
if(!CTX.mesh.constrained_bgmesh)
(*itp)->lc() = L;
(*itp)->lcBGM() = L;
++itp;
}
}
double OLDminL=1.E22,OLDmaxL=0;
const double MINE_ = 0.7, MAXE_=1.4;
while (1)
{
// we count the number of local mesh modifs.
int nb_split=0;
int nb_smooth=0;
int nb_collaps=0;
int nb_swap=0;
// 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();
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)
{
(*it)->p1->config_modified = false;
(*it)->p2->config_modified = false;
double lone = NewGetLc ( *it,gf);
maxL = std::max(maxL,lone);
minL = std::min(minL,lone);
}
int result = edgeSwapTestQuality(*it,5);
if (result >= 0)
if(edgeSwapTestDelaunay(*it,gf) || result > 0)
if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
nb_swap++;
++it;
}
}
}
if (OLDminL == minL && OLDmaxL == maxL)break;
OLDminL = minL;OLDmaxL = maxL;
if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT)))break;
NN1 = m.edges.size();
NN2 = 0;
it = m.edges.begin();
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;
......@@ -380,6 +338,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
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,
......@@ -393,32 +352,14 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
++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)
void collapseEdgePass ( GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
{
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();
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1)
{
if (NN2++ >= NN1)break;
......@@ -437,43 +378,129 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
++it;
}
}
// swap edges that provide a better configuration
NN1 = m.edges.size();
NN2 = 0;
it = m.edges.begin();
while (1)
void smoothVertexPass ( GFace *gf, BDS_Mesh &m, int &nb_smooth)
{
if (NN2++ >= NN1)break;
if (!(*it)->deleted)
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
{
int result = edgeSwapTestQuality(*it,5);
if (result >= 0)
if(edgeSwapTestDelaunay(*it,gf) || result > 0)
if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
nb_swap++;
++it;
if(m.smooth_point_centroid(*itp,gf))
nb_smooth ++;
++itp;
}
}
// smooth resulting mesh
void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
{
int IT =0;
// printf("lc (1,1) = %g\n",Attractor::lc(1,1,0));
int MAXNP = m.MAXPOINTNUMBER;
// computecharacteristic lengths using 1D mesh edge spacing
// those lengths will propagate in the 2D mesh.
if (0 && NIT > 0)
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
{
if(m.smooth_point_centroid(*itp,gf))
nb_smooth ++;
std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
double L = 1.e22;
while(it!=ite){
double l = (*it)->length();
if (l<L && (*it)->g && (*it)->g->classif_degree == 1)L=l;
++it;
}
(*itp)->lc() = L;
(*itp)->lcBGM() = L;
++itp;
}
}
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)
{
// we count the number of local mesh modifs.
int nb_split=0;
int nb_smooth=0;
int nb_collaps=0;
int nb_swap=0;
// 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();
while (1)
{
if (NN2++ >= NN1)break;
if (!(*it)->deleted)
{
(*it)->p1->config_modified = false;
(*it)->p2->config_modified = false;
double lone = NewGetLc ( *it,gf);
maxL = std::max(maxL,lone);
minL = std::min(minL,lone);
}
++it;
}
if (OLDminL == minL && OLDmaxL == maxL)break;
OLDminL = minL;OLDmaxL = maxL;
if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT)))break;
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
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;
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());
......@@ -1190,18 +1231,21 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
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
// {
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());
......
Mesh/meshGFaceDelaunayInsertion.cpp
+ 2
− 2
View file @ 0628eba0
// $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()],
......
benchmarks/2d/Square-Attr2.geo
+ 1
− 1
View file @ 0628eba0
......@@ -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;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment