From 2033212d5fc57a65f550c7a82ce80d41dbd67845 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 20 Aug 2013 13:51:34 +0000
Subject: [PATCH] 3D Boundary layer mesh is in better shape IMPORTANT : 3D
Delaunay Meshing now acts on connected volumes separately, allowing
to mesh regions that are not connected topologically but that may
intersect geometrically
---
Geo/GEdge.cpp | 18 +-
Geo/GEdge.h | 3 +
Geo/GFace.h | 7 +
Geo/GModel.h | 3 -
Geo/GModelIO_OCC.cpp | 2 +
Geo/GVertex.cpp | 29 ++
Geo/GVertex.h | 6 +
Geo/MFace.h | 2 +
Geo/boundaryLayersData.cpp | 677 +++++++++++++++++++++++---
Geo/boundaryLayersData.h | 116 ++++-
Mesh/CMakeLists.txt | 2 +-
Mesh/Field.cpp | 102 +++-
Mesh/Field.h | 5 +
Mesh/Generator.cpp | 35 +-
Mesh/meshGEdge.cpp | 25 -
Mesh/meshGFace.cpp | 74 +--
Mesh/meshGFaceDelaunayInsertion.cpp | 34 +-
Mesh/meshGFaceDelaunayInsertion.h | 1 +
Mesh/meshGFaceOptimize.cpp | 1 +
Mesh/meshGRegion.cpp | 608 ++++++++++++++++++++++-
Mesh/meshGRegion.h | 2 +
Mesh/meshGRegionBoundaryLayers.cpp | 212 --------
Mesh/meshGRegionDelaunayInsertion.cpp | 80 +--
Mesh/meshGRegionDelaunayInsertion.h | 4 +-
24 files changed, 1619 insertions(+), 429 deletions(-)
delete mode 100644 Mesh/meshGRegionBoundaryLayers.cpp
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 26245d295e..5ac4cd5134 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -97,7 +97,8 @@ void GEdge::resetMeshAttributes()
void GEdge::addFace(GFace *e)
{
- l_faces.push_back(e);
+ if (std::find(l_faces.begin(), l_faces.end(), e) == l_faces.end())
+ l_faces.push_back(e);
}
void GEdge::delFace(GFace *e)
@@ -485,3 +486,18 @@ void GEdge::replaceEndingPoints(GVertex *replOfv0, GVertex *replOfv1)
v1 = replOfv1;
}
}
+
+// regions that bound this entity or that this entity bounds.
+std::list<GRegion*> GEdge::regions() const
+{
+ std::list<GFace*> _faces = faces();
+ std::list<GFace*>::const_iterator it = _faces.begin();
+ std::set<GRegion*> _r;
+ for ( ; it != _faces.end() ; ++it){
+ std::list<GRegion*> temp = (*it)->regions();
+ _r.insert (temp.begin(), temp.end());
+ }
+ std::list<GRegion*> ret;
+ ret.insert (ret.begin(), _r.begin(), _r.end());
+ return ret;
+}
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index b43736925b..b6bb1c77d5 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -72,6 +72,9 @@ class GEdge : public GEntity {
// get the oriented bounding box
virtual SOrientedBoundingBox getOBB();
+ // regions that are boundedby this entity
+ virtual std::list<GRegion*> regions() const;
+
// faces that this entity bounds
virtual std::list<GFace*> faces() const { return l_faces; }
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 0c37ed4a7d..afc8765300 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -17,6 +17,7 @@
#include "SVector3.h"
#include "Pair.h"
#include "Numeric.h"
+#include "boundaryLayersData.h"
class MElement;
class MTriangle;
@@ -49,6 +50,7 @@ class GFace : public GEntity
// replace edges (for gluing) for specific modelers, we have to
// re-create internal data
virtual void replaceEdgesInternal(std::list<GEdge*> &){}
+ BoundaryLayerColumns _columns;
public: // this will become protected or private
std::list<GEdgeLoop> edgeLoops;
@@ -332,6 +334,11 @@ class GFace : public GEntity
void addTriangle(MTriangle *t){ triangles.push_back(t); }
void addQuadrangle(MQuadrangle *q){ quadrangles.push_back(q); }
void addPolygon(MPolygon *p){ polygons.push_back(p); }
+
+ // get the boundary layer columns
+ BoundaryLayerColumns *getColumns () {return &_columns;}
+
+
};
#endif
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 98dd7bf120..c21bc530ea 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -17,7 +17,6 @@
#include "GRegion.h"
#include "SPoint3.h"
#include "SBoundingBox3d.h"
-#include "boundaryLayersData.h"
template <class scalar> class simpleFunction;
@@ -280,8 +279,6 @@ class GModel
std::vector<GEdge*> bindingsGetEdges();
std::vector<GVertex*> bindingsGetVertices();
- // get the boundary layer columns
- BoundaryLayerColumns *getColumns () {return &_columns;}
// add/remove an entity in the model
void add(GRegion *r) { regions.insert(r); }
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 0181a5ab05..42bb88ce9e 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -727,12 +727,14 @@ int OCC_Internals::getGTagOfOCCSolidByNativePtr(GModel *model, TopoDS_Solid toFi
return 0;
}
+
void OCC_Internals::buildGModel(GModel *model)
{
// building geom vertices
int numv = model->getMaxElementaryNumber(0) + 1;
for(int i = 1; i <= vmap.Extent(); i++){
TopoDS_Vertex vertex = TopoDS::Vertex(vmap(i));
+
if(!getOCCVertexByNativePtr(model, vertex)){
model->add(new OCCVertex(model, numv, vertex));
numv++;
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index b7063057d3..42106a8e5d 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -92,3 +92,32 @@ bool GVertex::isOnSeam(const GFace *gf) const
}
return false;
}
+
+// faces that bound this entity or that this entity bounds.
+std::list<GFace*> GVertex::faces() const
+{
+ std::list<GEdge*>::const_iterator it = l_edges.begin();
+ std::set<GFace*> _f;
+ for ( ; it != l_edges.end() ; ++it){
+ std::list<GFace*> temp = (*it)->faces();
+ _f.insert (temp.begin(), temp.end());
+ }
+ std::list<GFace*> ret;
+ ret.insert (ret.begin(), _f.begin(), _f.end());
+ return ret;
+}
+
+// regions that bound this entity or that this entity bounds.
+std::list<GRegion*> GVertex::regions() const
+{
+ std::list<GFace*> _faces = faces();
+ std::list<GFace*>::const_iterator it = _faces.begin();
+ std::set<GRegion*> _r;
+ for ( ; it != _faces.end() ; ++it){
+ std::list<GRegion*> temp = (*it)->regions();
+ _r.insert (temp.begin(), temp.end());
+ }
+ std::list<GRegion*> ret;
+ ret.insert (ret.begin(), _r.begin(), _r.end());
+ return ret;
+}
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index 5332fa4db8..151a59a2ee 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -45,9 +45,15 @@ class GVertex : public GEntity
void addEdge(GEdge *e);
void delEdge(GEdge *e);
+ // regions that bound this entity or that this entity bounds.
+ virtual std::list<GRegion*> regions() const;
+
// get the edges that this vertex bounds
virtual std::list<GEdge*> edges() const{ return l_edges; }
+ // faces that bound this entity or that this entity bounds.
+ virtual std::list<GFace*> faces() const;
+
// get the dimension of the vertex (0)
virtual int dim() const { return 0; }
diff --git a/Geo/MFace.h b/Geo/MFace.h
index dd16c1c02b..6ad37f7cba 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -131,6 +131,8 @@ struct Equal_Face : public std::binary_function<MFace, MFace, bool> {
struct Less_Face : public std::binary_function<MFace, MFace, bool> {
bool operator()(const MFace &f1, const MFace &f2) const
{
+ if (f1.getNumVertices() != f2.getNumVertices())
+ return f1.getNumVertices() < f2.getNumVertices();
for(int i = 0; i < f1.getNumVertices(); i++) {
if(f1.getSortedVertex(i) < f2.getSortedVertex(i)) return true;
if(f1.getSortedVertex(i) > f2.getSortedVertex(i)) return false;
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index a9db02fd49..7c1a63370c 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -9,11 +9,41 @@
#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 "Field.h"
#include "OS.h"
+const int FANSIZE__ = 4;
+
+
+/*
+ ^ ni
+ |
+ |
+ +-----------------+
+ bi /
+ bj /
+ /\
+ / \ nj
+ / Z
+ +
+*/
+
+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)
{
SVector3 ez (0,0,1);
@@ -93,6 +123,66 @@ void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3])
metric[2] = res[1][1];
}
+const BoundaryLayerData & BoundaryLayerColumns::getColumn(MVertex *v, MFace f)
+{
+ int N = getNbColumns(v) ;
+ if (N == 1) return getColumn(v, 0);
+ if (isOnWedge (v)){
+ GFace *gf = _inverse_classification[f];
+ BoundaryLayerFanWedge3d w = getWedge(v);
+ if (w.isLeft(gf))return getColumn(v, 0);
+ if (w.isRight(gf))return getColumn(v, N-1);
+ Msg::Error("Strange behavior for a wedge");
+ }
+ if (isCorner (v)){
+ GFace *gf = _inverse_classification[f];
+ BoundaryLayerFanCorner3d c = getCorner(v);
+ int k = 0;
+ for (unsigned int i=0;i<c._gf.size();i++){
+ if (c._gf[i] == gf){
+ return getColumn(v, k);
+ }
+ k += (c._fanSize - 1);
+ }
+ }
+ 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)
{
Equal_Edge aaa;
@@ -188,11 +278,11 @@ edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side)
return error2;
}
-/*
-
-
+// FIXME
+static bool pointInFace (GFace *gf, double u, double v) {
+ return true;
+}
-*/
static void treat2Connections (GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2,
double _treshold, BoundaryLayerColumns *_columns,
std::vector<SVector3> &_dirs, bool test = false)
@@ -209,45 +299,83 @@ static void treat2Connections (GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2
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 (test) printf("angle %12.5E\n", 180*angle/M_PI);
+ // if (test)
+ // printf("angle %12.5E\n", 180*angle/M_PI);
if (angle < _treshold /*&& angle > - _treshold*/){
SVector3 x = N1[SIDE]*1.01+N2[SIDE];
x.normalize();
_dirs.push_back(x);
}
else if (angle >= _treshold){
- int fanSize = angle / _treshold;
- // if the angle is greater than PI, than reverse the sense
- double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x());
- double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x());
- double AMAX = std::max(alpha1,alpha2);
- double AMIN = std::min(alpha1,alpha2);
- MEdge ee[2];
- if (alpha1 > alpha2){
- ee[0] = e2;ee[1] = e1;
+
+ if (USEFANS__){
+ int fanSize = FANSIZE__; //angle / _treshold;
+ // if the angle is greater than PI, than reverse the sense
+ double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x());
+ double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x());
+ double AMAX = std::max(alpha1,alpha2);
+ double AMIN = std::min(alpha1,alpha2);
+ MEdge ee[2];
+ if (alpha1 > alpha2){
+ ee[0] = e2;ee[1] = e1;
+ }
+ else {
+ ee[0] = e1;ee[1] = e2;
+ }
+ if ( AMAX - AMIN >= M_PI){
+ double temp = AMAX;
+ AMAX = AMIN + 2*M_PI;
+ AMIN = temp;
+ MEdge eee0 = ee[0];
+ ee[0] = ee[1];ee[1] = eee0;
+ }
+ _columns->addFan (_myVert,ee[0],ee[1],true);
+ for (int i=-1; i<=fanSize; i++){
+ double t = (double)(i+1)/ (fanSize+1);
+ double alpha = t * AMAX + (1.-t)* AMIN;
+ SVector3 x (cos(alpha),sin(alpha),0);
+ x.normalize();
+ _dirs.push_back(x);
+ }
}
else {
- ee[0] = e1;ee[1] = e2;
- }
- if ( AMAX - AMIN >= M_PI){
- double temp = AMAX;
- AMAX = AMIN + 2*M_PI;
- AMIN = temp;
- MEdge eee0 = ee[0];
- ee[0] = ee[1];ee[1] = eee0;
- }
- _columns->addFan (_myVert,ee[0],ee[1],true);
- for (int i=-1; i<=fanSize; i++){
- double t = (double)(i+1)/ (fanSize+1);
- double alpha = t * AMAX + (1.-t)* AMIN;
- SVector3 x (cos(alpha),sin(alpha),0);
- x.normalize();
- _dirs.push_back(x);
+ _dirs.push_back(N1[SIDE]);
+ _dirs.push_back(N2[SIDE]);
}
}
}
}
}
+// static void treat2Connections (GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2,
+// double _treshold, BoundaryLayerColumns *_columns,
+// std::vector<SVector3> &_dirs, bool test = false)
+// {
+// std::vector<SVector3> N1,N2;
+// for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
+// _columns->_normals.lower_bound(e1);
+// itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second);
+// for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
+// _columns->_normals.lower_bound(e2);
+// itm != _columns->_normals.upper_bound(e2); ++itm) N2.push_back(itm->second);
+// if (test) printf("%d %d\n", (int)N1.size(), (int)N2.size());
+// 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 (test)
+// // printf("angle %12.5E\n", 180*angle/M_PI);
+// if (angle < _treshold /*&& angle > - _treshold*/){
+// SVector3 x = N1[SIDE]*1.01+N2[SIDE];
+// x.normalize();
+// _dirs.push_back(x);
+// }
+// else if (angle >= _treshold){
+// _dirs.push_back(N1[SIDE]);
+// _dirs.push_back(N2[SIDE]);
+// }
+// }
+// }
+// }
static void treat3Connections (GFace *gf, MVertex *_myVert, MEdge &e1,
MEdge &e2, MEdge &e3, double _treshold,
@@ -299,40 +427,50 @@ static void treat3Connections (GFace *gf, MVertex *_myVert, MEdge &e1,
_dirs.push_back(x2);
}
-
-bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
-{
-#if !defined(HAVE_ANN) || !defined(HAVE_MESH)
- return false;
-#else
- //// GET THE FIELD THAT DEFINES THE DISTANCE FUNCTION
- FieldManager *fields = gf->model()->getFields();
+BoundaryLayerField* getBLField (GModel *gm) {
+ FieldManager *fields = gm->getFields();
if(fields->getBoundaryLayerField() <= 0){
return 0;
}
Field *bl_field = fields->get(fields->getBoundaryLayerField());
- BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField*> (bl_field);
-
- if (!blf)return false;
-
- double _treshold = blf->fan_angle * M_PI / 180 ;
-
- // BoundaryLayerColumns * _columns = new BoundaryLayerColumns;
+ return dynamic_cast<BoundaryLayerField*> (bl_field);
+}
- // assume a field that i) gives us the closest point on one of the BL components
- // ii) Gives us mesh sizes in the 3 directions.
+static bool isEdgeOfFaceBL (GFace *gf,
+ GEdge *ge,
+ BoundaryLayerField *blf){
+ if (blf->isEdgeBL (ge->tag()))return true;
+ /*
+ std::list<GFace*> faces = ge->faces();
+ for (std::list<GFace*>::iterator it = faces.begin();
+ it != faces.end() ; ++it){
+ if ((*it) == gf)return false;
+ }
+ for (std::list<GFace*>::iterator it = faces.begin();
+ it != faces.end() ; ++it){
+ if (blf->isFaceBL ((*it)->tag()))return true;
+ }
+ */
+ return false;
+}
- ///// build vertex to vertex connexions
+static void getEdgesData (GFace *gf,
+ BoundaryLayerField *blf,
+ BoundaryLayerColumns *_columns,
+ std::set<MVertex*> &_vertices,
+ std::set<MEdge,Less_Edge> &allEdges,
+ std::multimap<MVertex*,MVertex*> &tangents)
+{
+ // get all model edges
std::list<GEdge*> edges = gf->edges();
std::list<GEdge*> embedded_edges = gf->embeddedEdges();
edges.insert(edges.begin(), embedded_edges.begin(),embedded_edges.end());
+
+ // iterate on model edges
std::list<GEdge*>::iterator ite = edges.begin();
- std::set<MVertex*> _vertices;
- std::set<MEdge,Less_Edge> allEdges;
- std::multimap<MVertex*,MVertex*> tangents;
while(ite != edges.end()){
// check if this edge generates a boundary layer
- if (blf->isEdgeBL ((*ite)->tag())){
+ if (isEdgeOfFaceBL (gf,*ite,blf)){
for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
MVertex *v1 = (*ite)->lines[i]->getVertex(0);
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
@@ -343,7 +481,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
_vertices.insert(v2);
}
}
- else {
+ else {
MVertex *v1 = (*ite)->lines[0]->getVertex(0);
MVertex *v2 = (*ite)->lines[0]->getVertex(1);
MVertex *v3 = (*ite)->lines[(*ite)->lines.size() - 1]->getVertex(1);
@@ -353,8 +491,13 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
}
++ite;
}
- if (!_vertices.size())return false;
+}
+static void getNormals (GFace *gf,
+ BoundaryLayerField *blf,
+ BoundaryLayerColumns *_columns,
+ std::set<MEdge,Less_Edge> &allEdges)
+{
// assume that the initial mesh has been created i.e. that there exist
// triangles inside the domain. Triangles are used to define
// exterior normals
@@ -365,13 +508,13 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
MVertex *v2 = gf->triangles[i]->getVertex(2);
reparamMeshEdgeOnFace(v0, v1, gf, p0, p1);
reparamMeshEdgeOnFace(v0, v2, gf, p0, p2);
-
+
MEdge me01(v0,v1);
if (allEdges.find(me01) != allEdges.end()){
SVector3 v01 = interiorNormal (p0,p1,p2);
_columns->_normals.insert(std::make_pair(me01,v01));
}
-
+
MEdge me02(v0,v2);
if (allEdges.find(me02) != allEdges.end()){
SVector3 v02 = interiorNormal (p0,p2,p1);
@@ -384,6 +527,69 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
_columns->_normals.insert(std::make_pair(me21,v21));
}
}
+}
+
+void addColumnAtTheEndOfTheBL (GEdge *ge,
+ GVertex *gv,
+ BoundaryLayerColumns* _columns,
+ BoundaryLayerField *blf)
+{
+ // printf("coucou %d\n",ge->tag());
+ if (!blf->isEdgeBL(ge->tag()))
+ {
+ GVertex *g0 = ge->getBeginVertex();
+ GVertex *g1 = ge->getEndVertex();
+ // printf("coucou 2 %d %d vs %d\n",g0->tag(),g1->tag(),gv->tag());
+ MVertex * v0 = g0->mesh_vertices[0];
+ MVertex * v1 = g1->mesh_vertices[0];
+ std::vector<MVertex*> invert;
+ 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));
+ }
+ 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);
+ }
+ else if (g1 == gv){
+ _columns->addColumn(t*-1.0, v1,invert,_metrics);
+ }
+ }
+}
+
+
+bool buildAdditionalPoints2D (GFace *gf)
+{
+#if !defined(HAVE_ANN) || !defined(HAVE_MESH)
+ return false;
+#else
+ BoundaryLayerColumns *_columns = gf->getColumns();
+
+ _columns->_normals.clear();
+ _columns->_non_manifold_edges.clear();
+ _columns->_data.clear();
+
+ //// GET THE FIELD THAT DEFINES THE DISTANCE FUNCTION
+ BoundaryLayerField *blf = getBLField (gf->model());
+
+ if (!blf)return false;
+
+ blf->setupFor2d(gf->tag());
+
+ double _treshold = blf->fan_angle * M_PI / 180 ;
+
+ std::set<MVertex*> _vertices;
+ std::set<MEdge,Less_Edge> allEdges;
+ std::multimap<MVertex*,MVertex*> tangents;
+
+ getEdgesData ( gf, blf, _columns, _vertices , allEdges , tangents );
+
+ if (!_vertices.size())return false;
+
+ getNormals ( gf, blf, _columns, allEdges);
// for all boundry points
for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
@@ -433,29 +639,29 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
for (std::multimap<MVertex*,MVertex*>::iterator itm =
tangents.lower_bound(*it);
itm != tangents.upper_bound(*it); ++itm) Ts.push_back(itm->second);
- // end of the BL
+ // end of the BL --> let's add a column that correspond to the
+ // model edge that lies after the end of teh BL
if (Ts.size() == 1){
- MEdge e2 (*it,Ts[0]);
- SPoint2 p0,p1;
- reparamMeshEdgeOnFace(*it,Ts[0], gf, p0, p1);
- dirEndOfBL = SVector3(p1.x()-p0.x(),p1.y()-p0.y(),0);
- dirEndOfBL.normalize();
- endOfTheBL = true;
- _columns->_normals.insert (std::make_pair(e2,dirEndOfBL));
- treat2Connections (gf, *it,e1,e2, _treshold, _columns, _dirs, true);
+ // printf("HERE WE ARE IN FACE %d %d\n",gf->tag(),Ts.size());
+ // printf("Classif dim %d %d\n",(*it)->onWhat()->dim(),Ts[0]->onWhat()->dim());
+ GEdge *ge = dynamic_cast<GEdge*>(Ts[0]->onWhat());
+ GVertex *gv = dynamic_cast<GVertex*>((*it)->onWhat());
+ if (ge && gv){
+ addColumnAtTheEndOfTheBL (ge,gv,_columns,blf);
+ }
}
else {
- Msg::Error("Impossible BL Configuration");
+ Msg::Error("Impossible BL Configuration -- One Edge -- Tscp.size() = %d",Ts.size());
}
}
else if (N1.size() == 2){
// printf("%g %g --> %g %g \n",e1.getVertex(0)->x(),e1.getVertex(0)->y(),
// e1.getVertex(1)->x(),e1.getVertex(1)->y());
- // printf("N1.size = %d %g %g %g %g\n",N1.size(),N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
+ // printf("N1.size = %d %g %g %g %g\n",N1.size(),N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
SPoint2 p0,p1;
reparamMeshEdgeOnFace(*it,_connections[0], gf, p0, p1);
- int fanSize = M_PI / _treshold;
+ int fanSize = FANSIZE__;//M_PI / _treshold;
double alpha1 = atan2(N1[0].y(),N1[0].x());
double alpha2 = atan2(N1[1].y(),N1[1].x());
double alpha3 = atan2(p1.y()-p0.y(),p1.x()-p0.x());
@@ -471,7 +677,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
AMIN = temp;
}
_columns->addFan (*it,e1,e1,true);
- // printf("%g %g --> %g %g\n",N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
+ // printf("%g %g --> %g %g\n",N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
for (int i=-1; i<=fanSize; i++){
double t = (double)(i+1)/ (fanSize+1);
double alpha = t * AMAX + (1.-t)* AMIN;
@@ -495,12 +701,13 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
// = e * (dX/du n_u + dX/dv n_v) //
// < ------------------------------- > //
- if (endOfTheBL){
+ /* if (endOfTheBL){
printf("%g %g %d %d %g\n", (*it)->x(), (*it)->y(), DIR, (int)_dirs.size(),
dot(n, dirEndOfBL));
}
+ */
if (endOfTheBL && dot(n,dirEndOfBL) > .99){
- printf( "coucou c'est moi\n");
+ // printf( "coucou c'est moi\n");
}
else {
MVertex *current = *it;
@@ -518,7 +725,7 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
SMetric3 m;
double metric[3];
double l;
- (*bl_field)(current->x(),current->y(), current->z(), m, current->onWhat());
+ (*blf)(current->x(),current->y(), current->z(), m, current->onWhat());
if (!catt){
catt = blf->current_closest;
_close = blf->_closest_point;
@@ -543,20 +750,19 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
_current->bl_data = new MVertexBoundaryLayerData;
current = _current;
_column.push_back(current);
- // printf("pnew %g %g new point %g %g n %g %g\n",pnew.x(),pnew.y(),gp.x(),gp.y(),n.x(),n.y());
_metrics.push_back(m);
- // const double l = n[0]*m(0,0) +;
if ((int)_column.size() > nbCol) break; // FIXME
p = pnew;
}
- // if (_dirs.size() > 1)printf("adding column with %d nodes\n",_column.size());
_columns->addColumn(n,*it, _column, _metrics);
}
}
}
// DEBUG STUFF
- FILE *f = Fopen ("test.pos","w");
+ char name[256];
+ sprintf(name,"points_face_%d.pos",gf->tag());
+ FILE *f = Fopen (name,"w");
fprintf(f,"View \"\" {\n");
for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
MVertex *v = *it;
@@ -577,3 +783,322 @@ bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
#endif
}
+double angle_0_180 (SVector3 &n1, SVector3 &n2){
+ double cosa = dot(n1,n2)/(n1.norm()*n2.norm());
+ // printf("NORMS = %12.5E %12.5E cosa = %22.15E a = %22.15E\n",n1.norm(),n2.norm(),cosa, acos(cosa));
+ if (cosa > 1.)cosa = 1.0;
+ if (cosa < -1.)cosa = -1.0;
+ return acos(cosa);
+}
+
+
+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 void 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,
+ double _treshold)
+{
+ SVector3 n[256];
+ SPoint3 c[256];
+ int count = 0;
+ GFace *gfs[256];
+
+
+ // generate datas per face;
+
+ 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){
+
+ n[count] += _normals[itm->second->getFace(0)];
+ c[count] = itm->second->getFace(0).barycenter();
+ }
+ gfs[count] = *it;
+ n[count].normalize();
+ count ++;
+ }
+
+ // we throw a column per set of faces that have normals that are sufficiently close
+
+ // printf("vertex %d %d faces\n",myV->getNum(),count);
+
+ std::set<int> done;
+ for (int i=0;i<count;i++){
+ if (done.find(i) == done.end()){
+ SVector3 n1 = n[i];
+ SPoint3 c1 = c[i];
+ SVector3 avg = n1;
+ std::vector<GFace*> joint;
+ joint.push_back(gfs[i]);
+ for (int j=i;j<count;j++){
+ if (done.find(j) == done.end()){
+ SVector3 n2 = n[j];
+ SPoint3 c2 = c[j];
+ double angle = angle_0_180 (n1,n2);
+ double sign = dot((n1-n2),(c1-c2));
+ if (!(angle > _treshold && sign > 0)){
+ joint.push_back(gfs[j]);
+ avg += n2;
+ done.insert(j);
+ }
+ }
+ }
+ if (joint.size()){
+ std::vector<MVertex*> _column;
+ std::vector<SMetric3> _metrics;
+ avg.normalize();
+ createBLPointsInDir (gr,myV,blf,avg,_column,_metrics);
+ _columns->addColumn(avg,myV, _column, _metrics, joint);
+ }
+ // printf("adding one column for %d faces\n",joint.size());
+ }
+ }
+}
+
+static bool createWedgeBetweenTwoFaces (bool testOnly,
+ MVertex *myV,
+ GFace *gf1, GFace *gf2,
+ std::multimap<GFace*,MTriangle*> & _faces,
+ std::map<MFace,SVector3,Less_Face> &_normals,
+ double _treshold,
+ std::vector<SVector3> &shoot)
+{
+ SVector3 n1,n2;
+ SPoint3 c1,c2;
+ for (std::multimap<GFace*,MTriangle*>::iterator itm =
+ _faces.lower_bound(gf1);
+ itm != _faces.upper_bound(gf1); ++itm){
+ n1 += _normals[itm->second->getFace(0)];
+ c1 = itm->second->getFace(0).barycenter();
+ }
+ for (std::multimap<GFace*,MTriangle*>::iterator itm =
+ _faces.lower_bound(gf2);
+ itm != _faces.upper_bound(gf2); ++itm){
+ n2 += _normals[itm->second->getFace(0)];
+ c2 = itm->second->getFace(0).barycenter();
+ }
+ n1.normalize();
+ n2.normalize();
+ // FIXME WRONG FOR INTERNAL CORNERS !!!
+ double angle = angle_0_180 (n1,n2);
+ double sign = dot((n1-n2),(c1-c2));
+ if (angle > _treshold && sign > 0){
+ if(testOnly)return true;
+ int fanSize = FANSIZE__; //angle / _treshold;
+ 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;
+ }
+ else {
+ if(testOnly)return false;
+ SVector3 n = n1+n2;
+ n.normalize();
+ shoot.push_back(n);
+ return false;
+ }
+}
+
+
+BoundaryLayerColumns * buildAdditionalPoints3D (GRegion *gr)
+{
+#if !defined(HAVE_ANN)
+ return 0;
+#else
+ BoundaryLayerField *blf = getBLField (gr->model());
+
+ if (!blf)return 0;
+
+ blf->setupFor3d();
+
+ double _treshold = blf->fan_angle * M_PI / 180 ;
+
+ BoundaryLayerColumns * _columns = new BoundaryLayerColumns;
+
+ 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;
+ 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));
+ _normals [(*itf)->triangles[i]->getFace(0)] = SVector3(0,0,0);
+ }
+ }
+ }
+ }
+ ++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));
+ }
+ }
+ }
+
+ // for all boundry points
+ for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
+ std::vector<MTriangle*> _connections;
+ std::vector<SVector3> _dirs, _allDirs;
+ std::list<GFace*> faces = (*it)->onWhat()->faces();
+
+ std::multimap<GFace*,MTriangle*> _faces;
+ std::set<GFace*> _allGFaces;
+ for (std::multimap<MVertex*,MTriangle*>::iterator itm =
+ _columns->_non_manifold_faces.lower_bound(*it);
+ itm != _columns->_non_manifold_faces.upper_bound(*it); ++itm){
+ _connections.push_back (itm->second);
+ _allDirs.push_back (_normals[itm->second->getFace(0)]);
+ GFace *gf = _gfaces[itm->second];
+ _faces.insert(std::make_pair(gf,itm->second));
+ _allGFaces.insert(gf);
+ }
+
+ bool onSymmetryPlane = 0;
+ if ((*it)->onWhat()->dim() != 2){
+ std::list<GFace*> faces = (*it)->onWhat()->faces();
+ if (faces.size() != _allGFaces.size()){
+ onSymmetryPlane = true;
+ }
+ }
+
+ if (onSymmetryPlane){
+ for ( std::list<GFace*>::iterator itf = faces.begin(); itf!= faces.end() ; ++itf){
+ BoundaryLayerColumns* _face_columns = (*itf)->getColumns();
+ int N = _face_columns->getNbColumns(*it);
+ if (N == 1){
+ std::vector<GFace*> _joint;
+ _joint.insert(_joint.begin(),_allGFaces.begin(),_allGFaces.end());
+ const BoundaryLayerData & c = _face_columns->getColumn(*it,0);
+ _columns->addColumn(_allDirs[0],*it, c._column, c._metrics, _joint);
+ }
+ else if (N > 1){
+ Msg::Error ("Impossible connexion between face and region BLs");
+ }
+ }
+ }
+ else
+ createColumnsBetweenFaces (gr,*it,blf,_columns,_allGFaces,_faces,_normals,_treshold);
+
+ }
+
+ // DEBUG STUFF
+
+ FILE *f = fopen ("test3D.pos","w");
+ 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 _columns;
+#endif
+ return 0;
+}
diff --git a/Geo/boundaryLayersData.h b/Geo/boundaryLayersData.h
index fce4b26f34..4da1d4ed2e 100644
--- a/Geo/boundaryLayersData.h
+++ b/Geo/boundaryLayersData.h
@@ -9,6 +9,7 @@
#include "SVector3.h"
#include "STensor3.h"
#include "MEdge.h"
+#include "MFace.h"
#include <map>
#include <set>
@@ -16,17 +17,27 @@ class Field;
class GFace;
class GRegion;
class MTriangle;
+class BoundaryLayerField;
+
+const int USEFANS__ = 1;
+
struct BoundaryLayerData
{
SVector3 _n;
std::vector<MVertex*> _column;
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){}
};
struct BoundaryLayerFan
@@ -37,6 +48,47 @@ struct BoundaryLayerFan
: _e1(e1),_e2(e2) , sense (s){}
};
+// wedge between 2 sets of continuous faces
+struct BoundaryLayerFanWedge3d
+{
+ std::vector<GFace*> _gf1, _gf2;
+ BoundaryLayerFanWedge3d(GFace *gf1=0, GFace *gf2=0)
+ {
+ if(gf1)_gf1.push_back(gf1);
+ if(gf2)_gf2.push_back(gf2);
+ }
+ BoundaryLayerFanWedge3d(std::vector<GFace*> &gf1,
+ std::vector<GFace*> &gf2):_gf1(gf1),_gf2(gf2)
+ {
+ }
+ bool isLeft (const GFace *gf) const {
+ for (unsigned int i=0;i<_gf1.size();i++)if (_gf1[i] == gf)return true;
+ return false;
+ }
+ bool isRight (const GFace *gf) const{
+ for (unsigned int i=0;i<_gf2.size();i++)if (_gf2[i] == gf)return true;
+ return false;
+ }
+ bool isLeft (const std::vector<GFace *> &gf) const {
+ for (unsigned int i=0;i<gf.size();i++)if (isLeft(gf[i]))return true;
+ return false;
+ }
+ bool isRight (const std::vector<GFace *> &gf) const {
+ for (unsigned int i=0;i<gf.size();i++)if (isRight(gf[i]))return true;
+ return false;
+ }
+
+};
+
+struct BoundaryLayerFanCorner3d
+{
+ int _fanSize;
+ std::vector<GFace *> _gf;
+ BoundaryLayerFanCorner3d() : _fanSize(0){}
+ BoundaryLayerFanCorner3d(int fanSize ,std::vector<GFace *> &gf)
+ : _fanSize(fanSize), _gf(gf){}
+};
+
struct edgeColumn {
const BoundaryLayerData &_c1, &_c2;
@@ -44,17 +96,35 @@ 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::multimap<MVertex*, BoundaryLayerData> _data;
std::map<MVertex*, BoundaryLayerFan> _fans;
+ std::map<MVertex*, BoundaryLayerFanWedge3d> _wedges;
+ std::map<MVertex*, BoundaryLayerFanCorner3d> _corners;
public:
+ std::map<MFace, GFace*, Less_Face> _inverse_classification;
+ std::multimap<MVertex*, BoundaryLayerData> _data;
size_t size () const {return _data.size();}
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;
iter begin() { return _data.begin(); }
iter end() { return _data.end(); }
iterf beginf() { return _fans.begin(); }
@@ -66,10 +136,48 @@ public:
{
_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)
+ {
+ _data.insert (std::make_pair(v,BoundaryLayerData(dir, _column,_metrics,_joint)));
+ }
inline void addFan(MVertex *v, MEdge e1, MEdge e2, bool s)
{
_fans.insert(std::make_pair(v,BoundaryLayerFan(e1,e2,s)));
}
+ inline void addWedge(MVertex *v, GFace *gf1, GFace *gf2)
+ {
+ _wedges.insert(std::make_pair(v,BoundaryLayerFanWedge3d(gf1,gf2)));
+ }
+ inline void addWedge(MVertex *v, std::vector<GFace *>&gf1, std::vector<GFace *> &gf2)
+ {
+ _wedges.insert(std::make_pair(v,BoundaryLayerFanWedge3d(gf1,gf2)));
+ }
+ inline void addCorner(MVertex *v, int fanSize, std::vector<GFace *> &gfs)
+ {
+ _corners.insert(std::make_pair(v,BoundaryLayerFanCorner3d(fanSize, gfs)));
+ }
+ inline bool isCorner (MVertex* v) const{
+ return _corners.find(v) != _corners.end();
+ }
+ inline bool isOnWedge (MVertex* v) const{
+ return _wedges.find(v) != _wedges.end();
+ }
+
+ inline BoundaryLayerFanWedge3d getWedge (MVertex* v) {
+ std::map<MVertex*, BoundaryLayerFanWedge3d>::iterator it = _wedges.find(v);
+ return it->second;
+ }
+ inline BoundaryLayerFanCorner3d getCorner (MVertex* v) {
+ std::map<MVertex*, BoundaryLayerFanCorner3d>::iterator it = _corners.find(v);
+ return it->second;
+ }
+
+
+ const BoundaryLayerData &getColumn(MVertex *v, MFace f);
+
inline const BoundaryLayerData &getColumn(MVertex *v, MEdge e)
{
std::map<MVertex*,BoundaryLayerFan>::const_iterator it = _fans.find(v);
@@ -87,6 +195,7 @@ public:
return error;
}
edgeColumn getColumns(MVertex *v1, MVertex *v2 , int side);
+ faceColumn getColumns(GFace *gf, MVertex *v1, MVertex *v2 , MVertex* v3, int side);
inline int getNbColumns(MVertex *v) { return _data.count(v); }
inline const BoundaryLayerData &getColumn(MVertex *v, int iColumn)
{
@@ -101,8 +210,9 @@ public:
void filterPoints();
};
-bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns * ) ;
-bool buildAdditionalPoints3D (GRegion *gr, BoundaryLayerColumns * ) ;
+bool buildAdditionalPoints2D (GFace *gf ) ;
+BoundaryLayerColumns * buildAdditionalPoints3D (GRegion *gr) ;
void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]);
+BoundaryLayerField* getBLField (GModel *gm);
#endif
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index fd4bc45550..bd9cde08ae 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -13,7 +13,7 @@ set(SRC
meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp
meshGFaceLloyd.cpp meshGFaceOptimize.cpp
meshGFaceQuadrilateralize.cpp
- meshGRegion.cpp meshGRegionBoundaryLayers.cpp
+ meshGRegion.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
meshGRegionLocalMeshMod.cpp meshGRegionMMG3D.cpp
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 696a6d247f..dd11918269 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1670,11 +1670,24 @@ class AttractorField : public Field
it != faces_id.end(); ++it) {
GFace *f = GModel::current()->getFaceByTag(*it);
if (f){
- SBoundingBox3d bb = f->bounds();
- SVector3 dd = bb.max() - bb.min();
- double maxDist = dd.norm() / n_nodes_by_edge ;
- f->fillPointCloud(maxDist, &points, &uvpoints);
- offset.push_back(points.size());
+
+ if (f->mesh_vertices.size()){
+ for (unsigned int i=0;i<f->mesh_vertices.size();i++){
+ MVertex *v = f->mesh_vertices[i];
+ double uu,vv;
+ v->getParameter(0,uu);
+ v->getParameter(1,vv);
+ points.push_back(SPoint3(v->x(),v->y(),v->z()));
+ uvpoints.push_back(SPoint2(uu,vv));
+ }
+ }
+ else {
+ SBoundingBox3d bb = f->bounds();
+ SVector3 dd = bb.max() - bb.min();
+ double maxDist = dd.norm() / n_nodes_by_edge ;
+ f->fillPointCloud(maxDist, &points, &uvpoints);
+ offset.push_back(points.size());
+ }
}
}
@@ -1823,6 +1836,85 @@ BoundaryLayerField::BoundaryLayerField()
(thickness, "Maximal thickness of the boundary layer");
}
+void BoundaryLayerField::removeAttractors()
+{
+ for (std::list<AttractorField *>::iterator it = _att_fields.begin();
+ it != _att_fields.end() ; ++it) delete *it;
+ _att_fields.clear();
+ update_needed = true;
+}
+
+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();
+ 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;
+ 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;
+ }
+ }
+ }
+ 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,edges_id.size());
+
+ removeAttractors();
+ }
+}
+
+void BoundaryLayerField::setupFor3d()
+{
+ faces_id = faces_id_saved;
+ removeAttractors();
+}
+
+
double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge)
{
diff --git a/Mesh/Field.h b/Mesh/Field.h
index 9cc0992dbc..7914f3c8cd 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -141,6 +141,7 @@ 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;
void operator() (AttractorField *cc, double dist, double x, double y, double z,
SMetric3 &metr, GEntity *ge);
public:
@@ -153,12 +154,16 @@ class BoundaryLayerField : public Field {
virtual const char *getName();
virtual std::string getDescription();
BoundaryLayerField();
+ ~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();}
bool isVertexBL (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 setupFor2d(int iF);
+ void setupFor3d();
+ void removeAttractors();
};
#endif
class FieldOptionString : public FieldOption
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 1e777e3104..5561081d18 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -4,6 +4,7 @@
// bugs and problems to the public mailing list <gmsh@geuz.org>.
#include <stdlib.h>
+#include <stack>
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "Numeric.h"
@@ -575,8 +576,36 @@ static void Mesh2D(GModel *m)
static void FindConnectedRegions(std::vector<GRegion*> &delaunay,
std::vector<std::vector<GRegion*> > &connected)
{
- // FIXME: need to split region vector into connected components here!
- connected.push_back(delaunay);
+ const unsigned int nbVolumes = delaunay.size();
+ if (!nbVolumes)return;
+ while (delaunay.size()){
+ std::set<GRegion*> oneDomain;
+ std::stack<GRegion*> _stack;
+ GRegion *r = delaunay[0];
+ _stack.push(r);
+ while(!_stack.empty()){
+ r = _stack.top();
+ _stack.pop();
+ oneDomain.insert(r);
+ std::list<GFace*> faces = r->faces();
+ for (std::list<GFace*> :: iterator it = faces.begin(); it != faces.end() ; ++it){
+ GFace *gf = *it;
+ GRegion *other = gf->getRegion(0) == r ? gf->getRegion(1) : gf->getRegion(0);
+ if (other != 0 && oneDomain.find(other) == oneDomain.end())
+ _stack.push (other);
+ }
+ }
+ std::vector<GRegion*> temp1,temp2;
+ for (unsigned int i=0;i<delaunay.size();i++){
+ r = delaunay[i];
+ if (oneDomain.find(r) == oneDomain.end())temp1.push_back(r);
+ else temp2.push_back(r);
+ }
+ connected.push_back(temp2);
+ delaunay=temp1;
+ }
+ Msg::Info("Delaunay Meshing %d volumes with %d connected components",
+ nbVolumes,connected.size());
}
static void Mesh3D(GModel *m)
@@ -606,6 +635,7 @@ static void Mesh3D(GModel *m)
FindConnectedRegions(delaunay, connected);
// remove quads elements for volumes that are recombined
+ // pragma OMP ICI ??
for(unsigned int i = 0; i < connected.size(); i++){
for(unsigned j=0;j<connected[i].size();j++){
GRegion *gr = connected[i][j];
@@ -617,6 +647,7 @@ static void Mesh3D(GModel *m)
}
}
+ // pragma OMP ICI ??
for(unsigned int i = 0; i < connected.size(); i++){
MeshDelaunayVolume(connected[i]);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index b790eb7ce6..c82d71fef2 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -492,31 +492,6 @@ void meshGEdge::operator() (GEdge *ge)
v0->z() = beg_p.z();
}
-#if defined(HAVE_ANN)
- if (blf && !blf->isEdgeBL(ge->tag()))
- {
- GVertex *g0 = ge->getBeginVertex();
- GVertex *g1 = ge->getEndVertex();
- BoundaryLayerColumns* _columns = ge->model()->getColumns();
- MVertex * v0 = g0->mesh_vertices[0];
- MVertex * v1 = g1->mesh_vertices[0];
- std::vector<MVertex*> invert;
- std::vector<SMetric3> _metrics;
- for(unsigned int i = 0; i < mesh_vertices.size() ; i++)
- {
- invert.push_back(mesh_vertices[mesh_vertices.size() - i - 1]);
- _metrics.push_back(SMetric3(1.0));
- }
- SVector3 t (v1->x()-v0->x(), v1->y()-v0->y(),v1->z()-v0->z());
- t.normalize();
- if (blf->isVertexBL(g0->tag())){
- _columns->addColumn(t, v0, mesh_vertices,_metrics);
- }
- if (blf->isVertexBL(g1->tag())){
- _columns->addColumn(t*-1.0, v1,invert,_metrics);
- }
- }
-#endif
ge->meshStatistics.status = GEdge::DONE;
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 0ff6ce8149..f82360f2ac 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -630,8 +630,8 @@ void filterOverlappingElements(int dim, std::vector<MElement*> &e,
void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
{
- BoundaryLayerColumns* _columns = gf->model()->getColumns();
- if (!buildAdditionalPoints2D (gf, _columns))return;
+ if (!buildAdditionalPoints2D (gf))return;
+ BoundaryLayerColumns* _columns = gf->getColumns();
std::set<MEdge,Less_Edge> bedges;
@@ -687,41 +687,43 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
++ite;
}
- for (BoundaryLayerColumns::iterf itf = _columns->beginf();
- itf != _columns->endf() ; ++itf){
- MVertex *v = itf->first;
- int nbCol = _columns->getNbColumns(v);
-
- for (int i=0;i<nbCol-1;i++){
- const BoundaryLayerData & c1 = _columns->getColumn(v,i);
- const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
- int N = std::min(c1._column.size(),c2._column.size());
- for (int l=0;l < N ;++l){
- MVertex *v11,*v12,*v21,*v22;
- v21 = c1._column[l];
- v22 = c2._column[l];
- if (l == 0){
- v11 = v;
- v12 = v;
- }
- else {
- v11 = c1._column[l-1];
- v12 = c2._column[l-1];
- }
- if (v11 != v12){
- blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
- fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
- v11->x(),v11->y(),v11->z(),
+ if (USEFANS__){
+ for (BoundaryLayerColumns::iterf itf = _columns->beginf();
+ itf != _columns->endf() ; ++itf){
+ MVertex *v = itf->first;
+ int nbCol = _columns->getNbColumns(v);
+
+ for (int i=0;i<nbCol-1;i++){
+ const BoundaryLayerData & c1 = _columns->getColumn(v,i);
+ const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
+ int N = std::min(c1._column.size(),c2._column.size());
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v21,*v22;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ if (l == 0){
+ v11 = v;
+ v12 = v;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ }
+ if (v11 != v12){
+ blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
+ fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
+ v11->x(),v11->y(),v11->z(),
v12->x(),v12->y(),v12->z(),
- v22->x(),v22->y(),v22->z(),
- v21->x(),v21->y(),v21->z());
- }
- else {
- blTris.push_back(new MTriangle(v,v22,v21));
- fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
- v->x(),v->y(),v->z(),
- v22->x(),v22->y(),v22->z(),
- v21->x(),v21->y(),v21->z());
+ v22->x(),v22->y(),v22->z(),
+ v21->x(),v21->y(),v21->z());
+ }
+ else {
+ blTris.push_back(new MTriangle(v,v22,v21));
+ fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
+ v->x(),v->y(),v->z(),
+ v22->x(),v22->y(),v22->z(),
+ v21->x(),v21->y(),v21->z());
+ }
}
}
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index aca7e4f7ff..19b59aa852 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1232,13 +1232,13 @@ void bowyerWatsonFrontal(GFace *gf,
int ITERATION = 0;
while (1){
++ITERATION;
- if(ITERATION % 10== 0 && CTX::instance()->mesh.saveAll){
- char name[245];
- sprintf(name,"delFrontal_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
- _printTris (name, AllTris.begin(), AllTris.end(), DATA,true);
- sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
- _printTris (name, ActiveTris.begin(), ActiveTris.end(), DATA,true);
- }
+ // if(ITERATION % 10== 0 && CTX::instance()->mesh.saveAll){
+ // char name[245];
+ // sprintf(name,"delFrontal_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
+ // _printTris (name, AllTris.begin(), AllTris.end(), DATA,true);
+ // sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ // _printTris (name, ActiveTris.begin(), ActiveTris.end(), DATA,true);
+ // }
/* if(ITER % 100== 0){
char name[245];
sprintf(name,"delfr2d%d-ITER%d.pos",gf->tag(),ITER);
@@ -1247,10 +1247,10 @@ void bowyerWatsonFrontal(GFace *gf,
// _printTris (name, ActiveTris, Us,Vs,false);
}
*/
+ // printf("active_tris.size = %d\n",ActiveTris.size());
if (!ActiveTris.size())break;
MTri3 *worst = (*ActiveTris.begin());
ActiveTris.erase(ActiveTris.begin());
- // printf("active_tris.size = %d\n",ActiveTris.size());
if (!worst->isDeleted() && isActive(worst, LIMIT_, active_edge) &&
worst->getRadius() > LIMIT_){
@@ -1285,7 +1285,9 @@ void bowyerWatsonFrontal(GFace *gf,
if(fields->getBoundaryLayerField() > 0){
Field *bl_field = fields->get(fields->getBoundaryLayerField());
blf = dynamic_cast<BoundaryLayerField*> (bl_field);
- if (blf && !blf->iRecombine)quadsToTriangles(gf,10000);
+ if (blf && !blf->iRecombine){
+ quadsToTriangles(gf,10000);
+ }
}
}
#endif
@@ -1467,13 +1469,13 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
int max_layers = quad ? 10000 : 4;
while (1){
ITERATION ++;
- if(ITERATION % 1== 0 && CTX::instance()->mesh.saveAll){
- char name[245];
- sprintf(name,"delInfinite_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
- _printTris (name, AllTris.begin(), AllTris.end(), DATA,true);
- sprintf(name,"delInfinite_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
- _printTris (name, ActiveTris.begin(), ActiveTris.end(),DATA,true);
- }
+ // if(ITERATION % 1== 0 && CTX::instance()->mesh.saveAll){
+ // char name[245];
+ // sprintf(name,"delInfinite_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
+ // _printTris (name, AllTris.begin(), AllTris.end(), DATA,true);
+ // sprintf(name,"delInfinite_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ // _printTris (name, ActiveTris.begin(), ActiveTris.end(),DATA,true);
+ // }
std::set<MTri3*,compareTri3Ptr> ActiveTrisNotInFront;
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index a07c34a0e2..ee07a19ad3 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -27,6 +27,7 @@ struct bidimMeshData
std::map<MVertex* , MVertex*>* equivalence;
std::map<MVertex*, SPoint2> * parametricCoordinates;
std::set<MEdge,Less_Edge> internalEdges; // embedded edges
+ // std::set<MVertex*> internalVertices; // embedded vertices
inline void addVertex (MVertex* mv, double u, double v, double size, double sizeBGM){
int index = Us.size();
if (mv->onWhat()->dim() == 2)mv->setIndex(index);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 4acba618ab..30181680c8 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -3516,6 +3516,7 @@ void recombineIntoQuads(GFace *gf,
void quadsToTriangles(GFace *gf, double minqual)
{
+
std::vector<MQuadrangle*> qds;
for (unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *q = gf->quadrangles[i];
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 57a631fbdf..00f641417d 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -8,6 +8,7 @@
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "meshGRegion.h"
+#include "meshGFace.h"
#include "meshGFaceOptimize.h"
#include "boundaryLayersData.h"
#include "meshGRegionDelaunayInsertion.h"
@@ -21,6 +22,8 @@
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MPyramid.h"
+#include "MPrism.h"
+#include "MHexahedron.h"
#include "BDS.h"
#include "OS.h"
#include "Context.h"
@@ -29,6 +32,7 @@
#include "simple3D.h"
#include "Levy3D.h"
#include "directions3D.h"
+#include "discreteFace.h"
#if defined(HAVE_ANN)
#include "ANN/ANN.h"
@@ -350,6 +354,8 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
// SHOULD be classified on the model face and get the right set of parametric
// coordinates.
+ const unsigned int initialSize = numberedV.size();
+
for(int i = numberedV.size(); i < out.numberofpoints; i++){
MVertex *v = new MVertex(out.pointlist[i * 3 + 0],
out.pointlist[i * 3 + 1],
@@ -363,6 +369,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
// Tetgen modifies both surface & edge mesh, so we need to re-create
// everything
+
gr->model()->destroyMeshCaches();
std::list<GFace*> faces = gr->faces();
for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
@@ -376,6 +383,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
gf->deleteVertexArrays();
}
+
// TODO: re-create 1D mesh
/*for(int i = 0; i < out.numberofedges; i++){
MVertex *v[2];
@@ -389,6 +397,8 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
// re-create the triangular meshes FIXME: this can lead to hanging nodes for
// non manifold geometries (single surface connected to volume)
for(int i = 0; i < out.numberoftrifaces; i++){
+ // printf("%d %d %d\n",i,out.numberoftrifaces,needParam);
+
if (out.trifacemarkerlist[i] > 0) {
MVertex *v[3];
v[0] = numberedV[out.trifacelist[i * 3 + 0] - 1];
@@ -438,7 +448,9 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
}
for(int j = 0; j < 3; j++){
- if(v[j]->onWhat()->dim() == 3){
+ if (out.trifacelist[i * 3 + j] - 1 >= initialSize){
+ printf("aaaaaaaaaaaaaaargh\n");
+ // if(v[j]->onWhat()->dim() == 3){
v[j]->onWhat()->mesh_vertices.erase
(std::find(v[j]->onWhat()->mesh_vertices.begin(),
v[j]->onWhat()->mesh_vertices.end(),
@@ -461,17 +473,20 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
else{
v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
}
-
+
gf->mesh_vertices.push_back(v1b);
numberedV[out.trifacelist[i * 3 + j] - 1] = v1b;
delete v[j];
v[j] = v1b;
}
}
+ // printf("new triangle %d %d %d\n",v[0]->onWhat()->dim(), v[1]->onWhat()->dim(), v[2]->onWhat()->dim());
MTriangle *t = new MTriangle(v[0], v[1], v[2]);
gf->triangles.push_back(t);
}// do not do anything with prismatic faces
}
+
+
for(int i = 0; i < out.numberoftetrahedra; i++){
MVertex *v1 = numberedV[out.tetrahedronlist[i * 4 + 0] - 1];
MVertex *v2 = numberedV[out.tetrahedronlist[i * 4 + 1] - 1];
@@ -483,9 +498,543 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
}
#endif
-static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
+static void addOrRemove(const MFace &f,
+ MElement *e,
+ std::map<MFace,MElement*,Less_Face> & bfaces)
{
- // buildAdditionalPoints3D (gr);
+ std::map<MFace,MElement*,Less_Face>::iterator it = bfaces.find(f);
+ if (it == bfaces.end())bfaces.insert(std::make_pair(f,e));
+ else bfaces.erase(it);
+}
+
+
+/*
+ here, we have a list of elements that actually do not form a mesh
+ --> a set boundary layer elements (prism, hexes, pyramids (and soon tets)
+ --> a set of tetrahedra that respect the external boundary of the BL mesh,
+ yet possiblty containing elements that are on the other side of the boundary
+ and therefore overlapping those elements. We have to extract the good ones !
+
+
+
+*/
+
+bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
+ std::vector<MTetrahedron*> &tets,
+ std::vector<MHexahedron*> &hexes,
+ std::vector<MPrism*> &prisms,
+ std::vector<MPyramid*> &pyramids)
+{
+ std::set<MElement*> all;
+ all.insert(hexes.begin(),hexes.end());
+ all.insert(prisms.begin(),prisms.end());
+ all.insert(pyramids.begin(),pyramids.end());
+ // start with one BL element !
+ MElement *FIRST = all.size() ? *(all.begin()) : 0;
+ if (!FIRST) return true;
+ all.insert(tets.begin(),tets.end());
+
+ // printf("coucou1 %d eleemnts\n",all.size());
+ fs_cont search;
+ buildFaceSearchStructure(gr->model(), search);
+
+ // create the graph face 2 elements for the region
+ std::map<MFace,std::pair<MElement*,MElement*> ,Less_Face> myGraph;
+
+ for (std::set<MElement*>::iterator it = all.begin(); it != all.end(); ++it){
+ MElement *t = *it;
+ const int nbf = t->getNumFaces();
+ for (int j=0;j<nbf;j++){
+ MFace f = t->getFace(j);
+ std::map<MFace,std::pair<MElement*,MElement*>,Less_Face>::iterator itf = myGraph.find(f);
+ if (itf == myGraph.end())myGraph.insert (std::make_pair (f, std::make_pair (t,(MElement*)0)));
+ else {
+ // what to do ??????
+ // two tets and one prism --> the prism should be
+ // geometrically on the other side of the
+ if (itf->second.second) {
+ MElement *prism=0, *t1=0, *t2=0;
+ if (itf->second.second->getType () == TYPE_PRI || itf->second.second->getType () == TYPE_PYR) {
+ prism = itf->second.second;
+ t1 = itf->second.first;
+ t2 = t;
+ }
+ else if (itf->second.first->getType () == TYPE_PRI || itf->second.first->getType () == TYPE_PYR) {
+ prism = itf->second.first;
+ t1 = itf->second.second;
+ t2 = t;
+ }
+ else if (t->getType () == TYPE_PRI || t->getType () == TYPE_PYR) {
+ prism = t;
+ t1 = itf->second.second;
+ t2 = itf->second.first;
+ }
+ else {
+ printf("types %d %d %d\n",t->getType () ,itf->second.first->getType (),itf->second.second->getType () );
+ }
+ if (!t1 || !t2 || !prism){
+ gr->model()->writeMSH("ooops.msh");
+ Msg::Error ("Wrong BL configuration");
+ return false;
+ }
+ SVector3 n = f.normal();
+ SPoint3 c = f.barycenter();
+ MVertex *v_prism,*v_t1,*v_t2;
+ for (int i=0;i<4;i++){
+ if (t1->getVertex(i) != f.getVertex(0) &&
+ t1->getVertex(i) != f.getVertex(1) &&
+ t1->getVertex(i) != f.getVertex(2))v_t1 = t1->getVertex(i);
+ if (t2->getVertex(i) != f.getVertex(0) &&
+ t2->getVertex(i) != f.getVertex(1) &&
+ t2->getVertex(i) != f.getVertex(2))v_t2 = t2->getVertex(i);
+ }
+ for (int i=0;i<prism->getNumVertices();i++){
+ if (prism->getVertex(i) != f.getVertex(0) &&
+ prism->getVertex(i) != f.getVertex(1) &&
+ prism->getVertex(i) != f.getVertex(2))v_prism = prism->getVertex(i);
+ }
+ SVector3 vf1 (v_t1->x() - c.x(),v_t1->y() - c.y(),v_t1->z() - c.z());
+ SVector3 vf2 (v_t2->x() - c.x(),v_t2->y() - c.y(),v_t2->z() - c.z());
+ SVector3 vfp (v_prism->x() - c.x(),v_prism->y() - c.y(),v_prism->z() - c.z());
+ // printf("%12.5E %12.5E%12.5E\n",dot(vf1,n),dot(vf2,n),dot(vfp,n));
+ if (dot (vf1,n) * dot (vfp,n) > 0){itf->second.first = prism;itf->second.second=t2; /*delete t1;*/}
+ else if (dot (vf2,n) * dot (vfp,n) > 0){itf->second.first = prism;itf->second.second=t1; /*delete t2;*/}
+ // else if (dot (vf2,vfp) > 0){itf->second.first = prism;itf->second.second=t2;}
+ else Msg::Fatal("Impossible Geom Config");
+ }
+ else itf->second.second = t;
+ }
+ }
+ }
+
+ std::stack<MElement*> myStack;
+ std::set<MElement*> connected;
+ std::set<GFace*> faces_bound;
+ myStack.push(FIRST);
+ while (!myStack.empty()){
+ FIRST = myStack.top();
+ myStack.pop();
+ connected.insert(FIRST);
+ for (int i=0;i<FIRST->getNumFaces();i++){
+ MFace f = FIRST->getFace(i);
+ GFace* gfound = findInFaceSearchStructure (f.getVertex(0),
+ f.getVertex(1),
+ f.getVertex(2),
+ search);
+ if (!gfound){
+ std::map<MFace,std::pair<MElement*,MElement*>,Less_Face>::iterator
+ itf = myGraph.find(f);
+ MElement *t_neigh = itf->second.first == FIRST ?
+ itf->second.second : itf->second.first;
+ if (connected.find(t_neigh) == connected.end())myStack.push(t_neigh);
+ }
+ else {
+ // if (faces_bound.find(gfound) == faces_bound.end())printf("face %d\n",gfound->tag());
+ faces_bound.insert(gfound);
+ }
+ }
+ }
+ // printf ("found a set of %d elements that are connected with %d bounding faces\n",connected.size(),faces_bound.size());
+ GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
+ // printf("REGION %d %d\n",myGRegion->tag(),gr->tag());
+ if (myGRegion == gr){
+ // printf("one region %d is found : %d elements\n",myGRegion->tag(),connected.size());
+ myGRegion->tetrahedra.clear();
+ for (std::set<MElement*>::iterator it = connected.begin(); it != connected.end(); ++it){
+ MElement *t = *it;
+ std::set<MVertex*> _mesh_vertices;
+ for (int i=0;i<t->getNumVertices();i++){
+ MVertex *_v = t->getVertex(i);
+ if (_v->onWhat() == gr){
+ _mesh_vertices.insert(_v);
+ }
+ }
+ // myGRegion->mesh_vertices.insert(myGRegion->mesh_vertices.end(),_mesh_vertices.begin(),_mesh_vertices.end());
+
+ if (t->getType() == TYPE_TET)
+ myGRegion->tetrahedra.push_back((MTetrahedron*)t);
+ else if (t->getType() == TYPE_HEX)
+ myGRegion->hexahedra.push_back((MHexahedron*)t);
+ else if (t->getType() == TYPE_PRI)
+ myGRegion->prisms.push_back((MPrism*)t);
+ else if (t->getType() == TYPE_PYR)
+ myGRegion->pyramids.push_back((MPyramid*)t);
+ }
+ }
+ else {
+ return false;
+ }
+ return true;
+}
+
+
+static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
+{
+ if (getBLField(gr->model())) insertVerticesInRegion(gr,-1);
+ BoundaryLayerColumns* _columns = buildAdditionalPoints3D (gr);
+ if (!_columns)return false;
+ std::map<MFace,MElement*,Less_Face> bfaces;
+
+ std::vector<MPrism*> blPrisms;
+ std::vector<MHexahedron*> blHexes;
+ std::vector<MPyramid*> blPyrs;
+ std::list<GFace*> faces = gr->faces();
+
+
+
+ std::list<GFace*> embedded_faces = gr->embeddedFaces();
+ faces.insert(faces.begin(), embedded_faces.begin(),embedded_faces.end());
+ FILE *ff2 = fopen ("tato3D.pos","w");
+ fprintf(ff2,"View \" \"{\n");
+ std::set<MVertex*> verts;
+
+ std::list<GFace*>::iterator itf = faces.begin();
+ while(itf != faces.end()){
+ for(unsigned int i = 0; i< (*itf)->triangles.size(); i++){
+ MVertex *v1 = (*itf)->triangles[i]->getVertex(0);
+ MVertex *v2 = (*itf)->triangles[i]->getVertex(1);
+ MVertex *v3 = (*itf)->triangles[i]->getVertex(2);
+ MFace dv(v1,v2,v3);
+ addOrRemove (dv,0,bfaces);
+ for (unsigned int SIDE = 0 ; SIDE < _columns->_normals3D.count(dv); SIDE ++){
+ faceColumn fc = _columns->getColumns(*itf,v1, v2, v3, SIDE);
+ const BoundaryLayerData & c1 = fc._c1;
+ const BoundaryLayerData & c2 = fc._c2;
+ const BoundaryLayerData & c3 = fc._c3;
+ int N = std::min(c1._column.size(),std::min(c2._column.size(),c3._column.size()));
+ // printf("%d Layers\n",N);
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v13,*v21,*v22,*v23;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ v23 = c3._column[l];
+ if (l == 0){
+ v11 = v1;
+ v12 = v2;
+ v13 = v3;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ v13 = c3._column[l-1];
+ }
+ // printf("coucoucouc %p %p %p %p %p %p\n",v11,v12,v13,v21,v22,v23);
+ blPrisms.push_back(new MPrism(v11,v12,v13,v21,v22,v23));
+ fprintf(ff2,"SI (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1,1,1};\n",
+ v11->x(),v11->y(),v11->z(),
+ v12->x(),v12->y(),v12->z(),
+ v13->x(),v13->y(),v13->z(),
+ v21->x(),v21->y(),v21->z(),
+ v22->x(),v22->y(),v22->z(),
+ v23->x(),v23->y(),v23->z());
+ // printf("done %d\n",l);
+ }
+ }
+ }
+ ++itf;
+ }
+
+ std::set<MEdge,Less_Edge> edges;
+ {
+ std::list<GEdge*> gedges = gr->edges();
+ for (std::list<GEdge*>::iterator it = gedges.begin(); it != gedges.end() ; ++it){
+ for (unsigned int i=0;i<(*it)->lines.size();++i){
+ edges.insert(MEdge((*it)->lines[i]->getVertex(0),(*it)->lines[i]->getVertex(1)));
+ }
+ }
+ }
+
+ // now treat the Wedges
+ // we have to know the two target GFaces that are concerned with a GEdge
+ std::set<MEdge>::iterator ite = edges.begin();
+ while(ite != edges.end()){
+ MEdge e = *ite;
+ MVertex *v1 = e.getVertex(0);
+ MVertex *v2 = e.getVertex(1);
+ bool onWedge1 = _columns->isOnWedge(v1);
+ bool onWedge2 = _columns->isOnWedge(v2);
+ bool onCorner1 = _columns->isCorner(v1);
+ bool onCorner2 = _columns->isCorner(v2);
+ if ((onWedge1 || onCorner1) && (onWedge2 || onCorner2)){
+ int N1=0,N2=0;
+ std::vector<GFace *>gfs1,gfs2;
+ BoundaryLayerFanWedge3d w1,w2;
+ BoundaryLayerFanCorner3d c1,c2;
+ if (onWedge1) {
+ N1 = _columns->getNbColumns(v1);
+ w1 = _columns->getWedge(v1);
+ gfs1 = w1._gf1;
+ gfs2 = w1._gf2;
+ }
+ else {
+ c1 = _columns->getCorner(v1);
+ N1 = c1._fanSize;
+ }
+ if (onWedge2) {
+ N2 = _columns->getNbColumns(v2);
+ w2 = _columns->getWedge(v2);
+ gfs1 = w2._gf1;
+ gfs2 = w2._gf2;
+ }
+ else {
+ c2 = _columns->getCorner(v2);
+ N2 = c2._fanSize;
+ }
+
+ // have to go from gf1 to gf2 to be aware of the sense!!!
+ if (N1 == N2){
+ for (unsigned int i=0;i<N1-1;i++){
+
+ unsigned int i11=i, i12=i+1, i21=i,i22=i+1;
+
+ if (onWedge1){
+ // if (w1.isLeft(gfs1)){i11=i;i12=i+1;}
+ // else if (w1.isRight(gfs1)){i11=N1-1-i;i12=N1-2-i;}
+ // printf("%d %d %d\n",w1.isLeft(gfs1),gfs1[0] == w1._gf1[0],gfs1.size());
+ // printf("%d %d %d\n",w1.isRight(gfs1),gfs1[0] == w1._gf2[0],gfs1.size());
+ if (gfs1[0] == w1._gf1[0]){i11=i;i12=i+1;}
+ else if (gfs1[0] == w1._gf2[0]){i11=N1-1-i;i12=N1-2-i;}
+ }
+ else {
+ /*
+
+
+
+ | 0 --> column 0
+ |
+ | fan with N1 columns
+ gf.size |
+ ---------- + ---------- 1 --> column N1-1
+ |
+ |
+ |
+ | 2 --> column 2*(N1-1)
+
+ 0 1 2 3 --> 0 -> 4-1
+ 3 4 5 6 --> 4-1 -> 2 (4-1)
+ 6 7 8 0
+
+ 1 --> 0 ==>
+
+ */
+ int K1=-1,K2=-1;
+ for (unsigned int s=0;s < c1._gf.size();s++){
+ if (w2.isLeft(c1._gf[s]))K1 = s;
+ if (w2.isRight(c1._gf[s]))K2 = s;
+ }
+ if (K1==-1)Msg::Error("K1 = -1");
+ if (K2==-1)Msg::Error("K2 = -1");
+ if (K1+1==K2){i11=K1*(N1-1)+i;i12=i11+1;}
+ else if (K2+1 == K1){i11=K1*(N1-1)-i;i12=i11-1;}
+ else if (K2 == 0 && K1 == c1._gf.size() - 1){i11=K1*(N1-1)+i;i12=i11+1;}
+ else if (K1 == 0 && K2 == c1._gf.size() - 1){i11=(K2+1)*(N1-1)-i;i12=i11-1;}
+ else Msg::Error("KROUPOUK 1 %d %d !",K1,K2);
+ if (i12 == (c1._gf.size()) * (N1-1))i12=0;
+ if (i11 == (c1._gf.size()) * (N1-1))i11=0;
+ }
+ if (onWedge2){
+ if (w2.isLeft(gfs1)){i21=i;i22=i+1;}
+ else if (w2.isRight(gfs1)){i21=N1-1-i;i22=N1-2-i;}
+ }
+ else {
+ int K1,K2;
+ for (unsigned int s=0;s<c2._gf.size();s++){
+ if (w1.isLeft(c2._gf[s]))K1 = s;
+ if (w1.isRight(c2._gf[s]))K2 = s;
+ }
+ if (K1+1 == K2){i21=K1*(N1-1)+i;i22=i21+1;}
+ else if (K2+1 == K1){i21=K1*(N1-1)-i;i22=i21-1;}
+ else if (K2 == 0 && K1 == c2._gf.size()-1 ){i21=K1*(N1-1)+i;i22=i21+1;}
+ else if (K1 == 0 && K2 == c2._gf.size()-1){i21=(K2+1)*(N1-1)-i;i22=i21-1;}
+ else Msg::Error("KROUPOUK 2 %d %d !",K1,K2);
+ if (i22 == (c2._gf.size()) * (N1-1))i22=0;
+ if (i21 == (c2._gf.size()) * (N1-1))i21=0;
+ }
+
+ BoundaryLayerData c11 = _columns->getColumn(v1,i11);
+ BoundaryLayerData c12 = _columns->getColumn(v1,i12);
+ BoundaryLayerData c21 = _columns->getColumn(v2,i21);
+ BoundaryLayerData c22 = _columns->getColumn(v2,i22);
+ int N = std::min(c11._column.size(),std::min(c12._column.size(), std::min(c21._column.size(),c22._column.size())));
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v13,*v14;
+ MVertex *v21,*v22,*v23,*v24;
+ v21 = c11._column[l];
+ v22 = c12._column[l];
+ v23 = c22._column[l];
+ v24 = c21._column[l];
+ if (l == 0){
+ v11 = v12 = v1;
+ v13 = v14 = v2;
+ }
+ else {
+ v11 = c11._column[l-1];
+ v12 = c12._column[l-1];
+ v13 = c22._column[l-1];
+ v14 = c21._column[l-1];
+ }
+
+ if (l == 0){
+ blPrisms.push_back(new MPrism(v12,v21,v22,v13,v24,v23));
+ fprintf(ff2,"SI (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){3,3,3,3,3,3};\n",
+ v12->x(),v12->y(),v12->z(),
+ v21->x(),v21->y(),v21->z(),
+ v22->x(),v22->y(),v22->z(),
+ v13->x(),v13->y(),v13->z(),
+ v24->x(),v24->y(),v24->z(),
+ v23->x(),v23->y(),v23->z());
+ }
+ else {
+ blHexes.push_back(new MHexahedron(v11,v12,v13,v14,v21,v22,v23,v24));
+ fprintf(ff2,"SH (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2,2,2,2,2,2};\n",
+ v11->x(),v11->y(),v11->z(),
+ v12->x(),v12->y(),v12->z(),
+ v13->x(),v13->y(),v13->z(),
+ v14->x(),v14->y(),v14->z(),
+ v21->x(),v21->y(),v21->z(),
+ v22->x(),v22->y(),v22->z(),
+ v23->x(),v23->y(),v23->z(),
+ v24->x(),v24->y(),v24->z());
+ }
+ }
+ }
+ }
+ else{
+ Msg::Error("cannot create 3D BL FAN %d %d -- %d %d %d %d",N1,N2,onWedge1,onWedge1,onCorner1,onCorner2);
+ }
+ }
+ ++ite;
+ }
+
+ fprintf(ff2,"};\n");
+ fclose(ff2);
+
+
+ for (unsigned int i = 0; i < blPrisms.size();i++){
+ for (unsigned int j=0;j<5;j++)
+ addOrRemove(blPrisms[i]->getFace(j),blPrisms[i],bfaces);
+ for (int j = 0; j < 6; j++)
+ if(blPrisms[i]->getVertex(j)->onWhat() == gr)verts.insert(blPrisms[i]->getVertex(j));
+ }
+ for (unsigned int i = 0; i < blHexes.size();i++){
+ for (unsigned int j=0;j<6;j++)
+ addOrRemove(blHexes[i]->getFace(j),blHexes[i],bfaces);
+ for (int j = 0; j < 8; j++)
+ if(blHexes[i]->getVertex(j)->onWhat() == gr)verts.insert(blHexes[i]->getVertex(j));
+ }
+
+ discreteFace *nf = new discreteFace(gr->model(), 444444);
+ gr->model()->add (nf);
+ std::list<GFace*> hop;
+ std::map<MFace,MElement*,Less_Face>::iterator it = bfaces.begin();
+
+ FILE *ff = fopen ("toto3D.pos","w");
+ fprintf(ff,"View \" \"{\n");
+ for (; it != bfaces.end(); ++it){
+ if (it->first.getNumVertices() == 3){
+ nf->triangles.push_back(new MTriangle (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2)));
+ fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n",
+ it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
+ it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
+ it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z());
+ }
+ else {
+
+ // we have a quad face --> create a pyramid
+
+ MElement *e = it->second;
+
+ // compute the center of the face;
+ SPoint3 center (0.25*(it->first.getVertex(0)->x()+it->first.getVertex(1)->x()+it->first.getVertex(2)->x()+it->first.getVertex(3)->x()),
+ 0.25*(it->first.getVertex(0)->y()+it->first.getVertex(1)->y()+it->first.getVertex(2)->y()+it->first.getVertex(3)->y()),
+ 0.25*(it->first.getVertex(0)->z()+it->first.getVertex(1)->z()+it->first.getVertex(2)->z()+it->first.getVertex(3)->z()));
+ // compute the center of the opposite face;
+ SPoint3 opposite (0,0,0);
+ int counter=0;
+ for (int i=0;i<e->getNumVertices();i++){
+ MVertex *vv = e->getVertex(i);
+ bool found = false;
+ for (int j=0;j<4;j++){
+ MVertex *ww = it->first.getVertex(j);
+ if (ww == vv)found = true;
+ }
+ if (!found){
+ counter ++;
+ opposite += SPoint3(vv->x(),vv->y(),vv->z());
+ }
+ }
+ // printf("counter = %d\n",counter);
+ opposite /= (double)counter;
+
+ SVector3 dir = center - opposite;
+ MTriangle temp (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2));
+ double D = temp.minEdge();
+ dir.normalize();
+ const double eps = D*1.e-2;
+ MVertex *newv = new MVertex (center.x() + eps * dir.x(),center.y() + eps * dir.y(), center.z() + eps * dir.z(), gr);
+
+ MPyramid *pyr = new MPyramid (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2),it->first.getVertex(3),newv);
+ verts.insert(newv);
+ blPyrs.push_back(pyr);
+
+ nf->triangles.push_back(new MTriangle (it->first.getVertex(0),it->first.getVertex(1),newv));
+ nf->triangles.push_back(new MTriangle (it->first.getVertex(1),it->first.getVertex(2),newv));
+ nf->triangles.push_back(new MTriangle (it->first.getVertex(2),it->first.getVertex(3),newv));
+ nf->triangles.push_back(new MTriangle (it->first.getVertex(3),it->first.getVertex(0),newv));
+
+ fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
+ it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
+ it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
+ newv->x(),newv->y(),newv->z());
+
+ fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
+ it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
+ it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(),
+ newv->x(),newv->y(),newv->z());
+
+ fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
+ it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(),
+ it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(),
+ newv->x(),newv->y(),newv->z());
+ fprintf(ff,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){2,2,2};\n",
+ it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z(),
+ it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
+ newv->x(),newv->y(),newv->z());
+
+ fprintf(ff,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){3,3,3,3};\n",
+ it->first.getVertex(0)->x(),it->first.getVertex(0)->y(),it->first.getVertex(0)->z(),
+ it->first.getVertex(1)->x(),it->first.getVertex(1)->y(),it->first.getVertex(1)->z(),
+ it->first.getVertex(2)->x(),it->first.getVertex(2)->y(),it->first.getVertex(2)->z(),
+ it->first.getVertex(3)->x(),it->first.getVertex(3)->y(),it->first.getVertex(3)->z());
+ }
+ }
+ fprintf(ff,"};\n");
+ fclose(ff);
+
+ printf("discrete face with %d %d elems\n",nf->triangles.size(),nf->quadrangles.size());
+ hop.push_back(nf);
+
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++) delete gr->tetrahedra[i];
+ gr->tetrahedra.clear();
+
+ gr->set(hop);
+ CreateAnEmptyVolumeMesh(gr);
+ printf("%d tets\n",gr->tetrahedra.size());
+ deMeshGFace _kill;
+ _kill (nf);
+ delete nf;
+ gr->model()->remove(nf);
+
+ gr->set(faces);
+ gr->mesh_vertices.insert(gr->mesh_vertices.begin(),verts.begin(),verts.end());
+
+ gr->model()->writeMSH("BL_start.msh");
+
+ AssociateElementsToModelRegionWithBoundaryLayers (gr, gr->tetrahedra , blHexes, blPrisms, blPyrs);
+
+ gr->model()->writeMSH("BL_start2.msh");
+
+ return true;
}
void _relocateVertex(MVertex *ver,
@@ -515,6 +1064,30 @@ void _relocateVertex(MVertex *ver,
}
}
+bool CreateAnEmptyVolumeMesh(GRegion *gr){
+ printf("creating an empty volume mesh\n");
+ splitQuadRecovery sqr;
+ tetgenio in, out;
+ std::vector<MVertex*> numberedV;
+ char opts[128];
+ buildTetgenStructure(gr, in, numberedV, sqr);
+ printf("tetgen structure created\n");
+ sprintf(opts, "-Ype%c",
+ (Msg::GetVerbosity() < 3) ? 'Q':
+ (Msg::GetVerbosity() > 6) ? 'V': '\0');
+ try{
+ tetrahedralize(opts, &in, &out);
+ }
+ catch (int error){
+ Msg::Error("Self intersecting surface mesh");
+ return false;
+ }
+ printf("creating an empty volume mesh done\n");
+ TransferTetgenMesh(gr, in, out, numberedV);
+ return true;
+}
+
+
void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
{
if(regions.empty()) return;
@@ -612,6 +1185,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
TransferTetgenMesh(gr, in, out, numberedV);
}
+
// sort triangles in all model faces in order to be able to search in vectors
std::list<GFace*>::iterator itf = allFaces.begin();
while(itf != allFaces.end()){
@@ -623,21 +1197,21 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
- modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr);
+ bool _BL = modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr);
// now do insertion of points
- if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
- bowyerWatsonFrontalLayers(gr, false);
- else if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_HEX)
- bowyerWatsonFrontalLayers(gr, true);
- else if(CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D){
- refineMeshMMG(gr);
- }
- else
- if(!Filler::get_nbr_new_vertices() && !LpSmoother::get_nbr_interior_vertices()){
- insertVerticesInRegion(gr);
- }
-
+ if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
+ bowyerWatsonFrontalLayers(gr, false);
+ else if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_HEX)
+ bowyerWatsonFrontalLayers(gr, true);
+ else if(CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D){
+ refineMeshMMG(gr);
+ }
+ else
+ if(!Filler::get_nbr_new_vertices() && !LpSmoother::get_nbr_interior_vertices()){
+ insertVerticesInRegion(gr,2000000000,!_BL);
+ }
+
//emi test frame field
// int NumSmooth = 10;//CTX::instance()->mesh.smoothCrossField
// std::cout << "NumSmooth = " << NumSmooth << std::endl;
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index 07b4079799..b0574663dd 100644
--- a/Mesh/meshGRegion.h
+++ b/Mesh/meshGRegion.h
@@ -6,6 +6,7 @@
#ifndef _MESH_GREGION_H_
#define _MESH_GREGION_H_
+#include <list>
#include <vector>
#include <map>
@@ -49,6 +50,7 @@ class deMeshGRegion {
};
void MeshDelaunayVolume(std::vector<GRegion*> &delaunay);
+bool CreateAnEmptyVolumeMesh(GRegion *gr);
int MeshTransfiniteVolume(GRegion *gr);
int SubdivideExtrudedMesh(GModel *m);
void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces);
diff --git a/Mesh/meshGRegionBoundaryLayers.cpp b/Mesh/meshGRegionBoundaryLayers.cpp
deleted file mode 100644
index db6f25e753..0000000000
--- a/Mesh/meshGRegionBoundaryLayers.cpp
+++ /dev/null
@@ -1,212 +0,0 @@
-// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@geuz.org>.
-
-#include "boundaryLayersData.h"
-#include "GModel.h"
-#include "GRegion.h"
-#include "MTriangle.h"
-#include "MTetrahedron.h"
-#include "MVertex.h"
-#include "Field.h"
-#include "OS.h"
-
-/*
- ^ ni
- |
- |
- +-----------------+
- bi /
- bj /
- /\
- / \ nj
- / Z
- +
-*/
-
-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);
-}
-
-BoundaryLayerColumns* buildAdditionalPoints3D_CAD_BASED (GRegion *gr)
-{
- return 0;
-}
-
-BoundaryLayerColumns* buildAdditionalPoints3D (GRegion *gr)
-{
-#if !defined(HAVE_ANN)
- return 0;
-#else
- FieldManager *fields = gr->model()->getFields();
- if(fields->getBoundaryLayerField() <= 0){
- return 0;
- }
- Field *bl_field = fields->get(fields->getBoundaryLayerField());
- BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField*> (bl_field);
-
- if (!blf)return 0;
-
- double _treshold = blf->fan_angle * M_PI / 180 ;
-
- BoundaryLayerColumns * _columns = new BoundaryLayerColumns;
-
- std::list<GFace*> faces = gr->faces();
- std::list<GFace*>::iterator itf = faces.begin();
- std::set<MVertex*> _vertices;
- std::map<MFace,SVector3,Less_Face> _normals;
-
- // 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++)
- 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));
- _normals [(*itf)->triangles[i]->getFace(0)] = SVector3(0,0,0);
- }
- }
- }
- ++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);
- }
- }
- }
-
- // for all boundry points
- for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
- std::vector<MTriangle*> _connections;
- std::vector<SVector3> _dirs, _allDirs;
- for (std::multimap<MVertex*,MTriangle*>::iterator itm =
- _columns->_non_manifold_faces.lower_bound(*it);
- itm != _columns->_non_manifold_faces.upper_bound(*it); ++itm){
- _connections.push_back (itm->second);
- _allDirs.push_back (_normals[itm->second->getFace(0)]);
- }
-
- // a list of normals is associated to a given vertex
- SPoint3 p ((*it)->x(),(*it)->y(),(*it)->z());
- std::vector<std::vector<int> > groupsOfDirections;
- for (unsigned int i=0;i<_connections.size();i++){
- SPoint3 bi = _connections[i]->barycenter();
- SVector3 ni = _allDirs[i];
- bool found = false;
- for (unsigned int j=0;j<groupsOfDirections.size();j++){
- for (unsigned int k=0;k<groupsOfDirections[j].size();k++){
- int iDir = groupsOfDirections[j][k];
- SPoint3 bj = _connections[iDir]->barycenter();
- SVector3 nj = _allDirs[iDir];
- double angle = solidAngle (ni,nj,bi,bj);
- // those two directions are
- if (angle < _treshold){
- found = true;
- groupsOfDirections[j].push_back(i);
- goto gotit;
- }
- }
- }
- gotit:
- if (!found){
- std::vector<int> newDir;
- newDir.push_back(i);
- groupsOfDirections.push_back(newDir);
- }
- }
-
-
- std::vector<SVector3> shoot;
- for (unsigned int j=0;j<groupsOfDirections.size();j++){
- SVector3 n;
- for (unsigned int k=0;k<groupsOfDirections[j].size();k++){
- int iDir = groupsOfDirections[j][k];
- n += _allDirs[iDir];
- }
- n.normalize();
- shoot.push_back(n);
- }
-
- // now create the BL points
- for (unsigned int DIR=0;DIR<shoot.size();DIR++){
- SVector3 n = shoot[DIR];
- MVertex *current = *it;
- SVector3 basis (current->x(),current->y(),current->z());
- double H = blf->hwall_n;
- std::vector<MVertex*> _column;
- std::vector<SMetric3> _metrics;
- double dist = H;
- while(dist < blf->thickness){
- SVector3 newp = basis + n * H;
- MVertex *_current = new MVertex (newp.x(),newp.y(),newp.z());
- _column.push_back(_current);
- H *= blf->ratio;
- dist += H;
- basis = newp;
- }
- _columns->addColumn(n,*it, _column, _metrics);
- }
- }
- // HERE WE SHOULD FILTER THE POINTS IN ORDER NOT TO HAVE POINTS THAT ARE TOO CLOSE
-
- // DEBUG STUFF
-
- FILE *f = Fopen ("test3D.pos","w");
- 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 _columns;
-#endif
- return 0;
-}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index d1455b602a..e044efcd5e 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -659,7 +659,7 @@ void non_recursive_classify(MTet4 *t, std::list<MTet4*> &theRegion,
while(!_stackounette.empty()){
t = _stackounette.top();
_stackounette.pop();
- if (!t) Msg::Error("a tet is not connected by a boundary face");
+ if (!t) Msg::Fatal("a tet is not connected by a boundary face");
if (!t->onWhat()) {
theRegion.push_back(t);
t->setOnWhat(bidon);
@@ -1141,7 +1141,7 @@ static void memoryCleanup(MTet4Factory &myFactory, std::set<MTet4*, compareTet4P
}
-void insertVerticesInRegion (GRegion *gr)
+void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
{
//printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
// sizeof(MTet4), sizeof(MTetrahedron), sizeof(MVertex));
@@ -1152,6 +1152,7 @@ void insertVerticesInRegion (GRegion *gr)
MTet4Factory myFactory(1600000);
std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
int NUM = 0;
+
{ // leave this in a block so the map gets deallocated directly
std::map<MVertex*, double> vSizesMap;
@@ -1183,32 +1184,40 @@ void insertVerticesInRegion (GRegion *gr)
// classify the tets on the right region
// Msg::Info("reclassifying %d tets", allTets.size());
- fs_cont search;
- buildFaceSearchStructure(gr->model(), search);
- for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end(); ++it){
- if(!(*it)->onWhat()){
- std::list<MTet4*> theRegion;
- std::set<GFace *> faces_bound;
- GRegion *bidon = (GRegion*)123;
- double _t1 = Cpu();
- Msg::Debug("start with a non classified tet");
- non_recursive_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
- double _t2 = Cpu();
- Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
- theRegion.size(), faces_bound.size(), _t2 - _t1);
- GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
- if(myGRegion){ // a geometrical region associated to the list of faces has been found
- Msg::Info("Found region %d", myGRegion->tag());
- for(std::list<MTet4*>::iterator it2 = theRegion.begin();
- it2 != theRegion.end(); ++it2) (*it2)->setOnWhat(myGRegion);
+ if (_classify) {
+ fs_cont search;
+ buildFaceSearchStructure(gr->model(), search);
+ for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if(!(*it)->onWhat()){
+ // printf("I'm in coucou\n");
+ std::list<MTet4*> theRegion;
+ std::set<GFace *> faces_bound;
+ GRegion *bidon = (GRegion*)123;
+ double _t1 = Cpu();
+ Msg::Debug("start with a non classified tet");
+ non_recursive_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
+ double _t2 = Cpu();
+ Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
+ theRegion.size(), faces_bound.size(), _t2 - _t1);
+ GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
+ if(myGRegion){ // a geometrical region associated to the list of faces has been found
+ Msg::Info("Found region %d", myGRegion->tag());
+ for(std::list<MTet4*>::iterator it2 = theRegion.begin();
+ it2 != theRegion.end(); ++it2) (*it2)->setOnWhat(myGRegion);
+ }
+ else // the tets are in the void
+ for(std::list<MTet4*>::iterator it2 = theRegion.begin();
+ it2 != theRegion.end(); ++it2)(*it2)->setDeleted(true);
}
- else // the tets are in the void
- for(std::list<MTet4*>::iterator it2 = theRegion.begin();
- it2 != theRegion.end(); ++it2)(*it2)->setDeleted(true);
}
+ search.clear();
+ }
+ else {
+ // FIXME ... too simple
+ for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end(); ++it)
+ (*it)->setOnWhat(gr);
}
- search.clear();
for(MTet4Factory::iterator it = allTets.begin(); it!=allTets.end(); ++it){
(*it)->setNeigh(0, 0);
@@ -1221,16 +1230,18 @@ void insertVerticesInRegion (GRegion *gr)
createAllEmbeddedFaces (gr, allEmbeddedFaces);
connectTets(allTets.begin(), allTets.end(),&allEmbeddedFaces);
Msg::Debug("All %d tets were connected", allTets.size());
-
+
// here the classification should be done
-
- int ITER = 0;
+
+ int ITER = 0, REALCOUNT = 0;
int NB_CORRECTION_OF_CAVITY = 0;
int COUNT_MISS_1 = 0;
int COUNT_MISS_2 = 0;
double t1 = Cpu();
while(1){
+ // break;
+ if (ITER > maxVert)break;
if(allTets.empty()){
Msg::Error("No tetrahedra in region %d %d", gr->tag(), allTets.size());
break;
@@ -1245,7 +1256,7 @@ void insertVerticesInRegion (GRegion *gr)
else{
if(ITER++ % 5000 == 0)
Msg::Info("%d points created - Worst tet radius is %g (PTS removed %d %d)",
- vSizes.size(), worst->getRadius(), COUNT_MISS_1,COUNT_MISS_2);
+ REALCOUNT, worst->getRadius(), COUNT_MISS_1,COUNT_MISS_2);
if(worst->getRadius() < 1) break;
double center[3];
double uvw[3];
@@ -1320,15 +1331,22 @@ void insertVerticesInRegion (GRegion *gr)
(*itc)->setDeleted(false);
delete v;
}
- else
+ else{
+ REALCOUNT++;
v->onWhat()->mesh_vertices.push_back(v);
+ }
}
else{
- // printf("coucou 2 %d\n",ITER);
- // printf("point outside %12.5E %12.5E %12.5E %12.5E %12.5E\n",VV,uvw[0], uvw[1], uvw[2],1-uvw[0]-uvw[1]-uvw[2]);
+ // printf("coucou 2 % cavity size %d\n",ITER,cavity.size());
+ // for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
+ // MTetrahedron *toto = (*itc)->tet();
+ // toto->xyz2uvw(center,uvw);
+ // printf("point outside %12.5E %12.5E %12.5E %12.5E\n",uvw[0], uvw[1], uvw[2],1-uvw[0]-uvw[1]-uvw[2]);
+ // }
myFactory.changeTetRadius(allTets.begin(), 0.0);
COUNT_MISS_2++;
for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc) (*itc)->setDeleted(false);
+ // if (cavity.size() > 10)printTets ("cavity.pos", cavity, true);
}
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index f74c034334..196c8a9159 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -173,8 +173,10 @@ class MTet4
void connectTets(std::list<MTet4*> &);
void connectTets(std::vector<MTet4*> &);
-void insertVerticesInRegion(GRegion *gr);
+void insertVerticesInRegion(GRegion *gr, int maxVert = 2000000000, bool _classify = true);
void bowyerWatsonFrontalLayers(GRegion *gr, bool hex);
+GRegion *getRegionFromBoundingFaces(GModel *model,
+ std::set<GFace *> &faces_bound);
class compareTet4Ptr
{
--
GitLab