Skip to content
Snippets Groups Projects
Commit 62457616 authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Added "CADDist" contribution to objective function in high-order mesh optimizer

parent aeea6a67
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
+ 29
− 55
View file @ 62457616
......@@ -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_ */
contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+ 1
− 0
View file @ 62457616
......@@ -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"
......
contrib/MeshOptimizer/MeshOptCommon.h
+ 1
− 0
View file @ 62457616
......@@ -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;
......
contrib/MeshOptimizer/MeshOptPatch.cpp
+ 60
− 59
View file @ 62457616
......@@ -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;
}
contrib/MeshOptimizer/MeshOptPatch.h
+ 2
− 4
View file @ 62457616
......@@ -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:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment