From 7d08e0731690dfb8e2d7edd0f8eb7e89c2977f97 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Sat, 10 Mar 2012 10:28:47 +0000
Subject: [PATCH] frontal mesher is now fully functional for compound faces
---
Geo/CMakeLists.txt | 1 +
Geo/GFaceCompound.cpp | 101 ++++++++++++++++++++++++++++
Geo/GFaceCompound.h | 2 +
Mesh/CMakeLists.txt | 2 +-
Mesh/meshGFace.cpp | 4 +-
Mesh/meshGFaceDelaunayInsertion.cpp | 75 +++++++++------------
Mesh/meshGFaceOptimize.cpp | 30 ++-------
Mesh/meshGFaceOptimize.h | 1 -
Numeric/Numeric.cpp | 3 +-
Numeric/fullMatrix.cpp | 8 +--
benchmarks/stl/artery.geo | 4 +-
benchmarks/stl/skull.geo | 4 +-
12 files changed, 155 insertions(+), 80 deletions(-)
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 32bbeeb13d..17466e538b 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -4,6 +4,7 @@
# bugs and problems to <gmsh@geuz.org>.
set(SRC
+ intersectCurveSurface.cpp
GEntity.cpp STensor3.cpp
GVertex.cpp GEdge.cpp GFace.cpp GRegion.cpp
GEdgeLoop.cpp GEdgeCompound.cpp GFaceCompound.cpp
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index cc4b6e3a65..78b3d77db3 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -11,6 +11,7 @@
#include "GmshDefines.h"
#include "GFaceCompound.h"
#include "GEdgeCompound.h"
+#include "intersectCurveSurface.h"
#if defined(HAVE_SOLVER)
@@ -2495,5 +2496,105 @@ void GFaceCompound::printStuff(int iNewton) const
// fclose(uvm);
}
+// useful for mesh generators ----------------------------------------
+// Intersection of a circle and a plane
+
+GPoint GFaceCompound::intersectionWithCircle (const SVector3 &n1, const SVector3 &n2, const SVector3 &p, const double &d, double uv[2]) const
+{
+ SVector3 n = crossprod(n1,n2);
+ n.normalize();
+ for (int i=-1;i<nbT;i++){
+ GFaceCompoundTriangle *ct;
+ double U,V;
+ if (i == -1) getTriangle(uv[0],uv[1], &ct, U,V);
+ else ct = &_gfct[i];
+ if (!ct) continue;
+ // we have two planes, defined with n1,n2 and t1,t2
+ // we have then a direction m = (n1 x n2) x (t1 x t2)
+ // and a point p that defines a straight line
+ // the point is situated in the half plane defined
+ // by direction n1 and point p (exclude the one on the
+ // other side)
+
+ SVector3 t1 = ct->v2 - ct->v1;
+ SVector3 t2 = ct->v3 - ct->v1;
+
+ // let us first compute point q to be the intersection
+ // of the plane of the triangle with the line L = p + t n1
+ // Compute w = t1 x t2 = (a,b,c) and write a point of the plabe as
+ // X(x,y,z) = ax + by + cz - (v1_x a + v1_y b + v1_z * c) = 0
+ // Then
+ // a (p_x + t n1_x) + a (p_y + t n1_y) + c (p_z + t n1_z) - (v1_x a + v1_y b + v1_z * c) = 0
+ // gives t
+
+ SVector3 w = crossprod(t1,t2);
+ double t = dot(ct->v1-p,w)/dot(n1,w);
+ SVector3 q = p + n1*t;
+
+ // consider the line that intersects both planes in its
+ // parametric form
+ // X(x,y,z) = q + t m
+ // We have now to intersect that line with the sphere
+ // (x-p_x)^2 + (y-p_y)^2 + (z-p_z)^2 = d^2
+ // Inserting the parametric form of the line into the sphere gives
+ // (t m_x + q_x - p_x)^2 + (t m_y + q_y - p_y)^2 + (t m_z + q_z - p_z)^2 = d^2
+ // t^2 (m_x^2 + m_y^2 + m_z^2) + 2 t (m_x (q_x - p_x) + m_y (q_y - p_z) + m_z (q_z - p_z)) + ((q_x - p_x)^2 + (q_y - p_y)^2 + (q_z - p_z)^2 - d^2) = 0
+ // t^2 m^2 + 2 t (m . (q-p)) + ((q-p)^2 - d^2) = 0
+ // Let us call ta and tb the two roots
+ // they correspond to two points s_1 = q + ta m and s_2 = q + tb m
+ // we choose the one that corresponds to (s_i - p) . n1 > 0
+ SVector3 m = crossprod(w,n);
+ const double a = dot(m,m);
+ const double b = 2*dot(m,q-p);
+ const double c = dot(q-p,q-p) - d*d;
+ const double delta = b*b-4*a*c;
+
+ if (delta >= 0){
+ const double ta = (-b + sqrt(delta)) / (2.*a);
+ const double tb = (-b - sqrt(delta)) / (2.*a);
+ SVector3 s1 = q + m * ta;
+ SVector3 s2 = q + m * tb;
+ SVector3 s;
+ if (dot(s1-p,n1) > 0){
+ s = s1;
+ }
+ else if (dot(s2-p,n1) > 0){
+ s = s2;
+ }
+ else continue;
+
+ // we have now to look if the point is inside the triangle
+ // s = v1 + u t1 + v t2
+ // we know the system has a solution because s is in the plane
+ // defined by v1 and v2 we solve then
+ // (s - v1) . t1 = u t1^2 + v t2 . t1
+ // (s - v1) . t2 = u t1 . t2 + v t2^2
+
+ double mat[2][2], b[2],uv[2];
+ mat[0][0] = dot(t1,t1);
+ mat[1][1] = dot(t2,t2);
+ mat[0][1] = mat[1][0] = dot(t1,t2);
+ b[0] = dot(s-ct->v1,t1) ;
+ b[1] = dot(s-ct->v1,t2) ;
+ sys2x2(mat,b,uv);
+ // check now if the point is inside the triangle
+ if (uv[0] >= -1.e-6 && uv[1] >= -1.e-6 &&
+ 1.-uv[0]-uv[1] >= -1.e-6 ) {
+ SVector3 pp =
+ ct->p1 * ( 1.-uv[0]-uv[1] ) +
+ ct->p2 *uv[0] +
+ ct->p3 *uv[1] ;
+ // GPoint ttt = point(pp.x(),pp.y());
+ // printf("%g %g %g vs %g %g %g\n",ttt.x(),ttt.y(),ttt.z(),s.x(),s.y(),s.z());
+ // printf("%d/%d\n",i,nbT);
+ return GPoint(s.x(),s.y(),s.z(),this,pp);
+ }
+ }
+ }
+ GPoint pp(0);
+ pp.setNoSuccess();
+ return pp;
+}
+
#endif
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 1424595bc0..3cd63b4d13 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -142,6 +142,8 @@ class GFaceCompound : public GFace {
int getNbSplit() const { return nbSplit; }
int allowPartition() const{ return _allowPartition; }
void setType(typeOfIsomorphism type){ _type=type;}
+ // useful for mesh generators ----------------------------------------
+ GPoint intersectionWithCircle (const SVector3 &n1, const SVector3 &n2, const SVector3 &p, const double &d, double uv[2]) const;
private:
mutable typeOfCompound _toc;
mutable typeOfMapping _mapping;
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 81d880f5e2..9118de87c7 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -7,7 +7,7 @@ set(SRC
CenterlineField.cpp
Generator.cpp
Field.cpp
- meshGEdge.cpp
+ meshGEdge.cpp
meshGEdgeExtruded.cpp
meshGFace.cpp
meshGFaceTransfinite.cpp meshGFaceExtruded.cpp
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index d1a9771d51..e39f254da1 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -521,7 +521,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
//avoid convergent errors
if (dv2.length() < 0.5 * dv.length())break;
- blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
+ blQuads.push_back(new MQuadrangle(v11,v21,v22,v12));
fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v11->x(),v11->y(),v11->z(),
v12->x(),v12->y(),v12->z(),
@@ -1141,7 +1141,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
!CTX::instance()->mesh.optimizeLloyd && !onlyInitialMesh)
- recombineIntoQuadsIterative(gf);
+ recombineIntoQuads(gf);
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 793f8fa035..a2d200abf4 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -20,6 +20,8 @@
#include "MQuadrangle.h"
#include "Field.h"
#include "GModel.h"
+#include "GFaceCompound.h"
+#include "intersectCurveSurface.h"
double LIMIT_ = 0.5 * sqrt(2.0) * 1;
@@ -269,8 +271,6 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double
double pc[3] =
{base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
- // compute the circumradius of the triangle
-
if (!metric){
if (radiusNorm == 2){
circumCenterXYZ(pa, pb, pc, center);
@@ -1072,30 +1072,6 @@ double optimalPointFrontal (GFace *gf,
and the edge length
*/
-struct intersectCurveSurfaceData
-{
- GFace *gf;
- SVector3 &n1,&n2;
- SVector3 &middle;
- double d;
- intersectCurveSurfaceData (GFace *_gf, SVector3 & _n1, SVector3 & _n2, SVector3 & _middle, double &_d)
- : gf(_gf),n1(_n1),n2(_n2),middle(_middle),d(_d)
- { }
-};
-
-static void intersectCircleSurface(fullVector<double> &uvt,
- fullVector<double> &res, void *_data){
- intersectCurveSurfaceData *data = (intersectCurveSurfaceData*)_data;
- GPoint gp = data->gf->point(uvt(0),uvt(1));
- SVector3 dir = (data->n1*cos(uvt(2))+data->n2*sin(uvt(2)))*data->d;
- SVector3 pp = data->middle + dir;
- res(0) = gp.x() - pp.x();
- res(1) = gp.y() - pp.y();
- res(2) = gp.z() - pp.z();
- // printf("f(%g %g %g) = %12.5E %12.5E %12.5E \n",uvt(0),uvt(1),uvt(2),res(0),res(1),res(2));
-}
-
-
void optimalPointFrontalB (GFace *gf,
MTri3* worst,
int active_edge,
@@ -1105,6 +1081,7 @@ void optimalPointFrontalB (GFace *gf,
std::vector<double> &vSizesBGM,
double newPoint[2],
double metric[3]){
+ static int missed = 1;
// as a starting point, let us use the "fast algo"
double d = optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
@@ -1127,19 +1104,31 @@ void optimalPointFrontalB (GFace *gf,
// we look for a point that is
// P = d * (n1 cos(t) + n2 sin(t)) that is on the surface
// so we have to find t, starting with t = 0
- fullVector<double> uvt(3);
- uvt(0) = newPoint[0];
- uvt(1) = newPoint[1];
- uvt(2) = 0.0;
- intersectCurveSurfaceData data (gf,n1,n2,middle,d);
- //printf("----------------------------\n");
- if(newton_fd(intersectCircleSurface, uvt, &data)){
- //printf("--- CONVERGED -- %12.5E -----------\n",d);
- newPoint[0] = uvt(0);
- newPoint[1] = uvt(1);
- return;
+
+ if (gf->geomType() == GEntity::CompoundSurface){
+ GFaceCompound *gfc = dynamic_cast<GFaceCompound*> (gf);
+ if (gfc){
+ GPoint gp = gfc->intersectionWithCircle(n1,n2,middle,d,newPoint);
+ if (gp.succeeded()){
+ newPoint[0] = gp.u();
+ newPoint[1] = gp.v();
+ return ;
+ }
+ }
+ }
+
+ double uvt[3] = {newPoint[0],newPoint[1],0.0};
+ curveFunctorCircle cc (n1,n2,middle,d);
+ surfaceFunctorGFace ss (gf);
+
+ if (intersectCurveSurface (cc,ss,uvt,d*1.e-8)){
+ newPoint[0] = uvt[0];
+ newPoint[1] = uvt[1];
}
- //printf("--- NOT CONVERGED ----------\n");
+ else {
+ Msg::Debug("--- Non optimal point found -----------");
+ }
+ // fclose(f);
}
void bowyerWatsonFrontal(GFace *gf)
@@ -1168,13 +1157,13 @@ void bowyerWatsonFrontal(GFace *gf)
// insert points
while (1){
- /*
- if(ITER % 1== 0){
+
+ /* if(ITER % 100== 0){
char name[245];
sprintf(name,"delfr2d%d-ITER%d.pos",gf->tag(),ITER);
- _printTris (name, AllTris, Us,Vs,false);
- sprintf(name,"delfr2dA%d-ITER%d.pos",gf->tag(),ITER);
- _printTris (name, ActiveTris, Us,Vs,false);
+ _printTris (name, AllTris, Us,Vs,true);
+ // sprintf(name,"delfr2dA%d-ITER%d.pos",gf->tag(),ITER);
+ // _printTris (name, ActiveTris, Us,Vs,false);
}
*/
if (!ActiveTris.size())break;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 0148e46973..1a70c55c42 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -860,6 +860,7 @@ void createRegularMesh (GFace *gf,
// create quads
for (int i=0;i<=N;i++){
for (int j=0;j<=M;j++){
+ //MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j+1][i],tab[j+1][i+1],tab[j][i+1]);
MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j][i+1],tab[j+1][i+1],tab[j+1][i]);
// printf("%p %p %p %p\n",tab[j][i],tab[j][i+1],tab[j+1][i+1],tab[j+1][i]);
q.push_back(qnew);
@@ -1011,7 +1012,11 @@ struct quadBlob {
int topo_convex_region = 1;
for (int i=0; i<bnodes.size();i++){
MVertex *v = bnodes[i];
- // if (v->onWhat()->dim() != 2)return 0;
+ // do not do anything in boundary layers !!!
+ if (v->onWhat()->dim() == 2){
+ MFaceVertex *fv = dynamic_cast<MFaceVertex*>(v);
+ if(fv && fv->bl_data) return 0;
+ }
//if (v->onWhat()->dim() == 2){
int topo_angle = topologicalAngle(v);
if (topo_angle < 0){
@@ -3006,29 +3011,6 @@ void quadsToTriangles(GFace *gf, double minqual)
gf->quadrangles = qds;
}
-void recombineIntoQuadsIterative(GFace *gf)
-{
- recombineIntoQuads(gf);
- // quadsToTriangles(gf, 1. - gf->meshAttributes.recombineAngle/90. );
- return;
- int COUNT = 0;
- while (1){
- quadsToTriangles(gf,-10000);
- {char NAME[245];sprintf(NAME,"iterT%d.msh",COUNT); gf->model()->writeMSH(NAME);}
- std::set<MTri3*,compareTri3Ptr> AllTris;
- std::vector<double> vSizes, vSizesBGM, Us, Vs;
- std::vector<SMetric3> vMetricsBGM;
- // buildMeshGenerationDataStructures
- // (gf, AllTris, vSizes, vSizesBGM, vMetricsBGM, Us, Vs);
- // int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
- // transferDataStructure(gf, AllTris, Us, Vs);
- {char NAME[245];sprintf(NAME,"iterTD%d.msh",COUNT); gf->model()->writeMSH(NAME);}
- recombineIntoQuads(gf,false,true);
- {char NAME[245];sprintf(NAME,"iter%d.msh",COUNT++); gf->model()->writeMSH(NAME);}
- if (COUNT == 5)break;
- }
-}
-
double Temporary::w1,Temporary::w2,Temporary::w3;
std::vector<SVector3> Temporary::gradients;
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 289f0b3351..2a44576ea5 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -97,7 +97,6 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
void recombineIntoQuads(GFace *gf,
bool topologicalOpti = true,
bool nodeRepositioning = true);
-void recombineIntoQuadsIterative(GFace *gf);
void quadsToTriangles(GFace *gf, double minqual);
struct swapquad{
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 1e5d1ef91a..1867dec2c3 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -680,7 +680,8 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
if (N == 1)
dx(0) = f(0) / J(0, 0);
else
- J.luSolve(f, dx);
+ if (!J.luSolve(f, dx))
+ return false;
for (int i = 0; i < N; i++)
x(i) -= relax * dx(i);
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 4b716034d4..52314e78d1 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -285,10 +285,10 @@ bool fullMatrix<double>::luSolve(const fullVector<double> &rhs, fullVector<doubl
F77NAME(dgesv)(&N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
delete [] ipiv;
if(info == 0) return true;
- if(info > 0)
- Msg::Debug("U(%d,%d)=0 in LU decomposition", info, info);
- else
- Msg::Debug("Wrong %d-th argument in LU decomposition", -info);
+ // if(info > 0)
+ // Msg::Error("U(%d,%d)=0 in LU decomposition", info, info);
+ // else
+ // Msg::Error("Wrong %d-th argument in LU decomposition", -info);
return false;
}
diff --git a/benchmarks/stl/artery.geo b/benchmarks/stl/artery.geo
index 2eabc40ad2..2f9d95a1c4 100644
--- a/benchmarks/stl/artery.geo
+++ b/benchmarks/stl/artery.geo
@@ -1,5 +1,5 @@
-Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
-Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
+Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
+Mesh.RemeshAlgorithm=2; //(0) nosplit (1) automatic (2) split metis
Mesh.CharacteristicLengthFactor=0.05;
diff --git a/benchmarks/stl/skull.geo b/benchmarks/stl/skull.geo
index ff18adfc99..3dee977428 100644
--- a/benchmarks/stl/skull.geo
+++ b/benchmarks/stl/skull.geo
@@ -1,8 +1,8 @@
Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
Mesh.RemeshAlgorithm=2; //(0) nosplit (1) automatic (2) split metis
-Mesh.CharacteristicLengthFactor=0.3;
-Mesh.Algorithm3D = 4; //Frontal (4) Delaunay(1)
+Mesh.CharacteristicLengthFactor=0.2;
+//Mesh.Algorithm3D = 4; //Frontal (4) Delaunay(1)
Merge "skullU.stl";
CreateTopology;
--
GitLab