From 3728935c5b2bf0f4259380ca3c42ecadb9eca482 Mon Sep 17 00:00:00 2001
From: Paul-Emile Bernard <paul-emile.bernard@uclouvain.be>
Date: Tue, 7 Feb 2012 09:29:10 +0000
Subject: [PATCH] Add a new metric computation for level set, simply based on
ls gradient and hessian eigendirections. meshAdapt was modified to receive
and intersect multiple metrics, preserving the most anisotropic metric. Also
added input and output functions to dgDofContainer, for the "NodeData"
format. Benchmark for anisotropic adaptation to level set iso-zero created in
benchmarks_gmsh/LSIsozeroAniso.
---
Geo/GModel.cpp | 17 +-
Geo/GModel.h | 19 +-
Geo/STensor3.cpp | 23 ++
Geo/STensor3.h | 3 +
Geo/gmshLevelset.cpp | 56 ++++
Geo/gmshLevelset.h | 16 ++
Mesh/meshMetric.cpp | 614 +++++++++++++++++++++++++++----------------
Mesh/meshMetric.h | 57 +++-
gmshpy/gmshGeo.i | 8 +
9 files changed, 571 insertions(+), 242 deletions(-)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e2752279f5..9d1bca0405 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -524,15 +524,13 @@ int GModel::mesh(int dimension)
#endif
}
-int GModel::adaptMesh(int technique, simpleFunction<double> *f,
- std::vector<double> parameters, bool meshAll)
+int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<double>* > f, std::vector<std::vector<double> > parameters, int niter, bool meshAll)
{
#if defined(HAVE_MESH)
if (getNumMeshElements() == 0) mesh(getDim());
int nbElemsOld = getNumMeshElements();
int nbElems;
- int niter = parameters.size() >=4 ? (int) parameters[3] : 3;
FieldManager *fields = getFields();
fields->reset();
@@ -544,7 +542,11 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f,
Msg::Info("-- adaptMesh (allDim) ITER =%d ", ITER);
fields->reset();
- fields->setBackgroundField(new meshMetric(this, technique, f, parameters));
+ meshMetric *metric = new meshMetric(this);
+ for (int imetric=0;imetric<technique.size();imetric++){;
+ metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
+ }
+ fields->setBackgroundField(metric);
// int id = fields->newId();
// (*fields)[id] = new meshMetric(this, technique, f, parameters);
// fields->background_field = id;
@@ -564,6 +566,7 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f,
char name[256];
sprintf(name, "meshAdapt-%d.msh", ITER);
writeMSH(name);
+// metric->exportInfo(name);
if (ITER++ >= niter) break;
if (ITER > 3 && fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
@@ -597,7 +600,11 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f,
if (elements.size() == 0)return -1;
fields->reset();
- fields->setBackgroundField(new meshMetric(this, technique, f, parameters));
+ meshMetric *metric = new meshMetric(this);
+ for (int imetric=0;imetric<technique.size();imetric++){;
+ metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
+ }
+ fields->setBackgroundField(metric);
// int id = fields->newId();
// (*fields)[id] = new meshMetric(this, technique, f, parameters);
// fields->background_field = id;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 2196de1c5d..c856ae18d3 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -376,23 +376,32 @@ class GModel
// mesh the model
int mesh(int dimension);
- // adapt the mesh anisotropically using a metric that is computed from a scalar function f(x,y,z).
- // One can either
+ // adapt the mesh anisotropically using metrics that are computed from a set of functions f(x,y,z).
// For all algorithms
// parameters[1] = lcmin (default : in global gmsh options CTX::instance()->mesh.lcMin)
// parameters[2] = lcmax (default : in global gmsh options CTX::instance()->mesh.lcMax)
- // parameters[3] = nb iterations
+ // niter is the maximum number of iterations
+ // Available algorithms:
// 1) Assume that the function is a levelset -> adapt using Coupez technique (technique = 1)
// parameters[0] = thickness of the interface (mandatory)
// 2) Assume that the function is a physical quantity -> adapt using the Hessain (technique = 2)
// parameters[0] = N, the final number of elements
// 3) A variant of 1) by P. Frey (= Coupez + takes curvature function into account)
// parameters[0] = thickness of the interface (mandatory)
- // The algorithm first generate a mesh if no one is available
+ // parameters[3] = the required minimum number of elements to represent a circle - used for curvature-based metric (default: = 15)
+ // 4) A variant (3), direct implementation in the metric eigendirections, assuming a level set (ls):
+ // - hmin is imposed in the ls gradient,
+ // - hmax is imposed in the two eigendirections of the ls hessian that are (almost ?) tangent to the iso-zero plane
+ // + the latter eigenvalues (1/hmax^2) are modified if necessary to capture the iso-zero curvature
+ // 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)
+ // 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
// we assume that boundaries are already adapted.
// This should be fixed.
- int adaptMesh (int technique, simpleFunction<double> *f, std::vector<double> parameters, bool meshAll=false);
+ int adaptMesh(std::vector<int> technique, std::vector<simpleFunction<double>*> f, std::vector<std::vector<double> > parameters, int niter, bool meshAll=false);
// make the mesh a high order mesh at order N
// linear is 1 if the high order points are not placed on the geometry of the model
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 7436bc7a51..0bbe5ee232 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -49,6 +49,29 @@ SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2)
return iv;
}
+// preserve orientation of the most anisotropic metric !!!
+SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2)
+{
+ fullMatrix<double> V1(3,3);
+ fullVector<double> S1(3);
+ m1.eig(V1,S1,true);
+ double lambda1_min = std::min(std::min(fabs(S1(0)),fabs(S1(1))),fabs(S1(2)));
+ double lambda1_max = std::max(std::max(fabs(S1(0)),fabs(S1(1))),fabs(S1(2)));
+ fullMatrix<double> V2(3,3);
+ fullVector<double> S2(3);
+ m2.eig(V2,S2,true);
+ double lambda2_min = std::min(std::min(fabs(S2(0)),fabs(S2(1))),fabs(S2(2)));
+ double lambda2_max = std::max(std::max(fabs(S2(0)),fabs(S2(1))),fabs(S2(2)));
+
+ double ratio1 = lambda1_min/lambda1_max;
+ double ratio2 = lambda2_min/lambda2_max;
+
+ if (ratio1<ratio2)
+ return intersection_conserveM1(m1,m2);
+ else
+ return intersection_conserveM1(m2,m1);
+}
+
// (1-t) * m1 + t * m2
SMetric3 interpolation (const SMetric3 &m1,
const SMetric3 &m2,
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index 6673cc5908..c99db1b89c 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -163,8 +163,11 @@ inline double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
b.z() * ( m(0,2) * a.x() + m(1,2) * a.y() + m(2,2) * a.z() ) ;
}
+// preserve orientation of m1
SMetric3 intersection_conserveM1 (const SMetric3 &m1,
const SMetric3 &m2);
+// preserve orientation of the most anisotropic metric
+SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2);
// compute the largest inscribed ellipsoid...
SMetric3 intersection (const SMetric3 &m1,
const SMetric3 &m2);
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index bed67951e9..f0515ed706 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -631,6 +631,62 @@ double gLevelsetQuadric::operator()(const double x, const double y, const double
+ 2. * A[1][2] * y * z + A[2][2] * z * z + B[0] * x + B[1] * y + B[2] * z + C);
}
+gLevelsetShamrock::gLevelsetShamrock(double _xmid, double _ymid, double _zmid, double _a, double _b, int _c, int tag) : gLevelsetPrimitive(tag), xmid(_xmid), ymid(_ymid), zmid(_zmid), a(_a), b(_b), c(_c) {
+ // creating the iso-zero
+ double angle = 0.;
+ double r;
+ while (angle<=2.*M_PI){
+ r = a+b*sin(c*angle);
+ iso_x.push_back(r*sin(angle)+xmid);
+ iso_y.push_back(r*cos(angle)+xmid);
+
+ angle += 2.*M_PI/1000.;
+ }
+}
+
+double gLevelsetShamrock::operator() (const double x, const double y, const double z) const {
+ // computing distance to pre-defined (sampled) iso-zero
+ double dx,dy,xi,yi,d;
+ dx = x-iso_x[0];
+ dy = y-iso_y[0];
+ double distance = sqrt(dx*dx+dy*dy);
+ std::vector<double>::const_iterator itx = iso_x.begin();
+ std::vector<double>::const_iterator itxen = iso_x.end();
+ std::vector<double>::const_iterator ity = iso_y.begin();
+ std::vector<double>::const_iterator itminx = iso_x.begin();
+ std::vector<double>::const_iterator itminy = iso_y.begin();
+ itx++;ity++;
+ for (;itx!=itxen;itx++,ity++){
+ xi = *itx;
+ yi = *ity;
+ dx = x-xi;
+ dy = y-yi;
+ d = sqrt(dx*dx+dy*dy);
+ if (distance>d){
+ distance = d;
+ itminx = itx;
+ itminy = ity;
+ }
+ }
+ // still need to find the LS sign... using vectorial prod with tangent vector
+ // if not, the LS gradient is discontinuous on iso-zero... oups !
+ SVector3 t,p;
+ p(0) = (*itminx) - x;
+ p(1) = (*itminy) - y;
+ if (itminx==(itxen-1)){
+ t(0) = iso_x[0] - (*itminx);
+ t(1) = iso_y[0] - (*itminy);
+ }
+ else{
+ t(0) = (*(itminx+1)) - (*itminx);
+ t(1) = (*(itminy+1)) - (*itminy);
+ }
+ double vectprod = (p(0)*t(1) - p(1)*t(0));
+ double sign = (vectprod < 0.) ? -1. : 1.;
+
+ return (sign*distance);
+}
+
gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r0, double _A, double _sigma, int tag) : gLevelsetPrimitive(tag) {
A = _A;
sigma = _sigma;
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 590e25d40b..2712cc2163 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -267,6 +267,22 @@ public:
int type() const { return UNKNOWN; }
};
+// creates the 2D (-approximate- signed distance !) level set
+// corresponding to the "shamrock-like" iso-zero from
+// Dobrzynski and Frey, "Anisotropic delaunay mesh adaptation for unsteady simulations",
+// 17th International Meshing Rountable (2008)(177–194)
+class gLevelsetShamrock: public gLevelsetPrimitive
+{
+ double xmid, ymid, zmid,a,b;
+ int c;
+ std::vector<double> iso_x,iso_y;
+public:
+ gLevelsetShamrock(double xmid, double ymid, double zmid, double a, double b, int c=3, int tag=1);
+ ~gLevelsetShamrock(){}
+ double operator () (const double x, const double y, const double z) const;
+ int type() const { return UNKNOWN; }
+};
+
class gLevelsetMathEval: public gLevelsetPrimitive
{
mathEvaluator *_expr;
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index eb690c3f92..766628d747 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -1,3 +1,5 @@
+
+
#include "meshMetric.h"
#include "meshGFaceOptimize.h"
#include "Context.h"
@@ -14,39 +16,16 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l
for (int j=0;j<lt[i]->getNumVertices();j++){
MVertex *v1 = lt[i]->getVertex(j);
if (v1 != v){
- std::vector<MElement*> <2 = adj[v1];
- stencil.insert(lt2.begin(),lt2.end());
+ std::vector<MElement*> <2 = adj[v1];
+ stencil.insert(lt2.begin(),lt2.end());
}
}
}
lt.clear();
lt.insert(lt.begin(),stencil.begin(),stencil.end());
}
-meshMetric::meshMetric(std::vector<MElement*> elements, int technique, simpleFunction<double> *fct, std::vector<double> parameters): _fct(fct){
-
- _dim = elements[0]->getDim();
- std::map<MElement*, MElement*> newP;
- std::map<MElement*, MElement*> newD;
-
- for (int i=0;i<elements.size();i++){
- MElement *e = elements[i];
- MElement *copy = e->copy(_vertexMap, newP, newD);
- _elements.push_back(copy);
- }
-
- _octree = new MElementOctree(_elements);
- hmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
- hmax = parameters.size() >=3 ? parameters[2] : CTX::instance()->mesh.lcMax;
-
- _E = parameters[0];
- _epsilon = parameters[0];
- _technique = (MetricComputationTechnique)technique;
-
- computeMetric();
-}
-
-meshMetric::meshMetric(GModel *gm, int technique, simpleFunction<double> *fct, std::vector<double> parameters): _fct(fct){
+meshMetric::meshMetric(GModel *gm){
_dim = gm->getDim();
std::map<MElement*, MElement*> newP;
std::map<MElement*, MElement*> newD;
@@ -54,40 +33,142 @@ meshMetric::meshMetric(GModel *gm, int technique, simpleFunction<double> *fct, s
if (_dim == 2){
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); ++fit){
for (int i=0;i<(*fit)->getNumMeshElements();i++){
- MElement *e = (*fit)->getMeshElement(i);
- MElement *copy = e->copy(_vertexMap, newP, newD);
- _elements.push_back(copy);
+ MElement *e = (*fit)->getMeshElement(i);
+ MElement *copy = e->copy(_vertexMap, newP, newD);
+ _elements.push_back(copy);
}
}
}
else if (_dim == 3){
for (GModel::riter rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit){
for (int i=0;i<(*rit)->getNumMeshElements();i++){
- MElement *e = (*rit)->getMeshElement(i);
- MElement *copy = e->copy(_vertexMap, newP, newD);
- _elements.push_back(copy);
+ MElement *e = (*rit)->getMeshElement(i);
+ MElement *copy = e->copy(_vertexMap, newP, newD);
+ _elements.push_back(copy);
}
}
}
+ _octree = new MElementOctree(_elements);
+}
+
+meshMetric::meshMetric(std::vector<MElement*> elements){
+ _dim = elements[0]->getDim();
+ std::map<MElement*, MElement*> newP;
+ std::map<MElement*, MElement*> newD;
+
+ for (int i=0;i<elements.size();i++){
+ MElement *e = elements[i];
+ MElement *copy = e->copy(_vertexMap, newP, newD);
+ _elements.push_back(copy);
+ }
_octree = new MElementOctree(_elements);
+}
+
+void meshMetric::addMetric(int technique, simpleFunction<double> *fct, std::vector<double> parameters){
+ needMetricUpdate=true;
+ _fct = fct;
hmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
hmax = parameters.size() >=3 ? parameters[2] : CTX::instance()->mesh.lcMax;
-
+
_E = parameters[0];
+ _E_moins = (parameters.size() >= 5) ? parameters[4] : -parameters[0];
_epsilon = parameters[0];
_technique = (MetricComputationTechnique)technique;
-
+ _Np = (parameters.size() >= 4) ? parameters[3] : 15.;
+
computeMetric();
+}
+
+void meshMetric::intersectMetrics(){
+ if (!setOfMetrics.size()){
+ std::cout << " meshMetric::intersectMetrics: Can't intersect metrics, no metric registered ! " << std::endl;
+ return;
+ }
+
+ v2t_cont adj;
+ buildVertexToElement (_elements,adj);
+ v2t_cont :: iterator it = adj.begin();
+ for (;it != adj.end();it++) {
+ MVertex *ver = it->first;
+ _nodalMetrics[ver] = setOfMetrics[0][ver];
+ _detMetric[ver] = setOfDetMetric[0][ver];
+ _nodalSizes[ver] = setOfSizes[0][ver];
+ for (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;
- printMetric("toto.pos");
- //exit(1);
+}
+void meshMetric::exportInfo(const char * fileendname){
+ if (needMetricUpdate) intersectMetrics();
+ std::stringstream sg,sm,sl;
+ sg << "meshmetric_gradients" << fileendname;
+ sm << "meshmetric_metric" << fileendname;
+ sl << "meshmetric_levelset" << fileendname;
+ std::ofstream out_grad(sg.str().c_str());
+ std::ofstream out_metric(sm.str().c_str());
+ std::ofstream out_ls(sl.str().c_str());
+ out_grad << "View \"ls_gradient\"{" << std::endl;
+ out_metric << "View \"metric\"{" << std::endl;
+ out_ls << "View \"ls\"{" << std::endl;
+ std::vector<MElement*>::iterator itelem = _elements.begin();
+ std::vector<MElement*>::iterator itelemen = _elements.end();
+ for (;itelem!=itelemen;itelem++){
+ MElement* e = *itelem;
+ out_metric << "TT(";
+ out_grad << "VT(";
+ out_ls << "ST(";
+ for ( int i = 0; i < e->getNumVertices(); i++) {
+ MVertex *ver = e->getVertex(i);
+ out_metric << ver->x() << "," << ver->y() << "," << ver->z();
+ out_grad << ver->x() << "," << ver->y() << "," << ver->z();
+ out_ls << ver->x() << "," << ver->y() << "," << ver->z();
+ if (i!=e->getNumVertices()-1){
+ out_metric << ",";
+ out_grad << ",";
+ out_ls << ",";
+ }
+ else{
+ out_metric << "){";
+ out_grad << "){";
+ out_ls << "){";
+ }
+ }
+ for ( int i = 0; i < e->getNumVertices(); i++) {
+ MVertex *ver = e->getVertex(i);
+ out_ls << vals[ver];
+ if ((i==(e->getNumVertices()-1))) out_ls << "};" << std::endl;
+ else out_ls << ",";
+ for (int k=0;k<3;k++){
+ out_grad << grads[ver](k);
+ if ((k==2)&&(i==(e->getNumVertices()-1))) out_grad << "};" << std::endl;
+ else out_grad << ",";
+ for (int l=0;l<3;l++){
+ out_metric << _nodalMetrics[ver](k,l);
+ if ((k==2)&&(l==2)&&(i==(e->getNumVertices()-1))) out_metric << "};" << std::endl;
+ else out_metric << ",";
+ }
+ }
+ }
+ }
+ out_grad << "};" << std::endl;
+ out_metric << "};" << std::endl;
+ out_ls << "};" << std::endl;
+ out_grad.close();
+ out_metric.close();
+ out_ls.close();
}
+
+
meshMetric::~meshMetric(){
if (_octree) delete _octree;
for (int i=0; i< _elements.size(); i++)
- delete _elements[i];
+ delete _elements[i];
}
void meshMetric::computeValues( v2t_cont adj){
@@ -116,7 +197,7 @@ void meshMetric::computeHessian( v2t_cont adj){
// increaseStencil(ver,adj,lt);
// }
// }
-
+
fullMatrix<double> A (lt.size(),DIM);
fullMatrix<double> AT (DIM,lt.size());
fullMatrix<double> ATA (DIM,DIM);
@@ -125,51 +206,51 @@ void meshMetric::computeHessian( v2t_cont adj){
fullVector<double> result (DIM);
fullMatrix<double> f (1,1);
for(int i = 0; i < lt.size(); i++) {
- MElement *e = lt[i];
- int npts; IntPt *pts;
- SPoint3 p;
- e->getIntegrationPoints(0,&npts,&pts);
- if (ITER == 0){
- SPoint3 p;
- e->pnt(pts[0].pt[0],pts[0].pt[1],pts[0].pt[2],p);
- b(i) = (*_fct)(p.x(),p.y(),p.z());
- }
- else {
- double value = 0.0;
- for (int j=0;j<e->getNumVertices();j++){
- SVector3 gr = grads[e->getVertex(j)];
- value += gr(ITER-1);
- }
- b(i) = value/(double)e->getNumVertices();
- }
- e->pnt(pts[0].pt[0],pts[0].pt[1],pts[0].pt[2],p);
- A(i,0) = 1.0;
- A(i,1) = p.x() - ver->x();
- A(i,2) = p.y() - ver->y();
- if (_dim == 3) A(i,3) = p.z() - ver->z();
+ MElement *e = lt[i];
+ int npts; IntPt *pts;
+ SPoint3 p;
+ e->getIntegrationPoints(0,&npts,&pts);
+ if (ITER == 0){
+ SPoint3 p;
+ e->pnt(pts[0].pt[0],pts[0].pt[1],pts[0].pt[2],p);
+ b(i) = (*_fct)(p.x(),p.y(),p.z());
+ }
+ else {
+ double value = 0.0;
+ for (int j=0;j<e->getNumVertices();j++){
+ SVector3 gr = grads[e->getVertex(j)];
+ value += gr(ITER-1);
+ }
+ b(i) = value/(double)e->getNumVertices();
+ }
+ e->pnt(pts[0].pt[0],pts[0].pt[1],pts[0].pt[2],p);
+ A(i,0) = 1.0;
+ A(i,1) = p.x() - ver->x();
+ A(i,2) = p.y() - ver->y();
+ if (_dim == 3) A(i,3) = p.z() - ver->z();
}
AT = A.transpose();
AT.mult(A,ATA);
AT.mult(b,ATb);
ATA.luSolve(ATb,result);
if (ITER == 0){
- double gr1 = result(1);
- double gr2 = result(2);
- double gr3 = (_dim==2) ? 0.0:result(3);
- double norm = sqrt(gr1*gr1+gr2*gr2+gr3*gr3);
- if (norm == 0.0 || _technique == meshMetric::HESSIAN) norm = 1.0;
- grads[ver] = SVector3(gr1/norm,gr2/norm,gr3/norm);
+ double gr1 = result(1);
+ double gr2 = result(2);
+ double gr3 = (_dim==2) ? 0.0:result(3);
+ double norm = sqrt(gr1*gr1+gr2*gr2+gr3*gr3);
+ if (norm == 0.0 || _technique == meshMetric::HESSIAN) norm = 1.0;
+ grads[ver] = SVector3(gr1/norm,gr2/norm,gr3/norm);
}
else
- dgrads[ITER-1][ver] = SVector3(result(1),result(2),_dim==2? 0.0:result(3));
+ dgrads[ITER-1][ver] = SVector3(result(1),result(2),_dim==2? 0.0:result(3));
++it;
}
if (_technique == meshMetric::LEVELSET) break;
}
-
+
}
void meshMetric::computeMetric(){
-
+
v2t_cont adj;
buildVertexToElement (_elements,adj);
@@ -178,38 +259,41 @@ void meshMetric::computeMetric(){
computeValues(adj);
computeHessian(adj);
+ int metricNumber = setOfMetrics.size();
+
v2t_cont :: iterator it = adj.begin();
while (it != adj.end()) {
- MVertex *ver = it->first;
- double dist = fabs(vals[ver]);
-
- SMetric3 H;
- SVector3 gradudx = dgrads[0][ver];
- SVector3 gradudy = dgrads[1][ver];
- SVector3 gradudz = dgrads[2][ver];
- SMetric3 hessian;
- hessian(0,0) = gradudx(0);
- hessian(1,1) = gradudy(1);
- hessian(2,2) = gradudz(2);
- hessian(1,0) = hessian(0,1) = 0.5 * (gradudx(1) +gradudy(0));
- hessian(2,0) = hessian(0,2) = 0.5 * (gradudx(2) +gradudz(0));
- hessian(2,1) = hessian(1,2) = 0.5 * (gradudy(2) +gradudz(1));
- if (_technique == meshMetric::HESSIAN){
+ MVertex *ver = it->first;
+ double signed_dist = vals[ver];
+ double dist = fabs(signed_dist);
+
+ SMetric3 H;
+ SVector3 gradudx = dgrads[0][ver];
+ SVector3 gradudy = dgrads[1][ver];
+ SVector3 gradudz = dgrads[2][ver];
+ SMetric3 hessian;
+ hessian(0,0) = gradudx(0);
+ hessian(1,1) = gradudy(1);
+ hessian(2,2) = gradudz(2);
+ hessian(1,0) = hessian(0,1) = 0.5 * (gradudx(1) +gradudy(0));
+ hessian(2,0) = hessian(0,2) = 0.5 * (gradudx(2) +gradudz(0));
+ hessian(2,1) = hessian(1,2) = 0.5 * (gradudy(2) +gradudz(1));
+ if (_technique == meshMetric::HESSIAN){
H = hessian;
}
- else if (_technique == meshMetric::LEVELSET ){
+ else if (_technique == meshMetric::LEVELSET ){
SVector3 gr = grads[ver];
SMetric3 hlevelset(1./(hmax*hmax));
double norm = gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2);
if (dist < _E && norm != 0.0){
- double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
- double C = 1./(h*h) -1./(hmax*hmax);
- hlevelset(0,0) += C*gr(0)*gr(0)/norm ;
- hlevelset(1,1) += C*gr(1)*gr(1)/norm ;
- hlevelset(2,2) += C*gr(2)*gr(2)/norm ;
- hlevelset(1,0) = hlevelset(0,1) = C*gr(1)*gr(0)/norm ;
- hlevelset(2,0) = hlevelset(0,2) = C*gr(2)*gr(0)/norm ;
- hlevelset(2,1) = hlevelset(1,2) = C*gr(2)*gr(1)/norm ;
+ double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
+ double C = 1./(h*h) -1./(hmax*hmax);
+ hlevelset(0,0) += C*gr(0)*gr(0)/norm ;
+ hlevelset(1,1) += C*gr(1)*gr(1)/norm ;
+ hlevelset(2,2) += C*gr(2)*gr(2)/norm ;
+ hlevelset(1,0) = hlevelset(0,1) = C*gr(1)*gr(0)/norm ;
+ hlevelset(2,0) = hlevelset(0,2) = C*gr(2)*gr(0)/norm ;
+ hlevelset(2,1) = hlevelset(1,2) = C*gr(2)*gr(1)/norm ;
}
// else if (dist > _E && dist < 2.*_E && norm != 0.0){
// double hmid = (hmax-hmin)/(1.*_E)*dist+(2.*hmin-1.*hmax);
@@ -223,90 +307,163 @@ void meshMetric::computeMetric(){
// }
H = hlevelset;
}
- //See paper Ducrot and Frey:
- //Anisotropic levelset adaptation for accurate interface capturing,
- //ijnmf, 2010
- else if (_technique == meshMetric::FREY ){
- SVector3 gr = grads[ver];
- SMetric3 hfrey(1./(hmax*hmax));
- double norm = gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2);
- if (dist < _E && norm != 0.0){
- double h = hmin*(hmax/hmin-1.0)*dist/_E + hmin;
- double C = 1./(h*h) -1./(hmax*hmax);
- double kappa = hessian(0,0)+hessian(1,1)+hessian(2,2);
- double Np = 15.0;
- double epsGeom = 4.0*3.14*3.14/(kappa*Np);
- hfrey(0,0) += C*gr(0)*gr(0)/(norm) + hessian(0,0)/epsGeom;
- hfrey(1,1) += C*gr(1)*gr(1)/(norm) + hessian(1,1)/epsGeom;
- hfrey(2,2) += C*gr(2)*gr(2)/(norm) + hessian(2,2)/epsGeom;
- hfrey(1,0) = hfrey(0,1) = C*gr(1)*gr(0)/(norm) + hessian(1,0)/epsGeom;
- hfrey(2,0) = hfrey(0,2) = C*gr(2)*gr(0)/(norm) + hessian(2,0)/epsGeom;
- hfrey(2,1) = hfrey(1,2) = C*gr(2)*gr(1)/(norm) + hessian(2,1)/epsGeom;
- // hfrey(0,0) += C*gr(0)*gr(0)/norm;
- // hfrey(1,1) += C*gr(1)*gr(1)/norm;
- // hfrey(2,2) += C*gr(2)*gr(2)/norm;
- // hfrey(1,0) = hfrey(0,1) = gr(1)*gr(0)/(norm) ;
- // hfrey(2,0) = hfrey(0,2) = gr(2)*gr(0)/(norm) ;
- // hfrey(2,1) = hfrey(1,2) = gr(2)*gr(1)/(norm) ;
- }
- // SMetric3 sss=hessian;
- // sss *= divEps;
- // sss(0,0) += 1/(hmax*hmax);
- // sss(1,1) += 1/(hmax*hmax);
- // sss(2,2) += 1/(hmax*hmax);
- // H = intersection(sss,hfrey);
- // if (dist < _E) H = intersection(sss,hfrey);
- // else H = hfrey;
- H = hfrey;
- }
- //See paper Hachem and Coupez:
- //Finite element solution to handle complex heat and fluid flows in
- //industrial furnaces using the immersed volume technique, ijnmf, 2010
-
- fullMatrix<double> V(3,3);
- fullVector<double> S(3);
- H.eig(V,S);
-
- double lambda1, lambda2, lambda3;
- // if (dist < _E && _technique == meshMetric::FREY){
- // fullMatrix<double> Vhess(3,3);
- // fullVector<double> Shess(3);
- // hessian.eig(Vhess,Shess);
- // double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
- // double lam1 = Shess(0);
- // double lam2 = Shess(1);
- // double lam3 = (_dim == 3)? Shess(2) : 1.;
- // lambda1 = lam1;
- // lambda2 = lam2/lambda1;
- // lambda3 = (_dim == 3)? lam3/lambda1: 1.0;
- // }
- // else{
+ //See paper Ducrot and Frey:
+ //Anisotropic levelset adaptation for accurate interface capturing,
+ //ijnmf, 2010
+ else if (_technique == meshMetric::FREY ){
+ SVector3 gr = grads[ver];
+ SMetric3 hfrey(1./(hmax*hmax));
+ double norm = gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2);
+ if (dist < _E && norm != 0.0){
+ double h = hmin*(hmax/hmin-1.0)*dist/_E + hmin;
+ double C = 1./(h*h) -1./(hmax*hmax);
+ double kappa = hessian(0,0)+hessian(1,1)+hessian(2,2);
+ double epsGeom = 4.0*3.14*3.14/(kappa*_Np*_Np);
+ hfrey(0,0) += C*gr(0)*gr(0)/(norm) + hessian(0,0)/epsGeom;
+ hfrey(1,1) += C*gr(1)*gr(1)/(norm) + hessian(1,1)/epsGeom;
+ hfrey(2,2) += C*gr(2)*gr(2)/(norm) + hessian(2,2)/epsGeom;
+ hfrey(1,0) = hfrey(0,1) = C*gr(1)*gr(0)/(norm) + hessian(1,0)/epsGeom;
+ hfrey(2,0) = hfrey(0,2) = C*gr(2)*gr(0)/(norm) + hessian(2,0)/epsGeom;
+ hfrey(2,1) = hfrey(1,2) = C*gr(2)*gr(1)/(norm) + hessian(2,1)/epsGeom;
+ // hfrey(0,0) += C*gr(0)*gr(0)/norm;
+ // hfrey(1,1) += C*gr(1)*gr(1)/norm;
+ // hfrey(2,2) += C*gr(2)*gr(2)/norm;
+ // hfrey(1,0) = hfrey(0,1) = gr(1)*gr(0)/(norm) ;
+ // hfrey(2,0) = hfrey(0,2) = gr(2)*gr(0)/(norm) ;
+ // hfrey(2,1) = hfrey(1,2) = gr(2)*gr(1)/(norm) ;
+ }
+ // SMetric3 sss=hessian;
+ // sss *= divEps;
+ // sss(0,0) += 1/(hmax*hmax);
+ // sss(1,1) += 1/(hmax*hmax);
+ // sss(2,2) += 1/(hmax*hmax);
+ // H = intersection(sss,hfrey);
+ // if (dist < _E) H = intersection(sss,hfrey);
+ // else H = hfrey;
+ H = hfrey;
+ }
+ else if (_technique == meshMetric::EIGENDIRECTIONS ){
+ double metric_value_hmax = 1./hmax/hmax;
+ SVector3 gr = grads[ver];
+ double norm = gr.normalize();
+ SMetric3 metric;
+
+ if (signed_dist < _E && signed_dist > _E_moins && norm != 0.0){
+ // first, compute the eigenvectors of the hessian, for a curvature-based metric
+ fullMatrix<double> eigvec_hessian(3,3);
+ fullVector<double> lambda_hessian(3);
+ hessian.eig(eigvec_hessian,lambda_hessian);
+ std::vector<SVector3> ti;
+ ti.push_back(SVector3(eigvec_hessian(0,0),eigvec_hessian(1,0),eigvec_hessian(2,0)));
+ ti.push_back(SVector3(eigvec_hessian(0,1),eigvec_hessian(1,1),eigvec_hessian(2,1)));
+ ti.push_back(SVector3(eigvec_hessian(0,2),eigvec_hessian(1,2),eigvec_hessian(2,2)));
+
+ // compute the required characteristic length relative to the curvature and the distance to iso
+ double maxkappa = std::max(std::max(fabs(lambda_hessian(0)),fabs(lambda_hessian(1))),fabs(lambda_hessian(2)));
+ // epsilon is set such that the characteristic element size in the tangent direction is h = 2pi R_min/_Np = sqrt(1/metric_maxmax) = sqrt(epsilon/kappa_max)
+ double un_sur_epsilon = maxkappa/((2.*3.14/_Np)*(2.*3.14/_Np));
+ // metric curvature-based eigenvalues
+ std::vector<double> eigenvals_curvature;
+ eigenvals_curvature.push_back(fabs(lambda_hessian(0))*un_sur_epsilon);
+ eigenvals_curvature.push_back(fabs(lambda_hessian(1))*un_sur_epsilon);
+ 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;
+
+ // 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
+ // need to find out which one of ti's is the closest to gr, to known which eigenvalue to discard...
+ std::vector<double> grad_dot_ti;
+ for (int i=0;i<3;i++)
+ grad_dot_ti.push_back(fabs(dot(gr,ti[i])));
+ std::vector<double>::iterator itbegin = grad_dot_ti.begin();
+ std::vector<double>::iterator itmax = std::max_element(itbegin,grad_dot_ti.end());
+ int grad_index = std::distance(itbegin,itmax);
+ std::vector<int> ti_index;
+ for (int i=0;i<3;i++)
+ if (i!=grad_index) ti_index.push_back(i);
+// std::cout << "gr tgr t1 t2 dots_tgr_gr_t1_t2 (" << gr(0) << "," << gr(1) << "," << gr(2) << ") (" << ti[grad_index](0) << "," << ti[grad_index](1) << "," << ti[grad_index](2) << ") (" << ti[ti_index[0]](0) << "," << ti[ti_index[0]](1) << "," << ti[ti_index[0]](2) << ") (" << ti[ti_index[1]](0) << "," << ti[ti_index[1]](1) << "," << ti[ti_index[1]](2) << ") " << fabs(dot(ti[grad_index],gr)) << " " << (dot(ti[grad_index],ti[ti_index[0]])) << " " << (dot(ti[grad_index],ti[ti_index[1]])) << std::endl;
+
+ // finally, creating the metric
+ std::vector<double> eigenvals;
+ eigenvals.push_back(std::min(std::max(eigenval_direction,metric_value_hmax),metric_value_hmin));// in gradient direction
+ eigenvals.push_back(std::min(std::max(eigenvals_curvature[ti_index[0]],metric_value_hmax),metric_value_hmin));
+ eigenvals.push_back(std::min(std::max(eigenvals_curvature[ti_index[1]],metric_value_hmax),metric_value_hmin));
+// eigenvals.push_back(std::min(std::max(metric_value_hmax,metric_value_hmax),metric_value_hmin));
+// eigenvals.push_back(std::min(std::max(metric_value_hmax,metric_value_hmax),metric_value_hmin));
+
+ metric = SMetric3(eigenvals[0],eigenvals[1],eigenvals[2],gr,ti[ti_index[0]],ti[ti_index[1]]);
+ setOfSizes[metricNumber].insert(std::make_pair(ver, std::min(std::min(1/sqrt(eigenvals[0]),1/sqrt(eigenvals[1])),1/sqrt(eigenvals[2]))));
+ }
+ else{// isotropic metric !
+ SMetric3 mymetric(metric_value_hmax);
+ metric = mymetric;
+ setOfSizes[metricNumber].insert(std::make_pair(ver, hmax));
+ }
+ _hessian[ver] = hessian;
+ setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
+ setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
+
+ it++;
+ }
+
+ if (_technique != meshMetric::EIGENDIRECTIONS ){
+
+ //See paper Hachem and Coupez:
+ //Finite element solution to handle complex heat and fluid flows in
+ //industrial furnaces using the immersed volume technique, ijnmf, 2010
+
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ H.eig(V,S);
+
+ double lambda1, lambda2, lambda3;
+ // if (dist < _E && _technique == meshMetric::FREY){
+ // fullMatrix<double> Vhess(3,3);
+ // fullVector<double> Shess(3);
+ // hessian.eig(Vhess,Shess);
+ // double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
+ // double lam1 = Shess(0);
+ // double lam2 = Shess(1);
+ // double lam3 = (_dim == 3)? Shess(2) : 1.;
+ // lambda1 = lam1;
+ // lambda2 = lam2/lambda1;
+ // lambda3 = (_dim == 3)? lam3/lambda1: 1.0;
+ // }
+ // else{
lambda1 = S(0);
lambda2 = S(1);
lambda3 = (_dim == 3)? S(2) : 1.;
- //}
+ //}
- if (_technique == meshMetric::HESSIAN || (dist < _E && _technique == meshMetric::FREY)){
- lambda1 = std::min(std::max(fabs(S(0))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
- lambda2 = std::min(std::max(fabs(S(1))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
- lambda3 = (_dim == 3) ? std::min(std::max(fabs(S(2))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin)) : 1.;
- }
+ if (_technique == meshMetric::HESSIAN || (dist < _E && _technique == meshMetric::FREY)){
+ lambda1 = std::min(std::max(fabs(S(0))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
+ lambda2 = std::min(std::max(fabs(S(1))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
+ lambda3 = (_dim == 3) ? std::min(std::max(fabs(S(2))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin)) : 1.;
+ }
- 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));
-
- SMetric3 metric(lambda1,lambda2,lambda3,t1,t2,t3);
+ 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));
- _hessian[ver] = hessian;
- _nodalMetrics[ver] = metric;
- _nodalSizes[ver] = std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3));
+ SMetric3 metric(lambda1,lambda2,lambda3,t1,t2,t3);
- _detMetric[ver] = sqrt(metric.determinant());
+ _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())));
- it++;
+ it++;
+ }
}
+
//Adapt epsilon
//-----------------
// dgDofContainer dofDet(groups, 1);
@@ -326,7 +483,7 @@ void meshMetric::computeMetric(){
// dgFunctionIntegrator integ(groups, dofDet.getFunction());
// integ.compute(res);
//fprintf("integ =%g \n", res(0,0))
-
+
//Smoothing
//-----------
//smoothMetric (sol);
@@ -335,6 +492,11 @@ void meshMetric::computeMetric(){
}
double meshMetric::operator() (double x, double y, double z, GEntity *ge) {
+ if (needMetricUpdate) intersectMetrics();
+ if (!setOfMetrics.size()){
+ std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
+ throw;
+ }
SPoint3 xyz(x,y,z), uvw;
double initialTol = MElement::getTolerance();
MElement::setTolerance(1.e-4);
@@ -354,6 +516,11 @@ 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 (!setOfMetrics.size()){
+ std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
+ throw;
+ }
metr = SMetric3(1.e-22);
SPoint3 xyz(x,y,z), uvw;
double initialTol = MElement::getTolerance();
@@ -378,25 +545,30 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
}
}
-void meshMetric::printMetric(const char* n) const{
+void meshMetric::printMetric(const char* n){
+ if (needMetricUpdate) intersectMetrics();
+ if (!setOfMetrics.size()){
+ std::cout << "meshMetric::printMetric : No metric defined ! " << std::endl;
+ throw;
+ }
FILE *f = fopen (n,"w");
fprintf(f,"View \"\"{\n");
- std::map<MVertex*,SMetric3 >::const_iterator it = _hessian.begin();
- //std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
+ //std::map<MVertex*,SMetric3 >::const_iterator it = _hessian.begin();
+ std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
//for (; it != _nodalMetrics.end(); ++it){
- for (; it != _hessian.end(); ++it){
- MVertex *v = it->first;
- SMetric3 h = it->second;
- double lapl = h(0,0)+h(1,1)+h(2,2);
- fprintf(f, "SP(%g,%g,%g){%g};\n", it->first->x(),it->first->y(),it->first->z(),lapl);
- // fprintf(f,"TP(%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
- // it->first->x(),it->first->y(),it->first->z(),
- // it->second(0,0),it->second(0,1),it->second(0,2),
- // it->second(1,0),it->second(1,1),it->second(1,2),
- // it->second(2,0),it->second(2,1),it->second(2,2)
- // );
+ for (; it != _hessian.end(); ++it){
+ MVertex *v = it->first;
+ SMetric3 h = it->second;
+ double lapl = h(0,0)+h(1,1)+h(2,2);
+ //fprintf(f, "SP(%g,%g,%g){%g};\n", it->first->x(),it->first->y(),it->first->z(),lapl);
+ fprintf(f,"TP(%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
+ it->first->x(),it->first->y(),it->first->z(),
+ h(0,0),h(0,1),h(0,2),
+ h(1,0),h(1,1),h(1,2),
+ h(2,0),h(2,1),h(2,2)
+ );
}
fprintf(f,"};\n");
@@ -432,43 +604,43 @@ double meshMetric::getLaplacian (MVertex *v) {
}*/
/*
- void meshMetric::smoothMetric (dgDofContainer *sol){
- return;
- dgGroupCollection *groups = sol->getGroups();
- std::set<MEdge,Less_Edge> edges;
- for (int i=0;i<groups->getNbElementGroups();i++){
- dgGroupOfElements *group = groups->getElementGroup(i);
- for (int j=0;j<group->getNbElements();j++){
- MElement *e = group->getElement(j);
- for (int k=0;k<e->getNumEdges();k++){
- edges.insert(e->getEdge(k));
- }
- }
- }
-
- std::set<MEdge,Less_Edge>::iterator it = edges.begin();
- double logRatio = log (CTX::instance()->mesh.smoothRatio);
- for ( ; it !=edges.end(); ++it){
- MEdge e = *it;
- std::map<MVertex*,SMetric3>::iterator it1 = _nodalMetrics.find(e.getVertex(0));
- std::map<MVertex*,SMetric3>::iterator it2 = _nodalMetrics.find(e.getVertex(1));
- SVector3 aij (e.getVertex(1)->x()-e.getVertex(0)->x(),
- e.getVertex(1)->y()-e.getVertex(0)->y(),
- e.getVertex(1)->z()-e.getVertex(0)->z());
- SMetric3 m1 = it1->second;
- double al1 = sqrt(dot(aij,m1,aij));
- SMetric3 m2 = it2->second;
- double al2 = sqrt(dot(aij,m2,aij));
- // it1->second.print("m1");
- // it2->second.print("m2");
- m1 *= 1./((1+logRatio*al1)*(1+logRatio*al1));
- m2 *= 1./((1+logRatio*al2)*(1+logRatio*al2));
- it1->second = intersection (it1->second,m2);
- it2->second = intersection (it2->second,m1);
- // it1->second.print("m1 after");
- // it2->second.print("m2 after");
- // getchar();
- }
- }
-*/
+ void meshMetric::smoothMetric (dgDofContainer *sol){
+ return;
+ dgGroupCollection *groups = sol->getGroups();
+ std::set<MEdge,Less_Edge> edges;
+ for (int i=0;i<groups->getNbElementGroups();i++){
+ dgGroupOfElements *group = groups->getElementGroup(i);
+ for (int j=0;j<group->getNbElements();j++){
+ MElement *e = group->getElement(j);
+ for (int k=0;k<e->getNumEdges();k++){
+ edges.insert(e->getEdge(k));
+ }
+ }
+ }
+
+ std::set<MEdge,Less_Edge>::iterator it = edges.begin();
+ double logRatio = log (CTX::instance()->mesh.smoothRatio);
+ for ( ; it !=edges.end(); ++it){
+ MEdge e = *it;
+ std::map<MVertex*,SMetric3>::iterator it1 = _nodalMetrics.find(e.getVertex(0));
+ std::map<MVertex*,SMetric3>::iterator it2 = _nodalMetrics.find(e.getVertex(1));
+ SVector3 aij (e.getVertex(1)->x()-e.getVertex(0)->x(),
+ e.getVertex(1)->y()-e.getVertex(0)->y(),
+ e.getVertex(1)->z()-e.getVertex(0)->z());
+ SMetric3 m1 = it1->second;
+ double al1 = sqrt(dot(aij,m1,aij));
+ SMetric3 m2 = it2->second;
+ double al2 = sqrt(dot(aij,m2,aij));
+// it1->second.print("m1");
+// it2->second.print("m2");
+m1 *= 1./((1+logRatio*al1)*(1+logRatio*al1));
+m2 *= 1./((1+logRatio*al2)*(1+logRatio*al2));
+it1->second = intersection (it1->second,m2);
+it2->second = intersection (it2->second,m1);
+// it1->second.print("m1 after");
+// it2->second.print("m2 after");
+// getchar();
+}
+}
+ */
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 9460b6cb3f..820b3a148e 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -1,6 +1,7 @@
#ifndef _MESH_METRIC_H_
#define _MESH_METRIC_H_
#include <map>
+#include <algorithm>
#include "STensor3.h"
#include "Field.h"
#include "meshGFaceOptimize.h"
@@ -13,10 +14,13 @@ class STensor3;
/**Anisotropic mesh size field based on a metric */
class meshMetric: public Field {
public:
- typedef enum {LEVELSET=1,HESSIAN=2, FREY=3} MetricComputationTechnique;
+ typedef enum {LEVELSET=1,HESSIAN=2, FREY=3, EIGENDIRECTIONS=4} MetricComputationTechnique;
private:
+ // intersect all metrics added in "setOfMetrics", preserve eigendirections of the "most anisotropic" metric
+ void intersectMetrics();
int _dim;
- double _epsilon, _E;
+ double _epsilon, _E, _E_moins, _Np;
+ bool needMetricUpdate;
meshMetric::MetricComputationTechnique _technique;
double hmin, hmax;
simpleFunction<double> *_fct;
@@ -28,16 +32,43 @@ class meshMetric: public Field {
std::map<MVertex*,double> vals;
std::map<MVertex*,SVector3> grads,dgrads[3];
- std::map<MVertex*,SMetric3> _nodalMetrics;
- std::map<MVertex*,double> _nodalSizes, _detMetric;
- std::map<MVertex*,SMetric3 > _hessian;
+ typedef std::map<MVertex*,SMetric3> nodalMetricTensor;
+ typedef std::map<MVertex*,double> nodalField;
+
+ nodalMetricTensor _nodalMetrics,_hessian;
+ nodalField _nodalSizes, _detMetric;
+
+ std::map<int,nodalMetricTensor> setOfMetrics;
+ std::map<int,nodalField> setOfSizes;
+ std::map<int,nodalField> setOfDetMetric;
+
public:
- meshMetric(std::vector<MElement*> elements, int technique,
- simpleFunction<double> *fct, std::vector<double> parameters);
- meshMetric(GModel *gm, int technique,
- simpleFunction<double> *fct, std::vector<double> parameters);
+ meshMetric(std::vector<MElement*> elements);
+ meshMetric(GModel *gm);
+
~meshMetric();
- inline SMetric3 metricAtVertex (MVertex* v) {return _nodalMetrics[v];}
+
+ // compute a new metric and add it to the set of metrics
+ // parameters[1] = lcmin (default : in global gmsh options CTX::instance()->mesh.lcMin)
+ // parameters[2] = lcmax (default : in global gmsh options CTX::instance()->mesh.lcMax)
+ // Available algorithms ("techniques"):
+ // 1: fct is a LS, metric based on Coupez technique
+ // parameters[0] = thickness of the interface (mandatory)
+ // 2: metric based on the hessian of fct
+ // parameters[0] = the final number of elements
+ // 3: fct is a LS, variant of (1) based on Frey technique (combines Coupez and curvature)
+ // parameters[0] = thickness of the interface (mandatory)
+ // parameters[3] = the required minimum number of elements to represent a circle - used for curvature-based metric (default: = 15)
+ // 4: fct is a LS, variant of (3), metric computed in LS eigendirections
+ // 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)
+ void addMetric(int technique, simpleFunction<double> *fct, std::vector<double> parameters);
+
+ inline SMetric3 metricAtVertex (MVertex* v) {
+ if (needMetricUpdate) intersectMetrics();
+ return _nodalMetrics[v];
+ }
void computeMetric() ;
void computeValues( v2t_cont adj);
void computeHessian( v2t_cont adj);
@@ -52,9 +83,13 @@ class meshMetric: public Field {
{
return "Anisotropic size field based on hessian of a given function";
}
+
+ // get metric at point(x,y,z) (previously computes intersection of metrics if not done yet)
virtual double operator() (double x, double y, double z, GEntity *ge=0) ;
virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
- void printMetric(const char* n) const;
+ void printMetric(const char* n);
+ // export pos files of fct, fct gradients (fct is the lattest fct passed to meshMetric !!) and resulting metric (intersection of all computed metrics)
+ void exportInfo(const char *fileendname);
};
#endif
diff --git a/gmshpy/gmshGeo.i b/gmshpy/gmshGeo.i
index 2fbcf2a963..a57ccb2dec 100644
--- a/gmshpy/gmshGeo.i
+++ b/gmshpy/gmshGeo.i
@@ -36,6 +36,7 @@
#include "SPoint2.h"
#include "SBoundingBox3d.h"
#include "Curvature.h"
+ #include "simpleFunction.h"
%}
namespace std {
@@ -54,6 +55,12 @@ namespace std {
}
%include "GmshConfig.h"
+%include "simpleFunction.h"
+%template(simpleFunctionDouble) simpleFunction<double>;
+%include std_vector.i
+namespace std {
+ %template(simpleFunctionDoubleVector) vector<simpleFunction<double>*, std::allocator<simpleFunction<double>*> >;
+}
%include "GModel.h"
%include "GPoint.h"
@@ -83,3 +90,4 @@ namespace std {
%include "Curvature.h"
%include "gmshLevelset.h"
+
--
GitLab