diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index bad1dd74008288dd212d79e4926b1b2fb0e40d8a..fabb515d453be6e03a9fb45c91503d592b04d3ac 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -532,6 +532,10 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.lcFromCurvature = 1;
i++;
}
+ else if(!strcmp(argv[i] + 1, "clcurviso")) {
+ CTX::instance()->mesh.lcFromCurvature = 2;
+ i++;
+ }
else if(!strcmp(argv[i] + 1, "smooth")) {
i++;
if(argv[i])
diff --git a/Common/Context.h b/Common/Context.h
index 8fcab24dfc6b0b6e9cfb04d70c19ed521efcfacb..39918459273b46f0f3300599f814d995accc10e5 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -42,6 +42,7 @@ struct contextMeshOptions {
int multiplePasses;
int cgnsImportOrder;
std::map<int,int> algo2d_per_face;
+ std::map<int,int> curvature_control_per_face;
int bunin;
int ignorePartBound;
};
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 6e8f8b42445bc3e0e8deb51adfe9fa97593505d1..bc24b569f5f53f170f2520dc639238ec7976f2b5 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -62,6 +62,19 @@ GFace::~GFace()
deleteMesh();
}
+int GFace::getCurvatureControlParameter () const
+{
+ std::map<int,int>::iterator it = CTX::instance()->mesh.curvature_control_per_face.find(tag());
+ return it == CTX::instance()->mesh.curvature_control_per_face.end() ?
+ CTX::instance()->mesh.minCircPoints : it->second ;
+}
+
+void GFace::setCurvatureControlParameter (int n)
+{
+ CTX::instance()->mesh.curvature_control_per_face[tag()] = n;
+}
+
+
int GFace::getMeshingAlgo () const
{
std::map<int,int>::iterator it = CTX::instance()->mesh.algo2d_per_face.find(tag());
@@ -1220,36 +1233,54 @@ int GFace::genusGeom()
return nSeams - single_seams.size();
}
-bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
+bool GFace::fillPointCloud(double maxDist,
+ std::vector<SPoint3> *points,
+ std::vector<SPoint2> *uvpoints,
std::vector<SVector3> *normals)
{
- if(!buildSTLTriangulation()){
- Msg::Error("No STL triangulation available to fill point cloud");
- return false;
- }
if(!points) return false;
- for(unsigned int i = 0; i < stl_triangles.size(); i += 3){
- SPoint2 &p0(stl_vertices[stl_triangles[i]]);
- SPoint2 &p1(stl_vertices[stl_triangles[i + 1]]);
- SPoint2 &p2(stl_vertices[stl_triangles[i + 2]]);
-
- GPoint gp0 = point(p0);
- GPoint gp1 = point(p1);
- GPoint gp2 = point(p2);
- double maxEdge = std::max(gp0.distance(gp1),
- std::max(gp1.distance(gp2), gp2.distance(gp0)));
- int N = (int)(maxEdge / maxDist);
- for(double u = 0.; u < 1.; u += 1. / N){
- for(double v = 0.; v < 1 - u; v += 1. / N){
- SPoint2 p = p0 * (1. - u - v) + p1 * u + p2 * v;
- GPoint gp(point(p));
- points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
- if(normals) normals->push_back(normal(p));
+ if (buildSTLTriangulation()){
+ for(unsigned int i = 0; i < stl_triangles.size(); i += 3){
+ SPoint2 &p0(stl_vertices[stl_triangles[i]]);
+ SPoint2 &p1(stl_vertices[stl_triangles[i + 1]]);
+ SPoint2 &p2(stl_vertices[stl_triangles[i + 2]]);
+
+ GPoint gp0 = point(p0);
+ GPoint gp1 = point(p1);
+ GPoint gp2 = point(p2);
+ double maxEdge = std::max(gp0.distance(gp1),
+ std::max(gp1.distance(gp2), gp2.distance(gp0)));
+ int N = (int)(maxEdge / maxDist);
+ for(double u = 0.; u < 1.; u += 1. / N){
+ for(double v = 0.; v < 1 - u; v += 1. / N){
+ SPoint2 p = p0 * (1. - u - v) + p1 * u + p2 * v;
+ GPoint gp(point(p));
+ points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
+ if (uvpoints)uvpoints->push_back(p);
+ if(normals) normals->push_back(normal(p));
+ }
}
}
}
+ else {
+ int N = 1000;//(int)(maxDX / maxDist);
+ Range<double> b1 = parBounds(0);
+ Range<double> b2 = parBounds(1);
+ for(int i = 0; i < N; i++) {
+ for(int j = 0; j < N; j++) {
+ double u = (double)i / (N - 1);
+ double v = (double)j / (N - 1);
+ double t1 = b1.low() + u * (b1.high() - b1.low());
+ double t2 = b2.low() + v * (b2.high() - b2.low());
+ GPoint gp = point(t1, t2);
+ points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
+ if (uvpoints)uvpoints->push_back(SPoint2(t1,t2));
+ if(normals) normals->push_back(normal(SPoint2(t1,t2)));
+ }
+ }
+ }
return true;
}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index f17eca36b71ce5214928db653838cdbc445aaf96..282b4bc92d483c0505288616371ae86c65a04345 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -250,7 +250,9 @@ class GFace : public GEntity
// add points (and optionally normals) in vectors so that two
// points are at most maxDist apart
- bool fillPointCloud(double maxDist, std::vector<SPoint3> *points,
+ bool fillPointCloud(double maxDist,
+ std::vector<SPoint3> *points,
+ std::vector<SPoint2> *uvpoints,
std::vector<SVector3> *normals=0);
// apply Lloyd's algorithm to the mesh
@@ -286,6 +288,8 @@ class GFace : public GEntity
int getMeshingAlgo () const;
void setMeshingAlgo (int) ;
+ int getCurvatureControlParameter () const;
+ void setCurvatureControlParameter (int) ;
struct {
mutable GEntity::MeshGenerationStatus status;
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 172f1e1b350540d50bd6341d3ffab872b57d2515..4eb29773b226f7dd83d201ef385099d86e7a91fd 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -14,6 +14,7 @@
#include "Context.h"
#if defined(HAVE_OCC)
+#include <TColStd_MapIteratorOfMapOfAsciiString.hxx>
#include <Standard_Version.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_ConicalSurface.hxx>
@@ -25,6 +26,8 @@
#include <Geom_Plane.hxx>
#include <gp_Pln.hxx>
#include <BRepMesh_FastDiscret.hxx>
+#include <BRepMesh_PDiscretRoot.hxx>
+#include <BRepMesh_DiscretFactory.hxx>
#include <IntTools_Context.hxx>
#include <BOPTools_Tools2D.hxx>
#include <BOPTools_Tools3D.hxx>
@@ -336,16 +339,37 @@ bool OCCFace::buildSTLTriangulation(bool force)
Bnd_Box aBox;
BRepBndLib::Add(s, aBox);
- BRepMesh_FastDiscret aMesher(0.1, 0.5, aBox, Standard_False, Standard_False,
- Standard_True, Standard_False);
-#if (OCC_VERSION_MAJOR == 6) && (OCC_VERSION_MINOR < 5)
- aMesher.Add(s);
-#else
- aMesher.Perform(s);
-#endif
+
+ Standard_Real aDiscret, aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
+ Standard_Real dX, dY, dZ, dMax, aCoeff, aAngle = .5;
+ aBox.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
+ dX=aXmax-aXmin;
+ dY=aYmax-aYmin;
+ dZ=aZmax-aZmin;
+ dMax=dX;
+ if (dY>dMax) {
+ dMax=dY;
+ }
+ if (dZ>dMax) {
+ dMax=dZ;
+ }
+ //
+ aCoeff=0.1;
+ aDiscret=aCoeff*dMax;
+
+ // std::cout << BRepMesh_DiscretFactory::Get().Names() << "\n";
+ // const TColStd_MapOfAsciiString &aaa = BRepMesh_DiscretFactory::Get().Names();
+ // TColStd_MapIteratorOfMapOfAsciiString It (aaa);
+ // for (; It.More(); It.Next()) {
+ // std::cout << It.Key() << "\n";
+ // }
+
+ BRepMesh_PDiscretRoot pAlgo = BRepMesh_DiscretFactory::Get().Discret(s, aDiscret, aAngle);
+ if (pAlgo) pAlgo->Perform();
TopLoc_Location loc;
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation(s, loc);
+
if(triangulation.IsNull() || !triangulation->HasUVNodes()){
if(triangulation.IsNull())
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 9793c2eca60372bbde04d97a3d43551e9e7c3cf6..ae0cf3e6cf069baadfc49b240ef98a5e143a4fb0 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -71,7 +71,7 @@ SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2
double ratio1 = lambda1_min/lambda1_max;
double ratio2 = lambda2_min/lambda2_max;
- if (ratio1<ratio2)
+ if (lambda1_min<lambda2_min)
return intersection_conserveM1(m1,m2);
else
return intersection_conserveM1(m2,m1);
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 742377dc654a543510d4e1c40e4a6099cc12f220..7a70f6094ceb8d8be5aae74cd59719badb05b6fc 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -35,10 +35,9 @@
// CTX::instance()->mesh.minCircPoints tells the minimum number of points per
// radius of curvature
-SMetric3 buildMetricTangentToCurve (SVector3 &t, double curvature, double &lambda)
+SMetric3 buildMetricTangentToCurve (SVector3 &t, double l_t, double l_n)
{
- lambda = 1.e22;
- if (curvature == 0.0)return SMetric3(1.e-22);
+ if (l_t == 0.0)return SMetric3(1.e-22);
SVector3 a;
if (t(0) <= t(1) && t(0) <= t(2)){
a = SVector3(1,0,0);
@@ -54,45 +53,59 @@ SMetric3 buildMetricTangentToCurve (SVector3 &t, double curvature, double &lambd
b.normalize();
c.normalize();
t.normalize();
- lambda = ((2 * M_PI) /( fabs(curvature) * CTX::instance()->mesh.minCircPoints ) );
-
- SMetric3 curvMetric (1./(lambda*lambda),1.e-10,1.e-10,t,b,c);
- //SMetric3 curvMetric (1./(lambda*lambda));
- return curvMetric;
+ SMetric3 Metric (1./(l_t*l_t),1./(l_n*l_n),1./(l_n*l_n),t,b,c);
+ return Metric;
+}
+
+SMetric3 buildMetricTangentToSurface (SVector3 &t1, SVector3 &t2, double l_t1, double l_t2, double l_n)
+{
+ t1.normalize();
+ t2.normalize();
+ SVector3 n = crossprod (t1,t2);
+ n.normalize();
+ SMetric3 Metric (1./(l_t1*l_t1),1./(l_t2*l_t2),1./(l_n*l_n),t1,t2,n);
+ return Metric;
}
-SMetric3 max_edge_curvature_metric(const GVertex *gv, double &length)
+
+SMetric3 max_edge_curvature_metric(const GVertex *gv)
{
SMetric3 val (1.e-12);
- length = 1.e22;
std::list<GEdge*> l_edges = gv->edges();
for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
ite != l_edges.end(); ++ite){
GEdge *_myGEdge = *ite;
Range<double> range = _myGEdge->parBounds(0);
SMetric3 cc;
- double l;
if (gv == _myGEdge->getBeginVertex()) {
SVector3 t = _myGEdge->firstDer(range.low());
t.normalize();
- cc = buildMetricTangentToCurve(t,_myGEdge->curvature(range.low()),l);
+ double l_t = ((2 * M_PI) /( fabs(_myGEdge->curvature(range.low()))
+ * CTX::instance()->mesh.minCircPoints ));
+ double l_n = 1.e12;
+ cc = buildMetricTangentToCurve(t,l_t,l_n);
}
else {
SVector3 t = _myGEdge->firstDer(range.high());
t.normalize();
- cc = buildMetricTangentToCurve(t,_myGEdge->curvature(range.high()),l);
+ double l_t = ((2 * M_PI) /( fabs(_myGEdge->curvature(range.high()))
+ * CTX::instance()->mesh.minCircPoints ));
+ double l_n = 1.e12;
+ cc = buildMetricTangentToCurve(t,l_t,l_n);
}
val = intersection(val,cc);
- length = std::min(length,l);
}
return val;
}
-SMetric3 max_edge_curvature_metric(const GEdge *ge, double u, double &l)
+SMetric3 max_edge_curvature_metric(const GEdge *ge, double u)
{
SVector3 t = ge->firstDer(u);
t.normalize();
- return buildMetricTangentToCurve(t,ge->curvature(u),l);
+ double l_t = ((2 * M_PI) /( fabs(ge->curvature(u))
+ * CTX::instance()->mesh.minCircPoints ));
+ double l_n = 1.e12;
+ return buildMetricTangentToCurve(t,l_t,l_n);
}
static double max_edge_curvature(const GVertex *gv)
@@ -144,35 +157,40 @@ static double max_surf_curvature(const GVertex *gv)
return val;
}
-static SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v)
+SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
+ bool surface_isotropic,
+ double d_normal ,
+ double d_tangent_max )
{
- if (gf->geomType() == GEntity::Plane)return SMetric3(1.e-6);
+ if (gf->geomType() == GEntity::Plane)return SMetric3(1.e-12);
double cmax, cmin;
SVector3 dirMax,dirMin;
cmax = gf->curvatures(SPoint2(u, v),&dirMax, &dirMin, &cmax,&cmin);
- if (cmin == 0)cmin =1.e-5;
- if (cmax == 0)cmax =1.e-5;
+ if (cmin == 0)cmin =1.e-12;
+ if (cmax == 0)cmax =1.e-12;
double lambda2 = ((2 * M_PI) /( fabs(cmax) * CTX::instance()->mesh.minCircPoints ) );
double lambda1 = ((2 * M_PI) /( fabs(cmin) * CTX::instance()->mesh.minCircPoints ) );
SVector3 Z = crossprod(dirMax,dirMin);
-
+ if (surface_isotropic) lambda2 = lambda1 = std::min(lambda2,lambda1);
+ dirMin.normalize();
+ dirMax.normalize();
+ Z.normalize();
lambda1 = std::max(lambda1, CTX::instance()->mesh.lcMin);
lambda2 = std::max(lambda2, CTX::instance()->mesh.lcMin);
lambda1 = std::min(lambda1, CTX::instance()->mesh.lcMax);
lambda2 = std::min(lambda2, CTX::instance()->mesh.lcMax);
- /* if (gf->tag() == 36 || gf->tag() == 62)
- printf("%g %g -- %g %g -- %g %g %g -- %g %g %g -- %g %g %g -- %g %g\n",u,v,cmax,cmin,
- dirMax.x(),dirMax.y(),dirMax.z(),
- dirMin.x(),dirMin.y(),dirMin.z(),
- Z.x(),Z.y(),Z.z(),
- lambda1,lambda2);
- */
- SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),1.e-5,
+ double lambda3 = std::min(d_normal, CTX::instance()->mesh.lcMax);
+ lambda3 = std::max(lambda3, CTX::instance()->mesh.lcMin);
+ lambda1 = std::min(lambda1, d_tangent_max);
+ lambda2 = std::min(lambda2, d_tangent_max);
+
+ SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),
+ 1./(lambda3*lambda3),
dirMin, dirMax, Z );
return curvMetric;
}
-static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
+static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u, bool iso_surf)
{
const GEdgeCompound* ptrCompoundEdge = dynamic_cast<const GEdgeCompound*>(ge);
if (ptrCompoundEdge)
@@ -200,30 +218,27 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
SMetric3 mesh_size(1.e-12);
std::list<GFace *> faces = ge->faces();
std::list<GFace *>::iterator it = faces.begin();
- int count = 0;
+ // we choose the metric eigenvectors to be the ones
+ // related to the edge ...
+ SMetric3 curvMetric = max_edge_curvature_metric(ge, u);
while(it != faces.end())
- {
- if (((*it)->geomType() != GEntity::CompoundSurface) &&
- ((*it)->geomType() != GEntity::DiscreteSurface)){
- SPoint2 par = ge->reparamOnFace((*it), u, 1);
- SMetric3 m = metric_based_on_surface_curvature (*it, par.x(), par.y());
- if (!count) mesh_size = m;
- else mesh_size = intersection(mesh_size, m);
- count++;
- }
- ++it;
- }
- double Crv = ge->curvature(u);
- double lambda = ((2 * M_PI) /( fabs(Crv) * CTX::instance()->mesh.minCircPoints ) );
- SMetric3 curvMetric (1./(lambda*lambda));
-
- return intersection_conserve_mostaniso(mesh_size,curvMetric);
+ {
+ if (((*it)->geomType() != GEntity::CompoundSurface) &&
+ ((*it)->geomType() != GEntity::DiscreteSurface)){
+ SPoint2 par = ge->reparamOnFace((*it), u, 1);
+ SMetric3 m = metric_based_on_surface_curvature (*it, par.x(), par.y(), iso_surf);
+ curvMetric = intersection_conserveM1(curvMetric,m);
+ }
+ ++it;
+ }
+
+ return curvMetric;
}
}
-static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
+static SMetric3 metric_based_on_surface_curvature(const GVertex *gv, bool iso_surf)
{
- SMetric3 mesh_size(1.e-5);
+ SMetric3 mesh_size(1.e-15);
std::list<GEdge*> l_edges = gv->edges();
for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
ite != l_edges.end(); ++ite){
@@ -239,11 +254,11 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
if (gv == _myGEdge->getBeginVertex())
mesh_size = intersection
(mesh_size,
- metric_based_on_surface_curvature(_myGEdge, bounds.low()));
+ metric_based_on_surface_curvature(_myGEdge, bounds.low(), iso_surf));
else
mesh_size = intersection
(mesh_size,
- metric_based_on_surface_curvature(_myGEdge, bounds.high()));
+ metric_based_on_surface_curvature(_myGEdge, bounds.high(), iso_surf));
}
}
@@ -262,15 +277,15 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
switch(ge->dim()){
case 0:
Crv = max_edge_curvature((const GVertex *)ge);
- Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
- //Crv = max_surf_curvature((const GVertex *)ge);
+ //Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
+ // Crv = max_surf_curvature((const GVertex *)ge);
break;
case 1:
{
GEdge *ged = (GEdge *)ge;
Crv = ged->curvature(U);
Crv = std::max(Crv, max_surf_curvature(ged, U));
- //Crv = max_surf_curvature(ged, U);
+ // Crv = max_surf_curvature(ged, U);
}
break;
case 2:
@@ -284,12 +299,14 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
return lc;
}
-static SMetric3 LC_MVertex_CURV_ANISO(GEntity *ge, double U, double V)
+SMetric3 LC_MVertex_CURV_ANISO(GEntity *ge, double U, double V)
{
+ bool iso_surf = CTX::instance()->mesh.lcFromCurvature == 2;
+
switch(ge->dim()){
- case 0: return metric_based_on_surface_curvature((const GVertex *)ge);
- case 1: return metric_based_on_surface_curvature((const GEdge *)ge, U);
- case 2: return metric_based_on_surface_curvature((const GFace *)ge, U, V);
+ case 0: return metric_based_on_surface_curvature((const GVertex *)ge, iso_surf);
+ case 1: return metric_based_on_surface_curvature((const GEdge *)ge, U, iso_surf);
+ case 2: return metric_based_on_surface_curvature((const GFace *)ge, U, V, iso_surf);
}
Msg::Error("Curvature control impossible to compute for a volume!");
return SMetric3();
@@ -382,19 +399,20 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
{
// default lc (mesh size == size of the model)
double l1 = CTX::instance()->lc;
+ const double max_lc = CTX::instance()->mesh.lcMax;
// lc from points
- double l2 = MAX_LC;
+ double l2 = max_lc;
if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
l2 = LC_MVertex_PNTS(ge, U, V);
// lc from curvature
- SMetric3 l3(1./(MAX_LC*MAX_LC));
+ SMetric3 l3(1./(max_lc*max_lc));
if(CTX::instance()->mesh.lcFromCurvature && ge->dim() < 3)
l3 = LC_MVertex_CURV_ANISO(ge, U, V);
// lc from fields
- SMetric3 l4(1./(MAX_LC*MAX_LC));
+ SMetric3 l4(1./(max_lc*max_lc));
FieldManager *fields = ge->model()->getFields();
if(fields->getBackgroundField() > 0){
Field *f = fields->get(fields->getBackgroundField());
@@ -426,8 +444,22 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
}
SMetric3 LC(1./(lc*lc));
- SMetric3 m = intersection(intersection (l4, LC),l3);
+ SMetric3 m = intersection_conserveM1(intersection_conserveM1 (l4, LC),l3);
//printf("%g %g %g %g %g %g\n",m(0,0),m(1,1),m(2,2),m(0,1),m(0,2),m(1,2));
+ {
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ m.eig(V,S,true);
+ if (S(0) < 0 || S(1) < 0 || S(2) < 0){
+ printf("%g %g %g\n",S(0),S(1),S(2));
+ l3.eig(V,S,true);
+ printf("%g %g %g\n",S(0),S(1),S(2));
+ l4.eig(V,S,true);
+ printf("%g %g %g\n",S(0),S(1),S(2));
+ }
+ }
+
+
return m;
// return lc * CTX::instance()->mesh.lcFactor;
}
@@ -510,7 +542,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
else {
std::map<MVertex*, MVertex*>::iterator itv2 = _2Dto3D.begin();
for ( ; itv2 != _2Dto3D.end(); ++itv2){
- _sizes[itv2->first] = MAX_LC;
+ _sizes[itv2->first] = CTX::instance()->mesh.lcMax;
}
}
// ensure that other criteria are fullfilled
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 8d5030c6856ea392b668e927e8364d37748384c1..e342b40d320bddeb5f842bb6dcd6e7a95225837f 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -79,11 +79,17 @@ class backgroundMesh : public simpleFunction<double>
MElementOctree* get_octree();
};
+SMetric3 buildMetricTangentToCurve (SVector3 &t, double l_t, double l_n);
+SMetric3 buildMetricTangentToSurface (SVector3 &t1, SVector3 &t2, double l_t1, double l_t2, double l_n);
double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double Z);
SMetric3 BGM_MeshMetric(GEntity *ge, double U, double V, double X, double Y, double Z);
bool Extend1dMeshIn2dSurfaces();
bool Extend2dMeshIn3dVolumes();
-SMetric3 max_edge_curvature_metric(const GVertex *gv, double &l);
+SMetric3 max_edge_curvature_metric(const GVertex *gv);
SMetric3 max_edge_curvature_metric(const GEdge *ge, double u, double &l);
+SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
+ bool surface_isotropic = false,
+ double d_normal = 1.e12,
+ double d_tangent_max = 1.e12);
#endif
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index eea3503d7be576375a718f8f26cd12d143ac84b1..ae3ead81827d3acb90233d6bf81c7fe4d0b2d8d2 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -35,7 +35,6 @@
#include "ANN/ANN.h"
#endif
-
Field::~Field()
{
for(std::map<std::string, FieldOption*>::iterator it = options.begin();
@@ -1649,8 +1648,33 @@ class AttractorField : public Field
annDeallocPts(zeronodes);
delete kdtree;
}
- int totpoints = nodes_id.size() + (n_nodes_by_edge-2) * edges_id.size() +
- n_nodes_by_edge * n_nodes_by_edge * faces_id.size();
+
+ std::vector<SPoint3> points;
+ std::vector<SPoint2> uvpoints;
+ std::vector<int> offset;
+ offset.push_back(0);
+ for(std::list<int>::iterator it = faces_id.begin();
+ it != faces_id.end(); ++it) {
+ Surface *s = FindSurface(*it);
+ if(!s) {
+ 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 ;
+ bool success = f->fillPointCloud(maxDist, &points, &uvpoints);
+ offset.push_back(points.size());
+ }
+ }
+ }
+
+
+ int totpoints =
+ nodes_id.size() +
+ (n_nodes_by_edge-2) * edges_id.size() +
+ ((points.size()) ? points.size() : n_nodes_by_edge * n_nodes_by_edge * faces_id.size());
+ printf("%d points found in points clouds (%d edges)\n",totpoints, edges_id.size());
+
if(totpoints){
zeronodes = annAllocPts(totpoints, 3);
_infos.resize(totpoints);
@@ -1676,8 +1700,7 @@ class AttractorField : public Field
for(std::list<int>::iterator it = edges_id.begin();
it != edges_id.end(); ++it) {
Curve *c = FindCurve(*it);
- GEdge *e = GModel::current()->getEdgeByTag(*it);
- if(c && !e) {
+ if(c) {
for(int i = 1; i < n_nodes_by_edge -1 ; i++) {
double u = (double)i / (n_nodes_by_edge - 1);
Vertex V = InterpolateCurve(c, u, 0);
@@ -1687,6 +1710,7 @@ class AttractorField : public Field
}
}
else {
+ GEdge *e = GModel::current()->getEdgeByTag(*it);
if(e) {
for(int i = 1; i < n_nodes_by_edge - 1; i++) {
double u = (double)i / (n_nodes_by_edge - 1);
@@ -1703,6 +1727,7 @@ class AttractorField : public Field
// This can lead to weird results as we generate attractors over
// the whole parametric plane (we should really use a mesh,
// e.g. a refined STL.)
+ int count = 0;
for(std::list<int>::iterator it = faces_id.begin();
it != faces_id.end(); ++it) {
Surface *s = FindSurface(*it);
@@ -1720,22 +1745,33 @@ class AttractorField : public Field
}
else {
GFace *f = GModel::current()->getFaceByTag(*it);
- if(f) {
- for(int i = 0; i < n_nodes_by_edge; i++) {
- for(int j = 0; j < n_nodes_by_edge; j++) {
- double u = (double)i / (n_nodes_by_edge - 1);
- double v = (double)j / (n_nodes_by_edge - 1);
- Range<double> b1 = ge->parBounds(0);
- Range<double> b2 = ge->parBounds(1);
- double t1 = b1.low() + u * (b1.high() - b1.low());
- double t2 = b2.low() + v * (b2.high() - b2.low());
- GPoint gp = f->point(t1, t2);
- getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
- zeronodes[k][1], zeronodes[k][2], f);
- _infos[k++] = AttractorInfo(*it,2,u,v);
- }
+ if(f) {
+ if (points.size()){
+ for(int j = offset[count]; j < offset[count+1];j++) {
+ zeronodes[k][0] = points[j].x();
+ zeronodes[k][1] = points[j].y();
+ zeronodes[k][2] = points[j].z();
+ _infos[k++] = AttractorInfo(*it,2,uvpoints[j].x(),uvpoints[j].y());
+ }
+ count++;
+ }
+ else{
+ for(int i = 0; i < n_nodes_by_edge; i++) {
+ for(int j = 0; j < n_nodes_by_edge; j++) {
+ double u = (double)i / (n_nodes_by_edge - 1);
+ double v = (double)j / (n_nodes_by_edge - 1);
+ Range<double> b1 = ge->parBounds(0);
+ Range<double> b2 = ge->parBounds(1);
+ double t1 = b1.low() + u * (b1.high() - b1.low());
+ double t2 = b2.low() + v * (b2.high() - b2.low());
+ GPoint gp = f->point(t1, t2);
+ getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
+ zeronodes[k][1], zeronodes[k][2], f);
+ _infos[k++] = AttractorInfo(*it,2,u,v);
+ }
+ }
}
- }
+ }
}
}
kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
@@ -1766,7 +1802,9 @@ BoundaryLayerField::BoundaryLayerField()
ratio = 1.1;
thickness = 1.e-2;
fan_angle = 30;
+ tgt_aniso_ratio = 1.e10;
iRecombine = 0;
+ iIntersect = 0;
options["NodesList"] = new FieldOptionList
(nodes_id, "Indices of nodes in the geometric model", &update_needed);
options["EdgesList"] = new FieldOptionList
@@ -1777,10 +1815,14 @@ BoundaryLayerField::BoundaryLayerField()
"layer is needed", &update_needed);
options["Quads"] = new FieldOptionInt
(iRecombine, "Generate recombined elements in the boundary layer");
+ options["IntersectMetrics"] = new FieldOptionInt
+ (iIntersect, "Intersect metrics of all faces");
options["hwall_n"] = new FieldOptionDouble
(hwall_n, "Mesh Size Normal to the The Wall");
options["fan_angle"] = new FieldOptionDouble
(fan_angle, "Threshold angle for creating a mesh fan in the boundary layer");
+ options["AnisoMax"] = new FieldOptionDouble
+ (tgt_aniso_ratio, "Threshold angle for creating a mesh fan in the boundary layer");
options["hwall_t"] = new FieldOptionDouble
(hwall_t, "Mesh Size Tangent to the Wall");
options["ratio"] = new FieldOptionDouble
@@ -1793,6 +1835,23 @@ BoundaryLayerField::BoundaryLayerField()
double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge)
{
+
+ if (update_needed){
+ for(std::list<int>::iterator it = nodes_id.begin();
+ it != nodes_id.end(); ++it) {
+ _att_fields.push_back(new AttractorField(0,*it,100000));
+ }
+ for(std::list<int>::iterator it = edges_id.begin();
+ it != edges_id.end(); ++it) {
+ _att_fields.push_back(new AttractorField(1,*it,300000));
+ }
+ for(std::list<int>::iterator it = faces_id.begin();
+ it != faces_id.end(); ++it) {
+ _att_fields.push_back(new AttractorField(2,*it,1200));
+ }
+ update_needed = false;
+ }
+
double dist = 1.e22;
AttractorField *cc;
for (std::list<AttractorField*>::iterator it = _att_fields.begin();
@@ -1831,122 +1890,59 @@ void BoundaryLayerField::operator() (AttractorField *cc, double dist,
const double ll2 = dist*(ratio-1) + hwall_t;
double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
- SVector3 t1,t2,t3;
- double L1,L2,L3;
+ lc_n = std::max(lc_n, CTX::instance()->mesh.lcMin);
+ lc_n = std::min(lc_n, CTX::instance()->mesh.lcMax);
+ lc_t = std::max(lc_t, CTX::instance()->mesh.lcMin);
+ lc_t = std::min(lc_t, CTX::instance()->mesh.lcMax);
std::pair<AttractorInfo,SPoint3> pp = cc->getAttractorInfo();
if (pp.first.dim ==0){
- L1 = lc_n;
- L2 = lc_n;
- L3 = lc_n;
GVertex *v = GModel::current()->getVertexByTag(pp.first.ent);
- t1 = SVector3(v->x() -x,v->y() -y,v->z() -z);
- t1.normalize();
- // printf("%p %g %g %g\n",v,t1.x(),t1.y(),t1.z());
- if (t1.norm() == 0)t1 = SVector3(1,0,0);
- if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
- t2 = SVector3(1,0,0);
- else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
- t2 = SVector3(0,1,0);
- else
- t2 = SVector3(0,0,1);
- t3 = crossprod(t1,t2);
- t3.normalize();
- t2 = crossprod(t3,t1);
- // printf("t1 %p %g %g %g\n",v,t1.x(),t1.y(),t1.z());
- // printf("t2 %p %g %g %g\n",v,t2.x(),t2.y(),t2.z());
- // printf("t3 %p %g %g %g\n",v,t3.x(),t3.y(),t3.z());
+ SVector3 t1 = SVector3(v->x() -x,v->y() -y,v->z() -z);
+ metr = buildMetricTangentToCurve(t1,lc_n,lc_n);
+ return;
}
else if (pp.first.dim ==1){
GEdge *e = GModel::current()->getEdgeByTag(pp.first.ent);
-
- // the tangent size at this point is the size of the
- // 1D mesh at this point !
- // hack : use curvature
- /*
- if(CTX::instance()->mesh.lcFromCurvature){
- double Crv = e->curvature(pp.first.u);
- double lc = Crv > 0 ? 2 * M_PI / Crv / CTX::instance()->mesh.minCircPoints : 1.e-22;
- const double ll2b = dist*(ratio-1) + lc;
- double lc_tb = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2b,hfar));
- lc_t = std::min(lc_t,lc_tb);
- }
- */
-
- if (dist < hwall_t){
- L1 = lc_t;
- L2 = lc_n;
- L3 = lc_n;
- t1 = e->firstDer(pp.first.u);
- t1.normalize();
- if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
- t2 = SVector3(1,0,0);
- else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
- t2 = SVector3(0,1,0);
- else
- t2 = SVector3(0,0,1);
- t3 = crossprod(t1,t2);
- t3.normalize();
- t2 = crossprod(t3,t1);
- // printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
+ if (dist < thickness){
+ SVector3 t1 = e->firstDer(pp.first.u);
+ double crv = e->curvature(pp.first.u);
+ const double b = lc_t;
+ const double h = lc_n;
+ double oneOverD2 = .5/(b*b) * (1. + sqrt (1. + ( 4.*crv*crv*b*b*b*b/ (h*h*CTX::instance()->mesh.smoothRatio*CTX::instance()->mesh.smoothRatio))));
+ metr = buildMetricTangentToCurve(t1,sqrt(1./oneOverD2),lc_n);
+ return;
}
else {
- L1 = lc_t;
- L2 = lc_n;
- L3 = lc_t;
GPoint p = e->point(pp.first.u);
- t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
- if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
- t1 = SVector3(1,0,0);
- else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
- t1 = SVector3(0,1,0);
- else
- t1 = SVector3(0,0,1);
- t2.normalize();
- t3 = crossprod(t1,t2);
- t3.normalize();
- t1 = crossprod(t3,t2);
+ SVector3 t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
+ metr = buildMetricTangentToCurve(t2,lc_t,lc_n);
+ return;
}
}
else {
GFace *gf = GModel::current()->getFaceByTag(pp.first.ent);
-
- if (dist < hwall_t){
- L1 = lc_n;
- L2 = lc_t;
- L3 = lc_t;
- t1 = gf->normal(SPoint2(pp.first.u,pp.first.v));
- t1.normalize();
- if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
- t2 = SVector3(1,0,0);
- else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
- t2 = SVector3(0,1,0);
- else
- t2 = SVector3(0,0,1);
- t3 = crossprod(t1,t2);
- t3.normalize();
- t2 = crossprod(t3,t1);
- // printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
+ if (dist < thickness){
+ double cmin, cmax;
+ SVector3 dirMax,dirMin;
+ cmax = gf->curvatures(SPoint2(pp.first.u,pp.first.v),&dirMax, &dirMin, &cmax,&cmin);
+ const double b = lc_t;
+ const double h = lc_n;
+ double oneOverD2_min = .5/(b*b) * (1. + sqrt (1. + ( 4.*cmin*cmin*b*b*b*b/ (h*h*CTX::instance()->mesh.smoothRatio*CTX::instance()->mesh.smoothRatio))));
+ double oneOverD2_max = .5/(b*b) * (1. + sqrt (1. + ( 4.*cmax*cmax*b*b*b*b/ (h*h*CTX::instance()->mesh.smoothRatio*CTX::instance()->mesh.smoothRatio))));
+ double dmin = sqrt(1./oneOverD2_min);
+ double dmax = sqrt(1./oneOverD2_max);
+ dmin = std::min(dmin,dmax*tgt_aniso_ratio);
+ metr = buildMetricTangentToSurface(dirMin,dirMax,dmin,dmax,lc_n);
+ return;
}
else {
- L1 = lc_t;
- L2 = lc_n;
- L3 = lc_t;
GPoint p = gf->point(SPoint2(pp.first.u,pp.first.v));
- t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
- if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
- t1 = SVector3(1,0,0);
- else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
- t1 = SVector3(0,1,0);
- else
- t1 = SVector3(0,0,1);
- t2.normalize();
- t3 = crossprod(t1,t2);
- t3.normalize();
- t1 = crossprod(t3,t2);
+ SVector3 t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
+ metr = buildMetricTangentToCurve(t2,lc_n,lc_t);
+ return;
}
}
- metr = SMetric3(1./(L1*L1), 1./(L2*L2), 1./(L3*L3), t1, t2, t3);
}
void BoundaryLayerField::operator() (double x, double y, double z,
@@ -1963,7 +1959,7 @@ void BoundaryLayerField::operator() (double x, double y, double z,
}
for(std::list<int>::iterator it = faces_id.begin();
it != faces_id.end(); ++it) {
- _att_fields.push_back(new AttractorField(2,*it,300));
+ _att_fields.push_back(new AttractorField(2,*it,1200));
}
update_needed = false;
}
@@ -1978,17 +1974,19 @@ void BoundaryLayerField::operator() (double x, double y, double z,
SPoint3 CLOSEST= (*it)->getAttractorInfo().second;
SMetric3 localMetric;
- (*this)(*it, cdist,x, y, z, localMetric, ge);
- hop.push_back(localMetric);
- //v = intersection(v,localMetric);
+ if (iIntersect){
+ (*this)(*it, cdist,x, y, z, localMetric, ge);
+ hop.push_back(localMetric);
+ }
if (cdist < current_distance){
+ if (!iIntersect)(*this)(*it, cdist,x, y, z, localMetric, ge);
current_distance = cdist;
current_closest = *it;
v = localMetric;
_closest_point = CLOSEST;
}
}
- // for (int i=0;i<hop.size();i++)v = intersection_conserveM1(v,hop[i]);
+ if (iIntersect)for (int i=0;i<hop.size();i++)v = intersection_conserve_mostaniso(v,hop[i]);
metr = v;
}
#endif
diff --git a/Mesh/Field.h b/Mesh/Field.h
index 8b562efc5c4aa7c34f60d54dd81274575d94f09d..aa7023987a85e3654714e1c2d2c290e366fe3db2 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -144,9 +144,9 @@ class BoundaryLayerField : public Field {
SMetric3 &metr, GEntity *ge);
public:
double hwall_n,hwall_t,ratio,hfar,thickness,fan_angle;
- double current_distance;
+ double current_distance, tgt_aniso_ratio;
SPoint3 _closest_point;
- int iRecombine;
+ int iRecombine, iIntersect;
AttractorField *current_closest;
virtual bool isotropic () const {return false;}
virtual const char *getName();
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index b8e0dfe5070265d345a6515ec38cca119256e626..82015f3824dd265cc7b02795f6524d1ced2429ee 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -121,7 +121,7 @@ static bool computeEquidistantParameters0(GEdge *ge, double u0, double uN, int N
return false;
}
// 1 = geodesics
-static int method_for_computing_intermediary_points = 0;
+static int method_for_computing_intermediary_points = 1;
static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
double *u, double underRelax){
if (method_for_computing_intermediary_points == 0) // use linear abscissa
@@ -1328,6 +1328,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// printJacobians(m, "smoothness.pos");
// m->writeMSH("SMOOTHED.msh");
// FIXME !!
+ checkHighOrderTriangles("Surface mesh", m, bad, worst);
if (!linear && CTX::instance()->mesh.smoothInternalEdges){
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
std::vector<MElement*> v;
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index b08328e5f91272fbdc324f1b843428b2dc9c0252..28fcbfe3cd966bed4bfe46b46c33307c2799766f 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -309,7 +309,7 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
double dx = dx0;
Msg::Debug(" dx0 = %12.5E", dx0);
int iter = 0;
- while(1){
+ while(0){
double dx2 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
Msg::Debug(" dx2 = %12.5E", dx2);
if (fabs(dx2 - dx) < 1.e-4 * dx0)break;
@@ -372,7 +372,7 @@ double highOrderTools::smooth_metric_(std::vector<MElement*> & v,
J23K33.gemm(J23, K33, 1, 0);
K22.gemm(J23K33, J32, 1, 0);
// K33.print("K33");
- // K22.print("K22");
+ // K22.print("K22");
J23K33.mult(D3, R2);
// R2.print("hopla");
for (int j = 0; j < n2; j++){
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 2863eb7cd6ecfb3739ff1789195b96d4d56639e9..064dd8ac2f83ebcc1fe92e205a1e7cc79b72c765 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -379,7 +379,7 @@ void meshGEdge::operator() (GEdge *ge)
SVector3 der = ge->firstDer(pt.t);
pt.xp = der.norm();
}
- a = smoothPrimitive(ge, CTX::instance()->mesh.smoothRatio, Points);
+ a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 8b12aef059497b165272c6ff96f145f9d7aecf31..12fc4e48e94f2ee4251beaa4ec3c7b06c4406b63 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1134,7 +1134,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf);
else {
- bowyerWatson(gf);
+ bowyerWatson(gf,15000);
meshGFaceBamg(gf);
}
laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 79647bf17a400e11403cc1488cce7e52ebaf2a64..1133b968c25678b41d14948366ee1bfac9495eca 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -802,6 +802,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
MTriangle *base = worst->tri();
Msg::Debug("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
+ // Msg::Info("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
+
AllTris.erase(it);
worst->forceRadius(-1);
AllTris.insert(worst);
@@ -834,7 +836,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
-void bowyerWatson(GFace *gf)
+void bowyerWatson(GFace *gf, int MAXPNT)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
std::vector<double> vSizes, vSizesBGM, Us, Vs;
@@ -871,7 +873,7 @@ void bowyerWatson(GFace *gf)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
vSizes.size(), worst->getRadius());
double center[2],metric[3],r2;
- if (worst->getRadius() < /*1.333333/(sqrt(3.0))*/0.5 * sqrt(2.0)) break;
+ if (worst->getRadius() < /*1.333333/(sqrt(3.0))*/0.5 * sqrt(2.0) || vSizes.size() > MAXPNT) break;
circUV(worst->tri(), Us, Vs, center, gf);
MTriangle *base = worst->tri();
double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index d5e1081a176d85b4d3911e0e81068d5a6fe21f38..8782668b315124cfee518ddc8e0827c97ac38b8a 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -138,7 +138,7 @@ void connectQuads(std::vector<MQua4*> &);
void connectTriangles(std::list<MTri3*> &);
void connectTriangles(std::vector<MTri3*> &);
void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
-void bowyerWatson(GFace *gf);
+void bowyerWatson(GFace *gf, int MAXPNT= 1000000000);
void bowyerWatsonFrontal(GFace *gf);
void bowyerWatsonFrontalLayers(GFace *gf, bool quad);
void buildBackGroundMesh (GFace *gf);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index f0a99fb19ece8a4b67d7ef4965415b154416cbb5..e5106e97b973d8cbcf38db26228f722d966e60d2 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1054,7 +1054,7 @@ struct quadBlob {
v2t_cont :: const_iterator it = adj.find(v);
int ss = temp.size();
std::vector<MElement*> elems = it->second;
-
+ /*
//EMI: try to orient BNodes the same as quad orientations
for (int j=0;j<elems.size();j++){
bool found =false;
@@ -1072,6 +1072,31 @@ struct quadBlob {
}
if (found) break;
}
+ */
+ for (int j=0;j<elems.size();j++){
+ bool found = false;
+ for (int i=0;i<4;i++){
+ MEdge e = it->second[j]->getEdge(i);
+ if (e.getVertex(0) == v &&
+ inBoundary(e.getVertex(1),bnodes) &&
+ !inBoundary(e.getVertex(1),temp)) {
+ v = e.getVertex(1);
+ temp.push_back(e.getVertex(1));
+ found = true;
+ break;
+ }
+ else if (e.getVertex(1) == v &&
+ inBoundary(e.getVertex(0),bnodes) &&
+ !inBoundary(e.getVertex(0),temp)) {
+ v = e.getVertex(0);
+ temp.push_back(e.getVertex(0));
+ found = true;
+ break;
+ }
+ }
+ if (found)break;
+ }
+
//JF this does not orient
// for (int j=0;j<elems.size();j++){
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 950c643777bd733a7c2a7c6046a3b92e1c118f64..325e598694fc59d0b3cf5b53b36670ef0dbe7611 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -526,11 +526,11 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D ||
CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD ||
CTX::instance()->mesh.algo2d == ALGO_2D_BAMG){
- sprintf(opts, "pY", (Msg::GetVerbosity() < 3) ? 'Q':
+ sprintf(opts, "-q1.5pY", (Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
}
else {
- sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q':
+ sprintf(opts, "-q3.5Ype%c", (Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
}
try{
@@ -594,11 +594,13 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
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)
+ else if(CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D){
refineMeshMMG(gr);
+ }
else
- if(!Filler::get_nbr_new_vertices() && !LpSmoother::get_nbr_interior_vertices())
+ if(!Filler::get_nbr_new_vertices() && !LpSmoother::get_nbr_interior_vertices()){
insertVerticesInRegion(gr);
+ }
#endif
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 154ca9e11675eddfcd00e3b5dbb3e281627ae7f5..8aef9689639f66f1fb8bcc9fa0b46ab4f618d567 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -199,7 +199,6 @@ bool insertVertex(MVertex *v,
while (it != shell.end()){
MTetrahedron *tr = new MTetrahedron(it->v[0], it->v[1], it->v[2], v);
- const double ONE_THIRD = 1./3.;
double lc = .25 * (vSizes[tr->getVertex(0)->getIndex()] +
vSizes[tr->getVertex(1)->getIndex()] +
vSizes[tr->getVertex(2)->getIndex()] +
@@ -259,6 +258,7 @@ bool insertVertex(MVertex *v,
return true;
}
else { // The cavity is NOT star shaped
+ // printf("hola %d %g %g %22.15E\n",onePointIsTooClose, oldVolume, newVolume, 100.*fabs(oldVolume - newVolume) / oldVolume);
for (unsigned int i = 0; i <shell.size(); i++) myFactory.Free(newTets[i]);
delete [] newTets;
ittet = cavity.begin();
@@ -746,6 +746,85 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
}
}
+
+static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
+ double circumcenter[3], double *xi, double *eta, double *zeta){
+ double xba, yba, zba, xca, yca, zca, xda, yda, zda;
+ double balength, calength, dalength;
+ double xcrosscd, ycrosscd, zcrosscd;
+ double xcrossdb, ycrossdb, zcrossdb;
+ double xcrossbc, ycrossbc, zcrossbc;
+ double denominator;
+ double xcirca, ycirca, zcirca;
+
+ /* Use coordinates relative to point `a' of the tetrahedron. */
+ xba = b[0] - a[0];
+ yba = b[1] - a[1];
+ zba = b[2] - a[2];
+ xca = c[0] - a[0];
+ yca = c[1] - a[1];
+ zca = c[2] - a[2];
+ xda = d[0] - a[0];
+ yda = d[1] - a[1];
+ zda = d[2] - a[2];
+ /* Squares of lengths of the edges incident to `a'. */
+ balength = xba * xba + yba * yba + zba * zba;
+ calength = xca * xca + yca * yca + zca * zca;
+ dalength = xda * xda + yda * yda + zda * zda;
+ /* Cross products of these edges. */
+ xcrosscd = yca * zda - yda * zca;
+ ycrosscd = zca * xda - zda * xca;
+ zcrosscd = xca * yda - xda * yca;
+ xcrossdb = yda * zba - yba * zda;
+ ycrossdb = zda * xba - zba * xda;
+ zcrossdb = xda * yba - xba * yda;
+ xcrossbc = yba * zca - yca * zba;
+ ycrossbc = zba * xca - zca * xba;
+ zcrossbc = xba * yca - xca * yba;
+
+ /* Calculate the denominator of the formulae. */
+ /* Use orient3d() from http://www.cs.cmu.edu/~quake/robust.html */
+ /* to ensure a correctly signed (and reasonably accurate) result, */
+ /* avoiding any possibility of division by zero. */
+ const double xxx = robustPredicates::orient3d(b, c, d, a);
+ denominator = 0.5 / xxx;
+
+ /* Calculate offset (from `a') of circumcenter. */
+ xcirca = (balength * xcrosscd + calength * xcrossdb + dalength * xcrossbc) *
+ denominator;
+ ycirca = (balength * ycrosscd + calength * ycrossdb + dalength * ycrossbc) *
+ denominator;
+ zcirca = (balength * zcrosscd + calength * zcrossdb + dalength * zcrossbc) *
+ denominator;
+ circumcenter[0] = xcirca + a[0];
+ circumcenter[1] = ycirca + a[1];
+ circumcenter[2] = zcirca + a[2];
+
+ /*
+ printf(" %g %g %g %g\n",
+ sqrt((a[0]-xcirca)*(a[0]-xcirca)+(a[1]-ycirca)*(a[1]-ycirca)+(a[2]-zcirca)*(a[2]-zcirca)),
+ sqrt((b[0]-xcirca)*(b[0]-xcirca)+(b[1]-ycirca)*(b[1]-ycirca)+(b[2]-zcirca)*(b[2]-zcirca)),
+ sqrt((c[0]-xcirca)*(c[0]-xcirca)+(c[1]-ycirca)*(c[1]-ycirca)+(c[2]-zcirca)*(c[2]-zcirca)),
+ sqrt((d[0]-xcirca)*(d[0]-xcirca)+(d[1]-ycirca)*(d[1]-ycirca)+(d[2]-zcirca)*(d[2]-zcirca)) );
+ */
+
+ if (xi != (double *) NULL) {
+ /* To interpolate a linear function at the circumcenter, define a */
+ /* coordinate system with a xi-axis directed from `a' to `b', */
+ /* an eta-axis directed from `a' to `c', and a zeta-axis directed */
+ /* from `a' to `d'. The values for xi, eta, and zeta are computed */
+ /* by Cramer's Rule for solving systems of linear equations. */
+ *xi = (xcirca * xcrosscd + ycirca * ycrosscd + zcirca * zcrosscd) *
+ (2.0 * denominator);
+ *eta = (xcirca * xcrossdb + ycirca * ycrossdb + zcirca * zcrossdb) *
+ (2.0 * denominator);
+ *zeta = (xcirca * xcrossbc + ycirca * ycrossbc + zcirca * zcrossbc) *
+ (2.0 * denominator);
+ }
+ return xxx;
+}
+
+
void insertVerticesInRegion (GRegion *gr)
{
//printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
@@ -836,10 +915,37 @@ void insertVerticesInRegion (GRegion *gr)
vSizes.size(), worst->getRadius());
if(worst->getRadius() < 1) break;
double center[3];
- worst->circumcenter(center);
double uvw[3];
- worst->tet()->xyz2uvw(center, uvw);
- if(worst->tet()->isInside(uvw[0], uvw[1], uvw[2])){
+ MTetrahedron *base = worst->tet();
+ double pa[3] = {base->getVertex(0)->x(),
+ base->getVertex(0)->y(),
+ base->getVertex(0)->z()};
+ double pb[3] = {base->getVertex(1)->x(),
+ base->getVertex(1)->y(),
+ base->getVertex(1)->z()};
+ double pc[3] = {base->getVertex(2)->x(),
+ base->getVertex(2)->y(),
+ base->getVertex(2)->z()};
+ double pd[3] = {base->getVertex(3)->x(),
+ base->getVertex(3)->y(),
+ base->getVertex(3)->z()};
+
+ // double UU,VV,WW;
+ tetcircumcenter(pa,pb,pc,pd, center,&uvw[0],&uvw[1],&uvw[2] );
+ // worst->tet()->xyz2uvw(center, uvw);
+ // printf("%12.5E %12.5E %12.5E -- %12.5E %12.5E %12.5E \n",
+ // UU,VV,WW,uvw[0],uvw[1],uvw[2]);
+ /*
+ if (!worst->tet()->isInside(uvw[0], uvw[1], uvw[2]) &&
+ worst->getRadius() > 20){
+ uvw[0] = uvw[1] = uvw[2] = 0.25;
+ center[0] = 0.25*(pa[0]+pb[0]+pc[0]+pd[0]);
+ center[1] = 0.25*(pa[1]+pb[1]+pc[1]+pd[1]);
+ center[2] = 0.25*(pa[2]+pb[2]+pc[2]+pd[2]);
+ }
+ */
+
+ if(worst->tet()->isInside(uvw[0], uvw[1], uvw[2])){
MVertex *v = new MVertex(center[0], center[1], center[2], worst->onWhat());
v->setIndex(NUM++);
double lc1 =
@@ -860,7 +966,8 @@ void insertVerticesInRegion (GRegion *gr)
v->onWhat()->mesh_vertices.push_back(v);
}
else{
- myFactory.changeTetRadius(allTets.begin(), 0.);
+ // 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]);
+ myFactory.changeTetRadius(allTets.begin(), 0.0);
}
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 331de8b3f209bab3527076fc3fb2ac586280d572..7b7c5d5fc0ef9a16389becaefc9e912e137e0bad 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -14,6 +14,7 @@
#include "Numeric.h"
#include "BackgroundMesh.h"
#include "qualityMeasures.h"
+#include "robustPredicates.h"
//#define _GMSH_PRE_ALLOCATE_STRATEGY_ 1
class GRegion;
@@ -148,7 +149,22 @@ class MTet4
{
return inCircumSphere(v->x(), v->y(), v->z());
}
- inline double getVolume() const { return base->getVolume(); }
+ inline double getVolume() const {
+
+ double pa[3] = {base->getVertex(0)->x(),
+ base->getVertex(0)->y(),
+ base->getVertex(0)->z()};
+ double pb[3] = {base->getVertex(1)->x(),
+ base->getVertex(1)->y(),
+ base->getVertex(1)->z()};
+ double pc[3] = {base->getVertex(2)->x(),
+ base->getVertex(2)->y(),
+ base->getVertex(2)->z()};
+ double pd[3] = {base->getVertex(3)->x(),
+ base->getVertex(3)->y(),
+ base->getVertex(3)->z()};
+ return fabs(robustPredicates::orient3d(pa, pb, pc, pd))/6.0;
+ }
inline void setDeleted(bool d)
{
deleted = d;
diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp
index 1a2cb0c8e5487efe95de076bdd013d5b001a9375..306d4fc8372d0c73f1d8be55bd7379f4736d5bcb 100644
--- a/Mesh/meshGRegionMMG3D.cpp
+++ b/Mesh/meshGRegionMMG3D.cpp
@@ -137,13 +137,13 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
if (itv != LCS.end()){
mmg2gmsh[(*it)->getNum()] = *it;
- if (CTX::instance()->mesh.lcExtendFromBoundary){
+ //if (CTX::instance()->mesh.lcExtendFromBoundary){
double LL = itv->second.first/itv->second.second;
SMetric3 l4(1./(LL*LL));
- SMetric3 MM = intersection (l4, m);
+ SMetric3 MM = intersection_conserve_mostaniso (l4, m);
m = MM;
//lc = std::min(LL,lc);
- }
+ // }
}
sol->met[count++] = m(0,0);
@@ -185,12 +185,12 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
}
-static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol)
+static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*> &mmg2gmsh)
{
std::list<GFace*> f = gr->faces();
- /*
+
std::map<MVertex*,std::pair<double,int> > LCS;
- if (CTX::instance()->mesh.lcExtendFromBoundary){
+ // if (CTX::instance()->mesh.lcExtendFromBoundary){
for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
for (int i=0;i<(*it)->triangles.size();i++){
MTriangle *t = (*it)->triangles[i];
@@ -208,26 +208,34 @@ static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol)
}
}
}
- }
- */
-
- int count = 1;
- for (int k=1 ; k<=mmg->np; k++){
- MMG_pPoint ppt = &mmg->point[k];
- if (ppt->tag & M_UNUSED) continue;
+ // }
+
- SMetric3 m = BGM_MeshMetric(gr, 0,0,ppt->c[0],ppt->c[1],ppt->c[2]);
- /*
- if (CTX::instance()->mesh.lcExtendFromBoundary){
- std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
+ int count = 1;
+ for (int k=1 ; k<=mmg->np; k++){
+ MMG_pPoint ppt = &mmg->point[k];
+ if (ppt->tag & M_UNUSED) continue;
+
+ SMetric3 m = BGM_MeshMetric(gr, 0,0,ppt->c[0],ppt->c[1],ppt->c[2]);
+
+ std::map<int,MVertex*>::iterator it = mmg2gmsh.find(k);
+
+ if (it != mmg2gmsh.end() && CTX::instance()->mesh.lcExtendFromBoundary){
+ std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(it->second);
if (itv != LCS.end()){
double LL = itv->second.first/itv->second.second;
SMetric3 l4(1./(LL*LL));
- SMetric3 MM = intersection (l4, m);
+ // printf("adding a size %g\n",LL);
+ SMetric3 MM = intersection_conserve_mostaniso (l4, m);
m = MM;
}
}
- */
+ if (m.determinant() < 1.e-30){
+ m(0,0) += 1.e-12;
+ m(1,1) += 1.e-12;
+ m(2,2) += 1.e-12;
+ }
+
sol->met[count++] = m(0,0);
sol->met[count++] = m(1,0);
sol->met[count++] = m(2,0);
@@ -273,10 +281,10 @@ void refineMeshMMG(GRegion *gr)
Msg::Info("MMG3D succeeded %d vertices %d tetrahedra",
mmg->np, mmg->ne);
// Here we should interact with BGM
- updateSizes(gr,mmg, sol);
+ updateSizes(gr,mmg, sol,mmg2gmsh);
int nTnow = mmg->ne;
- if (fabs((double)(nTnow - nT)) < 0.05 * nT) break;
+ if (fabs((double)(nTnow - nT)) < 0.03 * nT) break;
}
//char test[] = "test.mesh";
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 5deca329eaf15e0c9b73a0ce19de5ef9a1b3e65c..e2b599d3aa16c4dd7cd457a61a7d89922b4b1e59 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -648,7 +648,7 @@ void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
fullVector<double> &x, void *data, double relax, double tolx)
{
- const int MAXIT = 50;
+ const int MAXIT = 10;
const double EPS = 1.e-4;
const int N = x.size();