From fa055b8a85c04216da092b3f76f8a875f9d43241 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 7 May 2013 14:38:44 +0000
Subject: [PATCH] modified boundary layers spallart almaras is kinda working
---
Geo/CMakeLists.txt | 2 +
Geo/GModel.h | 9 +
Geo/OCCFace.cpp | 26 ++-
Geo/STensor3.cpp | 2 +
.../boundaryLayersData.cpp | 177 ++++++++++++++--
Geo/boundaryLayersData.h | 107 ++++++++++
Geo/intersectCurveSurface.cpp | 45 ++--
Geo/intersectCurveSurface.h | 1 +
Mesh/CMakeLists.txt | 2 +-
Mesh/Field.cpp | 36 ++++
Mesh/Field.h | 3 +
Mesh/meshGEdge.cpp | 62 +++++-
Mesh/meshGFace.cpp | 8 +-
Mesh/meshGFaceBoundaryLayers.h | 194 ------------------
Mesh/meshGRegion.cpp | 6 +-
Mesh/meshGRegionBoundaryLayers.cpp | 2 +-
Mesh/meshGRegionDelaunayInsertion.cpp | 23 ++-
Mesh/surfaceFiller.cpp | 37 +++-
Mesh/yamakawa.cpp | 8 +-
Numeric/Numeric.cpp | 4 +
Numeric/simpleFunction.h | 2 +-
benchmarks/3d/submarine/submarine_simple.geo | 21 +-
wrappers/gmshpy/gmshGeo.i | 2 +
23 files changed, 498 insertions(+), 281 deletions(-)
rename Mesh/meshGFaceBoundaryLayers.cpp => Geo/boundaryLayersData.cpp (68%)
create mode 100644 Geo/boundaryLayersData.h
delete mode 100644 Mesh/meshGFaceBoundaryLayers.h
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index f257938696..f241e458cd 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -4,6 +4,8 @@
# bugs and problems to the public mailing list <gmsh@geuz.org>.
set(SRC
+ boundaryLayersData.cpp
+ distanceToMesh.cpp
intersectCurveSurface.cpp
GEntity.cpp STensor3.cpp
GVertex.cpp GEdge.cpp GFace.cpp GRegion.cpp
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 9b6e59444f..5f6dbf1d73 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -17,6 +17,7 @@
#include "GRegion.h"
#include "SPoint3.h"
#include "SBoundingBox3d.h"
+#include "boundaryLayersData.h"
template <class scalar> class simpleFunction;
class FM_Internals;
@@ -151,10 +152,15 @@ class GModel
std::set<int> meshPartitions;
int partitionSize[2];
+ // boundary layer columns i.e. list of vertices that form columns
+ // in boundary layers
+ BoundaryLayerColumns _columns;
+
public:
GModel(std::string name="");
virtual ~GModel();
+
// the static list of all loaded models
static std::vector<GModel*> list;
@@ -266,6 +272,9 @@ 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); }
void add(GFace *f) { faces.insert(f); }
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index ebe176dd46..a3dfb968c3 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -163,8 +163,30 @@ void OCCFace::secondDer(const SPoint2 ¶m,
GPoint OCCFace::point(double par1, double par2) const
{
double pp[2] = {par1, par2};
- gp_Pnt val = occface->Value(par1, par2);
- return GPoint(val.X(), val.Y(), val.Z(), this, pp);
+ double umin2, umax2, vmin2, vmax2;
+ ShapeAnalysis::GetFaceUVBounds(s, umin2, umax2, vmin2, vmax2);
+
+ double du = umax2 - umin2;
+ double dv = vmax2 - vmin2;
+
+ if (par1 > (umax2+.1*du) || par1 < (umin2-.1*du) ||
+ par2 > (vmax2+.1*dv) || par2 < (vmin2-.1*dv))
+ {
+ GPoint p(0,0,0, this, pp);
+ p.setNoSuccess();
+ return p;
+ }
+
+ try{
+ gp_Pnt val;
+ val = occface->Value(par1, par2);
+ return GPoint(val.X(), val.Y(), val.Z(), this, pp);
+ }
+ catch(Standard_OutOfRange){
+ GPoint p(0,0,0, this, pp);
+ p.setNoSuccess();
+ return p;
+ }
}
GPoint OCCFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index a282ad4f04..63ba5ed059 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -51,6 +51,8 @@ SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
// preserve orientation of m1 !!!
SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2)
{
+ // we should do
+ // return intersection (m1,m2);
fullMatrix<double> V(3,3);
fullVector<double> S(3);
m1.eig(V,S,true);
diff --git a/Mesh/meshGFaceBoundaryLayers.cpp b/Geo/boundaryLayersData.cpp
similarity index 68%
rename from Mesh/meshGFaceBoundaryLayers.cpp
rename to Geo/boundaryLayersData.cpp
index fc8850ef30..6b01a2c6f7 100644
--- a/Mesh/meshGFaceBoundaryLayers.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -9,7 +9,7 @@
#include "MLine.h"
#include "MTriangle.h"
#include "MEdge.h"
-#include "meshGFaceBoundaryLayers.h"
+#include "boundaryLayersData.h"
#include "Field.h"
SVector3 interiorNormal (SPoint2 p1, SPoint2 p2, SPoint2 p3)
@@ -90,11 +90,105 @@ void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3])
metric[2] = res[1][1];
}
-BoundaryLayerColumns* buildAdditionalPoints2D (GFace *gf)
+edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side)
+{
+ Equal_Edge aaa;
+ MEdge e(v1, v2);
+ std::map<MVertex*,BoundaryLayerFan>::const_iterator it1 = _fans.find(v1);
+ std::map<MVertex*,BoundaryLayerFan>::const_iterator it2 = _fans.find(v2);
+ int N1 = getNbColumns(v1) ;
+ int N2 = getNbColumns(v2) ;
+
+
+ int nbSides = _normals.count(e);
+
+ // if (nbSides != 1)printf("I'm here %d sides\n",nbSides);
+ // Standard case, only two extruded columns from the two vertices
+ if (N1 == 1 && N2 == 1) return edgeColumn(getColumn(v1,0),getColumn(v2,0));
+ // one fan on
+ if (nbSides == 1){
+ if (it1 != _fans.end() && it2 == _fans.end() ){
+ if (aaa(it1->second._e1,e))
+ return edgeColumn(getColumn (v1,0),getColumn(v2,0));
+ else
+ return edgeColumn(getColumn (v1,N1-1),getColumn(v2,0));
+ }
+ if (it2 != _fans.end() && it1 == _fans.end() ){
+ if (aaa(it2->second._e1,e))
+ return edgeColumn(getColumn (v1,0),getColumn(v2,0));
+ else
+ return edgeColumn(getColumn (v1,0),getColumn(v2,N2-1));
+ }
+ if (it2 != _fans.end() && it1 != _fans.end() ){
+ int c1, c2;
+ if (aaa(it1->second._e1,e))
+ c1 = 0;
+ else
+ c1 = N1-1;
+ if (aaa(it2->second._e1,e))
+ c2 = 0;
+ else
+ c2 = N2-1;
+ return edgeColumn(getColumn (v1,c1),getColumn(v2,c2));
+ }
+ // fan on the right
+ if (N1 == 1 || N2 == 2){
+ const BoundaryLayerData & c10 = getColumn(v1,0);
+ const BoundaryLayerData & c20 = getColumn(v2,0);
+ const BoundaryLayerData & c21 = getColumn(v2,1);
+ if (dot(c10._n,c20._n) > dot(c10._n,c21._n) ) return edgeColumn(c10,c20);
+ else return edgeColumn(c10,c21);
+ }
+ // fan on the left
+ if (N1 == 2 || N2 == 1){
+ const BoundaryLayerData & c10 = getColumn(v1,0);
+ const BoundaryLayerData & c11 = getColumn(v1,1);
+ const BoundaryLayerData & c20 = getColumn(v2,0);
+ if (dot(c10._n,c20._n) > dot(c11._n,c20._n) ) return edgeColumn(c10,c20);
+ else return edgeColumn(c11,c20);
+ }
+
+ // Msg::Error ("Impossible Boundary Layer Configuration : one side and no fans %d %d",N1,N2);
+ // FIXME WRONG
+ return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) :
+ edgeColumn (getColumn (v2,0),getColumn(v2,0));
+ }
+ else if (nbSides == 2){
+ int i1=0,i2=1,j1=0,j2=1;
+ if (it1 != _fans.end()){
+ i1 = aaa(it1->second._e1,e) ? 0 : N1-1;
+ i2 = !aaa(it1->second._e1,e) ? 0 : N1-1;
+ }
+ if (it2 != _fans.end()){
+ j1 = aaa(it2->second._e1,e) ? 0 : N2-1;
+ j2 = !aaa(it2->second._e1,e) ? 0 : N2-1;
+ }
+ const BoundaryLayerData & c10 = getColumn(v1,i1);
+ const BoundaryLayerData & c11 = getColumn(v1,i2);
+ const BoundaryLayerData & c20 = getColumn(v2,j1);
+ const BoundaryLayerData & c21 = getColumn(v2,j2);
+
+ if (side == 0){
+ if (dot(c10._n,c20._n) > dot(c10._n,c21._n) ) return edgeColumn(c10,c20);
+ else return edgeColumn(c10,c21);
+ }
+ if (side == 1){
+ if (dot(c11._n,c20._n) > dot(c11._n,c21._n) ) return edgeColumn(c11,c20);
+ else return edgeColumn(c11,c21);
+ }
+ }
+
+ Msg::Error("Not yet Done in BoundaryLayerData nbSides = %d, ",nbSides );
+ static BoundaryLayerData error;
+ static edgeColumn error2(error, error);
+ return error2;
+}
+
+
+bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns *_columns)
{
- // return 0;
#if !defined(HAVE_ANN)
- return 0;
+ return false;
#else
FieldManager *fields = gf->model()->getFields();
@@ -104,11 +198,11 @@ BoundaryLayerColumns* buildAdditionalPoints2D (GFace *gf)
Field *bl_field = fields->get(fields->getBoundaryLayerField());
BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField*> (bl_field);
- if (!blf)return 0;
+ if (!blf)return false;
double _treshold = blf->fan_angle * M_PI / 180 ;
- BoundaryLayerColumns * _columns = new BoundaryLayerColumns;
+ // BoundaryLayerColumns * _columns = new BoundaryLayerColumns;
// 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.
@@ -120,18 +214,19 @@ BoundaryLayerColumns* buildAdditionalPoints2D (GFace *gf)
std::list<GEdge*>::iterator ite = edges.begin();
std::set<MVertex*> _vertices;
while(ite != edges.end()){
- for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
- MVertex *v1 = (*ite)->lines[i]->getVertex(0);
- MVertex *v2 = (*ite)->lines[i]->getVertex(1);
- _columns->_non_manifold_edges.insert(std::make_pair(v1,v2));
- _columns->_non_manifold_edges.insert(std::make_pair(v2,v1));
- _vertices.insert(v1);
- _vertices.insert(v2);
+ if (blf->isEdgeBL ((*ite)->tag())){
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ _columns->_non_manifold_edges.insert(std::make_pair(v1,v2));
+ _columns->_non_manifold_edges.insert(std::make_pair(v2,v1));
+ _vertices.insert(v1);
+ _vertices.insert(v2);
+ }
}
++ite;
}
-
- // printf("%d boundary points\n",_vertices.size());
+ if (!_vertices.size())return false;
// assume that the initial mesh has been created i.e. that there exist
// triangles inside the domain. Triangles are used to define
@@ -215,12 +310,12 @@ BoundaryLayerColumns* buildAdditionalPoints2D (GFace *gf)
_dirs.push_back(x1);
x2.normalize();
_dirs.push_back(x2);
- printf("%g %g vs %g %g\n",N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
- printf("%g %g vs %g %g\n",N2[0].x(),N2[0].y(),N3[0].x(),N3[0].y());
- printf("%g %g vs %g %g\n",x1.x(),x1.y(),x2.x(),x2.y());
+ // printf("%g %g vs %g %g\n",N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
+ // printf("%g %g vs %g %g\n",N2[0].x(),N2[0].y(),N3[0].x(),N3[0].y());
+ // printf("%g %g vs %g %g\n",x1.x(),x1.y(),x2.x(),x2.y());
}
// STANDARD CASE, one vertex connected to two neighboring vertices
- if (_connections.size() == 2){
+ else if (_connections.size() == 2){
MEdge e1 (*it,_connections[0]);
MEdge e2 (*it,_connections[1]);
//LL = 0.5 * (e1.length() + e2.length());
@@ -288,6 +383,48 @@ BoundaryLayerColumns* buildAdditionalPoints2D (GFace *gf)
}
}
}
+ }
+ else if (_connections.size() == 1){
+ MEdge e1 (*it,_connections[0]);
+ SPoint2 p0,p1;
+ reparamMeshEdgeOnFace(*it,_connections[0], gf, p0, p1);
+ std::vector<SVector3> N1;
+ 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);
+ // two sides at that point : end point of one embedded edge
+ // the fan angle is equal to PI
+
+ if (N1.size() == 2){
+ int fanSize = M_PI / _treshold;
+ 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());
+
+ 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());
+ double AMAX = std::max(alpha1,alpha2);
+ double AMIN = std::min(alpha1,alpha2);
+ if (alpha3 > AMAX){
+ AMIN += M_PI;
+ AMAX += M_PI;
+ }
+ if ( AMAX - AMIN >= M_PI){
+ double temp = AMAX;
+ AMAX = AMIN + 2*M_PI;
+ AMIN = temp;
+ }
+ _columns->addFan (*it,e1,e1,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);
+ printf("%d %g %g %g\n",i,x.x(),x.y(),alpha);
+ x.normalize();
+ _dirs.push_back(x);
+ }
+ }
}
// if (_dirs.size() > 1)printf("%d directions\n",_dirs.size());
@@ -389,7 +526,7 @@ BoundaryLayerColumns* buildAdditionalPoints2D (GFace *gf)
// END OF DEBUG STUFF
- return _columns;
+ return 1;
#endif
}
diff --git a/Geo/boundaryLayersData.h b/Geo/boundaryLayersData.h
new file mode 100644
index 0000000000..be5f0d13b4
--- /dev/null
+++ b/Geo/boundaryLayersData.h
@@ -0,0 +1,107 @@
+// 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>.
+
+#ifndef _BNDRYLRDATA_
+#define _BNDRYLRDATA_
+
+#include "SVector3.h"
+#include "STensor3.h"
+#include "MEdge.h"
+#include <map>
+#include <set>
+
+class Field;
+class GFace;
+class GRegion;
+class MTriangle;
+
+struct BoundaryLayerData
+{
+ SVector3 _n;
+ std::vector<MVertex*> _column;
+ std::vector<SMetric3> _metrics;
+ BoundaryLayerData(){}
+ BoundaryLayerData(const SVector3 & dir,
+ std::vector<MVertex*> column,
+ std::vector<SMetric3> metrics)
+ : _n(dir), _column(column), _metrics(metrics){}
+};
+
+struct BoundaryLayerFan
+{
+ MEdge _e1, _e2;
+ bool sense;
+ BoundaryLayerFan(MEdge e1, MEdge e2 , bool s = true)
+ : _e1(e1),_e2(e2) , sense (s){}
+};
+
+struct edgeColumn {
+ const BoundaryLayerData &_c1, &_c2;
+ edgeColumn(const BoundaryLayerData &c1, const BoundaryLayerData &c2)
+ : _c1(c1), _c2(c2){}
+};
+
+class BoundaryLayerColumns
+{
+ std::multimap<MVertex*, BoundaryLayerData> _data;
+ std::map<MVertex*, BoundaryLayerFan> _fans;
+public:
+ 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;
+ 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 addFan(MVertex *v, MEdge e1, MEdge e2, bool s)
+ {
+ _fans.insert(std::make_pair(v,BoundaryLayerFan(e1,e2,s)));
+ }
+ inline const BoundaryLayerData &getColumn(MVertex *v, MEdge e)
+ {
+ std::map<MVertex*,BoundaryLayerFan>::const_iterator it = _fans.find(v);
+ int N = getNbColumns(v) ;
+ if (N == 1) return getColumn(v, 0);
+ Equal_Edge aaa;
+ if (it != _fans.end()){
+ if (aaa(it->second._e1, e))
+ return getColumn(v, 0);
+ else
+ return getColumn(v, N-1);
+ }
+ Msg::Error("Cannot handle embedded lines in boundary layers");
+ static BoundaryLayerData error;
+ return error;
+ }
+ edgeColumn getColumns(MVertex *v1, MVertex *v2 , int side);
+ inline int getNbColumns(MVertex *v) { return _data.count(v); }
+ inline const BoundaryLayerData &getColumn(MVertex *v, int iColumn)
+ {
+ int count = 0;
+ for(std::multimap<MVertex*,BoundaryLayerData>::iterator itm = _data.lower_bound(v);
+ itm != _data.upper_bound(v); ++itm){
+ if (count++ == iColumn) return itm->second;
+ }
+ static BoundaryLayerData error;
+ return error;
+ }
+ void filterPoints();
+};
+
+bool buildAdditionalPoints2D (GFace *gf, BoundaryLayerColumns * ) ;
+bool buildAdditionalPoints3D (GRegion *gr, BoundaryLayerColumns * ) ;
+void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]);
+
+#endif
diff --git a/Geo/intersectCurveSurface.cpp b/Geo/intersectCurveSurface.cpp
index 7ac1ed2584..c125bef4ab 100644
--- a/Geo/intersectCurveSurface.cpp
+++ b/Geo/intersectCurveSurface.cpp
@@ -15,22 +15,30 @@ struct intersectCurveSurfaceData
{ }
bool apply (double newPoint[3]){
- fullVector<double> uvt(3);
- uvt(0) = newPoint[0];
- uvt(1) = newPoint[1];
- uvt(2) = newPoint[2];
- fullVector<double> res(3);
- _kaboom(uvt,res,this);
- // printf("start with %12.5E\n",res.norm());
- if (res.norm() < epsilon)return true;
-
- if(newton_fd(_kaboom, uvt, this)){
- // printf("--- CONVERGED -----------\n");
- newPoint[0] = uvt(0);
- newPoint[1] = uvt(1);
- newPoint[2] = uvt(2);
- return true;
- }
+ try {
+ fullVector<double> uvt(3);
+ uvt(0) = newPoint[0];
+ uvt(1) = newPoint[1];
+ uvt(2) = newPoint[2];
+ fullVector<double> res(3);
+ _kaboom(uvt,res,this);
+ // printf("start with %12.5E\n",res.norm());
+ if (res.norm() < epsilon)return true;
+
+
+ if(newton_fd(_kaboom, uvt, this)){
+ // printf("--- CONVERGED -----------\n");
+ newPoint[0] = uvt(0);
+ newPoint[1] = uvt(1);
+ newPoint[2] = uvt(2);
+ // printf("newton done\n");
+ return true;
+ }
+ }
+ catch (...){
+ // printf("intersect curve surface failed !\n");
+ }
+ // printf("newton failsed\n");
return false;
}
};
@@ -38,12 +46,15 @@ struct intersectCurveSurfaceData
void _kaboom(fullVector<double> &uvt,
fullVector<double> &res, void *_data){
intersectCurveSurfaceData *data = (intersectCurveSurfaceData*)_data;
+ // printf("coucuo1 %g %g\n",uvt(0),uvt(1));
SPoint3 s = data->s(uvt(0),uvt(1));
+ // printf("coucuo1\n");
SPoint3 c = data->c(uvt(2));
+ // printf("coucuo1\n");
res(0) = s.x() - c.x();
res(1) = s.y() - c.y();
res(2) = s.z() - c.z();
- // printf("%g %g %g vs %g %g %g \n",s.x(),s.y(),s.z(),c.x(),c.y(),c.z());
+ // printf("%g %g %g vs %g %g %g \n",s.x(),s.y(),s.z(),c.x(),c.y(),c.z());
}
int intersectCurveSurface (curveFunctor &c, surfaceFunctor & s, double uvt[3], double epsilon){
diff --git a/Geo/intersectCurveSurface.h b/Geo/intersectCurveSurface.h
index b3c14a7375..4c9fb5d247 100644
--- a/Geo/intersectCurveSurface.h
+++ b/Geo/intersectCurveSurface.h
@@ -35,6 +35,7 @@ public :
surfaceFunctorGFace(const GFace *_gf) : gf(_gf) {}
virtual SPoint3 operator () (double u, double v) const {
GPoint gp = gf->point(u,v);
+ if (! gp.succeeded())throw gf;
return SPoint3(gp.x(),gp.y(),gp.z());
}
};
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index b50d4d34d6..5d365bec22 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -12,7 +12,7 @@ set(SRC
meshGFaceTransfinite.cpp meshGFaceExtruded.cpp
meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp
meshGFaceLloyd.cpp meshGFaceOptimize.cpp
- meshGFaceQuadrilateralize.cpp meshGFaceBoundaryLayers.cpp
+ meshGFaceQuadrilateralize.cpp
meshGRegion.cpp meshGRegionBoundaryLayers.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 0f4ea7073a..b6ed27696d 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1859,6 +1859,42 @@ double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge
return std::min (hfar,lc);
}
+// assume that the closest point is one of the model vertices
+void BoundaryLayerField::computeFor1dMesh (double x, double y, double z,
+ SMetric3 &metr)
+{
+ double xpk,ypk,zpk;
+ double distk = 1.e22;
+ for(std::list<int>::iterator it = nodes_id.begin();
+ it != nodes_id.end(); ++it) {
+ GVertex *v = GModel::current()->getVertexByTag(*it);
+ double xp = v->x();
+ double yp = v->y();
+ double zp = v->z();
+ const double dist = sqrt ((x-xp)*(x-xp)+
+ (y-yp)*(y-yp)+
+ (z-zp)*(z-zp));
+ if (dist < distk){
+ distk = dist;
+ xpk=xp;
+ ypk=yp;
+ zpk=zp;
+ }
+ }
+
+ const double ll1 = (distk*(ratio-1) + hwall_n) / (1. + 0.5 * (ratio - 1));
+ //const double ll1 = (distk*(ratio-1) + hwall_n) / (1.);
+ double lc_n = std::min(ll1,hfar);
+
+ if (distk > thickness) lc_n = hfar;
+ lc_n = std::max(lc_n, CTX::instance()->mesh.lcMin);
+ lc_n = std::min(lc_n, CTX::instance()->mesh.lcMax);
+ SVector3 t1= SVector3(x-xpk,y-ypk,z-zpk);
+ t1.normalize();
+ metr = buildMetricTangentToCurve(t1,lc_n,lc_n);
+}
+
+
void BoundaryLayerField::operator() (AttractorField *cc, double dist,
double x, double y, double z,
SMetric3 &metr, GEntity *ge)
diff --git a/Mesh/Field.h b/Mesh/Field.h
index 4ebfe745f7..743267a7b7 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -156,6 +156,9 @@ class BoundaryLayerField : public Field {
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);
};
#endif
class FieldOptionString : public FieldOption
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 3bcd86a9d3..65d0172e13 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -9,6 +9,7 @@
#include "GEdge.h"
#include "MLine.h"
#include "BackgroundMesh.h"
+#include "boundaryLayersData.h"
#include "Numeric.h"
#include "GmshMessage.h"
#include "Context.h"
@@ -96,6 +97,16 @@ static double F_Lc(GEdge *ge, double t)
static double F_Lc_aniso(GEdge *ge, double t)
{
+#if defined(HAVE_ANN)
+ FieldManager *fields = ge->model()->getFields();
+ BoundaryLayerField *blf = 0;
+ Field *bl_field = fields->get(fields->getBoundaryLayerField());
+ blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+#else
+ bool blf = false;
+#endif
+
+
GPoint p = ge->point(t);
SMetric3 lc_here;
@@ -110,6 +121,14 @@ static double F_Lc_aniso(GEdge *ge, double t)
else
lc_here = BGM_MeshMetric(ge, t, 0, p.x(), p.y(), p.z());
+#if defined(HAVE_ANN)
+ if (blf && !blf->isEdgeBL(ge->tag())){
+ SMetric3 lc_bgm;
+ blf->computeFor1dMesh ( p.x(), p.y(), p.z() , lc_bgm );
+ lc_here = intersection_conserveM1 (lc_here, lc_bgm );
+ }
+#endif
+
SVector3 der = ge->firstDer(t);
double lSquared = dot(der, lc_here, der);
return sqrt(lSquared);
@@ -294,15 +313,17 @@ static void printFandPrimitive(int tag, std::vector<IntPoint> &Points)
}
*/
+void createBoundaryLayerPointsForGEdge (GEdge *ge, double &t_begin, double &t_end)
+{
+}
+
void meshGEdge::operator() (GEdge *ge)
{
#if defined(HAVE_ANN)
FieldManager *fields = ge->model()->getFields();
BoundaryLayerField *blf = 0;
- if(fields->getBackgroundField() > 0){
- Field *bl_field = fields->get(fields->getBoundaryLayerField());
- blf = dynamic_cast<BoundaryLayerField*> (bl_field);
- }
+ Field *bl_field = fields->get(fields->getBoundaryLayerField());
+ blf = dynamic_cast<BoundaryLayerField*> (bl_field);
#else
bool blf = false;
#endif
@@ -372,9 +393,7 @@ void meshGEdge::operator() (GEdge *ge)
N = ge->meshAttributes.nbPointsTransfinite;
}
else{
- if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG
- // || CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS
- || blf){
+ if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf){
a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
CTX::instance()->mesh.lcIntegrationPrecision);
}
@@ -391,8 +410,7 @@ void meshGEdge::operator() (GEdge *ge)
}
a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.99));
-
-}
+ }
// force odd number of points if blossom is used for recombination
if(ge->meshAttributes.method != MESH_TRANSFINITE &&
@@ -477,6 +495,32 @@ void meshGEdge::operator() (GEdge *ge)
v0->z() = beg_p.z();
}
+ 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);
+ // printf ("coucou %d %d %g %g %d\n",ge->tag(),g1->tag(),t.x(),t.y(),mesh_vertices.size());
+ }
+ if (blf->isVertexBL(g1->tag())){
+ _columns->addColumn(t*-1.0, v1,invert,_metrics);
+ // printf ("coucou %d %d\n",ge->tag(),g0->tag());
+
+ }
+ }
ge->meshStatistics.status = GEdge::DONE;
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e2d9778a58..4f463cceb5 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -43,7 +43,7 @@
#include "Context.h"
#include "multiscalePartition.h"
#include "meshGFaceLloyd.h"
-#include "meshGFaceBoundaryLayers.h"
+#include "boundaryLayersData.h"
inline double myAngle(const SVector3 &a, const SVector3 &b, const SVector3 &d)
{
@@ -622,7 +622,7 @@ void filterOverlappingElements(int dim, std::vector<MElement*> &e,
}
}
if (intersection){
- printf("intersection found\n");
+ // printf("intersection found\n");
einter.push_back(el);
}
else {
@@ -633,9 +633,9 @@ void filterOverlappingElements(int dim, std::vector<MElement*> &e,
void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
{
- BoundaryLayerColumns *_columns = buildAdditionalPoints2D (gf);
+ BoundaryLayerColumns* _columns = gf->model()->getColumns();
+ if (!buildAdditionalPoints2D (gf, _columns))return;
- if (!_columns)return;
std::set<MEdge,Less_Edge> bedges;
diff --git a/Mesh/meshGFaceBoundaryLayers.h b/Mesh/meshGFaceBoundaryLayers.h
deleted file mode 100644
index e6b4917ccd..0000000000
--- a/Mesh/meshGFaceBoundaryLayers.h
+++ /dev/null
@@ -1,194 +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>.
-
-#ifndef _MESHGFACE_BNDRYLR_
-#define _MESHGFACE_BNDRYLR_
-
-#include "SVector3.h"
-#include "STensor3.h"
-#include "MEdge.h"
-#include <map>
-#include <set>
-
-class Field;
-class GFace;
-class GRegion;
-class MTriangle;
-
-struct BoundaryLayerData
-{
- SVector3 _n;
- std::vector<MVertex*> _column;
- std::vector<SMetric3> _metrics;
- BoundaryLayerData(){}
- BoundaryLayerData(const SVector3 & dir,
- std::vector<MVertex*> column,
- std::vector<SMetric3> metrics)
- : _n(dir), _column(column), _metrics(metrics){}
-};
-
-struct BoundaryLayerFan
-{
- MEdge _e1, _e2;
- bool sense;
- BoundaryLayerFan(MEdge e1, MEdge e2 , bool s = true)
- : _e1(e1),_e2(e2) , sense (s){}
-};
-
-struct edgeColumn {
- const BoundaryLayerData &_c1, &_c2;
- edgeColumn(const BoundaryLayerData &c1, const BoundaryLayerData &c2)
- : _c1(c1), _c2(c2){}
-};
-
-class BoundaryLayerColumns
-{
- std::multimap<MVertex*, BoundaryLayerData> _data;
- std::map<MVertex*, BoundaryLayerFan> _fans;
-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;
- 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 addFan(MVertex *v, MEdge e1, MEdge e2, bool s)
- {
- _fans.insert(std::make_pair(v,BoundaryLayerFan(e1,e2,s)));
- }
- inline const BoundaryLayerData &getColumn(MVertex *v, MEdge e)
- {
- std::map<MVertex*,BoundaryLayerFan>::const_iterator it = _fans.find(v);
- int N = getNbColumns(v) ;
- if (N == 1) return getColumn(v, 0);
- Equal_Edge aaa;
- if (it != _fans.end()){
- if (aaa(it->second._e1, e))
- return getColumn(v, 0);
- else
- return getColumn(v, N-1);
- }
- Msg::Error("Cannot handle embedded lines in boundary layers");
- static BoundaryLayerData error;
- return error;
- }
- inline edgeColumn getColumns(MVertex *v1, MVertex *v2 , int side)
- {
- Equal_Edge aaa;
- MEdge e(v1, v2);
- std::map<MVertex*,BoundaryLayerFan>::const_iterator it1 = _fans.find(v1);
- std::map<MVertex*,BoundaryLayerFan>::const_iterator it2 = _fans.find(v2);
- int N1 = getNbColumns(v1) ;
- int N2 = getNbColumns(v2) ;
-
- int nbSides = _normals.count(e);
-
- // if (nbSides != 1)printf("I'm here %d sides\n",nbSides);
- // Standard case, only two extruded columns from the two vertices
- if (N1 == 1 && N2 == 1) return edgeColumn(getColumn(v1,0),getColumn(v2,0));
- // one fan on
- if (nbSides == 1){
- if (it1 != _fans.end() && it2 == _fans.end() ){
- if (aaa(it1->second._e1,e))
- return edgeColumn(getColumn (v1,0),getColumn(v2,0));
- else
- return edgeColumn(getColumn (v1,N1-1),getColumn(v2,0));
- }
- if (it2 != _fans.end() && it1 == _fans.end() ){
- if (aaa(it2->second._e1,e))
- return edgeColumn(getColumn (v1,0),getColumn(v2,0));
- else
- return edgeColumn(getColumn (v1,0),getColumn(v2,N2-1));
- }
- if (it2 != _fans.end() && it1 != _fans.end() ){
- int c1, c2;
- if (aaa(it1->second._e1,e))
- c1 = 0;
- else
- c1 = N1-1;
- if (aaa(it2->second._e1,e))
- c2 = 0;
- else
- c2 = N2-1;
- return edgeColumn(getColumn (v1,c1),getColumn(v2,c2));
- }
-
- if (N1 == 1 || N2 == 2){
- const BoundaryLayerData & c10 = getColumn(v1,0);
- const BoundaryLayerData & c20 = getColumn(v2,0);
- const BoundaryLayerData & c21 = getColumn(v2,1);
- if (dot(c10._n,c20._n) > dot(c10._n,c21._n) ) return edgeColumn(c10,c20);
- else return edgeColumn(c10,c21);
- }
- if (N1 == 2 || N2 == 1){
- const BoundaryLayerData & c10 = getColumn(v1,0);
- const BoundaryLayerData & c11 = getColumn(v1,1);
- const BoundaryLayerData & c20 = getColumn(v2,0);
- if (dot(c10._n,c20._n) > dot(c11._n,c20._n) ) return edgeColumn(c10,c20);
- else return edgeColumn(c11,c20);
- }
-
- Msg::Error ("Impossible Boundary Layer Configuration : one side and no fans %d %d",N1,N2);
- // FIXME WRONG
- return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) :
- edgeColumn (getColumn (v2,0),getColumn(v2,0));
- }
- else if (nbSides == 2){
- if (it1 == _fans.end() && it2 == _fans.end() ){
- if (N1 != 2 || N2 != 2){
- Msg::Error ("Impossible Boundary Layer Configuration");
- // FIXME WRONG
- return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) :
- edgeColumn (getColumn (v2,0),getColumn(v2,0));
- }
- const BoundaryLayerData & c10 = getColumn(v1,0);
- const BoundaryLayerData & c11 = getColumn(v1,1);
- const BoundaryLayerData & c20 = getColumn(v2,0);
- const BoundaryLayerData & c21 = getColumn(v2,1);
- if (side == 0){
- if (dot(c10._n,c20._n) > dot(c10._n,c21._n) ) return edgeColumn(c10,c20);
- else return edgeColumn(c10,c21);
- }
- if (side == 1){
- if (dot(c11._n,c20._n) > dot(c11._n,c21._n) ) return edgeColumn(c11,c20);
- else return edgeColumn(c11,c21);
- }
- }
- }
-
- Msg::Error("Not yet Done in BoundaryLayerData");
- static BoundaryLayerData error;
- static edgeColumn error2(error, error);
- return error2;
- }
- inline int getNbColumns(MVertex *v) { return _data.count(v); }
- inline const BoundaryLayerData &getColumn(MVertex *v, int iColumn)
- {
- int count = 0;
- for(std::multimap<MVertex*,BoundaryLayerData>::iterator itm = _data.lower_bound(v);
- itm != _data.upper_bound(v); ++itm){
- if (count++ == iColumn) return itm->second;
- }
- static BoundaryLayerData error;
- return error;
- }
- void filterPoints();
-};
-
-BoundaryLayerColumns * buildAdditionalPoints2D (GFace *gf) ;
-BoundaryLayerColumns * buildAdditionalPoints3D (GRegion *gr) ;
-void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]);
-
-#endif
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 4e78b3f24e..3d6203a167 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -9,7 +9,7 @@
#include "GmshMessage.h"
#include "meshGRegion.h"
#include "meshGFaceOptimize.h"
-#include "meshGFaceBoundaryLayers.h"
+#include "boundaryLayersData.h"
#include "meshGRegionDelaunayInsertion.h"
#include "GModel.h"
#include "GRegion.h"
@@ -484,7 +484,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
{
- buildAdditionalPoints3D (gr);
+ // buildAdditionalPoints3D (gr);
}
void _relocateVertex(MVertex *ver,
@@ -959,7 +959,7 @@ void meshGRegion::operator() (GRegion *gr)
}
// replace discreteFaces by their compounds
- if(gr->geomType() == GEntity::CompoundVolume){
+ if( 1 || gr->geomType() == GEntity::CompoundVolume){
std::set<GFace*> mySet;
std::list<GFace*>::iterator it = faces.begin();
while(it != faces.end()){
diff --git a/Mesh/meshGRegionBoundaryLayers.cpp b/Mesh/meshGRegionBoundaryLayers.cpp
index 315b74d70e..2310bcd56e 100644
--- a/Mesh/meshGRegionBoundaryLayers.cpp
+++ b/Mesh/meshGRegionBoundaryLayers.cpp
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
-#include "meshGFaceBoundaryLayers.h"
+#include "boundaryLayersData.h"
#include "GModel.h"
#include "GRegion.h"
#include "MTriangle.h"
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index f7ebac6393..77f2469970 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -402,19 +402,19 @@ bool insertVertexB(std::list<faceXtet> &shell,
while (it != shell.end()){
MTetrahedron *tr = new MTetrahedron(it->getVertex(0), it->getVertex(1), it->getVertex(2), v);
- double lc = .25 * (vSizes[tr->getVertex(0)->getIndex()] +
- vSizes[tr->getVertex(1)->getIndex()] +
- vSizes[tr->getVertex(2)->getIndex()] +
- vSizes[tr->getVertex(3)->getIndex()]);
- double lcBGM = .25 * (vSizesBGM[tr->getVertex(0)->getIndex()] +
- vSizesBGM[tr->getVertex(1)->getIndex()] +
- vSizesBGM[tr->getVertex(2)->getIndex()] +
- vSizesBGM[tr->getVertex(3)->getIndex()]);
- double LL = std::min(lc, lcBGM);
+ // double lc = .25 * (vSizes[tr->getVertex(0)->getIndex()] +
+ // vSizes[tr->getVertex(1)->getIndex()] +
+ // vSizes[tr->getVertex(2)->getIndex()] +
+ // vSizes[tr->getVertex(3)->getIndex()]);
+ // double lcBGM = .25 * (vSizesBGM[tr->getVertex(0)->getIndex()] +
+ // vSizesBGM[tr->getVertex(1)->getIndex()] +
+ // vSizesBGM[tr->getVertex(2)->getIndex()] +
+ // vSizesBGM[tr->getVertex(3)->getIndex()]);
+ // double LL = std::min(lc, lcBGM);
MTet4 *t4 = myFactory.Create(tr, vSizes, vSizesBGM);
t4->setOnWhat(t->onWhat());
-
+ /*
double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
(it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
(it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
@@ -426,6 +426,7 @@ bool insertVertexB(std::list<faceXtet> &shell,
(it->v[2]->z() - v->z()) * (it->v[2]->z() - v->z()));
if (d1 < LL * .25 || d2 < LL * .25 || d3 < LL * .25) onePointIsTooClose = true;
+ */
newTets[k++] = t4;
// all new tets are pushed front in order to ba able to destroy
// them if the cavity is not star shaped around the new vertex.
@@ -1282,7 +1283,7 @@ void insertVerticesInRegion (GRegion *gr)
// Normally, a tet mesh contains about 6 times more tets than
// vertices. This allows to clean up the set of tets when lots of
// deleted ones are present in the mesh
- if(allTets.size() > 7 * vSizes.size() && ITER > 20000){
+ if(allTets.size() > 7 * vSizes.size() && ITER > 1000){
memoryCleanup(myFactory, allTets);
}
}
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index f547098ba3..c1c7a170a1 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -181,6 +181,10 @@ bool compute4neighbors (GFace *gf, // the surface
SPoint2 newP[4][NUMDIR], // look into other directions
SMetric3 &metricField, FILE *crossf = 0) // the mesh metric
{
+
+ Range<double> rangeU = gf->parBounds(0);
+ Range<double> rangeV = gf->parBounds(1);
+
// we assume that v is on surface gf
// get the parameter of the point on the surface
@@ -247,14 +251,17 @@ bool compute4neighbors (GFace *gf, // the surface
// t1 . s2 = a E + b N --> solve the 2 x 2 system
// and get covariant coordinates a and b
double rhs1[2] = {dot(t1,s1),dot(t1,s2)}, covar1[2];
+ bool singular = false;
if (!sys2x2(metric,rhs1,covar1)){
- Msg::Info("Argh");
+ Msg::Info("Argh surface %d %g %g %g -- %g %g %g -- %g %g",gf->tag(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z(),size_1,size_2);
covar1[1] = 1.0; covar1[0] = 0.0;
+ singular = true;
}
double rhs2[2] = {dot(t2,s1),dot(t2,s2)}, covar2[2];
if (!sys2x2(metric,rhs2,covar2)){
- Msg::Info("Argh");
+ Msg::Info("Argh surface %d %g %g %g -- %g %g %g",gf->tag(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z());
covar2[0] = 1.0; covar2[1] = 0.0;
+ singular = true;
}
// transform the sizes with respect to the metric
@@ -263,18 +270,22 @@ bool compute4neighbors (GFace *gf, // the surface
// of size1 in direction v, it should be sqrt(v^T M v) * size1
double l1 = sqrt(covar1[0]*covar1[0]+covar1[1]*covar1[1]);
double l2 = sqrt(covar2[0]*covar2[0]+covar2[1]*covar2[1]);
-
covar1[0] /= l1;covar1[1] /= l1;
covar2[0] /= l2;covar2[1] /= l2;
-
- const double size_param_1 = size_1 / sqrt ( M*covar1[0]*covar1[0]+
+ double size_param_1 = size_1 / sqrt ( M*covar1[0]*covar1[0]+
2*E*covar1[1]*covar1[0]+
N*covar1[1]*covar1[1]);
- const double size_param_2 = size_2 / sqrt ( M*covar2[0]*covar2[0]+
- 2*E*covar2[1]*covar2[0]+
+ double size_param_2 = size_2 / sqrt ( M*covar2[0]*covar2[0]+
+ 2*E*covar2[1]*covar2[0]+
N*covar2[1]*covar2[1]);
+ if (singular){
+ size_param_1 = size_param_2 = std::min (size_param_1,size_param_2);
+ }
+
+ // printf("%12.5E %12.5E\n", size_param_1, size_param_2);
+
// if (v_center->onWhat() != gf && gf->tag() == 3)
// printf("M = (%g %g %g) L = %g %g LP = %g %g\n",metricField(0,0),metricField(1,1),metricField(0,1),l1,l2,size_param_1,size_param_2);
@@ -312,6 +323,10 @@ bool compute4neighbors (GFace *gf, // the surface
// surface
double ERR[4];
for (int i=0;i<4;i++){ //
+ // if (newPoint[i][0] < rangeU.low())newPoint[i][0] = rangeU.low();
+ // if (newPoint[i][0] > rangeU.high())newPoint[i][0] = rangeU.high();
+ // if (newPoint[i][1] < rangeV.low())newPoint[i][1] = rangeV.low();
+ // if (newPoint[i][1] > rangeV.high())newPoint[i][1] = rangeV.high();
GPoint pp = gf->point(SPoint2(newPoint[i][0], newPoint[i][1]));
double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
(pp.y() - v_center->y())*(pp.y() - v_center->y()) +
@@ -320,7 +335,7 @@ bool compute4neighbors (GFace *gf, // the surface
// printf("L = %12.5E D = %12.5E ERR = %12.5E\n",L,D,100*fabs(D-L)/(D+L));
}
- if (0 && goNonLinear){//---------------------------------------------------//
+ if (1 && goNonLinear){//---------------------------------------------------//
surfaceFunctorGFace ss (gf); //
SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1*(1.0),t2*(1.0)}; //
for (int i=0;i<4;i++){ //
@@ -330,7 +345,7 @@ bool compute4neighbors (GFace *gf, // the surface
curveFunctorCircle cf (dirs[i],n,
SVector3(v_center->x(),v_center->y(),v_center->z()),
L);
- if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-5)){ //
+ if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-3)){ //
GPoint pp = gf->point(SPoint2(uvt[0],uvt[1]));
double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
(pp.y() - v_center->y())*(pp.y() - v_center->y()) +
@@ -346,10 +361,12 @@ bool compute4neighbors (GFace *gf, // the surface
// printf("%12.5E vs %12.5E : %12.5E %12.5E vs %12.5E %12.5E \n",ERR[i],newErr,newPoint[i][0],newPoint[i][1],uvt[0],uvt[1]);
newPoint[i][0] = uvt[0]; //
newPoint[i][1] = uvt[1]; //
- } //
+ } //
+ // printf("OK\n");
}
else{
Msg::Debug("Cannot put a new point on Surface %d",gf->tag());
+ // printf("NOT OK\n");
}
}
} //
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 598fe22729..5a1afbe259 100755
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -338,17 +338,17 @@ void Recombinator::execute(GRegion* gr){
build_tuples(gr);
init_markings(gr);
+ Msg::Info("Building Connectivity...");
build_vertex_to_vertices(gr);
build_vertex_to_elements(gr);
- printf("connectivity\n");
potential.clear();
patern1(gr);
- printf("patern no. 1\n");
+ Msg::Info("Hex-merging patern nb. 1...");
patern2(gr);
- printf("patern no. 2\n");
+ Msg::Info("Hex-merging patern nb. 2...");
patern3(gr);
- printf("patern no. 3\n");
+ Msg::Info("Hex-merging patern nb. 3...");
std::sort(potential.begin(),potential.end());
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 1a124c83c5..0dd9bc3f03 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -684,10 +684,13 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
func(x, feps, data);
for (int i = 0; i < N; i++){
J(i, j) = (feps(i) - f(i)) / h;
+ // printf("J(%d,%d) = %12.5E \n",i,j,J(i,j));
}
x(j) -= h;
}
+ // printf("done\n");
+
if (N == 1)
dx(0) = f(0) / J(0, 0);
else
@@ -697,6 +700,7 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
for (int i = 0; i < N; i++)
x(i) -= relax * dx(i);
+ // printf("dx = %12.5E\n",dx.norm());
if(dx.norm() < tolx) return true;
}
return false;
diff --git a/Numeric/simpleFunction.h b/Numeric/simpleFunction.h
index 61754fc94a..313e2139f7 100644
--- a/Numeric/simpleFunction.h
+++ b/Numeric/simpleFunction.h
@@ -15,7 +15,7 @@ class simpleFunction {
scalar _val;
bool _hasDerivatives;
public :
- simpleFunction(scalar val=0.0) : _val(val), _hasDerivatives(false){}
+ simpleFunction(scalar val =0.0) : _val(val), _hasDerivatives(false){}
virtual ~simpleFunction(){}
virtual bool hasDerivatives() {return _hasDerivatives;};
virtual scalar operator () (double x, double y, double z) const { return _val; }
diff --git a/benchmarks/3d/submarine/submarine_simple.geo b/benchmarks/3d/submarine/submarine_simple.geo
index 2873ebf3a9..d00489aa7e 100644
--- a/benchmarks/3d/submarine/submarine_simple.geo
+++ b/benchmarks/3d/submarine/submarine_simple.geo
@@ -1,11 +1,24 @@
-lc = 200;
+lc = 25;
Mesh.CharacteristicLengthMin = lc;
Mesh.CharacteristicLengthMax = lc;
+Mesh.Smoothing = 0;
-Merge "submarine.stp";
+Merge "submarine.brep";
+
+
+Compound Surface(1034) = {12};
Physical Surface(1) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
-Physical Surface(2) = {12};
+Physical Surface(2) = {1034};
Physical Surface(3) = {11};
-Physical Volume(4) = {1};
+MeshAlgorithm Surface {1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,1034} 9;
+MeshAlgorithm Surface {12} 5;
+Mesh.Recombine3DAll = 1;
+Mesh.Algorithm3D = 9;
+Delete {Volume {1};}
+//Compound Volume(1035) = {1};
+
+Surface Loop(1000) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1034};
+Volume(1001) = {1000};
+Physical Volume(1) = {1001};
diff --git a/wrappers/gmshpy/gmshGeo.i b/wrappers/gmshpy/gmshGeo.i
index 6390afe435..512247447d 100644
--- a/wrappers/gmshpy/gmshGeo.i
+++ b/wrappers/gmshpy/gmshGeo.i
@@ -37,6 +37,7 @@
#include "SBoundingBox3d.h"
#include "Curvature.h"
#include "simpleFunction.h"
+ #include "distanceToMesh.h"
#include "GeomMeshMatcher.h"
%}
@@ -95,6 +96,7 @@ namespace std {
%include "SBoundingBox3d.h"
%include "Curvature.h"
%include "gmshLevelset.h"
+%include "distanceToMesh.h"
%include "GeomMeshMatcher.h"
--
GitLab