diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 81e03d992f15fdd55a18a16174cec322f071ddc4..a4499b78692cbcaa95c1dc78ccfe872c3b162d1f 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -51,23 +51,25 @@ static void statistics_histogram_cb(Fl_Widget *w, void *data)
new PView("Disto", "# Elements", x, y);
}
else{
- std::vector<GEntity*> entities;
- GModel::current()->getEntities(entities);
+
+ std::vector<GEntity*> entities_;
+ GModel::current()->getEntities(entities_);
std::map<int, std::vector<double> > d;
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < 2) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- if(name == "Gamma3D")
- d[e->getNum()].push_back(e->gammaShapeMeasure());
- else if(name == "Eta3D")
- d[e->getNum()].push_back(e->etaShapeMeasure());
- else if(name == "Rho3D")
- d[e->getNum()].push_back(e->rhoShapeMeasure());
- else
- d[e->getNum()].push_back(e->distoShapeMeasure());
+ for(unsigned int i = 0; i < entities_.size(); i++){
+ if(entities_[i]->dim() < 2) continue;
+ for(unsigned int j = 0; j < entities_[i]->getNumMeshElements(); j++){
+ MElement *e = entities_[i]->getMeshElement(j);
+ if(name == "Gamma3D")
+ d[e->getNum()].push_back(e->gammaShapeMeasure());
+ else if(name == "Eta3D")
+ d[e->getNum()].push_back(e->etaShapeMeasure());
+ else if(name == "Rho3D")
+ d[e->getNum()].push_back(e->rhoShapeMeasure());
+ else
+ d[e->getNum()].push_back(e->distoShapeMeasure());
}
}
+
name.resize(name.size() - 2);
new PView(name, "ElementData", GModel::current(), d);
}
@@ -186,7 +188,7 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
void statisticsWindow::compute(bool elementQuality)
{
//emi hack - MINIMUM ANGLES
- /*double minAngle = 120.0;
+ double minAngle = 120.0;
double meanAngle = 0.0;
int count = 0;
std::vector<GEntity*> entities;
@@ -194,17 +196,97 @@ void statisticsWindow::compute(bool elementQuality)
std::map<int, std::vector<double> > d;
for(unsigned int i = 0; i < entities.size(); i++){
if(entities[i]->dim() < 2) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- double angle = e->angleShapeMeasure();
- minAngle = std::min(minAngle, angle);
- meanAngle += angle;
- count++;
- }
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ double angle = e->angleShapeMeasure();
+ minAngle = std::min(minAngle, angle);
+ meanAngle += angle;
+ count++;
+ }
}
meanAngle = meanAngle / count;
- printf("Angles = min=%g av=%g \n", minAngle, meanAngle);*/
- //hack emi
+ printf("Angles = min=%g av=%g \n", minAngle, meanAngle);
+ // Emi hack - MESH DEGREE VERTICES
+ // std::vector<GEntity*> entities;
+ std::set<MEdge, Less_Edge> edges;
+ GModel::current()->getEntities(entities);
+ std::map<MVertex*, int > vert2Deg;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ // printf("entity dim =%d tag=%d \n", entities[i]->dim() , entities[i]->tag());
+ if(entities[i]->dim() < 2 ) continue;
+ if(entities[i]->tag() < 100) continue;
+ // printf("continuing \n");
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ for(unsigned int k = 0; k < e->getNumEdges(); k++){
+ edges.insert(e->getEdge(k));
+ }
+ for(unsigned int k = 0; k < e->getNumVertices(); k++){
+ MVertex *v = e->getVertex(k);
+ if (v->onWhat()->dim() < 2) continue;
+ std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
+ if (it == vert2Deg.end()) {
+ vert2Deg.insert(std::make_pair(v,1));
+ }
+ else{
+ int nbE = it->second+1;
+ it->second = nbE;
+ }
+ }
+ }
+ }
+ int dMin = 10;
+ int dMax = 0;
+ int d4 = 0;
+ int nbElems = vert2Deg.size();
+ std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
+ for(; itmap !=vert2Deg.end(); itmap++){
+ MVertex *v = itmap->first;
+ int nbE = itmap->second;
+ dMin = std::min(nbE, dMin);
+ dMax = std::max(nbE, dMax);
+ if (nbE == 4) d4 += 1;
+ }
+ if (nbElems > 0)
+ printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
+
+ FieldManager *fields = GModel::current()->getFields();
+ Field *f = fields->get(fields->background_field);
+ int nbEdges = edges.size();
+ printf("nb edges =%d \n", nbEdges);
+ system("rm qualEdges.txt");
+ FILE *fp = fopen("qualEdges.txt", "w");
+ std::vector<int> qualE;
+ int nbS = 50;
+ qualE.resize(nbS);
+ if(fields->background_field > 0){
+ printf("found field \n");
+ std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ double sum = 0;
+ for (; it !=edges.end();++it){
+ MVertex *v0 = it->getVertex(0);
+ MVertex *v1 = it->getVertex(1);
+ double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
+ (v0->y()-v1->y())*(v0->y()-v1->y())+
+ (v0->z()-v1->z())*(v0->z()-v1->z()));
+ double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat());
+ double el = l/lf;
+ int index = (int) ceil(el*nbS*0.5);
+ qualE[index]+= 1;
+ double e = (l>lf) ? lf/l : l/lf;
+ sum += e - 1.0;
+ }
+ double tau = exp ((1./edges.size()) * sum);
+ printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
+
+ double ibegin = 2./(2*nbS);
+ double inext = 2./nbS;
+ for (int i= 0; i< qualE.size(); i++){
+ fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ }
+
+ }
+ fclose(fp);
//Emi hack - MESH DEGREE VERTICES
#if 0
std::vector<GEntity*> entities;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index e62438a7ccfa877f37cedbf40e147295e4af6894..280e5554d85b8c491e682567d812436704da6def 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -450,6 +450,7 @@ void GFaceCompound::fillNeumannBCS() const
bool GFaceCompound::trivial() const
{
+ return false;
if(_compound.size() == 1 &&
(*(_compound.begin()))->getNativeType() == GEntity::OpenCascadeModel &&
(*(_compound.begin()))->geomType() != GEntity::DiscreteSurface &&
@@ -1763,9 +1764,24 @@ void GFaceCompound::buildOct() const
_gfct[count].p2 = it1->second;
_gfct[count].p3 = it2->second;
if((*it)->geomType() != GEntity::DiscreteSurface){
- reparamMeshVertexOnFace(t->getVertex(0), *it, _gfct[count].gfp1);
- reparamMeshVertexOnFace(t->getVertex(1), *it, _gfct[count].gfp2);
- reparamMeshVertexOnFace(t->getVertex(2), *it, _gfct[count].gfp3);
+ // take care of the seam !!!!
+ if (t->getVertex(0)->onWhat()->dim() == 2){
+ reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(1),*it, _gfct[count].gfp1, _gfct[count].gfp2);
+ reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(2),*it, _gfct[count].gfp1, _gfct[count].gfp3);
+ }
+ else if (t->getVertex(1)->onWhat()->dim() == 2){
+ reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(0),*it, _gfct[count].gfp2, _gfct[count].gfp1);
+ reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(2),*it, _gfct[count].gfp2, _gfct[count].gfp3);
+ }
+ else if (t->getVertex(2)->onWhat()->dim() == 2){
+ reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(0),*it, _gfct[count].gfp3, _gfct[count].gfp1);
+ reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(1),*it, _gfct[count].gfp3, _gfct[count].gfp2);
+ }
+ else {
+ reparamMeshVertexOnFace(t->getVertex(0), *it, _gfct[count].gfp1);
+ reparamMeshVertexOnFace(t->getVertex(1), *it, _gfct[count].gfp2);
+ reparamMeshVertexOnFace(t->getVertex(2), *it, _gfct[count].gfp3);
+ }
}
_gfct[count].v1 = SPoint3(t->getVertex(0)->x(), t->getVertex(0)->y(),
t->getVertex(0)->z());
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index fd431b0e54f2e5a43d5f6412c29c8c703b2aecef..b59e0e4422456faea67090c587141b70c89721e0 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -324,10 +324,10 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
else if (p1.size() == 1 && p2.size() == 2){
double d1 =
(p1[0].x() - p2[0].x()) * (p1[0].x() - p2[0].x()) +
- (p1[0].x() - p2[0].y()) * (p1[0].y() - p2[0].y());
+ (p1[0].y() - p2[0].y()) * (p1[0].y() - p2[0].y());
double d2 =
(p1[0].x() - p2[1].x()) * (p1[0].x() - p2[1].x()) +
- (p1[0].x() - p2[1].y()) * (p1[0].y() - p2[1].y());
+ (p1[0].y() - p2[1].y()) * (p1[0].y() - p2[1].y());
param1 = p1[0];
param2 = d2 < d1 ? p2[1] : p2[0];
return true;
@@ -335,10 +335,10 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
else if (p2.size() == 1 && p1.size() == 2){
double d1 =
(p2[0].x() - p1[0].x()) * (p2[0].x() - p1[0].x()) +
- (p2[0].x() - p1[0].y()) * (p2[0].y() - p1[0].y());
+ (p2[0].y() - p1[0].y()) * (p2[0].y() - p1[0].y());
double d2 =
(p2[0].x() - p1[1].x()) * (p2[0].x() - p1[1].x()) +
- (p2[0].x() - p1[1].y()) * (p2[0].y() - p1[1].y());
+ (p2[0].y() - p1[1].y()) * (p2[0].y() - p1[1].y());
param1 = d2 < d1 ? p1[1] : p1[0];
param2 = p2[0];
return true;
@@ -346,6 +346,9 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
else if(p1.size() > 1 && p2.size() > 1){
param1 = p1[0];
param2 = p2[0];
+
+ printf("NO WAY : TWO VERTICES ON THE SEAM, CANNOT CHOOSE\n");
+
// shout, both vertices are on seams
return false;
}
@@ -357,6 +360,8 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
}
}
+
+
bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 ¶m,
bool onSurface)
{
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index a931a7a9d0f820057797e98103bbc4ebf69dfb47..3e6e52c1746cb365715fcbcf2fd76e8dd81c82c8 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -145,10 +145,10 @@ class MFaceVertex : public MVertex{
MVertexBoundaryLayerData* bl_data;
MFaceVertex(double x, double y, double z, GEntity *ge, double u, double v, int num = 0)
- : MVertex(x, y, z, ge, num), _u(u), _v(v)
+ : MVertex(x, y, z, ge, num), _u(u), _v(v),bl_data(0)
{
}
- virtual ~MFaceVertex(){}
+ virtual ~MFaceVertex(){if(bl_data) delete bl_data;}
virtual bool getParameter(int i, double &par) const { par = (i ? _v : _u); return true; }
virtual bool setParameter(int i, double par)
{
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index fdbdedea78f9fee43aa9c0ebf276a0348d03722e..ebf4db718a2022a9ec4b4df38ed4323756901263 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -53,7 +53,7 @@ void OCCFace::setup()
}
else{
l_wire.push_back(e);
- Msg::Debug("Edge %d ori %d", e->tag(), edge.Orientation());
+ Msg::Debug("Edge %d (%d --> %d) ori %d", e->tag(), e->getBeginVertex()->tag(), e->getEndVertex()->tag(), edge.Orientation());
e->addFace(this);
if(!e->is3D()){
OCCEdge *occe = (OCCEdge*)e;
@@ -95,17 +95,12 @@ void OCCFace::setup()
umax += fabs(du) / 100.0;
vmax += fabs(dv) / 100.0;
occface = BRep_Tool::Surface(s);
- // specific parametrization for a sphere.
- _isSphere = isSphere(_radius,_center);
- if (_isSphere){
- for (std::list<GEdge*>::iterator it = l_edges.begin();
- it != l_edges.end(); ++it){
- GEdge *ge = *it;
- if (ge->isSeam(this)){
-
- }
- }
- }
+
+ // std::list<GEdge*>::const_iterator it = l_edges.begin();
+ // for (; it != l_edges.end(); ++it){
+ // printf("edge %d : %d %d iseam %d \n",(*it)->tag(),(*it)->getBeginVertex()->tag(),(*it)->getEndVertex()->tag(),(*it)->isSeam(this));
+ // }
+
}
Range<double> OCCFace::parBounds(int i) const
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 20e6594644ce449cdf369c6304b0b9d4e7c19445..7436bc7a510317408082580da8388b3a4d9d37fa 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -33,6 +33,22 @@ SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
return iv;
}
+// preserve orientation of m1 !!!
+SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2)
+{
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ m1.eig(V,S,true);
+ SVector3 v0(V(0,0),V(1,0),V(2,0));
+ SVector3 v1(V(0,1),V(1,1),V(2,1));
+ SVector3 v2(V(0,2),V(1,2),V(2,2));
+ double l0 = std::max(dot(v0,m1,v0),dot(v0,m2,v0));
+ double l1 = std::max(dot(v1,m1,v1),dot(v1,m2,v1));
+ double l2 = std::max(dot(v2,m1,v2),dot(v2,m2,v2));
+ SMetric3 iv(l0,l1,l2,v0,v1,v2);
+ return iv;
+}
+
// (1-t) * m1 + t * m2
SMetric3 interpolation (const SMetric3 &m1,
const SMetric3 &m2,
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index 3562f71e904800419d5a308000f2257e994d4444..a039a045fb4884f684c3a77e049f99c4df92a567 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -163,6 +163,8 @@ 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() ) ;
}
+SMetric3 intersection_conserveM1 (const SMetric3 &m1,
+ const SMetric3 &m2);
// compute the largest inscribed ellipsoid...
SMetric3 intersection (const SMetric3 &m1,
const SMetric3 &m2);
diff --git a/Graphics/drawGeom.cpp b/Graphics/drawGeom.cpp
index b6269dad7d2c401bec76344d6e0317d81726370b..f4634cffe6f078e559012be70ce151274e16fb43 100644
--- a/Graphics/drawGeom.cpp
+++ b/Graphics/drawGeom.cpp
@@ -421,7 +421,7 @@ class drawGRegion {
void operator () (GRegion *r)
{
if(!r->getVisibility()) return;
- if(r->geomType() == GEntity::DiscreteVolume) return;
+ // if(r->geomType() == GEntity::DiscreteVolume) return;
bool select = (_ctx->render_mode == drawContext::GMSH_SELECT &&
r->model() == GModel::current());
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 1f4dadee434042d7f657efec96d95192aedbd95c..4c3ceb2e0a5948503c0c8831f1d17fb6f322c0c5 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -457,7 +457,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
_sizes[itv2->first] = MAX_LC;
}
}
- // ensure that other criteria are fullfilled
+ // ensure that other criteria are fullfilled
updateSizes(_gf);
// compute optimal mesh orientations
@@ -717,6 +717,29 @@ void backgroundMesh::updateSizes(GFace *_gf)
itv->second = std::max(itv->second, CTX::instance()->mesh.lcMin);
itv->second = std::min(itv->second, CTX::instance()->mesh.lcMax);
}
+ // do not allow large variations in the size field
+ // INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 43, 1143–1165 (1998)
+ // MESH GRADATION CONTROL, BOROUCHAKI, HECHT, FREY
+ std::set<MEdge,Less_Edge> edges;
+ for (int i=0;i < _triangles.size();i++){
+ for (int j=0;j< _triangles[i]->getNumEdges();j++){
+ edges.insert(_triangles[i]->getEdge(j));
+ }
+ }
+ const double _beta = 1.3;
+ for (int i=0;i<0;i++){
+ std::set<MEdge,Less_Edge>::iterator it = edges.begin();
+ for ( ; it != edges.end(); ++it){
+ MVertex *v0 = it->getVertex(0);
+ MVertex *v1 = it->getVertex(1);
+ MVertex *V0 = _2Dto3D[v0];
+ MVertex *V1 = _2Dto3D[v1];
+ std::map<MVertex*,double>::iterator s0 = _sizes.find(V0);
+ std::map<MVertex*,double>::iterator s1 = _sizes.find(V1);
+ if (s0->second < s1->second)s1->second = std::min(s1->second,_beta*s0->second);
+ else s0->second = std::min(s0->second,_beta*s1->second);
+ }
+ }
}
double backgroundMesh::operator() (double u, double v, double w) const
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 329b6105876adb66a5129ddb032aee71dd79d002..1e6e62b777fe5bf6947b018b0189abbc02a1f037 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -12,7 +12,7 @@ set(SRC
meshGFaceTransfinite.cpp meshGFaceExtruded.cpp
meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp
meshGFaceLloyd.cpp meshGFaceOptimize.cpp
- meshGFaceQuadrilateralize.cpp
+ meshGFaceQuadrilateralize.cpp meshGFaceBoundaryLayers.cpp
meshGRegion.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 7f2a5304b064706002a782e4244bb1115b4eefad..88205789574eb00c72415ebf418c0ba0012c29c7 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1422,7 +1422,7 @@ class AttractorAnisoCurveField : public Field {
it != edges_id.end(); ++it) {
Curve *c = FindCurve(*it);
if(c) {
- for(int i = 0; i < n_nodes_by_edge; i++) {
+ for(int i = 1; i < n_nodes_by_edge-1; i++) {
double u = (double)i / (n_nodes_by_edge - 1);
Vertex V = InterpolateCurve(c, u, 0);
zeronodes[k][0] = V.Pos.X;
@@ -1437,7 +1437,7 @@ class AttractorAnisoCurveField : public Field {
else {
GEdge *e = GModel::current()->getEdgeByTag(*it);
if(e) {
- for(int i = 0; i < n_nodes_by_edge; i++) {
+ for(int i = 1; i < n_nodes_by_edge-1; i++) {
double u = (double)i / (n_nodes_by_edge - 1);
Range<double> b = e->parBounds(0);
double t = b.low() + u * (b.high() - b.low());
@@ -1493,6 +1493,16 @@ class AttractorField : public Field
Field *_xField, *_yField, *_zField;
int n_nodes_by_edge;
public:
+ AttractorField(int dim, int tag, int nbe) : kdtree(0), zeronodes(0), n_nodes_by_edge(nbe)
+ {
+ index = new ANNidx[1];
+ dist = new ANNdist[1];
+ if (dim == 0) nodes_id.push_back(tag);
+ else if (dim == 1) edges_id.push_back(tag);
+ else if (dim == 2) faces_id.push_back(tag);
+ _xFieldId = _yFieldId = _zFieldId = -1;
+ update_needed = true;
+ }
AttractorField() : kdtree(0), zeronodes(0)
{
index = new ANNidx[1];
@@ -1555,7 +1565,7 @@ class AttractorField : public Field
annDeallocPts(zeronodes);
delete kdtree;
}
- int totpoints = nodes_id.size() + n_nodes_by_edge * edges_id.size() +
+ int totpoints = nodes_id.size() + (n_nodes_by_edge-2) * edges_id.size() +
n_nodes_by_edge * n_nodes_by_edge * faces_id.size();
if(totpoints){
zeronodes = annAllocPts(totpoints, 3);
@@ -1567,22 +1577,22 @@ class AttractorField : public Field
Vertex *v = FindPoint(*it);
if(v) {
getCoord(v->Pos.X, v->Pos.Y, v->Pos.Z, zeronodes[k][0], zeronodes[k][1], zeronodes[k][2]);
- k++;
+ _infos[k++] = AttractorInfo(*it,0,0,0);
}
else {
GVertex *gv = GModel::current()->getVertexByTag(*it);
if(gv) {
getCoord(gv->x(), gv->y(), gv->z(), zeronodes[k][0], zeronodes[k][1], zeronodes[k][2], gv);
- k++;
+ _infos[k++] = AttractorInfo(*it,0,0,0);
}
}
}
for(std::list<int>::iterator it = edges_id.begin();
it != edges_id.end(); ++it) {
Curve *c = FindCurve(*it);
- GEdge *e = GModel::current()->getEdgeByTag(*it);
+ GEdge *e = GModel::current()->getEdgeByTag(*it);
if(c && !e) {
- for(int i = 0; i < n_nodes_by_edge; i++) {
+ for(int i = 1; i < n_nodes_by_edge -1 ; i++) {
double u = (double)i / (n_nodes_by_edge - 1);
Vertex V = InterpolateCurve(c, u, 0);
getCoord(V.Pos.X, V.Pos.Y, V.Pos.Z, zeronodes[k][0], zeronodes[k][1], zeronodes[k][2]);
@@ -1591,7 +1601,7 @@ class AttractorField : public Field
}
else {
if(e) {
- for(int i = 0; i < n_nodes_by_edge; i++) {
+ for(int i = 1; i < n_nodes_by_edge - 1; i++) {
double u = (double)i / (n_nodes_by_edge - 1);
Range<double> b = e->parBounds(0);
double t = b.low() + u * (b.high() - b.low());
@@ -1648,203 +1658,235 @@ class AttractorField : public Field
}
};
-class BoundaryLayerField : public Field {
- int iField;
- double hwall_n,hwall_t,ratio,hfar;
-public:
- virtual bool isotropic () const {return false;}
- virtual const char *getName()
- {
- return "BoundaryLayer";
- }
- virtual std::string getDescription()
- {
- return "hwall * ratio^(dist/hwall)";
- }
- BoundaryLayerField()
- {
- iField = 0;
- hwall_n = .1;
- hwall_t = .5;
- hfar = 1;
- ratio = 1.1;
- options["IField"] = new FieldOptionInt
- (iField, "Index of the field that contains the distance function");
- options["hwall_n"] = new FieldOptionDouble
- (hwall_n, "Mesh Size Normal to the The Wall");
- options["hwall_t"] = new FieldOptionDouble
- (hwall_t, "Mesh Size Tangent to the Wall");
- options["ratio"] = new FieldOptionDouble
- (ratio, "Size Ratio Between Two Successive Layers");
- options["hfar"] = new FieldOptionDouble
- (hfar, "Element size far from the wall");
- }
- virtual double operator() (double x, double y, double z, GEntity *ge=0)
- {
- Field *field = GModel::current()->getFields()->get(iField);
- if(!field || iField == id) {
- return MAX_LC;
+const char *BoundaryLayerField ::getName()
+{
+ return "BoundaryLayer";
+}
+std::string BoundaryLayerField :: getDescription()
+{
+ return "hwall * ratio^(dist/hwall)";
+}
+BoundaryLayerField :: BoundaryLayerField()
+{
+ hwall_n = .1;
+ hwall_t = .5;
+ hfar = 1;
+ ratio = 1.1;
+ thickness = 1.e-2;
+ options["NodesList"] = new FieldOptionList
+ (nodes_id, "Indices of nodes in the geometric model", &update_needed);
+ options["EdgesList"] = new FieldOptionList
+ (edges_id, "Indices of curves in the geometric model for which a boundary layer is needed", &update_needed);
+ // options["IField"] = new FieldOptionInt
+ // (iField, "Index of the field that contains the distance function");
+ options["hwall_n"] = new FieldOptionDouble
+ (hwall_n, "Mesh Size Normal to the The Wall");
+ options["hwall_t"] = new FieldOptionDouble
+ (hwall_t, "Mesh Size Tangent to the Wall");
+ options["ratio"] = new FieldOptionDouble
+ (ratio, "Size Ratio Between Two Successive Layers");
+ options["hfar"] = new FieldOptionDouble
+ (hfar, "Element size far from the wall");
+ options["thickness"] = new FieldOptionDouble
+ (thickness, "Maximal thickness of the boundary layer");
+}
+
+double BoundaryLayerField :: operator() (double x, double y, double z, GEntity *ge)
+{
+ double dist = 1.e22;
+ AttractorField *cc;
+ for (std::list<AttractorField*>::iterator it = _att_fields.begin();
+ it != _att_fields.end(); ++it){
+ double cdist = (*(*it)) (x, y, z);
+ if (cdist < dist){
+ cc = *it;
+ dist = cdist;
}
- const double dist = (*field) (x, y, z);
- const double lc = dist*(ratio-1) + hwall_n;
- // double lc = hwall * pow (ratio, dist / hwall);
- return std::min (hfar,lc);
}
- virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
- {
- Field *field = GModel::current()->getFields()->get(iField);
- if(!field || iField == id) {
- metr(0,0) = 1/(MAX_LC*MAX_LC);
- metr(1,1) = 1/(MAX_LC*MAX_LC);
- metr(2,2) = 1/(MAX_LC*MAX_LC);
- metr(0,1) = metr(0,2) = metr(1,2) = 0;
- return;
- }
- const double dist = (*field) (x, y, z);
-
- // dist = hwall -> lc = hwall * ratio
- // dist = hwall (1+ratio) -> lc = hwall ratio ^ 2
- // dist = hwall (1+ratio+ratio^2) -> lc = hwall ratio ^ 3
- // dist = hwall (1+ratio+ratio^2+...+ratio^{m-1}) = (ratio^{m} - 1)/(ratio-1) -> lc = hwall ratio ^ m
- // -> find m
- // dist/hwall = (ratio^{m} - 1)/(ratio-1)
- // (dist/hwall)*(ratio-1) + 1 = ratio^{m}
- // lc = dist*(ratio-1) + hwall
- const double ll1 = dist*(ratio-1) + hwall_n;
- double lc_n = std::min(ll1,hfar);
- const double ll2 = dist*(ratio-1) + hwall_t;
- double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
-
- SVector3 t1,t2,t3;
- double L1,L2,L3;
-
- AttractorField *cc = dynamic_cast<AttractorField*> (field);
- if (cc){
- std::pair<AttractorInfo,SPoint3> pp = cc->getAttractorInfo();
- if (pp.first.dim ==1){
- GEdge *e = GModel::current()->getEdgeByTag(pp.first.ent);
-
- // the tangent size at this point is the size of the
- // 1D mesh at this point !
- // hack : use curvature
- if(CTX::instance()->mesh.lcFromCurvature){
- double Crv = e->curvature(pp.first.u);
- double lc = Crv > 0 ? 2 * M_PI / Crv / CTX::instance()->mesh.minCircPoints : 1.e-22;
- const double ll2b = dist*(ratio-1) + lc;
- double lc_tb = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2b,hfar));
- lc_t = std::min(lc_t,lc_tb);
- }
+ current_distance = dist;
+ // const double dist = (*field) (x, y, z);
+ // current_distance = dist;
+ const double lc = dist*(ratio-1) + hwall_t;
+ // double lc = hwall * pow (ratio, dist / hwall);
+ // printf("coucou\n");
+ return std::min (hfar,lc);
+}
+void BoundaryLayerField :: operator() (AttractorField *cc, double dist, double x, double y, double z, SMetric3 &metr, GEntity *ge)
+{
+ // dist = hwall -> lc = hwall * ratio
+ // dist = hwall (1+ratio) -> lc = hwall ratio ^ 2
+ // dist = hwall (1+ratio+ratio^2) -> lc = hwall ratio ^ 3
+ // dist = hwall (1+ratio+ratio^2+...+ratio^{m-1}) = (ratio^{m} - 1)/(ratio-1) -> lc = hwall ratio ^ m
+ // -> find m
+ // dist/hwall = (ratio^{m} - 1)/(ratio-1)
+ // (dist/hwall)*(ratio-1) + 1 = ratio^{m}
+ // lc = dist*(ratio-1) + hwall
- if (dist < hwall_t){
- L1 = lc_t;
- L2 = lc_n;
- L3 = lc_n;
- t1 = e->firstDer(pp.first.u);
- t1.normalize();
- if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
- t2 = SVector3(1,0,0);
- else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
- t2 = SVector3(0,1,0);
- else
- t2 = SVector3(0,0,1);
- t3 = crossprod(t1,t2);
- t3.normalize();
- t2 = crossprod(t3,t1);
- // printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
- }
- else {
- L1 = lc_t;
- L2 = lc_n;
- L3 = lc_t;
- GPoint p = e->point(pp.first.u);
- t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
- if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
- t1 = SVector3(1,0,0);
- else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
- t1 = SVector3(0,1,0);
- else
- t1 = SVector3(0,0,1);
- t2.normalize();
- t3 = crossprod(t1,t2);
- t3.normalize();
- t1 = crossprod(t3,t2);
- }
- }
- else {
- GFace *gf = GModel::current()->getFaceByTag(pp.first.ent);
-
- if (dist < hwall_t){
- L1 = lc_n;
- L2 = lc_t;
- L3 = lc_t;
- t1 = gf->normal(SPoint2(pp.first.u,pp.first.v));
- t1.normalize();
- if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
- t2 = SVector3(1,0,0);
- else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
- t2 = SVector3(0,1,0);
- else
- t2 = SVector3(0,0,1);
- t3 = crossprod(t1,t2);
- t3.normalize();
- t2 = crossprod(t3,t1);
- // printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
- }
- else {
- L1 = lc_t;
- L2 = lc_n;
- L3 = lc_t;
- GPoint p = gf->point(SPoint2(pp.first.u,pp.first.v));
- t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
- if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
- t1 = SVector3(1,0,0);
- else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
- t1 = SVector3(0,1,0);
- else
- t1 = SVector3(0,0,1);
- t2.normalize();
- t3 = crossprod(t1,t2);
- t3.normalize();
- t1 = crossprod(t3,t2);
- }
+ const double ll1 = dist*(ratio-1) + hwall_n;
+ double lc_n = std::min(ll1,hfar);
+ const double ll2 = dist*(ratio-1) + hwall_t;
+ double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
+
+ SVector3 t1,t2,t3;
+ double L1,L2,L3;
+
+ std::pair<AttractorInfo,SPoint3> pp = cc->getAttractorInfo();
+ if (pp.first.dim ==0){
+ L1 = lc_n;
+ L2 = lc_n;
+ L3 = lc_n;
+ GVertex *v = GModel::current()->getVertexByTag(pp.first.ent);
+ t1 = SVector3(v->x() -x,v->y() -y,v->z() -z);
+ t1.normalize();
+ // printf("%p %g %g %g\n",v,t1.x(),t1.y(),t1.z());
+ if (t1.norm() == 0)t1 = SVector3(1,0,0);
+ if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
+ t2 = SVector3(1,0,0);
+ else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
+ t2 = SVector3(0,1,0);
+ else
+ t2 = SVector3(0,0,1);
+ t3 = crossprod(t1,t2);
+ t3.normalize();
+ t2 = crossprod(t3,t1);
+ // printf("t1 %p %g %g %g\n",v,t1.x(),t1.y(),t1.z());
+ // printf("t2 %p %g %g %g\n",v,t2.x(),t2.y(),t2.z());
+ // printf("t3 %p %g %g %g\n",v,t3.x(),t3.y(),t3.z());
+ }
+ else if (pp.first.dim ==1){
+ GEdge *e = GModel::current()->getEdgeByTag(pp.first.ent);
+
+ // the tangent size at this point is the size of the
+ // 1D mesh at this point !
+ // hack : use curvature
+ /*
+ if(CTX::instance()->mesh.lcFromCurvature){
+ double Crv = e->curvature(pp.first.u);
+ double lc = Crv > 0 ? 2 * M_PI / Crv / CTX::instance()->mesh.minCircPoints : 1.e-22;
+ const double ll2b = dist*(ratio-1) + lc;
+ double lc_tb = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2b,hfar));
+ lc_t = std::min(lc_t,lc_tb);
}
+ */
+
+ if (dist < hwall_t){
+ L1 = lc_t;
+ L2 = lc_n;
+ L3 = lc_n;
+ t1 = e->firstDer(pp.first.u);
+ t1.normalize();
+ if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
+ t2 = SVector3(1,0,0);
+ else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
+ t2 = SVector3(0,1,0);
+ else
+ t2 = SVector3(0,0,1);
+ t3 = crossprod(t1,t2);
+ t3.normalize();
+ t2 = crossprod(t3,t1);
+ // printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
}
- else{
- double delta = std::min(CTX::instance()->lc / 1e2, dist);
- double gx =
- ((*field) (x + delta / 2, y, z) -
- (*field) (x - delta / 2, y, z)) / delta;
- double gy =
- ((*field) (x, y + delta / 2, z) -
- (*field) (x, y - delta / 2, z)) / delta;
- double gz =
- ((*field) (x, y, z + delta / 2) -
- (*field) (x, y, z - delta / 2)) / delta;
- SVector3 g = SVector3 (gx,gy,gz);
- g.normalize();
-
- if (fabs(g.x()) < fabs(g.y()) && fabs(g.x()) < fabs(g.z()))
- t1 = SVector3(1,0,0);
- else if (fabs(g.y()) < fabs(g.x()) && fabs(g.y()) < fabs(g.z()))
+ else {
+ L1 = lc_t;
+ L2 = lc_n;
+ L3 = lc_t;
+ GPoint p = e->point(pp.first.u);
+ t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
+ if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
+ t1 = SVector3(1,0,0);
+ else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
t1 = SVector3(0,1,0);
else
t1 = SVector3(0,0,1);
-
- t2 = crossprod(g,t1);
t2.normalize();
- t1 = crossprod(t2,g);
- t3 = g;
+ t3 = crossprod(t1,t2);
+ t3.normalize();
+ t1 = crossprod(t3,t2);
+ }
+ }
+ else {
+ GFace *gf = GModel::current()->getFaceByTag(pp.first.ent);
+
+ if (dist < hwall_t){
L1 = lc_n;
L2 = lc_t;
L3 = lc_t;
+ t1 = gf->normal(SPoint2(pp.first.u,pp.first.v));
+ t1.normalize();
+ if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
+ t2 = SVector3(1,0,0);
+ else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
+ t2 = SVector3(0,1,0);
+ else
+ t2 = SVector3(0,0,1);
+ t3 = crossprod(t1,t2);
+ t3.normalize();
+ t2 = crossprod(t3,t1);
+ // printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
+ }
+ else {
+ L1 = lc_t;
+ L2 = lc_n;
+ L3 = lc_t;
+ GPoint p = gf->point(SPoint2(pp.first.u,pp.first.v));
+ t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
+ if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
+ t1 = SVector3(1,0,0);
+ else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
+ t1 = SVector3(0,1,0);
+ else
+ t1 = SVector3(0,0,1);
+ t2.normalize();
+ t3 = crossprod(t1,t2);
+ t3.normalize();
+ t1 = crossprod(t3,t2);
+ }
+ }
+ metr = SMetric3(1./(L1*L1),
+ 1./(L2*L2),
+ 1./(L3*L3),
+ t1,t2,t3);
+}
+
+void BoundaryLayerField :: operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge)
+{
+
+ if (update_needed){
+ for(std::list<int>::iterator it = nodes_id.begin();
+ it != nodes_id.end(); ++it) {
+ _att_fields.push_back(new AttractorField(0,*it,100000));
+ }
+ for(std::list<int>::iterator it = edges_id.begin();
+ it != edges_id.end(); ++it) {
+ _att_fields.push_back(new AttractorField(1,*it,300000));
}
- metr = SMetric3(1./(L1*L1),
- 1./(L2*L2),
- 1./(L3*L3),
- t1,t2,t3);
+ update_needed = false;
}
-};
+
+ current_distance = 1.e22;
+ current_closest = 0;
+ SMetric3 v (1./MAX_LC);
+ std::vector<SMetric3> hop;
+ for (std::list<AttractorField*>::iterator it = _att_fields.begin();
+ it != _att_fields.end(); ++it){
+ double cdist = (*(*it)) (x, y, z);
+ SPoint3 CLOSEST= (*it)->getAttractorInfo().second;
+
+ SMetric3 localMetric;
+ (*this)(*it, cdist,x, y, z, localMetric, ge);
+ hop.push_back(localMetric);
+ //v = intersection(v,localMetric);
+ if (cdist < current_distance){
+ current_distance = cdist;
+ current_closest = *it;
+ v = localMetric;
+ _closest_point = CLOSEST;
+ }
+ }
+ // for (int i=0;i<hop.size();i++)v = intersection_conserveM1(v,hop[i]);
+ metr = v;
+}
#endif
template<class F> class FieldFactoryT : public FieldFactory {
diff --git a/Mesh/Field.h b/Mesh/Field.h
index fdb070d201ff3e5122d30700810b72f1de6ca91b..dba2d5a1505001a028a0e1857bfa2d08bb3752d0 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -101,4 +101,29 @@ class FieldManager : public std::map<int, Field*> {
void setBackgroundMesh(int iView);
};
+// Boundary Layer Field
+// specific field that allow to build boundary layers
+// is used both for anisotropic meshing and BL extrusion
+
+#if defined(HAVE_ANN)
+class AttractorField;// : public Field;
+
+class BoundaryLayerField : public Field {
+ std::list<AttractorField *> _att_fields;
+ std::list<int> nodes_id, edges_id, faces_id;
+ void operator() (AttractorField *cc, double dist, double x, double y, double z, SMetric3 &metr, GEntity *ge);
+public:
+ double hwall_n,hwall_t,ratio,hfar,thickness;
+ double current_distance;
+ SPoint3 _closest_point;
+ AttractorField *current_closest;
+ virtual bool isotropic () const {return false;}
+ virtual const char *getName();
+ virtual std::string getDescription();
+ BoundaryLayerField();
+ 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);
+};
+#endif
+
#endif
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index a89fe66a2b5430c940c45111fec087f2dad6e1ee..a2abd8b623b3955990d99e6f3930e57184a0c897 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -474,11 +474,14 @@ static void Mesh2D(GModel *m)
// lloyd optimization
if (CTX::instance()->mesh.optimizeLloyd){
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
- if((*it)->geomType()==GEntity::CompoundSurface){
+ if((*it)->geomType()==GEntity::CompoundSurface || (*it)->geomType()==GEntity::Plane){
smoothing smm(CTX::instance()->mesh.optimizeLloyd,6);
+ m->writeMSH("beforeLLoyd.msh");
smm.optimize_face(*it);
- int rec = (CTX::instance()->mesh.recombineAll || (*it)->meshAttributes.recombine);
- if(rec) recombineIntoQuads(*it);
+ int rec = 1;//(CTX::instance()->mesh.recombineAll || (*it)->meshAttributes.recombine);
+ m->writeMSH("afterLLoyd.msh");
+ if(rec){printf("recombine\n"); recombineIntoQuads(*it);}
+ m->writeMSH("afterRecombine.msh");
}
}
/*
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index aaa02290e3f338529e610513fa639ca980f91a36..4e854c560f533795c4542e80f4cd654dc507c292 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1236,9 +1236,9 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// printJacobians(m, "smoothness.pos");
// m->writeMSH("SMOOTHED.msh");
-
+ // FIXME !!
if (!linear){
- hot.ensureMinimumDistorsion(0.1);
+ //hot.ensureMinimumDistorsion(0.1);
checkHighOrderTriangles("Final mesh", m, bad, worst);
checkHighOrderTetrahedron("Final mesh", m, bad, worst);
}
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index adb99fec89ac752a3313cb4b183b024bba1788b1..31797c86685930945a37cc0b4785855e7582d653 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -309,7 +309,8 @@ void meshGEdge::operator() (GEdge *ge)
N = ge->meshAttributes.nbPointsTransfinite;
}
else{
- if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG)
+ if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG /*|| 1 || // FIXME
+ CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD */)
a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
CTX::instance()->mesh.lcIntegrationPrecision);
else
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 393818c585818546ea0863b5f34416c1dc415f6d..74774d09576fe7adb2cf5a8734fca1535121cb4f 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -40,6 +40,7 @@
#include "Context.h"
#include "multiscalePartition.h"
#include "meshGFaceLloyd.h"
+#include "meshGFaceBoundaryLayers.h"
static void copyMesh(GFace *source, GFace *target)
{
@@ -390,19 +391,256 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
return true;
}
+void BDS2GMSH ( BDS_Mesh *m, GFace *gf, std::map<BDS_Point*, MVertex*> &recoverMap){
+ {
+ std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
+ while (itp != m->points.end()){
+ BDS_Point *p = *itp;
+ if(recoverMap.find(p) == recoverMap.end()){
+ MVertex *v = new MFaceVertex
+ (p->X, p->Y, p->Z, gf, m->scalingU * p->u, m->scalingV * p->v);
+ recoverMap[p] = v;
+ gf->mesh_vertices.push_back(v);
+ }
+ ++itp;
+ }
+ }
+ {
+ std::list<BDS_Face*>::iterator itt = m->triangles.begin();
+ while (itt != m->triangles.end()){
+ BDS_Face *t = *itt;
+ if(!t->deleted){
+ BDS_Point *n[4];
+ t->getNodes(n);
+ MVertex *v1 = recoverMap[n[0]];
+ MVertex *v2 = recoverMap[n[1]];
+ MVertex *v3 = recoverMap[n[2]];
+ if(!n[3]){
+ // when a singular point is present, degenerated triangles
+ // may be created, for example on a sphere that contains one
+ // pole
+ if(v1 != v2 && v1 != v3 && v2 != v3)
+ gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ }
+ else{
+ MVertex *v4 = recoverMap[n[3]];
+ gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
+ }
+ }
+ ++itt;
+ }
+ }
+}
+
+bool meshGenerator(GFace *gf, int RECUR_ITER,
+ bool repairSelfIntersecting1dMesh,
+ bool onlyInitialMesh,
+ bool debug = true,
+ std::list<GEdge*> *replacement_edges = 0);
+
+static void addOrRemove ( MVertex *v1, MVertex *v2, std::set<MEdge,Less_Edge> & bedges){
+ MEdge e(v1,v2);
+ std::set<MEdge,Less_Edge>::iterator it = bedges.find(e);
+ if (it == bedges.end())bedges.insert(e);
+ else bedges.erase(it);
+}
+
+void modifyInitialMeshForTakingIntoAccountBoundaryLayers (GFace *gf){
+
+ BoundaryLayerColumns *_columns = buidAdditionalPoints2D (gf, M_PI/6.);
+
+ if (!_columns)return;
+
+ std::set<MEdge,Less_Edge> bedges;
+
+ std::vector<MQuadrangle*> blQuads;
+ std::vector<MTriangle*> blTris;
+ std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*>::iterator ite = edges.begin();
+ FILE *ff2 = fopen ("tato.pos","w");
+ fprintf(ff2,"View \" \"{\n");
+ std::set<MVertex*> verts;
+ while(ite != edges.end()){
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ MEdge dv(v1,v2);
+ addOrRemove(v1,v2,bedges);
+ int nbCol1 = _columns->getNbColumns(v1);
+ int nbCol2 = _columns->getNbColumns(v2);
+ if (nbCol1 > 0 && nbCol2 > 0){
+ const BoundaryLayerData & c1 = _columns->getColumn(v1,MEdge(v1,v2));
+ const BoundaryLayerData & c2 = _columns->getColumn(v2,MEdge(v1,v2));
+ int N = std::min(c1._column.size(),c2._column.size());
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v21,*v22;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ if (l == 0){
+ v11 = v1;
+ v12 = v2;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ }
+ MEdge dv2 (v21,v22);
+
+ //avoid convergent errors
+ if (dv2.length() < 0.5 * dv.length())break;
+ // printf("quadrangle generated\n");
+ blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
+ //blTris.push_back(new MTriangle(v11,v12,v22));
+ // blTris.push_back(new MTriangle(v11,v22,v21));
+ addOrRemove(v11,v12,bedges);
+ addOrRemove(v12,v22,bedges);
+ addOrRemove(v22,v21,bedges);
+ addOrRemove(v21,v11,bedges);
+ if(v11->onWhat() == gf)verts.insert(v11);
+ if(v21->onWhat() == gf)verts.insert(v21);
+ if(v12->onWhat() == gf)verts.insert(v12);
+ if(v22->onWhat() == gf)verts.insert(v22);
+ fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
+ v11->x(),v11->y(),v11->z(),
+ v12->x(),v12->y(),v12->z(),
+ v22->x(),v22->y(),v22->z(),
+ v21->x(),v21->y(),v21->z());
+ }
+ int M = std::max(c1._column.size(),c2._column.size());
+
+ /*
+ if (M>N) M = N+1;
+ // close with triangles
+ for (int l=N-1;l < M-1 ;++l){
+ MVertex *v11,*v12,*v21,*v22;
+ v11 = c1._column[l>=c1._column.size() ? c1._column.size() -1 : l];
+ v12 = c2._column[l>=c2._column.size() ? c2._column.size() -1 : l];
+ v21 = c1._column[(l+1)>=c1._column.size() ? c1._column.size() -1 : l+1];
+ v22 = c2._column[(l+1)>=c2._column.size() ? c2._column.size() -1 : l+1];
+
+ MTriangle *nt = (v11 == v21) ? new MTriangle(v11,v12,v22) : new MTriangle(v11,v12,v21) ;
+ blTris.push_back(nt);
+ v11 = nt->getVertex(0);
+ v12 = nt->getVertex(1);
+ v21 = nt->getVertex(2);
+
+ addOrRemove(v11,v12,bedges);
+ addOrRemove(v12,v21,bedges);
+ addOrRemove(v21,v11,bedges);
+ if(v11->onWhat() == gf)verts.insert(v11);
+ if(v21->onWhat() == gf)verts.insert(v21);
+ if(v12->onWhat() == gf)verts.insert(v12);
+ fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n",
+ v11->x(),v11->y(),v11->z(),
+ v12->x(),v12->y(),v12->z(),
+ v21->x(),v21->y(),v21->z());
+ }
+ */
+ }
+ }
+ ++ite;
+ }
+
+ if (1){
+ for (BoundaryLayerColumns::iterf itf = _columns->beginf();
+ itf != _columns->endf() ; ++itf){
+ MVertex *v = itf->first;
+ int nbCol = _columns->getNbColumns(v);
+
+ for (int i=0;i<nbCol-1;i++){
+ // printf("permut %d %d\n",permut[i],permut[i+1]);
+ const BoundaryLayerData & c1 = _columns->getColumn(v,i);
+ const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
+ int N = std::min(c1._column.size(),c2._column.size());
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v21,*v22;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ if (l == 0){
+ v11 = v;
+ v12 = v;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ }
+ // printf("quadrangle generated\n");
+ if (v11 != v12){
+ addOrRemove(v11,v12,bedges);
+ addOrRemove(v12,v22,bedges);
+ addOrRemove(v22,v21,bedges);
+ addOrRemove(v21,v11,bedges);
+ blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
+ fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
+ v11->x(),v11->y(),v11->z(),
+ v12->x(),v12->y(),v12->z(),
+ v22->x(),v22->y(),v22->z(),
+ v21->x(),v21->y(),v21->z());
+ }
+ else {
+ addOrRemove(v,v22,bedges);
+ addOrRemove(v22,v21,bedges);
+ addOrRemove(v21,v,bedges);
+ blTris.push_back(new MTriangle(v,v22,v21));
+ fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
+ v->x(),v->y(),v->z(),
+ v22->x(),v22->y(),v22->z(),
+ v21->x(),v21->y(),v21->z());
+ }
+ if(v11->onWhat() == gf)verts.insert(v11);
+ if(v21->onWhat() == gf)verts.insert(v21);
+ if(v12->onWhat() == gf)verts.insert(v12);
+ if(v22->onWhat() == gf)verts.insert(v22);
+ }
+ }
+ }
+ }
+ fprintf(ff2,"};\n");
+ fclose(ff2);
+
+ discreteEdge ne (gf->model(), 444444,0,
+ (*edges.begin())->getEndVertex());
+ std::list<GEdge*> hop;
+ std::set<MEdge,Less_Edge>::iterator it = bedges.begin();
+
+ FILE *ff = fopen ("toto.pos","w");
+ fprintf(ff,"View \" \"{\n");
+ for (; it != bedges.end(); ++it){
+ ne.lines.push_back(new MLine (it->getVertex(0),it->getVertex(1)));
+ fprintf(ff,"SL (%g,%g,%g,%g,%g,%g){1,1};\n",
+ it->getVertex(0)->x(),it->getVertex(0)->y(),it->getVertex(0)->z(),
+ it->getVertex(1)->x(),it->getVertex(1)->y(),it->getVertex(1)->z());
+ }
+ fprintf(ff,"};\n");
+ fclose(ff);
+
+ hop.push_back(&ne);
+
+ deMeshGFace kil_;
+ kil_(gf);
+
+ meshGenerator(gf, 0, 0, true , true, &hop);
+
+ gf->quadrangles = blQuads;
+ gf->triangles.insert(gf->triangles.begin(),blTris.begin(),blTris.end());
+ gf->mesh_vertices.insert(gf->mesh_vertices.begin(),verts.begin(),verts.end());
+}
+
+
// Builds An initial triangular mesh that respects the boundaries of
// the domain, including embedded points and surfaces
-static bool meshGenerator(GFace *gf, int RECUR_ITER,
- bool repairSelfIntersecting1dMesh,
- bool onlyInitialMesh,
- bool debug = true)
+bool meshGenerator(GFace *gf, int RECUR_ITER,
+ bool repairSelfIntersecting1dMesh,
+ bool onlyInitialMesh,
+ bool debug,
+ std::list<GEdge*> *replacement_edges)
{
BDS_GeomEntity CLASS_F(1, 2);
BDS_GeomEntity CLASS_EXTERIOR(1, 3);
std::map<BDS_Point*, MVertex*> recoverMap;
std::map<MVertex*, BDS_Point*> recoverMapInv;
- std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*> edges = replacement_edges ? *replacement_edges : gf->edges();
std::list<int> dir = gf->edgeOrientations();
// replace edges by their compounds
@@ -494,6 +732,8 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
}
all_vertices.clear();
+ // here check if some boundary layer nodes should be added
+
// compute the bounding box in parametric space
SVector3 dd(bbox.max(), bbox.min());
double LC2D = norm(dd);
@@ -578,8 +818,6 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
++ite;
}
- Msg::Debug("Recovering %d mesh Edges (%d not recovered)", edgesToRecover.size(),
- edgesNotRecovered.size());
// effectively recover the medge
ite = edges.begin();
@@ -594,6 +832,9 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
++ite;
}
+ Msg::Debug("Recovering %d mesh Edges (%d not recovered)", edgesToRecover.size(),
+ edgesNotRecovered.size());
+
if(edgesNotRecovered.size()){
std::ostringstream sstream;
for(std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
@@ -601,7 +842,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
sstream << " " << itr->ge->tag();
Msg::Warning(":-( There are %d intersections in the 1D mesh (curves%s)",
edgesNotRecovered.size(), sstream.str().c_str());
- Msg::Warning("8-| Gmsh splits those edges and tries again");
+ if (repairSelfIntersecting1dMesh) Msg::Warning("8-| Gmsh splits those edges and tries again");
if(debug){
char name[245];
@@ -621,6 +862,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
int p1 = itr->p1;
int p2 = itr->p2;
int tag = itr->ge->tag();
+ printf("%d %d %d\n",p1,p2,tag);
_error[3 * I + 0] = p1;
_error[3 * I + 1] = p2;
_error[3 * I + 2] = tag;
@@ -633,7 +875,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
delete m;
if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
return meshGenerator
- (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh, debug);
+ (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh, debug,replacement_edges);
return false;
}
@@ -718,25 +960,27 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
}
// compute characteristic lengths at vertices
- Msg::Debug("Computing mesh size field at mesh vertices %d",
- edgesToRecover.size());
- for(int i = 0; i < doc.numPoints; i++){
- BDS_Point *pp = (BDS_Point*)doc.points[i].data;
- std::map<BDS_Point*, MVertex*>::iterator itv = recoverMap.find(pp);
- if(itv != recoverMap.end()){
- MVertex *here = itv->second;
- GEntity *ge = here->onWhat();
- if(ge->dim() == 0){
- pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
- }
- else if(ge->dim() == 1){
- double u;
- here->getParameter(0, u);
- pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
- }
- else
- pp->lcBGM() = MAX_LC;
- pp->lc() = pp->lcBGM();
+ if (!onlyInitialMesh){
+ Msg::Debug("Computing mesh size field at mesh vertices %d",
+ edgesToRecover.size());
+ for(int i = 0; i < doc.numPoints; i++){
+ BDS_Point *pp = (BDS_Point*)doc.points[i].data;
+ std::map<BDS_Point*, MVertex*>::iterator itv = recoverMap.find(pp);
+ if(itv != recoverMap.end()){
+ MVertex *here = itv->second;
+ GEntity *ge = here->onWhat();
+ if(ge->dim() == 0){
+ pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
+ }
+ else if(ge->dim() == 1){
+ double u;
+ here->getParameter(0, u);
+ pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
+ }
+ else
+ pp->lcBGM() = MAX_LC;
+ pp->lc() = pp->lcBGM();
+ }
}
}
}
@@ -780,8 +1024,8 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
sprintf(name, "surface%d-recovered-param.pos", gf->tag());
outputScalarField(m->triangles, name, 1);
}
-
- {
+
+ if(1){
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end()){
BDS_Face *t = *itt;
@@ -803,7 +1047,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
++itt;
}
}
-
+
if (Msg::GetVerbosity() == 10){
GEdge *ge = new discreteEdge(gf->model(), 1000, 0, 0);
MElementOctree octree(gf->model());
@@ -819,69 +1063,45 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
gf->triangles.clear();
gf->quadrangles.clear();
- int nb_swap;
- //outputScalarField(m->triangles, "beforeswop.pos",1);
- Msg::Debug("Delaunizing the initial mesh");
- delaunayizeBDS(gf, *m, nb_swap);
- //outputScalarField(m->triangles, "afterswop.pos",0);
- Msg::Debug("Starting to add internal points");
+ {
+ int nb_swap;
+ Msg::Debug("Delaunizing the initial mesh");
+ delaunayizeBDS(gf, *m, nb_swap);
+ }
+ Msg::Debug("Starting to add internal points");
// start mesh generation
if(!algoDelaunay2D(gf) && !onlyInitialMesh){
+ if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) {
+ printf("coucou here !!!\n");
+ backgroundMesh::unset();
+ buildBackGroundMesh (gf);
+ }
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true,
&recoverMapInv);
optimizeMeshBDS(gf, *m, 2);
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, false,
&recoverMapInv);
optimizeMeshBDS(gf, *m, 2);
+ if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) {
+ backgroundMesh::unset();
+ }
}
+ /*
computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index,
gf->meshStatistics.longest_edge_length,
gf->meshStatistics.smallest_edge_length,
gf->meshStatistics.nbEdge,
gf->meshStatistics.nbGoodLength);
+ */
gf->meshStatistics.status = GFace::DONE;
// fill the small gmsh structures
- {
- std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
- while (itp != m->points.end()){
- BDS_Point *p = *itp;
- if(recoverMap.find(p) == recoverMap.end()){
- MVertex *v = new MFaceVertex
- (p->X, p->Y, p->Z, gf, m->scalingU * p->u, m->scalingV * p->v);
- recoverMap[p] = v;
- gf->mesh_vertices.push_back(v);
- }
- ++itp;
- }
- }
- {
- std::list<BDS_Face*>::iterator itt = m->triangles.begin();
- while (itt != m->triangles.end()){
- BDS_Face *t = *itt;
- if(!t->deleted){
- BDS_Point *n[4];
- t->getNodes(n);
- MVertex *v1 = recoverMap[n[0]];
- MVertex *v2 = recoverMap[n[1]];
- MVertex *v3 = recoverMap[n[2]];
- if(!n[3]){
- // when a singular point is present, degenerated triangles
- // may be created, for example on a sphere that contains one
- // pole
- if(v1 != v2 && v1 != v3 && v2 != v3)
- gf->triangles.push_back(new MTriangle(v1, v2, v3));
- }
- else{
- MVertex *v4 = recoverMap[n[3]];
- gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
- }
- }
- ++itt;
- }
- }
+ BDS2GMSH ( m, gf, recoverMap);
+ // BOUNDARY LAYER
+ if (!onlyInitialMesh)modifyInitialMeshForTakingIntoAccountBoundaryLayers (gf);
+
// the delaunay algo is based directly on internal gmsh structures
// BDS mesh is passed in order not to recompute local coordinates of
// vertices
@@ -1415,17 +1635,27 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// start mesh generation for periodic face
if(!algoDelaunay2D(gf)){
+ // need for a BGM for cross field
+ if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) {
+ // printf("coucou here !!!\n");
+ // backgroundMesh::unset();
+ // buildBackGroundMesh (gf);
+ }
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true);
optimizeMeshBDS(gf, *m, 2);
refineMeshBDS(gf, *m, -CTX::instance()->mesh.refineSteps, false);
optimizeMeshBDS(gf, *m, 2, &recoverMap);
// compute mesh statistics
+ /*
computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index,
gf->meshStatistics.longest_edge_length,
gf->meshStatistics.smallest_edge_length,
gf->meshStatistics.nbEdge,
- gf->meshStatistics.nbGoodLength);
+ gf->meshStatistics.nbGoodLength);*/
gf->meshStatistics.status = GFace::DONE;
+ if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) {
+ // backgroundMesh::unset();
+ }
}
// fill the small gmsh structures
@@ -1554,6 +1784,9 @@ void meshGFace::operator() (GFace *gf)
}
const char *algo = "Unknown";
+
+
+
switch(CTX::instance()->mesh.algo2d){
case ALGO_2D_MESHADAPT : algo = "MeshAdapt"; break;
case ALGO_2D_FRONTAL : algo = "Frontal"; break;
@@ -1566,6 +1799,11 @@ void meshGFace::operator() (GFace *gf)
break;
}
+ if (!algoDelaunay2D(gf)){
+ algo = "MeshAdapt";
+ }
+
+
Msg::Info("Meshing surface %d (%s, %s)", gf->tag(), gf->getTypeString().c_str(), algo);
// compute loops on the fly (indices indicate start and end points
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 61ca4979b3a8f29b1019e450a03a255a0953f31f..1173d253194928e0ed5221b902c79f3d7526d723 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -28,19 +28,34 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
const double dx = p1->X - p2->X;
const double dy = p1->Y - p2->Y;
const double dz = p1->Z - p2->Z;
+ // FIXME SUPER HACK
+ //const double l = std::max(fabs(dx),std::max(fabs(dy),fabs(dz)));
const double l = sqrt(dx * dx + dy * dy + dz * dz);
return l;
}
+extern double lengthInfniteNorm(const double p[2], const double q[2],
+ const double quadAngle);
+
inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
{
+ // FIXME SUPER HACK
+ /*
+ if (CTX::instance()->mesh.recombineAll || f->meshAttributes.recombine){
+ double quadAngle = backgroundMesh::current()->getAngle (0.5 * (p1->u + p2->u) * SCALINGU, 0.5 * (p1->v + p2->v) * SCALINGV,0);
+ const double a [2] = {p1->u,p1->v};
+ const double b [2] = {p2->u,p2->v};
+ return lengthInfniteNorm(a,b, quadAngle);
+ }
+ */
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
0.5 * (p1->v + p2->v) * SCALINGV));
if (!GP.succeeded())
return computeEdgeLinearLength(p1,p2);
+
const double dx1 = p1->X - GP.x();
const double dy1 = p1->Y - GP.y();
const double dz1 = p1->Z - GP.z();
@@ -143,6 +158,7 @@ inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
double SCALINGU, double SCALINGV)
{
+ // FIXME !!!
if (f->geomType() == GEntity::Plane)
return e->length();
else
@@ -532,7 +548,7 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
mid->u = U;
mid->v = V;
- if (backgroundMesh::current()){
+ if (backgroundMesh::current() && 0){
mid->lc() = mid->lcBGM() =
backgroundMesh::current()->operator()
((coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
@@ -613,7 +629,8 @@ void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
{
- // return;
+ // FIXME SUPER HACK
+ // return;
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while(itp != m.points.end()){
if(m.smooth_point_centroid(*itp, gf,q))
@@ -676,7 +693,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
double t_spl = 0, t_sw = 0,t_col = 0,t_sm = 0;
- const double MINE_ = 0.67, MAXE_ = 1.4;
+ const double MINE_ = 0.6666, MAXE_ = 1.4;
while (1){
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 1c09478da741d8001a79496d4c222b8f908d086b..ffda6e0a222dc01432de3f4bd39e90dabcc13bd8 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -16,10 +16,12 @@
#include "Numeric.h"
#include "STensor3.h"
#include "Context.h"
+#include "meshGFaceBoundaryLayers.h"
+#include "MLine.h"
+#include "MQuadrangle.h"
double LIMIT_ = 0.5 * sqrt(2.0) * 1;
-static const bool _experimental_anisotropic_blues_band_ = false;
int MTri3::radiusNorm = 2;
// This function controls the frontal algorithm
@@ -493,7 +495,7 @@ double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
return s * 0.5;
}
-bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
+bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
std::set<MTri3*, compareTri3Ptr> &allTets,
std::set<MTri3*, compareTri3Ptr> *activeTets,
std::vector<double> &vSizes,
@@ -513,6 +515,7 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
}
else{
recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);
+ // printf("RECURSIVELY FIND A CAVITY %d %d \n",shell.size(),cavity.size());
}
// check that volume is conserved
@@ -543,17 +546,7 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
double LL = Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM;
MTri3 *t4;
- if (_experimental_anisotropic_blues_band_){
- SMetric3 metrBGM = interpolation(vMetricsBGM[t->getVertex(0)->getIndex()],
- vMetricsBGM[t->getVertex(1)->getIndex()],
- vMetricsBGM[t->getVertex(2)->getIndex()],
- ONE_THIRD, ONE_THIRD);
-
- t4 = new MTri3(t, LL,&metrBGM);
- }
- else{
- t4 = new MTri3(t, LL,0,&Us,&Vs,gf);
- }
+ t4 = new MTri3(t, LL,0,&Us,&Vs,gf);
double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
(it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
@@ -563,8 +556,10 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
(it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
const double MID[3] = {0.5*(it->v[0]->x()+it->v[1]->x()),0.5*(it->v[0]->y()+it->v[1]->y()),0.5*(it->v[0]->z()+it->v[1]->z())};
double d3 = sqrt((MID[0] - v->x()) * (MID[0] - v->x()) + (MID[1] - v->y()) * (MID[1] - v->y()) + (MID[2] - v->z()) * (MID[2] - v->z()));
- if (d1 < LL * .25 || d2 < LL * .25 || d3 < LL * .25) onePointIsTooClose = true;
-
+ if ((d1 < LL * .25 || d2 < LL * .25 || d3 < LL * .25) && !force) {
+ // printf("TOO CLOSE %g %g %g %g %g %g\n",d1,d2,d3,LL,lc,lcBGM);
+ onePointIsTooClose = true;
+ }
// if (t4->getRadius () < LIMIT_ / 2) onePointIsTooClose = true;
newTris[k++] = t4;
@@ -576,14 +571,15 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
if (otherSide)
new_cavity.push_back(otherSide);
double ss = fabs(getSurfUV(t4->tri(), Us, Vs));
- if (ss < 1.e-25) ss = 1256172121;
+ if (ss < 1.e-25) ss = 1.e22;
newVolume += ss;
++it;
}
- if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() > 3 &&
+ if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() >= 3 &&
!onePointIsTooClose){
connectTris(new_cavity.begin(), new_cavity.end());
allTets.insert(newTris, newTris + shell.size());
+ // printf("shell.size() = %d triangles.size() = %d \n",shell.size(),allTets.size());
if (activeTets){
for (std::list<MTri3*>::iterator i = new_cavity.begin(); i != new_cavity.end(); ++i){
int active_edge;
@@ -598,6 +594,8 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
}
// The cavity is NOT star shaped
else{
+ // printf("cavity not star shaped %22.15E vs %22.15E\n",oldVolume,newVolume);
+ // printf("shell.size() = %d triangles.size() = %d \n",shell.size(),allTets.size());
for (unsigned int i = 0; i < shell.size(); i++) delete newTris[i];
delete [] newTris;
ittet = cavity.begin();
@@ -717,6 +715,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
MTri3 *ptin = 0;
bool inside = false;
double uv[2];
+ // FIXME !!!
if (MTri3::radiusNorm == 2){
inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
if (inside)ptin = worst;
@@ -736,6 +735,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
else {
ptin = search4Triangle (worst, center, Us, Vs, AllTris);
if (ptin)inside = true;
+ // if (!ptin)printf("strange : %g %g seems to be out of the domain for face %d\n",center[0],center[1],gf->tag());
}
// if (!ptin)ptin = worst;
@@ -766,20 +766,13 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
Vs.push_back(center[1]);
if(!p.succeeded() || !insertVertex
- (gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
+ (false, gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric) ) {
MTriangle *base = worst->tri();
- /*
- Msg::Info("Point %g %g of (%g %g , %g %g , %g %g) cannot be inserted",
- center[0], center[1],
- Us[base->getVertex(0)->getIndex()],
- Vs[base->getVertex(0)->getIndex()],
- Us[base->getVertex(1)->getIndex()],
- Vs[base->getVertex(1)->getIndex()],
- Us[base->getVertex(2)->getIndex()],
- Vs[base->getVertex(2)->getIndex()]);
- */
+
+ // Msg::Info("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
+
AllTris.erase(it);
worst->forceRadius(-1);
AllTris.insert(worst);
@@ -787,13 +780,14 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
return false;
}
else {
+ // printf("ouf ! %d\n",AllTris.size());
gf->mesh_vertices.push_back(v);
return true;
}
}
else {
MTriangle *base = worst->tri();
- /*
+ /*
Msg::Info("Point %g %g is outside (%g %g , %g %g , %g %g)",
center[0], center[1],
Us[base->getVertex(0)->getIndex()],
@@ -810,6 +804,137 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
+// initially adds some boundary layer points onto the mesh
+void addBoundaryLayerPoints (GFace *gf){
+
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ std::vector<double> vSizes, vSizesBGM, Us, Vs;
+ std::vector<SMetric3> vMetricsBGM;
+
+ BoundaryLayerColumns *_columns = buidAdditionalPoints2D (gf, M_PI/6.);
+
+ buildMeshGenerationDataStructures
+ (gf, AllTris, vSizes, vSizesBGM, vMetricsBGM, Us, Vs);
+
+ for (BoundaryLayerColumns::iter it = _columns->begin(); it != _columns->end(); ++it){
+ for (int i=0;i<it->second._column.size();++i){
+ while(1){
+ MTri3 *worst = *AllTris.begin();
+ if (worst->isDeleted()){
+ delete worst->tri();
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ }
+ else {
+ break;
+ }
+ }
+ MVertex *v = it->second._column[i];
+ SMetric3 mmm = it->second._metrics[i];
+
+ double center[2],metric[3] = {1.,0,1.};
+ v->getParameter(0,center[0]);
+ v->getParameter(1,center[1]);
+ buildMeshMetric (gf,center,mmm,metric);
+ v->setIndex(Us.size());
+ vSizesBGM.push_back(vSizesBGM[((it->first))->getIndex()]);
+ vSizes.push_back(vSizes[((it->first))->getIndex()]);
+ Us.push_back(center[0]);
+ Vs.push_back(center[1]);
+ MTri3 *ptin = search4Triangle (*AllTris.begin(), center, Us, Vs, AllTris);
+ if (ptin){
+ bool ok = insertVertex
+ (true, gf, v, center, ptin, AllTris,0, vSizes, vSizesBGM,vMetricsBGM,
+ Us, Vs, metric);
+ if (ok)gf->mesh_vertices.push_back(v);
+ else printf("failed !\n");
+ }
+ else {
+ printf("boundary layer point outside the domain ;-)\n");
+ }
+ }
+ }
+
+ // Create all edges of the triangular mesh
+ std::set<MEdge,Less_Edge> allEdges;
+ {
+ std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
+ for (;it != AllTris.end();++it ){
+ MTri3 *worst = *it;
+ if (!worst->isDeleted()){
+ allEdges.insert(MEdge((*it)->tri()->getVertex(0),(*it)->tri()->getVertex(1)));
+ allEdges.insert(MEdge((*it)->tri()->getVertex(0),(*it)->tri()->getVertex(2)));
+ allEdges.insert(MEdge((*it)->tri()->getVertex(1),(*it)->tri()->getVertex(2)));
+ }
+ }
+ }
+ // printf("%d edges generated\n",allEdges.size());
+
+ // look for existing cavities and fill them up with quads
+ std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*>::iterator ite = edges.begin();
+ while(ite != edges.end()){
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ int nbCol1 = _columns->getNbColumns(v1);
+ int nbCol2 = _columns->getNbColumns(v2);
+ if (nbCol1 >= 0 && nbCol2 >= 0){
+ const BoundaryLayerData & c1 = _columns->getColumn(v1,MEdge(v1,v2));
+ const BoundaryLayerData & c2 = _columns->getColumn(v2,MEdge(v1,v2));
+ int N = std::min(c1._column.size(),c2._column.size());
+
+ // Look for the highest edge on the column
+ int highestHorizontal = -1;
+ int highestVertical = -1;
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v21,*v22;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ if (l == 0){
+ v11 = v1;
+ v12 = v2;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ }
+ if (allEdges.find(MEdge(v11,v21)) != allEdges.end() && allEdges.find(MEdge(v12,v22)) != allEdges.end()){
+ highestVertical = l;
+ }
+ else {
+ break;
+ }
+
+ if (allEdges.find(MEdge(v21,v22))!= allEdges.end()){
+ highestHorizontal = l;
+ }
+ }
+ // printf("%d vs %d %d\n",N,highestHorizontal,highestVertical);
+ N = std::min(highestHorizontal,highestVertical);
+ //N--;
+ for (int l=0;l <= N ;++l){
+ MVertex *v11,*v12,*v21,*v22;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ if (l == 0){
+ v11 = v1;
+ v12 = v2;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ }
+ // printf("quadrangle generated\n");
+ gf->quadrangles.push_back(new MQuadrangle(v11,v12,v22,v21));
+ }
+ }
+ }
+ ++ite;
+ }
+ transferDataStructure(gf, AllTris, Us, Vs);
+}
+
void bowyerWatson(GFace *gf)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
@@ -856,15 +981,7 @@ void bowyerWatson(GFace *gf)
(Vs[base->getVertex(0)->getIndex()] +
Vs[base->getVertex(1)->getIndex()] +
Vs[base->getVertex(2)->getIndex()]) / 3.};
- if (_experimental_anisotropic_blues_band_){
- SMetric3 m = interpolation (vMetricsBGM[base->getVertex(0)->getIndex()],
- vMetricsBGM[base->getVertex(1)->getIndex()],
- vMetricsBGM[base->getVertex(2)->getIndex()],
- pa[0],pa[1]);
- buildMetric(gf, pa, m, metric);
- }
- else
- buildMetric(gf, pa, metric);
+ buildMetric(gf, pa, metric);
circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
insertAPoint(gf, AllTris.begin(), center, metric, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris);
@@ -881,7 +998,7 @@ void bowyerWatson(GFace *gf)
The point isertion is done on the Voronoi Edges
*/
-static double lengthInfniteNorm(const double p[2], const double q[2],
+double lengthInfniteNorm(const double p[2], const double q[2],
const double quadAngle){
double xp = p[0] * cos(quadAngle) + p[1] * sin(quadAngle);
double yp = -p[0] * sin(quadAngle) + p[1] * cos(quadAngle);
@@ -1078,6 +1195,7 @@ void bowyerWatsonFrontal(GFace *gf)
// sprintf(name,"frontal%d-param.pos", gf->tag());
// _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
+ quadsToTriangles(gf,10000);
}
void optimalPointFrontalQuad (GFace *gf,
@@ -1163,14 +1281,134 @@ void optimalPointFrontalQuad (GFace *gf,
}
}
+void optimalPointFrontalQuadAniso (GFace *gf,
+ MTri3* worst,
+ int active_edge,
+ std::vector<double> &Us,
+ std::vector<double> &Vs,
+ std::vector<double> &vSizes,
+ std::vector<double> &vSizesBGM,
+ double newPoint[2],
+ double metric[3]){
+ MTriangle *base = worst->tri();
+ int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
+ int ip2 = active_edge;
+ int ip3 = (active_edge+1)%3;
+
+ double P[2] = {Us[base->getVertex(ip1)->getIndex()],
+ Vs[base->getVertex(ip1)->getIndex()]};
+ double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
+ Vs[base->getVertex(ip2)->getIndex()]};
+ double O[2] = {Us[base->getVertex(ip3)->getIndex()],
+ Vs[base->getVertex(ip3)->getIndex()]};
+ double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
+
+ // compute background mesh data
+ double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0);
+ double center[2];
+ circumCenterInfinite (base, quadAngle,Us,Vs,center);
+
+ // rotate the points with respect to the angle
+ double XP1 = 0.5*(Q[0] - P[0]);
+ double YP1 = 0.5*(Q[1] - P[1]);
-void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
-{
+ double DX = 0.5*(Q[0] + P[0]);
+ double DY = 0.5*(Q[0] + P[0]);
- std::set<MTri3*,compareTri3Ptr> AllTris;
- std::set<MTri3*,compareTri3Ptr> ActiveTris;
- std::vector<double> vSizes, vSizesBGM, Us, Vs;
- std::vector<SMetric3> vMetricsBGM;
+ double xp = XP1 * cos(quadAngle) + YP1 * sin(quadAngle);
+ double yp = -XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
+
+ double X4 = O[0] -DX;
+ double Y4 = O[1] -DY;
+
+ double x4 = X4 * cos(quadAngle) + Y4 * sin(quadAngle);
+ double y4 = -X4 * sin(quadAngle) + Y4 * cos(quadAngle);
+
+ // ensure xp > yp
+ bool exchange = false;
+ if (fabs(xp) < fabs(yp)){
+ double temp = xp;
+ xp = yp;
+ yp = temp;
+ temp = x4;
+ x4 = y4;
+ y4 = temp;
+ exchange = true;
+ }
+
+ buildMetric(gf, midpoint, metric);
+ double RATIO = 1./pow(metric[0]*metric[2]-metric[1]*metric[1],0.25);
+ const double rhoM1 = 0.5 * RATIO *
+ (vSizes[base->getVertex(ip1)->getIndex()] +
+ vSizes[base->getVertex(ip2)->getIndex()] );
+ const double rhoM2 = 0.5 * RATIO *
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ vSizesBGM[base->getVertex(ip2)->getIndex()] );// * RATIO;
+ const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
+
+
+ double npx,npy;
+
+ const double L_edge = lengthInfniteNorm(P,Q,quadAngle);
+
+ const double XP[2]={fabs(xp),fabs(yp)};
+ if (xp*yp > 0){
+ double xe[2] = { rhoM - fabs(xp) + fabs(yp), rhoM };
+ double xc[2] = {0.5*(x4 - fabs(xp)), 0.5*(x4 + fabs(xp)) - fabs(yp)};
+ double xl[2] = {rhoM, rhoM+fabs(xp)-fabs(yp)};
+ const double R_Active = lengthInfniteNorm(xc,XP,0.0); // alerady rotated !
+ const double L_ec = lengthInfniteNorm(xe,xc,0.0); // alerady rotated !
+ const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
+ double t = ( L_edge - rhoM)/(L_edge - R_Active);
+ t = std::max(1.,t);
+ t = std::min(-L_el/L_ec,t);
+ npx = xe[0] + t*(xc[0]-xe[0]);
+ npy = xe[1] + t*(xc[1]-xe[1]);
+ }
+ else {
+ double xe[2] = { rhoM - xp + yp, rhoM };
+ double xc[2] = {0.5*(x4 - xp), 0.5*(x4 + xp) - yp};
+ double xl[2] = {rhoM, rhoM+xp-yp};
+ const double R_Active = lengthInfniteNorm(xc,XP,0.0); // alerady rotated !
+ const double L_ec = lengthInfniteNorm(xe,xc,0.0); // alerady rotated !
+ const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
+ const double XP[2]={xp,yp};
+ double t = ( L_edge - rhoM)/(L_edge - R_Active);
+ t = std::max(1.,t);
+ t = std::min(-L_el/L_ec,t);
+ npx = xe[0] + t*(xc[0]-xe[0]);
+ npy = xe[1] + t*(xc[1]-xe[1]);
+ }
+
+ /*
+ if (xp*yp > 0){
+ npx = - fabs(xp)*factor;
+ npy = fabs(xp)*(1.+factor) - fabs(yp);
+ }
+ else {
+ npx = fabs(xp) * factor;
+ npy = (1.+factor)*fabs(xp) - fabs(yp);
+ }
+ */
+ if (exchange){
+ double temp = npx;
+ npx = npy;
+ npy = temp;
+ }
+
+
+ newPoint[0] = midpoint[0] + cos(quadAngle) * npx - sin(quadAngle) * npy;
+ newPoint[1] = midpoint[1] + sin(quadAngle) * npx + cos(quadAngle) * npy;
+
+ if ((midpoint[0] - newPoint[0])*(midpoint[0] - O[0]) +
+ (midpoint[1] - newPoint[1])*(midpoint[1] - O[1]) < 0){
+ newPoint[0] = midpoint[0] - cos(quadAngle) * npx + sin(quadAngle) * npy;
+ newPoint[1] = midpoint[1] - sin(quadAngle) * npx - cos(quadAngle) * npy;
+ }
+}
+
+
+void buildBackGroundMesh (GFace *gf) {
if (!backgroundMesh::current()) {
std::vector<MTriangle*> TR;
@@ -1179,10 +1417,16 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
gf->triangles[i]->getVertex(1),
gf->triangles[i]->getVertex(2)));
}
+ // avoid computing curvatures on the fly : only on the
+ // BGM computes once curvatures at each node
+ // Disable curvature control
int CurvControl = CTX::instance()->mesh.lcFromCurvature;
CTX::instance()->mesh.lcFromCurvature = 0;
+ // Do a background mesh
bowyerWatson(gf);
+ // Re-enable curv control if asked
CTX::instance()->mesh.lcFromCurvature = CurvControl;
+ // apply this to the BGM
backgroundMesh::set(gf);
char name[256];
if (CTX::instance()->mesh.saveAll){
@@ -1191,12 +1435,22 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
sprintf(name,"cross-%d.pos",gf->tag());
backgroundMesh::current()->print(name,gf,1);
}
- // FIXME DELETE CURRENT MESH
- gf->triangles = TR;
- }
+ }
+}
+
+void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
+{
+
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ std::set<MTri3*,compareTri3Ptr> ActiveTris;
+ std::vector<double> vSizes, vSizesBGM, Us, Vs;
+ std::vector<SMetric3> vMetricsBGM;
+
+ buildBackGroundMesh (gf);
if (quad){
LIMIT_ = sqrt(2.) * .99;
+ //LIMIT_ = 4./3.;//sqrt(2.) * .99;
MTri3::radiusNorm =-1;
}
@@ -1234,6 +1488,8 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
std::set<MTri3*,compareTri3Ptr> ActiveTrisNotInFront;
+ // printf("%d active triangles\n",ActiveTris.size());
+
while (1){
if (!ActiveTris.size())break;
@@ -1252,6 +1508,8 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
// for (active_edge = 0 ; active_edge < 0 ; active_edge ++){
// if (active_edges[active_edge])break;
// }
+ // Msg::Info("%7d points created -- Worst tri infinite radius is %8.3f -- front size %6d",
+ // vSizes.size(), worst->getRadius(),_front.size());
if(ITER++ % 5000 == 0)
Msg::Debug("%7d points created -- Worst tri infinite radius is %8.3f -- front size %6d",
vSizes.size(), worst->getRadius(),_front.size());
@@ -1262,12 +1520,13 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
else optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
+ // printf("start INSERT A POINT %g %g \n",newPoint[0],newPoint[1]);
insertAPoint(gf, AllTris.end(), newPoint, 0, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris, &ActiveTris, worst);
// else if (!worst->isDeleted() && worst->getRadius() > LIMIT_){
// ActiveTrisNotInFront.insert(worst);
// }
-
+ // printf("-----------------> size %d\n",AllTris.size());
/*
if(ITER % 1== 0){
@@ -1295,11 +1554,11 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
if (!ActiveTris.size())break;
}
- // char name[245];
+ char name[245];
// sprintf(name,"frontal%d-real.pos", gf->tag());
// _printTris (name, AllTris, Us, Vs,false);
- // sprintf(name,"frontal%d-param.pos", gf->tag());
- // _printTris (name, AllTris, Us, Vs,true);
+ // sprintf(name,"frontal%d-param.pos", gf->tag());
+ // _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
MTri3::radiusNorm = 2;
LIMIT_ = 0.5 * sqrt(2.0) * 1;
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 0890d9888c29358cd8a79d92f88da6472a523dda..ce4f189839760d02d98d1b418a7ab787ee69ee1c 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -97,6 +97,8 @@ void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
void bowyerWatson(GFace *gf);
void bowyerWatsonFrontal(GFace *gf);
void bowyerWatsonFrontalLayers(GFace *gf, bool quad);
+void buildBackGroundMesh (GFace *gf);
+void addBoundaryLayerPoints (GFace *gf);
struct edgeXface
{
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e69abe59048b625404201aeae8e7f5cd8a1b4e56..6c32f154564b24a5d70fe3cf8fce11884af0f296 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -846,6 +846,8 @@ static void _relocateVertexConstrained (GFace *gf,
const std::vector<MElement*> <){
if( ver->onWhat()->dim() == 2){
+ MFaceVertex *fv = dynamic_cast<MFaceVertex*>(ver);
+ if (fv->bl_data)return;
std::map<MVertex*,p1p2p3> ring;
double initu,initv;
@@ -923,8 +925,13 @@ static void _relocateVertexConstrained (GFace *gf,
static void _relocateVertex (GFace *gf,
MVertex *ver,
const std::vector<MElement*> <){
+
double R; SPoint3 c; bool isSphere = gf->isSphere(R,c);
if( ver->onWhat()->dim() == 2){
+
+ MFaceVertex *fv = dynamic_cast<MFaceVertex*>(ver);
+ if (fv->bl_data)return;
+
double initu,initv;
ver->getParameter(0, initu);
ver->getParameter(1, initv);
@@ -957,17 +964,18 @@ static void _relocateVertex (GFace *gf,
}
bool success = false;
double factor = 1.0;
+ SPoint2 newp;
for (int i=0;i<10;i++){
- SPoint2 newp = after + before*(1.-factor);
+ newp = after + before*(1.-factor);
success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,before,newp);
if (success)break;
factor *= 0.5;
}
- if (success){ // ne devrait-on pas utiliser newp ici au lieu de after ??? (vu la boucle precedente)
- ver->setParameter(0, after.x());
- ver->setParameter(1, after.y());
- GPoint pt = gf->point(after);
+ if (success){
+ ver->setParameter(0, newp.x());
+ ver->setParameter(1, newp.y());
+ GPoint pt = gf->point(newp);
if(pt.succeeded()){
ver->x() = pt.x();
ver->y() = pt.y();
@@ -1876,8 +1884,6 @@ void recombineIntoQuads(GFace *gf,
if (saveAll) gf->model()->writeMSH("before.msh");
int success = _recombineIntoQuads(gf, 0);
- // gf->addLayersOfQuads(1, 0);
-
if (saveAll) gf->model()->writeMSH("raw.msh");
if(haveParam && nodeRepositioning)
laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
@@ -1895,7 +1901,7 @@ void recombineIntoQuads(GFace *gf,
if(x && saveAll){ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);}
int y = removeDiamonds(gf);
if(y && saveAll){ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME); }
- laplaceSmoothing(gf);
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
int z = 0; //edgeSwapQuadsForBetterQuality(gf);
if(z && saveAll){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
if (!(x+y+z)) break;
@@ -1922,7 +1928,7 @@ void recombineIntoQuads(GFace *gf,
Msg::Info("Simple recombination algorithm completed (%g s)", t2 - t1);
}
-void quadsToTriangles(GFace *gf, double minqual = -10000.)
+void quadsToTriangles(GFace *gf, double minqual)
{
std::vector<MQuadrangle*> qds;
for (int i=0;i<gf->quadrangles.size();i++){
@@ -1968,7 +1974,7 @@ void recombineIntoQuadsIterative(GFace *gf)
return;
int COUNT = 0;
while (1){
- quadsToTriangles(gf);
+ quadsToTriangles(gf,-10000);
{char NAME[245];sprintf(NAME,"iterT%d.msh",COUNT); gf->model()->writeMSH(NAME);}
std::set<MTri3*,compareTri3Ptr> AllTris;
std::vector<double> vSizes, vSizesBGM, Us, Vs;
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 55b2d34e371aef4b2dcfe9d4c7a66723f4a8af78..5aec057914fe62794aa8407e0835a6409ccfd038 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -98,6 +98,7 @@ void recombineIntoQuads(GFace *gf,
bool topologicalOpti = true,
bool nodeRepositioning = true);
void recombineIntoQuadsIterative(GFace *gf);
+void quadsToTriangles(GFace *gf, double minqual);
struct swapquad{
int v[4];
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index ea84394a4b6b90f6323db0f9dd5e52fa830257d3..57717553ddf46e3d686ef3688a77268594e15039 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -987,8 +987,12 @@ void assignPartitionBoundary(GModel *model, MFace &me,
}
else ppe = *it;
// to finish !!
- ppe->triangles.push_back(new MTriangle (me.getVertex(0),me.getVertex(1),
- me.getVertex(2)));
+ if (me.getNumVertices() == 3)
+ ppe->triangles.push_back(new MTriangle (me.getVertex(0),me.getVertex(1),
+ me.getVertex(2)));
+ else
+ ppe->quadrangles.push_back(new MQuadrangle (me.getVertex(0),me.getVertex(1),
+ me.getVertex(2),me.getVertex(3)));
}
void assignPartitionBoundary(GModel *model,
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 18fe5f57aa5f5e87c7359fbe14ff4a97c12ae9b5..e39cd9e5d8aa2849c659d69e9317841324a9d9d9 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -323,8 +323,14 @@ double mesh_functional_distorsion_pN(MElement *t)
fullVector<double> Bi( jac->matrixLag2Bez.size1() );
jac->matrixLag2Bez.mult(Ji,Bi);
-
- // jac->matrixLag2Bez.print("Lag2Bez");
+ /*
+ jac->matrixLag2Bez.print("Lag2Bez");
+
+ jac->points.print("Points");
+ t->getFunctionSpace(t->getPolynomialOrder())->points.print("lagrangianNodes");
+ t->getFunctionSpace(t->getPolynomialOrder())->monomials.print("Monomials");
+ t->getFunctionSpace(t->getPolynomialOrder())->coefficients.print("Coefficients");
+ */
return *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
}
@@ -337,10 +343,10 @@ double qmDistorsionOfMapping (MTriangle *e)
else if (e->getPolynomialOrder() == 2) {
// const double exact = mesh_functional_distorsion_p2_exact(e);
const double bezier= mesh_functional_distorsion_pN(e);
- if (bezier < .1){
- const double bezier_refined= mesh_functional_distorsion_p2_bezier_refined(e);
- return bezier_refined;
- }
+ // if (bezier < .1){
+ // const double bezier_refined= mesh_functional_distorsion_p2_bezier_refined(e);
+ // return bezier_refined;
+ // }
/*
if (exact < .99 || exact > 1.01){
FILE *f = fopen ("statistics.dat","a");
@@ -527,7 +533,7 @@ double qmQuadrangleAngles (MQuadrangle *e) {
double c;
prosca(v1,v2,&c);
- printf("Youhou %g %g\n",c,acos(c)*180/M_PI);
+ // printf("Youhou %g %g\n",c,acos(c)*180/M_PI);
double x = fabs(acos(c))-M_PI/2;
double quality = (atan(a*(x+M_PI/4)) + atan(a*(2*M_PI/4 - (x+M_PI/4))))/den;
worst_quality = std::min(worst_quality, quality);
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 76110a308b8488b3ea62723fbf57ea21c9d356d2..1285646b83c413981754002955ca77e2e0ca444d 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -125,12 +125,12 @@ class fullVector
// printing and file treatment
void print(const char *name="") const
{
- printf("Printing vector %s:\n", name);
- printf(" ");
+ printf("double %s[%d]=\n", name,size());
+ printf("{ ");
for(int I = 0; I < size(); I++){
printf("%12.5E ", (*this)(I));
}
- printf("\n");
+ printf("};\n");
}
void binarySave (FILE *f) const
{
@@ -493,16 +493,19 @@ class fullMatrix
void print(const std::string name = "", const std::string format = "%12.5E ") const
{
- printf("Printing matrix %s:\n", name.c_str());
int ni = size1();
int nj = size2();
+ printf("double %s [ %d ][ %d ]= { \n", name.c_str(),ni,nj);
for(int I = 0; I < ni; I++){
- printf(" ");
+ printf("{ ");
for(int J = 0; J < nj; J++){
printf(format.c_str(), (*this)(I, J));
+ if (J!=nj-1)printf(",");
}
- printf("\n");
+ if (I!=ni-1)printf("},\n");
+ else printf("}\n");
}
+ printf("};\n");
}
// specific functions for dgshell
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index ae03f491d401334d85e947404be575961b0be502..3557184ccd3a0384577d4be076856068a1de398c 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -50,9 +50,9 @@ class helmholtzTerm : public femTerm<scalar> {
MElement *e = se->getMeshElement();
// compute integration rule
- const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
- 2 * (e->getPolynomialOrder() - 1);
- // const int integrationOrder = 2 * e->getPolynomialOrder() + 1;
+ // const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
+ //2 * (e->getPolynomialOrder() - 1);
+ const int integrationOrder = 2 * e->getPolynomialOrder() + 1;
int npts; IntPt *GP;
e->getIntegrationPoints(integrationOrder, &npts, &GP);
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index 41be12dc9bb454ffecbbb04f8c6042a61865eadc..b61af396b53d9e75a2fc658bddd30b61687f766c 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -78,5 +78,5 @@ Physical Line(25) = {12,13,14,15,16,17,18,19,20,10};
Physical Surface(26) = {24};
Physical Line(27) = {10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 14, 13, 12, 11};
Physical Line(28) = {17, 16, 20, 19, 18, 15};
-//Recombine Surface {24, 22};
-//Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
+Recombine Surface {24, 22};
+Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index 6d1b45d0d4653271f039e6b3dda161886ff8126d..432e233608e2c0bcb2e444669d82104cfdf70608 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -227,3 +227,14 @@ Plane Surface(11) = {9,10};
//Physical Surface(11)={11};
//Point(9999) = {0.6,0,0,1};
+Field[2] = BoundaryLayer;
+Field[2].NodesList = {1};
+//Field[2].EdgesList = {1,2,3,4};
+Field[2].EdgesList = {1,2,3,4};
+Field[2].hfar = 1.5;
+Field[2].hwall_n = 0.0001;
+Field[2].hwall_t = 0.03;
+Field[2].ratio = 1.3;
+Field[2].thickness = .05;
+Background Field = 2;
+
diff --git a/benchmarks/2d/wing-splines.geo b/benchmarks/2d/wing-splines.geo
index 53c3aa4b630c4d1d28c4365c07cf72f930df2d5b..1bb082ea43c2a4f79d8b56ec895b13fe1cacd5ca 100644
--- a/benchmarks/2d/wing-splines.geo
+++ b/benchmarks/2d/wing-splines.geo
@@ -1,5 +1,6 @@
+Mesh.CharacteristicLengthExtendFromBoundary= 0;
-scale = .01 ;
+scale = 1 ;
lc_wing = 0.05 * scale ;
lc_box = 30 * scale ;
@@ -604,5 +605,16 @@ Circle(408) = {4355,4351,4354};
Circle(409) = {4354,4351,4353};
Circle(410) = {4353,4351,4352};
*/
-Recombine Surface {406};
-Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
+//Recombine Surface {406};
+//Mesh.RecombinationAlgorithm=1;
+
+Field[2] = BoundaryLayer;
+Field[2].EdgesList = {12,1,2,3,4,5,6,11,10};
+Field[2].NodesList = {3847,3846,4006,4007,4008,4222,4221};
+Field[2].hfar = 20.1;
+Field[2].hwall_n = 1e-04;
+Field[2].hwall_t = 14.e-3;
+Field[2].ratio = 1.2;
+//Field[2].thickness = .01;
+Field[2].thickness = .0075;
+Background Field = 2;
diff --git a/benchmarks/step/linkrods.geo b/benchmarks/step/linkrods.geo
index 0ec28e0765c668d8cf31ef563d1e1e40cbf8b783..4770c2f0de162c91d70713d11749e44a61c0907b 100644
--- a/benchmarks/step/linkrods.geo
+++ b/benchmarks/step/linkrods.geo
@@ -1,12 +1,13 @@
-Mesh.CharacteristicLengthFactor=0.9;
+Mesh.CharacteristicLengthFactor=.4;
-Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
-Mesh.CharacteristicLengthMin = 0.05; //-clmin
+//Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
+Mesh.CharacteristicLengthMin = 0.25; //-clmin
Mesh.CharacteristicLengthMax = 4.0; //-clmax
-Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
+//Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
Mesh.MinimumCirclePoints=15; //default=7
-Mesh.CharacteristicLengthExtendFromBoundary=0;
+//Mesh.CharacteristicLengthExtendFromBoundary=0;
Merge "linkrods.step";
+
diff --git a/benchmarks/stl/mobilette.geo b/benchmarks/stl/mobilette.geo
index b5242b188bdd69fea7ec7dfb159b2e8ad2425335..e754da25569b598f8603b15ac76e945ae585e4e7 100644
--- a/benchmarks/stl/mobilette.geo
+++ b/benchmarks/stl/mobilette.geo
@@ -1,4 +1,4 @@
-//Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg)
+Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg)
Mesh.CharacteristicLengthMin=1.5;
Mesh.CharacteristicLengthMax=2.5;
Mesh.RemeshAlgorithm=1;
@@ -27,8 +27,8 @@ EndFor
Surface Loop(1) = {s : s + #ss[]-1};
Volume(1) = {1};
-//Mesh.RecombineAll=1;
-//Mesh.RecombinationAlgorithm=1;
+Mesh.RecombineAll=1;
+Mesh.RecombinationAlgorithm=1;
Physical Surface(1) = {s : s + #ss[]-1};
Physical Volume(1) = 1;