diff --git a/CMakeLists.txt b/CMakeLists.txt
index d23017c98d70f1501841195ce9dfa2a7195936cd..b0341479028e829cb389b2c7dc87385435b05048 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -148,6 +148,8 @@ set(GMSH_API
contrib/HighOrderMeshOptimizer/OptHOM.h contrib/HighOrderMeshOptimizer/OptHomMesh.h
contrib/HighOrderMeshOptimizer/OptHomRun.h contrib/HighOrderMeshOptimizer/ParamCoord.h
contrib/HighOrderMeshOptimizer/OptHomFastCurving.h contrib/HighOrderMeshOptimizer/SuperEl.h
+ contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
+ contrib/HighOrderMeshOptimizer/OptHomCADDist.h
contrib/HighOrderMeshOptimizer/OptHomObjContribScaledJac.h
contrib/HighOrderMeshOptimizer/OptHomObjContribMetricMin.h
contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 1b2f0307c019a5d2fbdb3a1e2ccf2fe840e5d2a3..3200e1bd2f70974a1dabd147b80b0ba48766f8ba 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -283,7 +283,7 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
for (unsigned int i=0;i<(*it)->triangles.size(); i++){
// compute the distance from the geometry to the mesh
if(doWeCompute[i]){
- const double d = MFaceGFaceDistance ( (*it)->triangles[i] , *it , &gsfT, &normalsToCAD);
+ const double d = MFaceGFaceDistanceOld((*it)->triangles[i], *it, &gsfT, &normalsToCAD);
dist[(*it)->triangles[i]] = d;
maxDistCAD = std::max(maxDistCAD,d);
distCAD += d * factor;
@@ -309,7 +309,7 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
v->setParameter(0,t_+eps);
v->setXYZ(gp.x(),gp.y(),gp.z());
const double distT = dist[t];
- double deriv = (MFaceGFaceDistance ( t , *it , &gsfT, &normalsToCAD) - distT) /eps;
+ double deriv = (MFaceGFaceDistanceOld(t, *it, &gsfT, &normalsToCAD) - distT) /eps;
v->setXYZ(pp.x(),pp.y(),pp.z());
v->setParameter(0,t_);
gradObj[index] += deriv * factor;
@@ -330,7 +330,7 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
GPoint gp = (*it)->point(uu+eps,vv);
v->setParameter(0,uu+eps);
v->setXYZ(gp.x(),gp.y(),gp.z());
- double deriv = (MFaceGFaceDistance ( t , *it, &gsfT, &normalsToCAD ) - distT) /eps;
+ double deriv = (MFaceGFaceDistanceOld(t, *it, &gsfT, &normalsToCAD) - distT) /eps;
v->setXYZ(pp.x(),pp.y(),pp.z());
v->setParameter(0,uu);
gradObj[index] += deriv * factor;
@@ -338,7 +338,7 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
gp = (*it)->point(uu,vv+eps);
v->setParameter(1,vv+eps);
v->setXYZ(gp.x(),gp.y(),gp.z());
- deriv = (MFaceGFaceDistance ( t , *it, &gsfT, &normalsToCAD ) - distT) /eps;
+ deriv = (MFaceGFaceDistanceOld(t, *it, &gsfT, &normalsToCAD) - distT) /eps;
v->setXYZ(pp.x(),pp.y(),pp.z());
v->setParameter(1,vv);
gradObj[index+1] += deriv * factor;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp b/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
index 4644d77435aa66665fc54eef96b489d76a0fc7ad..5e3d16b0076d037f17173129a13f6c44cd3d9b16 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
@@ -3,12 +3,16 @@
#include "MLine.h"
#include "MTriangle.h"
#include "GModel.h"
+
+#include "GEdge.h"
+#include "JacobianBasis.h"
+#include "BasisFactory.h"
#include "OptHomCADDist.h"
-double MFaceGFaceDistance(MTriangle *t, GFace *gf,
- std::vector<std::vector<SVector3> > *gsfT,
- std::map<MVertex*,SVector3> *normalsToCAD) {
+double MFaceGFaceDistanceOld(MTriangle *t, GFace *gf,
+ std::vector<std::vector<SVector3> > *gsfT,
+ std::map<MVertex*,SVector3> *normalsToCAD) {
const double h = t->maxEdge();
double jac[3][3];
double distFace = 0.0;
@@ -48,7 +52,7 @@ double MFaceGFaceDistance(MTriangle *t, GFace *gf,
}
-double MLineGEdgeDistance (MLine *l, GEdge *ge, FILE *f) {
+double MLineGEdgeDistanceOld(MLine *l, GEdge *ge, FILE *f) {
const nodalBasis &elbasis = *l->getFunctionSpace();
const double h = .25*0.5*distance (l->getVertex (0), l->getVertex (1) ) / (l->getNumVertices()-1);
double jac[3][3];
@@ -84,6 +88,93 @@ double MLineGEdgeDistance (MLine *l, GEdge *ge, FILE *f) {
}
+double distToCAD1D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &tanCAD, double edLength)
+{
+ const int nV = nodesXYZ.size1();
+ const double h = .25*0.5*edLength/(nV-1);
+ fullMatrix<double> dxyzdX(nV, 3);
+ gb->getGradientsFromNodes(nodesXYZ, &dxyzdX, 0, 0);
+ double dist = 0.;
+ for (int i=0; i<nV; i++) {
+ SVector3 tanMesh(dxyzdX(i, 0), dxyzdX(i, 1), dxyzdX(i, 2));
+ tanMesh.normalize();
+ SVector3 diff = (dot(tanCAD[i], tanMesh) > 0) ?
+ tanCAD[i] - tanMesh : tanCAD[i] + tanMesh;
+ dist += diff.norm();
+ }
+ return h*dist;
+}
+
+
+double distToCAD2D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &normCAD)
+{
+ const int nV = nodesXYZ.size1();
+ fullMatrix<double> dxyzdX(nV, 3),dxyzdY(nV, 3);
+ gb->getGradientsFromNodes(nodesXYZ, &dxyzdX, &dxyzdY, 0);
+ double dist = 0.;
+ for (int i=0; i<nV; i++) {
+ const double nz = dxyzdX(i, 0) * dxyzdY(i, 1) - dxyzdX(i, 1) * dxyzdY(i, 0);
+ const double ny = -dxyzdX(i, 0) * dxyzdY(i, 2) + dxyzdX(i, 2) * dxyzdY(i, 0);
+ const double nx = dxyzdX(i, 1) * dxyzdY(i, 2) - dxyzdX(i, 2) * dxyzdY(i, 1);
+ SVector3 normMesh(nx, ny, nz);
+ normMesh.normalize();
+ SVector3 diff = (dot(normCAD[i], normMesh) > 0) ?
+ normCAD[i] - normMesh : normCAD[i] + normMesh;
+ dist += diff.norm();
+ }
+ return dist;
+}
+
+
+double MLineGEdgeDistance(MLine *l, GEdge *ge)
+{
+ const int nV = l->getNumVertices();
+ const GradientBasis *gb = BasisFactory::getGradientBasis(FuncSpaceData(l));
+ const double edLength = l->getLength();
+
+ // Coordinates of vertices
+ fullMatrix<double> nodesXYZ(nV, 3);
+ l->getNodesCoord(nodesXYZ);
+
+ // Tangent to CAD at vertices
+ std::vector<SVector3> tanCAD(nV);
+ for (int i=0; i<nV; i++) {
+ double tCAD;
+ reparamMeshVertexOnEdge(l->getVertex(i), ge, tCAD);
+ tanCAD[i] = ge->firstDer(tCAD);
+ tanCAD[i].normalize();
+ }
+
+ // Compute distance
+ return distToCAD1D(gb, nodesXYZ, tanCAD, edLength);
+}
+
+
+double MFaceGFaceDistance(MElement *el, GFace *gf)
+{
+ const int nV = el->getNumVertices();
+ const GradientBasis *gb = BasisFactory::getGradientBasis(FuncSpaceData(el));
+
+ // Coordinates of vertices
+ fullMatrix<double> nodesXYZ(nV, 3);
+ el->getNodesCoord(nodesXYZ);
+
+ // Normal to CAD at vertices
+ std::vector<SVector3> normCAD(nV);
+ for (int i=0; i<nV; i++) {
+ SPoint2 pCAD;
+ reparamMeshVertexOnFace(el->getVertex(i), gf, pCAD);
+ normCAD[i] = gf->normal(pCAD);
+ normCAD[i].normalize();
+ }
+
+ // Compute distance
+ return distToCAD2D(gb, nodesXYZ, normCAD);
+}
+
+
void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,double> &distances){
std::map<MEdge,double,Less_Edge> dist2Edge;
@@ -134,35 +225,22 @@ void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,doub
double distanceToGeometry(GModel *gm)
{
-
- FILE *f = fopen("toto.pos","w");
- fprintf(f,"View \"\"{\n");
-
double Obj = 0.0;
- for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
- if ((*it)->geomType() == GEntity::Line)continue;
- for (unsigned int i=0;i<(*it)->lines.size(); i++){
- //Obj += MLineGEdgeDistance ( (*it)->lines[i] , *it,f );
- Obj = std::max(MLineGEdgeDistance ( (*it)->lines[i] , *it, f ),Obj);
- }
+ for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
+ if ((*it)->geomType() == GEntity::Line) continue;
+ for (unsigned int i=0;i<(*it)->lines.size(); i++)
+ Obj = std::max(MLineGEdgeDistance((*it)->lines[i], *it), Obj);
}
+ printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
- printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
-
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
- if ((*it)->geomType() == GEntity::Plane)continue;
- // printf("FACE %d with %d triangles\n",(*it)->tag(),(*it)->triangles.size());
- for (unsigned int i=0;i<(*it)->triangles.size(); i++){
- //Obj += MFaceGFaceDistance ( (*it)->triangles[i] , *it );
- Obj = std::max(Obj,MFaceGFaceDistance ( (*it)->triangles[i] , *it ));
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
+ if ((*it)->geomType() == GEntity::Plane) continue;
+ for (unsigned int i=0;i<(*it)->triangles.size(); i++) {
+ Obj = std::max(Obj,MFaceGFaceDistance( (*it)->triangles[i] , *it ));
}
}
-
printf("DISTANCE TO GEOMETRY : 1D AND 2D PART %22.15E\n",Obj);
- fprintf(f,"};\n");
- fclose(f);
+
return Obj;
}
-
-
diff --git a/contrib/HighOrderMeshOptimizer/OptHomCADDist.h b/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
index 465335331580b57c2b921c46cdcb6f30d2689031..d91deb6f96d689f90426c076090aa2e2fe962092 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
@@ -4,6 +4,7 @@
#define OPTHOMCADDIST_H_
+class GradientBasis;
class GModel;
class MElement;
class MLine;
@@ -11,10 +12,17 @@ class MTriangle;
class GEdge;
class GFace;
-double MFaceGFaceDistance(MTriangle *t, GFace *gf,
- std::vector<std::vector<SVector3> > *gsfT=0,
- std::map<MVertex*,SVector3> *normalsToCAD=0);
-double MLineGEdgeDistance (MLine *l, GEdge *ge, FILE *f = 0);
+double MFaceGFaceDistanceOld(MTriangle *t, GFace *gf,
+ std::vector<std::vector<SVector3> > *gsfT=0,
+ std::map<MVertex*,SVector3> *normalsToCAD=0);
+double MLineGEdgeDistanceOld(MLine *l, GEdge *ge, FILE *f = 0);
+
+double distToCAD1D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &tanCAD, double edLength);
+double distToCAD2D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &normCAD);
+double MLineGEdgeDistance (MLine *l, GEdge *ge);
+double MFaceGFaceDistance(MElement *el, GFace *gf);
void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*, double> &distances);
double distanceToGeometry(GModel *gm);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
index aee3baed41b9dfe77e91acf13527713bd8e5ca1f..7c31c4178e78af8372780b2da3bc1c67167d8fe0 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
@@ -10,7 +10,7 @@ template<class FuncType>
class ObjContribCADDist : public ObjContrib, public FuncType
{
public:
- ObjContribCADDist(double weight, double geomTol);
+ ObjContribCADDist(double weight, double refDist);
virtual ~ObjContribCADDist() {}
virtual ObjContrib *copy() const;
virtual void initialize(Patch *mesh);
@@ -24,14 +24,14 @@ public:
protected:
Patch *_mesh;
double _weight;
- double _geomTol;
+ double _refDist;
};
template<class FuncType>
-ObjContribCADDist<FuncType>::ObjContribCADDist(double weight, double geomTol) :
+ObjContribCADDist<FuncType>::ObjContribCADDist(double weight, double refDist) :
ObjContrib("CADDist", FuncType::getNamePrefix()+"CADDist"),
- _mesh(0), _weight(weight), _geomTol(geomTol)
+ _mesh(0), _weight(weight), _refDist(refDist)
{
}
@@ -47,7 +47,7 @@ template<class FuncType>
void ObjContribCADDist<FuncType>::initialize(Patch *mesh)
{
_mesh = mesh;
-
+ _mesh->initScaledCADDist(_refDist);
updateMinMax();
FuncType::initialize(_min, _max);
}
@@ -59,16 +59,26 @@ bool ObjContribCADDist<FuncType>::addContrib(double &Obj, alglib::real_1d_array
_min = BIGVAL;
_max = -BIGVAL;
- std::vector<double> gradF;
- for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ const int bndDim = _mesh->dim()-1;
+
+ for (int iBndEl = 0; iBndEl < _mesh->nBndEl(); iBndEl++) {
+ const int nVEl = _mesh->nNodBndEl(iBndEl);
double f;
- if (_mesh->bndDistAndGradients(iEl, f, gradF, _geomTol)) {
- _min = std::min(_min, f);
- _max = std::max(_max, f);
- Obj += FuncType::compute(f) * _weight;
- const double dFact = _weight * FuncType::computeDiff(f);
- for (size_t i = 0; i < _mesh->nPCEl(iEl); ++i)
- gradObj[_mesh->indPCEl(iEl, i)] += gradF[i] * dFact;
+ std::vector<double> gradF(nVEl*bndDim);
+ _mesh->scaledCADDistAndGradients(iBndEl, f, gradF);
+ _min = std::min(_min, f);
+ _max = std::max(_max, f);
+ Obj += FuncType::compute(f) * _weight;
+ const double dFact = _weight * FuncType::computeDiff(f);
+ for (int i=0; i<nVEl; i++) {
+ const int iFVi = _mesh->bndEl2FV(iBndEl, i);
+ if (iFVi >= 0) { // Skip if not free vertex
+ if (bndDim == 1) gradObj[_mesh->indPCFV(iFVi, 0)] += gradF[i] * dFact; // 2D
+ else { // 3D
+ gradObj[_mesh->indPCFV(iFVi, 0)] += gradF[2*i] * dFact;
+ gradObj[_mesh->indPCFV(iFVi, 1)] += gradF[2*i+1] * dFact;
+ }
+ }
}
}
@@ -82,13 +92,15 @@ void ObjContribCADDist<FuncType>::updateMinMax()
_min = BIGVAL;
_max = -BIGVAL;
- std::vector<double> dumGradF;
- for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ const int bndDim = _mesh->dim()-1;
+
+ for (int iBndEl = 0; iBndEl < _mesh->nBndEl(); iBndEl++) {
+ const int nVEl = _mesh->nNodBndEl(iBndEl);
double f;
- if (_mesh->bndDistAndGradients(iEl, f, dumGradF, _geomTol)) {
- _min = std::min(_min, f);
- _max = std::max(_max, f);
- }
+ std::vector<double> dumGradF(nVEl*bndDim);
+ _mesh->scaledCADDistAndGradients(iBndEl, f, dumGradF);
+ _min = std::min(_min, f);
+ _max = std::max(_max, f);
}
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h
index 944cac46bbb253c870d49a365e58fcacbb9f4e47..8031ba314c51d0a50c622cc07b47608d381271aa 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h
@@ -30,7 +30,7 @@ protected:
template<class FuncType>
ObjContribCADDistOld<FuncType>::ObjContribCADDistOld(double weight, double geomTol) :
- ObjContrib("CADDist", FuncType::getNamePrefix()+"CADDist"),
+ ObjContrib("CADDistOld", FuncType::getNamePrefix()+"CADDistOld"),
_mesh(0), _weight(weight), _geomTol(geomTol)
{
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 69e612606ea0cd6944008cb1aedc45f801bbba7f..0cd51f319212a06a09d96bc917c769f3ad2b992c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -691,6 +691,7 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
#include "GEntity.h"
//#include "MElement.h"
//#include "OptHomRun.h"
+#include "OptHomCADDist.h"
#include "MeshOptCommon.h"
#include "MeshOptObjContribFunc.h"
#include "MeshOptObjContrib.h"
@@ -706,6 +707,7 @@ struct HOPatchDefParameters : public MeshOptPatchDef
HOPatchDefParameters(const OptHomParameters &p);
virtual ~HOPatchDefParameters() {}
virtual double elBadness(MElement *el, GEntity* gEnt) const;
+ virtual double bndElBadness(MElement *el, GEntity* gEnt) const;
virtual double maxDistance(MElement *el) const;
virtual int inPatch(const SPoint3 &badBary,
double limDist, MElement *el,
@@ -714,7 +716,7 @@ private:
double jacMin, jacMax;
double distanceFactor;
bool optCAD;
- double optCADDistMax, optCADWeight, discrTolerance;
+ double optCADDistMax, optCADWeight;
};
@@ -737,7 +739,6 @@ HOPatchDefParameters::HOPatchDefParameters(const OptHomParameters &p)
optCAD = p.optCAD;
optCADDistMax = p.optCADDistMax;
optCADWeight = p.optCADWeight;
- discrTolerance = p.discrTolerance;
}
@@ -746,11 +747,23 @@ double HOPatchDefParameters::elBadness(MElement *el, GEntity* gEnt) const
double jmin, jmax;
el->scaledJacRange(jmin, jmax);
double badness = std::min(jmin-jacMin, 0.) + std::min(jacMax-jmax, 0.);
+ return badness;
+}
+
+
+double HOPatchDefParameters::bndElBadness(MElement *el, GEntity* gEnt) const
+{
if (optCAD) {
- const double dist = computeBndDist(el, 2, fabs(discrTolerance));
- badness += optCADWeight*std::min(optCADDistMax-dist, 0.);
+ if (el->getType() == TYPE_LIN) { // 2D
+ if (gEnt->geomType() != GEntity::Line) // Straight geometric line -> no distance
+ return MLineGEdgeDistance(static_cast<MLine*>(el), gEnt->cast2Edge());
+ }
+ else { // 3D
+ if (gEnt->geomType() != GEntity::Plane) // Straight geometric plance -> no distance
+ return MFaceGFaceDistance(el, gEnt->cast2Face());
+ }
}
- return badness;
+ return 1.;
}
@@ -778,6 +791,7 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
par.fixBndNodes = p.fixBndNodes;
par.useGeomForPatches = false;
par.useGeomForOpt = false;
+ par.useBoundaries = p.optCAD;
HOPatchDefParameters patchDef(p);
par.patchDef = &patchDef;
par.displayInterv = 30;
@@ -789,8 +803,8 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
minJacBarFunc.setTarget(p.BARRIER_MIN, 1.);
ObjContribScaledJac<ObjContribFuncBarrierFixMinMovMax> minMaxJacBarFunc(1.);
minMaxJacBarFunc.setTarget(p.BARRIER_MAX, 1.);
- ObjContribCADDist<ObjContribFuncSimpleTargetMax> CADDistFunc(p.optCADWeight, p.discrTolerance);
- CADDistFunc.setTarget(p.optCADDistMax);
+ ObjContribCADDist<ObjContribFuncSimpleTargetMax> CADDistFunc(p.optCADWeight, p.optCADDistMax);
+ CADDistFunc.setTarget(1.);
MeshOptPass minJacPass;
minJacPass.maxParamUpdates = p.optPassMax;
diff --git a/contrib/MeshOptimizer/MeshOpt.cpp b/contrib/MeshOptimizer/MeshOpt.cpp
index 261a8611f64072142fa116fad56c34090792f753..a35c3125ee64da15a598f9475b07825374068324 100644
--- a/contrib/MeshOptimizer/MeshOpt.cpp
+++ b/contrib/MeshOptimizer/MeshOpt.cpp
@@ -65,9 +65,10 @@ void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *MOInst)
MeshOpt::MeshOpt(const std::map<MElement*,GEntity*> &element2entity,
+ const std::map<MElement*, GEntity*> &bndEl2Ent,
const std::set<MElement*> &els, std::set<MVertex*> &toFix,
- const MeshOptParameters &par) :
- patch(element2entity, els, toFix, par.fixBndNodes), _verbose(0),
+ const std::set<MElement*> &bndEls, const MeshOptParameters &par) :
+ patch(element2entity, bndEl2Ent, els, toFix, bndEls, par.fixBndNodes), _verbose(0),
_iPass(0), _iter(0), _intervDisplay(0), _initObj(0)
{
_allObjFunc.resize(par.pass.size());
diff --git a/contrib/MeshOptimizer/MeshOpt.h b/contrib/MeshOptimizer/MeshOpt.h
index 844d0ca61e54b4e43438277f9a6dde991a6bf380..9c7a21c8085f34b9cf5cb4631d9352ed35305422 100644
--- a/contrib/MeshOptimizer/MeshOpt.h
+++ b/contrib/MeshOptimizer/MeshOpt.h
@@ -48,8 +48,9 @@ class MeshOpt
public:
Patch patch;
MeshOpt(const std::map<MElement*, GEntity*> &element2entity,
+ const std::map<MElement*, GEntity*> &bndEl2Ent,
const std::set<MElement*> &els, std::set<MVertex*> &toFix,
- const MeshOptParameters &par);
+ const std::set<MElement*> &bndEls, const MeshOptParameters &par);
~MeshOpt();
int optimize(const MeshOptParameters &par);
void updateMesh(const alglib::real_1d_array &x);
diff --git a/contrib/MeshOptimizer/MeshOptCommon.h b/contrib/MeshOptimizer/MeshOptCommon.h
index 151d115d1cf00264932ffc9bb0d99488df90476b..eddbcc042fbed34c85bd0f9b6f03a492b22079ec 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.h
+++ b/contrib/MeshOptimizer/MeshOptCommon.h
@@ -55,6 +55,8 @@ public:
virtual ~MeshOptPatchDef() {}
virtual double elBadness(MElement *el, // Determine "badness" of a given element (for patch creation)
GEntity* gEnt) const = 0;
+ virtual double bndElBadness(MElement *el, // Determine "badness" of a given boundary element (for patch creation)
+ GEntity* gEnt) const = 0;
virtual double maxDistance(MElement *el) const = 0; // Compute max. distance to a given bad element for elements in patch
virtual int inPatch(const SPoint3 &badBary, // Determine whether a given element should be included in the patch around a...
double limDist, // ... given bad element barycenter, with a limit distance if needed. Output: ...
@@ -78,6 +80,7 @@ struct MeshOptParameters { // Parameters controllin
bool onlyVisible ; // Apply optimization to visible entities ONLY
bool fixBndNodes; // If points can move on boundaries
bool useGeomForPatches, useGeomForOpt; // Whether to use info from CAD for creation of patches and for optimization
+ bool useBoundaries; // Whether to use boundary elements
MeshOptPatchDef *patchDef;
std::vector<MeshOptPass> pass;
int displayInterv; // Sampling rate in opt. iterations for display
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index b6d457cd91b1fd49dbf6c91f240b0febd2f45495..0f657ddf023bcc13e72badf983a200173d1d9d4e 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -36,13 +36,15 @@
#include "BasisFactory.h"
#include "CondNumBasis.h"
#include "OptHomIntegralBoundaryDist.h"
+#include "OptHomCADDist.h"
#include "qualityMeasures.h"
#include "MeshOptPatch.h"
Patch::Patch(const std::map<MElement*, GEntity*> &element2entity,
- const std::set<MElement*> &els, std::set<MVertex*> &toFix,
- bool fixBndNodes) :
+ const std::map<MElement*, GEntity*> &bndEl2Ent,
+ const std::set<MElement*> &els, std::set<MVertex*> &toFix,
+ const std::set<MElement*> &bndEls, bool fixBndNodes) :
_typeLengthScale(LS_NONE), _invLengthScaleSq(0.)
{
@@ -95,6 +97,39 @@ Patch::Patch(const std::map<MElement*, GEntity*> &element2entity,
if (nonGeoMove) Msg::Warning("Some vertices will be moved along local lines "
"or planes, they may not remain on the exact geometry");
+ // Initialize boundary elements and related connectivity if required
+ if (!bndEls.empty()) {
+ int nBndElts = bndEls.size();
+ _bndEl.resize(nBndElts);
+ _bndEl2Ent.resize(nBndElts);
+ _bndEl2V.resize(nBndElts);
+ _bndEl2FV.resize(nBndElts);
+ int iBndEl = 0;
+ bool unknownVert = false;
+ for(std::set<MElement*>::iterator it = bndEls.begin();
+ it != bndEls.end(); ++it, ++iBndEl) {
+ MElement* bndEl = *it;
+ _bndEl[iBndEl] = bndEl;
+ std::map<MElement*,GEntity*>::const_iterator itBndEl2Ent = bndEl2Ent.find(bndEl);
+ _bndEl2Ent[iBndEl] = (itBndEl2Ent == bndEl2Ent.end()) ? 0 : itBndEl2Ent->second;
+ int nBndElVerts = bndEl->getNumVertices();
+ _bndEl2V[iBndEl].resize(nBndElVerts);
+ _bndEl2FV[iBndEl].resize(nBndElVerts);
+ for (int iVBndEl=0; iVBndEl<nBndElVerts; iVBndEl++) {
+ MVertex *vert = bndEl->getVertex(iVBndEl);
+ std::vector<MVertex*>::iterator itV = std::find(_vert.begin(), _vert.end(), vert);
+ if (itV == _vert.end()) unknownVert = true;
+ else _bndEl2V[iBndEl][iVBndEl] = std::distance(_vert.begin(), itV);
+ std::vector<MVertex*>::iterator itFV = std::find(_freeVert.begin(),
+ _freeVert.end(), vert);
+ if (itFV == _freeVert.end()) _bndEl2FV[iBndEl][iVBndEl] = -1;
+ else _bndEl2FV[iBndEl][iVBndEl] = std::distance(_freeVert.begin(), itFV);
+ }
+ }
+ if (unknownVert) Msg::Error("Unknown vertices in boundary element "
+ "at patch initialization");
+ }
+
// Initial coordinates
_ixyz.resize(nVert());
for (int iV = 0; iV < nVert(); iV++) _ixyz[iV] = _vert[iV]->point();
@@ -309,6 +344,12 @@ void Patch::initMetricMin()
}
+void Patch::initScaledCADDist(double refCADDist)
+{
+ _invRefCADDist = 1./refCADDist;
+}
+
+
void Patch::calcNormalEl2D(int iEl, NormalScaling scaling,
fullMatrix<double> &elNorm, bool ideal)
{
@@ -524,6 +565,107 @@ bool Patch::bndDistAndGradients(int iEl, double &f, std::vector<double> &gradF,
}
+void Patch::scaledCADDistAndGradients(int iBndEl, double &scaledDist,
+ std::vector<double> &gradScaledDist)
+{
+ const std::vector<int> &iV = _bndEl2V[iBndEl], &iFV = _bndEl2FV[iBndEl];
+ const int nV = iV.size();
+ const GradientBasis *gb = BasisFactory::getGradientBasis(FuncSpaceData(_bndEl[iBndEl]));
+
+ // Coordinates of nodes
+ fullMatrix<double> nodesXYZ(nV, 3);
+ for (int i = 0; i < nV; i++) {
+ const int &iVi = iV[i];
+ nodesXYZ(i,0) = _xyz[iVi].x();
+ nodesXYZ(i,1) = _xyz[iVi].y();
+ nodesXYZ(i,2) = _xyz[iVi].z();
+ }
+
+ // Compute distance and gradients (CAUTION: returns gradients w.r.t. vertices, not free vertices)
+ if (_dim == 2) { // 2D
+ const GEdge *ge = _bndEl2Ent[iBndEl]->cast2Edge();
+ const Range<double> parBounds = ge->parBounds(0);
+ const double eps = 1.e-6 * (parBounds.high()-parBounds.low());
+ std::vector<SVector3> tanCAD(nV);
+ for (int i=0; i<nV; i++) {
+ const int &iVi = iV[i], &iFVi = iFV[i];
+ MVertex* &vert = _vert[iVi];
+ double tCAD;
+ if (iFVi >= 0) // If free vertex, ...
+ tCAD = _uvw[iFVi].x(); // ... get stored param. coord. (can be only line).
+ else { // Otherwise, get param. coord. from CAD.
+ if(ge->getBeginVertex() && ge->getBeginVertex()->mesh_vertices[0] == vert)
+ tCAD = parBounds.low();
+ else if(ge->getEndVertex() && ge->getEndVertex()->mesh_vertices[0] == vert)
+ tCAD = parBounds.high();
+ else
+ tCAD = ge->parFromPoint(_xyz[iVi]);
+ }
+ tanCAD[i] = ge->firstDer(tCAD); // Compute tangent at vertex
+ tanCAD[i].normalize(); // Normalize tangent
+ }
+ const double edLength = _xyz[iV[0]].distance(_xyz[iV[1]]); // Get edge length
+ scaledDist = _invRefCADDist * distToCAD1D(gb, nodesXYZ, tanCAD, edLength);
+ for (int i=0; i<nV; i++) {
+ const int &iFVi = iFV[i];
+ if (iFVi < 0) continue; // Skip if not free vertex
+ const double xS = nodesXYZ(i, 0), yS = nodesXYZ(i, 1), zS = nodesXYZ(i, 2); // Save coord. of perturbed node for FD
+ const SVector3 tanCADS = tanCAD[i]; // Save tangent to CAD at perturbed node
+ const double tCAD = _uvw[iFVi].x() + eps; // New param. coord. of perturbed node
+ GPoint gp = ge->point(tCAD); // New coord. of perturbed node
+ nodesXYZ(i, 0) = gp.x(); nodesXYZ(i, 1) = gp.y(); nodesXYZ(i, 2) = gp.z();
+ tanCAD[i] = ge->firstDer(tCAD); // New tangent to CAD at perturbed node
+ tanCAD[i].normalize(); // Normalize new tangent
+ double sDistDiff = _invRefCADDist * distToCAD1D(gb, nodesXYZ, tanCAD, edLength); // Compute distance with perturbed node
+ gradScaledDist[i] = (sDistDiff-scaledDist) / eps; // Compute gradient
+ nodesXYZ(i, 0) = xS; nodesXYZ(i, 1) = yS; nodesXYZ(i, 2) = zS; // Restore coord. of perturbed node
+ tanCAD[i] = tanCADS; // Restore tan. to CAD at perturbed node
+ }
+ }
+ else { // 3D
+ const GFace *gf = _bndEl2Ent[iBndEl]->cast2Face();
+ const Range<double> parBounds0 = gf->parBounds(0), parBounds1 = gf->parBounds(1);
+ const double eps0 = 1.e-6 * (parBounds0.high()-parBounds0.low());
+ const double eps1 = 1.e-6 * (parBounds1.high()-parBounds1.low());
+ std::vector<SVector3> normCAD(nV);
+ for (int i=0; i<nV; i++) {
+ const int &iVi = iV[i], &iFVi = iFV[i];
+ MVertex* &vert = _vert[iVi];
+ SPoint2 pCAD;
+ if ((iFVi >= 0) && (vert->onWhat() == gf)) // If free vertex and on surface, ...
+ pCAD = SPoint2(_uvw[iFVi].x(), _uvw[iFVi].y()); // ... get stored param. coord.
+ else
+ reparamMeshVertexOnFace(vert, gf, pCAD); // Otherwise, get param. coord. from CAD.
+ normCAD[i] = gf->normal(pCAD); // Compute normal at vertex
+ normCAD[i].normalize(); // Normalize normal
+ }
+ scaledDist = _invRefCADDist * distToCAD2D(gb, nodesXYZ, normCAD);
+ for (int i=0; i<nV; i++) {
+ const int &iFVi = iFV[i];
+ if (iFVi < 0) continue; // Skip if not free vertex
+ const double xS = nodesXYZ(i, 0), yS = nodesXYZ(i, 1), zS = nodesXYZ(i, 2); // Save coord. of perturbed node for FD
+ const SVector3 tanCADS = normCAD[i]; // Save normal to CAD at perturbed node
+ const SPoint2 pCAD0 = SPoint2(_uvw[iFVi].x()+eps0, _uvw[iFVi].y()); // New param. coord. of perturbed node in 1st dir.
+ GPoint gp0 = gf->point(pCAD0); // New coord. of perturbed node in 1st dir.
+ nodesXYZ(i, 0) = gp0.x(); nodesXYZ(i, 1) = gp0.y(); nodesXYZ(i, 2) = gp0.z();
+ normCAD[i] = gf->normal(pCAD0); // New normal to CAD at perturbed node in 1st dir.
+ normCAD[i].normalize(); // Normalize new normal
+ double sDistDiff0 = _invRefCADDist * distToCAD2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 1st dir.
+ gradScaledDist[2*i] = (sDistDiff0-scaledDist) / eps0; // Compute gradient in 1st dir.
+ const SPoint2 pCAD1 = SPoint2(_uvw[iFVi].x(), _uvw[iFVi].y()+eps1); // New param. coord. of perturbed node in 2nd dir.
+ GPoint gp1 = gf->point(pCAD1); // New coord. of perturbed node in 2nd dir.
+ nodesXYZ(i, 0) = gp1.x(); nodesXYZ(i, 1) = gp1.y(); nodesXYZ(i, 2) = gp1.z();
+ normCAD[i] = gf->normal(pCAD1); // New normal to CAD at perturbed node in 2nd dir.
+ normCAD[i].normalize(); // Normalize new normal
+ double sDistDiff1 = _invRefCADDist * distToCAD2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 2nd dir.
+ gradScaledDist[2*i+1] = (sDistDiff1-scaledDist) / eps0; // Compute gradient in 2nd dir.
+ nodesXYZ(i, 0) = xS; nodesXYZ(i, 1) = yS; nodesXYZ(i, 2) = zS; // Restore coord. of perturbed node
+ normCAD[i] = tanCADS; // Restore tan. to CAD at perturbed node
+ }
+ }
+}
+
+
void Patch::initIdealJac()
{
// Initialize _nBezEl
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index 6b1c0d16b0aead1816d8468306dfda0086e1cbb0..e0d7fb1bc98647a8be321705d12d204a2565ab99 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/contrib/MeshOptimizer/MeshOptPatch.h
@@ -47,7 +47,9 @@ class Patch
{
public:
Patch(const std::map<MElement*,GEntity*> &element2entity,
- const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBndNodes);
+ const std::map<MElement*, GEntity*> &bndEl2Ent,
+ const std::set<MElement*> &els, std::set<MVertex*> &toFix,
+ const std::set<MElement*> &bndEls, bool fixBndNodes);
// Mesh entities and variables
inline const int &dim() { return _dim; }
@@ -85,11 +87,17 @@ public:
// High-order: scaled Jacobian and metric measures, distance to CAD
inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; }
+ inline int nBndEl() { return _bndEl.size(); }
+ inline int nNodBndEl(int iBndEl) { return _bndEl2V[iBndEl].size(); }
+ inline const int &bndEl2FV(int iBndEl, int i) { return _bndEl2FV[iBndEl][i]; }
void initScaledJac();
void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
void initMetricMin();
void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
- bool bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps);
+ bool bndDistAndGradients(int iEl, double &f, std::vector<double> &gradF, double eps);
+ void initScaledCADDist(double refCADDist);
+ void scaledCADDistAndGradients(int iBndEl, double &scaledDist,
+ std::vector<double> &gradScaledDist);
// Mesh quality
inline const int &nIJacEl(int iEl) { return _nIJacEl[iEl]; }
@@ -147,6 +155,10 @@ private:
std::vector<int> _nBezEl; // Number of Bezier poly. for an el.
std::vector<fullMatrix<double> > _JacNormEl; // Normals to 2D elements for Jacobian regularization and scaling
std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements
+ std::vector<MElement*> _bndEl; // Boundary elements
+ std::vector<std::vector<int> > _bndEl2V, _bndEl2FV; // Vertices & corresponding free vertices on the boundary elements
+ std::vector<GEntity*> _bndEl2Ent; // Geometric entities corresponding to the boundary elements
+ double _invRefCADDist;
void calcNormalEl2D(int iEl, NormalScaling scaling,
fullMatrix<double> &elNorm, bool ideal);
diff --git a/contrib/MeshOptimizer/MeshOptimizer.cpp b/contrib/MeshOptimizer/MeshOptimizer.cpp
index f52165efd7323afd704a5971b548251daaf538ef..3365908d13e72f816483ce3abf214adc2ab6b0ae 100644
--- a/contrib/MeshOptimizer/MeshOptimizer.cpp
+++ b/contrib/MeshOptimizer/MeshOptimizer.cpp
@@ -55,6 +55,7 @@ typedef std::set<MElement*> elSet;
typedef elSet::iterator elSetIter;
typedef std::set<MVertex*> vertSet;
typedef std::map<MElement*, GEntity*> elEntMap;
+typedef std::map<MElement*, MElement*> elElMap;
typedef std::map<MVertex*, elVec> vertElVecMap;
typedef std::map<MElement*, elSet> elElSetMap;
@@ -178,6 +179,96 @@ void calcElement2Entity(GEntity *entity, elEntMap &element2entity)
}
+MElement *getFaceInBndElements(const MFace &f, std::list<GFace*> gFaces)
+{
+ for (std::list<GFace*>::iterator itGF = gFaces.begin(); itGF != gFaces.end(); itGF++) {
+ if (f.getNumVertices() == 3) {
+ std::vector<MTriangle*> &tris = (*itGF)->triangles;
+ for (int iEl=0; iEl<tris.size(); iEl++)
+ if (tris[iEl]->getFace(0) == f) return tris[iEl];
+ }
+ else {
+ std::vector<MQuadrangle*> &quads = (*itGF)->quadrangles;
+ for (int iEl=0; iEl<quads.size(); iEl++)
+ if (quads[iEl]->getFace(0) == f) return quads[iEl];
+ }
+ }
+ return 0;
+}
+
+
+MElement *getEdgeInBndElements(const MEdge &e, std::list<GEdge*> gEdges)
+{
+ for (std::list<GEdge*>::iterator itGE = gEdges.begin(); itGE != gEdges.end(); itGE++) {
+ std::vector<MLine*> &lines = (*itGE)->lines;
+ for (int iEl=0; iEl<lines.size(); iEl++)
+ if (lines[iEl]->getEdge(0) == e) return lines[iEl];
+ }
+ return 0;
+}
+
+
+void calcBndInfo(GEntity *entity, elElMap &el2BndEl, elEntMap &bndEl2Ent)
+{
+ typedef std::list<GFace*> GFaceList;
+ typedef std::list<GEdge*> GEdgeList;
+
+ if (entity->dim() == 3) { // 3D
+
+ // Fill boundary element -> GEntity connectivity
+ GFaceList gFaces = entity->faces();
+ for (GFaceList::iterator itGF = gFaces.begin(); itGF != gFaces.end(); itGF++) {
+ std::vector<MTriangle*> &tris = (*itGF)->triangles;
+ for (int i=0; i<tris.size(); i++)
+ bndEl2Ent.insert(std::pair<MElement*, GEntity*>(tris[i], *itGF));
+ std::vector<MQuadrangle*> &quads = (*itGF)->quadrangles;
+ for (int i=0; i<quads.size(); i++)
+ bndEl2Ent.insert(std::pair<MElement*, GEntity*>(quads[i], *itGF));
+ }
+
+ // Fill element -> boundary element connectivity
+ for (int iEl = 0; iEl < entity->getNumMeshElements(); iEl++) {
+ MElement *el = entity->getMeshElement(iEl);
+ int nBndVert = 0; // Compute nb. of bnd. vertices in element
+ for (int iV=0; iV<el->getNumPrimaryVertices(); iV++)
+ if (el->getVertex(iV)->onWhat() != entity) nBndVert++;
+ if (nBndVert >= 3) // If more than 3 primary vert. on bnd., look for bnd. face(s)
+ for (int iF=0; iF<el->getNumFaces(); iF++) {
+ MElement *bndEl = getFaceInBndElements(el->getFace(iF), gFaces);
+ if (bndEl != 0)
+ el2BndEl.insert(std::pair<MElement*, MElement*>(el, bndEl));
+ }
+ }
+
+ }
+ else { // 2D
+
+ // Fill boundary element -> GEntity connectivity
+ GEdgeList gEdges = entity->edges();
+ for (GEdgeList::iterator itGE = gEdges.begin(); itGE != gEdges.end(); itGE++) {
+ std::vector<MLine*> &lines = (*itGE)->lines;
+ for (int i=0; i<lines.size(); i++)
+ bndEl2Ent.insert(std::pair<MElement*, GEntity*>(lines[i], *itGE));
+ }
+
+ // Fill element -> boundary element connectivity
+ for (int iEl = 0; iEl < entity->getNumMeshElements(); iEl++) {
+ MElement *el = entity->getMeshElement(iEl);
+ int nBndVert = 0; // Compute nb. of bnd. vertices in element
+ for (int iV=0; iV<el->getNumPrimaryVertices(); iV++)
+ if (el->getVertex(iV)->onWhat() != entity) nBndVert++;
+ if (nBndVert >= 2) // If more than 2 primary vert. on bnd., look for bnd. edge(s)
+ for (int iE=0; iE<el->getNumEdges(); iE++) {
+ MElement *bndEl = getEdgeInBndElements(el->getEdge(iE), gEdges);
+ if (bndEl != 0)
+ el2BndEl.insert(std::pair<MElement*, MElement*>(el, bndEl));
+ }
+ }
+
+ }
+}
+
+
std::vector<elSetVertSetPair> getDisjointPatches(const vertElVecMap &vertex2elements,
const elEntMap &element2entity,
const elSet &badElements,
@@ -253,8 +344,25 @@ std::vector<elSetVertSetPair> getDisjointPatches(const vertElVecMap &vertex2elem
}
+// Get (bad) boundary elements adjacent to patch
+void getAdjacentBndElts(const elElMap &el2BndEl, const elEntMap &bndEl2Ent,
+ const elSet &elts, elSet &bndElts, MeshOptParameters &par)
+{
+ for (elSetIter itEl=elts.begin(); itEl!=elts.end(); itEl++) {
+ elElMap::const_iterator itBndEl = el2BndEl.find(*itEl);
+ if (itBndEl != el2BndEl.end()) {
+ MElement* bndEl = itBndEl->second;
+ elEntMap::const_iterator itEnt = bndEl2Ent.find(bndEl);
+ if (par.patchDef->bndElBadness(bndEl, itEnt->second) < 0.) bndElts.insert(bndEl);
+ }
+ }
+}
+
+
void optimizeDisjointPatches(const vertElVecMap &vertex2elements,
const elEntMap &element2entity,
+ const elElMap &el2BndEl,
+ const elEntMap &bndEl2Ent,
elSet &badasses, MeshOptParameters &par)
{
par.success = 1;
@@ -262,16 +370,27 @@ void optimizeDisjointPatches(const vertElVecMap &vertex2elements,
const elEntMap &e2ePatch = par.useGeomForPatches ? element2entity : elEntMap();
const elEntMap &e2eOpt = par.useGeomForOpt ? element2entity : elEntMap();
+ // Get patches
std::vector<elSetVertSetPair> toOptimize = getDisjointPatches(vertex2elements,
e2ePatch, badasses, par);
+ // Get boundary elements adjacent to patch if required
+ std::vector<elSet> bndElts;
+ bndElts.resize(toOptimize.size());
+ if (!el2BndEl.empty()) {
+ for (int iPatch = 0; iPatch < toOptimize.size(); ++iPatch)
+ getAdjacentBndElts(el2BndEl, bndEl2Ent, toOptimize[iPatch].first, bndElts[iPatch], par);
+ }
+
for (int iPatch = 0; iPatch < toOptimize.size(); ++iPatch) {
// Initialize optimization and output if asked
if (par.verbose > 1)
- Msg::Info("Optimizing patch %i/%i composed of %i elements", iPatch,
- toOptimize.size()-1, toOptimize[iPatch].first.size());
- MeshOpt opt(e2eOpt, toOptimize[iPatch].first, toOptimize[iPatch].second, par);
+ Msg::Info("Optimizing patch %i/%i composed of %i elements, "
+ "%i boundary elements", iPatch, toOptimize.size()-1,
+ toOptimize[iPatch].first.size(), bndElts[iPatch].size());
+ MeshOpt opt(e2eOpt, bndEl2Ent, toOptimize[iPatch].first,
+ toOptimize[iPatch].second, bndElts[iPatch], par);
if (par.verbose > 3) {
std::ostringstream ossI1;
ossI1 << "initial_patch-" << iPatch << ".msh";
@@ -327,6 +446,8 @@ MElement *getWorstElement(elSet &badElts,
void optimizeOneByOne(const vertElVecMap &vertex2elements,
const elEntMap &element2entity,
+ const elElMap &el2BndEl,
+ const elEntMap &bndEl2Ent,
elSet badElts, MeshOptParameters &par)
{
par.success = 1;
@@ -367,11 +488,17 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
badElts.begin(),badElts.end(),
std::inserter(toOptimize, toOptimize.end()));
+ // Get boundary elements adjacent to patch if required
+ elSet bndElts;
+ if (!el2BndEl.empty()) {
+ getAdjacentBndElts(el2BndEl, bndEl2Ent, toOptimize, bndElts, par);
+ }
+
// Initialize optimization and output if asked
if (par.verbose > 1)
Msg::Info("Optimizing patch %i (max. %i remaining) composed of %4d elements",
iBadEl, badElts.size(), toOptimize.size());
- MeshOpt opt(e2eOpt, toOptimize, toFix, par);
+ MeshOpt opt(e2eOpt, bndEl2Ent, toOptimize, toFix, bndElts, par);
if (par.verbose > 3) {
std::ostringstream ossI1;
ossI1 << "initial_patch-" << iBadEl << ".msh";
@@ -442,7 +569,8 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par)
gm->getEntities(entities);
vertElVecMap vertex2elements;
- elEntMap element2entity;
+ elEntMap element2entity, bndEl2Ent;
+ elElMap el2BndEl;
elSet badElts;
for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
GEntity* &entity = entities[iEnt];
@@ -453,18 +581,28 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par)
calcVertex2Elements(par.dim, entity, vertex2elements);
if ((par.useGeomForPatches) || (par.useGeomForOpt))
calcElement2Entity(entity, element2entity);
- for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
- double jmin, jmax;
+ if (par.useBoundaries) calcBndInfo(entity, el2BndEl, bndEl2Ent);
+ for (int iEl = 0; iEl < entity->getNumMeshElements(); iEl++) { // Detect bad elements
MElement *el = entity->getMeshElement(iEl);
- if ((el->getDim() == par.dim) && (par.patchDef->elBadness(el, entity) < 0.))
- badElts.insert(el);
+ if (el->getDim() == par.dim) {
+ if (par.patchDef->elBadness(el, entity) < 0.) badElts.insert(el);
+ else if (par.useBoundaries) {
+ elElMap::iterator bndElIt = el2BndEl.find(el);
+ if (bndElIt != el2BndEl.end()) {
+ MElement* &bndEl = bndElIt->second;
+ if (par.patchDef->bndElBadness(bndEl, bndEl2Ent[bndEl])) badElts.insert(el);
+ }
+ }
+ }
}
}
if (par.patchDef->strategy == MeshOptPatchDef::STRAT_DISJOINT)
- optimizeDisjointPatches(vertex2elements, element2entity, badElts, par);
+ optimizeDisjointPatches(vertex2elements, element2entity,
+ el2BndEl, bndEl2Ent, badElts, par);
else if (par.patchDef->strategy == MeshOptPatchDef::STRAT_ONEBYONE)
- optimizeOneByOne(vertex2elements, element2entity, badElts, par);
+ optimizeOneByOne(vertex2elements, element2entity,
+ el2BndEl, bndEl2Ent, badElts, par);
else
Msg::Error("Unknown strategy %d for mesh optimization", par.patchDef->strategy);
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
index 0faff2097a02671a7b2fd1df8968f47fe20f94bb..faf6a8452b5e21bdb1c9f625b390fcbde2ef4935 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
@@ -23,6 +23,7 @@ struct QualPatchDefParameters : public MeshOptPatchDef
QualPatchDefParameters(const MeshQualOptParameters &p);
virtual ~QualPatchDefParameters() {}
virtual double elBadness(MElement *el, GEntity* gEnt) const;
+ virtual double bndElBadness(MElement *el, GEntity* gEnt) const { return 1.; }
virtual double maxDistance(MElement *el) const;
virtual int inPatch(const SPoint3 &badBary, double limDist,
MElement *el, GEntity* gEnt) const;
@@ -126,6 +127,7 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
par.fixBndNodes = p.fixBndNodes;
par.useGeomForPatches = p.excludeBL;
par.useGeomForOpt = false;
+ par.useBoundaries = false;
QualPatchDefParameters patchDef(p);
par.patchDef = &patchDef;
par.displayInterv = 20;