diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index cd7580629bb81122cc81330c18746b347fd7178f..d1d93fdce441b510a799e938aff1a172c9af8aac 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -548,7 +548,7 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
fields->reset();
meshMetric *metric = new meshMetric(this);
- for (int imetric=0;imetric<technique.size();imetric++){;
+ for (int imetric=0;imetric<technique.size();imetric++){
metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
}
fields->setBackgroundField(metric);
@@ -561,9 +561,9 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
opt_mesh_algo3d(0, GMSH_SET, 7.0); //mmg3d
opt_mesh_lc_from_points(0, GMSH_SET, 0.0); //do not mesh lines with lc
- std::for_each(firstEdge(), lastEdge(), deMeshGEdge());
- std::for_each(firstFace(), lastFace(), deMeshGFace());
std::for_each(firstRegion(), lastRegion(), deMeshGRegion());
+ std::for_each(firstFace(), lastFace(), deMeshGFace());
+ std::for_each(firstEdge(), lastEdge(), deMeshGEdge());
GenerateMesh(this, getDim());
nbElems = getNumMeshElements();
@@ -571,7 +571,7 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
char name[256];
sprintf(name, "meshAdapt-%d.msh", ITER);
writeMSH(name);
-// metric->exportInfo(name);
+ //metric->exportInfo(name);
if (ITER++ >= niter) break;
if (ITER > 3 && fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
@@ -586,27 +586,27 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
std::vector<MElement*> elements;
if (getDim() == 2){
- for (fiter fit = firstFace(); fit != lastFace(); ++fit){
- if ((*fit)->quadrangles.size())return -1;
- for (unsigned i=0;i<(*fit)->triangles.size();i++){
- elements.push_back((*fit)->triangles[i]);
- }
- }
+ for (fiter fit = firstFace(); fit != lastFace(); ++fit){
+ if ((*fit)->quadrangles.size())return -1;
+ for (unsigned i=0;i<(*fit)->triangles.size();i++){
+ elements.push_back((*fit)->triangles[i]);
+ }
+ }
}
else if (getDim() == 3){
- for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
- if ((*rit)->hexahedra.size())return -1;
- for (unsigned i=0;i<(*rit)->tetrahedra.size();i++){
- elements.push_back((*rit)->tetrahedra[i]);
- }
- }
+ for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
+ if ((*rit)->hexahedra.size())return -1;
+ for (unsigned i=0;i<(*rit)->tetrahedra.size();i++){
+ elements.push_back((*rit)->tetrahedra[i]);
+ }
+ }
}
if (elements.size() == 0)return -1;
fields->reset();
meshMetric *metric = new meshMetric(this);
- for (int imetric=0;imetric<technique.size();imetric++){;
+ for (int imetric=0;imetric<technique.size();imetric++){
metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
}
fields->setBackgroundField(metric);
@@ -615,21 +615,21 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
// fields->background_field = id;
if (getDim() == 2){
- for (fiter fit = firstFace(); fit != lastFace(); ++fit){
- if((*fit)->geomType() != GEntity::DiscreteSurface){
- meshGFaceBamg(*fit);
- laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing);
- }
- if(_octree) delete _octree;
- _octree = 0;
- }
+ for (fiter fit = firstFace(); fit != lastFace(); ++fit){
+ if((*fit)->geomType() != GEntity::DiscreteSurface){
+ meshGFaceBamg(*fit);
+ laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing);
+ }
+ if(_octree) delete _octree;
+ _octree = 0;
+ }
}
else if (getDim() == 3){
- for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
- refineMeshMMG(*rit);
- if(_octree) delete _octree;
- _octree = 0;
- }
+ for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
+ refineMeshMMG(*rit);
+ if(_octree) delete _octree;
+ _octree = 0;
+ }
}
nbElems = getNumMeshElements();
diff --git a/Geo/GModel.h b/Geo/GModel.h
index c856ae18d315562deb5429f131ebeaa9f79aab82..3afcc3a34a4fd2537280cb806dafcca85f286b39 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -396,6 +396,7 @@ class GModel
// parameters[0] = thickness of the interface in the positive ls direction (mandatory)
// parameters[4] = thickness of the interface in the negative ls direction (=parameters[0] if not specified)
// parameters[3] = the required minimum number of elements to represent a circle - used for curvature-based metric (default: = 15)
+ // 5) Same as 4, except that the transition in band E uses linear interpolation of h, instead of linear interpolation of metric
// The algorithm first generate a mesh if no one is available
// In this first attempt, only the highest dimensional mesh is adapted, which is ok if
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 766628d747dfcc5049a6c570e6c1328a5604a8d1..ce465eb00d98b256f2d4ccd86cee330deb47a3f0 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -73,6 +73,7 @@ void meshMetric::addMetric(int technique, simpleFunction<double> *fct, std::vect
_E = parameters[0];
_E_moins = (parameters.size() >= 5) ? parameters[4] : -parameters[0];
+ if (_E_moins>0.) _E_moins *= -1.;
_epsilon = parameters[0];
_technique = (MetricComputationTechnique)technique;
_Np = (parameters.size() >= 4) ? parameters[3] : 15.;
@@ -342,7 +343,7 @@ void meshMetric::computeMetric(){
// else H = hfrey;
H = hfrey;
}
- else if (_technique == meshMetric::EIGENDIRECTIONS ){
+ else if ((_technique == meshMetric::EIGENDIRECTIONS )||(_technique == meshMetric::EIGENDIRECTIONS_LINEARINTERP_H )){
double metric_value_hmax = 1./hmax/hmax;
SVector3 gr = grads[ver];
double norm = gr.normalize();
@@ -369,12 +370,16 @@ void meshMetric::computeMetric(){
eigenvals_curvature.push_back(fabs(lambda_hessian(2))*un_sur_epsilon);
// metric distance-based eigenvalue, in gradient direction
double metric_value_hmin = 1./hmin/hmin;
- // 1/ linear interpolation between h_min and h_max...
- //double h_dist = hmin + ((hmax-hmin)/_E)*dist;// the charcteristic element size in the normal direction - linear interp between hmin and hmax
- //eigenval_direction = 1./h_dist/h_dist;
- // 2/ ... or linear interpolation between 1/h_min^2 and 1/h_max^2
- double maximum_distance = (signed_dist>0.) ? _E : fabs(_E_moins);
- double eigenval_direction = metric_value_hmin + ((metric_value_hmax-metric_value_hmin)/maximum_distance)*dist;
+ double eigenval_direction;
+ if (_technique==meshMetric::EIGENDIRECTIONS_LINEARINTERP_H){
+ double h_dist = hmin + ((hmax-hmin)/_E)*dist;// the charcteristic element size in the normal direction - linear interp between hmin and hmax
+ eigenval_direction = 1./h_dist/h_dist;
+ }
+ else if(_technique==meshMetric::EIGENDIRECTIONS){
+ // ... or linear interpolation between 1/h_min^2 and 1/h_max^2
+ double maximum_distance = (signed_dist>0.) ? _E : fabs(_E_moins);
+ eigenval_direction = metric_value_hmin + ((metric_value_hmax-metric_value_hmin)/maximum_distance)*dist;
+ }
// now, consider only three eigenvectors (i.e. (grad(LS), u1, u2) with u1 and u2 orthogonal vectors in the tangent plane to the LS)
// and imposing directly eigenvalues and directions, instead of intersecting the two metrics based on direction and curvature
@@ -413,7 +418,7 @@ void meshMetric::computeMetric(){
it++;
}
- if (_technique != meshMetric::EIGENDIRECTIONS ){
+ if (_technique != meshMetric::EIGENDIRECTIONS && _technique!=meshMetric::EIGENDIRECTIONS_LINEARINTERP_H){
//See paper Hachem and Coupez:
//Finite element solution to handle complex heat and fluid flows in
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 820b3a148edfe2f21b257e7e75a928bf389a3459..238338df20296714385b9b6a3772f93057852dbb 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -14,7 +14,7 @@ class STensor3;
/**Anisotropic mesh size field based on a metric */
class meshMetric: public Field {
public:
- typedef enum {LEVELSET=1,HESSIAN=2, FREY=3, EIGENDIRECTIONS=4} MetricComputationTechnique;
+ typedef enum {LEVELSET=1,HESSIAN=2, FREY=3, EIGENDIRECTIONS=4, EIGENDIRECTIONS_LINEARINTERP_H=5} MetricComputationTechnique;
private:
// intersect all metrics added in "setOfMetrics", preserve eigendirections of the "most anisotropic" metric
void intersectMetrics();
@@ -63,6 +63,7 @@ class meshMetric: public Field {
// parameters[0] = thickness of the interface in the positive ls direction (mandatory)
// parameters[4] = thickness of the interface in the negative ls direction (default: =parameters[0] if not specified)
// parameters[3] = the required minimum number of elements to represent a circle - used for curvature-based metric (default: = 15)
+ // 5: same as 4, except that the transition in band E uses linear interpolation of h, instead of linear interpolation of metric
void addMetric(int technique, simpleFunction<double> *fct, std::vector<double> parameters);
inline SMetric3 metricAtVertex (MVertex* v) {