diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e2752279f5d69ab43c507ff5d019ff9275616008..9d1bca04050127bb482df659fdf0da22e57f702a 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 2196de1c5db60f61b3f05e909350fc8aa9f0c432..c856ae18d315562deb5429f131ebeaa9f79aab82 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 7436bc7a510317408082580da8388b3a4d9d37fa..0bbe5ee232c39c2a1764f711391a77b8d92eb117 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 6673cc5908b506b7c94f31e36653d862c1f72bce..c99db1b89c3870ec043bfc528fc13748fc087a4b 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 bed67951e91762c7c9cf88c62bb919783e1f4bba..f0515ed706f0e1a1807b4e7cf1b76ae1a6387bbf 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 590e25d40b5f5d0b16a5408775119712ddf940f1..2712cc21631e590e1b7f789083e73733d84fb352 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 eb690c3f925e5bc644a0bf627d14752190bc2db7..766628d747dfcc5049a6c570e6c1328a5604a8d1 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 9460b6cb3fdeb8e7423d0e40e1010ff7cdf35911..820b3a148edfe2f21b257e7e75a928bf389a3459 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 2fbcf2a9639e04ddc91b7e8cdf8c716998dcd263..a57ccb2dec132ebe00ccf7c4fe1d8aa5a088a25f 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"
+