diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 9c4e4f92a2acd7f543e8bd0a3c4847d18b1ab85e..b4faa37a2157880aebf3c1e5c82816389f94a727 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -541,7 +541,8 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
GenerateMesh(this, getDim());
nbElems = getNumMeshElements();
- //if (fields) fields->deleteField(1);
+ writeMSH("meshAdapt.msh");
+
fields->reset();
if (++ITER >= niter) break;
if (ITER > 5 && fabs((double)(nbElems - nbElemsOld)) < 0.005 * nbElemsOld) break;
@@ -549,6 +550,7 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
int id = fields->newId();
(*fields)[id] = new meshMetric(this, technique, f, parameters);;
fields->background_field = id;
+
std::for_each(firstEdge(), lastEdge(), deMeshGEdge());
std::for_each(firstFace(), lastFace(), deMeshGFace());
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 1bddd8249a1fb89d8896283e75ea907a4d1be233..ba6203e8c33d61381003684ee66bdfc06788d037 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1722,10 +1722,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
else if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY ||
CTX::instance()->mesh.algo2d == ALGO_2D_AUTO)
bowyerWatson(gf);
- else {
- printf("in bamg \n");
+ else
meshGFaceBamg(gf);
- }
laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
}
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index 6809d494e3de90be5ed8226a4797840921e8f996..a4463d38075a1e03acdaa03ce7096ab40b2498db 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -17,6 +17,7 @@
#include <list>
#include <map>
#include "BackgroundMesh.h"
+#include "meshGFaceDelaunayInsertion.h"
#if defined(HAVE_BAMG)
@@ -71,8 +72,10 @@ void meshGFaceBamg(GFace *gf){
std::list<GEdge*> edges = gf->edges();
std::set<GEdge*> mySet;
std::list<GEdge*>::iterator it = edges.begin();
+ bool hasCompounds = false;
while(it != edges.end()){
if((*it)->getCompound()){
+ hasCompounds = true;
GEdge *gec = (GEdge*)(*it)->getCompound();
mySet.insert(gec);
}
@@ -90,6 +93,12 @@ void meshGFaceBamg(GFace *gf){
bcVertex.insert((*it)->lines[i]->getVertex(1));
}
}
+ //if has compound edges use a frontal mesher to remesh compound surface
+ if (hasCompounds){
+ printf("tris compound = %d \n", gf->triangles.size());
+ bowyerWatsonFrontal(gf);
+ printf("after tris compound = %d \n", gf->triangles.size());
+ }
//fill mesh data fo bamg (bamgVertices, bamgTriangles, bamgBoundary)
std::set<MVertex*> all;
@@ -230,23 +239,6 @@ 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/meshMetric.cpp b/Mesh/meshMetric.cpp
index 1a8f0825806bd47e947609665ee5ea5adc20aaf3..4a65f0d6f2bb0c06e6e6fb51159a712d5cf988f2 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -22,11 +22,32 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l
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){
_dim = gm->getDim();
- std::map<int, MVertex*> vertexMap;
std::map<MElement*, MElement*> newP;
std::map<MElement*, MElement*> newD;
@@ -34,7 +55,7 @@ meshMetric::meshMetric(GModel *gm, int technique, simpleFunction<double> *fct, s
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);
+ MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
}
}
@@ -43,7 +64,7 @@ meshMetric::meshMetric(GModel *gm, int technique, simpleFunction<double> *fct, s
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);
+ MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
}
}
@@ -58,7 +79,9 @@ meshMetric::meshMetric(GModel *gm, int technique, simpleFunction<double> *fct, s
_technique = (MetricComputationTechnique)technique;
computeMetric();
+ //printf("computing metric \n");
//printMetric("toto.pos");
+ //exit(1);
}
meshMetric::~meshMetric(){
@@ -85,7 +108,7 @@ void meshMetric::computeHessian( v2t_cont adj){
while (it != adj.end()) {
std::vector<MElement*> lt = it->second;
MVertex *ver = it->first;
- while (lt.size() < 10) increaseStencil(ver,adj,lt);
+ while (lt.size() < 12) increaseStencil(ver,adj,lt); //<7
// if ( ver->onWhat()->dim() < _dim ){
// while (lt.size() < 12){
// increaseStencil(ver,adj,lt);
@@ -133,7 +156,7 @@ void meshMetric::computeHessian( v2t_cont adj){
dgrads[ITER-1][ver] = SVector3(result(1),result(2),_dim==2? 0.0:result(3));
++it;
}
- if (_technique != meshMetric::HESSIAN)break;
+ if (_technique == meshMetric::LEVELSET) break;
}
}
@@ -142,7 +165,7 @@ void meshMetric::computeMetric(){
v2t_cont adj;
buildVertexToElement (_elements,adj);
- //printf("%d elements are considered in the meshMetric \n",_elements.size());
+ //printf("%d elements are considered in the meshMetric \n",(int)_elements.size());
computeValues(adj);
computeHessian(adj);
@@ -251,21 +274,6 @@ void meshMetric::computeMetric(){
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));
@@ -300,12 +308,11 @@ void meshMetric::computeMetric(){
//smoothMetric (sol);
//curvatureContributionToMetric();
-
putOnNewView();
}
-double meshMetric::operator() (double x, double y, double z, GEntity *ge){
+double meshMetric::operator() (double x, double y, double z, GEntity *ge) {
SPoint3 xyz(x,y,z), uvw;
double initialTol = MElement::getTolerance();
MElement::setTolerance(1.e-4);
@@ -324,7 +331,7 @@ double meshMetric::operator() (double x, double y, double z, GEntity *ge){
return 1.e22;
}
-void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge){
+void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge) {
metr = SMetric3(1.e-22);
SPoint3 xyz(x,y,z), uvw;
double initialTol = MElement::getTolerance();
@@ -360,30 +367,35 @@ void meshMetric::printMetric(const char* n) const{
FILE *f = fopen (n,"w");
fprintf(f,"View \"\"{\n");
- std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
+
//std::map<MVertex*,fullMatrix<double> >::const_iterator it = _hessian.begin();
+ std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
for (; it != _nodalMetrics.end(); ++it){
//for (; it != _hessian.end(); ++it){
- 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)
- );
+ MVertex *v = it->first;
+ // double lapl = getLaplacian(v);
+ // 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)
+ );
+
}
fprintf(f,"};\n");
fclose (f);
}
-double meshMetric::getLaplacian (MVertex *v){
-
- std::map<MVertex*, fullMatrix<double> >::iterator it = _hessian.find(v);
+double meshMetric::getLaplacian (MVertex *v) {
+ MVertex *vNew = _vertexMap[v->getNum()];
+ std::map<MVertex*, fullMatrix<double> >::const_iterator it = _hessian.find(vNew);
fullMatrix<double> h = it->second;
double laplace = h(0,0)+h(1,1)+h(2,2);
return laplace;
}
-/*void dgMeshMetric::curvatureContributionToMetric (){
+/*void meshMetric::curvatureContributionToMetric (){
std::map<MVertex*,SMetric3>::iterator it = _nodalMetrics.begin();
std::map<MVertex*,double>::iterator its = _nodalSizes.begin();
for (; it != _nodalMetrics.end(); ++it,++its){
@@ -405,7 +417,7 @@ double meshMetric::getLaplacian (MVertex *v){
}*/
/*
- void dgMeshMetric::smoothMetric (dgDofContainer *sol){
+ void meshMetric::smoothMetric (dgDofContainer *sol){
return;
dgGroupCollection *groups = sol->getGroups();
std::set<MEdge,Less_Edge> edges;
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 84b63e5d3c194f0f9280a4a529ea54ab69b1bb99..7c72d941b67e409874582ec74509b514a58097c3 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -22,6 +22,7 @@ class meshMetric: public Field {
std::vector<MElement*> _elements;
MElementOctree *_octree;
+ std::map<int, MVertex*> _vertexMap;
std::map<MVertex*,double> vals;
std::map<MVertex*,SVector3> grads,dgrads[3];
@@ -30,11 +31,13 @@ class meshMetric: public Field {
std::map<MVertex*,double> _nodalSizes, _detMetric;
std::map<MVertex*,fullMatrix<double> > _hessian;
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();
inline SMetric3 metricAtVertex (MVertex* v) {return _nodalMetrics[v];}
- void computeMetric();
+ void computeMetric() ;
void computeValues( v2t_cont adj);
void computeHessian( v2t_cont adj);
void setAsBackgroundMesh (GModel *gm);
@@ -49,7 +52,7 @@ class meshMetric: public Field {
{
return "Anisotropic size field based on hessian of a given function";
}
- virtual double operator() (double x, double y, double z, GEntity *ge=0);
+ 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;
diff --git a/contrib/bamg/bamg-gmsh.cpp b/contrib/bamg/bamg-gmsh.cpp
index 4ac8de89f9028cfe0f3d640f2f436266ab2a056c..999226ceed6bb579042f7fd45f0e3fba718f214e 100644
--- a/contrib/bamg/bamg-gmsh.cpp
+++ b/contrib/bamg/bamg-gmsh.cpp
@@ -569,7 +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);
+ //nTh->Write("bamg.mesh",bamg::Triangles::AutoMesh);
Mesh2 * g= bamg2msh(nTh,true);
delete nTh;
diff --git a/gmshpy/gmshMesh.i b/gmshpy/gmshMesh.i
index cba7a004aef6c80daebb0d25a25d42941f6bcc13..7e57df9349ef047230c21fc5912bb43380e48290 100644
--- a/gmshpy/gmshMesh.i
+++ b/gmshpy/gmshMesh.i
@@ -15,6 +15,11 @@
#include "meshMetric.h"
%}
+%include std_vector.i
+namespace std {
+ %template(DoubleVector) vector<double, std::allocator<double> >;
+}
+
%include "GmshConfig.h"
%include "Generator.h"
%include "highOrderTools.h"