From c0952f45229eb352602ce1595a28a0331be338b2 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 27 May 2011 07:48:49 +0000
Subject: [PATCH] New algorithm for quad meshing (delquad)
---
Common/CommandLine.cpp | 4 +-
Common/GmshDefines.h | 1 +
Common/Options.cpp | 3 +
Fltk/optionWindow.cpp | 2 +
Geo/MElementOctree.cpp | 4 +-
Mesh/BackgroundMesh.cpp | 13 +-
Mesh/Field.cpp | 6 -
Mesh/HighOrder.cpp | 2 +-
Mesh/meshGFace.cpp | 8 +-
Mesh/meshGFaceBDS.cpp | 18 +-
Mesh/meshGFaceDelaunayInsertion.cpp | 526 ++++++++++++++++++++--------
Mesh/meshGFaceDelaunayInsertion.h | 3 +
benchmarks/2d/Square-01.geo | 2 +-
benchmarks/2d/conge.geo | 5 +-
benchmarks/2d/recombine.geo | 3 +-
benchmarks/2d/slot.geo | 2 +-
benchmarks/2d/wing-splines.geo | 6 +-
17 files changed, 428 insertions(+), 180 deletions(-)
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 2863b339c3..12b07737ea 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -66,7 +66,7 @@ void PrintUsage(const char *name)
Msg::Direct(" -bin Use binary format when available");
Msg::Direct(" -parametric Save vertices with their parametric coordinates");
Msg::Direct(" -numsubedges Set the number of subdivisions when displaying high order elements");
- Msg::Direct(" -algo string Select mesh algorithm (meshadapt, del2d, front2d, del3d, front3d)");
+ Msg::Direct(" -algo string Select mesh algorithm (meshadapt, del2d, front2d, delquad, del3d, front3d)");
Msg::Direct(" -smooth int Set number of mesh smoothing steps");
Msg::Direct(" -order int Set mesh order (1, ..., 5)");
Msg::Direct(" -optimize[_netgen] Optimize quality of tetrahedral elements");
@@ -540,6 +540,8 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.algo2d = ALGO_2D_MESHADAPT_OLD;
else if(!strncmp(argv[i], "del2d", 5) || !strncmp(argv[i], "tri", 3))
CTX::instance()->mesh.algo2d = ALGO_2D_DELAUNAY;
+ else if(!strncmp(argv[i], "delquad", 5) || !strncmp(argv[i], "tri", 3))
+ CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL_QUAD;
else if(!strncmp(argv[i], "front2d", 7) || !strncmp(argv[i], "frontal", 7))
CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL;
else if(!strncmp(argv[i], "bamg",4))
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 4325ce96d0..7e0d859c57 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -181,6 +181,7 @@
#define ALGO_2D_DELAUNAY 5
#define ALGO_2D_FRONTAL 6
#define ALGO_2D_BAMG 7
+#define ALGO_2D_FRONTAL_QUAD 8
// 3D meshing algorithms (numbers should not be changed)
#define ALGO_3D_DELAUNAY 1
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 9cd4ce697c..82435f710a 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5667,6 +5667,9 @@ double opt_mesh_algo2d(OPT_ARGS_NUM)
case ALGO_2D_FRONTAL:
FlGui::instance()->options->mesh.choice[2]->value(3);
break;
+ case ALGO_2D_FRONTAL_QUAD:
+ FlGui::instance()->options->mesh.choice[2]->value(4);
+ break;
case ALGO_2D_AUTO:
default:
FlGui::instance()->options->mesh.choice[2]->value(0);
diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index 2573a9dbd8..2cedcb5e5c 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -476,6 +476,7 @@ static void mesh_options_ok_cb(Fl_Widget *w, void *data)
(o->mesh.choice[2]->value() == 1) ? ALGO_2D_MESHADAPT :
(o->mesh.choice[2]->value() == 2) ? ALGO_2D_DELAUNAY :
(o->mesh.choice[2]->value() == 3) ? ALGO_2D_FRONTAL :
+ (o->mesh.choice[2]->value() == 4) ? ALGO_2D_FRONTAL_QUAD :
ALGO_2D_AUTO);
opt_mesh_algo3d(0, GMSH_SET,
(o->mesh.choice[3]->value() == 0) ? ALGO_3D_DELAUNAY :
@@ -2051,6 +2052,7 @@ optionWindow::optionWindow(int deltaFontSize)
{"MeshAdapt", 0, 0, 0},
{"Delaunay", 0, 0, 0},
{"Frontal", 0, 0, 0},
+ {"Delaunay for quads", 0, 0, 0},
{0}
};
static Fl_Menu_Item menu_3d_algo[] = {
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index d212da204d..391437909a 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -78,12 +78,12 @@ MElementOctree::MElementOctree(GModel *m) : _gm(m)
Octree_Arrange(_octree);
}
-MElementOctree::MElementOctree(std::vector<MElement*> &v)
+MElementOctree::MElementOctree(std::vector<MElement*> &v) : _gm(0)
{
SBoundingBox3d bb;
for (unsigned int i = 0; i < v.size(); i++){
for(int j = 0; j < v[i]->getNumVertices(); j++){
- if (!_gm) _gm = v[i]->getVertex(j)->onWhat()->model();
+ // if (!_gm) _gm = v[i]->getVertex(j)->onWhat()->model();
bb += SPoint3(v[i]->getVertex(j)->x(),
v[i]->getVertex(j)->y(),
v[i]->getVertex(j)->z());
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index e6188996b8..a38167eb9a 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -369,6 +369,12 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
// take the minimum, then constrain by lcMin and lcMax
double lc = std::min(l1, l2);
+
+ if (backgroundMesh::current()){
+ const double lcBG = backgroundMesh::current()->operator() (U,V,0);
+ lc = std::min(lc,lcBG);
+ }
+
lc = std::max(lc, CTX::instance()->mesh.lcMin);
lc = std::min(lc, CTX::instance()->mesh.lcMax);
@@ -381,7 +387,10 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
SMetric3 LC(1./(lc*lc));
SMetric3 m = intersection(intersection (l4, LC),l3);
// printf("%g %g %g %g %g %g\n",m(0,0),m(1,1),m(2,2),m(0,1),m(0,2),m(1,2));
- return m;
+
+
+
+ return m;
// return lc * CTX::instance()->mesh.lcFactor;
}
@@ -448,7 +457,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
}
}
// ensure that other criteria are fullfilled
- updateSizes(_gf);
+ // updateSizes(_gf);
// compute optimal mesh orientations
propagatecrossField(_gf);
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index cc79e9695c..de40cee7b4 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1715,12 +1715,6 @@ public:
const double ll2 = dist*(ratio-1) + hwall_t;
double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
- if (backgroundMesh::current()){
- const double lcBG = backgroundMesh::current()->operator() (x,y,z);
- lc_n = std::min(lc_n, lcBG);
- lc_t = std::min(lc_t, lcBG);
- }
-
SVector3 t1,t2,t3;
double L1,L2,L3;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index aaa02290e3..8372d0cd8c 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1237,7 +1237,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// printJacobians(m, "smoothness.pos");
// m->writeMSH("SMOOTHED.msh");
- if (!linear){
+ if (0 && !linear){
hot.ensureMinimumDistorsion(0.1);
checkHighOrderTriangles("Final mesh", m, bad, worst);
checkHighOrderTetrahedron("Final mesh", m, bad, worst);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 828d00a090..c848f82ab1 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -300,7 +300,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 ||
+ CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD)
return true;
if(CTX::instance()->mesh.algo2d == ALGO_2D_AUTO && gf->geomType() == GEntity::Plane)
@@ -884,6 +885,8 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
if(algoDelaunay2D(gf)){
if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
bowyerWatsonFrontal(gf);
+ else if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD)
+ bowyerWatsonFrontalQuad(gf);
else if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY ||
CTX::instance()->mesh.algo2d == ALGO_2D_AUTO)
bowyerWatson(gf);
@@ -1473,6 +1476,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
if(algoDelaunay2D(gf)){
if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
bowyerWatsonFrontal(gf);
+ else if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD)
+ bowyerWatsonFrontalQuad(gf);
else if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY ||
CTX::instance()->mesh.algo2d == ALGO_2D_AUTO)
bowyerWatson(gf);
@@ -1544,6 +1549,7 @@ void meshGFace::operator() (GFace *gf)
switch(CTX::instance()->mesh.algo2d){
case ALGO_2D_MESHADAPT : algo = "MeshAdapt"; break;
case ALGO_2D_FRONTAL : algo = "Frontal"; break;
+ case ALGO_2D_FRONTAL_QUAD : algo = "Frontal Quad"; break;
case ALGO_2D_DELAUNAY : algo = "Delaunay"; break;
case ALGO_2D_MESHADAPT_OLD : algo = "MeshAdapt (old)"; break;
case ALGO_2D_BAMG : algo = "Bamg"; break;
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index ef3e7c514c..7066b5f2fe 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -518,27 +518,14 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
std::sort(edges.begin(), edges.end());
- double RADIUS;
- SPoint3 CENTER;
- bool isSphere = gf->isSphere(RADIUS,CENTER);
-
for (unsigned int i = 0; i < edges.size(); ++i){
BDS_Edge *e = edges[i].second;
if (!e->deleted){
const double coord = 0.5;
- //const double coord = computeEdgeMiddleCoord((*it)->p1, (*it)->p2, gf,
- // m.scalingU, m.scalingV);
BDS_Point *mid ;
-
double U, V;
- if (0 && isSphere){
- midpointsphere(gf,e->p1->u,e->p1->v,e->p2->u,e->p2->v,U,V,
- CENTER,RADIUS);
- }
- else{
- U = coord * e->p1->u + (1 - coord) * e->p2->u;
- V = coord * e->p1->v + (1 - coord) * e->p2->v;
- }
+ U = coord * e->p1->u + (1 - coord) * e->p2->u;
+ V = coord * e->p1->v + (1 - coord) * e->p2->v;
GPoint gpp = gf->point(m.scalingU*U,m.scalingV*V);
if (gpp.succeeded()){
@@ -558,7 +545,6 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
(coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
(coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
mid->X,mid->Y,mid->Z);
- // mid->lc() = exp(0.5 * (log(e->p1->lc()) + log(e->p2->lc())));
mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
}
if(!m.split_edge(e, mid)) m.del_point(mid);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 12b11684f8..8f2cc8b99a 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -16,21 +16,66 @@
#include "Numeric.h"
#include "STensor3.h"
-const double LIMIT_ = 0.5 * sqrt(2.0) * 1;
+double LIMIT_ = 0.5 * sqrt(2.0) * 1;
static const bool _experimental_anisotropic_blues_band_ = false;
+int MTri3::radiusNorm = 2;
// This function controls the frontal algorithm
+static bool isBoundary(MTri3 *t, double limit_, int &active){
+ if (t->isDeleted()) return false;
+ for (active = 0; active < 3; active++){
+ MTri3 *neigh = t->getNeigh(active);
+ if (!neigh) return true;
+ }
+ return false;
+}
+
static bool isActive(MTri3 *t, double limit_, int &active)
{
if (t->isDeleted()) return false;
for (active = 0; active < 3; active++){
MTri3 *neigh = t->getNeigh(active);
- if (!neigh || neigh->getRadius() < limit_) return true;
+ if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
+ return true;
+ }
+ }
+ return false;
+}
+
+static bool isActive(MTri3 *t, double limit_, int &i, std::set<MEdge,Less_Edge> *front)
+{
+ if (t->isDeleted()) return false;
+ for (i = 0; i < 3; i++){
+ MTri3 *neigh = t->getNeigh(i);
+ if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
+ int ip1 = i - 1 < 0 ? 2 : i - 1;
+ int ip2 = i;
+ MEdge me (t->tri()->getVertex(ip1),t->tri()->getVertex(ip2));
+ if(front->find(me) != front->end()){
+ return true;
+ }
+ }
}
return false;
}
+
+static void updateActiveEdges(MTri3 *t, double limit_, std::set<MEdge,Less_Edge> &front)
+{
+ if (t->isDeleted()) return;
+ for (int active = 0; active < 3; active++){
+ MTri3 *neigh = t->getNeigh(active);
+ if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
+ int ip1 = active - 1 < 0 ? 2 : active - 1;
+ int ip2 = active;
+ MEdge me (t->tri()->getVertex(ip1),t->tri()->getVertex(ip2));
+ front.insert(me);
+ }
+ }
+}
+
+
bool circumCenterMetricInTriangle(MTriangle *base, const double *metric,
const std::vector<double> &Us,
const std::vector<double> &Vs)
@@ -221,12 +266,29 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric)
// compute the circumradius of the triangle
if (!metric){
- circumCenterXYZ(pa, pb, pc, center);
- const double dx = base->getVertex(0)->x() - center[0];
- const double dy = base->getVertex(0)->y() - center[1];
- const double dz = base->getVertex(0)->z() - center[2];
- circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
- circum_radius /= lc;
+ if (radiusNorm == 2){
+ circumCenterXYZ(pa, pb, pc, center);
+ const double dx = base->getVertex(0)->x() - center[0];
+ const double dy = base->getVertex(0)->y() - center[1];
+ const double dz = base->getVertex(0)->z() - center[2];
+ circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
+ circum_radius /= lc;
+ }
+ else {
+ double quadAngle = 0.0;//backgroundMesh::current() ? backgroundMesh::current()->getAngle ((pa[0]+pb[0]+pc[0])/3.0,(pa[1]+pb[1]+pc[1])/3.0,0) : 0.0;
+ double x0 = base->getVertex(0)->x() * cos(quadAngle) + base->getVertex(0)->y() * sin(quadAngle);
+ double y0 = -base->getVertex(0)->x() * sin(quadAngle) + base->getVertex(0)->y() * cos(quadAngle);
+ double x1 = base->getVertex(1)->x() * cos(quadAngle) + base->getVertex(1)->y() * sin(quadAngle);
+ double y1 = -base->getVertex(1)->x() * sin(quadAngle) + base->getVertex(1)->y() * cos(quadAngle);
+ double x2 = base->getVertex(2)->x() * cos(quadAngle) + base->getVertex(2)->y() * sin(quadAngle);
+ double y2 = -base->getVertex(2)->x() * sin(quadAngle) + base->getVertex(2)->y() * cos(quadAngle);
+ double xmax = std::max(std::max(x0,x1),x2);
+ double ymax = std::max(std::max(y0,y1),y2);
+ double xmin = std::min(std::min(x0,x1),x2);
+ double ymin = std::min(std::min(y0,y1),y2);
+ circum_radius = std::max(xmax-xmin,ymax-ymin);
+ circum_radius /= lc;
+ }
}
else{
double uv[2];
@@ -234,6 +296,7 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric)
}
}
+
int MTri3::inCircumCircle(const double *p) const
{
double pa[3] =
@@ -473,7 +536,7 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
else{
t4 = new MTri3(t, LL);
}
-
+
double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
(it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
(it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
@@ -572,7 +635,22 @@ void _printTris(char *name, std::set<MTri3*, compareTri3Ptr> &AllTris,
fclose (ff);
}
-static void insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it,
+static MTri3* search4Triangle (MTri3 *t, double pt[2],
+ std::vector<double> &Us, std::vector<double> &Vs,
+ std::set<MTri3*,compareTri3Ptr> &AllTris) {
+
+ double uv[2];
+ bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
+ if (inside) return starting_point;
+ while (1){
+ for (int i=0;i<3;i++){
+ MTri3 *tn = t->getNeigh(0);
+
+ }
+ }
+}
+
+static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it,
double center[2], double metric[3],
std::vector<double> &Us, std::vector<double> &Vs,
std::vector<double> &vSizes,
@@ -586,14 +664,13 @@ static void insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
it = AllTris.find(worst);
if (worst != *it){
Msg::Error("Could not insert point");
- return;
+ return false;
}
}
else worst = *it;
MTri3 *ptin = 0;
double uv[2];
- //MTriangle *base = worst->tri();
bool inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
if (inside)ptin = worst;
if (!inside && worst->getNeigh(0)){
@@ -608,7 +685,22 @@ static void insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
inside |= invMapUV(worst->getNeigh(2)->tri(), center, Us, Vs, uv, 1.e-8);
if (inside)ptin = worst->getNeigh(2);
}
- if (inside) {
+
+ // TEST :
+ if (MTri3::radiusNorm == -1 && !ptin){
+ for(std::set<MTri3*,compareTri3Ptr>::iterator itx = AllTris.begin(); itx != AllTris.end();++itx){
+ if (!(*itx)->isDeleted()){
+ inside = invMapUV((*itx)->tri(), center, Us, Vs, uv, 1.e-8);
+ if (inside){
+ ptin = *it;
+ break;
+ }
+ }
+ }
+ }
+
+ // if (!ptin)ptin = worst;
+ if ( inside) {
// we use here local coordinates as real coordinates
// x,y and z will be computed hereafter
// Msg::Info("Point is inside");
@@ -637,61 +729,48 @@ static void insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
if(!p.succeeded() || !insertVertex
(gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric) ) {
- // Msg::Debug("2D Delaunay : a cavity is not star shaped");
+
+ MTriangle *base = worst->tri();
+ /*
+ Msg::Info("Point %g %g of (%g %g , %g %g , %g %g) cannot be inserted",
+ center[0], center[1],
+ Us[base->getVertex(0)->getIndex()],
+ Vs[base->getVertex(0)->getIndex()],
+ Us[base->getVertex(1)->getIndex()],
+ Vs[base->getVertex(1)->getIndex()],
+ Us[base->getVertex(2)->getIndex()],
+ Vs[base->getVertex(2)->getIndex()]);
+ */
AllTris.erase(it);
worst->forceRadius(-1);
AllTris.insert(worst);
delete v;
+ return false;
}
else {
gf->mesh_vertices.push_back(v);
+ return true;
}
}
else {
- /* Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",
- center[0], center[1],
- Us[base->getVertex(0)->getIndex()],
- Vs[base->getVertex(0)->getIndex()],
- Us[base->getVertex(1)->getIndex()],
- Vs[base->getVertex(1)->getIndex()],
- Us[base->getVertex(2)->getIndex()],
- Vs[base->getVertex(2)->getIndex()],
- metric[0], metric[1], metric[2]);*/
+ MTriangle *base = worst->tri();
+ /*
+ Msg::Info("Point %g %g is outside (%g %g , %g %g , %g %g)",
+ center[0], center[1],
+ Us[base->getVertex(0)->getIndex()],
+ Vs[base->getVertex(0)->getIndex()],
+ Us[base->getVertex(1)->getIndex()],
+ Vs[base->getVertex(1)->getIndex()],
+ Us[base->getVertex(2)->getIndex()],
+ Vs[base->getVertex(2)->getIndex()]);
+ */
AllTris.erase(it);
worst->forceRadius(0);
AllTris.insert(worst);
+ return false;
}
}
-
-static void insertManyPoints(GFace *gf,
- std::list<SPoint2> &points,
- std::vector<double> &Us,
- std::vector<double> &Vs,
- std::vector<double> &vSizes,
- std::vector<double> &vSizesBGM,
- std::vector<SMetric3> &vMetricsBGM,
- std::set<MTri3*,compareTri3Ptr> &AllTris,
- std::set<MTri3*,compareTri3Ptr> *ActiveTris = 0){
-
- // a first implementation : greeeeedy algorithm
- for (std::list<SPoint2>::iterator itp = points.begin(); itp != points.end() ; ++itp){
- std::set<MTri3*,compareTri3Ptr> :: iterator it = AllTris.begin();
- double metric[3];
- double pa[2] = {itp->x(),itp->y()};
- buildMetric(gf, pa, metric);
- for (; it != AllTris.end() ; ++it){
- int found = inCircumCircleAniso(gf, (*it)->tri(), pa, metric, Us, Vs);
- if (found){
- insertAPoint(gf, it, pa, metric, Us, Vs, vSizes, vSizesBGM, vMetricsBGM,
- AllTris, ActiveTris);
- break;
- }
- }
- }
-}
-
-
void bowyerWatson(GFace *gf)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
@@ -763,6 +842,43 @@ void bowyerWatson(GFace *gf)
The point isertion is done on the Voronoi Edges
*/
+static double lengthInfniteNorm(const double p[2], const double q[2],
+ const double quadAngle){
+ double xp = p[0] * cos(quadAngle) - p[1] * sin(quadAngle);
+ double yp = p[0] * sin(quadAngle) + p[1] * cos(quadAngle);
+ double xq = q[0] * cos(quadAngle) - q[1] * sin(quadAngle);
+ double yq = q[0] * sin(quadAngle) + q[1] * cos(quadAngle);
+ double xmax = std::max(xp,xq);
+ double xmin = std::min(xp,xq);
+ double ymax = std::max(yp,yq);
+ double ymin = std::min(yp,yq);
+ return std::max(xmax-xmin,ymax-ymin);
+}
+
+void circumCenterInfinite (MTriangle *base, double quadAngle,
+ const std::vector<double> &Us,
+ const std::vector<double> &Vs, double *x) {
+ double pa[2] = {Us[base->getVertex(0)->getIndex()],
+ Vs[base->getVertex(0)->getIndex()]};
+ double pb[2] = {Us[base->getVertex(1)->getIndex()],
+ Vs[base->getVertex(1)->getIndex()]};
+ double pc[2] = {Us[base->getVertex(2)->getIndex()],
+ Vs[base->getVertex(2)->getIndex()]};
+ double xa = pa[0] * cos(quadAngle) + pa[1] * sin(quadAngle);
+ double ya = -pa[0] * sin(quadAngle) + pa[1] * cos(quadAngle);
+ double xb = pb[0] * cos(quadAngle) + pb[1] * sin(quadAngle);
+ double yb = -pb[0] * sin(quadAngle) + pb[1] * cos(quadAngle);
+ double xc = pc[0] * cos(quadAngle) + pc[1] * sin(quadAngle);
+ double yc = -pc[0] * sin(quadAngle) + pc[1] * cos(quadAngle);
+ double xmax = std::max(std::max(xa,xb),xc);
+ double ymax = std::max(std::max(ya,yb),yc);
+ double xmin = std::min(std::min(xa,xb),xc);
+ double ymin = std::min(std::min(ya,yb),yc);
+ x[0] = 0.5 * (xmax-xmin);
+ x[1] = 0.5 * (ymax-ymin);
+}
+
+
static double lengthMetric(const double p[2], const double q[2],
const double metric[3])
{
@@ -912,105 +1028,227 @@ void bowyerWatsonFrontal(GFace *gf)
transferDataStructure(gf, AllTris, Us, Vs);
}
-// give it a try : add one quad layer on the
- /*
-void addOneLayerOnContour(GFace *gf, GVertex *gv){
-, int nbLayers, double hplus, double factor){
- // for each vertex
- std::map<MVertex*,std::vector<MVertex*> >layers;
- std::vector<MQuadrangle*> newQuads;
- std::vector<MTriangle*> newTris;
-
- std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin();
- for (; it != gf->edgeLoops.end(); ++it){
- bool found = false;
- std::list<GEdge*> ed;
- for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
- if (it2->ge->getBeginVertex() == gv || it2->ge->getEndVertex() == gv) {
- found = true;
- }
- ed.push_back(it2->ge);
+
+void bowyerWatsonFrontalQuad(GFace *gf)
+{
+
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ std::set<MTri3*,compareTri3Ptr> ActiveTris;
+ std::vector<double> vSizes, vSizesBGM, Us, Vs;
+ std::vector<SMetric3> vMetricsBGM;
+
+ if (!backgroundMesh::current()) {
+ std::vector<MTriangle*> TR;
+ for(int i=0;i<gf->triangles.size();i++){
+ TR.push_back(new MTriangle(gf->triangles[i]->getVertex(0),
+ gf->triangles[i]->getVertex(1),
+ gf->triangles[i]->getVertex(2)));
}
- // we found an edge loop with the GVertex that was specified
- if (found){
- // compute model vertices that will produce fans
- for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
- GEdgeLoop::iter it3 = it2; ++it3;
- GVertex *gv = it2->getEndVertex();
- GEdgeSigned *before,*after = *it2;
- if (it3 == it->end()){
- before = *(it->begin());
+ bowyerWatson(gf);
+ backgroundMesh::set(gf);
+ char name[256];
+ sprintf(name,"bgm-%d.pos",gf->tag());
+ backgroundMesh::current()->print(name,gf);
+ sprintf(name,"cross-%d.pos",gf->tag());
+ backgroundMesh::current()->print(name,gf,1);
+ // FIXME DELETE CURRENT MESH
+ gf->triangles = TR;
+ }
+
+ LIMIT_ = sqrt(2)*.99;
+ // LIMIT_ = 1.7;
+ MTri3::radiusNorm = -1;
+
+
+ buildMeshGenerationDataStructures
+ (gf, AllTris, vSizes, vSizesBGM, vMetricsBGM,Us, Vs);
+
+ // delaunise the initial mesh
+ int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
+ Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
+
+ int ITER = 0, active_edge;
+ // compute active triangle
+ std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
+ std::set<MEdge,Less_Edge> _front;
+ for ( ; it!=AllTris.end();++it){
+ if(isActive(*it,LIMIT_,active_edge)){
+ ActiveTris.insert(*it);
+ updateActiveEdges(*it, LIMIT_, _front);
+ }
+ else if ((*it)->getRadius() < LIMIT_)break;
+ }
+
+ // insert points
+ int ITERATION = 1;
+ while (1){
+ ITERATION ++;
+ if(ITERATION % 1== 0){
+ char name[245];
+ sprintf(name,"denInfinite_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
+ _printTris (name, AllTris, Us,Vs,true);
+ sprintf(name,"denInfinite_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ _printTris (name, ActiveTris, Us,Vs,true);
+ }
+
+ std::set<MTri3*,compareTri3Ptr> ActiveTrisNotInFront;
+
+ while (1){
+
+ if (!ActiveTris.size())break;
+
+ std::set<MTri3*,compareTri3Ptr>::iterator WORST_ITER = ActiveTris.begin();
+
+ MTri3 *worst = (*WORST_ITER);
+ ActiveTris.erase(WORST_ITER);
+ if (!worst->isDeleted() && isActive(worst, LIMIT_, active_edge,&_front) &&
+ worst->getRadius() > LIMIT_){
+ // for (active_edge = 0 ; active_edge < 0 ; active_edge ++){
+ // if (active_edges[active_edge])break;
+ // }
+ if(ITER++ % 5000 == 0)
+ Msg::Debug("%7d points created -- Worst tri infinite radius is %8.3f",
+ vSizes.size(), worst->getRadius());
+
+ // compute the middle point of the edge
+ MTriangle *base = worst->tri();
+ int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
+ int ip2 = active_edge;
+ int ip3 = (active_edge+1)%3;
+
+ double P[2] = {Us[base->getVertex(ip1)->getIndex()],
+ Vs[base->getVertex(ip1)->getIndex()]};
+ double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
+ Vs[base->getVertex(ip2)->getIndex()]};
+ double O[2] = {Us[base->getVertex(ip3)->getIndex()],
+ Vs[base->getVertex(ip3)->getIndex()]};
+ double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
+
+ // compute background mesh data
+ double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0);
+ // quadAngle = 35*M_PI/180;
+ //double quadAngle = backgroundMesh::current()->getAngle (0.5*(base->getVertex(ip1)->x()+base->getVertex(ip2)->x()),
+ // 0.5*(base->getVertex(ip1)->y()+base->getVertex(ip2)->y()),
+ // 0.5*(base->getVertex(ip1)->x()+base->getVertex(ip2)->x()));
+ // printf("quadAngle = %12.5E\n",quadAngle*180/M_PI);
+ // double meshSize = backgroundMesh::current()->operator()(midpoint[0],midpoint[1],0);
+ //double quadAngle = 0;
+ double center[2];
+ circumCenterInfinite (base, quadAngle,Us,Vs,center);
+
+ // rotate the points with respect to the angle
+ double XP1 = 0.5*(Q[0] - P[0]);
+ double YP1 = 0.5*(Q[1] - P[1]);
+ double xp = XP1 * cos(quadAngle) - YP1 * sin(quadAngle);
+ double yp = XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
+ // ensure xp > yp
+ bool exchange = false;
+ if (fabs(xp) < fabs(yp)){
+ double temp = xp;
+ xp = yp;
+ yp = temp;
+ exchange = true;
}
- else{
- before = *it2;
+
+
+
+ const double p = 0.5 * lengthInfniteNorm(P, Q, quadAngle);
+ const double q = lengthInfniteNorm(center, midpoint, quadAngle);
+ const double rhoM1 = 0.5 *
+ (vSizes[base->getVertex(ip1)->getIndex()] +
+ vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
+ const double rhoM2 = 0.5 *
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ vSizesBGM[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
+ const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
+
+ const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
+ // const double rhoM_hat = std::max(rhoM, p);
+ // assume rhoM = L/\sqrt{3}
+ // assume that p = L/2
+ // d = L/\sqrt{3} + \sqrt{L/3 - L/4} = L/\sqrt{3} + L/2\sqrt{3} = \sqrt{3} L / 2 ... OK
+ const double factor = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) /(sqrt(3)*p);
+ // printf("factor = %g\n",factor);
+
+ double npx,npy;
+ if (xp*yp > 0){
+ npx = - fabs(xp)*factor;
+ npy = fabs(xp)*(1.+factor) - fabs(yp);
}
- }
-
- for (std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); ++it){
- GEdge *ge = *it;
- for (int i=0;i<ge->lines.size();i++){
- SPoint2 p[2];
- reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),gf,p[0],p[1]);
- MVertex *vd[2];
- for (int j=0;j<2;j++){
- MVertex *v = ge->lines[i]->getVertex(j);
- std::map<MVertex*,MVertex*>::iterator itv = duplicates.find(v);
- if (itv == duplicates.end()){
- vd[j] = new MFaceVertex(v->x(),v->y(),v->z(),gf,p[j].x(),p[j].y());
- duplicates[v] = vd[j];
- gf->mesh_vertices.push_back(vd[j]);
- }
- else
- vd[j] = itv->second;
- }
- newQuads.push_back(new MQuadrangle(ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),vd[1],vd[0]));
+ else {
+ npx = fabs(xp) * factor;
+ npy = (1.+factor)*fabs(xp) - fabs(yp);
}
- }
- for (int i=0;i<gf->quadrangles.size();i++){
- MQuadrangle *q = gf->quadrangles[i];
- MVertex *vs[4];
- for (int j=0;j<4;j++){
- MVertex *v = q->getVertex(j);
- std::map<MVertex*,MVertex*>::iterator itv = duplicates.find(v);
- if (itv == duplicates.end()){
- vs[j] = v;
- }
- else{
- vs[j] = itv->second;
+ if (exchange){
+ double temp = npx;
+ npx = npy;
+ npy = temp;
+ }
+
+
+ double newPoint[2] = {midpoint[0] + cos(quadAngle) * npx + sin(quadAngle) * npy,
+ midpoint[1] - sin(quadAngle) * npx + cos(quadAngle) * npy};
+ /*
+ printf("exchange %d\n",exchange);
+ printf("P %g %g\n",P[0],P[1]);
+ printf("Q %g %g\n",Q[0],Q[1]);
+ printf("midpoint %g %g\n",midpoint[0],midpoint[1]);
+ printf("xp yp %g %g\n",xp,yp);
+ printf("O %g %g\n",O[0],O[1]);
+ printf("dx %g %g\n",npx,npy);
+ */
+ if ((midpoint[0] - newPoint[0])*(midpoint[0] - O[0]) +
+ (midpoint[1] - newPoint[1])*(midpoint[1] - O[1]) < 0){
+ newPoint[0] = midpoint[0] - cos(quadAngle) * npx - sin(quadAngle) * npy;
+ newPoint[1] = midpoint[1] + sin(quadAngle) * npx - cos(quadAngle) * npy;
+
+ // printf("wrong sense %g \n",(midpoint[0] - newPoint[0])*(midpoint[0] - O[0]) +
+ // (midpoint[1] - newPoint[1])*(midpoint[1] - O[1]));
+ }
+
+ // printf("new %g %g\n",newPoint[0],newPoint[1]);
+
+
+ insertAPoint(gf, AllTris.end(), newPoint, 0, Us, Vs, vSizes,
+ vSizesBGM, vMetricsBGM, AllTris, &ActiveTris, worst);
+ // else if (!worst->isDeleted() && worst->getRadius() > LIMIT_){
+ // ActiveTrisNotInFront.insert(worst);
+ // }
+
+ /*
+ if(ITER % 100== 0){
+ char name[245];
+ sprintf(name,"frontal%d-ITER%d.pos",gf->tag(),ITER);
+ _printTris (name, AllTris, Us,Vs,false);
}
- }
- newQuads.push_back(new MQuadrangle(vs[0],vs[1],vs[2],vs[3]));
- delete q;
+ */
}
- for (int i=0;i<gf->triangles.size();i++){
- MTriangle *t = gf->triangles[i];
- MVertex *vs[3];
- for (int j=0;j<3;j++){
- MVertex *v = t->getVertex(j);
- std::map<MVertex*,MVertex*>::iterator itv = duplicates.find(v);
- if (itv == duplicates.end()){
- vs[j] = v;
- }
- else{
- vs[j] = itv->second;
- }
- }
- newTris.push_back(new MTriangle(vs[0],vs[1],vs[2]));
- delete t;
+ else if (!worst->isDeleted() && worst->getRadius() > LIMIT_){
+ ActiveTrisNotInFront.insert(worst);
}
-
- gf->triangles = newTris;
- gf->quadrangles = newQuads;
}
+ _front.clear();
+ it = ActiveTrisNotInFront.begin();
+ //it = AllTris.begin();
+ for ( ; it!=ActiveTrisNotInFront.end();++it){
+ // for ( ; it!=AllTris.end();++it){
+ if((*it)->getRadius() > LIMIT_ && isActive(*it,LIMIT_,active_edge)){
+ ActiveTris.insert(*it);
+ updateActiveEdges(*it, LIMIT_, _front);
+ }
+ }
+ Msg::Info("%d active tris %d front edges %d not in front",ActiveTris.size(),_front.size(),ActiveTrisNotInFront.size());
+ if (!ActiveTris.size())break;
}
-}
-*/
+
+ // char name[245];
+ // sprintf(name,"frontal%d-real.pos", gf->tag());
+ // _printTris (name, AllTris, Us, Vs,false);
+ // sprintf(name,"frontal%d-param.pos", gf->tag());
+ // _printTris (name, AllTris, Us, Vs,true);
+ transferDataStructure(gf, AllTris, Us, Vs);
+ MTri3::radiusNorm = 2;
+ LIMIT_ = 0.5 * sqrt(2.0) * 1;
+ backgroundMesh::unset();
+}
-void addBoundaryLayers(GFace *gf)
-{
- if (backgroundMesh::current()){
- }
- // first compute the cross field if it is not computed yet
- // start from a selection of edges and create points in the boundary layer
- // connect everybody with delaunay
-}
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index da5af89261..3dd85c5d03 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -38,12 +38,14 @@ bool invMapUV(MTriangle *t, double *p,
class MTri3
{
+ protected :
bool deleted;
double circum_radius;
MTriangle *base;
MTri3 *neigh[3];
public :
+ static int radiusNorm; // 2 is euclidian norm, -1 is infinite norm
bool isDeleted() const { return deleted; }
void forceRadius(double r) { circum_radius = r; }
double getRadius() const { return circum_radius; }
@@ -94,6 +96,7 @@ void connectTriangles(std::vector<MTri3*> &);
void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
void bowyerWatson(GFace *gf);
void bowyerWatsonFrontal(GFace *gf);
+void bowyerWatsonFrontalQuad(GFace *gf);
struct edgeXface
{
diff --git a/benchmarks/2d/Square-01.geo b/benchmarks/2d/Square-01.geo
index e6aab03099..ba31ef8cd4 100644
--- a/benchmarks/2d/Square-01.geo
+++ b/benchmarks/2d/Square-01.geo
@@ -2,7 +2,7 @@ fact = 100;
lc = .1 * fact;
Point(1) = {0.0,0.0,0,lc};
Point(2) = {1* fact,0.0,0,lc};
-Point(3) = {1* fact,1* fact,0,lc*.1};
+Point(3) = {1* fact,1* fact,0,lc*.5};
Point(4) = {0,1* fact,0,lc};
Line(1) = {3,2};
Line(2) = {2,1};
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index 4ff78396df..60dd0ebd3b 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -21,7 +21,7 @@ Point(1) = { -e1-e2, 0.0 , 0.0 , Lc1};
Point(2) = { -e1-e2, h1 , 0.0 , Lc1};
Point(3) = { -e3-r , h1 , 0.0 , Lc2};
Point(4) = { -e3-r , h1+r , 0.0 , Lc2};
-Point(5) = { -e3 , h1+r , 0.0 , Lc2};
+Point(5) = { -e3 , h1+r , 0.0 , Lc2*.1};
Point(6) = { -e3 , h1+h2, 0.0 , Lc1};
Point(7) = { e3 , h1+h2, 0.0 , Lc1};
Point(8) = { e3 , h1+r , 0.0 , Lc2};
@@ -78,4 +78,5 @@ Physical Line(25) = {12,13,14,15,16,17,18,19,20,10};
Physical Surface(26) = {22,24};
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};
-//Recombine Surface {24, 22};
+Recombine Surface {24, 22};
+Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
diff --git a/benchmarks/2d/recombine.geo b/benchmarks/2d/recombine.geo
index 083f744b62..b0afab1e8b 100644
--- a/benchmarks/2d/recombine.geo
+++ b/benchmarks/2d/recombine.geo
@@ -9,4 +9,5 @@ Line(3) = {3,4} ;
Line(4) = {4,1} ;
Line Loop(5) = {4,1,-2,3} ;
Plane Surface(6) = {5} ;
-Recombine Surface{6};
+//Recombine Surface{6};
+//Mesh.RecombinationAlgorithm=1;
diff --git a/benchmarks/2d/slot.geo b/benchmarks/2d/slot.geo
index 7e4d787585..bc70926c25 100644
--- a/benchmarks/2d/slot.geo
+++ b/benchmarks/2d/slot.geo
@@ -1,4 +1,4 @@
-Mesh.SubdivisionAlgorithm = 3;
+Mesh.RecombinationAlgorithm = 1;
Point(1) = {0, 0, 0, 9.999999999999999e-05};
Point(2) = {0.00125, 0.045983013167908, 0, 0.0002};
Point(3) = {-0.00125, 0.045983013167908, 0, 0.0002};
diff --git a/benchmarks/2d/wing-splines.geo b/benchmarks/2d/wing-splines.geo
index 49d30ddc7b..53c3aa4b63 100644
--- a/benchmarks/2d/wing-splines.geo
+++ b/benchmarks/2d/wing-splines.geo
@@ -2,7 +2,7 @@
scale = .01 ;
lc_wing = 0.05 * scale ;
-lc_box = 10 * scale ;
+lc_box = 30 * scale ;
Point(3895) = {1.177410e-02*scale,-2.768003e-03*scale,0,lc_wing};
Point(3897) = {1.081196e-02*scale,-3.794565e-03*scale,0,lc_wing};
@@ -524,7 +524,7 @@ Spline(1) = {3846,3981,4001,3942,3940,3892,3943,3895,
3897,3896,3968,3995,4003,3857,3856,3860,3861,3863,
3864,3865,3866,3867,3869,3870,3871,3872,3977,
3877,3876,3878,3934,3873,3874,3875,3935,3880,3879,
- 3881,3936,3882,3883,3885,3884,1218,3933,3996,3989,
+ 3881,3936,3882,3883,3885,3884,3933,3996,3989,
3990,3978,3979,3974,3973,3963,3948,3904,3903,
3946,3902,3901,3900,3908,3907,3951,3950,3847};
Spline(2) = {3847,3949,3952,3905,3906,3909,3969,3970,3997,3998,
@@ -604,3 +604,5 @@ Circle(408) = {4355,4351,4354};
Circle(409) = {4354,4351,4353};
Circle(410) = {4353,4351,4352};
*/
+Recombine Surface {406};
+Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
--
GitLab