Skip to content
Snippets Groups Projects
Commit f39f118d authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

pp

parent 82ad0b39
No related branches found
No related tags found
No related merge requests found
...@@ -32,13 +32,9 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2) ...@@ -32,13 +32,9 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
return l; return l;
} }
inline double computeEdgeLinearLength(BDS_Point *p1, inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
BDS_Point *p2, double SCALINGU, double SCALINGV)
GFace *f,
double SCALINGU,
double SCALINGV)
{ {
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU, GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
0.5 * (p1->v + p2->v) * SCALINGV)); 0.5 * (p1->v + p2->v) * SCALINGV));
...@@ -56,11 +52,8 @@ inline double computeEdgeLinearLength(BDS_Point *p1, ...@@ -56,11 +52,8 @@ inline double computeEdgeLinearLength(BDS_Point *p1,
return l1 + l2; return l1 + l2;
} }
inline double computeEdgeLinearLength_new(BDS_Point *p1, inline double computeEdgeLinearLength_new(BDS_Point *p1, BDS_Point *p2, GFace *f,
BDS_Point *p2, double SCALINGU, double SCALINGV)
GFace *f,
double SCALINGU,
double SCALINGV)
{ {
const int nbSb = 10; const int nbSb = 10;
GPoint GP[nbSb-1]; GPoint GP[nbSb-1];
...@@ -96,13 +89,8 @@ inline double computeEdgeLinearLength_new(BDS_Point *p1, ...@@ -96,13 +89,8 @@ inline double computeEdgeLinearLength_new(BDS_Point *p1,
return l; return l;
} }
inline double computeEdgeMiddleCoord_new(BDS_Point *p1, BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
inline double computeEdgeMiddleCoord_new(BDS_Point *p1,
BDS_Point *p2,
GFace *f,
double SCALINGU,
double SCALINGV)
{ {
const int nbSb = 3; const int nbSb = 3;
double L = computeEdgeLinearLength(p1,p2); double L = computeEdgeLinearLength(p1,p2);
...@@ -128,11 +116,8 @@ inline double computeEdgeMiddleCoord_new(BDS_Point *p1, ...@@ -128,11 +116,8 @@ inline double computeEdgeMiddleCoord_new(BDS_Point *p1,
return 0.5; return 0.5;
} }
inline double computeEdgeMiddleCoord(BDS_Point *p1, inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
BDS_Point *p2, double SCALINGU, double SCALINGV)
GFace *f,
double SCALINGU,
double SCALINGV)
{ {
if (f->geomType() == GEntity::Plane) if (f->geomType() == GEntity::Plane)
return 0.5; return 0.5;
...@@ -155,10 +140,8 @@ inline double computeEdgeMiddleCoord(BDS_Point *p1, ...@@ -155,10 +140,8 @@ inline double computeEdgeMiddleCoord(BDS_Point *p1,
return 0.25 * (3 * l2 - l1) / l2; return 0.25 * (3 * l2 - l1) / l2;
} }
inline double computeEdgeLinearLength(BDS_Edge *e, inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
GFace *f, double SCALINGU, double SCALINGV)
double SCALINGU,
double SCALINGV)
{ {
if (f->geomType() == GEntity::Plane) if (f->geomType() == GEntity::Plane)
return e->length(); return e->length();
...@@ -173,15 +156,13 @@ double NewGetLc(BDS_Point *p) ...@@ -173,15 +156,13 @@ double NewGetLc(BDS_Point *p)
} }
static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f, static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV){ double SCALINGU, double SCALINGV)
{
double l1 = NewGetLc(p1); double l1 = NewGetLc(p1);
double l2 = NewGetLc(p2); double l2 = NewGetLc(p2);
double l = 0.5 * (l1 + l2); double l = 0.5 * (l1 + l2);
// printf(" %g %g -- %g %g\n",SCALINGU,SCALINGV,l1,l2); if(CTX::instance()->mesh.lcFromCurvature){
if(CTX::instance()->mesh.lcFromCurvature)
{
// GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU, // GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
// 0.5 * (p1->v + p2->v) * SCALINGV)); // 0.5 * (p1->v + p2->v) * SCALINGV));
// double l3 = BGM_MeshSize(f,GP.u(),GP.v(),GP.x(),GP.y(),GP.z()); // double l3 = BGM_MeshSize(f,GP.u(),GP.v(),GP.x(),GP.y(),GP.z());
...@@ -282,7 +263,8 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m) ...@@ -282,7 +263,8 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
BDS_Point *p21, *p22, *p23; BDS_Point *p21, *p22, *p23;
BDS_Point *p31, *p32, *p33; BDS_Point *p31, *p32, *p33;
BDS_Point *p41, *p42, *p43; BDS_Point *p41, *p42, *p43;
swap_config(e, &p11, &p12, &p13, &p21, &p22, &p23, &p31, &p32, &p33, &p41, &p42, &p43); swap_config(e, &p11, &p12, &p13, &p21, &p22, &p23, &p31, &p32, &p33, &p41,
&p42, &p43);
// First, evaluate what we gain in element quality if the // First, evaluate what we gain in element quality if the
// swap is performed // swap is performed
...@@ -365,7 +347,6 @@ bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf) ...@@ -365,7 +347,6 @@ bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
return result > 0.; return result > 0.;
} }
bool edgeSwapTestHighOrder(BDS_Edge *e,GFace *gf) bool edgeSwapTestHighOrder(BDS_Edge *e,GFace *gf)
{ {
// must evaluate the swap with the perspectve of // must evaluate the swap with the perspectve of
...@@ -375,7 +356,6 @@ bool edgeSwapTestHighOrder(BDS_Edge *e,GFace *gf) ...@@ -375,7 +356,6 @@ bool edgeSwapTestHighOrder(BDS_Edge *e,GFace *gf)
return false; return false;
} }
bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &configs) bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &configs)
{ {
BDS_Point *op[2]; BDS_Point *op[2];
...@@ -508,10 +488,12 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split) ...@@ -508,10 +488,12 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
// m.scalingU, m.scalingV); // m.scalingU, m.scalingV);
BDS_Point *mid; BDS_Point *mid;
GPoint gpp = gf->point(m.scalingU*(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u), GPoint gpp = gf->point
(m.scalingU*(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u),
m.scalingV*(coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)); m.scalingV*(coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v));
if (gpp.succeeded()){ if (gpp.succeeded()){
mid = m.add_point(++m.MAXPOINTNUMBER, mid = m.add_point
(++m.MAXPOINTNUMBER,
coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u, coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v, gf); coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v, gf);
mid->lcBGM() = BGM_MeshSize mid->lcBGM() = BGM_MeshSize
...@@ -529,15 +511,10 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split) ...@@ -529,15 +511,10 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
} }
} }
/* // A test for spheres
A test for spheres static void midpointsphere(GFace *gf, double u1, double v1, double u2,
*/ double v2, double &u12, double &v12, double r)
{
static void midpointsphere (GFace *gf,
double u1,
double v1,
double u2,
double v2, double &u12, double &v12, double r){
GPoint p1 = gf->point(u1, v1); GPoint p1 = gf->point(u1, v1);
GPoint p2 = gf->point(u2, v2); GPoint p2 = gf->point(u2, v2);
double guess [2] = {0.5 * (u1 + u2), 0.5 * (v1 + v2)}; double guess [2] = {0.5 * (u1 + u2), 0.5 * (v1 + v2)};
...@@ -562,10 +539,10 @@ static void midpointsphere (GFace *gf, ...@@ -562,10 +539,10 @@ static void midpointsphere (GFace *gf,
} }
return; return;
printf("%g %g -- %g %g -- %g %g -- %g %g\n",u1,v1,u2,v2,u12,v12,0.5*(u1+u2),0.5*(v1+v2)); printf("%g %g -- %g %g -- %g %g -- %g %g\n",
u1, v1, u2, v2, u12, v12, 0.5 * (u1 + u2), 0.5 * (v1 + v2));
} }
void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split) void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
{ {
std::list<BDS_Edge*>::iterator it = m.edges.begin(); std::list<BDS_Edge*>::iterator it = m.edges.begin();
...@@ -689,14 +666,11 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q) ...@@ -689,14 +666,11 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const bool computeNodalSizeField) const bool computeNodalSizeField)
{ {
int IT = 0; int IT = 0;
int MAXNP = m.MAXPOINTNUMBER; int MAXNP = m.MAXPOINTNUMBER;
// classify correctly the embedded vertices // classify correctly the embedded vertices use a negative model
// use a negative model face number to avoid // face number to avoid mesh motion
// mesh motion
std::list<GVertex*> emb_vertx = gf->embeddedVertices(); std::list<GVertex*> emb_vertx = gf->embeddedVertices();
std::list<GVertex*>::iterator itvx = emb_vertx.begin(); std::list<GVertex*>::iterator itvx = emb_vertx.begin();
while(itvx != emb_vertx.end()){ while(itvx != emb_vertx.end()){
...@@ -709,10 +683,8 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, ...@@ -709,10 +683,8 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
++itvx; ++itvx;
} }
// If asked, compute nodal size field using 1D Mesh
// IF ASKED , compute nodal size field using 1D Mesh
if (computeNodalSizeField){ if (computeNodalSizeField){
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin(); std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end()){ while (itp != m.points.end()){
std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin(); std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
...@@ -766,7 +738,6 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, ...@@ -766,7 +738,6 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
++it; ++it;
} }
if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT))) break; if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT))) break;
double maxE = MAXE_; double maxE = MAXE_;
double minE = MINE_; double minE = MINE_;
...@@ -801,7 +772,6 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, ...@@ -801,7 +772,6 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
} }
double t_total = t_spl + t_sw + t_col + t_sm; double t_total = t_spl + t_sw + t_col + t_sm;
if(!t_total) t_total = 1.e-6; if(!t_total) t_total = 1.e-6;
Msg::Debug(" ---------------------------------------"); Msg::Debug(" ---------------------------------------");
...@@ -820,7 +790,7 @@ void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2) ...@@ -820,7 +790,7 @@ void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2)
{ {
} }
void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_map, void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recoverMap,
std::set<BDS_Edge*> &toSplit) std::set<BDS_Edge*> &toSplit)
{ {
// first look for degenerated vertices // first look for degenerated vertices
...@@ -829,10 +799,10 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma ...@@ -829,10 +799,10 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma
while (it != m.edges.end()){ while (it != m.edges.end()){
BDS_Edge *e = *it; BDS_Edge *e = *it;
if (!e->deleted && e->numfaces() == 1){ if (!e->deleted && e->numfaces() == 1){
std::map<BDS_Point*,MVertex*>::iterator itp1 = recover_map->find(e->p1); std::map<BDS_Point*,MVertex*>::iterator itp1 = recoverMap->find(e->p1);
std::map<BDS_Point*,MVertex*>::iterator itp2 = recover_map->find(e->p2); std::map<BDS_Point*,MVertex*>::iterator itp2 = recoverMap->find(e->p2);
if (itp1 != recover_map->end() && if (itp1 != recoverMap->end() &&
itp2 != recover_map->end() && itp2 != recoverMap->end() &&
itp1->second == itp2->second){ itp1->second == itp2->second){
degenerated.insert(itp1->second); degenerated.insert(itp1->second);
} }
...@@ -846,12 +816,12 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma ...@@ -846,12 +816,12 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma
while (it != m.edges.end()){ while (it != m.edges.end()){
BDS_Edge *e = *it; BDS_Edge *e = *it;
if (!e->deleted && e->numfaces() == 2){ if (!e->deleted && e->numfaces() == 2){
std::map<BDS_Point*,MVertex*>::iterator itp1 = recover_map->find(e->p1); std::map<BDS_Point*,MVertex*>::iterator itp1 = recoverMap->find(e->p1);
std::map<BDS_Point*,MVertex*>::iterator itp2 = recover_map->find(e->p2); std::map<BDS_Point*,MVertex*>::iterator itp2 = recoverMap->find(e->p2);
if (itp1 != recover_map->end() && if (itp1 != recoverMap->end() &&
itp2 != recover_map->end() && itp2 != recoverMap->end() &&
itp1->second == itp2->second) toSplit.insert(e); itp1->second == itp2->second) toSplit.insert(e);
else if (itp1 != recover_map->end() && itp2 == recover_map->end()){ else if (itp1 != recoverMap->end() && itp2 == recoverMap->end()){
std::pair<MVertex*,BDS_Point*> a ( itp1->second, e->p2 ); std::pair<MVertex*,BDS_Point*> a ( itp1->second, e->p2 );
std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf = std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf =
touchPeriodic.find(a); touchPeriodic.find(a);
...@@ -860,7 +830,7 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma ...@@ -860,7 +830,7 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma
toSplit.insert(e); toSplit.insert(itf->second); toSplit.insert(e); toSplit.insert(itf->second);
} }
} }
else if (itp1 == recover_map->end() && itp2 != recover_map->end()){ else if (itp1 == recoverMap->end() && itp2 != recoverMap->end()){
std::pair<MVertex*,BDS_Point*> a(itp2->second, e->p1); std::pair<MVertex*,BDS_Point*> a(itp2->second, e->p1);
std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf = std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf =
touchPeriodic.find (a); touchPeriodic.find (a);
...@@ -882,10 +852,10 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma ...@@ -882,10 +852,10 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*,MVertex*> *recover_ma
// if p1 p2 exists and it is internal, split it // if p1 p2 exists and it is internal, split it
int gmshSolveInvalidPeriodic(GFace *gf, BDS_Mesh &m, int gmshSolveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
std::map<BDS_Point*,MVertex*> *recover_map) std::map<BDS_Point*,MVertex*> *recoverMap)
{ {
std::set<BDS_Edge*> toSplit; std::set<BDS_Edge*> toSplit;
invalidEdgesPeriodic(m, recover_map, toSplit); invalidEdgesPeriodic(m, recoverMap, toSplit);
std::set<BDS_Edge*>::iterator ite = toSplit.begin(); std::set<BDS_Edge*>::iterator ite = toSplit.begin();
for (;ite !=toSplit.end(); ++ite){ for (;ite !=toSplit.end(); ++ite){
...@@ -911,7 +881,7 @@ int gmshSolveInvalidPeriodic(GFace *gf, BDS_Mesh &m, ...@@ -911,7 +881,7 @@ int gmshSolveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
} }
void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
std::map<BDS_Point*,MVertex*> *recover_map=0) std::map<BDS_Point*,MVertex*> *recoverMap=0)
{ {
int nb_swap; int nb_swap;
gmshDelaunayizeBDS(gf, m, nb_swap); gmshDelaunayizeBDS(gf, m, nb_swap);
...@@ -934,8 +904,8 @@ void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, ...@@ -934,8 +904,8 @@ void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
} }
} }
if (recover_map){ if (recoverMap){
while(gmshSolveInvalidPeriodic(gf, m, recover_map)){ while(gmshSolveInvalidPeriodic(gf, m, recoverMap)){
} }
} }
} }
......
...@@ -23,7 +23,7 @@ void computeElementShapes(GFace *gf, BDS_Mesh &m, double &worst, double &avg, ...@@ -23,7 +23,7 @@ void computeElementShapes(GFace *gf, BDS_Mesh &m, double &worst, double &avg,
void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const bool computeNodalSizeField); const bool computeNodalSizeField);
void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
std::map<BDS_Point*,MVertex*> *recover_map=0); std::map<BDS_Point*,MVertex*> *recoverMap=0);
void gmshDelaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap); void gmshDelaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap);
void gmshCollapseSmallEdges(GModel &gm); void gmshCollapseSmallEdges(GModel &gm);
BDS_Mesh *gmsh2BDS(std::list<GFace*> &l); BDS_Mesh *gmsh2BDS(std::list<GFace*> &l);
......
...@@ -18,20 +18,18 @@ class BDS_Mesh; ...@@ -18,20 +18,18 @@ class BDS_Mesh;
class BDS_Point; class BDS_Point;
void buildMetric(GFace *gf, double *uv, double *metric); void buildMetric(GFace *gf, double *uv, double *metric);
int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3, double *p4, int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
double *metric); double *p4, double *metric);
int inCircumCircleAniso(GFace *gf, MTriangle *base, const double *uv, const double *metric, int inCircumCircleAniso(GFace *gf, MTriangle *base, const double *uv,
const std::vector<double> &Us, const std::vector<double> &Vs); const double *metric, const std::vector<double> &Us,
void circumCenterMetric(double *pa, double *pb, double *pc, const std::vector<double> &Vs);
const double *metric, void circumCenterMetric(double *pa, double *pb, double *pc, const double *metric,
double *x, double &Radius2); double *x, double &Radius2);
void circumCenterMetric(MTriangle *base, void circumCenterMetric(MTriangle *base, const double *metric,
const double *metric,
const std::vector<double> &Us, const std::vector<double> &Us,
const std::vector<double> &Vs, const std::vector<double> &Vs,
double *x, double &Radius2); double *x, double &Radius2);
bool circumCenterMetricInTriangle(MTriangle *base, bool circumCenterMetricInTriangle(MTriangle *base, const double *metric,
const double *metric,
const std::vector<double> &Us, const std::vector<double> &Us,
const std::vector<double> &Vs); const std::vector<double> &Vs);
bool invMapUV(MTriangle *t, double *p, bool invMapUV(MTriangle *t, double *p,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment