Skip to content
Snippets Groups Projects
Commit 92bb2f31 authored by Emilie Marchandise's avatar Emilie Marchandise
Browse files

work on adaptMesh

parent 099a1bf4
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/GModel.cpp
+ 11
− 6
View file @ 92bb2f31
...@@ -524,7 +524,10 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub ...@@ -524,7 +524,10 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
meshMetric *mm; meshMetric *mm;
int ITER = 0; int ITER = 0;
int nbElemsOld = 0;
int nbElems;
while(1){ while(1){
Msg::Info("-- adaptMesh ITER =%d ", ITER);
std::vector<MElement*> elements; std::vector<MElement*> elements;
if (getDim() == 2){ if (getDim() == 2){
...@@ -544,7 +547,8 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub ...@@ -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 lcmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
double lcmax = parameters.size() >=3 ? parameters[2] : CTX::instance()->mesh.lcMax; 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 ...@@ -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]); mm = new meshMetric (elements, f, meshMetric::HESSIAN,lcmin,lcmax,0, parameters[0]);
break; break;
case 3 : 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; break;
default : Msg::Error("Unknown Adaptive Strategy");return -1; default : Msg::Error("Unknown Adaptive Strategy");return -1;
} }
...@@ -567,10 +571,8 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub ...@@ -567,10 +571,8 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
if (getDim() == 2){ if (getDim() == 2){
for (fiter fit = firstFace(); fit != lastFace(); ++fit){ for (fiter fit = firstFace(); fit != lastFace(); ++fit){
if((*fit)->geomType() != GEntity::DiscreteSurface){ if((*fit)->geomType() != GEntity::DiscreteSurface){
opt_mesh_lc_from_points(0, GMSH_SET, 0);
meshGFaceBamg(*fit); meshGFaceBamg(*fit);
laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing); //laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing);
} }
if(_octree) delete _octree; if(_octree) delete _octree;
_octree = 0; _octree = 0;
...@@ -584,7 +586,10 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub ...@@ -584,7 +586,10 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
} }
} }
delete mm; delete mm;
if (++ITER > niter) break; if (++ITER >= niter) break;
if (fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
nbElemsOld = nbElems;
} }
......
Geo/GModel.h
+ 1
− 1
View file @ 92bb2f31
...@@ -373,7 +373,7 @@ class GModel ...@@ -373,7 +373,7 @@ class GModel
// parameters[0] = thickness of the interface (mandatory) // parameters[0] = thickness of the interface (mandatory)
// 2) Assume that the function is a physical quantity -> adapt using the Hessain (technique = 2) // 2) Assume that the function is a physical quantity -> adapt using the Hessain (technique = 2)
// parameters[0] = N, the final number of elements // 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) // parameters[0] = thickness of the interface (mandatory)
// The algorithm first generate a mesh if no one is available // The algorithm first generate a mesh if no one is available
......
Geo/gmshLevelset.cpp
+ 30
− 24
View file @ 92bb2f31
...@@ -626,30 +626,36 @@ double gLevelsetQuadric::operator()(const double x, const double y, const double ...@@ -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); + 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) { gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r0, double _A, double _sigma, int tag) : gLevelsetPrimitive(tag) {
// } A = _A;
sigma = _sigma;
// double gLevelsetPopcorn::operator() (const double x, const double y, const double z) const { r0 = _r0;
// double r0 = 0.25; xc = _xc;
// double r = sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5)); yc = _yc;
// double val = r - r0; zc = _zc;
// 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 gLevelsetPopcorn::operator() (const double x, const double y, const double z) const {
// double zk = 1.0; double s2 = (sigma)*(sigma);
// 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); 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);
// for (int k = 0; k< 5; k++){ double val = r - r0;
// xk = r0/sqrt(5)*(2*cos(2*((k-5)-1)*pi/5)); for (int k = 0; k< 5; k++){
// yk = r0/sqrt(5)*(2*sin(2*((k-5)-1)*pi/5)); double xk = r0/sqrt(5.0)*(2.*cos(2*k*M_PI/5.0));
// zk = 1.0; double yk = r0/sqrt(5.0)*(2.*sin(2*k*M_PI/5.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); 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 += -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); for (int k = 5; k< 10; k++){
// return val; 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) { gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitive(tag) {
std::vector<std::string> expressions(1, f); std::vector<std::string> expressions(1, f);
......
Geo/gmshLevelset.h
+ 13
− 0
View file @ 92bb2f31
...@@ -250,6 +250,19 @@ public: ...@@ -250,6 +250,19 @@ public:
int type() const {return QUADRIC;} 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 class gLevelsetMathEval: public gLevelsetPrimitive
{ {
mathEvaluator *_expr; mathEvaluator *_expr;
......
Mesh/meshGFaceBamg.cpp
+ 26
− 8
View file @ 92bb2f31
...@@ -25,6 +25,7 @@ long verbosity = 0; ...@@ -25,6 +25,7 @@ long verbosity = 0;
#include "Mesh2.h" #include "Mesh2.h"
Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bool); Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bool);
static void computeMeshMetricsForBamg(GFace *gf, int numV, static void computeMeshMetricsForBamg(GFace *gf, int numV,
Vertex2 *bamgVertices, Vertex2 *bamgVertices,
double *mm11, double *mm12, double *mm22) double *mm11, double *mm12, double *mm22)
...@@ -120,30 +121,30 @@ void meshGFaceBamg(GFace *gf){ ...@@ -120,30 +121,30 @@ void meshGFaceBamg(GFace *gf){
std::list<GEdge*> edges = gf->edges(); std::list<GEdge*> edges = gf->edges();
int numEdges = 0; int numEdges = 0;
for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){ 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]; Seg *bamgBoundary = new Seg[numEdges];
int count = 0; int count = 0;
for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){ for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
for (unsigned int i = 0; i < (*it)->lines.size(); ++i){ for (unsigned int i = 0; i < (*it)->lines.size(); ++i){
int nodes [2] = {(*it)->lines[i]->getVertex(0)->getIndex(), 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].init (bamgVertices, nodes, (*it)->tag());
bamgBoundary[count++].lab = count; bamgBoundary[count++].lab = count;
} }
} }
// Seg *bamgBoundary = NULL;
// numEdges = 0;
Mesh2 *bamgMesh = new Mesh2 ( all.size(), gf->triangles.size(), numEdges, Mesh2 *bamgMesh = new Mesh2 ( all.size(), gf->triangles.size(), numEdges,
bamgVertices, bamgTriangles, bamgBoundary); bamgVertices, bamgTriangles, bamgBoundary);
//*************** refine loop here
Mesh2 *refinedBamgMesh = 0; Mesh2 *refinedBamgMesh = 0;
int iterMax = 11; int iterMax = 11;
for (int k= 0; k < iterMax; k++){ for (int k= 0; k < iterMax; k++){
int nbVert = bamgMesh->nv;// all.size(); int nbVert = bamgMesh->nv;
double *mm11 = new double[nbVert]; double *mm11 = new double[nbVert];
double *mm12 = new double[nbVert]; double *mm12 = new double[nbVert];
...@@ -152,7 +153,8 @@ void meshGFaceBamg(GFace *gf){ ...@@ -152,7 +153,8 @@ void meshGFaceBamg(GFace *gf){
for (int i=0;i<256;i++)args[i] = -1.1e100; for (int i=0;i<256;i++)args[i] = -1.1e100;
args[16] = CTX::instance()->mesh.anisoMax; args[16] = CTX::instance()->mesh.anisoMax;
args[ 7] = CTX::instance()->mesh.smoothRatio; 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{ try{
refinedBamgMesh = Bamg(bamgMesh, args, mm11, mm12, mm22, false); refinedBamgMesh = Bamg(bamgMesh, args, mm11, mm12, mm22, false);
...@@ -174,8 +176,7 @@ void meshGFaceBamg(GFace *gf){ ...@@ -174,8 +176,7 @@ void meshGFaceBamg(GFace *gf){
bamgMesh = refinedBamgMesh; bamgMesh = refinedBamgMesh;
if (fabs((double)(nTnow - nT)) < 0.01 * nT) break; if (fabs((double)(nTnow - nT)) < 0.01 * nT) break;
} }
//*************** end for loop refine
std::map<int,MVertex*> yetAnother; std::map<int,MVertex*> yetAnother;
for (int i = 0; i < refinedBamgMesh->nv; i++){ for (int i = 0; i < refinedBamgMesh->nv; i++){
Vertex2 &v = refinedBamgMesh->vertices[i]; Vertex2 &v = refinedBamgMesh->vertices[i];
...@@ -203,6 +204,23 @@ void meshGFaceBamg(GFace *gf){ ...@@ -203,6 +204,23 @@ void meshGFaceBamg(GFace *gf){
yetAnother[(*refinedBamgMesh)(v2)], yetAnother[(*refinedBamgMesh)(v2)],
yetAnother[(*refinedBamgMesh)(v3)])); 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; if (refinedBamgMesh) delete refinedBamgMesh;
} }
......
Mesh/meshGRegionMMG3D.cpp
+ 10
− 3
View file @ 92bb2f31
...@@ -262,7 +262,10 @@ void refineMeshMMG(GRegion *gr) ...@@ -262,7 +262,10 @@ void refineMeshMMG(GRegion *gr)
std::map<int,MVertex*> mmg2gmsh; std::map<int,MVertex*> mmg2gmsh;
gmsh2MMG (gr, mmg, sol,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 verb_mmg = (Msg::GetVerbosity() > 9) ? -1 : 0;
int opt[9] = {1,0,64,0,0,0, verb_mmg , 0,0}; int opt[9] = {1,0,64,0,0,0, verb_mmg , 0,0};
Msg::Debug("-------- GMSH LAUNCHES MMG3D ---------------"); Msg::Debug("-------- GMSH LAUNCHES MMG3D ---------------");
...@@ -270,9 +273,13 @@ void refineMeshMMG(GRegion *gr) ...@@ -270,9 +273,13 @@ void refineMeshMMG(GRegion *gr)
Msg::Debug("-------- MG3D TERMINATED -------------------"); Msg::Debug("-------- MG3D TERMINATED -------------------");
// Here we should interact with BGM // Here we should interact with BGM
updateSizes(gr,mmg, sol); 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(); gr->deleteVertexArrays();
for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i]; for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i];
......
Mesh/meshMetric.cpp
+ 71
− 46
View file @ 92bb2f31
...@@ -22,14 +22,12 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l ...@@ -22,14 +22,12 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l
lt.insert(lt.begin(),stencil.begin(),stencil.end()); lt.insert(lt.begin(),stencil.begin(),stencil.end());
} }
meshMetric::meshMetric ( std::vector<MElement*> &elements, simpleFunction<double> *fct, meshMetric::MetricComputationTechnique technique, meshMetric::meshMetric ( std::vector<MElement*> &elements, simpleFunction<double> *fct, meshMetric::MetricComputationTechnique technique, double lcmin, double lcmax, double E, double eps)
double lcmin, double lcmax,
double E, double eps)
: _technique(technique), _epsilon(eps), _E(E), _fct(fct), hmax(lcmax),hmin(lcmin) : _technique(technique), _epsilon(eps), _E(E), _fct(fct), hmax(lcmax),hmin(lcmin)
{ {
_dim = elements[0]->getDim(); _dim = elements[0]->getDim();
computeMetric(elements); computeMetric(elements);
printMetric("toto.pos"); //printMetric("toto.pos");
} }
...@@ -110,7 +108,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){ ...@@ -110,7 +108,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
v2t_cont adj; v2t_cont adj;
buildVertexToElement (e,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); computeValues(adj);
computeHessian(adj); computeHessian(adj);
...@@ -118,6 +116,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){ ...@@ -118,6 +116,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
v2t_cont :: iterator it = adj.begin(); v2t_cont :: iterator it = adj.begin();
while (it != adj.end()) { while (it != adj.end()) {
MVertex *ver = it->first; MVertex *ver = it->first;
double dist = fabs(vals[ver]);
SMetric3 H; SMetric3 H;
SVector3 gradudx = dgrads[0][ver]; SVector3 gradudx = dgrads[0][ver];
...@@ -134,48 +133,60 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){ ...@@ -134,48 +133,60 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
if (_technique == meshMetric::HESSIAN){ if (_technique == meshMetric::HESSIAN){
H.setMat(hessian); H.setMat(hessian);
} }
else if (_technique == meshMetric::FREY){ //See paper Ducrot and Frey:
SVector3 gr = grads[ver]; //Anisotropic levelset adaptation for accurate interface capturing,
fullMatrix<double> hfrey(hessian); //ijnmf, 2010
hfrey(0,0) += gr(0)*gr(0)/(hmin*hmin); else if (_technique == meshMetric::FREY ){
hfrey(1,1) += gr(1)*gr(1)/(hmin*hmin); SVector3 gr = grads[ver];
hfrey(2,2) += gr(2)*gr(2)/(hmin*hmin); fullMatrix<double> hfrey(3,3);
hfrey(1,0) += gr(1)*gr(0)/(hmin*hmin); hfrey(0,0) = 1./(hmax*hmax);
hfrey(2,0) += gr(2)*gr(0)/(hmin*hmin); hfrey(1,1) = 1./(hmax*hmax);
hfrey(2,1) += gr(2)*gr(1)/(hmin*hmin); hfrey(2,2) = 1./(hmax*hmax);
hfrey(0,1) = hfrey(1,0) ; double divEps = 1./0.01;
hfrey(0,2) = hfrey(2,0); double norm = gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2);
hfrey(1,2) = hfrey(2,1); 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); H.setMat(hfrey);
} }
else if (_technique == meshMetric::LEVELSET){ //See paper Hachem and Coupez:
double dist = fabs(vals[ver]); //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]; SVector3 gr = grads[ver];
fullMatrix<double> hcoupez(3,3); fullMatrix<double> hlevelset(3,3);
hcoupez(0,0) = 1./(hmax*hmax); hlevelset(0,0) = 1./(hmax*hmax);
hcoupez(1,1) = 1./(hmax*hmax); hlevelset(1,1) = 1./(hmax*hmax);
hcoupez(2,2) = 1./(hmax*hmax); hlevelset(2,2) = 1./(hmax*hmax);
double norm = sqrt(gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2)); double norm = gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2);
if (dist < _E && norm != 0.0){ if (dist < _E && norm != 0.0){
double C = 1./(hmin*hmin) -1./(hmax*hmax); double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
hcoupez(0,0) += C*gr(0)*gr(0)/norm ; double C = 1./(h*h) -1./(hmax*hmax);
hcoupez(1,1) += C*gr(1)*gr(1)/norm ; hlevelset(0,0) += C*gr(0)*gr(0)/norm ;
hcoupez(2,2) += C*gr(2)*gr(2)/norm ; hlevelset(1,1) += C*gr(1)*gr(1)/norm ;
hcoupez(1,0) = hcoupez(0,1) = C*gr(1)*gr(0)/norm; hlevelset(2,2) += C*gr(2)*gr(2)/norm ;
hcoupez(2,0) = hcoupez(0,2) = C*gr(2)*gr(0)/norm; hlevelset(1,0) = hlevelset(0,1) = C*gr(1)*gr(0)/norm ;
hcoupez(2,1) = hcoupez(1,2) = C*gr(2)*gr(1)/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);
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;
} }
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); fullMatrix<double> V(3,3);
...@@ -196,6 +207,22 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){ ...@@ -196,6 +207,22 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
SVector3 t3 (V(0,2),V(1,2),V(2,2)); SVector3 t3 (V(0,2),V(1,2),V(2,2));
SMetric3 metric(lambda1,lambda2,lambda3,t1,t2,t3); 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; _hessian[ver] = hessian;
_nodalMetrics[ver] = metric; _nodalMetrics[ver] = metric;
_nodalSizes[ver] = std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3)); _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){ ...@@ -230,12 +257,11 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
//smoothMetric (sol); //smoothMetric (sol);
//curvatureContributionToMetric(); //curvatureContributionToMetric();
putOnNewView(); putOnNewView();
} }
double meshMetric::operator() (double x, double y, double z, GEntity *ge){ double meshMetric::operator() (double x, double y, double z, GEntity *ge){
GModel *gm = _nodalMetrics.begin()->first->onWhat()->model(); GModel *gm = _nodalMetrics.begin()->first->onWhat()->model();
SPoint3 xyz(x,y,z),uvw; SPoint3 xyz(x,y,z),uvw;
...@@ -319,7 +345,6 @@ double meshMetric::getLaplacian (MVertex *v){ ...@@ -319,7 +345,6 @@ double meshMetric::getLaplacian (MVertex *v){
return laplace; return laplace;
} }
/*void dgMeshMetric::curvatureContributionToMetric (){ /*void dgMeshMetric::curvatureContributionToMetric (){
std::map<MVertex*,SMetric3>::iterator it = _nodalMetrics.begin(); std::map<MVertex*,SMetric3>::iterator it = _nodalMetrics.begin();
std::map<MVertex*,double>::iterator its = _nodalSizes.begin(); std::map<MVertex*,double>::iterator its = _nodalSizes.begin();
......
contrib/bamg/bamg-gmsh.cpp
+ 6
− 4
View file @ 92bb2f31
...@@ -414,10 +414,11 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo ...@@ -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()); oTh = msh2bamg(*Thh,cutoffradian,nbdfv,ndfv,nbdfe,ndfe,reqedges,reqedges.N());
} }
else */ else */
oTh = msh2bamg(*Thh,cutoffradian,reqedges,reqedges.N());
oTh = msh2bamg(*Thh,cutoffradian,reqedges,reqedges.N());
Triangles &Th(*oTh); Triangles &Th(*oTh);
// printf("COUCOUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU\n"); // printf("COUCOUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU\n");
// Th.Write("toto.mesh",bamg::Triangles::AutoMesh); // Th.Write("toto.mesh",bamg::Triangles::AutoMesh);
...@@ -537,7 +538,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo ...@@ -537,7 +538,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo
Triangles* nTh = 0; Triangles* nTh = 0;
// Th.Write("toto.msh",bamg::Triangles::AutoMesh); // Th.Write("toto.msh",bamg::Triangles::AutoMesh);
if (initialMesh){ if (initialMesh){
nTh= new Triangles(nbsx,Th.Gh); nTh= new Triangles(nbsx,Th.Gh);
...@@ -546,7 +547,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo ...@@ -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= new Triangles(nbsx,Th,KeepBackVertices); // Adaption is here
} }
// nTh->Write("tata.mesh",bamg::Triangles::AutoMesh); //nTh->Write("tata.mesh",bamg::Triangles::AutoMesh);
if (split) if (split)
...@@ -568,6 +569,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo ...@@ -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++) for (int iv=0;iv < Th.nbv;iv++)
Th[iv].m = M; Th[iv].m = M;
//nTh->Write("bamg.mesh",bamg::Triangles::AutoMesh);
Mesh2 * g= bamg2msh(nTh,true); Mesh2 * g= bamg2msh(nTh,true);
delete nTh; delete nTh;
......
contrib/mmg3d/build/sources/mmg3dlib.c
+ 3
− 4
View file @ 92bb2f31
...@@ -256,8 +256,8 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) { ...@@ -256,8 +256,8 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) {
short imprim; short imprim;
int k,iadr,i,jj,kk,ii,im; int k,iadr,i,jj,kk,ii,im;
double lambda[3],v[3][3],*mold,*m; double lambda[3],v[3][3],*mold,*m;
fprintf(stdout," %s \n", M_STR); //fprintf(stdout," %s \n", M_STR);
fprintf(stdout," -- START MMG3d \n") ; 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] < 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," Copyright (c) LJLL/IMB, 2010\n");
if (opt[6] < -10) fprintf(stdout," %s\n",COMPIL); 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) { ...@@ -469,7 +469,7 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) {
MMG_prilen(mesh,sol); MMG_prilen(mesh,sol);
MMG_ratio(mesh,sol,NULL); MMG_ratio(mesh,sol,NULL);
} }
fprintf(stdout," -- END MMG3D \n"); fprintf(stdout," -- END MMG3D (%d ELEMS)\n", mesh->ne);
if ( alert ) if ( alert )
fprintf(stdout,"\n ## WARNING: INCOMPLETE MESH %d , %d\n", fprintf(stdout,"\n ## WARNING: INCOMPLETE MESH %d , %d\n",
mesh->np,mesh->ne); mesh->np,mesh->ne);
...@@ -487,7 +487,6 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) { ...@@ -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," NUMBER OF GIVEN ELEMENTS %8d\n",mesh->nefixe);
fprintf(stdout," TOTAL NUMBER OF VERTICES %8d\n",mesh->np); 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 TRIANGLES %8d\n",mesh->nt);
fprintf(stdout," TOTAL NUMBER OF ELEMENTS %8d\n",mesh->ne);
} }
if ( MMG_imprim ) fprintf(stdout," -- SAVING COMPLETED\n"); if ( MMG_imprim ) fprintf(stdout," -- SAVING COMPLETED\n");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment