diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index 32d5403ac6246dd824dfe346ce8475fad0bc771a..494383ed2d6dbe6a140bdbd00c9d7c02b8eab5ea 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -3,32 +3,26 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
-#include <stack>
+#include "boundaryLayersData.h"
#include "GmshConfig.h"
#include "GModel.h"
-#include "GFace.h"
-#include "MVertex.h"
#include "MLine.h"
#include "MTriangle.h"
-#include "MTetrahedron.h"
-#include "MHexahedron.h"
-#include "MPrism.h"
#include "MEdge.h"
-#include "boundaryLayersData.h"
#include "OS.h"
-#include "BackgroundMesh.h"
#if !defined(HAVE_MESH) || !defined(HAVE_ANN)
BoundaryLayerField* getBLField(GModel *gm){ return 0; }
bool buildAdditionalPoints2D (GFace *gf ) { return false; }
bool buildAdditionalPoints3D (GRegion *gr) { return false; }
-void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]) {}
-faceColumn BoundaryLayerColumns::getColumns(GFace *gf, MVertex *v1, MVertex *v2,
+//void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]) {}
+/*faceColumn BoundaryLayerColumns::getColumns(GFace *gf, MVertex *v1, MVertex *v2,
MVertex *v3, int side)
{
return faceColumn(BoundaryLayerData(),BoundaryLayerData(),BoundaryLayerData());
}
+*/
edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side)
{
return edgeColumn(BoundaryLayerData(),BoundaryLayerData());
@@ -54,18 +48,6 @@ const int FANSIZE__ = 4;
+
*/
-/*
-static double solidAngle (SVector3 &ni, SVector3 &nj,
- SPoint3 &bi, SPoint3 &bj)
-{
- double cosa = dot (ni, nj);
- SVector3 bibj = bj - bi;
- SVector3 sina = crossprod ( ni , nj );
- double a = atan2(sina.norm(),cosa);
- double sign = dot (nj, bibj);
- return sign > 0 ? fabs (a) : -fabs(a);
-}
-*/
SVector3 interiorNormal(SPoint2 p1, SPoint2 p2, SPoint2 p3)
{
@@ -78,124 +60,6 @@ SVector3 interiorNormal(SPoint2 p1, SPoint2 p2, SPoint2 p3)
return n*(-1.);
}
-double computeAngle(GFace *gf, const MEdge &e1, const MEdge &e2,
- SVector3 &n1, SVector3 &n2)
-{
- double cosa = dot(n1,n2);
- SPoint2 p0,p1,p2;
- MVertex *v11 = e1.getVertex(0);
- MVertex *v12 = e1.getVertex(1);
- MVertex *v21 = e2.getVertex(0);
- MVertex *v22 = e2.getVertex(1);
- MVertex *v0,*v1,*v2;
- if (v11 == v21){
- v0 = v12 ; v1 = v11 ; v2 = v22;
- }
- else if (v11 == v22){
- v0 = v12 ; v1 = v11 ; v2 = v21;
- }
- else if (v12 == v21){
- v0 = v11 ; v1 = v12 ; v2 = v22;
- }
- else if (v12 == v22){
- v0 = v11 ; v1 = v12 ; v2 = v21;
- }
- else throw;
-
- reparamMeshEdgeOnFace(v0, v1, gf, p0, p1);
- reparamMeshEdgeOnFace(v0, v2, gf, p0, p2);
-
- SVector3 t1 (p1.x()-p0.x(),p1.y()-p0.y(),0);
- SVector3 t2 (p2.x()-p1.x(),p2.y()-p1.y(),0);
- t1.normalize();
- t2.normalize();
- SVector3 n = crossprod(t1,t2);
-
- double sign = dot(t1,n2);
-
- double a = atan2 (n.z(),cosa);
- a = sign > 0 ? fabs(a) : -fabs(a);
-
- // printf("a = %12.5e cos %12.5E sin %12.5E %g %g vs %g %g\n",
- // a,cosa,n.z(),n1.x(),n1.y(),n2.x(),n2.y());
- return a;
-}
-
-void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3])
-{
- Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
-
- double res[2][2];
-
- double M[2][3] = {{der.first().x(),der.first().y(),der.first().z()},
- {der.second().x(),der.second().y(),der.second().z()}};
-
-
- for (int i=0;i<2;i++){
- for (int l=0;l<2;l++){
- res[i][l] = 0;
- for (int j=0;j<3;j++){
- for (int k=0;k<3;k++){
- res[i][l] += M[i][j]*m(j,k)*M[l][k];
- }
- }
- }
- }
- metric[0] = res[0][0];
- metric[1] = res[1][0];
- metric[2] = res[1][1];
-}
-
-const BoundaryLayerData & BoundaryLayerColumns::getColumn(MVertex *v, MFace f) const
-{
- int N = getNbColumns(v) ;
- if (N == 1) return getColumn(v, 0);
- std::map<MFace, GFace*, Less_Face>::const_iterator it = _inverse_classification.find(f);
- if (it != _inverse_classification.end()) {
- GFace *gf = it->second;
- for (int i=0;i<N;i++){
- const BoundaryLayerData & c = getColumn(v, i);
- if (std::find(c._joint.begin(),c._joint.end(),gf) != c._joint.end())return c;
- }
- }
- static BoundaryLayerData error;
- return error;
-}
-
-
-faceColumn BoundaryLayerColumns::getColumns(GFace *gf, MVertex *v1, MVertex *v2,
- MVertex *v3, int side)
-{
- // printf("%d %d %d for vertex face %d\n",getNbColumns(v1),getNbColumns(v3),getNbColumns(v3),
- // gf->tag());
- int i1=-1,i2=-1,i3=-1;
- for (int i=0;i<getNbColumns(v1);i++){
- const BoundaryLayerData &d1 = getColumn(v1,i);
- if (std::find(d1._joint.begin(),d1._joint.end(),gf) != d1._joint.end()){
- i1 = i;
- // printf("1 Yeah column %d among %d\n",i,d1._joint.size());
- break;
- }
- }
- for (int i=0;i<getNbColumns(v2);i++){
- const BoundaryLayerData &d2 = getColumn(v2,i);
- if (std::find(d2._joint.begin(),d2._joint.end(),gf) != d2._joint.end()){
- i2 = i;
- // printf("2 Yeah column %d among %d\n",i,d2._joint.size());
- break;
- }
- }
- for (int i=0;i<getNbColumns(v3);i++){
- const BoundaryLayerData &d3 = getColumn(v3,i);
- if (std::find(d3._joint.begin(),d3._joint.end(),gf) != d3._joint.end()){
- i3 = i;
- // printf("3 Yeah column %d among %d\n",i,d3._joint.size());
- break;
- }
- }
-
- return faceColumn(getColumn(v1,i1), getColumn(v2,i2), getColumn(v3,i3));
-}
edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side)
{
@@ -312,8 +176,6 @@ static void treat2Connections(GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2,
itm != _columns->_normals.upper_bound(e2); ++itm) N2.push_back(itm->second);
if (N1.size() == N2.size()){
for (unsigned int SIDE = 0; SIDE < N1.size() ; SIDE++){
- // IF THE ANGLE IS GREATER THAN THRESHOLD, ADD DIRECTIONS !!
- // double angle = computeAngle (gf,e1,e2,N1[SIDE],N2[SIDE]);
if (!fan){
SVector3 x = N1[SIDE]*1.01+N2[SIDE];
x.normalize();
@@ -321,7 +183,7 @@ static void treat2Connections(GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2,
}
else if (fan){
- printf("fan \n");
+ // printf("fan \n");
int fanSize = FANSIZE__;
// if the angle is greater than PI, than reverse the sense
@@ -526,31 +388,27 @@ static void addColumnAtTheEndOfTheBL(GEdge *ge,
MVertex * v0 = g0->mesh_vertices[0];
MVertex * v1 = g1->mesh_vertices[0];
std::vector<MVertex*> invert;
- std::vector<SMetric3> _metrics;
+ // std::vector<SMetric3> _metrics;
for(unsigned int i = 0; i < ge->mesh_vertices.size() ; i++){
invert.push_back(ge->mesh_vertices[ge->mesh_vertices.size() - i - 1]);
- _metrics.push_back(SMetric3(1.0));
+ // _metrics.push_back(SMetric3(1.0));
}
SVector3 t (v1->x()-v0->x(), v1->y()-v0->y(),v1->z()-v0->z());
t.normalize();
if (g0 == gv){
- _columns->addColumn(t, v0, ge->mesh_vertices,_metrics);
+ _columns->addColumn(t, v0, ge->mesh_vertices/*,_metrics*/);
}
else if (g1 == gv){
- _columns->addColumn(t*-1.0, v1,invert,_metrics);
+ _columns->addColumn(t*-1.0, v1,invert/*,_metrics*/);
}
}
}
bool buildAdditionalPoints2D(GFace *gf)
{
- // printf("coucou1\n");
-
+
BoundaryLayerColumns *_columns = gf->getColumns();
-
- _columns->_normals.clear();
- _columns->_non_manifold_edges.clear();
- _columns->_data.clear();
+ _columns->clearData();
// GET THE FIELD THAT DEFINES THE DISTANCE FUNCTION
BoundaryLayerField *blf = getBLField (gf->model());
@@ -673,7 +531,6 @@ bool buildAdditionalPoints2D(GFace *gf)
// now create the BL points
for (unsigned int DIR=0;DIR<_dirs.size();DIR++){
- SPoint2 p;
SVector3 n = _dirs[DIR];
// < ------------------------------- > //
@@ -690,51 +547,23 @@ bool buildAdditionalPoints2D(GFace *gf)
// printf( "coucou c'est moi\n");
}
else {
- MVertex *current = *it;
- reparamMeshVertexOnFace(current,gf,p);
-
- int nbCol = 100;
+ MVertex *first = *it;
std::vector<MVertex*> _column;
- std::vector<SMetric3> _metrics;
- // printf("start with point %g %g (%g %g)\n",current->x(),current->y(),p.x(),p.y());
- AttractorField *catt = 0;
- SPoint3 _close;
- //double _current_distance = 0.;
+ SPoint2 par = gf->parFromPoint(SPoint3(first->x(),first->y(),first->z()));
+ double L = blf->hwall_n ;
while(1){
-
- SMetric3 m;
- double metric[3];
- double l;
- (*blf)(current->x(),current->y(), current->z(), m, current->onWhat());
- if (!catt){
- catt = blf->current_closest;
- _close = blf->_closest_point;
- //_current_distance = blf -> current_distance;
- }
- SPoint2 poffset (p.x() + 1.e-12 * n.x(),
- p.y() + 1.e-12 * n.y());
- buildMeshMetric(gf, poffset, m, metric);
- const double l2 = n.x()*n.x()*metric[0] + 2*n.x()*n.y()*metric[1] + n.y()*n.y()*metric[2] ;
- l = 1./sqrt(l2);
- if (l >= blf->hfar){
- break;
- }
- if (blf -> current_distance > blf->thickness) break;
- catt = blf->current_closest;
- _close = blf->_closest_point;
- //_current_distance = blf -> current_distance;
- SPoint2 pnew (p.x() + l * n.x(),
- p.y() + l * n.y());
- GPoint gp = gf->point (pnew);
- MFaceVertex *_current = new MFaceVertex (gp.x(),gp.y(),gp.z(),gf,pnew.x(),pnew.y());
+ // printf("L = %g\n",L);
+ if (L > blf->thickness) break;
+ SPoint2 pnew (par.x() + L * n.x(),
+ par.y() + L * n.y());
+ GPoint pp = gf->point(pnew);
+ MFaceVertex *_current = new MFaceVertex (pp.x(),pp.y(),pp.z(),gf,pnew.x(),pnew.y());
_current->bl_data = new MVertexBoundaryLayerData;
- current = _current;
- _column.push_back(current);
- _metrics.push_back(m);
- if ((int)_column.size() > nbCol) break; // FIXME
- p = pnew;
+ _column.push_back(_current);
+ int ith = _column.size() ;
+ L+= blf->hwall_n * pow (blf->ratio, ith);
}
- _columns->addColumn(n,*it, _column, _metrics);
+ _columns->addColumn(n,*it, _column /*,_metrics*/);
}
}
}
@@ -763,594 +592,7 @@ bool buildAdditionalPoints2D(GFace *gf)
return 1;
}
-static double angle_0_180(SVector3 &n1, SVector3 &n2)
-{
- double cosa = dot(n1,n2)/(n1.norm()*n2.norm());
- if (cosa > 1.) cosa = 1.0;
- if (cosa < -1.) cosa = -1.0;
- return acos(cosa);
-}
-
-static void createBLPointsInDir(GRegion *gr,
- MVertex *current,
- BoundaryLayerField *blf,
- SVector3 & n,
- std::vector<MVertex*> &_column,
- std::vector<SMetric3> &_metrics)
-{
- SVector3 basis (current->x(),current->y(),current->z());
- double H = blf->hwall_n;
- double dist = H;
- while(dist < blf->thickness){
- SVector3 newp = basis + n * H;
- MVertex *_current = new MVertex (newp.x(),newp.y(),newp.z(),gr);
- // gr->mesh_vertices.push_back(_current);
- _column.push_back(_current);
- H *= blf->ratio;
- dist += H;
- basis = newp;
- }
-}
-
-
-static bool createWedgeBetweenTwoFaces(MVertex *myV,
- SVector3 n1, SVector3 n2,
- std::vector<SVector3> &shoot)
-{
- double angle = angle_0_180 (n1,n2);
- int fanSize = FANSIZE__;
- for (int i=-1; i<=fanSize; i++){
-
- double ti = (double)(i+1)/ (fanSize+1);
- double angle_t = ti * angle;
- double cosA = cos (angle_t);
- double cosAlpha = dot(n1,n2);
-
- const double A = (1.- 2.*cosA*cosA) + cosAlpha*cosAlpha - 2 * cosAlpha*(1.-cosA*cosA);
- const double B = -2*(1.-cosA*cosA)*(1-cosAlpha);
- const double C = 1.-cosA*cosA;
- double DELTA = B*B-4*A*C;
- double t = 0.0;
- if (A == 0.0){
- t = -C / B;
- }
- else if (C != 0.0){
- if (DELTA < 0){
- Msg::Error("this should not happen DELTA = %12.5E",DELTA);
- DELTA = 0;
- }
- const double t1 = (-B+sqrt(DELTA))/(2.*A);
- const double t2 = (-B-sqrt(DELTA))/(2.*A);
-
- SVector3 x1 (n1*(1.-t1) + n2 * t2);
- SVector3 x2 (n1*(1.-t2) + n2 * t2);
- double a1 = angle_0_180 (n1,x1);
- double a2 = angle_0_180 (n1,x2);
- if (fabs(a1 - angle_t) < fabs(a2 - angle_t))t = t1;
- else t = t2;
- }
- SVector3 x (n1*(1.-t) + n2 * t);
- x.normalize();
- shoot.push_back(x);
- }
- return true;
-}
-
-
-
-/*
- a "fanned" edge is one that opens a fan of element between the 2 faces
- it connects in the volume.
-
- fan start (end)
- --> at a symmetry plane or at an inflow/outflow bdry. This implies
- that the fan has a trace on this boundary
- --> at a corner where at least 3 incident model edges are "fanned"
-
- fan can also be doing a closed loop
-
- when looking at a model vertex, one can then have
- --> At least 3 fanned edges : this is a corner
- --> Two fanned edges : this is the passage of a fan
- --> A model vertex that has one single incident fanned edge is forbidden
-
- The algo that allows to build the topology of fans
-
- 1) Compute fanned edges : use threshold angle between adjacent faces as
- user parameter
- 2) If one model vertex is connected to one single fanned edge
- -) If the model vertex is in the symmetry plane --> OK
- -) If the number of incident edges is 2 --> add the second one as a fan
- -) If more that 1 other edge is not fanned --> ask the user to choose
- -) Repeat 2) until everything is OK
- 3) Create sets of model faces that are NOT connected by fans and therefore that
- are considered as equivalent
-*/
-
-// return the edges connected to a vertex for this volume
-
-struct fanTopology {
- std::set<GVertex*> corners;
- std::set<GVertex*> endings;
- std::set<GEdge*> fans;
- std::vector< std::set<GFace *> > groups;
- bool isFan (GEdge* ge) const {return fans.find(ge) != fans.end();}
- int getGroupNum (GFace *gf) const
- {
- for (unsigned int i=0; i< groups.size(); i++){
- if (std::find(groups[i].begin(),groups[i].end(),gf) != groups[i].end())return i;
- }
- return -1;
- }
- fanTopology (GRegion * gr, const std::set<GEdge*> &detectedFans, const std::set<GVertex*> &sym, BoundaryLayerField *blf );
-};
-
-static std::list<GEdge*> edges (GVertex* gv, GRegion* gr) {
- std::list<GEdge*> er = gr->edges();
- std::list<GEdge*> ev = gv->edges();
- std::list<GEdge*> e ;
-
- for (std::list<GEdge*>::iterator ite = ev.begin() ; ite !=ev.end() ; ite++){
- if (std::find (er.begin(), er.end(), *ite) != er.end())e.push_back(*ite);
- }
- return e;
-}
-
-static GFace *otherSideOfEdge (GRegion * gr, GFace *gf, GEdge *ge)
-{
- std::list<GFace*> fe = ge->faces();
- std::list<GFace*> fr = gr->faces();
- for (std::list<GFace*>::iterator it = fe.begin(); it != fe.end(); ++it){
- if ((*it) != gf && std::find(fr.begin(), fr.end(), *it) != fr.end())
- return *it;
- }
- // printf("ouille autre cote de face %d par arerte %d dans region %d foireux\n",gf->tag(),ge->tag(),gr->tag());
- return 0;
-}
-fanTopology :: fanTopology (GRegion * gr, const std::set<GEdge*> &detectedFans, const std::set<GVertex*> &sym, BoundaryLayerField *blf )
-{
- fans = detectedFans;
- endings = sym;
- std::list<GVertex*> vs = gr->vertices();
- while (1){
- unsigned int nb = fans.size();
- for (std::list<GVertex*>::iterator it = vs.begin() ; it != vs.end(); it++){
- GVertex *gv = *it;
- std::list<GEdge*> e = edges (gv, gr);
- std::list<GEdge*> fe ;
- for (std::list<GEdge*>::iterator ite = e.begin() ; ite !=e.end() ; ite++)
- if (isFan(*ite))fe.push_back(*ite);
- if (fe.size() == 1){
- if (endings.find(gv) != endings.end()) continue;
- else if (e.size() == 2)fans.insert(e.begin(), e.end());
- else Msg::Fatal("Ensure that it is possible to find closed loops in BL fans !");
- }
- }
- if (nb == fans.size())break;
- }
- for (std::list<GVertex*>::iterator it = vs.begin() ; it != vs.end(); it++){
- GVertex *gv = *it;
- std::list<GEdge*> e = edges (gv, gr);
- std::list<GEdge*> fe ;
- for (std::list<GEdge*>::iterator ite = e.begin() ; ite !=e.end() ; ite++)
- if (isFan(*ite))fe.push_back(*ite);
- // printf("%d fans connected to vert %d\n",fe.size(),gv->tag(), e.size());
- if (fe.size() > 2)corners.insert(gv);
- }
-
- std::list<GFace*> fs = gr->faces();
-
- while (!fs.empty()) {
- GFace *gf = *fs.begin();
- std::set<GFace *> groupOfFaces;
- groupOfFaces.insert(gf);
- std::stack<GFace*> _stack;
- _stack.push(gf);
- // printf("pushing %d\n",gf->tag());
- while(!_stack.empty()){
- gf = _stack.top();
- _stack.pop();
- std::list<GEdge*> ed = gf->edges();
- for (std::list<GEdge*>::iterator ite = ed.begin() ; ite !=ed.end() ; ite++){
- if (!isFan(*ite) && !blf->isEdgeBLSaved((*ite)->tag())){
- GFace *otherSide = otherSideOfEdge (gr,gf,*ite);
- if (otherSide && groupOfFaces.find(otherSide) == groupOfFaces.end()){
- groupOfFaces.insert(otherSide);
- _stack.push(otherSide);
- // printf("pushing %d through %d\n",otherSide->tag(),(*ite)->tag());
- }
- }
- }
- }
- groups.push_back(groupOfFaces);
- std::list<GFace*> fsb;
- for (std::list<GFace*>::iterator it = fs.begin() ; it != fs.end() ; it++){
- if (groupOfFaces.find(*it) == groupOfFaces.end())fsb.push_back(*it);
- }
- fs = fsb;
- }
-}
-
-// This is the tricky part
-// We have to find a vector N that has the following properties
-// V = (x,y,z) / (x^2+y^2+z^2)^{1/2} is a point on the unit sphere
-// the n[i]'s are points on the unit sphere
-// V maximizes min_i (V * n[i])
-// This means I'd like to find point V that is
-
-/*
-static void filterVectors(std::vector<SVector3> &n)
-{
- std::vector<SVector3> temp;
- temp.push_back(n[0]);
- for (unsigned int i = 1 ; i<n.size() ; i++){
- bool found = false;
- for (unsigned int j = 0 ; j<temp.size() ; j++){
- double d = dot(n[i],temp[j]);
- if (d < 0.98)found = true;
- }
- if (found) temp.push_back(n[i]);
- }
- n = temp;
-}
-*/
-
-static SVector3 computeBestNormal(std::vector<SVector3> &n)
-{
- // filterVectors (n);
- SVector3 V;
- if (n.size() == 1)V = n[0];
- else if (n.size() == 2)V = n[0]+n[1];
- else if (n.size() == 3)circumCenterXYZ(n[0],n[1],n[2],V);
- else {
- Msg::Warning("suboptimal choice for exterior normal: %d vectors to average",n.size());
- for (unsigned int i = 0 ; i<n.size() ; i++)V+=n[i];
- }
- V.normalize();
- return V;
-}
-
-
-static int createColumnsBetweenFaces(GRegion *gr,
- MVertex *myV,
- BoundaryLayerField *blf,
- BoundaryLayerColumns *_columns,
- std::set<GFace*> &_gfaces,
- std::multimap<GFace*,MTriangle*> & _faces,
- std::map<MFace,SVector3,Less_Face> &_normals,
- fanTopology &ft)
-{
- SVector3 n[256];
- int count = 0;
- GFace *gfs[256];
-
- // we compute normals per face using the surface mesh (we could try to use the
- // model but I doubt it'd be better, especially if the mesh is coarse
- for( std::set<GFace*> ::iterator it = _gfaces.begin() ; it != _gfaces.end(); ++it){
- for (std::multimap<GFace*,MTriangle*>::iterator itm =
- _faces.lower_bound(*it);
- itm != _faces.upper_bound(*it); ++itm){
-
- SVector3 N = _normals[itm->second->getFace(0)];
- n[count] += N;
- }
- gfs[count] = *it;
- n[count].normalize();
- count ++;
- }
-
- // the topology of the fans is known
- // we look over all faces and look in which
- // group it lies in
-
- std::vector< std::vector<GFace*> > joints;
-
- std::multimap<int, GFace*> lGroup;
- std::map<GFace*, int> inv;
- std::set<int> gs;
- for (int i=0;i<count;i++) {
- inv[gfs[i]] = i;
- int iGroup = ft.getGroupNum (gfs[i]);
- gs.insert(iGroup);
- lGroup.insert(std::make_pair(iGroup,gfs[i]));
- }
-
-
- for (std::set<int>::iterator it = gs.begin(); it != gs.end() ; ++it){
- std::pair<std::multimap<int, GFace*>::iterator,
- std::multimap<int, GFace*>::iterator> range = lGroup.equal_range(*it);
- std::vector<GFace*> joint;
- for (std::multimap<int, GFace*>::iterator itm = range.first ;
- itm != range.second ; itm++)
- joint.push_back(itm->second);
- joints.push_back(joint);
- // SVector3 avg (0,0,0);
- std::vector<SVector3> ns;
-
- for (unsigned int i=0;i<joint.size(); i++){
- ns.push_back( n[inv[joint[i]]] );
- // avg += n[inv[joint[i]]];
- }
- SVector3 avg = computeBestNormal(ns);
-
- std::vector<MVertex*> _column;
- std::vector<SMetric3> _metrics;
- avg.normalize();
- createBLPointsInDir (gr,myV,blf,avg,_column,_metrics);
- _columns->addColumn(avg,myV, _column, _metrics, joint);
- }
-
- // // DEBUG STUFF
- // if (myV->onWhat()->dim() == 0 && myV->onWhat()->tag() == 6){
- // printf("%d %d\n",gs.size(),joints.size());
- // }
-
- // create wedges
- if (joints.size() > 1){
- for (unsigned int I = 0 ; I < joints.size() ; I++){
- const BoundaryLayerData & c0 = _columns->getColumn(myV, I);
- for (unsigned int J = I+1 ; J < joints.size() ; J++){
- const BoundaryLayerData & c1 = _columns->getColumn(myV, J);
- std::vector<SVector3> shoot;
- createWedgeBetweenTwoFaces(myV,c0._n,c1._n,shoot);
- for (unsigned int i=1;i<shoot.size()-1;i++){
- std::vector<MVertex*> _column;
- std::vector<SMetric3> _metrics;
- createBLPointsInDir (gr,myV,blf,shoot[i],_column,_metrics);
- _columns->addColumn(shoot[i] , myV, _column, _metrics);
- }
- }
- }
- }
- // we have a corner : in a first step, only add one single
- // in the average direction
- if (joints.size() > 2){
- }
-
- return joints.size();
-}
-
-static void createColumnsOnSymmetryPlane(MVertex *myV,
- BoundaryLayerColumns *_columns,
- std::set<GFace*> &_allGFaces,
- std::list<GFace*> &faces,
- fanTopology &ft)
-{
- // get all model faces for that vertex that have boundary layer columns attached
- for ( std::list<GFace*>::iterator itf = faces.begin(); itf!= faces.end() ; ++itf){
- BoundaryLayerColumns* _face_columns = (*itf)->getColumns();
- int N = _face_columns->getNbColumns(myV);
- if (N == 1){
- std::vector<GFace*> _joint;
- _joint.insert(_joint.begin(),_allGFaces.begin(),_allGFaces.end());
- const BoundaryLayerData & c = _face_columns->getColumn(myV,0);
- _columns->addColumn(c._n,myV, c._column, c._metrics, _joint);
- }
- else if (N > 1){
- if (_allGFaces.size() != 2){
- Msg::Fatal("cannot solve such a strange stuff in the BL");
- }
- // printf("%d columns\n",N);
- std::set<GFace*>::iterator itff = _allGFaces.begin();
- GFace *g1 = *itff ; ++itff; GFace *g2 = *itff;
- bool sense = true;
- std::vector<GFace*> _joint;
-
- const BoundaryLayerFan *fan = _face_columns->getFan(myV);
-
- if (fan){
- MVertex *v11 = fan->_e1.getVertex(0);
- MVertex *v12 = fan->_e1.getVertex(1);
- std::list<GEdge*> l1 = g1->edges();
- std::list<GEdge*> l2 = g2->edges();
- if (v11 == myV){
- if (v12->onWhat()->dim() == 1){
- GEdge *ge1 = (GEdge*)v12->onWhat();
- if (std::find(l1.begin(),l1.end(),ge1) != l1.end())sense = fan->sense;
- else if (std::find(l2.begin(),l2.end(),ge1) != l2.end())sense = !fan->sense;
- }
- else
- Msg::Error("Cannot choose between directions in a BL (dim = %d)",
- v12->onWhat()->dim());
- }
- else {
- if (v11->onWhat()->dim() == 1){
- GEdge *ge1 = (GEdge*)v11->onWhat();
- if (std::find(l1.begin(),l1.end(),ge1) != l1.end())sense = fan->sense;
- else if (std::find(l2.begin(),l2.end(),ge1) != l2.end())sense = !fan->sense;
- }
- else Msg::Error("Cannot choose between directions in a BL");
- }
- }
- else{
- Msg::Error("No fan on the outgoing BL");
- }
- _joint.push_back(g1);
- const BoundaryLayerData & c0 = _face_columns->getColumn(myV,sense ? 0 : N-1);
- _columns->addColumn(c0._n,myV, c0._column, c0._metrics,_joint);
- _joint.clear();
- _joint.push_back(g2);
- const BoundaryLayerData & cN = _face_columns->getColumn(myV,sense ? N-1 : 0);
- _columns->addColumn(cN._n,myV, cN._column, cN._metrics,_joint);
- if (sense){
- for (int k=1;k<N-1;k++){
- const BoundaryLayerData & c = _face_columns->getColumn(myV,k);
- _columns->addColumn(c._n,myV, c._column, c._metrics);
- }
- }
- else {
- for (int k=N-2;k>0;k--){
- const BoundaryLayerData & c = _face_columns->getColumn(myV,k);
- _columns->addColumn(c._n,myV, c._column, c._metrics);
- }
- }
- }
- }
-}
-
-static bool preprocessVertex (MVertex *v,
- std::map<MTriangle*,GFace*> &_gfaces,
- BoundaryLayerColumns * _columns,
- std::list<GFace*> &faces,
- std::multimap<GFace*,MTriangle*> &_faces,
- std::set<GFace*> &_allGFaces){
- faces = v->onWhat()->faces();
- for (std::multimap<MVertex*,MTriangle*>::iterator itm =
- _columns->_non_manifold_faces.lower_bound(v);
- itm != _columns->_non_manifold_faces.upper_bound(v); ++itm){
- GFace *gf = _gfaces[itm->second];
- _faces.insert(std::make_pair(gf,itm->second));
- _allGFaces.insert(gf);
- }
-
- bool onSymmetryPlane = false;
- if (v->onWhat()->dim() != 2)
- if (faces.size() != _allGFaces.size())
- onSymmetryPlane = true;
- return onSymmetryPlane;
-}
-
-bool buildAdditionalPoints3D(GRegion *gr)
-{
- BoundaryLayerField *blf = getBLField (gr->model());
-
- if (!blf) return false;
-
- blf->setupFor3d();
-
- BoundaryLayerColumns *_columns = gr->getColumns();
-
- std::list<GFace*> faces = gr->faces();
- std::list<GFace*>::iterator itf = faces.begin();
- std::set<MVertex*> _vertices;
- std::map<MFace,SVector3,Less_Face> _normals;
- std::map<MTriangle*,GFace*> _gfaces;
-
- // filter vertices : belong to BL and are classified on FACES
- while(itf != faces.end()){
- if (blf->isFaceBL((*itf)->tag())){
- // printf("FACE %d is a boundary layer face %d triangles\n",(*itf)->tag(),
- // (int)(*itf)->triangles.size());
- for(unsigned int i = 0; i< (*itf)->triangles.size(); i++){
- _gfaces[(*itf)->triangles[i]] = *itf;
- _columns->_inverse_classification [(*itf)->triangles[i]->getFace(0)] = *itf;
- _normals [(*itf)->triangles[i]->getFace(0)] = SVector3(0,0,0);
- for(unsigned int j = 0; j< 3; j++){
- if ((*itf)->triangles[i]->getVertex(j)->onWhat()->dim() != 3){
- _columns->_non_manifold_faces.insert
- (std::make_pair((*itf)->triangles[i]->getVertex(j),(*itf)->triangles[i]));
- _vertices.insert((*itf)->triangles[i]->getVertex(j));
- }
- }
- }
- }
- ++itf;
- }
- // printf("%d vertices \n", (int)_vertices.size());
-
- // assume that the initial mesh has been created i.e. that there exist
- // tetrahedra inside the domain. Tetrahedra are used to define
- // exterior normals
- for (unsigned int i=0;i<gr->tetrahedra.size();i++){
- for (int j=0;j<4;j++){
- MFace f = gr->tetrahedra[i]->getFace(j);
- std::map<MFace,SVector3,Less_Face>::iterator it = _normals.find(f);
- if (it != _normals.end()){
- MVertex *v0 = f.getVertex(0);
- MVertex *v1 = f.getVertex(1);
- MVertex *v2 = f.getVertex(2);
- MVertex *v3 = 0;
- for (int k=0;k<4;k++){
- if (gr->tetrahedra[i]->getVertex(k) != v0 &&
- gr->tetrahedra[i]->getVertex(k) != v1 &&
- gr->tetrahedra[i]->getVertex(k) != v2 ){
- v3 = gr->tetrahedra[i]->getVertex(k);
- }
- }
- SVector3 v01 (v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z());
- SVector3 v02 (v2->x()-v0->x(),v2->y()-v0->y(),v2->z()-v0->z());
- SVector3 v03 (v3->x()-v0->x(),v3->y()-v0->y(),v3->z()-v0->z());
- SVector3 n = crossprod (v01,v02);
- double sign = dot(n,v03);
- n.normalize();
- if (sign > 0)it->second = n;
- else it->second = n*(-1.0);
- if (_columns->_normals3D.find(it->first) != _columns->_normals3D.end())
- printf("aaaaaaaarghhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh\n");
- _columns->_normals3D.insert(std::make_pair(it->first,it->second));
- }
- }
- }
-
- // look for all the wedges
- std::set<GEdge*> fans;
- std::set<GVertex*> fanNodes;
- for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
- std::list<GFace*> faces;
- std::multimap<GFace*,MTriangle*> _faces;
- std::set<GFace*> _allGFaces;
- preprocessVertex (*it, _gfaces, _columns, faces, _faces, _allGFaces);
- if ((*it)->onWhat()->dim() == 0){
- GVertex *gv = static_cast<GVertex*>((*it)->onWhat());
- if (gv && blf->isFanNode(gv->tag())){
- fanNodes.insert(gv);
- }
- }
- else if ((*it)->onWhat()->dim() == 1){
- GEdge *ge = static_cast<GEdge*>((*it)->onWhat());
- if (ge && blf->isFan(ge->tag())){
- fans.insert(ge);
- }
- }
- }
-
- // printf("STARTING ANALYSIS OF TOPOLOGY (%d fan nodes %d fan edges\n",fanNodes.size(),fans.size());
-
- fanTopology ft (gr, fans, fanNodes, blf);
-
- Msg::Info("BL Topology for region %d : %d FANS %d CORNERS %d GROUPS %d ENDINGS",
- gr->tag(),ft.fans.size(),ft.corners.size(),ft.groups.size(),ft.endings.size());
-
- // for all boundry points
- for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
- std::list<GFace*> faces;
- std::multimap<GFace*,MTriangle*> _faces;
- std::set<GFace*> _allGFaces;
- bool onSymmetryPlane = preprocessVertex (*it, _gfaces, _columns, faces, _faces, _allGFaces);
-
- if (onSymmetryPlane){
- // we have a 3D boundary layer that connects a 2D boundary layer
- // this is the tricky case
- createColumnsOnSymmetryPlane(*it,_columns, _allGFaces, faces, ft);
- }
- else {
- // internal columns (inside the volume)
- createColumnsBetweenFaces (gr,*it,blf,_columns,_allGFaces,_faces,_normals,ft);
- }
- }
-
- // DEBUG STUFF
- FILE *f = Fopen("POINTS3D.pos","w");
- if(f){
- fprintf(f,"View \"\" {\n");
- for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
- MVertex *v = *it;
- for (int i=0;i<_columns->getNbColumns(v);i++){
- const BoundaryLayerData &data = _columns->getColumn(v,i);
- for (unsigned int j=0;j<data._column.size();j++){
- MVertex *blv = data._column[j];
- fprintf(f,"SP(%g,%g,%g){%d};\n",blv->x(),blv->y(),blv->z(),v->getNum());
- }
- }
- }
- fprintf(f,"};\n");
- fclose (f);
- }
- // END OF DEBUG STUFF
-
- return true;
-}
#endif
diff --git a/Geo/boundaryLayersData.h b/Geo/boundaryLayersData.h
index 419dba893a6bf644f0b64c4dd3a45fe885d4bfc4..76fc336e13ac368936e1ec4b53ed6df0da005103 100644
--- a/Geo/boundaryLayersData.h
+++ b/Geo/boundaryLayersData.h
@@ -8,6 +8,7 @@
#include "SVector3.h"
#include "STensor3.h"
+#include "MElement.h"
#include "MEdge.h"
#include "MFace.h"
#include <map>
@@ -24,18 +25,13 @@ struct BoundaryLayerData
{
SVector3 _n;
std::vector<MVertex*> _column;
- std::vector<SMetric3> _metrics;
- std::vector<GFace*> _joint;
+ // std::vector<SMetric3> _metrics;
+ // std::vector<GFace*> _joint;
BoundaryLayerData(){}
BoundaryLayerData(const SVector3 & dir,
- std::vector<MVertex*> column,
- std::vector<SMetric3> metrics)
- : _n(dir), _column(column), _metrics(metrics){}
- BoundaryLayerData(const SVector3 & dir,
- std::vector<MVertex*> column,
- std::vector<SMetric3> metrics,
- std::vector<GFace*> joint)
- : _n(dir), _column(column), _metrics(metrics),_joint(joint){}
+ std::vector<MVertex*> column/*,
+ std::vector<SMetric3> metrics*/)
+ : _n(dir), _column(column) /*,_metrics(metrics)*/{}
};
struct BoundaryLayerFan
@@ -51,20 +47,6 @@ struct edgeColumn {
: _c1(c1), _c2(c2){}
};
-struct faceColumn {
- const BoundaryLayerData &_c1, &_c2, &_c3, &_c4;
- faceColumn(const BoundaryLayerData &c1,
- const BoundaryLayerData &c2,
- const BoundaryLayerData &c3)
- : _c1(c1), _c2(c2), _c3(c3), _c4(c3){}
- faceColumn(const BoundaryLayerData &c1,
- const BoundaryLayerData &c2,
- const BoundaryLayerData &c4,
- const BoundaryLayerData &c3)
- : _c1(c1), _c2(c2), _c3(c3), _c4(c4){}
-};
-
-
class BoundaryLayerColumns
{
std::map<MVertex*, BoundaryLayerFan> _fans;
@@ -79,26 +61,30 @@ public:
typedef std::multimap<MVertex*,BoundaryLayerData>::iterator iter;
typedef std::map<MVertex*, BoundaryLayerFan>::iterator iterf;
std::multimap<MVertex*, MVertex*> _non_manifold_edges;
- std::multimap<MVertex*,MTriangle*> _non_manifold_faces;
std::multimap<MEdge, SVector3, Less_Edge> _normals;
- std::multimap<MFace, SVector3, Less_Face> _normals3D;
+ void clearData () {
+ _toFirst.clear();
+ _elemColumns.clear();
+ _inverse_classification.clear();
+ _data.clear();
+ _normals.clear();
+ _non_manifold_edges.clear();
+ _elemColumns.clear();
+ _fans.clear();
+ }
+
iter begin() { return _data.begin(); }
iter end() { return _data.end(); }
iterf beginf() { return _fans.begin(); }
iterf endf() { return _fans.end(); }
BoundaryLayerColumns (){}
inline void addColumn(const SVector3 &dir, MVertex* v,
- std::vector<MVertex*> _column,
- std::vector<SMetric3> _metrics)
- {
- _data.insert (std::make_pair(v,BoundaryLayerData(dir, _column,_metrics)));
- }
- inline void addColumn(const SVector3 &dir, MVertex* v,
- std::vector<MVertex*> _column,
- std::vector<SMetric3> _metrics,
- std::vector<GFace*> _joint)
+ std::vector<MVertex*> _column//,
+ // std::vector<SMetric3> _metrics,
+ // std::vector<GFace*> _joint
+ )
{
- _data.insert (std::make_pair(v,BoundaryLayerData(dir, _column,_metrics,_joint)));
+ _data.insert (std::make_pair(v,BoundaryLayerData(dir, _column/*,_metrics,_joint*/)));
}
inline void addFan(MVertex *v, MEdge e1, MEdge e2, bool s)
{
@@ -113,7 +99,7 @@ public:
return 0;
}
- const BoundaryLayerData &getColumn(MVertex *v, MFace f) const;
+ // const BoundaryLayerData &getColumn(MVertex *v, MFace f) const;
inline const BoundaryLayerData &getColumn(MVertex *v, MEdge e) const
{
@@ -132,7 +118,7 @@ public:
return error;
}
edgeColumn getColumns(MVertex *v1, MVertex *v2 , int side);
- faceColumn getColumns(GFace *gf, MVertex *v1, MVertex *v2 , MVertex* v3, int side);
+ // faceColumn getColumns(GFace *gf, MVertex *v1, MVertex *v2 , MVertex* v3, int side);
inline int getNbColumns(MVertex *v) const { return _data.count(v); }
inline const BoundaryLayerData &getColumn(MVertex *v, int iColumn) const
{
diff --git a/Mesh/BackgroundMeshTools.cpp b/Mesh/BackgroundMeshTools.cpp
index fe6a091a6a816ef38b4ff036915cfb65d19e4049..4795e374f00a55fa40102b2f25aa8523003d738c 100644
--- a/Mesh/BackgroundMeshTools.cpp
+++ b/Mesh/BackgroundMeshTools.cpp
@@ -66,6 +66,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
GEdge *ged = (GEdge *)ge;
Crv = ged->curvature(U);
Crv = std::max(Crv, max_surf_curvature(ged, U));
+ // printf("%g %d\n",Crv, CTX::instance()->mesh.minCircPoints);
// Crv = max_surf_curvature(ged, U);
}
break;
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 3927929f18976f44bc5cd1f38b30c73088be9099..6a20fa24df4aa2368feec545ae504344fdcb627c 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -2199,7 +2199,7 @@ std::string BoundaryLayerField::getDescription()
BoundaryLayerField::BoundaryLayerField()
{
hwall_n = .1;
- hwall_t = .5;
+ // hwall_t = .5;
hfar = 1;
ratio = 1.1;
thickness = 1.e-2;
@@ -2207,20 +2207,15 @@ BoundaryLayerField::BoundaryLayerField()
tgt_aniso_ratio = 1.e10;
iRecombine = 0;
iIntersect = 0;
- options["NodesList"] = new FieldOptionList
- (nodes_id, "Indices of nodes in the geometric model", &update_needed);
options["EdgesList"] = new FieldOptionList
(edges_id, "Indices of curves in the geometric model for which a boundary "
"layer is needed", &update_needed);
- options["FacesList"] = new FieldOptionList
- (faces_id, "Indices of faces in the geometric model for which a boundary "
- "layer is needed", &update_needed);
- options["FansList"] = new FieldOptionList
- (fans_id, "Indices of edges in the geometric model for which a fan "
- "is created", &update_needed);
options["FanNodesList"] = new FieldOptionList
(fan_nodes_id, "Indices of vertices in the geometric model for which a fan "
"is created", &update_needed);
+ options["NodesList"] = new FieldOptionList
+ (nodes_id, "Indices of vertices in the geometric model for which a BL "
+ "ends", &update_needed);
options["Quads"] = new FieldOptionInt
(iRecombine, "Generate recombined elements in the boundary layer");
options["IntersectMetrics"] = new FieldOptionInt
@@ -2231,8 +2226,8 @@ BoundaryLayerField::BoundaryLayerField()
// (fan_angle, "Threshold angle for creating a mesh fan in the boundary layer");
options["AnisoMax"] = new FieldOptionDouble
(tgt_aniso_ratio, "Threshold angle for creating a mesh fan in the boundary layer");
- options["hwall_t"] = new FieldOptionDouble
- (hwall_t, "Mesh Size Tangent to the Wall");
+ // options["hwall_t"] = new FieldOptionDouble
+ // (hwall_t, "Mesh Size Tangent to the Wall");
options["ratio"] = new FieldOptionDouble
(ratio, "Size Ratio Between Two Successive Layers");
options["hfar"] = new FieldOptionDouble
@@ -2250,18 +2245,13 @@ void BoundaryLayerField::removeAttractors()
}
void BoundaryLayerField::setupFor1d(int iE) {
-
- if (faces_id_saved.empty() &&
- edges_id_saved.empty() &&
- faces_id_saved.empty() ){
- faces_id_saved = faces_id;
+ if (edges_id_saved.empty() ){
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();
@@ -2289,17 +2279,13 @@ void BoundaryLayerField::setupFor2d(int iF)
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;
+ if (edges_id_saved.empty()){
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());
@@ -2313,9 +2299,7 @@ void BoundaryLayerField::setupFor2d(int iF)
std::list<GEdge*> embedded_edges = gf->embeddedEdges();
ed.insert(ed.begin(), embedded_edges.begin(),embedded_edges.end());
-
-
- // printf("face %d has %d edges\n",iF,ed.size());
+ // 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;
@@ -2327,22 +2311,9 @@ void BoundaryLayerField::setupFor2d(int iF)
if (found){
std::list<GFace*> fc = (*it)->faces();
// one only face --> 2D --> BL
- if (fc.size() == 1) isIn = true;
+ 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 {
- 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;
- }
+ Msg::Error ("Only 2D Boundary Layers are supported (edge %d is adjacet to %d faces\n",iE,fc.size());
}
}
if (isIn){
@@ -2355,12 +2326,6 @@ void BoundaryLayerField::setupFor2d(int iF)
removeAttractors();
}
-void BoundaryLayerField::setupFor3d()
-{
- faces_id = faces_id_saved;
- removeAttractors();
-}
-
double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge)
{
if (update_needed){
@@ -2372,10 +2337,6 @@ double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge
it != edges_id.end(); ++it) {
_att_fields.push_back(new AttractorField(1,*it,300000));
}
- for(std::list<int>::iterator it = faces_id.begin();
- it != faces_id.end(); ++it) {
- _att_fields.push_back(new AttractorField(2,*it,1200));
- }
update_needed = false;
}
@@ -2442,6 +2403,9 @@ void BoundaryLayerField::operator() (AttractorField *cc, double dist,
double x, double y, double z,
SMetric3 &metr, GEntity *ge)
{
+
+ // printf("WHAT THE FUCK\n");
+
// dist = hwall -> lc = hwall * ratio
// dist = hwall (1+ratio) -> lc = hwall ratio ^ 2
// dist = hwall (1+ratio+ratio^2) -> lc = hwall ratio ^ 3
@@ -2454,8 +2418,8 @@ void BoundaryLayerField::operator() (AttractorField *cc, double dist,
const double ll1 = dist*(ratio-1) + hwall_n;
double lc_n = std::min(ll1,hfar);
- const double ll2 = dist*(ratio-1) + hwall_t;
- double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
+ // const double ll2 = dist*(ratio-1) + hwall_t;
+ double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, hfar/*std::min(ll2,hfar)*/);
lc_n = std::max(lc_n, CTX::instance()->mesh.lcMin);
lc_n = std::min(lc_n, CTX::instance()->mesh.lcMax);
@@ -2535,10 +2499,6 @@ void BoundaryLayerField::operator() (double x, double y, double z,
it != edges_id.end(); ++it) {
_att_fields.push_back(new AttractorField(1,*it,10000));
}
- for(std::list<int>::iterator it = faces_id.begin();
- it != faces_id.end(); ++it) {
- _att_fields.push_back(new AttractorField(2,*it,1200));
- }
update_needed = false;
}
diff --git a/Mesh/Field.h b/Mesh/Field.h
index bbc1ed23cafcf7f9bb51dcfc6a9efc9b5e335189..df99fad2a5124b44be33e8e5178f49ed51744901 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -139,12 +139,12 @@ class AttractorField;
class BoundaryLayerField : public Field {
private:
std::list<AttractorField *> _att_fields;
- std::list<int> nodes_id, edges_id, faces_id;
- std::list<int> faces_id_saved, edges_id_saved, nodes_id_saved, fans_id, fan_nodes_id;
+ std::list<int> nodes_id, edges_id;
+ std::list<int> edges_id_saved, nodes_id_saved, fan_nodes_id;
void operator() (AttractorField *cc, double dist, double x, double y, double z,
SMetric3 &metr, GEntity *ge);
public:
- double hwall_n,hwall_t,ratio,hfar,thickness,fan_angle;
+ double hwall_n,/*hwall_t,*/ratio,hfar,thickness,fan_angle;
double current_distance, tgt_aniso_ratio;
SPoint3 _closest_point;
int iRecombine, iIntersect;
@@ -156,10 +156,6 @@ class BoundaryLayerField : public Field {
~BoundaryLayerField() {removeAttractors();}
virtual double operator() (double x, double y, double z, GEntity *ge=0);
virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
- bool isFaceBL (int iF) const
- {
- return std::find(faces_id.begin(),faces_id.end(),iF) != faces_id.end();
- }
bool isEdgeBL (int iE) const
{
return std::find(edges_id.begin(),edges_id.end(),iE) != edges_id.end();
@@ -169,22 +165,17 @@ class BoundaryLayerField : public Field {
return std::find(edges_id_saved.begin(),edges_id_saved.end(),iE)
!= edges_id_saved.end();
}
- bool isFan (int iE) const
- {
- return std::find(fans_id.begin(),fans_id.end(),iE) != fans_id.end();
- }
bool isFanNode (int iV) const
{
return std::find(fan_nodes_id.begin(),fan_nodes_id.end(),iV) != fan_nodes_id.end();
}
- bool isVertexBL (int iV) const
+ bool isEndNode (int iV) const
{
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/filterElements.cpp b/Mesh/filterElements.cpp
index 380838fcadd317e72300195cef16e2f3b51c3c80..5990306d09dd039acb56e153973925dc94b8a42c 100644
--- a/Mesh/filterElements.cpp
+++ b/Mesh/filterElements.cpp
@@ -8,14 +8,12 @@
#include <set>
#include "GmshDefines.h"
#include "filterElements.h"
-#include "fullMatrix.h"
-#include "GaussIntegration.h"
#include "MElement.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
-#include "MPrism.h"
-#include "MHexahedron.h"
+#include "MLine.h"
#include "rtree.h"
+#include "robustPredicates.h"
void MElementBB(void *a, double *min, double *max);
int MElementInEle(void *a, double *x);
@@ -32,89 +30,113 @@ struct MElement_Wrapper
}
};
-static bool inBB (double *mi, double *ma, double *x){
- if (x[0] < mi[0])return false;
- if (x[1] < mi[1])return false;
- if (x[2] < mi[2])return false;
- if (x[0] > ma[0])return false;
- if (x[1] > ma[1])return false;
- if (x[2] > ma[2])return false;
- return true;
+
+/** OVERLAP TEST IN 2D
+ only use predicate orient2D :-)
+
+2 cases
+
+1-->
+
+ *---------*
+ | |
+ | *--------*
+ | | | |
+ *---------*
+ | | |
+ *--------*
+
+2-->
+
+**/
+
+inline double orientationTest (double a[2], double b[2], double c[2]){
+ double s = -robustPredicates::orient2d(a,b,c);
+ return s > 0 ? 1.0 : s < 0 ? -1.0 : 0.0;
}
-bool rtree_callback(MElement *e1,void* pe2)
-{
- MElement_Wrapper *wrapper = static_cast<MElement_Wrapper*>(pe2);
- MElement *e2 = wrapper->_e;
+inline double orientationTest (MVertex *va, MVertex *vb, MVertex *vc){
+ double a[2] = {va->x(),va->y()};
+ double b[2] = {vb->x(),vb->y()};
+ double c[2] = {vc->x(),vc->y()};
+ return orientationTest (a,b,c);
+}
- if (std::binary_search(wrapper->_notOverlap.begin(),wrapper->_notOverlap.end(),e1))
- return true;
+inline double orientationTest (SVector3 &va, SVector3 &vb, SVector3 &vc){
+ double a[2] = {va.x(),va.y()};
+ double b[2] = {vb.x(),vb.y()};
+ double c[2] = {vc.x(),vc.y()};
+ return orientationTest (a,b,c);
+}
- for (int i=0;i<e1->getNumVertices();i++){
- for (int j=0;j<e2->getNumVertices();j++){
- if(e1->getVertex(i) == e2->getVertex(j))return true;
- }
- }
- double min2[3],max2[3];
- MElementBB(e2,min2,max2);
+bool intersectEdge2d(const MEdge &ed1, const MEdge &ed2) {
- for (int i=0;i<e1->getNumVertices();i++){
- SPoint3 xyz (e1->getVertex(i)->x(),e1->getVertex(i)->y(),e1->getVertex(i)->z());
- if (inBB (min2,max2,xyz)){
- if (MElementInEle(e2,xyz)){
- wrapper->_overlap = true;
- return false;
- }
- }
- }
+ /* SVector3 a1(ed1.getVertex(0)->x(), ed1.getVertex(0)->y(), 0);
+ SVector3 a2(ed1.getVertex(1)->x(), ed1.getVertex(1)->y(), 0);
- double min1[3],max1[3];
- MElementBB(e1,min1,max1);
- for (int i=0;i<e2->getNumVertices();i++){
- SPoint3 xyz (e2->getVertex(i)->x(),e2->getVertex(i)->y(),e2->getVertex(i)->z());
- if (inBB(min1,max1,xyz)){
- if (MElementInEle(e1,xyz)){
- wrapper->_overlap = true;
- return false;
- }
- }
- }
+ SVector3 b1(ed2.getVertex(0)->x(), ed2.getVertex(0)->y(), 0);
+ SVector3 b2(ed2.getVertex(1)->x(), ed2.getVertex(1)->y(), 0);
- // return false;
+ SVector3 a = (a1+a2)*0.5;
+ SVector3 b = (b1+b2)*0.5;
- fullMatrix<double> pts1,pts2;
- fullVector<double> weights;
+ SVector3 a1l = a + (a2-a1)*0.55;
+ SVector3 a2l = a - (a2-a1)*0.55;
+
+ SVector3 b1l = b + (b2-b1)*0.55;
+ SVector3 b2l = b - (b2-b1)*0.55;
- gaussIntegration::get (e1->getType(),1,pts1,weights);
- gaussIntegration::get (e2->getType(),1,pts2,weights);
+ if (orientationTest (a1l, a2l, b1l) * orientationTest (a1l, a2l, b2l) <= 0 &&
+ orientationTest (b1l, b2l, a1l) * orientationTest (b1l, b2l, a2l) <= 0 )
+ return true;
+ */
+ if (ed1.getVertex(0) == ed2.getVertex(0) ||
+ ed1.getVertex(0) == ed2.getVertex(1) ||
+ ed1.getVertex(1) == ed2.getVertex(0) ||
+ ed1.getVertex(1) == ed2.getVertex(1) ) return false;
+
+ if ((orientationTest (ed1.getVertex(0),ed1.getVertex(1),ed2.getVertex(0)) *
+ orientationTest (ed1.getVertex(0),ed1.getVertex(1),ed2.getVertex(1)) <= 0) &&
+ (orientationTest (ed2.getVertex(0),ed2.getVertex(1),ed1.getVertex(0)) *
+ orientationTest (ed2.getVertex(0),ed2.getVertex(1),ed1.getVertex(1)) <= 0))
+ return true;
+ return false;
+}
- for (int i=0;i<pts1.size1();i++){
- double u = pts1(i,0);
- double v = pts1(i,1);
- double w = pts1(i,2);
- SPoint3 xyz;
- e1->pnt(u,v,w,xyz);
- if (inBB(min2,max2,xyz)){
- if (MElementInEle(e2,xyz)){
- wrapper->_overlap = true;
- return false;
+bool overlap2D (MElement *e1, MElement *e2) {
+ for (int i=0;i<e1->getNumEdges();i++){
+ MEdge ed1 = e1->getEdge (i);
+ for (int j=0;j<e2->getNumEdges();j++){
+ MEdge ed2 = e2->getEdge (j);
+ if (intersectEdge2d(ed1,ed2)){
+ // printf("apero time \n");
+ return true;
}
}
}
- for (int i=0;i<pts2.size1();i++){
- double u = pts2(i,0);
- double v = pts2(i,1);
- double w = pts2(i,2);
- SPoint3 xyz;
- e2->pnt(u,v,w,xyz);
- if (inBB(min1,max1,xyz)){
- if (MElementInEle(e1,xyz)){
- wrapper->_overlap = true;
- return false;
- }
- }
+ return false;
+}
+
+bool rtree_callback(MElement *e1,void* pe2)
+{
+ MElement_Wrapper *wrapper = static_cast<MElement_Wrapper*>(pe2);
+ MElement *e2 = wrapper->_e;
+
+ if (std::binary_search(wrapper->_notOverlap.begin(),wrapper->_notOverlap.end(),e1))
+ return true;
+
+ // for (int i=0;i<e1->getNumVertices();i++){
+ // for (int j=0;j<e2->getNumVertices();j++){
+ // if(e1->getVertex(i) == e2->getVertex(j))return true;
+ // }
+ // }
+
+ if (e1->getDim() <= 2 && e2->getDim() <= 2) {
+ wrapper->_overlap = overlap2D (e1,e2);
+ return !wrapper->_overlap;
}
+ Msg::Error("overlapping of elements in 3D not done yet");
return true;
}
@@ -133,75 +155,47 @@ void filterColumns(std::vector<MElement*> &elem,
for (std::map<MElement*,std::vector<MElement*> >::iterator it = _elemColumns.begin();
it !=_elemColumns.end() ; ++it){
const std::vector<MElement*> &c = it->second;
- unsigned int MAX = c.size() - 1;
+ unsigned int MAX = c.size();
+ // printf("size of column %d\n",c.size());
for (unsigned int i=0;i<c.size(); i++){
if (!std::binary_search(elem.begin(),elem.end(),c[i])){
MAX = i - 1;
break;
}
}
+ if (!MAX)MAX=1;
// printf("MAX = %d c = %d\n",MAX,c.size());
for (unsigned int i=0;i<MAX;i++){
toKeep.push_back(c[i]);
}
- for (unsigned int i=MAX;i<c.size();i++){
+ // for (unsigned int i=MAX;i<c.size();i++){
/// FIXME !!!
// delete c[i];
- }
+ // }
}
- printf("%d --> %d\n", (int)elem.size(), (int)toKeep.size());
+ // printf("%d --> %d\n", (int)elem.size(), (int)toKeep.size());
elem = toKeep;
}
-void filterOverlappingElements (std::vector<MTriangle*> &blTris,
- std::vector<MQuadrangle*>&blQuads,
- std::map<MElement*,std::vector<MElement*> > &_elemColumns,
- std::map<MElement*,MElement*> &_toFirst)
-{
- std::vector<MElement*> vvv;
- vvv.insert(vvv.begin(),blTris.begin(),blTris.end());
- vvv.insert(vvv.begin(),blQuads.begin(),blQuads.end());
- Less_Partition lp;
- std::sort(vvv.begin(),vvv.end(), lp);
- filterOverlappingElements (vvv,_elemColumns,_toFirst);
- filterColumns (vvv,_elemColumns);
- blTris.clear();
- blQuads.clear();
- for (unsigned int i=0;i<vvv.size();i++){
- if (vvv[i]->getType() == TYPE_TRI)blTris.push_back((MTriangle*)vvv[i]);
- else if (vvv[i]->getType() == TYPE_QUA)blQuads.push_back((MQuadrangle*)vvv[i]);
- }
-}
-
-void filterOverlappingElements (std::vector<MPrism*> &blPrisms,
- std::vector<MHexahedron*>&blHexes,
- std::map<MElement*,std::vector<MElement*> > &_elemColumns,
- std::map<MElement*,MElement*> &_toFirst)
-{
- printf("filtering !!\n");
- std::vector<MElement*> vvv;
- vvv.insert(vvv.begin(),blPrisms.begin(),blPrisms.end());
- vvv.insert(vvv.begin(),blHexes.begin(),blHexes.end());
- Less_Partition lp;
- std::sort(vvv.begin(),vvv.end(), lp);
- filterOverlappingElements (vvv,_elemColumns,_toFirst);
- filterColumns (vvv,_elemColumns);
- blPrisms.clear();
- blHexes.clear();
- for (unsigned int i=0;i<vvv.size();i++){
- if (vvv[i]->getType() == TYPE_PRI)blPrisms.push_back((MPrism*)vvv[i]);
- else if (vvv[i]->getType() == TYPE_HEX)blHexes.push_back((MHexahedron*)vvv[i]);
- }
-}
-void filterOverlappingElements (std::vector<MElement*> &els,
- std::map<MElement*,std::vector<MElement*> > &_elemColumns,
- std::map<MElement*,MElement*> &_toFirst )
+static void filterOverlappingElements (std::vector<MLine*> &lines,
+ std::vector<MElement*> &els,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst )
{
std::vector<MElement*> newEls;
RTree<MElement*,double,3,double> rtree;
+
+ for (unsigned int i=0;i<lines.size();i++){
+ MElement *e = lines[i];
+ double _min[3],_max[3];
+ MElementBB(e,_min,_max);
+ rtree.Insert(_min,_max,e);
+ }
+
+
for (unsigned int i=0;i<els.size();i++){
MElement *e = els[i];
double _min[3],_max[3];
@@ -210,7 +204,7 @@ void filterOverlappingElements (std::vector<MElement*> &els,
MElement_Wrapper w(e,_elemColumns[first]);
rtree.Search(_min,_max,rtree_callback,&w);
if (w._overlap){
- delete e;
+ // delete e;
}
else {
rtree.Insert(_min,_max,e);
@@ -220,3 +214,28 @@ void filterOverlappingElements (std::vector<MElement*> &els,
els = newEls;
}
+// WE SHOULD ADD THE BOUNDARY OF THE DOMAIN IN ORDER TO AVOID
+// ELEMENTS THAT ARE OUTSIDE THE DOMAIN --> FIXME
+
+void filterOverlappingElements (std::vector<MLine*> &bdry,
+ std::vector<MTriangle*> &blTris,
+ std::vector<MQuadrangle*>&blQuads,
+ std::map<MElement*,std::vector<MElement*> > &_elemColumns,
+ std::map<MElement*,MElement*> &_toFirst)
+{
+ std::vector<MElement*> vvv;
+ vvv.insert(vvv.begin(),blTris.begin(),blTris.end());
+ vvv.insert(vvv.begin(),blQuads.begin(),blQuads.end());
+ Less_Partition lp;
+ std::sort(vvv.begin(),vvv.end(), lp);
+ filterOverlappingElements (bdry,vvv,_elemColumns,_toFirst);
+ filterColumns (vvv,_elemColumns);
+ blTris.clear();
+ blQuads.clear();
+ for (unsigned int i=0;i<vvv.size();i++){
+ if (vvv[i]->getType() == TYPE_TRI)blTris.push_back((MTriangle*)vvv[i]);
+ else if (vvv[i]->getType() == TYPE_QUA)blQuads.push_back((MQuadrangle*)vvv[i]);
+ }
+
+
+}
diff --git a/Mesh/filterElements.h b/Mesh/filterElements.h
index 60f16ff06a381eeb15843efdfe12a0c770265c0a..bd18e2937c3ab9194522d5bbd0b6c0fc5240c9e6 100644
--- a/Mesh/filterElements.h
+++ b/Mesh/filterElements.h
@@ -10,19 +10,12 @@
#include <vector>
class MElement;
-class MPrism;
-class MHexahedron;
class MTriangle;
class MQuadrangle;
+class MLine;
-void filterOverlappingElements(std::vector<MElement*> &,
- std::map<MElement*,std::vector<MElement*> > &_elemColumns,
- std::map<MElement*,MElement*> &_toFirst);
-void filterOverlappingElements(std::vector<MPrism*> &blPrisms,
- std::vector<MHexahedron*>&blHexes,
- std::map<MElement*,std::vector<MElement*> > &_elemColumns,
- std::map<MElement*,MElement*> &_toFirst);
-void filterOverlappingElements(std::vector<MTriangle*> &blTris,
+void filterOverlappingElements(std::vector<MLine*> &_lines,
+ std::vector<MTriangle*> &blTris,
std::vector<MQuadrangle*>&blQuads,
std::map<MElement*,std::vector<MElement*> > &_elemColumns,
std::map<MElement*,MElement*> &_toFirst);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 068d23662c7394d172e0d45ba4664998a1972a0e..5330284cc749a191ef843c6ac05fbeef2ffd7397 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -78,34 +78,36 @@ static double smoothPrimitive(GEdge *ge, double alpha,
static double F_LcB(GEdge *ge, double t)
{
- BoundaryLayerField *blf = 0;
+ /* 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 = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
- if (blf){
+ /* 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;
@@ -119,13 +121,13 @@ static double F_Lc(GEdge *ge, double t)
lc_here = BGM_MeshSize(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
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);
@@ -445,7 +447,24 @@ static void filterPoints (GEdge*ge, int nMinimumPoints)
}
}
-void meshGEdge::operator() (GEdge *ge)
+static void createPoints ( GVertex *gv, GEdge *ge, BoundaryLayerField *blf ,
+ std::vector<MVertex*>& v, const SVector3 &dir){
+ double L = blf->hwall_n;
+ while (1){
+ if (L > blf->thickness) break;
+ SPoint3 p (gv->x() + dir.x() * L, gv->y() + dir.y() * L, 0.0);
+ v.push_back(new MEdgeVertex (p.x(), p.y(), p.z(), ge, ge->parFromPoint(p), blf->hfar));
+ int ith = v.size() ;
+ L+= blf->hwall_n * pow (blf->ratio, ith);
+ // printf("parameter %g length %g\n",ge->parFromPoint(p),L);
+ }
+}
+
+void addBoundaryLayerPoints (GEdge *ge,
+ double &t_begin, // may change the left parameter of the interval
+ double &t_end, // may change the right parameter of the interval
+ std::vector<MVertex*> &_addBegin, // additional points @ left
+ std::vector<MVertex*> &_addEnd) // additional points @ right
{
BoundaryLayerField *blf = 0;
#if defined(HAVE_ANN)
@@ -454,7 +473,46 @@ void meshGEdge::operator() (GEdge *ge)
blf = dynamic_cast<BoundaryLayerField*> (bl_field);
if (blf) blf->setupFor1d(ge->tag());
#endif
+ if (!blf) return;
+ if (blf->isEdgeBL(ge->tag()))return;
+ SVector3 dir ( ge->getEndVertex()->x() - ge->getBeginVertex()->x(),
+ ge->getEndVertex()->y() - ge->getBeginVertex()->y(),
+ ge->getEndVertex()->z() - ge->getBeginVertex()->z());
+ dir.normalize();
+ GVertex *gvb = ge->getBeginVertex();
+ GVertex *gve = ge->getEndVertex();
+ if (blf->isEndNode(gvb->tag())){
+ if (ge->geomType() != GEntity::Line){
+ Msg::Error ("Boundary layer end point %d should lie on a straight line", gvb->tag());
+ return;
+ }
+ createPoints (gvb, ge, blf, _addBegin, dir);
+ if (!_addBegin.empty())_addBegin[_addBegin.size()-1]->getParameter(0,t_begin);
+ }
+ if (blf->isEndNode(gve->tag())){
+ if (ge->geomType() != GEntity::Line){
+ Msg::Error ("Boundary layer end point %d should lie on a straight line", gve->tag());
+ return;
+ }
+ createPoints (gve, ge, blf, _addEnd, dir * -1.0);
+ if (!_addEnd.empty())_addEnd[_addEnd.size()-1]->getParameter(0,t_end);
+ }
+
+ // printf("Edge %d ยง %d %d points added (%g %g)-\n", ge->tag(), _addBegin.size(),_addEnd.size(),t_begin,t_end);
+}
+void meshGEdge::operator() (GEdge *ge)
+{
+ /*
+ 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);
+ if (blf) blf->setupFor1d(ge->tag());
+#endif
+ */
+
ge->model()->setCurrentMeshEntity(ge);
// if(ge->geomType() == GEntity::DiscreteCurve) return;
@@ -488,6 +546,10 @@ void meshGEdge::operator() (GEdge *ge)
double t_begin = bounds.low();
double t_end = bounds.high();
+ // if a BL is ending at one of the ends, then create specific points
+ std::vector<MVertex*> _addBegin, _addEnd;
+ addBoundaryLayerPoints (ge, t_begin, t_end, _addBegin, _addEnd);
+
// first compute the length of the curve by integrating one
double length;
std::vector<IntPoint> Points;
@@ -541,7 +603,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);
filterMinimumN = ge->minimumMeshSegments() + 1;
N = std::max(filterMinimumN, (int)(a + 1.99));
}
@@ -620,8 +682,20 @@ void meshGEdge::operator() (GEdge *ge)
mesh_vertices.resize(NUMP - 1);
}
+ // Boundary Layer points are added
+ {
+ std::vector<MVertex*> vv;
+ vv.insert(vv.end(), _addBegin.begin(), _addBegin.end());
+ vv.insert(vv.end(), mesh_vertices.begin(), mesh_vertices.end());
+ for (unsigned int i=0;i<_addEnd.size();i++)
+ vv.push_back(_addEnd[_addEnd.size()-1-i]);
+ // vv.insert(vv.end(), _addEnd.rend(), _addEnd.rbegin());
+ mesh_vertices = vv;
+ }
+
// printf("%ld ----> ", ge->mesh_vertices.size());
- filterPoints (ge, filterMinimumN - 2);
+ if (_addBegin.empty() && _addEnd.empty())
+ filterPoints (ge, filterMinimumN - 2);
// printf("%ld \n", ge->mesh_vertices.size());
for(unsigned int i = 0; i < mesh_vertices.size() + 1; i++){
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 4a90b7c24c09082f7dc404ec38d650dd730ca609..450f36d3836982e10707ac7213119613410e8fd2 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -662,8 +662,12 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
FILE *ff2 = Fopen ("tato.pos","w");
if(ff2) fprintf(ff2,"View \" \"{\n");
std::set<MVertex*> verts;
+
+ std::vector<MLine*> _lines;
+
while(ite != edges.end()){
for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ _lines.push_back((*ite)->lines[i]);
MVertex *v1 = (*ite)->lines[i]->getVertex(0);
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
MEdge dv(v1,v2);
@@ -691,6 +695,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
//avoid convergent errors
if (dv2.length() < 0.03 * dv.length())break;
MQuadrangle *qq = new MQuadrangle(v11,v21,v22,v12);
+ qq->setPartition (l+1);
myCol.push_back(qq);
blQuads.push_back(qq);
if(ff2)
@@ -732,6 +737,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
}
if (v11 != v12){
MQuadrangle *qq = new MQuadrangle(v11,v12,v22,v21);
+ qq->setPartition (l+1);
myCol.push_back(qq);
blQuads.push_back(qq);
if(ff2)
@@ -743,6 +749,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
}
else {
MTriangle *qq = new MTriangle(v,v22,v21);
+ qq->setPartition (l+1);
myCol.push_back(qq);
blTris.push_back(qq);
if(ff2)
@@ -762,7 +769,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
fclose(ff2);
}
- // filterOverlappingElements (blTris,blQuads,_columns->_elemColumns,_columns->_toFirst);
+ filterOverlappingElements (_lines, blTris,blQuads,_columns->_elemColumns,_columns->_toFirst);
for (unsigned int i = 0; i < blQuads.size();i++){
addOrRemove(blQuads[i]->getVertex(0),blQuads[i]->getVertex(1),bedges, removed);
@@ -803,7 +810,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
deMeshGFace kil_;
kil_(gf);
meshGenerator(gf, 0, 0, true , false, &hop);
-
+
gf->quadrangles = blQuads;
gf->triangles.insert(gf->triangles.begin(),blTris.begin(),blTris.end());
gf->mesh_vertices.insert(gf->mesh_vertices.begin(),verts.begin(),verts.end());
@@ -1543,8 +1550,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// BDS mesh is passed in order not to recompute local coordinates of
// vertices
if(algoDelaunay2D(gf) && !onlyInitialMesh){
- if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL)
+ if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL){
bowyerWatsonFrontal(gf);
+ }
else if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD){
bowyerWatsonFrontalLayers(gf,true);
}
diff --git a/Mesh/meshRefine.cpp b/Mesh/meshRefine.cpp
index aa5dbe3f9d01a80ab3c4e670b7a8fe86252e7882..8ebafac03340947d66d3204e8ea323e03b2620c8 100644
--- a/Mesh/meshRefine.cpp
+++ b/Mesh/meshRefine.cpp
@@ -449,8 +449,8 @@ static void Subdivide(GRegion *gr, bool splitIntoHexas, faceContainer &faceVerti
void RefineMesh(GModel *m, bool linear, bool splitIntoQuads, bool splitIntoHexas)
{
- // splitIntoQuads = true;
- // splitIntoHexas = true;
+ splitIntoQuads = true;
+ splitIntoHexas = true;
Msg::StatusBar(true, "Refining mesh...");
double t1 = Cpu();
diff --git a/benchmarks/2d/BL1.geo b/benchmarks/2d/BL1.geo
new file mode 100755
index 0000000000000000000000000000000000000000..6030cc9c59a386e65f686756c820b23ad87356c3
--- /dev/null
+++ b/benchmarks/2d/BL1.geo
@@ -0,0 +1,54 @@
+Mesh.CharacteristicLengthFromCurvature=1;
+Mesh.MinimumCirclePoints=20;
+l = .2;
+//+
+Point(1) = {-1.1, -0.3, 0, l};
+//+
+Point(2) = {-0.3, -0, 0, l};
+//+
+Point(3) = {1, -0.3, 0, l};
+//+
+Point(4) = {-0.3, -0.6, 0, l};
+//+
+Point(5) = {0.7, -0.2, 0, l};
+//+
+Point(6) = {0.7, -0.4, 0, l};
+//+
+Spline(1) = {3, 5, 2, 1, 4, 6, 3};
+//+
+Point(7) = {-1.4, 0.5, 0, l};
+//+
+Point(8) = {-1.4, -1.1, 0, l};
+//+
+Point(9) = {1.4, -1.1, 0, l};
+//+
+Point(10) = {1.4, 0.5, 0, l};
+//+
+Line(2) = {7, 10};
+//+
+Line(3) = {9, 10};
+//+
+Line(4) = {9, 8};
+//+
+Line(5) = {8, 7};
+//+
+Line Loop(6) = {2, -3, 4, 5};
+//+
+Line Loop(7) = {1};
+//+
+Plane Surface(8) = {6, 7};
+//+
+Field[1] = BoundaryLayer;
+//+
+Field[1].EdgesList = {1};
+//+
+Field[1].hfar = 0.1;
+//+
+Field[1].hwall_n = 0.01;
+//+
+Field[1].thickness = 0.1;
+//+
+BoundaryLayer Field = 1;
+
+//+
+Field[1].ratio = 1.4;
diff --git a/benchmarks/2d/BL2.geo b/benchmarks/2d/BL2.geo
new file mode 100644
index 0000000000000000000000000000000000000000..49e1815d23799a7573f49c06dc2ebcee135993ca
--- /dev/null
+++ b/benchmarks/2d/BL2.geo
@@ -0,0 +1,63 @@
+Mesh.CharacteristicLengthFromCurvature=1;
+Mesh.MinimumCirclePoints=20;
+l = .2;
+//+
+Point(1) = {-1.1, -0.3, 0, l};
+//+
+Point(2) = {-0.3, -0, 0, l};
+//+
+Point(3) = {1, -0.3, 0, l};
+//+
+Point(4) = {-0.3, -0.6, 0, l};
+//+
+Point(5) = {0.7, -0.2, 0, l};
+//+
+Point(6) = {0.7, -0.4, 0, l};
+//+
+Spline(1) = {1, 4, 6, 3};
+Spline(11) = {3 , 5 , 2 , 1};
+//+
+Point(7) = {-1.4, .5, 0, l};
+//+
+Point(8) = {-1.4, -1.1, 0, l};
+//+
+Point(9) = {3.4, -1.1, 0, l};
+//+
+Point(10) = {3.4, .5, 0, l};
+//+
+Line(2) = {7, 10};
+//+
+Line(3) = {9, 10};
+//+
+Line(4) = {9, 8};
+//+
+Line(5) = {8, 7};
+//+
+Line Loop(6) = {2, -3, 4, 5};
+//+
+Line Loop(7) = {1, 11};
+//+
+Plane Surface(8) = {6, 7};
+//+
+Field[1] = BoundaryLayer;
+//+
+Field[1].EdgesList = {1, 11, 9};
+//+
+Field[1].hfar = 0.1;
+//+
+Field[1].hwall_n = 0.01;
+//+
+Field[1].thickness = 0.1;
+//+
+BoundaryLayer Field = 1;
+
+//+
+Field[1].ratio = 1.4;
+//+
+Point(11) = {2.2, -0.3, 0, l};
+//+
+Line(9) = {3, 11};
+Line {9} In Surface {8};//+
+Field[1].FanNodesList = {1};
+//+
+Field[1].Quads = 1;
diff --git a/benchmarks/2d/BL3.geo b/benchmarks/2d/BL3.geo
new file mode 100644
index 0000000000000000000000000000000000000000..71c6a54732ef134310f3192cbb0948bff2bb1721
--- /dev/null
+++ b/benchmarks/2d/BL3.geo
@@ -0,0 +1,84 @@
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthMin = 0.1;
+Mesh.CharacteristicLengthMax = 0.1;
+//+
+Point(1) = {-1.7, -0.3, 0, .1};
+//+
+Point(2) = {-1.2, -0, 0, .1};
+//+
+Point(3) = {-0.5, -0.3, 0, .1};
+//+
+Point(4) = {-0.4, -0.7, 0, .1};
+//+
+Point(5) = {-0.8, -0.8, 0, .1};
+//+
+Point(6) = {-1.2, -0.5, 0, .1};
+//+
+Point(7) = {-1.8, -0.5, 0, .1};
+//+
+Point(8) = {-1.1, 0.6, 0, .1};
+//+
+Point(9) = {-0.5, 0.7, 0, .1};
+//+
+Point(10) = {-0.1, 0.2, 0, .1};
+//+
+Point(11) = {-0.3, -0.1, 0, .1};
+//+
+Point(12) = {-0.7, 0.2, 0, .1};
+//+
+Point(13) = {-1.1, 0.4, 0, 1.0};
+//+
+Point(14) = {-1.5, 0.5, 0, .1};
+//+
+Point(15) = {-2.3, 1.1, 0, .1};
+//+
+Point(16) = {-2.3, -1, 0, .1};
+//+
+Point(17) = {0.7, -1, 0, .1};
+//+
+Point(18) = {0.5, 1.1, 0, .1};
+//+
+Line(1) = {15, 16};
+//+
+Line(2) = {16, 17};
+//+
+Line(3) = {17, 18};
+//+
+Line(4) = {18, 15};
+//+
+Spline(5) = {14, 8, 9, 10, 11, 12, 13, 14};
+//+
+Spline(6) = {2, 3, 4, 5, 6, 7, 1, 2};
+//+
+Field[1] = BoundaryLayer;
+//+
+Field[1].EdgesList = {5, 6, 1};
+//+
+Field[1].hfar = 0.1;
+//+
+Field[1].hwall_n = 0.001;
+//+
+Field[1].ratio = 1.4;
+//+
+Field[1].thickness = 0.38;
+BoundaryLayer Field = 1;
+//+
+Line Loop(7) = {1, 2, 3, 4};
+//+
+Line Loop(8) = {5};
+//+
+Line Loop(9) = {6};
+//+
+Plane Surface(10) = {7, 8, 9};
+//+
+Characteristic Length {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14} = 0.1;
+//+
+Field[1].FanNodesList = {14};
+//+
+Characteristic Length {17, 16, 5, 4, 6, 7, 3, 1, 11, 2, 10, 12, 13, 14, 8, 9, 15, 18} = 0.1;
+//+
+Field[1].Quads = 1;
+//+
+//Field[1].thickness = 1;
+//+
+Field[1].NodesList = {15, 16};
diff --git a/benchmarks/3d/hex.geo b/benchmarks/3d/hex.geo
index a838d0ff2d5e56d4bbc774f154f4bb518fdc954c..35c542e57b029f9e33db3acbf0df19cf4c9a2a5c 100644
--- a/benchmarks/3d/hex.geo
+++ b/benchmarks/3d/hex.geo
@@ -4,7 +4,7 @@ lc = 0.3;
// example of a purely hexahedral mesh using only transfinite
// mesh constraints
z=1;
-deform=0.0;
+deform=0;
Point(1) = {-2-deform,0,0,lc};
Point(2) = {-1,0,0,lc};
@@ -40,8 +40,8 @@ Line Loop(24) = {6,-12,3,10};
Ruled Surface(25) = {24};
Surface Loop(1) = {17,-25,-23,-21,19,15};
Volume(1) = {1};
-Transfinite Line{1:4,6:9} = 20;
-Transfinite Line{10:13} = 20;
+Transfinite Line{1:4,6:9} = 2;
+Transfinite Line{10:13} = 2;
Transfinite Surface {15} = {1,2,3,4};
Transfinite Surface {17} = {5,6,7,8};
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 34c59b566c715c494a9199681e3f1942090e5e51..f71e14d54207d47cb6ed638f2d376061f6eb3178 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -308,9 +308,10 @@ bool getTopPrimVertices(std::map<MVertex*, std::vector<MElement *> > &vertex2ele
topPrimVert = std::vector<MVertex*>(nbVert);
for(int i = 0; i < nbVert; i++) {
if (blc && (blc->size() > 1)) { // Is there BL data?
- const BoundaryLayerData &c = blc->getColumn(baseVert[i], baseBnd);
- if (c._column.size() == 0) return false; // Give up element if column not found
- topPrimVert[i] = *(c._column.end()-2);
+ return false;
+ // const BoundaryLayerData &c = blc->getColumn(baseVert[i], baseBnd);
+ // if (c._column.size() == 0) return false; // Give up element if column not found
+ // topPrimVert[i] = *(c._column.end()-2);
}
else { // No BL data, look for columns of vertices
std::map<MVertex*, SVector3>::const_iterator itNorm = normVert.find(baseVert[i]);
@@ -553,12 +554,12 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
if (blf->isEdgeBL((*itEd)->tag())) // Already skip model edge if no BL there
entities.insert(std::pair<GEntity*,GEntity*>(entity, *itEd));
}
- else if (p.dim == 3) { // "Domain" region?
- std::list<GFace*> faces = entity->faces();
- for (std::list<GFace*>::iterator itF = faces.begin(); itF != faces.end(); itF++) // Loop over model boundary faces
- if (blf->isFaceBL((*itF)->tag())) // Already skip model face if no BL there
- entities.insert(std::pair<GEntity*,GEntity*>(entity, *itF));
- }
+ // else if (p.dim == 3) { // "Domain" region?
+ // std::list<GFace*> faces = entity->faces();
+ // for (std::list<GFace*>::iterator itF = faces.begin(); itF != faces.end(); itF++) // Loop over model boundary faces
+ // if (blf->isFaceBL((*itF)->tag())) // Already skip model face if no BL there
+ // entities.insert(std::pair<GEntity*,GEntity*>(entity, *itF));
+ // }
}
}
else { // No BL field