diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 49b262a1e14848c8ce1f36125822a1a2fde6cffe..f77eeb18a37fa925e27d7b28ea047a7b54a00ae2 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -1439,6 +1439,64 @@ void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs)
}
+//========================================================================================================
+
+void Curvature::edgeNodalValuesAndDirections(MLine* edge, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs)
+{
+ std::vector<MVertex*> LineVertices;
+ //LineVertices.resize(2);
+ edge->getEdgeVertices(0,LineVertices);//triangle->getVertex(0);
+
+ int V0 = 0;
+ int V1 = 0;
+
+ //LineVertices[0] is a pointer the first vertex of the edge
+ //LineVertices[1] is a pointer the second vertex of the edge
+
+
+ std::map<int,int>::iterator vertexIterator;
+ vertexIterator = _VertexToInt.find( LineVertices[0]->getNum() );
+ if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << LineVertices[0]->getNum() << " in _VertextToInt !" << std::endl;
+
+ vertexIterator = _VertexToInt.find( LineVertices[1]->getNum() );
+ if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << LineVertices[1]->getNum() << " in _VertextToInt !" << std::endl;
+
+
+ if (isAbs){
+ dMax[0] = _pdir1[V0];
+ dMax[1] = _pdir1[V1];
+
+ dMin[0] = _pdir2[V0];
+ dMin[1] = _pdir2[V1];
+
+ cMax[0] = std::abs(_curv1[V0]);
+ cMax[1] = std::abs(_curv1[V1]);
+
+ cMin[0] = std::abs(_curv2[V0]);
+ cMin[1] = std::abs(_curv2[V1]);
+
+ }
+ else{
+
+ dMax[0] = _pdir1[V0];
+ dMax[1] = _pdir1[V1];
+
+ dMin[0] = _pdir2[V0];
+ dMin[1] = _pdir2[V1];
+
+ cMax[0] = _curv1[V0];
+ cMax[1] = _curv1[V1];
+
+ cMin[0] = _curv2[V0];
+ cMin[1] = _curv2[V1];
+ }
+}
+
+
//========================================================================================================
double Curvature::getAtVertex(const MVertex *v) const {
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index a4674916ac030a458932fe495557f268103982de..cda44bdcc993693e3feb7aea277945c347012bc7 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -209,10 +209,13 @@ public:
void computeCurvature_RBF();
void triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs=0);
+
void triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs=0);
void edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs=0);
+ void edgeNodalValuesAndDirections(MLine* edge, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs=0);
+
void writeToMshFile( const std::string & filename);
void writeToPosFile( const std::string & filename);
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index a5e44e97bd978d65c1d5b0a8b64c0731e7cbced1..6b79a31df1d92c4f1004cd07f9546ca77455b124 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -10,7 +10,9 @@
#include <stdlib.h>
#include "GmshConfig.h"
#include "GEdgeCompound.h"
+#include "discreteEdge.h"
#include "Numeric.h"
+#include "Curvature.h"
GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound,
std::vector<int> &orientation)
@@ -241,6 +243,51 @@ double GEdgeCompound::curvature(double par) const
return _compound[iEdge]->curvature(tLoc);
}
+double GEdgeCompound::curvatures(const double par, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const
+{
+ double tLoc = -1.0;
+ int iEdge;
+ getLocalParameter(par,iEdge,tLoc);
+
+ if( _compound[iEdge]->geomType() == GEntity::DiscreteCurve)
+ {
+ Curvature& curvature = Curvature::getInstance();
+ if( !Curvature::valueAlreadyComputed() )
+ {
+ Msg::Info("Need to compute discrete curvature for anisotropic remesh (in GFace)");
+ Curvature::typeOfCurvature type = Curvature::RUSIN; //RBF
+ curvature.computeCurvature(model(), type);
+ }
+ double tLocMLine;
+ int iMLine;
+ discreteEdge* ptrEdge = dynamic_cast<discreteEdge*>(_compound[iEdge]);
+ ptrEdge->getLocalParameter(tLoc, iMLine, tLocMLine);
+
+ //vector of MLines
+ MLine* mline = ptrEdge->lines[iMLine];
+
+ double cMin[2];
+ double cMax[2];
+ SVector3 dMin[2];
+ SVector3 dMax[2];
+ curvature.edgeNodalValuesAndDirections(mline, dMax, dMin, cMax, cMin, 0);
+
+ *dirMax = (1-tLocMLine)*dMax[0] + tLocMLine*dMax[1];
+ *dirMin = (1-tLocMLine)*dMin[0] + tLocMLine*dMin[1];
+ *curvMax = (1-tLocMLine)*cMax[0] + tLocMLine*cMax[1];
+ *curvMin = (1-tLocMLine)*cMin[0] + tLocMLine*cMin[1];
+
+ return *curvMax;
+
+ }
+ else
+ {
+ Msg::Error("Case of CAD Geometry, don't know what to do here...");
+ }
+ return 0.0;
+}
+
GPoint GEdgeCompound::point(double par) const
{
double tLoc;
diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index d30fb8fd83d986bc1a2a5512e20672870a458eee..9562588528b6d9473604a20dfbc28699cb986d7d 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -32,6 +32,8 @@ class GEdgeCompound : public GEdge {
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return 0; }
virtual double curvature(double t) const;
+ virtual double curvatures(const double par, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const;
virtual int minimumMeshSegments() const;
virtual int minimumDrawSegments() const;
std::vector<GEdge*> getCompounds() const { return _compound; }
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 87039f8aa84634b11d906e66017cdfef27cfddca..fd8047b94e792b31a6d007d639b5b3b40d32cf77 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -20,6 +20,7 @@
#include "Geo.h"
#include "OS.h"
#include "Curvature.h"
+#include "GEdgeCompound.h"
discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
: GEdge(model, num, _v0, _v1)
@@ -513,6 +514,18 @@ double discreteEdge::curvature(double par) const
}
+double discreteEdge::curvatures(const double par, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const
+{
+ if (getCompound()){
+ return getCompound()->curvatures(par, dirMax, dirMin, curvMax, curvMin);
+ }
+ else{
+ Msg::Error("Cannot evaluate curvatures and curvature directions on discrete edge");
+ return false;
+ }
+}
+
Range<double> discreteEdge::parBounds(int i) const
{
return Range<double>(0, lines.size());
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index a136f8c8cd3a2ace5ad4d7307f3386b23ff149f3..67fc1f7573700bec095a68843d9bb7c1a7e1425b 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -23,6 +23,7 @@ class discreteEdge : public GEdge {
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
virtual double curvature(double par) const;
+ virtual double curvatures(const double par, SVector3 *dirMax, SVector3 *dirMin, double *curvMax, double *curvMin) const;
virtual bool haveParametrization(){ return !_pars.empty(); }
virtual Range<double> parBounds(int) const;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index a846a8f1a03bf63edba7794e29f4eef6341b8637..02e05831927324c14e18439a994e1a1e574f9825 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -174,11 +174,35 @@ static SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, dou
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;
- while(it != faces.end()){
+ const GEdgeCompound* ptrCompoundEdge = dynamic_cast<const GEdgeCompound*>(ge);
+ if (ptrCompoundEdge)
+ {
+ double cmax, cmin;
+ SVector3 dirMax,dirMin;
+ cmax = ptrCompoundEdge->curvatures(u,&dirMax, &dirMin, &cmax,&cmin);
+ if (cmin == 0)cmin =1.e-5;
+ if (cmax == 0)cmax =1.e-5;
+ 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);
+
+ 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);
+
+ SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),1.e-5,dirMin, dirMax, Z );
+ return curvMetric;
+ }
+
+ else
+ {
+ SMetric3 mesh_size(1.e-12);
+ std::list<GFace *> faces = ge->faces();
+ std::list<GFace *>::iterator it = faces.begin();
+ int count = 0;
+ while(it != faces.end())
+ {
if (((*it)->geomType() != GEntity::CompoundSurface) &&
((*it)->geomType() != GEntity::DiscreteSurface)){
SPoint2 par = ge->reparamOnFace((*it), u, 1);
@@ -193,42 +217,8 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
double lambda = ((2 * M_PI) /( fabs(Crv) * CTX::instance()->mesh.minCircPoints ) );
SMetric3 curvMetric (1./(lambda*lambda));
- return intersection(mesh_size,curvMetric);
-}
-
-static SMetric3 metric_based_on_surface_curvature_ES(const GEdge *ge, double u)
-{
- // If the curve is a compound curve, find the corresponding discrete curve. Then loop over all
- // neighboring (discrete) surfaces and compute the metric as the intersection of the surface metrics
-
- SMetric3 mesh_size(1.e-12);
- int iDiscEdge;
- double tLoc;
- if (ge->geomType() == GEntity::CompoundCurve)
- {
- const GEdgeCompound* ptrEdgeCompound = dynamic_cast<const GEdgeCompound*>(ge);
- ptrEdgeCompound->getLocalParameter(u, iDiscEdge, tLoc);
- std::vector <GEdge*> DiscreteEdge = ptrEdgeCompound->getCompounds(); //vector of pointers to GEdge
- GEdge* ptrDiscEdge = DiscreteEdge[iDiscEdge]; //Single pointer to GEdge
-
- std::list<GFace *> faces = ptrDiscEdge->faces();
- std::list<GFace *>::iterator it = faces.begin();
- int count = 0;
- while(it != faces.end())
- {
-
-// SMetric3 m = metric_based_on_surface_curvature(*it, u, v); //This doesn't work!!!!
- 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;
- }
- }
-
- return mesh_size;
+ return intersection_conserve_mostaniso(mesh_size,curvMetric);
+ }
}
static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
@@ -239,77 +229,27 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
ite != l_edges.end(); ++ite){
GEdge *_myGEdge = *ite;
Range<double> bounds = _myGEdge->parBounds(0);
- if (gv == _myGEdge->getBeginVertex())
- mesh_size = intersection
- (mesh_size,
- metric_based_on_surface_curvature(_myGEdge, bounds.low()));
- else
- mesh_size = intersection
- (mesh_size,
- metric_based_on_surface_curvature(_myGEdge, bounds.high()));
- }
- return mesh_size;
-}
-static SMetric3 metric_based_on_surface_curvature_ES(const GVertex *gv)
-{
- SMetric3 mesh_size(1.e-5);
- std::list <GFace*>::const_iterator itf;
-
- int iDiscEdge;
- double tLoc;
-
- std::list<GEdge*> l_edges = gv->edges();
- std::list<GFace*> l_faces;
- for (std::list<GEdge*>::const_iterator ite = l_edges.begin(); ite != l_edges.end(); ++ite)
- {
- if ((*ite)->geomType() == GEntity::CompoundCurve)
- {
- const GEdgeCompound* ptrEdgeCompound = dynamic_cast<const GEdgeCompound*>(*ite);
- Range<double> bounds = (*ite)->parBounds(0);
-
- if (gv == (*ite)->getBeginVertex())
- {
- ptrEdgeCompound->getLocalParameter(bounds.low(), iDiscEdge, tLoc);
- }
- else
- {
- ptrEdgeCompound->getLocalParameter(bounds.high(), iDiscEdge, tLoc);
- }
-
- std::vector <GEdge*> DiscreteEdge = ptrEdgeCompound->getCompounds(); //vector of pointers to GEdge
- GEdge* ptrDiscEdge = DiscreteEdge[iDiscEdge]; //Single pointer to GEdge
-
-
- std::list<GFace *> neighFaces = ptrDiscEdge->faces();
- for (itf = neighFaces.begin(); itf != neighFaces.end(); ++itf)
- {
- l_faces.push_back(*itf);
- }
- }
- }
-
- l_faces.sort();
- l_faces.unique();
-
- itf = l_faces.begin();
- int count = 0;
- while(itf != l_faces.end())
- {
- SPoint2 par = gv->reparamOnFace((*itf), 1);
- SMetric3 m = metric_based_on_surface_curvature (*itf, par.x(), par.y());
- if (!count) mesh_size = m;
- else mesh_size = intersection(mesh_size, m);
- count++;
+ // ES: Added extra if condition to use the code below only with compund curves
+ // This is because we want to call the function
+ // metric_based_on_surface_curvature(const GEdge *ge, double u) for the case when
+ // ge is a compound edge
+ if (_myGEdge->geomType() == GEntity::CompoundCurve)
+ {
+ if (gv == _myGEdge->getBeginVertex())
+ mesh_size = intersection
+ (mesh_size,
+ metric_based_on_surface_curvature(_myGEdge, bounds.low()));
+ else
+ mesh_size = intersection
+ (mesh_size,
+ metric_based_on_surface_curvature(_myGEdge, bounds.high()));
+ }
- ++itf;
}
-
return mesh_size;
}
-
-
// the mesh vertex is classified on a model vertex. we compute the
// maximum of the curvature of model faces surrounding this point if
// it is classified on a model edge, we do the same for all model