From 3544d2082ac741e1e90f162e97779de48fe224fa Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 14 Jun 2016 14:05:50 +0000
Subject: [PATCH] 2D boundary layers now are ablt take into account axis of
symmetry
---
Mesh/Field.cpp | 159 ++++++++++++++++++++++++++++-----------------
Mesh/Field.h | 1 +
Mesh/meshGEdge.cpp | 42 ++++++++++--
3 files changed, 138 insertions(+), 64 deletions(-)
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 3707d5d710..6f32a72df4 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -2249,70 +2249,105 @@ void BoundaryLayerField::removeAttractors()
update_needed = true;
}
+void BoundaryLayerField::setupFor1d(int iE) {
+
+ if (faces_id_saved.empty() &&
+ edges_id_saved.empty() &&
+ faces_id_saved.empty() ){
+ faces_id_saved = faces_id;
+ edges_id_saved = edges_id;
+ nodes_id_saved = nodes_id;
+ }
+
+ nodes_id.clear();
+ edges_id.clear();
+ faces_id.clear();
+
+ bool found = std::find(edges_id_saved.begin(), edges_id_saved.end(), iE) !=
+ edges_id_saved.end();
+
+ if (!found) {
+ GEdge *ge = GModel::current()->getEdgeByTag(iE);
+ GVertex *gv0 = ge->getBeginVertex();
+ found = std::find(nodes_id_saved.begin(), nodes_id_saved.end(), gv0->tag()) !=
+ nodes_id_saved.end();
+ if (found)nodes_id.push_back(gv0->tag());
+ GVertex *gv1 = ge->getEndVertex();
+ found = std::find(nodes_id_saved.begin(), nodes_id_saved.end(), gv1->tag()) !=
+ nodes_id_saved.end();
+ if (found)nodes_id.push_back(gv1->tag());
+ }
+ // printf("edge %d %d nodes added\n",iE,nodes_id.size());
+ // getchar();
+ removeAttractors();
+}
+
+
void BoundaryLayerField::setupFor2d(int iF)
{
- if (1 || faces_id.size()){
- /* preprocess data in the following way
- remove GFaces from the attarctors (only used in 2D)
- for edges and vertices
- */
- if ( !faces_id_saved.size()){
- faces_id_saved = faces_id;
- edges_id_saved = edges_id;
- nodes_id_saved = nodes_id;
- }
- nodes_id.clear();
- edges_id.clear();
- faces_id.clear();
-
- // printf("have %d %d\n",faces_id_saved.size(),edges_id_saved.size());
-
- /// FIXME :
- /// NOT REALLY A NICE WAY TO DO IT (VERY AD HOC)
- /// THIS COULD BE PART OF THE INPUT
- /// OR (better) CHANGE THE PHILOSOPHY
-
- GFace *gf = GModel::current()->getFaceByTag(iF);
- std::list<GEdge*> ed = gf->edges();
- // printf("face %d has %d edges\n",iF,ed.size());
- for (std::list<GEdge*>::iterator it = ed.begin();
+ /* preprocess data in the following way
+ remove GFaces from the attarctors (only used in 2D)
+ for edges and vertices
+ */
+ if (faces_id_saved.empty() &&
+ edges_id_saved.empty() &&
+ faces_id_saved.empty() ){
+ faces_id_saved = faces_id;
+ edges_id_saved = edges_id;
+ nodes_id_saved = nodes_id;
+ }
+
+ nodes_id.clear();
+ edges_id.clear();
+ faces_id.clear();
+
+ // printf("have %d %d\n",faces_id_saved.size(),edges_id_saved.size());
+
+ /// FIXME :
+ /// NOT REALLY A NICE WAY TO DO IT (VERY AD HOC)
+ /// THIS COULD BE PART OF THE INPUT
+ /// OR (better) CHANGE THE PHILOSOPHY
+
+ GFace *gf = GModel::current()->getFaceByTag(iF);
+ std::list<GEdge*> ed = gf->edges();
+ // printf("face %d has %d edges\n",iF,ed.size());
+ for (std::list<GEdge*>::iterator it = ed.begin();
it != ed.end() ; ++it){
- bool isIn = false;
- int iE = (*it)->tag();
- bool found = std::find(edges_id_saved.begin(), edges_id_saved.end(), iE) !=
- edges_id_saved.end();
- // printf("edges %d found %d\n",iE,found);
- // this edge is a BL Edge
- if (found){
- std::list<GFace*> fc = (*it)->faces();
- // one only face --> 2D --> BL
- if (fc.size() == 1) isIn = true;
+ bool isIn = false;
+ int iE = (*it)->tag();
+ bool found = std::find(edges_id_saved.begin(), edges_id_saved.end(), iE) !=
+ edges_id_saved.end();
+ // printf("edges %d found %d\n",iE,found);
+ // this edge is a BL Edge
+ if (found){
+ std::list<GFace*> fc = (*it)->faces();
+ // one only face --> 2D --> BL
+ if (fc.size() == 1) isIn = true;
+ else {
+ // more than one face and
+ std::list<GFace*>::iterator itf = fc.begin();
+ bool found_this = std::find(faces_id_saved.begin(), faces_id_saved.end(), iF) !=
+ faces_id_saved.end();
+ if (!found_this)isIn = true;
else {
- // more than one face and
- std::list<GFace*>::iterator itf = fc.begin();
- bool found_this = std::find(faces_id_saved.begin(), faces_id_saved.end(), iF) !=
- faces_id_saved.end();
- if (!found_this)isIn = true;
- else {
- bool foundAll = true;
- for ( ; itf != fc.end() ; ++itf){
- int iF2 = (*itf)->tag();
- foundAll &= std::find(faces_id_saved.begin(), faces_id_saved.end(), iF2) !=
- faces_id_saved.end();
- }
- if(foundAll) isIn = true;
+ bool foundAll = true;
+ for ( ; itf != fc.end() ; ++itf){
+ int iF2 = (*itf)->tag();
+ foundAll &= std::find(faces_id_saved.begin(), faces_id_saved.end(), iF2) !=
+ faces_id_saved.end();
}
+ if(foundAll) isIn = true;
}
}
- if (isIn){
- edges_id.push_back(iE);
- nodes_id.push_back ((*it)->getBeginVertex()->tag());
- nodes_id.push_back ((*it)->getEndVertex()->tag());
- }
}
- // printf("face %d %d BL Edges\n", iF, (int)edges_id.size());
- removeAttractors();
+ if (isIn){
+ edges_id.push_back(iE);
+ nodes_id.push_back ((*it)->getBeginVertex()->tag());
+ nodes_id.push_back ((*it)->getEndVertex()->tag());
+ }
}
+ // printf("face %d %d BL Edges\n", iF, (int)edges_id.size());
+ removeAttractors();
}
void BoundaryLayerField::setupFor3d()
@@ -2339,21 +2374,27 @@ double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge
update_needed = false;
}
+
double dist = 1.e22;
- //AttractorField *cc;
+ if (_att_fields.empty())return dist;
+ // AttractorField *cc;
for (std::list<AttractorField*>::iterator it = _att_fields.begin();
it != _att_fields.end(); ++it){
double cdist = (*(*it)) (x, y, z);
if (cdist < dist){
- //cc = *it;
+ // cc = *it;
dist = cdist;
}
}
+
+ if (dist > thickness*ratio) return 1.e22;
current_distance = dist;
// const double dist = (*field) (x, y, z);
// current_distance = dist;
- const double lc = dist*(ratio-1) + hwall_t;
- // double lc = hwall * pow (ratio, dist / hwall);
+ double lc = dist*(ratio-1) + hwall_n;
+
+ // double lc = hwall_n;
+ // double lc = hwall_n * pow (ratio, dist / hwall_t);
return std::min (hfar,lc);
}
diff --git a/Mesh/Field.h b/Mesh/Field.h
index ee4f902e39..bbc1ed23ca 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -182,6 +182,7 @@ class BoundaryLayerField : public Field {
return std::find(nodes_id.begin(),nodes_id.end(),iV) != nodes_id.end();
}
void computeFor1dMesh(double x, double y, double z, SMetric3 &metr);
+ void setupFor1d(int iE);
void setupFor2d(int iF);
void setupFor3d();
void removeAttractors();
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 752793b9b4..89e8da3718 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -78,12 +78,34 @@ static double smoothPrimitive(GEdge *ge, double alpha,
static double F_LcB(GEdge *ge, double t)
{
+ BoundaryLayerField *blf = 0;
+#if defined(HAVE_ANN)
+ FieldManager *fields = ge->model()->getFields();
+ Field *bl_field = fields->get(fields->getBoundaryLayerField());
+ blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+#endif
+
GPoint p = ge->point(t);
- return BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
+ double lc = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
+
+ if (blf){
+ double lc2 = (*blf)( p.x(), p.y(), p.z() , ge);
+ // printf("p %g %g lc %g\n",p.x(),p.y(),lc2);
+ lc = std::min(lc, lc2);
+ }
+
+ return lc;
}
static double F_Lc(GEdge *ge, double t)
{
+ BoundaryLayerField *blf = 0;
+#if defined(HAVE_ANN)
+ FieldManager *fields = ge->model()->getFields();
+ Field *bl_field = fields->get(fields->getBoundaryLayerField());
+ blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+#endif
+
GPoint p = ge->point(t);
double lc_here;
@@ -98,6 +120,12 @@ static double F_Lc(GEdge *ge, double t)
else
lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
+ if (blf){
+ double lc2 = (*blf)( p.x(), p.y(), p.z() , ge);
+ // printf("p %g %g lc %g\n",p.x(),p.y(),lc2);
+ lc_here = std::min(lc_here, lc2);
+ }
+
SVector3 der = ge->firstDer(t);
const double d = norm(der);
@@ -115,6 +143,8 @@ static double F_Lc_aniso(GEdge *ge, double t)
bool blf = false;
#endif
+ printf("coucou\n");
+
GPoint p = ge->point(t);
SMetric3 lc_here;
@@ -397,11 +427,12 @@ static void filterPoints (GEdge*ge) {
void meshGEdge::operator() (GEdge *ge)
{
+ BoundaryLayerField *blf = 0;
#if defined(HAVE_ANN)
FieldManager *fields = ge->model()->getFields();
- BoundaryLayerField *blf = 0;
Field *bl_field = fields->get(fields->getBoundaryLayerField());
blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+ if (blf)blf->setupFor1d(ge->tag());
#else
bool blf = false;
#endif
@@ -476,7 +507,7 @@ void meshGEdge::operator() (GEdge *ge)
N /= CTX::instance()->mesh.lcFactor;
}
else{
- if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf){
+ if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG/* || blf*/){
a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
CTX::instance()->mesh.lcIntegrationPrecision);
}
@@ -491,7 +522,7 @@ void meshGEdge::operator() (GEdge *ge)
SVector3 der = ge->firstDer(pt.t);
pt.xp = der.norm();
}
- a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
+ // a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.99));
}
@@ -517,7 +548,8 @@ void meshGEdge::operator() (GEdge *ge)
}
}
- // printFandPrimitive(ge->tag(),Points);
+
+ //printFandPrimitive(ge->tag(),Points);
// if the curve is periodic and if the begin vertex is identical to
// the end vertex and if this vertex has only one model curve
--
GitLab