diff --git a/CMakeLists.txt b/CMakeLists.txt
index 24f6b03a09eff41cb2d9fc4bf431d0698bf65f1d..93bdf2cb985a6d7fc94e17bec6002661f95befa6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -102,7 +102,7 @@ set(GMSH_API
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h
Mesh/meshGFaceElliptic.h Mesh/meshPartition.h Mesh/meshGFaceDelaunayInsertion.h
Mesh/simple3D.h Mesh/meshPartitionOptions.h Mesh/directions3D.h Mesh/yamakawa.h
- Mesh/Voronoi3D.h Mesh/Levy3D.h Mesh/periodical.h
+ Mesh/Voronoi3D.h Mesh/Levy3D.h Mesh/periodical.h Mesh/meshMetric.h
Numeric/mathEvaluator.h
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Solver/crossConfTerm.h Solver/orthogonalTerm.h
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index e52fd40b5c1b7665d7c14cc98e43016e2c96bd54..59de16224e6306d9e8f5d3f0cd1bb2f7a7b7ce99 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -105,13 +105,13 @@ void meshMetric::intersectMetrics()
for (;it != _adj.end();it++) {
MVertex *ver = it->first;
_nodalMetrics[ver] = setOfMetrics[0][ver];
- _detMetric[ver] = setOfDetMetric[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]);
}
- _detMetric[ver] = sqrt(_nodalMetrics[ver].determinant());
+ // _detMetric[ver] = sqrt(_nodalMetrics[ver].determinant());
}
needMetricUpdate=false;
@@ -356,16 +356,20 @@ void meshMetric::computeMetricLevelSet()
_hessian[ver] = hessian;
setOfSizes[metricNumber].insert(std::make_pair(ver, std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3))));
setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
- setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
+ // setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
}
}
-void meshMetric::computeMetricHessian()
+void meshMetric::computeMetricHessian( )
{
int metricNumber = setOfMetrics.size();
+ double _epsilonP = 1.;
+ double hminP = 1.e-12;
+ double hmaxP = 1.e+12;
+
for (v2t_cont::iterator it=_adj.begin(); it!=_adj.end(); it++) {
MVertex *ver = it->first;
@@ -381,10 +385,11 @@ void meshMetric::computeMetricHessian()
fullMatrix<double> V(3,3);
fullVector<double> S(3);
H.eig(V,S);
+
- double lambda1 = std::min(std::max(fabs(S(0))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
- double lambda2 = std::min(std::max(fabs(S(1))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
- double lambda3 = (_dim == 3) ? std::min(std::max(fabs(S(2))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin)) : 1.;
+ double lambda1 = std::min(std::max(fabs(S(0))/_epsilonP,1./(hmaxP*hmaxP)),1./(hminP*hminP));
+ double lambda2 = std::min(std::max(fabs(S(1))/_epsilonP,1./(hmaxP*hmaxP)),1./(hminP*hminP));
+ double lambda3 = (_dim == 3) ? std::min(std::max(fabs(S(2))/_epsilonP,1./(hmaxP*hmaxP)),1./(hminP*hminP)) : 1.;
SVector3 t1 (V(0,0),V(1,0),V(2,0));
SVector3 t2 (V(0,1),V(1,1),V(2,1));
@@ -395,10 +400,12 @@ void meshMetric::computeMetricHessian()
_hessian[ver] = H;
setOfSizes[metricNumber].insert(std::make_pair(ver, std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3))));
setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
- setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
+ // setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
}
+ scaleMetric(_epsilon, setOfMetrics[metricNumber]);
+
}
void meshMetric::computeMetricFrey()
@@ -459,7 +466,7 @@ void meshMetric::computeMetricFrey()
_hessian[ver] = hessian;
setOfSizes[metricNumber].insert(std::make_pair(ver, std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3))));
setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
- setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
+ // setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
}
}
@@ -529,7 +536,7 @@ void meshMetric::computeMetricEigenDir()
}
_hessian[ver] = hessian;
setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
- setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
+ // setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
}
@@ -557,7 +564,7 @@ void meshMetric::computeMetricIsoLinInterp()
_hessian[ver] = H;
setOfSizes[metricNumber].insert(std::make_pair(ver, lambda));
setOfMetrics[metricNumber].insert(std::make_pair(ver,H));
- setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(H.determinant())));
+ // setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(H.determinant())));
}
@@ -623,13 +630,76 @@ void meshMetric::computeMetricScaledHessian()
SMetric3 metric(l1,l2,l3,*itT1,*itT2,*itT3);
setOfSizes[metricNumber].insert(std::make_pair(ver, 1./sqrt(lMax)));
setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
- setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
+ // setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
itL1++; itL2++; if (_dim == 3) itL3++;
itT1++; itT2++; itT3++;
}
}
+
+// this function scales the mesh metric in order
+// to reach a target number of elements
+// We know that the number of elements in the final
+// mesh will be (assuming M_e the metric at centroid of element e)
+// N = \sum_e \sqrt {\det (M_e)} V_e
+// where V_e is the volume of e
+// assuming that N_{target} = K N, we have
+// K N = K \sum_e \sqrt {\det (M_e)} V_e
+// = \sum_e \sqrt {K^2 \det (M_e)} V_e
+// = \sum_e \sqrt {\det (K^{2/d} M_e)} V_e
+// where d is the dimension of the problem.
+// This means that the metric should be scaled by K^{2/d} where
+// K is N_target / N
+
+void meshMetric::scaleMetric( int nbElementsTarget,
+ nodalMetricTensor &nmt )
+{
+ // compute N
+ double N = 0;
+ for (unsigned int i=0;i<_elements.size();i++){
+ MElement *e = _elements[i];
+ SMetric3 m1 = nmt[e->getVertex(0)];
+ SMetric3 m2 = nmt[e->getVertex(1)];
+ SMetric3 m3 = nmt[e->getVertex(2)];
+ if (_dim == 2){
+ SMetric3 m = interpolation(m1,m2,m3,0.3333,0.3333);
+ N += sqrt(m.determinant()) * e->getVolume() * 3.0;
+ // printf("%12.5E %12.5E\n",m.determinant(),e->getVolume());
+ }
+ else{
+ SMetric3 m4 = nmt[e->getVertex(3)];
+ SMetric3 m = interpolation(m1,m2,m3,m4,0.25,0.25,0.25);
+ N += sqrt(m.determinant()) * e->getVolume() * 4.0;
+ }
+ }
+ double scale = pow ((double)nbElementsTarget/N,2.0/_dim);
+ // printf("%d elements --- %d element target --- %12.5E elements with the present metric\n",
+ // _elements.size(),nbElementsTarget,N);
+ // getchar();
+ for (nodalMetricTensor::iterator it = nmt.begin(); it != nmt.end() ; ++it){
+ if (_dim == 3){
+ it->second *= scale;
+ }
+ else {
+ it->second(0,0) *= scale;
+ it->second(1,0) *= scale;
+ it->second(1,1) *= scale;
+ }
+ SMetric3 &m = it->second;
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ m.eig(V,S);
+ S(0) = std::min(std::max(S(0),1/(hmax*hmax)),1/(hmin*hmin));
+ S(1) = std::min(std::max(S(1),1/(hmax*hmax)),1/(hmin*hmin));
+ S(2) = std::min(std::max(S(2),1/(hmax*hmax)),1/(hmin*hmin));
+ SVector3 t1 (V(0,0),V(1,0),V(2,0));
+ SVector3 t2 (V(0,1),V(1,1),V(2,1));
+ SVector3 t3 (V(0,2),V(1,2),V(2,2));
+ m = SMetric3(S(0),S(1),S(2),t1,t2,t3);
+ }
+}
+
void meshMetric::computeMetric()
{
//printf("%d elements are considered in the meshMetric \n",(int)_elements.size());
@@ -676,7 +746,9 @@ 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();
+ if (needMetricUpdate) {
+ intersectMetrics();
+ }
if (!setOfMetrics.size()){
std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
throw;
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index b2eed0f41b3de24fba8da470d8c1b796d685e991..3b0d014ac10487239fa017a544e0b133e13195ac 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -6,6 +6,7 @@
#ifndef _MESH_METRIC_H_
#define _MESH_METRIC_H_
+
#include <map>
#include <algorithm>
#include "STensor3.h"
@@ -41,15 +42,16 @@ class meshMetric: public Field {
std::map<MVertex*,double> vals;
std::map<MVertex*,SVector3> grads,dgrads[3];
+ public:
typedef std::map<MVertex*,SMetric3> nodalMetricTensor;
typedef std::map<MVertex*,double> nodalField;
-
+ private:
nodalMetricTensor _nodalMetrics,_hessian;
nodalField _nodalSizes, _detMetric;
std::map<int,nodalMetricTensor> setOfMetrics;
std::map<int,nodalField> setOfSizes;
- std::map<int,nodalField> setOfDetMetric;
+ // std::map<int,nodalField> setOfDetMetric;
public:
meshMetric(std::vector<MElement*> elements);
@@ -81,6 +83,10 @@ class meshMetric: public Field {
if (needMetricUpdate) intersectMetrics();
return _nodalMetrics[v];
}
+ // this function scales the mesh metric in order
+ // to reach a target number of elements
+ void scaleMetric( int nbElementsTarget,
+ nodalMetricTensor &nmt );
void computeMetric();
void computeMetricLevelSet();
diff --git a/benchmarks/python/adaptivity.py b/benchmarks/python/adaptivity.py
index e91fe8d861fec4080c29a9f78d34635e01b99cc4..7d624af1e5efaf12eb1223f2968fd3469c5f480e 100644
--- a/benchmarks/python/adaptivity.py
+++ b/benchmarks/python/adaptivity.py
@@ -19,13 +19,13 @@ def Adaptation ( gm , gm_adapt) :
print "hello 3", myMap.size()
pv = PView ('xcarre','NodeData', gm, myMap)
pve = PViewEvaluator (pv)
- lcMin = 0.00001
+ lcMin = 0.001
lcMax = .1
- eps = 1.e-2
- nbIter = 4
+ NBTARGET = 13000
+ nbIter = 2
allDim = 1
- gm_adapt.adaptMesh([2], [pve], [[eps, lcMin, lcMax]], nbIter, allDim)
+ gm_adapt.adaptMesh([2], [pve], [[NBTARGET, lcMin, lcMax]], nbIter, allDim)
print "hello 4"
gm = GModel()
@@ -38,6 +38,11 @@ gm_adapt.load('square.geo')
gm_adapt2 = GModel()
gm_adapt2.load('square.geo')
+gm_adapt3 = GModel()
+gm_adapt3.load('square.geo')
+
Adaptation ( gm , gm_adapt2)
print "------------------------------------------------------------------------"
Adaptation ( gm_adapt2 , gm_adapt)
+print "------------------------------------------------------------------------"
+Adaptation ( gm_adapt , gm_adapt3)