diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 9f8f4f3aa38b8dd87704f40938534e6b9d6595dc..18245e2c76efd08abac2df234fb160654065c783 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -524,7 +524,10 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
meshMetric *mm;
int ITER = 0;
+ int nbElemsOld = 0;
+ int nbElems;
while(1){
+ Msg::Info("-- adaptMesh ITER =%d ", ITER);
std::vector<MElement*> elements;
if (getDim() == 2){
@@ -544,7 +547,8 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
}
}
- if (elements.size() == 0)return -1;
+ nbElems = elements.size();
+ if (nbElems == 0)return -1;
double lcmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
double lcmax = parameters.size() >=3 ? parameters[2] : CTX::instance()->mesh.lcMax;
@@ -558,7 +562,7 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
mm = new meshMetric (elements, f, meshMetric::HESSIAN,lcmin,lcmax,0, parameters[0]);
break;
case 3 :
- mm = new meshMetric (elements, f, meshMetric::FREY,lcmin,lcmax,parameters[0], 0);
+ mm = new meshMetric (elements, f, meshMetric::FREY,lcmin,lcmax,parameters[0], 0.01);
break;
default : Msg::Error("Unknown Adaptive Strategy");return -1;
}
@@ -567,10 +571,8 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
if (getDim() == 2){
for (fiter fit = firstFace(); fit != lastFace(); ++fit){
if((*fit)->geomType() != GEntity::DiscreteSurface){
- opt_mesh_lc_from_points(0, GMSH_SET, 0);
-
meshGFaceBamg(*fit);
- laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing);
+ //laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing);
}
if(_octree) delete _octree;
_octree = 0;
@@ -584,7 +586,10 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
}
}
delete mm;
- if (++ITER > niter) break;
+ if (++ITER >= niter) break;
+ if (fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
+
+ nbElemsOld = nbElems;
}
diff --git a/Geo/GModel.h b/Geo/GModel.h
index e03481e504bea758a1a1edcf8ea6b8e51cef4968..e345f19f020a7ca7a44874291c7b751625a02e00 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -373,7 +373,7 @@ class GModel
// 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
+ // 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
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 3070425df9475a36cb096cd11c77f6b329898814..67dc0a9f1c3bf6d87eb962389ee9d1761b7235cb 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -626,30 +626,36 @@ 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);
}
-// gLevelsetPopcorn::gLevelsetPopcorn(int tag) : gLevelsetPrimitive(tag) {
-// }
-
-// double gLevelsetPopcorn::operator() (const double x, const double y, const double z) const {
-// double r0 = 0.25;
-// double r = sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5));
-// double val = r - r0;
-// for (int k = 0; k< 5; k++){
-// double xk = r0/sqrt(5)*(2.*cos(2*k*pi/5));
-// double yk = r0/sqrt(5)*(2.*sin(2*k*pi/5));
-// double zk = 1.0;
-// val += -4.*exp(((x-0.5-xk)*(x-0.5-xk)+(y-0.5-yk)*(y-0.5-yk)+(z-0.5-zk)*(z-0.5-zk))/0.01);
-// }
-// for (int k = 0; k< 5; k++){
-// xk = r0/sqrt(5)*(2*cos(2*((k-5)-1)*pi/5));
-// yk = r0/sqrt(5)*(2*sin(2*((k-5)-1)*pi/5));
-// zk = 1.0;
-// val += 4.*exp(((x-0.5-xk)*(x-0.5-xk)+(y-0.5-yk)*(y-0.5-yk)+(z-0.5-zk)*(z-0.5-zk))/0.01);
-// }
-// val += -4.*exp(((x-0.5-xk)*(x-0.5-xk)+(y-0.5-yk)*(y-0.5-yk)+(z-0.5-zk)*(z-0.5-zk))/0.01);
-// val += -4.*math.exp(((x-0.5)**2+(y-0.5)**2+(z-0.5+r0)**2)/0.01);
-// return val;
-// }
-
+gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r0, double _A, double _sigma, int tag) : gLevelsetPrimitive(tag) {
+ A = _A;
+ sigma = _sigma;
+ r0 = _r0;
+ xc = _xc;
+ yc = _yc;
+ zc = _zc;
+}
+
+double gLevelsetPopcorn::operator() (const double x, const double y, const double z) const {
+ double s2 = (sigma)*(sigma);
+ double r = sqrt((x-xc)*(x-xc)+(y-yc)*(y-yc)+(z-zc)*(z-zc));
+ //printf("z=%g zc=%g r=%g \n", z, zc, r);
+ double val = r - r0;
+ for (int k = 0; k< 5; k++){
+ double xk = r0/sqrt(5.0)*(2.*cos(2*k*M_PI/5.0));
+ double yk = r0/sqrt(5.0)*(2.*sin(2*k*M_PI/5.0));
+ double zk = r0/sqrt(5.0);
+ val -= A*exp(-((x-xc-xk)*(x-xc-xk)+(y-yc-yk)*(y-yc-yk)+(z-zc-zk)*(z-zc-zk))/s2);
+ }
+ for (int k = 5; k< 10; k++){
+ double xk = r0/sqrt(5.0)*(2.*cos((2.*(k-5.)-1.)*M_PI/5.0));
+ double yk = r0/sqrt(5.0)*(2.*sin((2.*(k-5.)-1.)*M_PI/5.0));
+ double zk = -r0/sqrt(5.0);
+ val -= A*exp(-((x-xc-xk)*(x-xc-xk)+(y-yc-yk)*(y-yc-yk)+(z-zc-zk)*(z-zc-zk))/s2);
+ }
+ val -= A*exp(-((x-xc)*(x-xc)+(y-yc)*(y-yc)+(z-zc-r0)*(z-zc-r0))/s2);
+ val -= A*exp(-((x-xc)*(x-xc)+(y-yc)*(y-yc)+(z-zc+r0)*(z-zc+r0))/s2);
+ return val;
+}
gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitive(tag) {
std::vector<std::string> expressions(1, f);
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 38ef8df4430d0c7e60318a31c85d976cc9ddcc0e..7ce2a06cfda44e94154069d8d57e861e5527d0f9 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -250,6 +250,19 @@ public:
int type() const {return QUADRIC;}
};
+class gLevelsetPopcorn: public gLevelsetPrimitive
+{
+ double A;
+ double sigma;
+ double r0;
+ double xc, yc, zc;
+public:
+ gLevelsetPopcorn(double xc, double yc, double zc, double r0, double A, double sigma, int tag=1);
+ ~gLevelsetPopcorn(){}
+ 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/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index 1b24017ea7fdfb83ebd022bc58d3f8eeed4854d5..eabd8c80fd83d6f5b3017d4b596aacb0c831e56c 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -25,6 +25,7 @@ long verbosity = 0;
#include "Mesh2.h"
Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bool);
+
static void computeMeshMetricsForBamg(GFace *gf, int numV,
Vertex2 *bamgVertices,
double *mm11, double *mm12, double *mm22)
@@ -120,30 +121,30 @@ void meshGFaceBamg(GFace *gf){
std::list<GEdge*> edges = gf->edges();
int numEdges = 0;
for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
- numEdges += (*it)->lines.size();
+ numEdges += (*it)->lines.size();
}
-
Seg *bamgBoundary = new Seg[numEdges];
int count = 0;
for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
for (unsigned int i = 0; i < (*it)->lines.size(); ++i){
int nodes [2] = {(*it)->lines[i]->getVertex(0)->getIndex(),
- (*it)->lines[i]->getVertex(1)->getIndex()};
+ (*it)->lines[i]->getVertex(1)->getIndex()};
bamgBoundary[count].init (bamgVertices, nodes, (*it)->tag());
bamgBoundary[count++].lab = count;
}
}
+ // Seg *bamgBoundary = NULL;
+ // numEdges = 0;
Mesh2 *bamgMesh = new Mesh2 ( all.size(), gf->triangles.size(), numEdges,
bamgVertices, bamgTriangles, bamgBoundary);
- //*************** refine loop here
Mesh2 *refinedBamgMesh = 0;
int iterMax = 11;
for (int k= 0; k < iterMax; k++){
- int nbVert = bamgMesh->nv;// all.size();
+ int nbVert = bamgMesh->nv;
double *mm11 = new double[nbVert];
double *mm12 = new double[nbVert];
@@ -152,7 +153,8 @@ void meshGFaceBamg(GFace *gf){
for (int i=0;i<256;i++)args[i] = -1.1e100;
args[16] = CTX::instance()->mesh.anisoMax;
args[ 7] = CTX::instance()->mesh.smoothRatio;
- computeMeshMetricsForBamg (gf, nbVert, bamgMesh->vertices , mm11,mm12,mm22); //bamgVertices
+ //args[ 21] = 90.0;//cutoffrad = 90 degree
+ computeMeshMetricsForBamg (gf, nbVert, bamgMesh->vertices , mm11,mm12,mm22);
try{
refinedBamgMesh = Bamg(bamgMesh, args, mm11, mm12, mm22, false);
@@ -174,8 +176,7 @@ void meshGFaceBamg(GFace *gf){
bamgMesh = refinedBamgMesh;
if (fabs((double)(nTnow - nT)) < 0.01 * nT) break;
}
- //*************** end for loop refine
-
+
std::map<int,MVertex*> yetAnother;
for (int i = 0; i < refinedBamgMesh->nv; i++){
Vertex2 &v = refinedBamgMesh->vertices[i];
@@ -203,6 +204,23 @@ void meshGFaceBamg(GFace *gf){
yetAnother[(*refinedBamgMesh)(v2)],
yetAnother[(*refinedBamgMesh)(v3)]));
}
+
+ //EMI PRINT TRIS
+ // FILE * fi = fopen("TRIS.pos","w");
+ // fprintf(fi, "View \"\"{\n");
+ // for( int i =0; i< gf->triangles.size(); i++){
+ // MTriangle *t = gf->triangles[i];
+ // fprintf(fi, "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,"
+ // "%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
+ // t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+ // t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+ // t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
+ // 1., 1., 1.);
+ // }
+ // fprintf(fi,"};\n");
+ // fclose(fi);
+ //END EMI PRINT TRIS
+
if (refinedBamgMesh) delete refinedBamgMesh;
}
diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp
index b664d1fb351d145ec5ba07683d043b96790a5fbb..4d9533e4de520e347b19fe120fe958fa182ae83e 100644
--- a/Mesh/meshGRegionMMG3D.cpp
+++ b/Mesh/meshGRegionMMG3D.cpp
@@ -262,7 +262,10 @@ void refineMeshMMG(GRegion *gr)
std::map<int,MVertex*> mmg2gmsh;
gmsh2MMG (gr, mmg, sol,mmg2gmsh);
- for (int ITER=0;ITER<2;ITER++){
+ int iterMax = 11;
+ for (int ITER=0;ITER<iterMax;ITER++){
+ int nT = mmg->ne;
+
int verb_mmg = (Msg::GetVerbosity() > 9) ? -1 : 0;
int opt[9] = {1,0,64,0,0,0, verb_mmg , 0,0};
Msg::Debug("-------- GMSH LAUNCHES MMG3D ---------------");
@@ -270,9 +273,13 @@ void refineMeshMMG(GRegion *gr)
Msg::Debug("-------- MG3D TERMINATED -------------------");
// Here we should interact with BGM
updateSizes(gr,mmg, sol);
+
+ int nTnow = mmg->ne;
+ if (fabs((double)(nTnow - nT)) < 0.05 * nT) break;
}
- //char test[] = "test.mesh";
- //MMG_saveMesh(mmg, test);
+
+ char test[] = "test.mesh";
+ MMG_saveMesh(mmg, test);
gr->deleteVertexArrays();
for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i];
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 830813f6fcbe7f4d9725b4b11db86752c615e45f..4282817ec506c1f7365fea4690089213d60a4fcb 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -22,14 +22,12 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l
lt.insert(lt.begin(),stencil.begin(),stencil.end());
}
-meshMetric::meshMetric ( std::vector<MElement*> &elements, simpleFunction<double> *fct, meshMetric::MetricComputationTechnique technique,
- double lcmin, double lcmax,
- double E, double eps)
+meshMetric::meshMetric ( std::vector<MElement*> &elements, simpleFunction<double> *fct, meshMetric::MetricComputationTechnique technique, double lcmin, double lcmax, double E, double eps)
: _technique(technique), _epsilon(eps), _E(E), _fct(fct), hmax(lcmax),hmin(lcmin)
{
_dim = elements[0]->getDim();
computeMetric(elements);
- printMetric("toto.pos");
+ //printMetric("toto.pos");
}
@@ -110,7 +108,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
v2t_cont adj;
buildVertexToElement (e,adj);
- printf("%d elements are considered in the metric \n",e.size());
+ //printf("%d elements are considered in the metric \n",e.size());
computeValues(adj);
computeHessian(adj);
@@ -118,6 +116,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
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];
@@ -134,48 +133,60 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
if (_technique == meshMetric::HESSIAN){
H.setMat(hessian);
}
- else if (_technique == meshMetric::FREY){
- SVector3 gr = grads[ver];
- fullMatrix<double> hfrey(hessian);
- hfrey(0,0) += gr(0)*gr(0)/(hmin*hmin);
- hfrey(1,1) += gr(1)*gr(1)/(hmin*hmin);
- hfrey(2,2) += gr(2)*gr(2)/(hmin*hmin);
- hfrey(1,0) += gr(1)*gr(0)/(hmin*hmin);
- hfrey(2,0) += gr(2)*gr(0)/(hmin*hmin);
- hfrey(2,1) += gr(2)*gr(1)/(hmin*hmin);
- hfrey(0,1) = hfrey(1,0) ;
- hfrey(0,2) = hfrey(2,0);
- hfrey(1,2) = hfrey(2,1);
+ //See paper Ducrot and Frey:
+ //Anisotropic levelset adaptation for accurate interface capturing,
+ //ijnmf, 2010
+ else if (_technique == meshMetric::FREY ){
+ SVector3 gr = grads[ver];
+ fullMatrix<double> hfrey(3,3);
+ hfrey(0,0) = 1./(hmax*hmax);
+ hfrey(1,1) = 1./(hmax*hmax);
+ hfrey(2,2) = 1./(hmax*hmax);
+ double divEps = 1./0.01;
+ 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);
+ hfrey(0,0) += C*gr(0)*gr(0)/(norm) + hessian(0,0)*divEps; //metric intersection ???
+ hfrey(1,1) += C*gr(1)*gr(1)/(norm) + hessian(1,1)*divEps;
+ hfrey(2,2) += C*gr(2)*gr(2)/(norm) + hessian(2,2)*divEps;
+ hfrey(1,0) = hfrey(0,1) = C*gr(1)*gr(0)/(norm) + hessian(1,0)*divEps;
+ hfrey(2,0) = hfrey(0,2) = C*gr(2)*gr(0)/(norm) + hessian(2,0)*divEps;
+ hfrey(2,1) = hfrey(1,2) = C*gr(2)*gr(1)/(norm) + hessian(2,1)*divEps;
+ }
H.setMat(hfrey);
- }
- else if (_technique == meshMetric::LEVELSET){
- double dist = fabs(vals[ver]);
+ }
+ //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
+ else if (_technique == meshMetric::LEVELSET ){
SVector3 gr = grads[ver];
- fullMatrix<double> hcoupez(3,3);
- hcoupez(0,0) = 1./(hmax*hmax);
- hcoupez(1,1) = 1./(hmax*hmax);
- hcoupez(2,2) = 1./(hmax*hmax);
- double norm = sqrt(gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2));
+ fullMatrix<double> hlevelset(3,3);
+ hlevelset(0,0) = 1./(hmax*hmax);
+ hlevelset(1,1) = 1./(hmax*hmax);
+ hlevelset(2,2) = 1./(hmax*hmax);
+ double norm = gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2);
if (dist < _E && norm != 0.0){
- double C = 1./(hmin*hmin) -1./(hmax*hmax);
- hcoupez(0,0) += C*gr(0)*gr(0)/norm ;
- hcoupez(1,1) += C*gr(1)*gr(1)/norm ;
- hcoupez(2,2) += C*gr(2)*gr(2)/norm ;
- hcoupez(1,0) = hcoupez(0,1) = C*gr(1)*gr(0)/norm;
- hcoupez(2,0) = hcoupez(0,2) = C*gr(2)*gr(0)/norm;
- hcoupez(2,1) = hcoupez(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);
- double C = 1./(hmid*hmid) -1./(hmax*hmax);
- hcoupez(0,0) += C*gr(0)*gr(0)/norm ;
- hcoupez(1,1) += C*gr(1)*gr(1)/norm ;
- hcoupez(2,2) += C*gr(2)*gr(2)/norm ;
- hcoupez(1,0) = hcoupez(0,1) = C*gr(1)*gr(0)/norm;
- hcoupez(2,0) = hcoupez(0,2) = C*gr(2)*gr(0)/norm;
- hcoupez(2,1) = hcoupez(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 ;
}
- H.setMat(hcoupez);
+ // else if (dist > _E && dist < 2.*_E && norm != 0.0){
+ // double hmid = (hmax-hmin)/(1.*_E)*dist+(2.*hmin-1.*hmax);
+ // double C = 1./(hmid*hmid) -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 ;
+ // }
+ H.setMat(hlevelset);
}
fullMatrix<double> V(3,3);
@@ -196,6 +207,22 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
SVector3 t3 (V(0,2),V(1,2),V(2,2));
SMetric3 metric(lambda1,lambda2,lambda3,t1,t2,t3);
+
+ // if (_technique == meshMetric::FREY && dist < _E) {
+ // SMetric3 Hess;
+ // fullMatrix<double> hessianFrey(hessian);
+ // hessianFrey.scale(1./0.1);
+ // Hess.setMat(hessianFrey);
+ // fullMatrix<double> Vh(3,3);
+ // fullVector<double> Sh(3);
+ // Hess.eig(Vh,Sh);
+ // lambda1 = std::min(std::max(fabs(Sh(0))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
+ // lambda2 = std::min(std::max(fabs(Sh(1))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
+ // lambda3 = (_dim == 3) ? std::min(std::max(fabs(Sh(2))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin)) : 1.;
+ // SMetric3 inter = intersection(metric, Hess);
+ // metric = inter;
+ // }
+
_hessian[ver] = hessian;
_nodalMetrics[ver] = metric;
_nodalSizes[ver] = std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3));
@@ -230,12 +257,11 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
//smoothMetric (sol);
//curvatureContributionToMetric();
+
putOnNewView();
}
-
-
double meshMetric::operator() (double x, double y, double z, GEntity *ge){
GModel *gm = _nodalMetrics.begin()->first->onWhat()->model();
SPoint3 xyz(x,y,z),uvw;
@@ -319,7 +345,6 @@ double meshMetric::getLaplacian (MVertex *v){
return laplace;
}
-
/*void dgMeshMetric::curvatureContributionToMetric (){
std::map<MVertex*,SMetric3>::iterator it = _nodalMetrics.begin();
std::map<MVertex*,double>::iterator its = _nodalSizes.begin();
diff --git a/contrib/bamg/bamg-gmsh.cpp b/contrib/bamg/bamg-gmsh.cpp
index 72da8563a938c911e62ef48d1bbb79c9ddd56714..999226ceed6bb579042f7fd45f0e3fba718f214e 100644
--- a/contrib/bamg/bamg-gmsh.cpp
+++ b/contrib/bamg/bamg-gmsh.cpp
@@ -414,10 +414,11 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo
oTh = msh2bamg(*Thh,cutoffradian,nbdfv,ndfv,nbdfe,ndfe,reqedges,reqedges.N());
}
else */
- oTh = msh2bamg(*Thh,cutoffradian,reqedges,reqedges.N());
+
+ oTh = msh2bamg(*Thh,cutoffradian,reqedges,reqedges.N());
Triangles &Th(*oTh);
- // printf("COUCOUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU\n");
+ // printf("COUCOUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU\n");
// Th.Write("toto.mesh",bamg::Triangles::AutoMesh);
@@ -537,7 +538,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo
Triangles* nTh = 0;
- // Th.Write("toto.msh",bamg::Triangles::AutoMesh);
+ // Th.Write("toto.msh",bamg::Triangles::AutoMesh);
if (initialMesh){
nTh= new Triangles(nbsx,Th.Gh);
@@ -546,7 +547,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo
nTh= new Triangles(nbsx,Th,KeepBackVertices); // Adaption is here
}
- // nTh->Write("tata.mesh",bamg::Triangles::AutoMesh);
+ //nTh->Write("tata.mesh",bamg::Triangles::AutoMesh);
if (split)
@@ -568,6 +569,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo
for (int iv=0;iv < Th.nbv;iv++)
Th[iv].m = M;
+ //nTh->Write("bamg.mesh",bamg::Triangles::AutoMesh);
Mesh2 * g= bamg2msh(nTh,true);
delete nTh;
diff --git a/contrib/mmg3d/build/sources/mmg3dlib.c b/contrib/mmg3d/build/sources/mmg3dlib.c
index 05615c533fec11d82edf903cc6657a9236bc7879..724671a0d55d65544bd0540c861dfc56896729fd 100644
--- a/contrib/mmg3d/build/sources/mmg3dlib.c
+++ b/contrib/mmg3d/build/sources/mmg3dlib.c
@@ -256,8 +256,8 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) {
short imprim;
int k,iadr,i,jj,kk,ii,im;
double lambda[3],v[3][3],*mold,*m;
- fprintf(stdout," %s \n", M_STR);
- fprintf(stdout," -- START MMG3d \n") ;
+ //fprintf(stdout," %s \n", M_STR);
+ fprintf(stdout," -- START MMG3d (%d ELEMS)\n", mesh->ne) ;
if (opt[6] < 0) fprintf(stdout," -- MMG3d, Release %s (%s) \n",M_VER,M_REL);
if (opt[6] < -10) fprintf(stdout," Copyright (c) LJLL/IMB, 2010\n");
if (opt[6] < -10) fprintf(stdout," %s\n",COMPIL);
@@ -469,7 +469,7 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) {
MMG_prilen(mesh,sol);
MMG_ratio(mesh,sol,NULL);
}
- fprintf(stdout," -- END MMG3D \n");
+ fprintf(stdout," -- END MMG3D (%d ELEMS)\n", mesh->ne);
if ( alert )
fprintf(stdout,"\n ## WARNING: INCOMPLETE MESH %d , %d\n",
mesh->np,mesh->ne);
@@ -487,7 +487,6 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) {
fprintf(stdout," NUMBER OF GIVEN ELEMENTS %8d\n",mesh->nefixe);
fprintf(stdout," TOTAL NUMBER OF VERTICES %8d\n",mesh->np);
fprintf(stdout," TOTAL NUMBER OF TRIANGLES %8d\n",mesh->nt);
- fprintf(stdout," TOTAL NUMBER OF ELEMENTS %8d\n",mesh->ne);
}
if ( MMG_imprim ) fprintf(stdout," -- SAVING COMPLETED\n");