diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index a1f6f2fa652ae534dd0c23a05827208fb0fc5407..a7a0f011e8b3632c17d93a61342237823cd8ff02 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -43,10 +43,10 @@ SMetric3 buildMetricTangentToCurve(SVector3 &t, double l_t, double l_n)
{
if (l_t == 0.0) return SMetric3(1.e-22);
SVector3 a;
- if (t(0) <= t(1) && t(0) <= t(2)){
+ if (fabs(t(0)) <= fabs(t(1)) && fabs(t(0)) <= fabs(t(2))){
a = SVector3(1,0,0);
}
- else if (t(1) <= t(0) && t(1) <= t(2)){
+ else if (fabs(t(1)) <= fabs(t(0)) && fabs(t(1)) <= fabs(t(2))){
a = SVector3(0,1,0);
}
else{
@@ -58,6 +58,7 @@ SMetric3 buildMetricTangentToCurve(SVector3 &t, double l_t, double l_n)
c.normalize();
t.normalize();
SMetric3 Metric (1./(l_t*l_t),1./(l_n*l_n),1./(l_n*l_n),t,b,c);
+ // printf("bmttc %g %g %g %g %g\n",l_t,l_n,Metric(0,0),Metric(0,1),Metric(1,1));
return Metric;
}
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index efb78c8a9eb975a91309f8ba3f4336b8834dfd68..acd71be9d3998923e3aea2c2ae42c2814f457e78 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1931,7 +1931,7 @@ void BoundaryLayerField::operator() (AttractorField *cc, double dist,
const double b = lc_t;
const double h = lc_n;
double oneOverD2 = .5/(b * b) *
- (1. + sqrt (1. + (4. * crv * crv * b * b * b * b / (h * h * beta * beta))));
+ (1. + sqrt (1. + (4. * crv * crv * b * b * b * b / (h * h * beta * beta))));
metr = buildMetricTangentToCurve(t1, sqrt(1. / oneOverD2), lc_n);
return;
}
@@ -2013,7 +2013,7 @@ void BoundaryLayerField::operator() (double x, double y, double z,
}
if (iIntersect)
for (unsigned int i = 0; i < hop.size(); i++)
- v = intersection_conserve_mostaniso(v, hop[i]);
+ v = intersection_conserveM1(v, hop[i]);
metr = v;
}
#endif
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index cc99a7c6b2c477e87c62c7ac17eb2ce77a1dd542..7ff2757035e1f5f0db0ad6663b21f41fb0e23089 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -369,7 +369,7 @@ void meshGEdge::operator() (GEdge *ge)
N = ge->meshAttributes.nbPointsTransfinite;
}
else{
- if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf)
+ if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS || blf)
a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
CTX::instance()->mesh.lcIntegrationPrecision);
else{
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 947e0a333745e8e540f836ddf6de2f068f4c2e5e..223eb34515c706e83c79fe8ac06be43b76b2f804 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -274,10 +274,10 @@ int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
const double a = metric[0];
const double b = metric[1];
const double d = metric[2];
- double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
- + (x[1] - uv[1]) * (x[1] - uv[1]) * d
- + 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
- return d2 < Radius2;
+ const double d0 = (x[0] - uv[0]);
+ const double d1 = (x[1] - uv[1]);
+ const double d3 = d0*d0*a + d1*d1*d + 2.0*d0*d1*b;
+ return d3 < Radius2;
}
int inCircumCircleAniso(GFace *gf, MTriangle *base,
@@ -292,7 +292,6 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base,
(Vs[base->getVertex(0)->getIndex()] +
Vs[base->getVertex(1)->getIndex()] +
Vs[base->getVertex(2)->getIndex()]) / 3.};
-
buildMetric(gf, pa, metric);
circumCenterMetric(base, metric, Us, Vs, x, Radius2);
@@ -300,11 +299,10 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base,
const double b = metric[1];
const double d = metric[2];
- double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
- + (x[1] - uv[1]) * (x[1] - uv[1]) * d
- + 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
-
- return d2 < Radius2;
+ const double d0 = (x[0] - uv[0]);
+ const double d1 = (x[1] - uv[1]);
+ const double d3 = d0*d0*a + d1*d1*d + 2.0*d0*d1*b;
+ return d3 < Radius2;
}
MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric,
@@ -558,7 +556,7 @@ bool invMapUV(MTriangle *t, double *p,
return false;
}
-double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
+inline double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
{
double u1 = Us[t->getVertex(0)->getIndex()];
double v1 = Vs[t->getVertex(0)->getIndex()];
@@ -572,6 +570,7 @@ double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
return s * 0.5;
}
+double __DT2;
bool insertVertexB (std::list<edgeXface> &shell,
std::list<MTri3*> &cavity,
bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
@@ -585,6 +584,8 @@ bool insertVertexB (std::list<edgeXface> &shell,
double *metric,
MTri3 **oneNewTriangle)
{
+ if (shell.size() <= 3 || shell.size() != cavity.size() + 2) return false;
+
std::list<MTri3*> new_cavity;
// check that volume is conserved
@@ -616,7 +617,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
MTri3 *t4;
t4 = new MTri3(t, LL,0,&Us,&Vs,gf);
- if (oneNewTriangle) *oneNewTriangle = t4;
+ if (oneNewTriangle) {force = true; *oneNewTriangle = t4;}
// double din = t->getInnerRadius();
double d1 = distance(it->v[0],v);
@@ -626,7 +627,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
// avoid angles that are too obtuse
double cosv = ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2));
- if ((d1 < LL * .0025 || d2 < LL * .0025 || cosv < -.9999999) && !force) {
+ if ((d1 < LL * .5 || d2 < LL * .5 || cosv < -.9) && !force) {
onePointIsTooClose = true;
}
@@ -643,15 +644,16 @@ bool insertVertexB (std::list<edgeXface> &shell,
newVolume += ss;
++it;
}
+
-
- if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() > 3 &&
- !onePointIsTooClose){
+ if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && !onePointIsTooClose){
connectTris_vector(new_cavity.begin(), new_cavity.end());
- // clock_t t1 = clock();
+ // printf("%d %d\n",shell.size(),cavity.size());
+ clock_t t1 = clock();
+ // 30 % of the time is spent here !!!
allTets.insert(newTris, newTris + shell.size());
// clock_t t2 = clock();
- // DT_INSERT_VERTEX += (double)(t2-t1)/CLOCKS_PER_SEC;
+ __DT2 += (double)(clock()-t1)/CLOCKS_PER_SEC;
if (activeTets){
for (std::list<MTri3*>::iterator i = new_cavity.begin(); i != new_cavity.end(); ++i){
int active_edge;
@@ -668,7 +670,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
// The cavity is NOT star shaped
else{
- // printf("(%g %g) %22.15E %22.15E %d %d %d %d %d\n",v->x(),v->y(),oldVolume, newVolume, newVolume,shell.size(), onePointIsTooClose, cavity.size(), new_cavity.size(),allTets.size());
+ // printf("(%g %g) %22.15E %22.15E %d %d %d %d %d\n",v->x(),v->y(),oldVolume, newVolume, newVolume,shell.size(), onePointIsTooClose, cavity.size(), new_cavity.size(),allTets.size());
// printf("12.5E 12.5E 12.5E 12.5E 12.5E\n",d1,d2,LL,cosv);
ittet = cavity.begin();
@@ -764,7 +766,7 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2],
}
if (!force)return 0; // FIXME: removing this leads to horrible performance
-
+
N_GLOBAL_SEARCH ++ ;
for(std::set<MTri3*,compareTri3Ptr>::iterator itx = AllTris.begin();
itx != AllTris.end();++itx){
@@ -778,6 +780,8 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2],
return 0;
}
+//double __DT1;
+
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,
@@ -806,7 +810,9 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
// if the point is able to break the bad triangle "worst"
if (1){
if (inCircumCircleAniso(gf, worst->tri(), center, metric, Us, Vs)){
+ // clock_t t1 = clock();
recurFindCavityAniso(gf, shell, cavity, metric, center, worst, Us, Vs);
+ // __DT1 += (double) (clock() - t1)/CLOCKS_PER_SEC ;
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;
@@ -816,14 +822,14 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
// else look for it
else {
- printf("we should never be here !!!\n");
+ // printf("cocuou\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) {
@@ -853,8 +859,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
if(!p.succeeded() || !insertVertexB
(shell, cavity,false, gf, v, center, ptin, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric, oneNewTriangle) ) {
- // Msg::Debug("Point %g %g cannot be inserted because %d",
- // center[0], center[1], p.succeeded() );
+ Msg::Debug("Point %g %g cannot be inserted because %d",
+ center[0], center[1], p.succeeded() );
// printf("Point %g %g cannot be inserted because %d",
// center[0], center[1], p.succeeded() );
AllTris.erase(it);
@@ -910,6 +916,8 @@ void bowyerWatson(GFace *gf, int MAXPNT)
}
int ITER = 0;
+ int NBDELETED = 0;
+ // double DT1 = 0 , DT2=0, DT3=0;
while (1){
// if(ITER % 1== 0){
// char name[245];
@@ -918,14 +926,20 @@ void bowyerWatson(GFace *gf, int MAXPNT)
// }
MTri3 *worst = *AllTris.begin();
if (worst->isDeleted()){
+ // clock_t t1 = clock();
delete worst->tri();
delete worst;
AllTris.erase(AllTris.begin());
+ NBDELETED ++;
+ // DT1 += (double) (clock() - t1)/CLOCKS_PER_SEC ;
}
else{
- if(ITER++ % 5000 == 0)
+ // clock_t t2 = clock();
+ if(ITER++ % 5000 == 0){
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
vSizes.size(), worst->getRadius());
+ printf("%d %d %d\n",vSizes.size(), AllTris.size(),NBDELETED);
+ }
double center[2],metric[3],r2;
if (worst->getRadius() < /*1.333333/(sqrt(3.0))*/0.5 * sqrt(2.0) ||
(int)vSizes.size() > MAXPNT) break;
@@ -939,10 +953,15 @@ void bowyerWatson(GFace *gf, int MAXPNT)
Vs[base->getVertex(2)->getIndex()]) / 3.};
buildMetric(gf, pa, metric);
circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
+ // DT2 += (double) (clock() - t2)/CLOCKS_PER_SEC ;
+ // clock_t t3 = clock() ;
insertAPoint(gf, AllTris.begin(), center, metric, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris);
+ // DT3 += (double) (clock() - t3)/CLOCKS_PER_SEC ;
}
}
+ // printf("%12.5E %12.5E %12.5E %12.5E %12.5E\n",DT1,DT2,DT3,__DT1,__DT2);
+ printf("%12.5E \n",__DT2);
#if defined(HAVE_ANN)
{
FieldManager *fields = gf->model()->getFields();
@@ -1176,6 +1195,7 @@ void optimalPointFrontalB (GFace *gf,
}
else {
Msg::Debug("--- Non optimal point found -----------");
+ // Msg::Info("--- Non optimal point found -----------");
}
}
@@ -1230,7 +1250,7 @@ void bowyerWatsonFrontal(GFace *gf)
worst->getRadius() > LIMIT_){
if(ITER++ % 5000 == 0)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
- vSizes.size(), worst->getRadius());
+ gf->mesh_vertices.size(), worst->getRadius());
double newPoint[2], metric[3];
//optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
optimalPointFrontalB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
@@ -1703,7 +1723,7 @@ void bowyerWatsonParallelograms(GFace *gf)
buildMetric(gf, newPoint, metric);
SMetric3 ANIZO_MESH = metrics[i];
- bool success = insertAPoint(gf, AllTris.begin(), newPoint, metric, Us, Vs, vSizes,
+ bool success = insertAPoint( gf, AllTris.begin(), newPoint, metric, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris, 0, oneNewTriangle, &oneNewTriangle);
if (!success) oneNewTriangle = 0;
// if (!success)printf("success %d %d\n",success,AllTris.size());
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 003341ca3facfff7454404f2ab1793bf1afaf37d..19ddefc7fd7267ffa31799d3e412915de74850db 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -49,7 +49,7 @@ class MTri3
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; }
+ inline double getRadius() const { return circum_radius; }
MTri3(MTriangle *t, double lc, SMetric3 *m = 0, const std::vector<double> *Us = 0, const std::vector<double> *Vs = 0, GFace *gf = 0);
inline MTriangle *tri() const { return base; }
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index 8c5d44206412d03a5b97a59ef068496ee6733c04..c531b3d8aeeea05232d222e9c67031bb0685e6ac 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -17,7 +17,7 @@
#include "BackgroundMesh.h"
#include "intersectCurveSurface.h"
-static const double FACTOR = .81;
+static const double FACTOR = .71;
static const int NUMDIR = 3;
static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.};
@@ -163,6 +163,8 @@ bool compute4neighbors (GFace *gf, // the surface
metricField = SMetric3(1./(L*L));
}
}
+
+ // printf("M = (%g %g %g)\n",metricField(0,0),metricField(1,1),metricField(0,1));
// get the unit normal at that point
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(midpoint[0],midpoint[1]));
@@ -220,7 +222,15 @@ bool compute4neighbors (GFace *gf, // the surface
2*E*covar2[1]*covar2[0]+
N*covar2[1]*covar2[1]);
- // printf("%12.5E %12.5E %g %g %g %g %g %g %g %g %g %g %g\n",l1,l2,size_1,size_2,size_param_1,size_param_2,M,N,E,s1.x(),s1.y(),s2.x(),s2.y());
+ /* printf("%12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %g %g %g %g %g %g %g %g %g %g %g\n",
+ M*covar1[0]*covar1[0]+
+ 2*E*covar1[1]*covar1[0]+
+ N*covar1[1]*covar1[1],
+ M*covar2[0]*covar2[0]+
+ 2*E*covar2[1]*covar2[0]+
+ N*covar2[1]*covar2[1]
+ ,covar1[0],covar1[1],covar2[0],covar2[1],l1,l2,size_1,size_2,size_param_1,size_param_2,M,N,E,s1.x(),s1.y(),s2.x(),s2.y());*/
+
// this is the rectangle in the parameter plane.
double r1 = 0*1.e-8*(double)rand() / RAND_MAX;
double r2 = 0*1.e-8*(double)rand() / RAND_MAX;
@@ -270,6 +280,8 @@ bool compute4neighbors (GFace *gf, // the surface
void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics){
#if defined(HAVE_RTREE)
+ // FILE *f = fopen ("parallelograms.pos","w");
+
// get all the boundary vertices
std::set<MVertex*> bnd_vertices;
for(unsigned int i=0;i<gf->getNumMeshElements();i++){
@@ -296,7 +308,9 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
vertices.push_back(sp);
double _min[2],_max[2];
sp->minmax(_min,_max);
- rtree.Insert(_min,_max,sp);
+ // printf("%g %g .. %g %g\n",_min[0],_min[1],_max[0],_max[1]);
+ rtree.Insert(_min,_max,sp);
+ // sp->print(f);
}
// printf("initially : %d vertices in the domain\n",vertices.size());
@@ -345,8 +359,8 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
metrics.push_back(vertices[i]->_meshMetric);
SPoint2 midpoint;
reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint);
- // fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",vertices[i]._v->x(),vertices[i]._v->y(),vertices[i]._v->z());
- fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",midpoint.x(),midpoint.y(),0.0);
+ fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",vertices[i]->_v->x(),vertices[i]->_v->y(),vertices[i]->_v->z());
+ //fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",midpoint.x(),midpoint.y(),0.0);
}
delete vertices[i];
}
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index 90867af8e248e1ab0e7a97fbc1966b938bb0ca9e..dafab665c4533c3d17554c149ec7df4e358488c5 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -235,7 +235,8 @@ Field[2].hfar = 1.5;
Field[2].hwall_n = 0.001;
Field[2].hwall_t = 0.01;
Field[2].ratio = 1.3;
-Field[2].thickness = .1;
+Field[2].thickness = .03;
+Field[2].IntersectMetrics = 1;
//Background Field = 2;
Field[1] = Box;