From 2462a6b2be8ebd97297eb84a3218ebb1419a4867 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 30 Apr 2014 09:42:00 +0000
Subject: [PATCH] added control on geometrical accuracy of high order elements
coherence now removes degenerate entities (curves, surfaces & volumes). Not
toroughly tested !!!
---
Common/Options.cpp | 3 +-
Common/StringUtils.cpp | 8 +-
Geo/CMakeLists.txt | 3 +-
Geo/GEdge.cpp | 60 +-
Geo/GEdge.h | 5 +-
Geo/GEdgeLoop.cpp | 4 +-
Geo/GModelIO_MSH2.cpp | 1 +
Geo/GModelIO_OCC.cpp | 35 ++
Geo/GModelIO_OCC.h | 1 +
Geo/Geo.cpp | 167 +++++-
Geo/Geo.h | 5 +
Geo/MElement.cpp | 16 +
Geo/MElement.h | 1 +
Geo/MPrism.cpp | 33 ++
Geo/MPrism.h | 1 +
Geo/MQuadrangle.cpp | 1 -
Geo/MVertexPositionSet.h | 1 +
Geo/boundaryLayersData.cpp | 110 ++--
Geo/gmshEdge.cpp | 9 +
Geo/gmshEdge.h | 1 +
Geo/gmshEdgeDiscretize.cpp | 196 +++++++
Mesh/Field.cpp | 2 +-
Mesh/HighOrder.cpp | 14 -
Mesh/HighOrder.h | 5 -
Mesh/meshGFace.cpp | 2 +-
Mesh/meshGRegion.cpp | 236 +++++---
Mesh/meshGRegion.h | 4 +-
Mesh/periodical.cpp | 5 +
Numeric/HilbertCurve.cpp | 4 +-
Plugin/CMakeLists.txt | 1 +
Plugin/PluginManager.cpp | 3 +
benchmarks/boundaryLayers/lann_ver1.geo | 1 +
contrib/HighOrderMeshOptimizer/CMakeLists.txt | 1 +
contrib/HighOrderMeshOptimizer/OptHOM.cpp | 46 +-
contrib/HighOrderMeshOptimizer/OptHOM.h | 13 +-
.../OptHomIntegralBoundaryDist.cpp | 526 ++++++++++++++++++
.../OptHomIntegralBoundaryDist.h | 56 ++
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp | 62 +++
contrib/HighOrderMeshOptimizer/OptHomMesh.h | 1 +
contrib/HighOrderMeshOptimizer/OptHomRun.cpp | 89 ++-
contrib/HighOrderMeshOptimizer/OptHomRun.h | 9 +-
41 files changed, 1513 insertions(+), 228 deletions(-)
create mode 100644 Geo/gmshEdgeDiscretize.cpp
create mode 100644 contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp
create mode 100644 contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 7b39a2c11d..6e6d7fca73 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -4418,8 +4418,9 @@ double opt_geometry_occ_fix_small_faces(OPT_ARGS_NUM)
double opt_geometry_occ_sew_faces(OPT_ARGS_NUM)
{
- if(action & GMSH_SET)
+ if(action & GMSH_SET){
CTX::instance()->geom.occSewFaces = val ? 1 : 0;
+ }
#if defined(HAVE_FLTK)
if(FlGui::available() && (action & GMSH_GUI)) {
FlGui::instance()->options->geo.butt[13]->value
diff --git a/Common/StringUtils.cpp b/Common/StringUtils.cpp
index f609237f6f..788cb4238f 100644
--- a/Common/StringUtils.cpp
+++ b/Common/StringUtils.cpp
@@ -96,13 +96,17 @@ std::vector<std::string> SplitFileName(const std::string &fileName)
int islash = (int)fileName.find_last_of("/\\");
if(idot == (int)std::string::npos) idot = -1;
if(islash == (int)std::string::npos) islash = -1;
- std::vector<std::string> s(3);
+ std::string s[3];
if(idot > 0)
s[2] = fileName.substr(idot);
if(islash > 0)
s[0] = fileName.substr(0, islash + 1);
s[1] = fileName.substr(s[0].size(), fileName.size() - s[0].size() - s[2].size());
- return s;
+ std::vector<std::string> ss;
+ ss.push_back(s[0]);
+ ss.push_back(s[1]);
+ ss.push_back(s[2]);
+ return ss;
}
std::string GetFileNameWithoutPath(const std::string &fileName)
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 2c56e9b75c..ed90152c76 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -5,6 +5,7 @@
set(SRC
boundaryLayersData.cpp
+ closestPoint.cpp
intersectCurveSurface.cpp
GEntity.cpp STensor3.cpp
GVertex.cpp GEdge.cpp GFace.cpp GRegion.cpp
@@ -12,7 +13,7 @@ set(SRC
GRegionCompound.cpp GRbf.cpp
gmshVertex.cpp gmshEdge.cpp gmshFace.cpp gmshRegion.cpp
gmshSurface.cpp
- OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp
+ OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp
discreteEdge.cpp discreteFace.cpp discreteRegion.cpp
fourierEdge.cpp fourierFace.cpp fourierProjectionFace.cpp
ACISVertex.cpp
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index ceffd557f2..66dd62cc9c 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -14,9 +14,10 @@
#include "MLine.h"
#include "GaussLegendre1D.h"
#include "Context.h"
+#include "closestPoint.h"
GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
- : GEntity(model, tag), _tooSmall(false), v0(_v0), v1(_v1), compound(0)
+ : GEntity(model, tag), _tooSmall(false), v0(_v0), v1(_v1), compound(0),_cp(0)
{
if(v0) v0->addEdge(this);
if(v1 && v1 != v0) v1->addEdge(this);
@@ -28,7 +29,7 @@ GEdge::~GEdge()
{
if(v0) v0->delEdge(this);
if(v1 && v1 != v0) v1->delEdge(this);
-
+ if (_cp)delete _cp;
deleteMesh();
}
@@ -376,49 +377,6 @@ GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
return point(t);
}
-bool GEdge::computeDistanceFromMeshToGeometry (double &d2, double &dmax)
-{
- d2 = 0.0; dmax = 0.0;
- if (geomType() == Line) return true;
- if (!lines.size())return false;
- IntPt *pts;
- int npts;
- lines[0]->getIntegrationPoints(2*lines[0]->getPolynomialOrder(), &npts, &pts);
-
- for (unsigned int i = 0; i < lines.size(); i++){
- MLine *l = lines[i];
- double t[256];
-
- for (int j=0; j< l->getNumVertices();j++){
- MVertex *v = l->getVertex(j);
- if (v->onWhat() == getBeginVertex()){
- t[j] = getLowerBound();
- }
- else if (v->onWhat() == getEndVertex()){
- t[j] = getUpperBound();
- }
- else {
- v->getParameter(0,t[j]);
- }
- }
- for (int j=0;j<npts;j++){
- SPoint3 p;
- l->pnt(pts[j].pt[0],0,0,p);
- double tinit = l->interpolate(t,pts[j].pt[0],0,0);
- GPoint pc = closestPoint(p, tinit);
- if (!pc.succeeded())continue;
- double dsq =
- (pc.x()-p.x())*(pc.x()-p.x()) +
- (pc.y()-p.y())*(pc.y()-p.y()) +
- (pc.z()-p.z())*(pc.z()-p.z());
- d2 += pts[i].weight * fabs(l->getJacobianDeterminant(pts[j].pt[0],0,0)) * dsq;
- dmax = std::max(dmax,sqrt(dsq));
- }
- }
- d2 = sqrt(d2);
- return true;
-}
-
double GEdge::parFromPoint(const SPoint3 &P) const
{
double t;
@@ -527,6 +485,18 @@ void GEdge::relocateMeshVertices()
}
}
+SPoint3 GEdge :: closestPoint (SPoint3 &p, double tolerance)
+{
+ if (!_cp || _cp->tol() != tolerance) {
+ if(_cp)printf("coucou %12.15E %22.15E \n",tolerance,_cp->tol());
+ else printf("coucou %12.5E \n",tolerance);
+ if (_cp) delete _cp;
+ _cp = new closestPointFinder (this, tolerance);
+ }
+ return (*_cp)(p);
+}
+
+
typedef struct {
SPoint3 p;
double t;
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index e08fd177cd..dbb584c557 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -22,13 +22,14 @@ class MElement;
class MLine;
class ExtrudeParams;
class GEdgeCompound;
+class closestPointFinder;
// A model edge.
class GEdge : public GEntity {
private:
double _length;
bool _tooSmall;
-
+ closestPointFinder *_cp;
protected:
GVertex *v0, *v1;
// FIXME: normals need to be mutable at the moment, because thay can
@@ -208,8 +209,8 @@ class GEdge : public GEntity {
void addLine(MLine *line){ lines.push_back(line); }
- bool computeDistanceFromMeshToGeometry (double &d2, double &dmax);
virtual void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
+ SPoint3 closestPoint (SPoint3 &p, double tolerance);
};
#endif
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 73ae2c4e23..6ff21f8646 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -153,10 +153,10 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
wire.push_front(degenerated[0]);
}
else if (degenerated.size() > 2){
- Msg::Error("More than two degenerated edges in one model face");
+ Msg::Error("More than two degenerated edges in one model face of an OCC model");
}
- while (!wire.empty()){
+ while (!wire.empty() ){
// printf("wire.size = %d\n",wire.size());
loopTheLoop(wire,loop,°eneratedToInsert);
// break;
diff --git a/Geo/GModelIO_MSH2.cpp b/Geo/GModelIO_MSH2.cpp
index 3d3b555e13..3b41ed2131 100644
--- a/Geo/GModelIO_MSH2.cpp
+++ b/Geo/GModelIO_MSH2.cpp
@@ -844,6 +844,7 @@ int GModel::_writeMSH2(const std::string &name, double version, bool binary,
bool saveAll, bool saveParametric, double scalingFactor,
int elementStartNum, int saveSinglePartition, bool multipleView)
{
+
FILE *fp;
if(multipleView)
fp = Fopen(name.c_str(), binary ? "ab" : "a");
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index f1b7bbd0ee..3632cf7769 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -35,6 +35,37 @@ void OCC_Internals::buildLists()
addShapeToLists(shape);
}
+void OCC_Internals::buildShapeFromGModel(GModel* gm){
+ somap.Clear();
+ shmap.Clear();
+ fmap.Clear();
+ wmap.Clear();
+ emap.Clear();
+ vmap.Clear();
+ for (GModel::riter it = gm->firstRegion(); it != gm->lastRegion() ; ++it){
+ if ((*it)->getNativeType() == GEntity::OpenCascadeModel){
+ OCCRegion *occ = static_cast<OCCRegion*> (*it);
+ if (occ)addShapeToLists (occ->getTopoDS_Shape());
+ }
+ }
+ for (GModel::fiter it = gm->firstFace(); it != gm->lastFace() ; ++it){
+ if ((*it)->getNativeType() == GEntity::OpenCascadeModel){
+ OCCFace *occ = static_cast<OCCFace*> (*it);
+ if(occ)addShapeToLists (occ->getTopoDS_Face ());
+ }
+ }
+ BRep_Builder B;
+ TopoDS_Compound C;
+ B.MakeCompound(C);
+ for(int i = 1; i <= vmap.Extent(); i++) B.Add(C, vmap(i));
+ for(int i = 1; i <= emap.Extent(); i++) B.Add(C, emap(i));
+ for(int i = 1; i <= wmap.Extent(); i++) B.Add(C, wmap(i));
+ for(int i = 1; i <= fmap.Extent(); i++) B.Add(C, fmap(i));
+ for(int i = 1; i <= shmap.Extent(); i++) B.Add(C, shmap(i));
+ for(int i = 1; i <= somap.Extent(); i++) B.Add(C, somap(i));
+ shape = C;
+}
+
void OCC_Internals::buildShapeFromLists(TopoDS_Shape _shape)
{
BRep_Builder B;
@@ -211,6 +242,7 @@ void OCC_Internals::healGeometry(double tolerance, bool fixdegenerated,
bool fixsmalledges, bool fixspotstripfaces,
bool sewfaces, bool makesolids, bool connect)
{
+
if(!fixdegenerated && !fixsmalledges && !fixspotstripfaces &&
!sewfaces && !makesolids && !connect) return;
@@ -1047,6 +1079,8 @@ int GModel::readOCCIGES(const std::string &fn)
int GModel::writeOCCBREP(const std::string &fn)
{
+ _occ_internals->buildShapeFromGModel(this);
+
if(!_occ_internals){
Msg::Error("No OpenCASCADE model found");
return 0;
@@ -1058,6 +1092,7 @@ int GModel::writeOCCBREP(const std::string &fn)
int GModel::writeOCCSTEP(const std::string &fn)
{
+ _occ_internals->buildShapeFromGModel(this);
if(!_occ_internals){
Msg::Error("No OpenCASCADE model found");
return 0;
diff --git a/Geo/GModelIO_OCC.h b/Geo/GModelIO_OCC.h
index f45d293528..445a9d5943 100644
--- a/Geo/GModelIO_OCC.h
+++ b/Geo/GModelIO_OCC.h
@@ -27,6 +27,7 @@ class OCC_Internals {
TopoDS_Shape getShape () { return shape; }
void buildLists();
void buildShapeFromLists(TopoDS_Shape _shape);
+ void buildShapeFromGModel(GModel*);
void addShapeToLists(TopoDS_Shape shape);
void healGeometry(double tolerance, bool fixdegenerated,
bool fixsmalledges, bool fixspotstripfaces,
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index d24fa1a701..9f5b2f1cb6 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2384,30 +2384,33 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
for(int i = 0; i < List_Nbr(tmp); i++) {
Vertex *V;
List_Read(tmp, i, &V);
- pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol);
+ Vertex *found = v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)];
+ // printf("v(%d) = %d\n",V->Num,found->Num);
}
List_Delete(tmp);
// replace points in curves
tmp = Tree2List(GModel::current()->getGEOInternals()->Curves);
+
for(int i = 0; i < List_Nbr(tmp); i++) {
Curve *c;
List_Read(tmp, i, &c);
// replace begin/end points
c->beg = v2V[pos.find(c->beg->Pos.X, c->beg->Pos.Y, c->beg->Pos.Z, tol)];
c->end = v2V[pos.find(c->end->Pos.X, c->end->Pos.Y, c->end->Pos.Z, tol)];
+
// replace control points
for(int j = 0; j < List_Nbr(c->Control_Points); j++) {
Vertex *V;
List_Read(c->Control_Points, j, &V);
List_Write(c->Control_Points, j,
- &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]);
+ &v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]);
}
// replace extrusion sources
if(c->Extrude && c->Extrude->geo.Mode == EXTRUDED_ENTITY){
Vertex *V = FindPoint(std::abs(c->Extrude->geo.Source));
if(V) c->Extrude->geo.Source =
- v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]->Num;
+ v2V[pos.find(V->Pos.X, V->Pos.Y, V->Pos.Z, tol)]->Num;
}
}
List_Delete(tmp);
@@ -2791,6 +2794,158 @@ static void ReplaceDuplicateCurves(std::map<int, int> * c_report = 0)
Tree_Delete(allNonDuplicatedCurves);
}
+/*
+ 1) Find duplicate points and replace in curves
+ 2) Find duplicate curves and replace in surfaces
+ 3) Find duplicate surfaces and replace in volumes
+
+--> some curves are degenerate (zero length)
+--> some surfaces are degenerate (zero surface)
+--> some volumes are degenerate (zero volume)
+
+*/
+
+static void RemoveDegenerateCurves(){
+
+ { // remove degenerate curves from surface generatrices
+ List_T *All = Tree2List(GModel::current()->getGEOInternals()->Surfaces);
+ for(int i = 0; i < List_Nbr(All); i++) {
+ Surface *s;
+ List_Read(All, i, &s);
+ List_T *ll = s->Generatrices;
+ s->Generatrices = List_Create(4, 1, sizeof(Curve *));
+ // List_Delete(s->GeneratricesByTag);
+ // s->GeneratricesByTag = List_Create(4, 1, sizeof(int));
+ for(int j = 0; j < List_Nbr(ll); j++) {
+ Curve *c;
+ List_Read(ll, j, &c);
+ if (!c->degenerate()){
+ List_Add(s->Generatrices, &c);
+ // List_Add(s->GeneratricesByTag, &c->Num);
+ }
+ }
+ if (List_Nbr(ll) != List_Nbr(s->Generatrices)) Msg::Info("Coherence : Surface %d goes from %d to %d boundary curves",s->Num,List_Nbr(ll),List_Nbr(s->Generatrices));
+ List_Delete(ll);
+ }
+ }
+
+ { // actually remove the curves
+ List_T *All = Tree2List(GModel::current()->getGEOInternals()->Curves);
+ for(int i = 0; i < List_Nbr(All); i++) {
+ Curve *c;
+ List_Read(All, i, &c);
+ if(c->degenerate()) {
+ DeleteCurve(c->Num);
+ // DeleteCurve(-c->Num);
+ }
+ }
+ }
+}
+
+static void RemoveDegenerateVolumes(){
+ List_T *Vols = Tree2List(GModel::current()->getGEOInternals()->Volumes);
+ for(int k = 0; k < List_Nbr(Vols); k++) {
+ Volume *v;
+ List_Read(Vols, k, &v);
+ std::set<int> unique;
+ int N = List_Nbr(v->Surfaces);
+ for(int j = 0; j < N; j++) {
+ Surface *s;
+ List_Read(v->Surfaces, j, &s);
+ std::set<int>::iterator it = unique.find(-s->Num);
+ if (it == unique.end())unique.insert(s->Num);
+ else unique.erase(it);
+ }
+ if (N-unique.size()) Msg::Info("Coherence : Removing %d seams on Volume %d",N-unique.size(),v->Num);
+
+ List_T *ll= v->Surfaces;
+ List_T *ll2= v->SurfacesOrientations;
+ v->Surfaces = List_Create(1, 2, sizeof(Surface *));
+ v->SurfacesOrientations = List_Create(1, 2, sizeof(int));
+ for(int j = 0; j < List_Nbr(ll); j++) {
+ Surface *s;
+ List_Read(ll, j, &s);
+ if (unique.find(s->Num) != unique.end()){
+ List_Add(v->Surfaces,&s);
+ List_Add(v->SurfacesOrientations, List_Pointer(ll2, j));
+ }
+ }
+ List_Delete(ll);
+ List_Delete(ll2);
+ if (List_Nbr(v->Surfaces) == 0){
+ Msg::Info("Coherence Volume %d is removed (degenerated)",v->Num);
+ DeleteVolume(v->Num);
+ }
+ }
+}
+static void RemoveDegenerateSurfaces(){
+ List_T *All = Tree2List(GModel::current()->getGEOInternals()->Surfaces);
+
+ for(int i = 0; i < List_Nbr(All); i++) {
+ Surface *s;
+ std::set<int> unique;
+ List_Read(All, i, &s);
+ int N = List_Nbr(s->Generatrices);
+ for(int j = 0; j < N; j++) {
+ Curve *c;
+ List_Read(s->Generatrices, j, &c);
+ std::set<int>::iterator it = unique.find(-c->Num);
+ if (it == unique.end())unique.insert(c->Num);
+ else unique.erase(it);
+ }
+
+ if (N-unique.size()) Msg::Info("Coherence : Removing %d seams on Surface %d",N-unique.size(),s->Num);
+
+ List_T *ll = s->Generatrices;
+ s->Generatrices = List_Create(4, 1, sizeof(Curve *));
+ // List_Delete(s->GeneratricesByTag);
+ // s->GeneratricesByTag = List_Create(4, 1, sizeof(int));
+ for(int j = 0; j < List_Nbr(ll); j++) {
+ Curve *c;
+ List_Read(ll, j, &c);
+ if (unique.find(c->Num) != unique.end()){
+ List_Add(s->Generatrices,&c);
+ // List_Add(s->GeneratricesByTag, &c->Num);
+ }
+ }
+ List_Delete(ll);
+
+ if(s->degenerate()) {
+ Msg::Info("Coherence Surface %d is removed (degenerated)",s->Num);
+ List_T *Vols = Tree2List(GModel::current()->getGEOInternals()->Volumes);
+ for(int k = 0; k < List_Nbr(Vols); k++) {
+ Volume *v;
+ List_Read(Vols, k, &v);
+ List_T *ll= v->Surfaces;
+ List_T *ll2= v->SurfacesOrientations;
+ v->Surfaces = List_Create(1, 2, sizeof(Surface *));
+ v->SurfacesOrientations = List_Create(1, 2, sizeof(int));
+ for(int j = 0; j < List_Nbr(ll); j++) {
+ if(compareSurface(List_Pointer(ll, j), &s)){
+ List_Add(v->Surfaces, List_Pointer(ll, j));
+ List_Add(v->SurfacesOrientations, List_Pointer(ll2, j));
+ }
+ }
+ List_Delete (ll);
+ List_Delete (ll2);
+ }
+ DeleteSurface(s->Num);
+ }
+ }
+}
+
+bool Surface::degenerate() const {
+ int N = List_Nbr(Generatrices);
+ int Nd = 0;
+ for(int i = 0; i < N; i++) {
+ Curve *c;
+ List_Read(Generatrices, i, &c);
+ if(!c->degenerate())Nd++;
+ }
+ return Nd == 0;
+}
+
+
static void ReplaceDuplicateSurfaces(std::map<int, int> *s_report = 0)
{
Surface *s, *s2, **ps, **ps2;
@@ -2918,9 +3073,13 @@ static void ReplaceAllDuplicates(std::vector<std::map<int, int> > &report)
if(report.size() >= 3 && report[2].size())
surface_report = &(report[2]);
- ReplaceDuplicatePoints(vertex_report);
+ ReplaceDuplicatePoints();
ReplaceDuplicateCurves(curve_report);
ReplaceDuplicateSurfaces(surface_report);
+
+ RemoveDegenerateCurves();
+ RemoveDegenerateSurfaces();
+ RemoveDegenerateVolumes();
}
void ReplaceAllDuplicates()
diff --git a/Geo/Geo.h b/Geo/Geo.h
index f11a0edfd0..024c7adb16 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -155,6 +155,10 @@ class Curve{
Color.type = 1;
Color.geom = Color.mesh = value;
}
+ bool degenerate() const{
+ if (beg == end && Typ == MSH_SEGM_LINE)return true;
+ return false;
+ }
};
class EdgeLoop{
@@ -214,6 +218,7 @@ class Surface{
}
}
}
+ bool degenerate() const;
};
class SurfaceLoop{
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index f277c07b45..49cdd33cc7 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1567,3 +1567,19 @@ MElement *MElementFactory::create(int num, int type, const std::vector<int> &dat
return element;
}
+
+double MElement::skewness() {
+ double minsk = 1.0;
+ for (int i=0;i<getNumFaces();i++){
+ MFace f = getFace(i);
+ if (f.getNumVertices() == 3){
+ // MTriangle t (f.getVertex(0),f.getVertex(1),f.getVertex(2));
+ // minsk = std::min(minsk, t.etaShapeMeasure ());
+ }
+ else if (f.getNumVertices() == 4){
+ MQuadrangle q (f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3));
+ minsk = std::min(minsk, q.etaShapeMeasure ());
+ }
+ }
+ return minsk;
+}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7df20c6ae0..1e7bbe3972 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -179,6 +179,7 @@ class MElement
double maxDistToStraight();
// get the quality measures
+ double skewness();
virtual double rhoShapeMeasure();
virtual double gammaShapeMeasure(){ return 0.; }
virtual double etaShapeMeasure(){ return 0.; }
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 36dcfc1d10..41c2b69ef4 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -516,3 +516,36 @@ void MPrismN::getFaceVertices(const int num, std::vector<MVertex*> &v) const
_addFaceNodes(num,_order,_vs,ind,v);
}
+
+static double scaled_jacobian(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
+ double val;
+ double l1,l2,l3;
+ SVector3 vec1,vec2,vec3;
+
+ vec1 = SVector3(b->x()-a->x(),b->y()-a->y(),b->z()-a->z());
+ vec2 = SVector3(c->x()-a->x(),c->y()-a->y(),c->z()-a->z());
+ vec3 = SVector3(d->x()-a->x(),d->y()-a->y(),d->z()-a->z());
+
+ l1 = vec1.norm();
+ l2 = vec2.norm();
+ l3 = vec3.norm();
+
+ val = dot(vec1,crossprod(vec2,vec3));
+ return fabs(val)/(l1*l2*l3);
+}
+
+double MPrism::gammaShapeMeasure() {
+ MVertex *a = getVertex(0),*b= getVertex(1),*c= getVertex(2);
+ MVertex *d = getVertex(3),*e= getVertex(4),*f= getVertex(5);
+ const double j [6] = {
+ scaled_jacobian(a,b,c,d),
+ scaled_jacobian(b,a,c,e),
+ scaled_jacobian(c,a,b,f),
+ scaled_jacobian(d,a,e,f),
+ scaled_jacobian(e,b,d,f),
+ scaled_jacobian(f,c,d,e)};
+ const double result = *std::min_element(j,j+6);
+ printf("%12.5E\n",result);
+ return result;
+}
+
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 1ac48c5d0d..df37180fcf 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -129,6 +129,7 @@ class MPrism : public MElement {
tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
}
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
+ virtual double gammaShapeMeasure();
virtual int getVolumeSign();
virtual void getNode(int num, double &u, double &v, double &w) const
{
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 95c6313636..56aa05f5dc 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -7,7 +7,6 @@
#include "MQuadrangle.h"
#include "GaussLegendre1D.h"
#include "Context.h"
-#include "qualityMeasures.h"
#include "Numeric.h"
#include "BasisFactory.h"
diff --git a/Geo/MVertexPositionSet.h b/Geo/MVertexPositionSet.h
index 27f379e635..d8f3fc2c42 100644
--- a/Geo/MVertexPositionSet.h
+++ b/Geo/MVertexPositionSet.h
@@ -63,6 +63,7 @@ class MVertexPositionSet{
return _vertices[_index[i]];
}
if(_index[0] >= 0 && sqrt(_dist[0]) < tolerance){
+ // printf("tol %g dist %g x %g %g %g\n",tolerance,_dist[0],x,y,z);
_vertices[_index[0]]->setIndex(-1);
return _vertices[_index[0]];
}
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index a96851a753..72432af3f7 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -321,41 +321,39 @@ static void treat2Connections(GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2,
}
else if (fan){
- if (USEFANS__){
- int fanSize = FANSIZE__;
- // if the angle is greater than PI, than reverse the sense
- double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x());
- double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x());
- double AMAX = std::max(alpha1,alpha2);
- double AMIN = std::min(alpha1,alpha2);
- MEdge ee[2];
- if (alpha1 > alpha2){
- ee[0] = e2;ee[1] = e1;
- }
- else {
- ee[0] = e1;ee[1] = e2;
- }
- if ( AMAX - AMIN >= M_PI){
- double temp = AMAX;
- AMAX = AMIN + 2*M_PI;
- AMIN = temp;
- MEdge eee0 = ee[0];
- ee[0] = ee[1];ee[1] = eee0;
- }
- _columns->addFan (_myVert,ee[0],ee[1],true);
- for (int i=-1; i<=fanSize; i++){
- double t = (double)(i+1)/ (fanSize+1);
- double alpha = t * AMAX + (1.-t)* AMIN;
- SVector3 x (cos(alpha),sin(alpha),0);
- x.normalize();
- _dirs.push_back(x);
- }
+ int fanSize = FANSIZE__;
+ // if the angle is greater than PI, than reverse the sense
+ double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x());
+ double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x());
+ double AMAX = std::max(alpha1,alpha2);
+ double AMIN = std::min(alpha1,alpha2);
+ MEdge ee[2];
+ if (alpha1 > alpha2){
+ ee[0] = e2;ee[1] = e1;
}
else {
- _dirs.push_back(N1[SIDE]);
- _dirs.push_back(N2[SIDE]);
+ ee[0] = e1;ee[1] = e2;
+ }
+ if ( AMAX - AMIN >= M_PI){
+ double temp = AMAX;
+ AMAX = AMIN + 2*M_PI;
+ AMIN = temp;
+ MEdge eee0 = ee[0];
+ ee[0] = ee[1];ee[1] = eee0;
+ }
+ _columns->addFan (_myVert,ee[0],ee[1],true);
+ for (int i=-1; i<=fanSize; i++){
+ double t = (double)(i+1)/ (fanSize+1);
+ double alpha = t * AMAX + (1.-t)* AMIN;
+ SVector3 x (cos(alpha),sin(alpha),0);
+ x.normalize();
+ _dirs.push_back(x);
}
}
+ else {
+ _dirs.push_back(N1[SIDE]);
+ _dirs.push_back(N2[SIDE]);
+ }
}
}
}
@@ -969,6 +967,41 @@ fanTopology :: fanTopology (GRegion * gr, const std::set<GEdge*> &detectedFans,
}
}
+// This is the tricky part
+// We have to find a vector N that has the following properties
+// V = (x,y,z) / (x^2+y^2+z^2)^{1/2} is a point on the unit sphere
+// the n[i]'s are points on the unit sphere
+// V maximizes min_i (V * n[i])
+// This means I'd like to find point V that is
+
+static void filterVectors(std::vector<SVector3> &n){
+ std::vector<SVector3> temp;
+ temp.push_back(n[0]);
+ for (unsigned int i = 1 ; i<n.size() ; i++){
+ bool found = false;
+ for (unsigned int j = 0 ; j<temp.size() ; j++){
+ double d = dot(n[i],temp[j]);
+ if (d < 0.98)found = true;
+ }
+ if (found) temp.push_back(n[i]);
+ }
+ n = temp;
+}
+
+static SVector3 computeBestNormal(std::vector<SVector3> &n){
+ // filterVectors (n);
+ SVector3 V;
+ if (n.size() == 1)V = n[0];
+ else if (n.size() == 2)V = n[0]+n[1];
+ else if (n.size() == 3)circumCenterXYZ(n[0],n[1],n[2],V);
+ else {
+ Msg::Warning("suboptimal choice for exterior normal: %d vectors to average",n.size());
+ for (unsigned int i = 0 ; i<n.size() ; i++)V+=n[i];
+ }
+ V.normalize();
+ return V;
+}
+
static int createColumnsBetweenFaces(GRegion *gr,
MVertex *myV,
@@ -980,7 +1013,6 @@ static int createColumnsBetweenFaces(GRegion *gr,
fanTopology &ft)
{
SVector3 n[256];
- SPoint3 c[256];
int count = 0;
GFace *gfs[256];
@@ -991,8 +1023,8 @@ static int createColumnsBetweenFaces(GRegion *gr,
_faces.lower_bound(*it);
itm != _faces.upper_bound(*it); ++itm){
- n[count] += _normals[itm->second->getFace(0)];
- c[count] = itm->second->getFace(0).barycenter();
+ SVector3 N = _normals[itm->second->getFace(0)];
+ n[count] += N;
}
gfs[count] = *it;
n[count].normalize();
@@ -1022,10 +1054,16 @@ static int createColumnsBetweenFaces(GRegion *gr,
for (std::multimap<int, GFace*>::iterator itm = range.first ; itm != range.second ; itm++)
joint.push_back(itm->second);
joints.push_back(joint);
- SVector3 avg (0,0,0);
+ // SVector3 avg (0,0,0);
+ std::vector<SVector3> ns;
+
for (unsigned int i=0;i<joint.size(); i++){
- avg += n[inv[joint[i]]];
+ ns.push_back( n[inv[joint[i]]] );
+ // avg += n[inv[joint[i]]];
}
+ SVector3 avg = computeBestNormal(ns);
+
+
std::vector<MVertex*> _column;
std::vector<SMetric3> _metrics;
avg.normalize();
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index d23172b307..7780d8a1ff 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -20,6 +20,15 @@ gmshEdge::gmshEdge(GModel *m, Curve *edge, GVertex *v1, GVertex *v2)
resetMeshAttributes();
}
+bool gmshEdge::degenerate(int dim) const {
+ if (c->beg == c->end && c->Typ == MSH_SEGM_LINE){
+ Msg::Info("Model Edge %d is degenerate",tag());
+ return true;
+ }
+ return false;
+}
+
+
void gmshEdge::resetMeshAttributes()
{
meshAttributes.method = c->Method;
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index 3c92afd14a..c66d33cac7 100644
--- a/Geo/gmshEdge.h
+++ b/Geo/gmshEdge.h
@@ -30,6 +30,7 @@ class gmshEdge : public GEdge {
virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
virtual void writeGEO(FILE *fp);
void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
+ virtual bool degenerate(int dim) const;
};
#endif
diff --git a/Geo/gmshEdgeDiscretize.cpp b/Geo/gmshEdgeDiscretize.cpp
new file mode 100644
index 0000000000..b85e429a10
--- /dev/null
+++ b/Geo/gmshEdgeDiscretize.cpp
@@ -0,0 +1,196 @@
+#include <cstdio>
+#include <cmath>
+#include <vector>
+#include "SPoint3.h"
+#include "SVector3.h"
+#include "GEdge.h"
+#include "gmshEdge.h"
+#include "Geo.h"
+
+
+class discreteList {
+ std::vector<std::pair<SPoint3, double> > _pts;
+ std::vector<int> _next;
+ public:
+ int insertPoint(int pos, const SPoint3 &pt, double t) {
+ _pts.push_back(std::make_pair(pt, t));
+ _next.push_back(_next[pos + 1]);
+ _next[pos + 1] = _pts.size() - 1;
+ return _pts.size() - 1;
+ }
+ void sort(std::vector<SPoint3> &spts, std::vector<double> &ts) {
+ spts.clear();
+ spts.reserve(_pts.size());
+ ts.clear();
+ ts.reserve(_pts.size());
+ for (int p = _next[0]; p != -1; p = _next[p + 1]) {
+ spts.push_back(_pts[p].first);
+ ts.push_back(_pts[p].second);
+ }
+ }
+ discreteList() {
+ _next.push_back(-1);
+ }
+};
+
+static void decasteljau(double tol, discreteList &discrete, int pos, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p3, double t0, double t3)
+{
+ SVector3 d30 = p3 - p0;
+ SVector3 d13 = p1 - p3;
+ SVector3 d23 = p2 - p3;
+ SVector3 d130 = crossprod(d13, d30);
+ SVector3 d230 = crossprod(d23, d30);
+ double d = std::max(dot(d130, d130), dot(d230, d230));
+ double l2 = dot(d30, d30);
+
+ if(d < tol * tol * l2) {
+ return;
+ }
+
+ SPoint3 p01((p0 + p1) * 0.5);
+ SPoint3 p12((p1 + p2) * 0.5);
+ SPoint3 p23((p2 + p3) * 0.5);
+ SPoint3 p012((p01 + p12) * 0.5);
+ SPoint3 p123((p12 + p23) * 0.5);
+ SPoint3 p0123((p012 + p123) * 0.5);
+ double t0123 = 0.5 * (t0 + t3);
+ int newpos = discrete.insertPoint(pos, p0123, t0123);
+
+ decasteljau(tol, discrete, pos, p0, p01, p012, p0123, t0, t0123);
+ decasteljau(tol, discrete, newpos, p0123, p123, p23, p3, t0123, t3);
+}
+
+static int discretizeBezier(double tol, discreteList &discrete, int pos, const SPoint3 pt[4], double t0, double t3, bool insertFirstPoint)
+{
+ if (insertFirstPoint)
+ pos = discrete.insertPoint(pos, pt[0], t0);
+ int newp = discrete.insertPoint(pos, pt[3], t3);
+ decasteljau(tol, discrete, pos, pt[0], pt[1], pt[2], pt[3], t0, t3);
+ return newp;
+}
+
+static int discretizeBSpline(double tol, discreteList &discrete, int pos, const SPoint3 pt[4], double t0, double t3, bool insertFirstPoint)
+{
+ SPoint3 bpt[4] = {
+ SPoint3((pt[0] + 4 * pt[1] + pt[2]) * (1./6.)),
+ SPoint3((2 * pt[1] + pt[2]) * (1./3.)),
+ SPoint3((pt[1] + 2 * pt[2]) * (1./3.)),
+ SPoint3((pt[1] + 4 * pt[2] + pt[3]) * (1./6.))
+ };
+ return discretizeBezier(tol, discrete, pos, bpt, t0, t3, insertFirstPoint);
+}
+
+static int discretizeCatmullRom(double tol, discreteList &discrete, int pos, const SPoint3 pt[4], double t0, double t3, bool insertFirstPoint)
+{
+ SPoint3 bpt[4] = {
+ pt[1],
+ SPoint3(( 6 * pt[1] + pt[2] - pt[0]) * (1./6.)),
+ SPoint3(( 6 * pt[2] - pt[3] + pt[1]) * (1./6.)),
+ pt[2]
+ };
+ return discretizeBezier(tol, discrete, pos, bpt, t0, t3, insertFirstPoint);
+}
+
+static inline SPoint3 curveGetPoint(Curve *c, int i)
+{
+ Vertex *v;
+ List_Read(c->Control_Points, i , &v);
+ return SPoint3(v->Pos.X, v->Pos.Y, v->Pos.Z);
+}
+
+static void discretizeCurve(Curve *c, double tol, std::vector<SPoint3> &pts, std::vector<double> &ts)
+{
+ discreteList discrete;
+ switch(c->Typ) {
+ case MSH_SEGM_LINE :
+ {
+ int NPt = List_Nbr(c->Control_Points);
+ pts.resize(NPt);
+ ts.resize(NPt);
+ for (int i = 0; i < NPt; ++i) {
+ pts[i]= curveGetPoint(c, i);
+ ts[i] = i / (double) (NPt - 1);
+ }
+ return;
+ }
+ case MSH_SEGM_BEZIER :
+ {
+ int back = -1;
+ int NbCurves = (List_Nbr(c->Control_Points) - 1) / 3;
+ for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
+ double t1 = (iCurve) / (double)(NbCurves);
+ double t2 = (iCurve+1) / (double)(NbCurves);
+ SPoint3 pt[4];
+ for(int i = 0; i < 4; i++) {
+ pt[i] = curveGetPoint(c, iCurve * 3 + i);
+ }
+ back = discretizeBezier(tol, discrete, back, pt, t1, t2, iCurve == 0);
+ }
+ break;
+ }
+ case MSH_SEGM_BSPLN:
+ {
+ int back = -1;
+ bool periodic = (c->end == c->beg);
+ int NbControlPoints = List_Nbr(c->Control_Points);
+ int NbCurves = NbControlPoints + (periodic ? -1 : 1);
+ SPoint3 pt[4];
+ for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
+ double t1 = (iCurve) / (double)(NbCurves);
+ double t2 = (iCurve+1) / (double)(NbCurves);
+ for(int i = 0; i < 4; i++) {
+ int k;
+ if (periodic) {
+ k = (iCurve - 1 + i) % (NbControlPoints - 1);
+ if (k < 0)
+ k += NbControlPoints - 1;
+ }
+ else {
+ k = std::max(0, std::min(iCurve - 2 + i, NbControlPoints -1));
+ }
+ pt[i] = curveGetPoint(c, k);
+ }
+ back = discretizeBSpline(tol, discrete, back, pt, t1, t2, iCurve == 0);
+ }
+ break;
+ }
+ case MSH_SEGM_SPLN:
+ {
+ int NbCurves = List_Nbr(c->Control_Points) - 1;
+ SPoint3 pt[4];
+ int back = -1;
+ for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
+ double t1 = (iCurve) / (double)(NbCurves);
+ double t2 = (iCurve+1) / (double)(NbCurves);
+ pt[1] = curveGetPoint(c, iCurve);
+ pt[2] = curveGetPoint(c, iCurve + 1);
+ if(iCurve == 0) {
+ if(c->beg == c->end)
+ pt[0] = curveGetPoint(c, NbCurves - 1);
+ else
+ pt[0] = SPoint3(pt[1] * 2 - pt[2]);
+ }
+ else
+ pt[0] = curveGetPoint(c, iCurve - 1);
+ if(iCurve == NbCurves - 1) {
+ if(c->beg == c->end)
+ pt[3] = curveGetPoint(c, 1);
+ else
+ pt[3] = SPoint3(2 * pt[2] - pt[1]);
+ }
+ else
+ pt[3] = curveGetPoint(c, iCurve + 2);
+ back = discretizeCatmullRom(tol, discrete, back, pt, t1, t2, iCurve == 0);
+ }
+ break;
+ }
+ default :
+ Msg::Fatal("not implemented");
+ }
+ discrete.sort(pts, ts);
+}
+
+void gmshEdge::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts)
+{
+ discretizeCurve(c, tol, dpts, ts);
+}
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index f6f5f3c19f..dccd862052 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1916,7 +1916,7 @@ void BoundaryLayerField::setupFor2d(int iF)
nodes_id.push_back ((*it)->getEndVertex()->tag());
}
}
- printf("face %d %d BL Edges\n", iF, (int)edges_id.size());
+ // printf("face %d %d BL Edges\n", iF, (int)edges_id.size());
removeAttractors();
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index e7d8f378b8..2a5a9887ef 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1519,20 +1519,6 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
Msg::StatusBar(true, "Done meshing order %d (%g s)", order, t2 - t1);
}
-void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t &dist)
-{
- for (GModel::eiter itEdge = m->firstEdge(); itEdge != m->lastEdge(); ++itEdge) {
- double d2,dmax;
- (*itEdge)->computeDistanceFromMeshToGeometry(d2, dmax);
- dist.d2[*itEdge] = d2;
- dist.d_max[*itEdge] = dmax;
- }
-
- for (GModel::fiter itFace = m->firstFace(); itFace != m->lastFace(); ++itFace) {
-
- }
-}
-
void SetHighOrderComplete(GModel *m, bool onlyVisible)
{
faceContainer faceVertices;
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index 7164d06985..5d2ccea882 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -31,11 +31,6 @@ void checkHighOrderTriangles(const char* cc, GModel *m,
void checkHighOrderTetrahedron(const char* cc, GModel *m,
std::vector<MElement*> &bad, double &minJGlob);
-struct distanceFromMeshToGeometry_t {
- std::map<GEntity*, double> d_max, d2;
-};
-
-void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t &dist);
void getMeshInfoForHighOrder(GModel *gm, int &meshOrder, bool &complete, bool &CAD);
#endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 35e2746c74..beae02b54c 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -708,7 +708,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
fprintf(ff2,"};\n");
fclose(ff2);
- filterOverlappingElements (blTris,blQuads,_columns->_elemColumns,_columns->_toFirst);
+ // filterOverlappingElements (blTris,blQuads,_columns->_elemColumns,_columns->_toFirst);
for (unsigned int i = 0; i < blQuads.size();i++){
addOrRemove(blQuads[i]->getVertex(0),blQuads[i]->getVertex(1),bedges);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 3c4e7f852e..74fe7e0f27 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -45,6 +45,8 @@ class splitQuadRecovery {
std::multimap<GEntity*, std::pair<MVertex*,MFace> >_data;
bool _empty;
public :
+ std::map<MFace, MVertex*, Less_Face>_invmap;
+ std::set<MFace, Less_Face>_toDelete;
splitQuadRecovery() : _empty(true) {}
bool empty(){ return _empty; }
void setEmpty(bool val){ _empty = val; }
@@ -65,43 +67,49 @@ class splitQuadRecovery {
(*it)->triangles.clear();
for (std::multimap<GEntity*, std::pair<MVertex*,MFace> >::iterator it2 =
_data.lower_bound(*it); it2 != _data.upper_bound(*it) ; ++it2){
- MVertex *v = it2->second.first;
- v->onWhat()->mesh_vertices.erase(std::find(v->onWhat()->mesh_vertices.begin(),
- v->onWhat()->mesh_vertices.end(), v));
const MFace &f = it2->second.second;
- std::set<MFace, Less_Face>::iterator itf0 = allFaces.find(MFace(f.getVertex(0),
- f.getVertex(1),v));
- std::set<MFace, Less_Face>::iterator itf1 = allFaces.find(MFace(f.getVertex(1),
- f.getVertex(2),v));
- std::set<MFace, Less_Face>::iterator itf2 = allFaces.find(MFace(f.getVertex(2),
- f.getVertex(3),v));
- std::set<MFace, Less_Face>::iterator itf3 = allFaces.find(MFace(f.getVertex(3),
- f.getVertex(0),v));
- if (itf0 != allFaces.end() && itf1 != allFaces.end() &&
- itf2 != allFaces.end() && itf3 != allFaces.end()){
- (*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0), f.getVertex(1),
- f.getVertex(2), f.getVertex(3)));
- allFaces.erase(*itf0);
- allFaces.erase(*itf1);
- allFaces.erase(*itf2);
- allFaces.erase(*itf3);
- // printf("some pyramids should be created %d regions\n", (*it)->numRegions());
- for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){
- if (iReg == 1) {
- Msg::Error("Cannot build pyramids on non manifold faces");
- v = new MVertex(v->x(), v->y(), v->z(), (*it)->getRegion(iReg));
- }
- else
- v->setEntity ((*it)->getRegion(iReg));
- (*it)->getRegion(iReg)->pyramids.push_back
- (new MPyramid(f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), v));
- (*it)->getRegion(iReg)->mesh_vertices.push_back(v);
- NBPY++;
+ MVertex *v = it2->second.first;
+ v->onWhat()->mesh_vertices.erase(std::find(v->onWhat()->mesh_vertices.begin(),
+ v->onWhat()->mesh_vertices.end(), v));
+ std::set<MFace, Less_Face>::iterator itf0 = allFaces.find(MFace(f.getVertex(0),
+ f.getVertex(1),v));
+ std::set<MFace, Less_Face>::iterator itf1 = allFaces.find(MFace(f.getVertex(1),
+ f.getVertex(2),v));
+ std::set<MFace, Less_Face>::iterator itf2 = allFaces.find(MFace(f.getVertex(2),
+ f.getVertex(3),v));
+ std::set<MFace, Less_Face>::iterator itf3 = allFaces.find(MFace(f.getVertex(3),
+ f.getVertex(0),v));
+ if (itf0 != allFaces.end() && itf1 != allFaces.end() &&
+ itf2 != allFaces.end() && itf3 != allFaces.end()){
+ (*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0), f.getVertex(1),
+ f.getVertex(2), f.getVertex(3)));
+ allFaces.erase(*itf0);
+ allFaces.erase(*itf1);
+ allFaces.erase(*itf2);
+ allFaces.erase(*itf3);
+ // printf("some pyramids should be created %d regions\n", (*it)->numRegions());
+ for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){
+ if (iReg == 1) {
+ Msg::Error("Cannot build pyramids on non manifold faces");
+ v = new MVertex(v->x(), v->y(), v->z(), (*it)->getRegion(iReg));
+ }
+ else
+ v->setEntity ((*it)->getRegion(iReg));
+ // A quad face connected to an hex or a primsm --> leave the quad face as is
+ if (_toDelete.find(f) == _toDelete.end()){
+ (*it)->getRegion(iReg)->pyramids.push_back
+ (new MPyramid(f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), v));
+ (*it)->getRegion(iReg)->mesh_vertices.push_back(v);
+ NBPY++;
+ }
+ else {
+ delete v;
+ }
+ }
}
- }
}
for (std::set<MFace, Less_Face>::iterator itf = allFaces.begin();
- itf != allFaces.end(); ++itf){
+ itf != allFaces.end(); ++itf){
(*it)->triangles.push_back
(new MTriangle(itf->getVertex(0), itf->getVertex(1), itf->getVertex(2)));
}
@@ -258,6 +266,7 @@ void getBoundingInfoAndSplitQuads(GRegion *gr,
(v0->y() + v1->y() + v2->y() + v3->y())*0.25,
(v0->z() + v1->z() + v2->z() + v3->z())*0.25,itx->second);
sqr.add(f,newv,itx->second);
+ sqr._invmap[f] = newv;
allBoundingFaces[MFace(v0,v1,newv)] = itx->second;
allBoundingFaces[MFace(v1,v2,newv)] = itx->second;
allBoundingFaces[MFace(v2,v3,newv)] = itx->second;
@@ -289,6 +298,15 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
std::map<MFace,GEntity*,Less_Face> allBoundingFaces;
getBoundingInfoAndSplitQuads(gr, allBoundingFaces, allBoundingVertices, sqr);
+ //// TEST
+ {
+ std::vector<MVertex*>ALL;
+ std::vector<MTetrahedron*> MESH;
+ ALL.insert(ALL.begin(),allBoundingVertices.begin(),allBoundingVertices.end());
+ // delaunayMeshIn3D (ALL,MESH);
+ // exit(1);
+ }
+
in.mesh_dim = 3;
in.firstnumber = 1;
in.numberofpoints = allBoundingVertices.size() + Filler::get_nbr_new_vertices() +
@@ -504,8 +522,20 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
static void addOrRemove(const MFace &f,
MElement *e,
- std::map<MFace,MElement*,Less_Face> & bfaces)
+ std::map<MFace,MElement*,Less_Face> & bfaces,
+ splitQuadRecovery &sqr)
{
+ {
+ std::map<MFace, MVertex*, Less_Face>::const_iterator it = sqr._invmap.find(f);
+ if (it != sqr._invmap.end()){
+ addOrRemove (MFace(it->second, f.getVertex(0),f.getVertex(1)),e,bfaces,sqr);
+ addOrRemove (MFace(it->second, f.getVertex(1),f.getVertex(2)),e,bfaces,sqr);
+ addOrRemove (MFace(it->second, f.getVertex(2),f.getVertex(3)),e,bfaces,sqr);
+ addOrRemove (MFace(it->second, f.getVertex(3),f.getVertex(0)),e,bfaces,sqr);
+ return;
+ }
+ }
+
std::map<MFace,MElement*,Less_Face>::iterator it = bfaces.find(f);
if (it == bfaces.end())bfaces.insert(std::make_pair(f,e));
else bfaces.erase(it);
@@ -523,11 +553,12 @@ static void addOrRemove(const MFace &f,
*/
-bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
- std::vector<MTetrahedron*> &tets,
- std::vector<MHexahedron*> &hexes,
- std::vector<MPrism*> &prisms,
- std::vector<MPyramid*> &pyramids)
+static bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
+ std::vector<MTetrahedron*> &tets,
+ std::vector<MHexahedron*> &hexes,
+ std::vector<MPrism*> &prisms,
+ std::vector<MPyramid*> &pyramids,
+ splitQuadRecovery & sqr)
{
std::set<MElement*> all;
all.insert(hexes.begin(),hexes.end());
@@ -555,7 +586,7 @@ bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
else {
// what to do ??????
// two tets and one prism --> the prism should be
- // geometrically on the other side of the
+ // geometrically on the other side of the tet
if (itf->second.second) {
MElement *prism=0, *t1=0, *t2=0;
if (itf->second.second->getType () == TYPE_PRI || itf->second.second->getType () == TYPE_PYR) {
@@ -621,15 +652,21 @@ bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
connected.insert(FIRST);
for (int i=0;i<FIRST->getNumFaces();i++){
MFace f = FIRST->getFace(i);
- GFace* gfound = findInFaceSearchStructure (f.getVertex(0),
- f.getVertex(1),
- f.getVertex(2),
- search);
+ std::map<MFace, MVertex*, Less_Face>::iterator it = sqr._invmap.find(f);
+ GFace* gfound = 0;
+ if (it != sqr._invmap.end()){
+ gfound = (GFace*)it->second->onWhat();
+ // one pyramid is useless because one element with a quad face impacts the
+ // boundary of the domain.
+ sqr._toDelete.insert(f);
+ }
+ else gfound = findInFaceSearchStructure (f,search);
if (!gfound){
std::map<MFace,std::pair<MElement*,MElement*>,Less_Face>::iterator
itf = myGraph.find(f);
MElement *t_neigh = itf->second.first == FIRST ?
itf->second.second : itf->second.first;
+ if (!t_neigh)printf("oulalalalalalalala %d vertices\n",f.getNumVertices());
if (connected.find(t_neigh) == connected.end())myStack.push(t_neigh);
}
else {
@@ -638,6 +675,7 @@ bool AssociateElementsToModelRegionWithBoundaryLayers (GRegion *gr,
}
}
}
+
// printf ("found a set of %d elements that are connected with %d bounding faces\n",connected.size(),faces_bound.size());
GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
// printf("REGION %d %d\n",myGRegion->tag(),gr->tag());
@@ -770,7 +808,8 @@ static int getWedge(BoundaryLayerColumns* _columns, MVertex *v1, MVertex *v2,
return fanSize + 2;
}
-static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
+
+static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr, splitQuadRecovery & sqr)
{
if (getBLField(gr->model())) insertVerticesInRegion(gr,-1);
if (!buildAdditionalPoints3D (gr)) return false;
@@ -786,20 +825,34 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
faces.insert(faces.begin(), embedded_faces.begin(),embedded_faces.end());
std::set<MVertex*> verts;
+ {
+ std::list<GFace*>::iterator itf = faces.begin();
+ while(itf != faces.end()){
+ for(unsigned int i = 0; i< (*itf)->getNumMeshElements(); i++){
+ MElement *e = (*itf)->getMeshElement(i);
+ addOrRemove (e->getFace(0),0,bfaces,sqr);
+ }
+ ++itf;
+ }
+ }
+
std::list<GFace*>::iterator itf = faces.begin();
while(itf != faces.end()){
for(unsigned int i = 0; i< (*itf)->triangles.size(); i++){
MVertex *v1 = (*itf)->triangles[i]->getVertex(0);
MVertex *v2 = (*itf)->triangles[i]->getVertex(1);
MVertex *v3 = (*itf)->triangles[i]->getVertex(2);
- MFace dv(v1,v2,v3);
- addOrRemove (dv,0,bfaces);
+ MFace dv (v1,v2,v3);
for (unsigned int SIDE = 0 ; SIDE < _columns->_normals3D.count(dv); SIDE ++){
faceColumn fc = _columns->getColumns(*itf,v1, v2, v3, SIDE);
const BoundaryLayerData & c1 = fc._c1;
const BoundaryLayerData & c2 = fc._c2;
const BoundaryLayerData & c3 = fc._c3;
int N = std::min(c1._column.size(),std::min(c2._column.size(),c3._column.size()));
+
+ // double distMax = getDistMax (v1, v2, v3, c1._n, c2._n, c3._n);
+ MFace f_low (v1,v2,v3);
+ SVector3 n_low = f_low.normal();
// printf("%d Layers\n",N);
std::vector<MElement*> myCol;
for (int l=0;l < N ;++l){
@@ -817,12 +870,19 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
v12 = c2._column[l-1];
v13 = c3._column[l-1];
}
- // printf("coucoucouc %p %p %p %p %p %p\n",v11,v12,v13,v21,v22,v23);
+ MFace f_up (v21,v22,v23);
+ SVector3 n_up = f_up.normal();
+ double dotProd = dot(n_up,n_low);
MPrism *prism = new MPrism(v11,v12,v13,v21,v22,v23);
- // store the layer the element belongs
- prism->setPartition(l+1);
- blPrisms.push_back(prism);
- myCol.push_back(prism);
+ if (dotProd > 0.2 && prism->skewness() > 0.1){
+ prism->setPartition(l+1);
+ blPrisms.push_back(prism);
+ myCol.push_back(prism);
+ }
+ else {
+ delete prism;
+ l = N+1;
+ }
}
if (!myCol.empty()){
for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
@@ -909,8 +969,10 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
}
++ite;
}
-
- // filterOverlappingElements (blPrisms,blHexes,_columns->_elemColumns,_columns->_toFirst);
+ // ------------------------------------------------------------------------------------
+ // FIXME : NOT 100 % CORRECT
+ // filterOverlappingElements (blPrisms,blHexes,_columns->_elemColumns,_columns->_toFirst);
+ // ------------------------------------------------------------------------------------
{
FILE *ff2 = fopen ("tato3D.pos","w");
fprintf(ff2,"View \" \"{\n");
@@ -926,13 +988,13 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
for (unsigned int i = 0; i < blPrisms.size();i++){
for (unsigned int j=0;j<5;j++)
- addOrRemove(blPrisms[i]->getFace(j),blPrisms[i],bfaces);
+ addOrRemove(blPrisms[i]->getFace(j),blPrisms[i],bfaces,sqr);
for (int j = 0; j < 6; j++)
if(blPrisms[i]->getVertex(j)->onWhat() == gr)verts.insert(blPrisms[i]->getVertex(j));
}
for (unsigned int i = 0; i < blHexes.size();i++){
for (unsigned int j=0;j<6;j++)
- addOrRemove(blHexes[i]->getFace(j),blHexes[i],bfaces);
+ addOrRemove(blHexes[i]->getFace(j),blHexes[i],bfaces, sqr);
for (int j = 0; j < 8; j++)
if(blHexes[i]->getVertex(j)->onWhat() == gr)verts.insert(blHexes[i]->getVertex(j));
}
@@ -1037,7 +1099,11 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
printf("%d tets\n", (int)gr->tetrahedra.size());
deMeshGFace _kill;
_kill (nf);
+ //<<<<<<< .mine
gr->model()->remove(nf);
+ // delete nf;
+ //=======
+ //>>>>>>> .r18259
delete nf;
gr->set(faces);
@@ -1045,8 +1111,7 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
gr->model()->writeMSH("BL_start.msh");
- AssociateElementsToModelRegionWithBoundaryLayers(gr, gr->tetrahedra, blHexes,
- blPrisms, blPyrs);
+ AssociateElementsToModelRegionWithBoundaryLayers (gr, gr->tetrahedra , blHexes, blPrisms, blPyrs, sqr);
gr->model()->writeMSH("BL_start2.msh");
@@ -1224,7 +1289,7 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
- bool _BL = modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr);
+ bool _BL = modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr,sqr);
// now do insertion of points
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
@@ -1274,7 +1339,7 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
}
// uncomment this to test the new code
-#define NEW_CODE
+//##define NEW_CODE
void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
{
@@ -1284,6 +1349,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
MeshDelaunayVolumeTetgen(regions);
return;
#endif
+ splitQuadRecovery sqr;
for(unsigned int i = 0; i < regions.size(); i++)
Msg::Info("Meshing volume %d (Delaunay)", regions[i]->tag());
@@ -1326,7 +1392,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
- bool _BL = modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr);
+ bool _BL = modifyInitialMeshForTakingIntoAccountBoundaryLayers(gr,sqr);
// now do insertion of points
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
@@ -1471,8 +1537,12 @@ void deMeshGRegion::operator() (GRegion *gr)
gr->deleteMesh();
}
+/// X_1 (1-u-v) + X_2 u + X_3 v = P_x + t N_x
+/// Y_1 (1-u-v) + Y_2 u + Y_3 v = P_y + t N_y
+/// Z_1 (1-u-v) + Z_2 u + Z_3 v = P_z + t N_z
+
int intersect_line_triangle(double X[3], double Y[3], double Z[3] ,
- double P[3], double N[3], double eps_prec)
+ double P[3], double N[3], const double eps_prec)
{
double mat[3][3], det;
double b[3], res[3];
@@ -1494,8 +1564,10 @@ int intersect_line_triangle(double X[3], double Y[3], double Z[3] ,
b[2] = P[2] - Z[0];
if(!sys3x3_with_tol(mat, b, res, &det))
- return 0;
-
+ {
+ return 0;
+ }
+ // printf("coucou %g %g %g\n",res[0],res[1],res[2]);
if(res[0] >= eps_prec && res[0] <= 1.0 - eps_prec &&
res[1] >= eps_prec && res[1] <= 1.0 - eps_prec &&
1 - res[0] - res[1] >= eps_prec && 1 - res[0] - res[1] <= 1.0 - eps_prec){
@@ -1509,6 +1581,7 @@ int intersect_line_triangle(double X[3], double Y[3], double Z[3] ,
return 0;
}
else{
+ printf("non robust stuff\n");
// the intersection is not robust, try another triangle
return -10000;
}
@@ -1709,6 +1782,7 @@ void optimizeMeshGRegionGmsh::operator() (GRegion *gr)
optimizeMesh(gr, QMTET_2);
}
+
bool buildFaceSearchStructure(GModel *model, fs_cont &search)
{
search.clear();
@@ -1723,13 +1797,9 @@ bool buildFaceSearchStructure(GModel *model, fs_cont &search)
std::set<GFace*>::iterator fit = faces_to_consider.begin();
while(fit != faces_to_consider.end()){
- for(unsigned int i = 0; i < (*fit)->triangles.size(); i++){
- MVertex *p1 = (*fit)->triangles[i]->getVertex(0);
- MVertex *p2 = (*fit)->triangles[i]->getVertex(1);
- MVertex *p3 = (*fit)->triangles[i]->getVertex(2);
- MVertex *p = std::min(p1, std::min(p2, p3));
- search.insert(std::pair<MVertex*, std::pair<MTriangle*, GFace*> >
- (p, std::pair<MTriangle*, GFace*>((*fit)->triangles[i], *fit)));
+ for(unsigned int i = 0; i < (*fit)->getNumMeshElements(); i++){
+ MFace ff = (*fit)->getMeshElement(i)->getFace(0);
+ search[ff] = *fit;
}
++fit;
}
@@ -1757,21 +1827,21 @@ bool buildEdgeSearchStructure(GModel *model, es_cont &search)
GFace *findInFaceSearchStructure(MVertex *p1, MVertex *p2, MVertex *p3,
const fs_cont &search)
{
- MVertex *p = std::min(p1, std::min(p2, p3));
+ MFace ff(p1,p2,p3);
+ fs_cont::const_iterator it = search.find(ff);
+ if (it == search.end())return 0;
+ return it->second;
+}
- for(fs_cont::const_iterator it = search.lower_bound(p);
- it != search.upper_bound(p);
- ++it){
- MTriangle *t = it->second.first;
- GFace *gf= it->second.second;
- if((t->getVertex(0) == p1 || t->getVertex(0) == p2 || t->getVertex(0) == p3) &&
- (t->getVertex(1) == p1 || t->getVertex(1) == p2 || t->getVertex(1) == p3) &&
- (t->getVertex(2) == p1 || t->getVertex(2) == p2 || t->getVertex(2) == p3))
- return gf;
- }
- return 0;
+GFace *findInFaceSearchStructure(const MFace &ff,
+ const fs_cont &search)
+{
+ fs_cont::const_iterator it = search.find(ff);
+ if (it == search.end())return 0;
+ return it->second;
}
+
GEdge *findInEdgeSearchStructure(MVertex *p1, MVertex *p2, const es_cont &search)
{
MVertex *p = std::min(p1, p2);
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index 2469a17219..5e7ab28944 100644
--- a/Mesh/meshGRegion.h
+++ b/Mesh/meshGRegion.h
@@ -9,6 +9,7 @@
#include <list>
#include <vector>
#include <map>
+#include "MFace.h"
class GModel;
class GRegion;
@@ -55,10 +56,11 @@ int MeshTransfiniteVolume(GRegion *gr);
int SubdivideExtrudedMesh(GModel *m);
void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces);
-typedef std::multimap<MVertex*, std::pair<MTriangle*, GFace*> > fs_cont ;
+typedef std::map<MFace,GFace*,Less_Face > fs_cont ;
typedef std::multimap<MVertex*, std::pair<MLine*, GEdge*> > es_cont ;
GFace* findInFaceSearchStructure(MVertex *p1, MVertex *p2, MVertex *p3,
const fs_cont &search);
+GFace* findInFaceSearchStructure(const MFace &f, const fs_cont &search);
GEdge* findInEdgeSearchStructure(MVertex *p1, MVertex *p2, const es_cont &search);
bool buildFaceSearchStructure(GModel *model, fs_cont &search);
bool buildEdgeSearchStructure(GModel *model, es_cont &search);
diff --git a/Mesh/periodical.cpp b/Mesh/periodical.cpp
index 6028889d46..9e80a7accc 100644
--- a/Mesh/periodical.cpp
+++ b/Mesh/periodical.cpp
@@ -185,8 +185,13 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
delta = 0;
container contA(min_x-delta,max_x+delta,min_y-delta,max_y+delta,min_z-delta,max_z+delta,6,6,6,true,true,true,vertices.size());
+ //container contA(min_x-delta,max_x+delta,min_y-delta,max_y+delta,min_z-delta,max_z+delta,6,6,6,false,false,false,vertices.size());
container_poly contB(min_x-delta,max_x+delta,min_y-delta,max_y+delta,min_z-delta,max_z+delta,6,6,6,true,true,true,vertices.size());
+ // wall_cylinder cyl(.5,.5,-3,0,0,6,.4);
+ // contA.add_wall(cyl);
+
+
for(i=0;i<vertices.size();i++){
if(radical==0){
contA.put(i,vertices[i].x(),vertices[i].y(),vertices[i].z());
diff --git a/Numeric/HilbertCurve.cpp b/Numeric/HilbertCurve.cpp
index 4987c5ce33..afd34460c2 100644
--- a/Numeric/HilbertCurve.cpp
+++ b/Numeric/HilbertCurve.cpp
@@ -25,7 +25,7 @@ struct HilbertSort
double BoundingBoxXmin, double BoundingBoxXmax,
double BoundingBoxYmin, double BoundingBoxYmax,
double BoundingBoxZmin, double BoundingBoxZmax, int depth);
- HilbertSort (int m = 0, int l=1) : maxDepth(m),Limit(l)
+ HilbertSort (int m = 0, int l=2) : maxDepth(m),Limit(l)
{
ComputeGrayCode(3);
}
@@ -54,7 +54,7 @@ struct HilbertSort
bbox *= 1.01;
MVertex**pv = &v[0];
int depth;
- MultiscaleSortHilbert(pv, v.size(), 128, 0.125,&depth);
+ MultiscaleSortHilbert(pv, v.size(), 64, 0.125,&depth);
}
};
diff --git a/Plugin/CMakeLists.txt b/Plugin/CMakeLists.txt
index 2545e44f89..e68d68267e 100644
--- a/Plugin/CMakeLists.txt
+++ b/Plugin/CMakeLists.txt
@@ -25,6 +25,7 @@ set(SRC
HomologyComputation.cpp HomologyPostProcessing.cpp
Distance.cpp ExtractEdges.cpp NearestNeighbor.cpp
AnalyseCurvedMesh.cpp
+ CurvedBndDist.cpp
FieldFromAmplitudePhase.cpp
Bubbles.cpp NearToFarField.cpp
DiscretizationError.cpp
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index a444dc6374..2aa18ba377 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -20,6 +20,7 @@
#include "CutBox.h"
#include "Skin.h"
#include "AnalyseCurvedMesh.h"
+#include "CurvedBndDist.h"
#include "MathEval.h"
#include "ExtractElements.h"
#include "SimplePartition.h"
@@ -185,6 +186,8 @@ void PluginManager::registerDefaultPlugins()
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("AnalyseCurvedMesh", GMSH_RegisterAnalyseCurvedMeshPlugin()));
#endif
+ allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
+ ("CurvedBndDist", GMSH_RegisterCurvedBndDistPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("ModifyComponent", GMSH_RegisterModifyComponentPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
diff --git a/benchmarks/boundaryLayers/lann_ver1.geo b/benchmarks/boundaryLayers/lann_ver1.geo
index 378850f880..0f223d9a8b 100644
--- a/benchmarks/boundaryLayers/lann_ver1.geo
+++ b/benchmarks/boundaryLayers/lann_ver1.geo
@@ -31,3 +31,4 @@ Characteristic Length {5, 6, 7, 3, 4} = 10;
Surface Loop(129) = {1,2,3,4,5,6};
Volume(1000) = {129};
+Field[1].Quads = 1;
diff --git a/contrib/HighOrderMeshOptimizer/CMakeLists.txt b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
index 3147c7291a..f7311c8442 100644
--- a/contrib/HighOrderMeshOptimizer/CMakeLists.txt
+++ b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SRC
OptHomMesh.cpp
OptHOM.cpp
OptHomRun.cpp
+ OptHomIntegralBoundaryDist.cpp
ParamCoord.cpp
SuperEl.cpp
OptHomElastic.cpp
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 05d88b3be1..fd0e6d1fdd 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -106,6 +106,29 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
return true;
}
+bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gradObj)
+{
+ maxDistCAD = 0.0;
+
+ std::vector<double> gradF;
+ double DISTANCE = 0.0;
+ for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
+ double f;
+ if (mesh.bndDistAndGradients(iEl, f, gradF, geomTol)) {
+ maxDistCAD = std::max(maxDistCAD,f);
+ DISTANCE += f;
+ Obj += f * factor;
+ for (size_t i = 0; i < mesh.nPCEl(iEl); ++i){
+ gradObj[mesh.indPCEl(iEl, i)] += gradF[i] * factor;
+ // printf("gradf[%d] = %12.5E\n",i,gradF[i]*factor);
+ }
+ }
+ }
+ // printf("DIST = %12.5E\n",DISTANCE);
+ return true;
+
+}
+
bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
{
@@ -132,7 +155,7 @@ bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
}
// Contribution of the vertex distance to the objective function value and
-// gradients (2D version)
+// gradients
bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj,
alglib::real_1d_array &gradObj)
{
@@ -169,7 +192,10 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj,
addDistObjGrad(lambda, lambda2, Obj, gradObj);
if(_optimizeMetricMin)
addMetricMinObjGrad(Obj, gradObj);
- if ((minJac > barrier_min) && (maxJac < barrier_max || !_optimizeBarrierMax)) {
+ if(_optimizeCAD)
+ addBndObjGrad(lambda3, Obj, gradObj);
+ // printf("maxDistCAD = %12.5E distMax = %12.5E Obj %12.5E\n",maxDistCAD,distance_max,Obj);
+ if ((minJac > barrier_min) && (maxJac < barrier_max || !_optimizeBarrierMax) && (maxDistCAD < distance_max|| !_optimizeCAD) ) {
Msg::Info("Reached %s (%g %g) requirements, setting null gradient",
_optimizeMetricMin ? "svd" : "jacobian", minJac, maxJac);
Obj = 0.;
@@ -268,11 +294,12 @@ void OptHOM::calcScale(alglib::real_1d_array &scale)
void OptHOM::OptimPass(alglib::real_1d_array &x,
const alglib::real_1d_array &initGradObj, int itMax)
{
+
static const double EPSG = 0.;
static const double EPSF = 0.;
static const double EPSX = 0.;
static int OPTMETHOD = 1;
-
+
Msg::Info("--- Optimization pass with initial jac. range (%g, %g), jacBar = %g",
minJac, maxJac, jacBar);
@@ -332,20 +359,24 @@ void OptHOM::OptimPass(alglib::real_1d_array &x,
}
}
-int OptHOM::optimize(double weightFixed, double weightFree, double b_min,
+int OptHOM::optimize(double weightFixed, double weightFree, double weightCAD, double b_min,
double b_max, bool optimizeMetricMin, int pInt,
- int itMax, int optPassMax)
+ int itMax, int optPassMax, int optCAD, double distanceMax, double tolerance)
{
barrier_min = b_min;
barrier_max = b_max;
+ distance_max = distanceMax;
progressInterv = pInt;
// powM = 4;
// powP = 3;
_optimizeMetricMin = optimizeMetricMin;
+ _optimizeCAD = optCAD;
// Set weights & length scale for non-dimensionalization
lambda = weightFixed;
lambda2 = weightFree;
+ lambda3 = weightCAD;
+ geomTol = tolerance;
std::vector<double> dSq(mesh.nEl());
mesh.distSqToStraight(dSq);
const double maxDSq = *max_element(dSq.begin(),dSq.end());
@@ -384,13 +415,14 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min,
int ITER = 0;
bool minJacOK = true;
- while (minJac < barrier_min) {
+ while (minJac < barrier_min || (maxDistCAD > distance_max && _optimizeCAD)) {
const double startMinJac = minJac;
OptimPass(x, gradObj, itMax);
recalcJacDist();
jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
+ if (_optimizeCAD) jacBar = std::min(jacBar,barrier_min);
if (ITER ++ > optPassMax) {
- minJacOK = (minJac > barrier_min);
+ minJacOK = (minJac > barrier_min && (maxDistCAD < distance_max || !_optimizeCAD));
break;
}
if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index 098d829132..7ca9db1572 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -50,8 +50,8 @@ public:
// are in the range; returns 0 if the mesh is valid (all jacobians positive,
// JMIN > 0) but JMIN < barrier_min || JMAX > barrier_max; returns -1 if the
// mesh is invalid : some jacobians cannot be made positive
- int optimize(double lambda, double lambda2, double barrier_min, double barrier_max,
- bool optimizeMetricMin, int pInt, int itMax, int optPassMax);
+ int optimize(double lambda, double lambda2, double lambda3, double barrier_min, double barrier_max,
+ bool optimizeMetricMin, int pInt, int itMax, int optPassMax, int optimizeCAD, double optCADDistMax, double tolerance);
void recalcJacDist();
inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
void updateMesh(const alglib::real_1d_array &x);
@@ -59,17 +59,20 @@ public:
alglib::real_1d_array &gradObj);
void printProgress(const alglib::real_1d_array &x, double Obj);
- double barrier_min, barrier_max;
+ double barrier_min, barrier_max, distance_max, geomTol;
private:
- double lambda, lambda2, jacBar, invLengthScaleSq;
+ double lambda, lambda2, lambda3, jacBar, invLengthScaleSq;
int iter, progressInterv; // Current iteration, interval of iterations for reporting
bool _optimizeMetricMin;
double initObj, initMaxDist, initAvgDist; // Values for reporting
- double minJac, maxJac, maxDist, avgDist; // Values for reporting
+ double minJac, maxJac, maxDist, maxDistCAD, avgDist; // Values for reporting
bool _optimizeBarrierMax; // false : only moving barrier min;
// true : fixed barrier min + moving barrier max
+ bool _optimizeCAD; // false : do not minimize the distance between mesh and CAD
+ // true : minimize the distance between mesh and CAD
bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
+ bool addBndObjGrad(double Fact, double &Obj, alglib::real_1d_array &gradObj);
bool addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj);
bool addDistObjGrad(double Fact, double Fact2, double &Obj,
alglib::real_1d_array &gradObj);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp
new file mode 100644
index 0000000000..16b4c29a14
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp
@@ -0,0 +1,526 @@
+#include "GEdge.h"
+#include "nodalBasis.h"
+#include "SVector3.h"
+#include <algorithm>
+#include <limits>
+#include "OptHomIntegralBoundaryDist.h"
+#include "discreteFrechetDistance.h"
+#include "MLine.h"
+
+
+parametricLineNodalBasis::parametricLineNodalBasis(const nodalBasis &basis, const std::vector<SPoint3> &xyz):
+ _basis(basis), _xyz(xyz) {};
+
+SPoint3 parametricLineNodalBasis::operator()(double xi) const
+{
+ std::vector<double> psi(_xyz.size());
+ SPoint3 p(0, 0, 0);
+ _basis.f(-1 + 2 * xi, 0, 0, &psi[0]);
+ for (size_t j = 0; j < psi.size(); ++j) {
+ p[0] += psi[j] * _xyz[j].x();
+ p[1] += psi[j] * _xyz[j].y();
+ p[2] += psi[j] * _xyz[j].z();
+ }
+ return p;
+}
+
+SVector3 parametricLineNodalBasis::derivative(double xi) const
+{
+ double dpsi[_xyz.size()][3];
+ SVector3 p(0, 0, 0);
+ _basis.df(-1 + 2 * xi, 0, 0, dpsi);
+ for (size_t j = 0; j < _xyz.size(); ++j) {
+ p[0] += dpsi[j][0] * _xyz[j].x();
+ p[1] += dpsi[j][0] * _xyz[j].y();
+ p[2] += dpsi[j][0] * _xyz[j].z();
+ }
+ return p;
+}
+
+SVector3 parametricLineNodalBasis::secondDerivative(double xi) const
+{
+ double ddpsi[_xyz.size()][3][3];
+ SVector3 p(0, 0, 0);
+ _basis.ddf(-1 + 2 * xi, 0, 0, ddpsi);
+ for (size_t j = 0; j < _xyz.size(); ++j) {
+ p[0] += ddpsi[j][0][0] * _xyz[j].x();
+ p[1] += ddpsi[j][0][0] * _xyz[j].y();
+ p[2] += ddpsi[j][0][0] * _xyz[j].z();
+ }
+ return p;
+}
+
+
+
+parametricLineGEdge::parametricLineGEdge(const GEdge *edge, double t0, double t1):
+ _edge(edge), _t0(t0), _t1(t1) {}
+
+SPoint3 parametricLineGEdge::operator()(double xi) const
+{
+ GPoint gp = _edge->point(_t0 + (_t1 - _t0) * xi);
+ return SPoint3 (gp.x(), gp.y(), gp.z());
+}
+SVector3 parametricLineGEdge::derivative(double xi) const
+{
+ return _edge->firstDer(_t0 + (_t1 - _t0) * xi);
+}
+
+SVector3 parametricLineGEdge::secondDerivative(double xi) const
+{
+ return _edge->secondDer(_t0 + (_t1 - _t0) * xi);
+}
+
+static void oversample (std::vector<SPoint3> &s, double tol){
+ std::vector<SPoint3> t;
+
+ for (unsigned int i=1;i<s.size();i++){
+ SPoint3 p0 = s[i-1];
+ SPoint3 p1 = s[i];
+ double d = p0.distance(p1);
+ int N = (int) (d / tol);
+ // printf("N = %d %g %g\n",N,d,tol);
+ t.push_back(p0);
+ for (int j=1;j<N;j++){
+ const double xi = (double) j/ N;
+ t.push_back(p0 + (p1-p0)*xi);
+ }
+ }
+ t.push_back(s[s.size()-1]);
+ s = t;
+}
+
+// FAST IMPLEMENTATION OF DISCRETE UNIDIRECTIONAL HAUSDORFF DISTANCE
+double computeBndDistH(GEdge *edge, std::vector<double> & params, // the model edge
+ const std::vector<MVertex*> &vs,
+ const nodalBasis &basis, const std::vector<SPoint3> &xyz,
+ const double tolerance) // the mesh edge
+{
+ if (edge->geomType() == GEntity::Line)return 0.0;
+ std::vector<SPoint3> dpts;
+ std::vector<double> ts;
+ std::vector<MVertex*> hov;
+ for (unsigned int i=2;i<vs.size();i++)hov.push_back(vs[i]);
+ MLineN l (vs[0],vs[1],hov);
+ l.discretize (tolerance,dpts,ts);
+ oversample(dpts,tolerance);
+ double maxDist = 0.0;
+ for (unsigned int i = 0; i<dpts.size(); i++){
+ maxDist = std::max (maxDist,dpts[i].distance(edge->closestPoint(dpts[i],tolerance)));
+ }
+ return maxDist;
+}
+
+SVector3 parametricLine::curvature(double xi) const
+{
+ SVector3 xp = derivative(xi);
+ SVector3 xpp = secondDerivative(xi);
+ const double nxp = xp.norm();
+ const double onxp = 1./nxp;
+ SVector3 c = (onxp*onxp*onxp)*(xpp*nxp-xp*dot(xp,xpp)*onxp);
+ return c;
+}
+
+double parametricLine::frechetDistance(const parametricLine &l,
+ SPoint3 &p1, SPoint3 &p2,
+ double tol) const
+{
+ std::vector<SPoint3> dpts1,dpts2;
+ std::vector<double> ts1,ts2;
+ discretize (dpts1,ts1,tol);
+ l.discretize (dpts2,ts2,tol);
+ // printf("discretizing gives %d %d points\n",dpts1.size(),dpts2.size());
+ oversample(dpts1,tol);
+ oversample(dpts2,tol);
+ // printf("after oversaplinf an discretizing gives %d %d points\n",dpts1.size(),dpts2.size());
+ return discreteFrechetDistance(dpts1,dpts2);
+}
+double parametricLine::hausdorffDistance(const parametricLine &l, SPoint3 &p1, SPoint3 &p2,
+ double tolerance) const
+{
+ std::vector<SPoint3> dpts1,dpts2;
+ std::vector<double> ts1,ts2;
+ discretize (dpts1,ts1,tolerance);
+ l.discretize (dpts2,ts2,tolerance);
+
+ // oversample(dpts1,tolerance);
+ // oversample(dpts2,tolerance);
+ // printf("coucou4 %d %d points\n",dpts1.size(),dpts2.size());
+ double h1 = 0.0;
+ int I1=0,J1=0;
+ int I2=0,J2=0;
+ for (unsigned int i = 0; i<dpts1.size();i++){
+ double hl = 1.e22;
+ int JLOC = 0;
+ for (unsigned int j = 0; j<dpts2.size();j++){
+ double H = dpts1[i].distance(dpts2[j]);
+ if (hl < H){
+ hl = H;
+ JLOC = j;
+ }
+ }
+ if (hl > h1){
+ h1 = hl;
+ J1 = JLOC;
+ I1 = i;
+ }
+ }
+ double h2 = 0.0;
+ for (unsigned int i = 0; i<dpts2.size();i++){
+ double hl = 1.e22;
+ int JLOC = 0;
+ for (unsigned int j = 0; j<dpts1.size();j++){
+ double H = dpts1[j].distance(dpts2[i]);
+ if (hl < H){
+ hl = H;
+ JLOC = j;
+ }
+ }
+ if (hl > h2){
+ h2 = hl;
+ J2 = JLOC;
+ I2 = i;
+ }
+ }
+ if (h1 > h2){
+ p1 = dpts1[I1];
+ p2 = dpts2[J1];
+ return h1;
+ }
+ else{
+ p1 = dpts2[I2];
+ p2 = dpts1[J2];
+ return h2;
+ }
+}
+
+
+// DISCRETE FRECHET DISTANCE
+double computeBndDistF(GEdge *edge,
+ std::vector<double> & params, // the model edge
+ const nodalBasis &basis,
+ const std::vector<SPoint3> &xyz,
+ const double tolerance) // the mesh edge
+{
+ if (edge->geomType() == GEntity::Line)return 0.0;
+ parametricLineGEdge l1 = parametricLineGEdge(edge, params[0],params[1]);
+ parametricLineNodalBasis l2 = parametricLineNodalBasis(basis, xyz);
+ SPoint3 p1,p2;
+ return l1.frechetDistance(l2,p1,p2,tolerance);
+}
+
+
+// GMSH's DISTANCE
+/*
+ */
+double computeBndDistGb(GEdge *edge, std::vector<double> & params, // the model edge
+ const nodalBasis &basis, const std::vector<SPoint3> &xyz, double tolerance) // the mesh edge
+{
+ parametricLineGEdge l1 = parametricLineGEdge(edge, params[0], params[1]);
+ parametricLineNodalBasis l2 = parametricLineNodalBasis(basis, xyz);
+ const unsigned int N = 20;
+ SPoint3 P1[N];
+ SPoint3 P2[N];
+ double D = 0.0;
+ for (unsigned int i=0;i<N;i++){
+ const double _x2 = (double)i / (N-1);
+ P1[i] = l1(_x2);
+ P2[i] = l2(_x2);
+ }
+ double L = 0.0;
+ for (unsigned int i=0;i<N-1;i++){
+ SPoint3 p11 = P1[i];
+ SPoint3 p12 = P1[i+1];
+ SPoint3 p21 = P2[i];
+ SPoint3 p22 = P2[i+1];
+ SVector3 v1 (p21,p11);
+ SVector3 v2 (p12,p11);
+ SVector3 v12 (p22,p11);
+ SVector3 vl (p22,p21);
+ SVector3 x1 = crossprod(v12,v1);
+ SVector3 x2 = crossprod(v12,v2);
+ D += 0.5*(x1.norm()+x2.norm());
+ L += vl.norm();
+ }
+ return D;
+}
+
+double computeBndDistG_(GEdge *edge, std::vector<double> & p, // the model edge
+ const nodalBasis &basis, const std::vector<SPoint3> &xyz,
+ const unsigned int N) // the mesh edge
+{
+ std::vector<int> o;
+ o.push_back(0);
+ for (unsigned int i=2; i < p.size();i++)o.push_back(i);
+ o.push_back(1);
+
+ // printf("computing diustance with tolerance %g\n",tolerac);
+
+
+
+ double D = 0.0;
+ const double U0 = basis.points(0,0);
+ const double U1 = basis.points(1,0);
+ for (int i = 0 ; i < basis.order ; i++) {
+ const double u0 = basis.points(o[i],0);
+ const double u1 = basis.points(o[i+1],0);
+
+ // U0 ----u0-----u1-----U1
+
+ const double t0 = p[o[i]];
+ const double t1 = p[o[i+1]];
+ parametricLineGEdge l1 = parametricLineGEdge(edge, t0, t1);
+ parametricLineNodalBasis l2 = parametricLineNodalBasis(basis, xyz);
+ SPoint3 P1[N];
+ SPoint3 P2[N];
+ for (unsigned int i=0;i<N;i++){
+ const double _x2 = (double)i / (N-1);
+ const double u = u0 + _x2 * (u1-u0);
+ // U0 + uu * (U1 - U0) = u
+ const double uu = (u - U0) / (U1-U0);
+ P1[i] = l1(_x2);
+ P2[i] = l2(uu);
+ }
+ // double L = 0.0;
+ for (unsigned int i=0;i<N-1;i++){
+ SPoint3 p11 = P1[i];
+ SPoint3 p12 = P1[i+1];
+ SPoint3 p21 = P2[i];
+ SPoint3 p22 = P2[i+1];
+ SVector3 v1 (p21,p11);
+ SVector3 v2 (p12,p11);
+ SVector3 v12 (p22,p11);
+ // SVector3 vl (p22,p21);
+ SVector3 x1 = crossprod(v12,v1);
+ SVector3 x2 = crossprod(v12,v2);
+ D += 0.5*(x1.norm()+x2.norm());
+ // L += vl.norm();
+ }
+ }
+ return D;
+}
+
+double computeBndDistG(GEdge *edge, std::vector<double> & p, // the model edge
+ const nodalBasis &basis, const std::vector<SPoint3> &xyz,
+ double tolerance) // the mesh edge
+{
+ int N = 4;
+ double d = computeBndDistG_(edge, p, basis, xyz, N);
+
+ // printf("GO !!\n");
+ while (1){
+ N *= 2;
+ double dp = computeBndDistG_(edge, p, basis, xyz, N);
+ // printf("%12.5E %12.5E %12.5E %12.5E\n",d,dp,fabs(d - dp),tolerance);
+ if (fabs(d - dp) < tolerance) // Richardson with assumed linear convergence ...
+ return dp;
+ d = dp;
+ }
+}
+
+void parametricLine::recur_discretize(const double &t1, const double &t2,
+ const SPoint3 &p1, const SPoint3 &p2,
+ std::vector<SPoint3> &dpts,
+ std::vector<double> &ts,
+ double Prec, int depth) const
+{
+ double t = 0.5 * (t2 + t1);
+ SPoint3 p = (*this)(t);
+ SVector3 dx (p,(p1+p2)*0.5);
+ // printf("%g %g -- %g %g dist %12.5E %12.5E\n",p1.x(),p1.y(),p2.x(),p2.y(),dx.norm(),Prec);
+ if ((depth > 20 && dx.norm() < Prec) || depth > 45){
+ dpts.push_back(p);
+ ts.push_back(t);
+ dpts.push_back(p2);
+ ts.push_back(t2);
+ }
+ else {
+ recur_discretize (t1,t,p1,p,dpts,ts,Prec,depth+1);
+ recur_discretize (t,t2,p,p2,dpts,ts,Prec,depth+1);
+ }
+}
+
+void parametricLine::discretize(std::vector<SPoint3> &dpts,
+ std::vector<double> &ts,
+ double Prec, double t0, double t1) const
+{
+ dpts.push_back((*this)(t0));
+ ts.push_back(t0);
+ recur_discretize (t0,t1,dpts[0],(*this)(t1),dpts,ts,Prec,0);
+
+ // printf("discretizing from %g to %g\n",t0,t1);
+ // for (unsigned int i=0;i<ts.size();i++){
+ // printf("%g ",ts[i]);
+ // }
+ // printf("\n");
+
+}
+
+double trapeze (SPoint3 &p1, SPoint3 &p2){
+ return (p2.x() - p1.x())*(p1.y()+p2.y())*0.5;
+}
+
+double computeDeviationOfTangents(GEdge *edge,
+ std::vector<double> &p, // parameters of mesh vertices on the model edge
+ const nodalBasis &basis,
+ const std::vector<SPoint3> &xyz) // the mesh edge
+{
+ // parametricLineGEdge l1 = parametricLineGEdge(edge,p[0],p[p.size()-1]);
+ parametricLineNodalBasis l2 = parametricLineNodalBasis(basis, xyz);
+ double deviation = 0;
+ double ddeviation = 0;
+ std::vector<int> o;
+ o.push_back(0);
+ for (unsigned int i=2; i < p.size();i++)o.push_back(i);
+ o.push_back(1);
+
+ SVector3 dx (xyz[xyz.size() - 1], xyz[0]);
+
+ for (unsigned int i=0; i<p.size();i++){
+ const double u = basis.points(o[i],0);
+ SVector3 xp = edge->firstDer (p[o[i]]);
+ SVector3 xpp = edge->secondDer (p[o[i]]);
+ const double nxp = xp.norm();
+ const double onxp = 1./nxp;
+ SVector3 c = (onxp*onxp*onxp)*(xpp*nxp-xp*dot(xp,xpp)*onxp);
+
+ SVector3 t_mesh_edge = l2.derivative(0.5*(1+u));
+ SVector3 c2 = l2.curvature(0.5*(1+u));
+ // GPoint p0 = edge->point(p[o[i]]);
+ // SPoint3 p1 = l2 (0.5*(1+u));
+ // printf("%g = %g %g vs %g %g\n",u,p0.x(),p0.y(),p1.x(),p1.y());
+ xp.normalize();
+ t_mesh_edge.normalize();
+ SVector3 diff1 = (dot(xp, t_mesh_edge) > 0) ? xp - t_mesh_edge : xp + t_mesh_edge;
+ SVector3 diff2 = (dot(c, c2) > 0) ? c - c2 : c + c2;
+ //printf("%g %g %g vs %g %g %g diff %g %g %g\n",c.x(),c.y(),c.z(),c2.x(),c2.y(),c2.z(),diff2.x(),diff2.y(),diff2.z());
+ // printf("%g %g %g vs %g %g %g val %g\n",t_model_edge.x(),t_model_edge.y(),t_model_edge.z(),
+ // t_mesh_edge.x(),t_mesh_edge.y(),t_mesh_edge.z(),c.norm());
+ // deviation = std::max(diff1.norm(),deviation);
+ // ddeviation = std::max(diff2.norm(),ddeviation);
+ deviation += diff1.norm();
+ ddeviation += diff2.norm();
+ }
+ const double h = dx.norm();
+ // printf ("%g %g\n",deviation * h,ddeviation * h * h * 0.5);
+ return deviation * h;// + ddeviation * h * h * 0.5;
+}
+
+double computeBndDistAccurateArea(GEdge *edge,
+ std::vector<double> &p, // parameters of mesh vertices on the model edge
+ const nodalBasis &basis,
+ const std::vector<SPoint3> &xyz,
+ double tolerance) // the mesh edge
+{
+ // assume mesh and CAD non intersecting except at mesh vertices
+ // compute distance for "order" sub-polygonal curves
+ double area = 0.0;
+ double length = 0.0;
+ std::vector<int> o;
+ o.push_back(0);
+ for (unsigned int i=2; i < p.size();i++)o.push_back(i);
+ o.push_back(1);
+
+ // printf("computing diustance with tolerance %g\n",tolerac);
+
+ for (int i = 0 ; i < basis.order ; i++) {
+ const double u0 = basis.points(o[i],0);
+ const double u1 = basis.points(o[i+1],0);
+ const double t0 = p[o[i]];
+ const double t1 = p[o[i+1]];
+ parametricLineGEdge l1 = parametricLineGEdge(edge, t0, t1);
+ parametricLineNodalBasis l2 = parametricLineNodalBasis(basis, xyz);
+ std::vector<SPoint3> dpts1,dpts2;
+ std::vector<double> ts1,ts2;
+ l1.discretize (dpts1,ts1,tolerance);
+ l2.discretize (dpts2,ts2,tolerance,0.5*(1+u0),0.5*(1+u1));
+ // printf("discretizing : %g %g and %g %g\n",u0,u1,t0,t1);
+ // simple 2D version
+ double arealocal = 0.0;
+ for (unsigned int j=1;j<dpts1.size(); j++) {
+ length += dpts1[j-1].distance (dpts1[j]);
+ arealocal += trapeze (dpts1[j-1],dpts1[j]);
+ }
+ for (unsigned int j=1;j<dpts2.size(); j++) {
+ arealocal -= trapeze (dpts2[j-1],dpts2[j]);
+ }
+ area += fabs(arealocal);
+ }
+ return area;
+}
+
+
+// INPUT FCT FOR OPTIMIZATION
+double computeBndDistAndGradient(GEdge *edge,
+ std::vector<double> ¶m, // parameters of mesh vertices on the model edge
+ const std::vector<MVertex*> &vs, // vertices
+ const nodalBasis &basis, // what is the FE basis of the edge
+ std::vector<SPoint3> &xyz, // real coordinates of mesh vertices on the model edge
+ std::vector<bool> &onEdge, // tell if a given vertex sits on the model edge and therefore can be movd
+ std::vector<double> &grad,
+ double tolerance)// model tolerance
+{
+ grad.resize(xyz.size());
+ double ref;
+ if (tolerance < 0 ) ref = computeDeviationOfTangents(edge, param,basis, xyz);
+ else ref = computeBndDistG(edge, param,basis, xyz, tolerance);
+ //double ref = computeBndDistAccurateArea(edge, param,basis, xyz, tolerance);
+ double delta = (edge->getUpperBound() - edge->getLowerBound()) * 1e-8;
+ for (size_t i = 0; i < xyz.size(); ++i) {
+ if (!onEdge[i]) {
+ grad[i] = 0;
+ continue;
+ }
+ double p = param[i];
+ double delta = 1e-6;
+ param[i] += delta;
+ xyz[i] = SPoint3(edge->position(param[i]));
+ if (tolerance > 0) grad[i] = (computeBndDistG(edge, param,basis, xyz,tolerance) - ref) / delta;
+ else grad[i] = (computeDeviationOfTangents(edge, param,basis, xyz) - ref) / delta;
+ param[i] = p;
+ xyz[i] = SPoint3(edge->position(param[i]));
+ }
+ return ref;
+}
+
+#include "MElement.h"
+#include "MVertex.h"
+#include "BasisFactory.h"
+double computeBndDist(MElement *element, int distanceDefinition, double tolerance)
+{
+ double dist = 0;
+ const nodalBasis &elbasis = *element->getFunctionSpace();
+ for (int iEdge = 0; iEdge < element->getNumEdges(); ++iEdge) {
+ int clId = elbasis.getClosureId(iEdge, 1);
+ const std::vector<int> &closure = elbasis.closures[clId];
+ std::vector<SPoint3> xyz;
+ GEdge *edge = NULL;
+ std::vector<MVertex *> vertices (closure.size());
+ for (size_t i = 0; i < closure.size(); ++i) {
+ MVertex *v = element->getVertex(closure[i]);
+ vertices[i] = v;
+ xyz.push_back(v->point());
+ if ((int)i >= 2 && v->onWhat() && v->onWhat()->dim() == 1) {
+ edge = v->onWhat()->cast2Edge();
+ }
+ }
+ if (edge){
+ std::vector<double> params(closure.size());
+ for (size_t i = 0; i < closure.size(); ++i) {
+ reparamMeshVertexOnEdge(element->getVertex(closure[i]), edge, params[i]);
+ }
+ if (distanceDefinition == 1)
+ dist = std::max(computeBndDistH(edge, params, vertices, *BasisFactory::getNodalBasis(elbasis.getClosureType(clId)), xyz,tolerance),dist);
+ else if (distanceDefinition == 2)
+ dist = std::max(computeBndDistG(edge, params, *BasisFactory::getNodalBasis(elbasis.getClosureType(clId)), xyz, tolerance),dist);
+ else if (distanceDefinition == 4)
+ dist = std::max(computeBndDistF(edge, params, *BasisFactory::getNodalBasis(elbasis.getClosureType(clId)), xyz, tolerance),dist);
+ else if (distanceDefinition == 5)
+ dist = std::max(computeBndDistAccurateArea(edge, params,*BasisFactory::getNodalBasis(elbasis.getClosureType(clId)), xyz, tolerance),dist);
+ else if (distanceDefinition == 6)
+ dist = std::max(computeDeviationOfTangents(edge, params,*BasisFactory::getNodalBasis(elbasis.getClosureType(clId)), xyz),dist);
+ else
+ Msg::Fatal("unknown distance definition %d. Choose 1 for Hausdorff and 2 for Area/Length 4 for Discrete Frechet",distanceDefinition);
+ }
+ }
+ return dist;
+}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
new file mode 100644
index 0000000000..a0e22b57fe
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
@@ -0,0 +1,56 @@
+#ifndef OPT_HOM_INTEGRAL_BOUNDARY_DIST_H
+#define OPT_HOM_INTEGRAL_BOUNDARY_DIST_H
+#include <vector>
+class GEdge;
+class nodalBasis;
+class SPoint3;
+class MElement;
+class MVertex;
+double computeBndDist(MElement *element, int distanceDefinition, double tolerance);
+double computeBndDistAndGradient(GEdge *edge,
+ std::vector<double> ¶m,
+ const std::vector<MVertex*> &vs,
+ const nodalBasis &basis, std::vector<SPoint3> &xyz,
+ std::vector<bool> &onEdge, std::vector<double> &grad, double tolerance);
+
+class parametricLine {
+ public :
+ virtual SPoint3 operator()(double xi) const = 0;
+ virtual SVector3 derivative(double xi) const = 0;
+ virtual SVector3 secondDerivative(double xi) const = 0;
+ SVector3 curvature(double xi) const;
+ void discretize(std::vector<SPoint3> &dpts,std::vector<double> &ts,double Prec,
+ double t0 = 0.0, double t1 = 1.0) const;
+ void recur_discretize(const double &t1, const double &t2,
+ const SPoint3 &p1, const SPoint3 &p2,
+ std::vector<SPoint3> &dpts,
+ std::vector<double> &ts,
+ double Prec, int depth) const;
+ double frechetDistance(const parametricLine &l, SPoint3 &p1, SPoint3 &p2, double tol = 1.e-6) const;
+ double hausdorffDistance(const parametricLine &l, SPoint3 &p1, SPoint3 &p2, double tol = 1.e-6) const;
+};
+
+class parametricLineNodalBasis : public parametricLine
+{
+ const nodalBasis &_basis;
+ const std::vector<SPoint3> &_xyz;
+ double _u0, _u1;
+ public :
+ parametricLineNodalBasis(const nodalBasis &basis, const std::vector<SPoint3> &xyz);
+ virtual SPoint3 operator()(double xi) const;
+ virtual SVector3 derivative(double xi) const;
+ virtual SVector3 secondDerivative(double xi) const;
+};
+
+class parametricLineGEdge : public parametricLine
+{
+ const GEdge *_edge;
+ double _t0, _t1;
+ public :
+ parametricLineGEdge(const GEdge *edge, double t0, double t1);
+ virtual SPoint3 operator()(double xi) const;
+ virtual SVector3 derivative(double xi) const;
+ virtual SVector3 secondDerivative(double xi) const;
+};
+
+#endif
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 3d10e1aca3..5c7614c6d2 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -34,6 +34,8 @@
#include "MTetrahedron.h"
#include "ParamCoord.h"
#include "OptHomMesh.h"
+#include "BasisFactory.h"
+#include "OptHomIntegralBoundaryDist.h"
Mesh::Mesh(const std::map<MElement*,GEntity*> &element2entity,
const std::set<MElement*> &els, std::set<MVertex*> &toFix,
@@ -292,6 +294,66 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda,
}
}
+bool Mesh::bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps)
+{
+ MElement *element = _el[iEl];
+ f = 0.;
+ // dommage ;-)
+ if (element->getDim() != 2)
+ return false;
+
+ int currentId = 0;
+ std::vector<int> vertex2param(element->getNumVertices());
+ for (size_t i = 0; i < element->getNumVertices(); ++i) {
+ if (_el2FV[iEl][i] >= 0) {
+ vertex2param[i] = currentId;
+ currentId += _nPCFV[_el2FV[iEl][i]];
+ }
+ else
+ vertex2param[i] = -1;
+ }
+ gradF.clear();
+ gradF.resize(currentId, 0.);
+
+ const nodalBasis &elbasis = *element->getFunctionSpace();
+ bool edgeFound = false;
+ for (int iEdge = 0; iEdge < element->getNumEdges(); ++iEdge) {
+ int clId = elbasis.getClosureId(iEdge, 1);
+ const std::vector<int> &closure = elbasis.closures[clId];
+ std::vector<MVertex *> vertices;
+ GEdge *edge = NULL;
+ for (size_t i = 0; i < closure.size(); ++i) {
+ MVertex *v = element->getVertex(closure[i]);
+ vertices.push_back(v);
+ // only valid in 2D
+ if ((int)i >= 2 && v->onWhat() && v->onWhat()->dim() == 1) {
+ edge = v->onWhat()->cast2Edge();
+ }
+ }
+ if (edge) {
+ edgeFound = true;
+ std::vector<double> localgrad;
+ std::vector<SPoint3> nodes(closure.size());
+ std::vector<double> params(closure.size());
+ std::vector<bool> onedge(closure.size());
+ for (size_t i = 0; i < closure.size(); ++i) {
+ nodes[i] = _xyz[_el2V[iEl][closure[i]]];
+ onedge[i] = element->getVertex(closure[i])->onWhat() == edge && _el2FV[iEl][closure[i]] >= 0;
+ if (onedge[i]) {
+ params[i] = _uvw[_el2FV[iEl][closure[i]]].x();
+ }else
+ reparamMeshVertexOnEdge(element->getVertex(closure[i]), edge, params[i]);
+ }
+ f += computeBndDistAndGradient(edge, params, vertices, *BasisFactory::getNodalBasis(elbasis.getClosureType(clId)), nodes, onedge, localgrad, eps);
+ for (size_t i = 0; i < closure.size(); ++i) {
+ if (onedge[i])
+ gradF[vertex2param[closure[i]]] += localgrad[i];
+ }
+ }
+ }
+ return edgeFound;
+}
+
void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
std::vector<double> &gSJ)
{
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 8a6d6bdb2e..f1f48c9fed 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -57,6 +57,7 @@ public:
void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
+ bool bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps);
inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; }
inline double distSq(int iFV);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 821a8ead4b..ab0a0f7921 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -47,6 +47,35 @@
#if defined(HAVE_BFGS)
+double distMaxStraight(MElement *el)
+{
+ const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
+ const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
+ int nV = lagrange->points.size1();
+ int nV1 = lagrange1->points.size1();
+ int dim = lagrange1->dimension;
+ SPoint3 sxyz[256];
+ for (int i = 0; i < nV1; ++i) {
+ sxyz[i] = el->getVertex(i)->point();
+ }
+ for (int i = nV1; i < nV; ++i) {
+ double f[256];
+ lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1),
+ dim < 3 ? 0 : lagrange->points(i, 2), f);
+ for (int j = 0; j < nV1; ++j)
+ sxyz[i] += sxyz[j] * f[j];
+ }
+
+ double maxdx = 0.0;
+ for (int iV = nV1; iV < nV; iV++) {
+ SVector3 d = el->getVertex(iV)->point()-sxyz[iV];
+ double dx = d.norm();
+ if (dx > maxdx) maxdx = dx;
+ }
+
+ return maxdx;
+}
+
void exportMeshToDassault(GModel *gm, const std::string &fn, int dim)
{
FILE *f = fopen(fn.c_str(),"w");
@@ -221,7 +250,8 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
std::vector<std::set<MElement*> > primBlobs;
primBlobs.reserve(badElements.size());
for (std::set<MElement*>::const_iterator it = badElements.begin(); it != badElements.end(); ++it) {
- const int minLayers = ((*it)->getDim() == 3) ? 1 : 0;
+ //const int minLayers = ((*it)->getDim() == 3) ? 1 : 0;
+ const int minLayers = 3;
primBlobs.push_back(getSurroundingBlob(*it, depth, vertex2elements,
distFactor, minLayers, optPrimSurfMesh));
}
@@ -308,12 +338,12 @@ static void optimizeConnectedBlobs
if (temp.mesh.nPC() == 0)
Msg::Info("Blob %i has no degree of freedom, skipping", i+1);
else
- success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
- p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
+ success = temp.optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN,
+ p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax, p.discrTolerance);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
Msg::Info("Jacobian optimization succeed, starting svd optimization");
- success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
- true, samples, p.itMax, p.optPassMax);
+ success = temp.optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
+ true, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax,p.discrTolerance);
}
double minJac, maxJac, distMaxBND, distAvgBND;
temp.recalcJacDist();
@@ -321,11 +351,11 @@ static void optimizeConnectedBlobs
p.minJac = std::min(p.minJac,minJac);
p.maxJac = std::max(p.maxJac,maxJac);
temp.mesh.updateGEntityPositions();
- if (success <= 0) {
+ //if (success <= 0) {
std::ostringstream ossI2;
ossI2 << "final_ITER_" << i << ".msh";
temp.mesh.writeMSH(ossI2.str().c_str());
- }
+ //}
//#pragma omp critical
p.SUCCESS = std::min(p.SUCCESS, success);
}
@@ -550,12 +580,12 @@ static void optimizeOneByOne
std::ostringstream ossI1;
ossI1 << "initial_blob-" << iBadEl << ".msh";
opt->mesh.writeMSH(ossI1.str().c_str());
- success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
- p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
+ success = opt->optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN,
+ p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax,p.discrTolerance);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
Msg::Info("Jacobian optimization succeed, starting svd optimization");
- success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
- p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
+ success = opt->optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN_METRIC,
+ p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax,p.discrTolerance);
}
// Measure min and max Jac., update mesh
@@ -594,6 +624,27 @@ static void optimizeOneByOne
#endif
+#include "OptHomIntegralBoundaryDist.h"
+double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition, double tolerance){
+ double maxd = 0.0;
+ double sum = 0.0;
+ int NUM;
+ for (int iEl = 0; iEl < ge->getNumMeshElements();iEl++) {
+ MElement *el = ge->getMeshElement(iEl);
+ if (ge->dim() == el->getDim()){
+ const double DISTE =computeBndDist(el,distanceDefinition, tolerance);
+ if (DISTE != 0.0){
+ NUM++;
+ // if(distanceDefinition == 1)printf("%d %12.5E\n",iEl,DISTE);
+ maxd = std::max(maxd,DISTE);
+ sum += DISTE;
+ }
+ }
+ }
+ if (distanceDefinition == 2)return sum/NUM;
+ return maxd;
+}
+
void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
{
#if defined(HAVE_BFGS)
@@ -610,6 +661,7 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
std::map<MVertex*, std::vector<MElement *> > vertex2elements;
std::map<MElement*,GEntity*> element2entity;
std::set<MElement*> badasses;
+ double maxdist = 0;
for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
GEntity* &entity = entities[iEnt];
if (entity->dim() != p.dim || (p.onlyVisible && !entity->getVisibility())) continue;
@@ -620,13 +672,22 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
double jmin, jmax;
MElement *el = entity->getMeshElement(iEl);
+ const double DISTE =computeBndDist(el,2,fabs(p.discrTolerance));
+ // printf("Element %d Distance %12.5E\n",iEl,DISTE);
+ maxdist = std::max(DISTE, maxdist);
if (el->getDim() == p.dim) {
- el->scaledJacRange(jmin, jmax, p.optPrimSurfMesh ? entity : 0);
- if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
- if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) badasses.insert(el);
+ if (p.optCAD && DISTE > p.optCADDistMax)
+ badasses.insert(el);
+
+ el->scaledJacRange(jmin, jmax, p.optPrimSurfMesh ? entity : 0);
+ if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
+ if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX){
+ badasses.insert(el);
+ }
}
}
}
+ printf("maxdist = %g badasses size = %lu\n", maxdist, badasses.size());
if (p.strategy == 0)
optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, false);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index e934f01ad3..f1d172f693 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -54,6 +54,10 @@ struct OptHomParameters {
int adaptBlobLayerFact; // Growth factor in number of layers for blob adaptation
double adaptBlobDistFact; // Growth factor in distance factor for blob adaptation
bool optPrimSurfMesh; // Enable optimisation of p1 surface meshes
+ bool optCAD;//Enable optimization of mesh vertices positions for geometrical fitting
+ double optCADWeight;//Weight
+ double optCADDistMax;//Maximum allowed distance from the CAD
+ double discrTolerance;
// OUTPUT ------>
int SUCCESS ; // 0 --> success , 1 --> Not converged
@@ -64,11 +68,14 @@ struct OptHomParameters {
: BARRIER_MIN_METRIC(-1.), BARRIER_MIN(0.1), BARRIER_MAX(2.0), weightFixed(1000.),
weightFree (1.), nbLayers (6) , dim(3) , itMax(300), onlyVisible(true),
distanceFactor(12), fixBndNodes(false), strategy(0), maxAdaptBlob(3),
- adaptBlobLayerFact(2.), adaptBlobDistFact(2.), optPrimSurfMesh(false)
+ adaptBlobLayerFact(2.), adaptBlobDistFact(2.), optPrimSurfMesh(false),optCAD(false),
+ optCADWeight(1000.),optCADDistMax(1.e22),discrTolerance(1.e-4)
{
}
};
void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p);
+// distanceDefinition 1) Hausdorff 2) Area/Length 3) Frechet (not done)
+double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition,double tolerance) ;
#endif
--
GitLab