From d3e7a3e539e8ea2f40a2c9f245ad93742e3238b4 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 8 Feb 2012 09:26:23 +0000
Subject: [PATCH] added a BL test case in validation suite; smoothing 1D mesh
is now (maybe) correct
---
Geo/GeoInterpolation.cpp | 2 +-
Mesh/BackgroundMesh.cpp | 3 +-
Mesh/Field.cpp | 2 +-
Mesh/HighOrder.cpp | 2 +-
Mesh/meshGEdge.cpp | 50 +++++++++++++++++------------
Mesh/meshGFaceBoundaryLayers.cpp | 6 ++--
Mesh/meshGFaceDelaunayInsertion.cpp | 3 +-
7 files changed, 40 insertions(+), 28 deletions(-)
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index ed1e23b88a..45ba40fb94 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -736,7 +736,7 @@ static Vertex InterpolateExtrudedSurface(Surface *s, double u, double v)
Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v)
{
if(derivee == 1) {
- double eps = 1.e-6;
+ double eps = 1.e-8;
Vertex D[4];
if(u_v == 1) {
if(u - eps < 0.0) {
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index e0c01c4ba1..fc8fffdf66 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -174,7 +174,7 @@ static SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, dou
static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
{
- SMetric3 mesh_size(1.e-05);
+ SMetric3 mesh_size(1.e-12);
std::list<GFace *> faces = ge->faces();
std::list<GFace *>::iterator it = faces.begin();
int count = 0;
@@ -192,6 +192,7 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
double Crv = ge->curvature(u);
double lambda = ((2 * M_PI) /( fabs(Crv) * CTX::instance()->mesh.minCircPoints ) );
SMetric3 curvMetric (1./(lambda*lambda));
+
return intersection(mesh_size,curvMetric);
}
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 192e9b1cb5..c7a5d197b2 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1281,7 +1281,7 @@ class MinAnisoField : public Field
else{
(*f) (x, y, z, ff, ge);
}
- v = intersection(v,ff);
+ v = intersection_conserve_mostaniso(v,ff);
}
}
metr = v;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index c7e17705f0..1a9049d7be 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -122,7 +122,7 @@ static bool computeEquidistantParameters0(GEdge *ge, double u0, double uN, int N
return false;
}
// 1 = geodesics
-static int method_for_computing_intermediary_points = 0;
+static int method_for_computing_intermediary_points = 1;
static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
double *u, double underRelax){
if (method_for_computing_intermediary_points == 0) // use linear abscissa
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 73c2a71ff8..bb79388edd 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -31,11 +31,8 @@ static double smoothPrimitive (GEdge *ge,
double alpha,
std::vector<IntPoint> &Points)
{
- for (int i=0; i< Points.size(); i++){
- IntPoint &pt = Points[i];
- SVector3 der = ge->firstDer(pt.t);
- pt.xp = der.norm();
- }
+ // printf("alpha = %g\n",alpha);
+
int ITER = 0;
while (1){
int count1 = 0;
@@ -45,28 +42,34 @@ static double smoothPrimitive (GEdge *ge,
// iterate forward and then backward
// convergence is usually very fast.
for (int i=1; i< Points.size(); i++){
- double hprime;
- hprime = ((Points[i].xp/Points[i].lc) - (Points[i-1].xp/Points[i-1].lc))/(Points[i].t - Points[i-1].t);
- if (hprime * Points[i].xp > alpha*1.01){
- double hnew = Points[i-1].xp / Points[i-1].lc + (Points[i].t - Points[i-1].t) * alpha / Points[i].xp;
+ double dh = (Points[i].xp/Points[i].lc - Points[i-1].xp/Points[i-1].lc);
+ double dt = Points[i].t - Points[i-1].t;
+ double dhdt = dh/dt;
+ if (dhdt / Points[i].xp > (alpha - 1.)*1.01){
+ double hnew = Points[i-1].xp / Points[i-1].lc + dt * (alpha-1) * Points[i].xp;
Points[i].lc = Points[i].xp / hnew;
count1 ++;
}
}
-
+
for (int i=Points.size()-1; i>0 ; i--){
- double hprime;
- hprime = ((Points[i].xp/Points[i].lc) - (Points[i-1].xp/Points[i-1].lc))/(Points[i].t - Points[i-1].t);
- if (-hprime * Points[i].xp > alpha*1.01){
- double hnew = Points[i].xp / Points[i].lc + (Points[i].t - Points[i-1].t) * alpha / Points[i-1].xp;
+ double dh = (Points[i-1].xp/Points[i-1].lc - Points[i].xp/Points[i].lc);
+ double dt = fabs(Points[i-1].t - Points[i].t);
+ double dhdt = dh/dt;
+ if (dhdt / Points[i-1].xp > (alpha-1.)*1.01){
+ double hnew = Points[i].xp / Points[i].lc + dt * (alpha-1) * Points[i].xp;
Points[i-1].lc = Points[i-1].xp / hnew;
count2 ++;
+ // dh = (Points[i-1].xp/Points[i-1].lc - Points[i].xp/Points[i].lc);
+ // dt = fabs(Points[i-1].t - Points[i].t);
+ // dhdt = dh/dt;
+ // printf("dhdt / Points[i-1].xp
}
}
-
+
++ITER;
- if (ITER > 1000){printf("OUUUUUH\n");break;}
+ if (ITER > 2000){break;}
if (!(count1+count2))break;
}
// recompute the primitive
@@ -226,7 +229,7 @@ static void RecursiveIntegration(GEdge *ge, IntPoint *from, IntPoint *to,
double val3 = trapezoidal(&P, to);
double err = fabs(val1 - val2 - val3);
- if(((err < Prec) && (*depth > 1)) || (*depth > 25)) {
+ if(((err < Prec) && (*depth > 3)) || (*depth > 25)) {
p1=Points.back();
P.p = p1.p + val2;
Points.push_back(P);
@@ -300,9 +303,11 @@ static void printFandPrimitive(int tag, std::vector<IntPoint> &Points)
sprintf(name, "line%d.dat", tag);
FILE *f = fopen(name, "w");
if(!f) return;
+ double l = 0;
for (unsigned int i = 0; i < Points.size(); i++){
const IntPoint &P = Points[i];
- fprintf(f, "%g %g %g %g\n", P.t, 1./P.lc, P.p,P.lc);
+ if (i) l +=(P.t - Points[i-1].t)*P.xp;
+ fprintf(f, "%g %g %g %g %g\n", P.t, P.xp/P.lc, P.p,P.lc, l);
}
fclose(f);
}
@@ -313,7 +318,7 @@ void meshGEdge::operator() (GEdge *ge)
FieldManager *fields = ge->model()->getFields();
BoundaryLayerField *blf = 0;
if(fields->getBackgroundField() > 0){
- Field *bl_field = fields->get(fields->getBackgroundField());
+ Field *bl_field = fields->get(fields->getBoundaryLayerField());
blf = dynamic_cast<BoundaryLayerField*> (bl_field);
}
#else
@@ -395,8 +400,13 @@ void meshGEdge::operator() (GEdge *ge)
CTX::instance()->mesh.lcIntegrationPrecision);
// FIXME JF : MAYBE WE SHOULD NOT ALWAYS SMOOTH THE 1D MESH SIZE FIELD ??
+ for (int i=0; i< Points.size(); i++){
+ IntPoint &pt = Points[i];
+ SVector3 der = ge->firstDer(pt.t);
+ pt.xp = der.norm();
+ }
//printFandPrimitive(ge->tag(), Points);
- // a = smoothPrimitive (ge,CTX::instance()->mesh.smoothRatio*CTX::instance()->mesh.smoothRatio,Points);
+ a = smoothPrimitive (ge,/*CTX::instance()->mesh.smoothRatio*/CTX::instance()->mesh.smoothRatio,Points);
// printf(" coucou %g\n",a);
//printFandPrimitive(ge->tag()+10000, Points);
diff --git a/Mesh/meshGFaceBoundaryLayers.cpp b/Mesh/meshGFaceBoundaryLayers.cpp
index 889cd5ffc0..f313d257d3 100644
--- a/Mesh/meshGFaceBoundaryLayers.cpp
+++ b/Mesh/meshGFaceBoundaryLayers.cpp
@@ -241,9 +241,9 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf, double _treshold)
}
else if (angle >= _treshold){
int fanSize = angle / _treshold;
- printf("ONE FAN HAS BEEN CREATED : %d %d %d %d ANGLE = %g | %g %g %g %g\n",e1.getVertex(0)->getNum(),
- e1.getVertex(1)->getNum(),e2.getVertex(0)->getNum(),e2.getVertex(1)->getNum(),
- angle/M_PI*180,N1[SIDE].x(),N1[SIDE].y(),N2[SIDE].x(),N2[SIDE].y());
+ // printf("ONE FAN HAS BEEN CREATED : %d %d %d %d ANGLE = %g | %g %g %g %g\n",e1.getVertex(0)->getNum(),
+ // e1.getVertex(1)->getNum(),e2.getVertex(0)->getNum(),e2.getVertex(1)->getNum(),
+ // angle/M_PI*180,N1[SIDE].x(),N1[SIDE].y(),N2[SIDE].x(),N2[SIDE].y());
// if the angle is greater than PI, than reverse the sense
double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x());
double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x());
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index f42fdde41d..0ffbe0d84f 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1129,7 +1129,8 @@ void bowyerWatsonFrontal(GFace *gf)
// sprintf(name,"frontal%d-param.pos", gf->tag());
// _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
- quadsToTriangles(gf,10000);
+ // in case of boundary layer meshing
+ // quadsToTriangles(gf,10000);
}
void optimalPointFrontalQuad (GFace *gf,
--
GitLab