diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 7096523b7166d8105f68b669b2b7ba0db546cbce..574e1f971f3eccfab6be7eb3054f6f38bc3bc873 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -22,8 +22,6 @@ void STensor3::print (const char *s) const
(*this)(2,0), (*this)(2,1), (*this)(2,2));
}
-//// ORIGINAL CODE FROM GMSH
-/*
SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
{
SMetric3 im1 = m1.invert();
@@ -37,48 +35,16 @@ SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
double l0 = std::max(dot(v0,m1,v0),dot(v0,m2,v0));
double l1 = std::max(dot(v1,m1,v1),dot(v1,m2,v1));
double l2 = std::max(dot(v2,m1,v2),dot(v2,m2,v2));
- SMetric3 iv(l0,l1,l2,v0,v1,v2);
-
- return iv;
-}
-*/
-
-//// MODIFIED INTERSECTION
-SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
-{
- SMetric3 im1 = m1.invert();
- fullMatrix<double> V(3,3);
- fullVector<double> S(3);
- im1 *= m2;
- im1.eig(V,S,true);
- SVector3 v0(V(0,0),V(1,0),V(2,0));
- SVector3 v1(V(0,1),V(1,1),V(2,1));
- SVector3 v2(V(0,2),V(1,2),V(2,2));
- double l0 = std::max(dot(v0,m1,v0),dot(v0,m2,v0));
- double l1 = std::max(dot(v1,m1,v1),dot(v1,m2,v1));
- double l2 = std::max(dot(v2,m1,v2),dot(v2,m2,v2));
-
-
- // Correction as explained in the PhD thesis of Frederic Alauzet p.16
- double max_eig = std::max(l0, std::max(l1,l2));
- double min_eig = std::min(l0, std::min(l1,l2));
- double range_eig = (max_eig - min_eig);
- if( range_eig < 1.0e-1) // Change this value if you think it's too big ...
- {
- if(max_eig <= 1.0)
- {
- return m1;
- }
-
- else
- {
- return m2;
- }
- }
+ // Correction from the PhD thesis of Frederic Alauzet p.16
+ // CAUTION: error in the thesis, cases alpha<=1 and alpha>1 are switched
+ static const double eps = 1.e-2; // Tolerance to detect triple eigenvalue (i.e. proportional metrics)
+ const double max_eig = std::max(l0, std::max(l1,l2));
+ const double min_eig = std::min(l0, std::min(l1,l2));
+ const double range_eig = fabs((max_eig-min_eig)/max_eig);
+ if (range_eig < eps) return (max_eig <= 1.) ? m2 : m1;
SMetric3 iv(l0,l1,l2,v0,v1,v2);
-
return iv;
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index e6ed5405035d89d41d7bfae931061cace31c9b0c..a4750c55f927294a1ab3377847b1b2bdb862c445 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -402,73 +402,46 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
// anisotropic version of the background field - for now, only works
// with bamg in 2D, work in progress
-
SMetric3 BGM_MeshMetric(GEntity *ge,
double U, double V,
double X, double Y, double Z)
{
- // 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;
- 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));
- 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));
- FieldManager *fields = ge->model()->getFields();
- if(fields->getBackgroundField() > 0){
- Field *f = fields->get(fields->getBackgroundField());
- if(f){
- if (!f->isotropic()){
- (*f)(X, Y, Z, l4,ge);
- }
- else{
- double L = (*f)(X, Y, Z, ge);
- l4 = SMetric3(1/(L*L));
- }
- }
- }
- // take the minimum, then constrain by lcMin and lcMax
- double lc = std::min(l1, l2);
+ // Metrics based on element size
+ // Element size = min. between default lc and lc from point (if applicable), constrained by lcMin and lcMax
+ double lc = CTX::instance()->lc;
+ if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2) lc = std::min(lc, LC_MVertex_PNTS(ge, U, V));
lc = std::max(lc, CTX::instance()->mesh.lcMin);
lc = std::min(lc, CTX::instance()->mesh.lcMax);
-
-
if(lc <= 0.){
Msg::Error("Wrong mesh element size lc = %g (lcmin = %g, lcmax = %g)",
lc, CTX::instance()->mesh.lcMin, CTX::instance()->mesh.lcMax);
- lc = l1;
+ lc = CTX::instance()->lc;
}
+ SMetric3 m0(1./(lc*lc));
- SMetric3 LC(1./(lc*lc));
- //SMetric3 m = intersection_conserveM1(intersection_conserveM1 (l4, LC),l3);
- SMetric3 m = intersection(intersection (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));
+ // Intersect with metrics from fields if applicable
+ FieldManager *fields = ge->model()->getFields();
+ SMetric3 m1 = m0;
+ if(fields->getBackgroundField() > 0){
+ Field *f = fields->get(fields->getBackgroundField());
+ if(f) {
+ SMetric3 l4;
+ if (!f->isotropic()) (*f)(X, Y, Z, l4, ge);
+ else {
+ const double L = (*f)(X, Y, Z, ge);
+ l4 = SMetric3(1/(L*L));
+ }
+ m1 = intersection(l4, m0);
}
}
+ // Intersect with metrics from curvature if applicable
+ SMetric3 m = (CTX::instance()->mesh.lcFromCurvature && ge->dim() < 3) ?
+ intersection(m1, LC_MVertex_CURV_ANISO(ge, U, V)) : m1;
return m;
- // return lc * CTX::instance()->mesh.lcFactor;
+
}
bool Extend1dMeshIn2dSurfaces()
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 58401607bcac1faf953b5ed8bd4cff15e4d342fd..d185e6ce9f4ed5b16f9a17f909ec31d3bcff3c22 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -81,7 +81,7 @@ meshMetric::meshMetric(std::vector<MElement*> elements)
void meshMetric::addMetric(int technique, simpleFunction<double> *fct,
std::vector<double> parameters)
{
- needMetricUpdate=true;
+ needMetricUpdate = true;
_fct = fct;
hmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
hmax = parameters.size() >=3 ? parameters[2] : CTX::instance()->mesh.lcMax;
@@ -96,7 +96,7 @@ void meshMetric::addMetric(int technique, simpleFunction<double> *fct,
computeMetric();
}
-void meshMetric::intersectMetrics()
+void meshMetric::updateMetrics()
{
if (!setOfMetrics.size()){
std::cout << " meshMetric::intersectMetrics: Can't intersect metrics, no metric registered ! " << std::endl;
@@ -109,10 +109,11 @@ void meshMetric::intersectMetrics()
_nodalMetrics[ver] = setOfMetrics[0][ver];
// _detMetric[ver] = setOfDetMetric[0][ver];
_nodalSizes[ver] = setOfSizes[0][ver];
- for (unsigned int i=1;i<setOfMetrics.size();i++){
- _nodalMetrics[ver] = intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]);
- _nodalSizes[ver] = std::min(_nodalSizes[ver],setOfSizes[i][ver]);
- }
+ if (setOfMetrics.size() > 1)
+ for (unsigned int i=1;i<setOfMetrics.size();i++){
+ _nodalMetrics[ver] = intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]);
+ _nodalSizes[ver] = std::min(_nodalSizes[ver],setOfSizes[i][ver]);
+ }
// _detMetric[ver] = sqrt(_nodalMetrics[ver].determinant());
}
needMetricUpdate=false;
@@ -121,7 +122,7 @@ void meshMetric::intersectMetrics()
void meshMetric::exportInfo(const char * fileendname)
{
- if (needMetricUpdate) intersectMetrics();
+ if (needMetricUpdate) updateMetrics();
std::stringstream sg,sm,sl,sh;
sg << "meshmetric_gradients_" << fileendname;
sm << "meshmetric_metric_" << fileendname;
@@ -762,7 +763,7 @@ void meshMetric::computeMetric()
double meshMetric::operator() (double x, double y, double z, GEntity *ge)
{
- if (needMetricUpdate) intersectMetrics();
+ if (needMetricUpdate) updateMetrics();
if (!setOfMetrics.size()){
std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
throw;
@@ -799,7 +800,7 @@ double meshMetric::operator() (double x, double y, double z, GEntity *ge)
void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge)
{
if (needMetricUpdate) {
- intersectMetrics();
+ updateMetrics();
}
if (!setOfMetrics.size()){
std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 5bf50e131249f625539818eba489695230940fc7..0fba36cd2f97bb5c687c4ca627f0bc52fcd0d225 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -26,7 +26,7 @@ class meshMetric: public Field {
ISOTROPIC_LINEARINTERP_H=6, SCALED_HESSIAN=7} MetricComputationTechnique;
private:
// intersect all metrics added in "setOfMetrics", preserve eigendirections of the "most anisotropic" metric
- void intersectMetrics();
+ void updateMetrics();
int _dim;
double _epsilon, _E, _E_moins, _Np;
bool needMetricUpdate;
@@ -80,7 +80,7 @@ class meshMetric: public Field {
void addMetric(int technique, simpleFunction<double> *fct, std::vector<double> parameters);
inline SMetric3 metricAtVertex (MVertex* v) {
- if (needMetricUpdate) intersectMetrics();
+ if (needMetricUpdate) updateMetrics();
return _nodalMetrics[v];
}
// this function scales the mesh metric in order