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

2D delaunay & frontal meshers are substantially faster

parent abe69c4b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/meshGFaceDelaunayInsertion.cpp
+ 83
− 40
View file @ 0cabe2d1
...@@ -572,30 +572,21 @@ double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs) ...@@ -572,30 +572,21 @@ double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
return s * 0.5; return s * 0.5;
} }
bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t, bool insertVertexB (std::list<edgeXface> &shell,
std::set<MTri3*, compareTri3Ptr> &allTets, std::list<MTri3*> &cavity,
std::set<MTri3*, compareTri3Ptr> *activeTets, bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
std::vector<double> &vSizes, std::set<MTri3*, compareTri3Ptr> &allTets,
std::vector<double> &vSizesBGM, std::set<MTri3*, compareTri3Ptr> *activeTets,
std::vector<SMetric3> &vMetricsBGM, std::vector<double> &vSizes,
std::vector<double> &Us, std::vector<double> &vSizesBGM,
std::vector<double> &Vs, std::vector<SMetric3> &vMetricsBGM,
double *metric, std::vector<double> &Us,
MTri3 **oneNewTriangle) std::vector<double> &Vs,
double *metric,
MTri3 **oneNewTriangle)
{ {
std::list<edgeXface> shell;
std::list<MTri3*> cavity;
std::list<MTri3*> new_cavity; std::list<MTri3*> new_cavity;
if (!metric){
double p[3] = {v->x(), v->y(), v->z()};
recurFindCavity(shell, cavity, p, param, t, Us, Vs);
}
else{
recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);
}
// check that volume is conserved // check that volume is conserved
double newVolume = 0; double newVolume = 0;
double oldVolume = 0; double oldVolume = 0;
...@@ -699,6 +690,41 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t, ...@@ -699,6 +690,41 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
} }
} }
bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
std::set<MTri3*, compareTri3Ptr> &allTets,
std::set<MTri3*, compareTri3Ptr> *activeTets,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
std::vector<SMetric3> &vMetricsBGM,
std::vector<double> &Us,
std::vector<double> &Vs,
double *metric,
MTri3 **oneNewTriangle)
{
std::list<edgeXface> shell;
std::list<MTri3*> cavity;
if (!metric){
double p[3] = {v->x(), v->y(), v->z()};
recurFindCavity(shell, cavity, p, param, t, Us, Vs);
}
else{
recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);
}
return insertVertexB(shell, cavity, force, gf, v, param , t,
allTets,
activeTets,
vSizes,
vSizesBGM,
vMetricsBGM,
Us,
Vs,
metric,
oneNewTriangle);
}
static MTri3* search4Triangle (MTri3 *t, double pt[2], static MTri3* search4Triangle (MTri3 *t, double pt[2],
std::vector<double> &Us, std::vector<double> &Vs, std::vector<double> &Us, std::vector<double> &Vs,
std::set<MTri3*,compareTri3Ptr> &AllTris, double uv[2], bool force = false) { std::set<MTri3*,compareTri3Ptr> &AllTris, double uv[2], bool force = false) {
...@@ -771,16 +797,36 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it ...@@ -771,16 +797,36 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
} }
} }
else worst = *it; else worst = *it;
MTri3 *ptin = 0; MTri3 *ptin = 0;
bool inside = false; std::list<edgeXface> shell;
std::list<MTri3*> cavity;
double uv[2]; double uv[2];
// printf("inserting %g %g\n",center[0],center[1]);
ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv, oneNewTriangle ? true : false);
if (ptin) inside = true;
if (inside) { // if the point is able to break the bad triangle "worst"
if (1){
if (inCircumCircleAniso(gf, worst->tri(), center, metric, Us, Vs)){
recurFindCavityAniso(gf, shell, cavity, metric, center, worst, Us, Vs);
for (std::list<MTri3*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
if (invMapUV((*itc)->tri(), center, Us, Vs, uv, 1.e-8)) {
ptin = *itc;
break;
}
}
}
// else look for it
else {
printf("we should never be here !!!\n");
ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv, oneNewTriangle ? true : false);
if (ptin) {
recurFindCavityAniso(gf, shell, cavity, metric, center, ptin, Us, Vs);
}
}
}
// ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv, oneNewTriangle ? true : false);
if (ptin) {
// we use here local coordinates as real coordinates // we use here local coordinates as real coordinates
// x,y and z will be computed hereafter // x,y and z will be computed hereafter
// Msg::Info("Point is inside"); // Msg::Info("Point is inside");
...@@ -789,16 +835,11 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it ...@@ -789,16 +835,11 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]); MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]);
v->setIndex(Us.size()); v->setIndex(Us.size());
double lc1,lc; double lc1,lc;
if (0 && backgroundMesh::current()){ lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getIndex()] +
lc1 = lc = uv[0] * vSizes[ptin->tri()->getVertex(1)->getIndex()] +
backgroundMesh::current()->operator()(center[0], center[1], 0.0); uv[1] * vSizes[ptin->tri()->getVertex(2)->getIndex()]);
} lc = BGM_MeshSize(gf, center[0], center[1], p.x(), p.y(), p.z());
else {
lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getIndex()] +
uv[0] * vSizes[ptin->tri()->getVertex(1)->getIndex()] +
uv[1] * vSizes[ptin->tri()->getVertex(2)->getIndex()]);
lc = BGM_MeshSize(gf, center[0], center[1], p.x(), p.y(), p.z());
}
//SMetric3 metr = BGM_MeshMetric(gf, center[0], center[1], p.x(), p.y(), p.z()); //SMetric3 metr = BGM_MeshMetric(gf, center[0], center[1], p.x(), p.y(), p.z());
// vMetricsBGM.push_back(metr); // vMetricsBGM.push_back(metr);
vSizesBGM.push_back(lc); vSizesBGM.push_back(lc);
...@@ -809,8 +850,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it ...@@ -809,8 +850,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
// clock_t t1 = clock(); // clock_t t1 = clock();
if(!p.succeeded() || !insertVertex if(!p.succeeded() || !insertVertexB
(false, gf, v, center, ptin, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM, (shell, cavity,false, gf, v, center, ptin, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric, oneNewTriangle) ) { Us, Vs, metric, oneNewTriangle) ) {
// Msg::Debug("Point %g %g cannot be inserted because %d", // Msg::Debug("Point %g %g cannot be inserted because %d",
// center[0], center[1], p.succeeded() ); // center[0], center[1], p.succeeded() );
...@@ -820,6 +861,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it ...@@ -820,6 +861,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
worst->forceRadius(-1); worst->forceRadius(-1);
AllTris.insert(worst); AllTris.insert(worst);
delete v; delete v;
for (std::list<MTri3*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc)(*itc)->setDeleted(false);
return false; return false;
} }
else { else {
...@@ -832,6 +874,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it ...@@ -832,6 +874,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
} }
else { else {
MTriangle *base = worst->tri(); MTriangle *base = worst->tri();
for (std::list<MTri3*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc)(*itc)->setDeleted(false);
Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g)", Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g)",
center[0], center[1], center[0], center[1],
Us[base->getVertex(0)->getIndex()], Us[base->getVertex(0)->getIndex()],
......
benchmarks/2d/conge.geo
+ 1
− 1
View file @ 0cabe2d1
...@@ -78,5 +78,5 @@ Physical Line(25) = {12,13,14,15,16,17,18,19,20,10}; ...@@ -78,5 +78,5 @@ Physical Line(25) = {12,13,14,15,16,17,18,19,20,10};
Physical Surface(26) = {24}; Physical Surface(26) = {24};
Physical Line(27) = {10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 14, 13, 12, 11}; Physical Line(27) = {10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 14, 13, 12, 11};
Physical Line(28) = {17, 16, 20, 19, 18, 15}; Physical Line(28) = {17, 16, 20, 19, 18, 15};
//Recombine Surface {24, 22}; Recombine Surface {24, 22};
Mesh.RecombinationAlgorithm=1; Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
benchmarks/2d/naca12_2d.geo
+ 7
− 6
View file @ 0cabe2d1
...@@ -228,15 +228,15 @@ Plane Surface(11) = {9,10}; ...@@ -228,15 +228,15 @@ Plane Surface(11) = {9,10};
//Point(9999) = {0.6,0,0,1}; //Point(9999) = {0.6,0,0,1};
Field[2] = BoundaryLayer; Field[2] = BoundaryLayer;
Field[2].NodesList = {1}; //Field[2].NodesList = {1};
//Field[2].EdgesList = {1,2,3,4}; //Field[2].EdgesList = {1,2,3,4};
Field[2].EdgesList = {1,2,3,4}; Field[2].EdgesList = {1,2,3,4};
Field[2].hfar = 1.5; Field[2].hfar = 1.5;
Field[2].hwall_n = 0.0001; Field[2].hwall_n = 0.001;
Field[2].hwall_t = 0.01; Field[2].hwall_t = 0.01;
Field[2].ratio = 1.3; Field[2].ratio = 1.3;
Field[2].thickness = .01; Field[2].thickness = .1;
Background Field = 2; //Background Field = 2;
Field[1] = Box; Field[1] = Box;
Field[1].VIn = 0.01; Field[1].VIn = 0.01;
...@@ -250,5 +250,6 @@ Field[1].ZMin = -1; ...@@ -250,5 +250,6 @@ Field[1].ZMin = -1;
Field[3] = MinAniso; Field[3] = MinAniso;
Field[3].FieldsList = {1, 2}; Field[3].FieldsList = {1, 2};
Background Field = 3; Background Field = 2;
BoundaryLayer Field = 2; //BoundaryLayer Field = 2;
//Recombine Surface {11};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment