From 624576160c0239a55a393af2c7d6428b0d565690 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Mon, 22 Sep 2014 08:45:01 +0000
Subject: [PATCH] Added "CADDist" contribution to objective function in
high-order mesh optimizer
---
.../OptHomObjContribCADDist.h | 84 +++++--------
contrib/HighOrderMeshOptimizer/OptHomRun.cpp | 1 +
contrib/MeshOptimizer/MeshOptCommon.h | 1 +
contrib/MeshOptimizer/MeshOptPatch.cpp | 119 +++++++++---------
contrib/MeshOptimizer/MeshOptPatch.h | 6 +-
5 files changed, 93 insertions(+), 118 deletions(-)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
index d4b8d1e2aa..40aeccd354 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
@@ -10,8 +10,9 @@ template<class FuncType>
class ObjContribCADDist : public ObjContrib, public FuncType
{
public:
- ObjContribCADDist(double weight);
+ ObjContribCADDist(double weight, double geomTol);
virtual ~ObjContribCADDist() {}
+ virtual ObjContrib *copy() const;
virtual void initialize(Patch *mesh);
virtual bool fail() { return false; }
virtual bool addContrib(double &Obj, alglib::real_1d_array &gradObj);
@@ -24,31 +25,30 @@ public:
protected:
Patch *_mesh;
double _weight;
- std::vector<int> _nBezEl; // Number of Bezier poly. and for an el.
- inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
- inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; }
- void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
+ double _geomTol;
};
template<class FuncType>
-ObjContribCADDist<FuncType>::ObjContribCADDist(double weight) :
- ObjContrib("MetricMin", FuncType::getNamePrefix()+"MetricMin"),
- _mesh(0), _weight(weight)
+ObjContribCADDist<FuncType>::ObjContribCADDist(double weight, double geomTol) :
+ ObjContrib("MetricMin", FuncType::getNamePrefix()+"CADDist"),
+ _mesh(0), _weight(weight), _geomTol(geomTol)
{
}
+template<class FuncType>
+ObjContrib *ObjContribCADDist<FuncType>::copy() const
+{
+ return new ObjContribCADDist<FuncType>(*this);
+}
+
+
template<class FuncType>
void ObjContribCADDist<FuncType>::initialize(Patch *mesh)
{
_mesh = mesh;
- // Initialize _nBezEl
- _nBezEl.resize(_mesh->nEl());
- for (int iEl=0; iEl<_mesh->nEl(); iEl++)
- _nBezEl[iEl] = _mesh->el(iEl)->getJacobianFuncSpace()->getNumJacNodes();
-
updateMinMax();
FuncType::initialize(_min, _max);
}
@@ -60,17 +60,16 @@ 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++) {
- std::vector<double> mM(_nBezEl[iEl]); // Min. of Metric
- std::vector<double> gMM(_nBezEl[iEl]*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
- metricMinAndGradients(iEl, mM, gMM);
- for (int l = 0; l < _nBezEl[iEl]; l++) { // Add contribution for each Bezier coeff.
- Obj += _weight * FuncType::compute(mM[l]);
- for (int iPC = 0; iPC < _mesh->nPCEl(iEl); iPC++)
- gradObj[_mesh->indPCEl(iEl,iPC)] +=
- _weight * FuncType::computeDiff(mM[l], gMM[indGSJ(iEl, l, iPC)]);
- _min = std::min(_min, mM[l]);
- _max = std::max(_max, mM[l]);
+ 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;
}
}
@@ -84,13 +83,12 @@ void ObjContribCADDist<FuncType>::updateMinMax()
_min = BIGVAL;
_max = -BIGVAL;
+ std::vector<double> dumGradF;
for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
- std::vector<double> mM(_nBezEl[iEl]); // Min. of Metric
- std::vector<double> dumGMM(_nBezEl[iEl]*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
- metricMinAndGradients(iEl, mM, dumGMM);
- for (int l = 0; l < _nBezEl[iEl]; l++) { // Check each Bezier coeff.
- _min = std::min(_min, mM[l]);
- _max = std::max(_max, mM[l]);
+ double f;
+ if (_mesh->bndDistAndGradients(iEl, f, dumGradF, _geomTol)) {
+ _min = std::min(_min, f);
+ _max = std::max(_max, f);
}
}
}
@@ -99,33 +97,9 @@ void ObjContribCADDist<FuncType>::updateMinMax()
template<class FuncType>
void ObjContribCADDist<FuncType>::updateResults(MeshOptResults &res) const
{
- res.minMetricMin = std::min(_min, res.minMetricMin);
- res.maxMetricMin = std::max(_max, res.maxMetricMin);
+ res.minCADDist = std::min(_min, res.minCADDist);
+ res.maxCADDist = std::max(_max, res.maxCADDist);
}
-//bool MeshOpt::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;
-//
-//}
-
-
#endif /* _OPTHOMOBJCONTRIBCADDIST_H_ */
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 52095e5862..cb47ed375b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -749,6 +749,7 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
#include "MeshOptObjContribScaledNodeDispSq.h"
#include "OptHomObjContribScaledJac.h"
#include "OptHomObjContribMetricMin.h"
+#include "OptHomObjContribCADDist.h"
#include "MeshOptimizer.h"
diff --git a/contrib/MeshOptimizer/MeshOptCommon.h b/contrib/MeshOptimizer/MeshOptCommon.h
index 7314c0ea5b..c2e2f2f366 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.h
+++ b/contrib/MeshOptimizer/MeshOptCommon.h
@@ -43,6 +43,7 @@ struct MeshOptResults { // Output of mesh optimizati
double minNodeDisp, maxNodeDisp; // Range of node displacement
double minScaledJac, maxScaledJac; // Range of Scaled Jacobians
double minMetricMin, maxMetricMin; // Range of min. of metric
+ double minCADDist, maxCADDist; // Range of distance to CAD
MeshOptResults();
private:
static const double BIGVAL;
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index d812dcf0a7..0577a25257 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -29,11 +29,12 @@
#include "GmshMessage.h"
#include "GRegion.h"
+#include "GFace.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "BasisFactory.h"
-//#include "OptHomIntegralBoundaryDist.h"
+#include "OptHomIntegralBoundaryDist.h"
#include "MeshOptPatch.h"
@@ -437,61 +438,61 @@ void Patch::metricMinAndGradients(int iEl, std::vector<double> &lambda,
}
-//bool Patch::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;
-//}
+bool Patch::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;
+}
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index 2207c839bb..209f173c0a 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/contrib/MeshOptimizer/MeshOptPatch.h
@@ -83,16 +83,14 @@ public:
double scaledNodeDispSq(int iFV);
void gradScaledNodeDispSq(int iFV, std::vector<double> &gDSq);
- // High-order: scaled Jacobian and metric measures
+ // 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; }
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);
-
- // TODO: Re-introduce distance to boundary
-// bool bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps);
+ bool bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps);
private:
--
GitLab