diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribMetricMin.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribMetricMin.h index e6d98681b64c31ed917b0499e83795bdd4665d30..3d44138e33fe4b2b16e95cbeb07a94ae2999ea5b 100644 --- a/contrib/HighOrderMeshOptimizer/OptHomObjContribMetricMin.h +++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribMetricMin.h @@ -25,10 +25,6 @@ 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); }; @@ -44,12 +40,7 @@ template<class FuncType> void ObjContribMetricMin<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(); - + _mesh->initScaledJac(); updateMinMax(); FuncType::initialize(_min, _max); } @@ -62,14 +53,14 @@ bool ObjContribMetricMin<FuncType>::addContrib(double &Obj, alglib::real_1d_arra _max = -BIGVAL; 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. + std::vector<double> mM(_mesh->nBezEl(iEl)); // Min. of Metric + std::vector<double> gMM(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min. + _mesh->metricMinAndGradients(iEl, mM, gMM); + for (int l = 0; l < _mesh->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)]); + _weight * FuncType::computeDiff(mM[l], gMM[_mesh->indGSJ(iEl, l, iPC)]); _min = std::min(_min, mM[l]); _max = std::max(_max, mM[l]); } @@ -86,10 +77,10 @@ void ObjContribMetricMin<FuncType>::updateMinMax() _max = -BIGVAL; 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. + std::vector<double> mM(_mesh->nBezEl(iEl)); // Min. of Metric + std::vector<double> dumGMM(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min. + _mesh->metricMinAndGradients(iEl, mM, dumGMM); + for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Check each Bezier coeff. _min = std::min(_min, mM[l]); _max = std::max(_max, mM[l]); } @@ -105,60 +96,4 @@ void ObjContribMetricMin<FuncType>::updateResults(MeshOptResults &res) const } -template<class FuncType> -void ObjContribMetricMin<FuncType>::metricMinAndGradients(int iEl, std::vector<double> &lambda, - std::vector<double> &gradLambda) -{ - const JacobianBasis *jacBasis = _mesh->el(iEl)->getJacobianFuncSpace(); - const int &numJacNodes = jacBasis->getNumJacNodes(); - const int &numMapNodes = jacBasis->getNumMapNodes(); - const int &numPrimMapNodes = jacBasis->getNumPrimMapNodes(); - fullVector<double> lambdaJ(numJacNodes), lambdaB(numJacNodes); - fullMatrix<double> gradLambdaJ(numJacNodes, 2*numMapNodes); - fullMatrix<double> gradLambdaB(numJacNodes, 2*numMapNodes); - - // Coordinates of nodes - fullMatrix<double> nodesXYZ(numMapNodes, 3), nodesXYZStraight(numPrimMapNodes, 3); - for (int i = 0; i < numMapNodes; i++) { - const int &iVi = _mesh->el2V(iEl, i); - nodesXYZ(i, 0) = _mesh->xyz(iVi).x(); - nodesXYZ(i, 1) = _mesh->xyz(iVi).y(); - nodesXYZ(i, 2) = _mesh->xyz(iVi).z(); - if (i < numPrimMapNodes) { - nodesXYZStraight(i, 0) = _mesh->ixyz(iVi).x(); - nodesXYZStraight(i, 1) = _mesh->ixyz(iVi).y(); - nodesXYZStraight(i, 2) = _mesh->ixyz(iVi).z(); - } - } - - jacBasis->getMetricMinAndGradients(nodesXYZ, nodesXYZStraight, lambdaJ, gradLambdaJ); - - lambdaB = lambdaJ; - gradLambdaB = gradLambdaJ; - - int iPC = 0; - std::vector<SPoint3> gXyzV(numJacNodes); - std::vector<SPoint3> gUvwV(numJacNodes); - for (int l = 0; l < numJacNodes; l++) { - lambda[l] = lambdaB(l); - } - for (int i = 0; i < numMapNodes; i++) { - const int &iFVi = _mesh->el2FV(iEl, i); - if (iFVi >= 0) { - for (int l = 0; l < numJacNodes; l++) { - gXyzV [l] = SPoint3(gradLambdaB(l, i), - gradLambdaB(l, i+numMapNodes),/*BDB(l,i+2*nbNod)*/ 0.); - } - _mesh->gXyz2gUvw(iFVi, gXyzV, gUvwV); - for (int l = 0; l < numJacNodes; l++) { - gradLambda[indGSJ(iEl, l, iPC)] = gUvwV[l][0]; - if (_mesh->nPCFV(iFVi) >= 2) gradLambda[indGSJ(iEl, l, iPC+1)] = gUvwV[l][1]; - if (_mesh->nPCFV(iFVi) == 3) gradLambda[indGSJ(iEl, l, iPC+2)] = gUvwV[l][2]; - } - iPC += _mesh->nPCFV(iFVi); - } - } -} - - #endif /* _OPTHOMOBJCONTRIBMETRICMIN_H_ */ diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribScaledJac.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribScaledJac.h index 6aa91d0b1810549911415ab3afa2c7d9a95d77db..dd260fc8c18796628a3244059f0e8a6000ecf5db 100644 --- a/contrib/HighOrderMeshOptimizer/OptHomObjContribScaledJac.h +++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribScaledJac.h @@ -25,14 +25,6 @@ public: protected: Patch *_mesh; double _weight; - std::vector<int> _nBezEl; // Number of Bezier poly. and for an el. - std::vector<fullMatrix<double> > _scaledNormEl; // Normals to 2D elements for Jacobian regularization and scaling - std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements -// inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; } - inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; } - void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ); -// void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl); - void calcScaledNormalEl2D(int iEl); }; @@ -48,26 +40,7 @@ template<class FuncType> void ObjContribScaledJac<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(); - - // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial - // Jacobians of 3D elements - if (_mesh->dim() == 2) { - _scaledNormEl.resize(_mesh->nEl()); -// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl); - for (int iEl = 0; iEl < _mesh->nEl(); iEl++) calcScaledNormalEl2D(iEl); - } - else { - _invStraightJac.resize(_mesh->nEl()); - double dumJac[3][3]; - for (int iEl = 0; iEl < _mesh->nEl(); iEl++) - _invStraightJac[iEl] = 1. / fabs(_mesh->el(iEl)->getPrimaryJacobian(0.,0.,0.,dumJac)); - } - + _mesh->initScaledJac(); updateMinMax(); FuncType::initialize(_min, _max); } @@ -80,14 +53,14 @@ bool ObjContribScaledJac<FuncType>::addContrib(double &Obj, alglib::real_1d_arra _max = -BIGVAL; for (int iEl = 0; iEl < _mesh->nEl(); iEl++) { - std::vector<double> sJ(_nBezEl[iEl]); // Scaled Jacobians - std::vector<double> gSJ(_nBezEl[iEl]*_mesh->nPCEl(iEl)); // Gradients of scaled Jacobians - scaledJacAndGradients(iEl, sJ, gSJ); - for (int l = 0; l < _nBezEl[iEl]; l++) { // Add contribution for each Bezier coeff. + std::vector<double> sJ(_mesh->nBezEl(iEl)); // Scaled Jacobians + std::vector<double> gSJ(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Gradients of scaled Jacobians + _mesh->scaledJacAndGradients(iEl, sJ, gSJ); + for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Add contribution for each Bezier coeff. Obj += _weight * FuncType::compute(sJ[l]); const double dfact = _weight * FuncType::computeDiff(sJ[l]); for (int iPC = 0; iPC < _mesh->nPCEl(iEl); iPC++) - gradObj[_mesh->indPCEl(iEl, iPC)] += dfact * gSJ[indGSJ(iEl, l, iPC)]; + gradObj[_mesh->indPCEl(iEl, iPC)] += dfact * gSJ[_mesh->indGSJ(iEl, l, iPC)]; _min = std::min(_min, sJ[l]); _max = std::max(_max, sJ[l]); } @@ -104,10 +77,10 @@ void ObjContribScaledJac<FuncType>::updateMinMax() _max = -BIGVAL; for (int iEl = 0; iEl < _mesh->nEl(); iEl++) { - std::vector<double> sJ(_nBezEl[iEl]); // Scaled Jacobians - std::vector<double> dumGSJ(_nBezEl[iEl]*_mesh->nPCEl(iEl)); // Gradients of scaled Jacobians - scaledJacAndGradients(iEl, sJ, dumGSJ); - for (int l = 0; l < _nBezEl[iEl]; l++) { // Check each Bezier coeff. + std::vector<double> sJ(_mesh->nBezEl(iEl)); // Scaled Jacobians + std::vector<double> dumGSJ(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Gradients of scaled Jacobians + _mesh->scaledJacAndGradients(iEl, sJ, dumGSJ); + for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Check each Bezier coeff. _min = std::min(_min, sJ[l]); _max = std::max(_max, sJ[l]); } @@ -123,97 +96,4 @@ void ObjContribScaledJac<FuncType>::updateResults(MeshOptResults &res) const } -// TODO: Re-introduce normal to geometry? -//void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl) -template<class FuncType> -void ObjContribScaledJac<FuncType>::calcScaledNormalEl2D(int iEl) -{ - const JacobianBasis *jac = _mesh->el(iEl)->getJacobianFuncSpace(); - - fullMatrix<double> primNodesXYZ(jac->getNumPrimMapNodes(), 3); -// SVector3 geoNorm(0.,0.,0.); -// std::map<MElement*,GEntity*>::const_iterator itEl2ent = element2entity.find(_el[iEl]); -// GEntity *ge = (itEl2ent == element2entity.end()) ? 0 : itEl2ent->second; -// const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization(); - for (int i=0; i<jac->getNumPrimMapNodes(); i++) { - const int &iV = _mesh->el2V(iEl, i); - primNodesXYZ(i,0) = _mesh->xyz(iV).x(); - primNodesXYZ(i,1) = _mesh->xyz(iV).y(); - primNodesXYZ(i,2) = _mesh->xyz(iV).z(); -// if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) { -// double u, v; -// _vert[iV]->getParameter(0,u); -// _vert[iV]->getParameter(1,v); -// geoNorm += ((GFace*)ge)->normal(SPoint2(u,v)); -// } - } -// if (hasGeoNorm && (geoNorm.normSq() == 0.)) { -// SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false); -// geoNorm = ((GFace*)ge)->normal(param); -// } - - fullMatrix<double> &elNorm = _scaledNormEl[iEl]; - elNorm.resize(1, 3); - const double norm = jac->getPrimNormal2D(primNodesXYZ, elNorm); - double factor = 1./norm; -// if (hasGeoNorm) { -// const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2); -// if (scal < 0.) factor = -factor; -// } - elNorm.scale(factor); // Re-scaling normal here is faster than an extra scaling operation on the Jacobian - -} - - -template<class FuncType> -void ObjContribScaledJac<FuncType>::scaledJacAndGradients(int iEl, std::vector<double> &sJ, - std::vector<double> &gSJ) -{ - const JacobianBasis *jacBasis = _mesh->el(iEl)->getJacobianFuncSpace(); - const int &numJacNodes = _nBezEl[iEl]; - const int &numMapNodes = _mesh->nNodEl(iEl); - fullMatrix<double> JDJ(numJacNodes, 3*numMapNodes+1), BDB(numJacNodes, 3*numMapNodes+1); - - // Coordinates of nodes - fullMatrix<double> nodesXYZ(numMapNodes, 3), normals(_mesh->dim(), 3); - for (int i = 0; i < numMapNodes; i++) { - const int &iVi = _mesh->el2V(iEl, i); - nodesXYZ(i, 0) = _mesh->xyz(iVi).x(); - nodesXYZ(i, 1) = _mesh->xyz(iVi).y(); - nodesXYZ(i, 2) = _mesh->xyz(iVi).z(); - } - - // Calculate Jacobian and gradients, scale if 3D (already scaled by - // regularization normals in 2D) - jacBasis->getSignedJacAndGradients(nodesXYZ, _scaledNormEl[iEl],JDJ); - if (_mesh->dim() == 3) JDJ.scale(_invStraightJac[iEl]); - - // Transform Jacobian and gradients from Lagrangian to Bezier basis - jacBasis->lag2Bez(JDJ,BDB); - - // Scaled jacobian - for (int l = 0; l < numJacNodes; l++) sJ[l] = BDB (l, 3*numMapNodes); - - // Gradients of the scaled jacobian - int iPC = 0; - std::vector<SPoint3> gXyzV(numJacNodes); - std::vector<SPoint3> gUvwV(numJacNodes); - for (int i = 0; i < numMapNodes; i++) { - const int &iFVi = _mesh->el2FV(iEl, i); - if (iFVi >= 0) { - for (int l = 0; l < numJacNodes; l++) - gXyzV[l] = SPoint3(BDB(l, i), BDB(l, i+numMapNodes), BDB(l, i+2*numMapNodes)); - _mesh->gXyz2gUvw(iFVi, gXyzV, gUvwV); - for (int l = 0; l < numJacNodes; l++) { - gSJ[indGSJ(iEl, l, iPC)] = gUvwV[l][0]; - if (_mesh->nPCFV(iFVi) >= 2) gSJ[indGSJ(iEl, l, iPC+1)] = gUvwV[l][1]; - if (_mesh->nPCFV(iFVi) == 3) gSJ[indGSJ(iEl, l, iPC+2)] = gUvwV[l][2]; - } - iPC += _mesh->nPCFV(iFVi); - } - } - -} - - #endif /* _OPTHOMOBJCONTRIBSCALEDJAC_H_ */ diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp index 877beae09aef80f69f0649c9fd797343de6a2342..c462487e0a8dc16d10f0a805eef71a55055b92e8 100644 --- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp +++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp @@ -753,10 +753,10 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p) #include "MeshOptimizer.h" -struct HOPatchParameters : public MeshOptParameters::PatchParameters +struct HOPatchDefParameters : public MeshOptParameters::PatchDefParameters { - HOPatchParameters(const OptHomParameters &p); - virtual ~HOPatchParameters() {} + HOPatchDefParameters(const OptHomParameters &p); + virtual ~HOPatchDefParameters() {} virtual double elBadness(const MElement *el); virtual double maxDistance(const MElement *el); private: @@ -765,7 +765,7 @@ private: }; -HOPatchParameters::HOPatchParameters(const OptHomParameters &p) +HOPatchDefParameters::HOPatchDefParameters(const OptHomParameters &p) { jacMin = p.BARRIER_MIN; jacMax = (p.BARRIER_MAX > 0.) ? p.BARRIER_MAX : 1.e300; @@ -784,14 +784,14 @@ HOPatchParameters::HOPatchParameters(const OptHomParameters &p) } -double HOPatchParameters::elBadness(const MElement *el) { +double HOPatchDefParameters::elBadness(const MElement *el) { double jmin, jmax; el->scaledJacRange(jmin, jmax); return std::min(jmin-jacMin, 0.) + std::min(jacMax-jmax, 0.); } -double HOPatchParameters::maxDistance(const MElement *el) { +double HOPatchDefParameters::maxDistance(const MElement *el) { return distanceFactor * el->maxDistToStraight(); } @@ -804,8 +804,8 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p) par.dim = p.dim; par.onlyVisible = p.onlyVisible; par.fixBndNodes = p.fixBndNodes; - HOPatchParameters patch(p); - par.patch = &patch; + HOPatchDefParameters patchDef(p); + par.patchDef = &patchDef; par.optDisplay = 30; par.verbose = 4; ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree); diff --git a/contrib/MeshOptimizer/MeshOptCommon.h b/contrib/MeshOptimizer/MeshOptCommon.h index 43137280d045b04369f936846fd940ed2346c4b8..b57db04592c10bc544acbd347a9c9e0db04e3727 100644 --- a/contrib/MeshOptimizer/MeshOptCommon.h +++ b/contrib/MeshOptimizer/MeshOptCommon.h @@ -56,7 +56,7 @@ struct MeshOptParameters { // Parameters controlling int optIterMax ; // Max. number of opt. iterations each time the barrier is moved int barrierIterMax ; // Max. number of times the barrier is moved }; - struct PatchParameters { + struct PatchDefParameters { int strategy; // Strategy: connected patches or adaptive one-by-one int minLayers, maxLayers; // Min. and max. nb. of layers around a bad element in patch union { @@ -73,7 +73,7 @@ struct MeshOptParameters { // Parameters controlling int dim ; // Which dimension to optimize bool onlyVisible ; // Apply optimization to visible entities ONLY bool fixBndNodes; // If points can move on boundaries - PatchParameters *patch; + PatchDefParameters *patchDef; std::vector<PassParameters> pass; int optDisplay; // Sampling rate in opt. iterations for display int verbose; // Level of information displayed and written to disk diff --git a/contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h b/contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h index b1829bb7ae0671e6929103b71f2c4a694d9cc631..2a452eaf4f43d2eb54d54eb119dc27e107cf397c 100644 --- a/contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h +++ b/contrib/MeshOptimizer/MeshOptObjContribScaledNodeDispSq.h @@ -25,9 +25,6 @@ public: protected: Patch *_mesh; double _weightFixed, _weightFree; - double _invLengthScaleSq; // Square inverse of a length for node displacement scaling - inline double scaledNodeDispSq(int iFV); - inline void gradScaledNodeDispSq(int iFV, std::vector<double> &gDSq); }; @@ -36,7 +33,7 @@ template<class FuncType> ObjContribScaledNodeDispSq<FuncType>::ObjContribScaledNodeDispSq(double weightFixed, double weightFree) : ObjContrib("ScaledNodeDispSq", FuncType::getNamePrefix()+"ScaledNodeDispSq"), - _mesh(0), _weightFixed(weightFixed), _weightFree(weightFree), _invLengthScaleSq(0.) + _mesh(0), _weightFixed(weightFixed), _weightFree(weightFree) { } @@ -45,22 +42,7 @@ template<class FuncType> void ObjContribScaledNodeDispSq<FuncType>::initialize(Patch *mesh) { _mesh = mesh; - - double maxDSq = 0.; - for (int iEl = 0; iEl < _mesh->nEl(); iEl++) { - const double d = _mesh->el(iEl)->maxDistToStraight(), dd = d*d; - if (dd > maxDSq) maxDSq = dd; - } - if (maxDSq < 1.e-10) { - double maxSSq = 0.; - for (int iEl = 0; iEl < _mesh->nEl(); iEl++) { - const double s = _mesh->el(iEl)->getOuterRadius(), ss = s*s; - if (ss > maxSSq) maxSSq = ss; - } - _invLengthScaleSq = 1./maxSSq; - } - else _invLengthScaleSq = 1./maxDSq; - + _mesh->initScaledNodeDispSq(); updateMinMax(); FuncType::initialize(_min, _max); } @@ -75,10 +57,10 @@ bool ObjContribScaledNodeDispSq<FuncType>::addContrib(double &Obj, for (int iFV = 0; iFV < _mesh->nFV(); iFV++) { const double Factor = _mesh->forced(iFV) ? _weightFixed : _weightFree; - const double dSq = scaledNodeDispSq(iFV); + const double dSq = _mesh->scaledNodeDispSq(iFV); Obj += Factor * FuncType::compute(dSq); std::vector<double> gDSq(_mesh->nPCFV(iFV)); - gradScaledNodeDispSq(iFV, gDSq); + _mesh->gradScaledNodeDispSq(iFV, gDSq); const double dfact = Factor * FuncType::computeDiff(dSq); for (int iPC = 0; iPC < _mesh->nPCFV(iFV); iPC++) gradObj[_mesh->indPCFV(iFV, iPC)] += dfact * gDSq[iPC]; @@ -97,7 +79,7 @@ void ObjContribScaledNodeDispSq<FuncType>::updateMinMax() _max = -BIGVAL; for (int iFV = 0; iFV < _mesh->nFV(); iFV++) { - const double dSq = scaledNodeDispSq(iFV); + const double dSq = _mesh->scaledNodeDispSq(iFV); _min = std::min(_min, dSq); _max = std::max(_max, dSq); } @@ -107,33 +89,10 @@ void ObjContribScaledNodeDispSq<FuncType>::updateMinMax() template<class FuncType> void ObjContribScaledNodeDispSq<FuncType>::updateResults(MeshOptResults &res) const { - const double scaleSq = 1./_invLengthScaleSq; + const double scaleSq = 1./_mesh->invLengthScaleSq(); res.minNodeDisp = std::min(sqrt(_min*scaleSq), res.minNodeDisp); res.maxNodeDisp = std::max(sqrt(_max*scaleSq), res.maxNodeDisp); } -template<class FuncType> -double ObjContribScaledNodeDispSq<FuncType>::scaledNodeDispSq(int iFV) -{ - const int &iV = _mesh->fv2V(iFV); - const SPoint3 d = _mesh->xyz(iV)-_mesh->ixyz(iV); - return (d[0]*d[0]+d[1]*d[1]+d[2]*d[2])*_invLengthScaleSq; -} - - -template<class FuncType> -void ObjContribScaledNodeDispSq<FuncType>::gradScaledNodeDispSq(int iFV, std::vector<double> &gDSq) -{ - const int &iV = _mesh->fv2V(iFV); - const SPoint3 gXyz = (_mesh->xyz(iV)-_mesh->ixyz(iV))*2.*_invLengthScaleSq; - SPoint3 gUvw; - _mesh->gXyz2gUvw(iFV, gXyz, gUvw); - - gDSq[0] = gUvw[0]; - if (_mesh->nPCFV(iFV) >= 2) gDSq[1] = gUvw[1]; - if (_mesh->nPCFV(iFV) == 3) gDSq[2] = gUvw[2]; -} - - #endif /* _MESHOPTOBJCONTRIBSCALEDNODEDISPSQ_H_ */ diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp index 8b9f9d585961eb180382899d0a2abc127c9e267d..d812dcf0a7bb96a682710c833046ca88eaf01c14 100644 --- a/contrib/MeshOptimizer/MeshOptPatch.cpp +++ b/contrib/MeshOptimizer/MeshOptPatch.cpp @@ -209,3 +209,289 @@ void Patch::writeMSH(const char *filename) fclose(f); } + +// TODO: allow user to choose type of scaling length +void Patch::initScaledNodeDispSq() +{ + if (_invLengthScaleSq == 0.) { + double maxDSq = 0.; + for (int iEl = 0; iEl < nEl(); iEl++) { + const double d = el(iEl)->maxDistToStraight(), dd = d*d; + if (dd > maxDSq) maxDSq = dd; + } + if (maxDSq < 1.e-10) { + double maxSSq = 0.; + for (int iEl = 0; iEl < nEl(); iEl++) { + const double s = el(iEl)->getOuterRadius(), ss = s*s; + if (ss > maxSSq) maxSSq = ss; + } + _invLengthScaleSq = 1./maxSSq; + } + else _invLengthScaleSq = 1./maxDSq; + } +} + + +double Patch::scaledNodeDispSq(int iFV) +{ + const int &iV = _fv2V[iFV]; + const SPoint3 d = _xyz[iV]-_ixyz[iV]; + return (d[0]*d[0]+d[1]*d[1]+d[2]*d[2])*_invLengthScaleSq; +} + + +void Patch::gradScaledNodeDispSq(int iFV, std::vector<double> &gDSq) +{ + const int &iV = _fv2V[iFV]; + const SPoint3 gXyz = (_xyz[iV]-_ixyz[iV])*2.*_invLengthScaleSq; + SPoint3 gUvw; + gXyz2gUvw(iFV, gXyz, gUvw); + + gDSq[0] = gUvw[0]; + if (_nPCFV[iFV] >= 2) gDSq[1] = gUvw[1]; + if (_nPCFV[iFV] == 3) gDSq[2] = gUvw[2]; +} + + +void Patch::initScaledJac() +{ + // Initialize _nBezEl + if (_nBezEl.empty()) { + _nBezEl.resize(nEl()); + for (int iEl=0; iEl<nEl(); iEl++) + _nBezEl[iEl] = _el[iEl]->getJacobianFuncSpace()->getNumJacNodes(); + } + + // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial + // Jacobians of 3D elements + if ((_dim == 2) && _scaledNormEl.empty()) { + _scaledNormEl.resize(nEl()); +// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl); + for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(iEl); + } + else if (_invStraightJac.empty()) { + _invStraightJac.resize(nEl(),1.); + double dumJac[3][3]; + for (int iEl = 0; iEl < nEl(); iEl++) + _invStraightJac[iEl] = 1. / fabs(_el[iEl]->getPrimaryJacobian(0.,0.,0.,dumJac)); + } +} + + +void Patch::initMetricMin() +{ + // Initialize _nBezEl + if (_nBezEl.empty()) { + _nBezEl.resize(nEl()); + for (int iEl=0; iEl<nEl(); iEl++) + _nBezEl[iEl] = _el[iEl]->getJacobianFuncSpace()->getNumJacNodes(); + } +} + + +// TODO: Re-introduce normal to geometry? +//void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl) +void Patch::calcScaledNormalEl2D(int iEl) +{ + const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace(); + + fullMatrix<double> primNodesXYZ(jac->getNumPrimMapNodes(),3); +// SVector3 geoNorm(0.,0.,0.); +// std::map<MElement*,GEntity*>::const_iterator itEl2ent = element2entity.find(_el[iEl]); +// GEntity *ge = (itEl2ent == element2entity.end()) ? 0 : itEl2ent->second; +// const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization(); + for (int i=0; i<jac->getNumPrimMapNodes(); i++) { + const int &iV = _el2V[iEl][i]; + primNodesXYZ(i,0) = _xyz[iV].x(); + primNodesXYZ(i,1) = _xyz[iV].y(); + primNodesXYZ(i,2) = _xyz[iV].z(); +// if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) { +// double u, v; +// _vert[iV]->getParameter(0,u); +// _vert[iV]->getParameter(1,v); +// geoNorm += ((GFace*)ge)->normal(SPoint2(u,v)); +// } + } +// if (hasGeoNorm && (geoNorm.normSq() == 0.)) { +// SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false); +// geoNorm = ((GFace*)ge)->normal(param); +// } + + fullMatrix<double> &elNorm = _scaledNormEl[iEl]; + elNorm.resize(1,3); + const double norm = jac->getPrimNormal2D(primNodesXYZ,elNorm); + double factor = 1./norm; +// if (hasGeoNorm) { +// const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2); +// if (scal < 0.) factor = -factor; +// } + elNorm.scale(factor); // Re-scaling normal here is faster than an extra scaling operation on the Jacobian + +} + + +void Patch::scaledJacAndGradients(int iEl, std::vector<double> &sJ, + std::vector<double> &gSJ) +{ + const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace(); + const int &numJacNodes = _nBezEl[iEl]; + const int &numMapNodes = _nNodEl[iEl]; + fullMatrix<double> JDJ(numJacNodes,3*numMapNodes+1), BDB(numJacNodes,3*numMapNodes+1); + + // Coordinates of nodes + fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3); + for (int i = 0; i < numMapNodes; i++) { + int &iVi = _el2V[iEl][i]; + nodesXYZ(i,0) = _xyz[iVi].x(); + nodesXYZ(i,1) = _xyz[iVi].y(); + nodesXYZ(i,2) = _xyz[iVi].z(); + } + + // Calculate Jacobian and gradients, scale if 3D (already scaled by + // regularization normals in 2D) + jacBasis->getSignedJacAndGradients(nodesXYZ,_scaledNormEl[iEl],JDJ); + if (_dim == 3) JDJ.scale(_invStraightJac[iEl]); + + // Transform Jacobian and gradients from Lagrangian to Bezier basis + jacBasis->lag2Bez(JDJ,BDB); + + // Scaled jacobian + for (int l = 0; l < numJacNodes; l++) sJ [l] = BDB (l,3*numMapNodes); + + // Gradients of the scaled jacobian + int iPC = 0; + std::vector<SPoint3> gXyzV(numJacNodes); + std::vector<SPoint3> gUvwV(numJacNodes); + for (int i = 0; i < numMapNodes; i++) { + int &iFVi = _el2FV[iEl][i]; + if (iFVi >= 0) { + for (int l = 0; l < numJacNodes; l++) + gXyzV [l] = SPoint3(BDB(l,i+0*numMapNodes), BDB(l,i+1*numMapNodes), + BDB(l,i+2*numMapNodes)); + _coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV); + for (int l = 0; l < numJacNodes; l++) { + gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0]; + if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1]; + if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2]; + } + iPC += _nPCFV[iFVi]; + } + } +} + + +void Patch::metricMinAndGradients(int iEl, std::vector<double> &lambda, + std::vector<double> &gradLambda) +{ + const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace(); + const int &numJacNodes = jacBasis->getNumJacNodes(); + const int &numMapNodes = jacBasis->getNumMapNodes(); + const int &numPrimMapNodes = jacBasis->getNumPrimMapNodes(); + fullVector<double> lambdaJ(numJacNodes), lambdaB(numJacNodes); + fullMatrix<double> gradLambdaJ(numJacNodes, 2 * numMapNodes); + fullMatrix<double> gradLambdaB(numJacNodes, 2 * numMapNodes); + + // Coordinates of nodes + fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZStraight(numPrimMapNodes,3); + for (int i = 0; i < numMapNodes; i++) { + int &iVi = _el2V[iEl][i]; + nodesXYZ(i,0) = _xyz[iVi].x(); + nodesXYZ(i,1) = _xyz[iVi].y(); + nodesXYZ(i,2) = _xyz[iVi].z(); + if (i < numPrimMapNodes) { + nodesXYZStraight(i,0) = _ixyz[iVi].x(); + nodesXYZStraight(i,1) = _ixyz[iVi].y(); + nodesXYZStraight(i,2) = _ixyz[iVi].z(); + } + } + + jacBasis->getMetricMinAndGradients(nodesXYZ,nodesXYZStraight,lambdaJ,gradLambdaJ); + + //l2b.mult(lambdaJ, lambdaB); + //l2b.mult(gradLambdaJ, gradLambdaB); + lambdaB = lambdaJ; + gradLambdaB = gradLambdaJ; + + int iPC = 0; + std::vector<SPoint3> gXyzV(numJacNodes); + std::vector<SPoint3> gUvwV(numJacNodes); + for (int l = 0; l < numJacNodes; l++) { + lambda[l] = lambdaB(l); + } + for (int i = 0; i < numMapNodes; i++) { + int &iFVi = _el2FV[iEl][i]; + if (iFVi >= 0) { + for (int l = 0; l < numJacNodes; l++) { + gXyzV [l] = SPoint3(gradLambdaB(l,i+0*numMapNodes), + gradLambdaB(l,i+1*numMapNodes),/*BDB(l,i+2*nbNod)*/ 0.); + } + _coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV); + for (int l = 0; l < numJacNodes; l++) { + gradLambda[indGSJ(iEl,l,iPC)] = gUvwV[l][0]; + if (_nPCFV[iFVi] >= 2) gradLambda[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1]; + if (_nPCFV[iFVi] == 3) gradLambda[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2]; + } + iPC += _nPCFV[iFVi]; + } + } +} + + +//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 2c1d43df71818411d4b32586be49a521f5bc1e52..2207c839bba3a458e5de564e53259cd5894879c4 100644 --- a/contrib/MeshOptimizer/MeshOptPatch.h +++ b/contrib/MeshOptimizer/MeshOptPatch.h @@ -77,6 +77,23 @@ public: void updateGEntityPositions(); void writeMSH(const char *filename); + // Node distance measure + void initScaledNodeDispSq(); + inline double invLengthScaleSq() { return _invLengthScaleSq; } + double scaledNodeDispSq(int iFV); + void gradScaledNodeDispSq(int iFV, std::vector<double> &gDSq); + + // High-order: scaled Jacobian and metric measures + 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); + private: // Mesh entities and variables @@ -113,6 +130,16 @@ private: { return indJB3DBase(_nNodEl[iEl],l,i,j,m); } + + // Node displacement + double _invLengthScaleSq; // Square inverse of a length for node displacement scaling + + // High-order: scaled Jacobian and metric measures + std::vector<int> _nBezEl; // Number of Bezier poly. and for an el. + std::vector<fullMatrix<double> > _scaledNormEl; // Normals to 2D elements for Jacobian regularization and scaling + std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements +// void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl); + void calcScaledNormalEl2D(int iEl); }; diff --git a/contrib/MeshOptimizer/MeshOptimizer.cpp b/contrib/MeshOptimizer/MeshOptimizer.cpp index 52f7a37cf21a1b06e9218ea0391cd24783f8c94a..f672741ad0eb1e2120570384ffac684788ac2987 100644 --- a/contrib/MeshOptimizer/MeshOptimizer.cpp +++ b/contrib/MeshOptimizer/MeshOptimizer.cpp @@ -189,9 +189,9 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn std::vector<std::set<MElement*> > primPatches; primPatches.reserve(badElements.size()); for (std::set<MElement*>::const_iterator it = badElements.begin(); it != badElements.end(); ++it) { - const double limDist = par.patch->maxDistance(*it); - primPatches.push_back(getSurroundingPatch(*it, par.patch->minLayers, - par.patch->maxLayers, limDist, vertex2elements)); + const double limDist = par.patchDef->maxDistance(*it); + primPatches.push_back(getSurroundingPatch(*it, par.patchDef->minLayers, + par.patchDef->maxLayers, limDist, vertex2elements)); } // Compute patch connectivity @@ -203,7 +203,7 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn for(std::set<MElement*>::iterator itEl = patch.begin(); itEl != patch.end(); ++itEl) { std::set<int> &patchInd = tags[*itEl]; if (!patchInd.empty() && (badElements.find(*itEl) != badElements.end() || - !par.patch->weakMerge)) { + !par.patchDef->weakMerge)) { for (std::set<int>::iterator itBS = patchInd.begin(); itBS != patchInd.end(); ++itBS) patchConnect[*itBS].insert(iB); patchConnect[iB].insert(patchInd.begin(), patchInd.end()); @@ -304,7 +304,7 @@ static MElement *getWorstElement(std::set<MElement*> &badElts, const MeshOptPara MElement *worstEl = 0; for (std::set<MElement*>::iterator it=badElts.begin(); it!=badElts.end(); it++) { - const double val = par.patch->elBadness(*it); + const double val = par.patchDef->elBadness(*it); if (val < worst) { worst = val; worstEl = *it; @@ -335,16 +335,16 @@ static void optimizeOneByOne badElts.erase(worstEl); // Initialize patch size to be adapted - int maxLayers = par.patch->maxLayers; + int maxLayers = par.patchDef->maxLayers; double distanceFactor = 1.; int success; // Patch adaptation loop - for (int iAdapt=0; iAdapt<par.patch->maxAdaptPatch; iAdapt++) { + for (int iAdapt=0; iAdapt<par.patchDef->maxAdaptPatch; iAdapt++) { // Set up patch - const double limDist = par.patch->maxDistance(worstEl); - std::set<MElement*> toOptimizePrim = getSurroundingPatch(worstEl, par.patch->minLayers, + const double limDist = par.patchDef->maxDistance(worstEl); + std::set<MElement*> toOptimizePrim = getSurroundingPatch(worstEl, par.patchDef->minLayers, maxLayers, limDist, vertex2elements); std::set<MVertex*> toFix = getAllBndVertices(toOptimizePrim, vertex2elements); std::set<MElement*> toOptimize; @@ -379,15 +379,15 @@ static void optimizeOneByOne } // If (partial) success, update mesh and break adaptation loop, otherwise adapt - if ((success > 0) || (iAdapt == par.patch->maxAdaptPatch-1)) { + if ((success > 0) || (iAdapt == par.patchDef->maxAdaptPatch-1)) { opt.updateResults(res); if (success >= 0) { opt.patch.updateGEntityPositions(); break; } else { - distanceFactor *= par.patch->distanceAdaptFact; - maxLayers *= par.patch->maxLayersAdaptFact; + distanceFactor *= par.patchDef->distanceAdaptFact; + maxLayers *= par.patchDef->maxLayersAdaptFact; if (par.verbose > 1) Msg::Info("Patch %i failed (adapt #%i), adapting with increased size", iBadEl, iAdapt); } @@ -442,16 +442,16 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par, MeshOptResults &res) // if (jmin < par.patchJac.min || jmax > par.patchJac.max) badElts.insert(el); // } // } - if ((el->getDim() == par.dim) && (par.patch->elBadness(el) < 0.)) badElts.insert(el); + if ((el->getDim() == par.dim) && (par.patchDef->elBadness(el) < 0.)) badElts.insert(el); } } - if (par.patch->strategy == MeshOptParameters::STRAT_CONNECTED) + if (par.patchDef->strategy == MeshOptParameters::STRAT_CONNECTED) optimizeConnectedPatches(vertex2elements, element2entity, badElts, par, res); - else if (par.patch->strategy == MeshOptParameters::STRAT_ONEBYONE) + else if (par.patchDef->strategy == MeshOptParameters::STRAT_ONEBYONE) optimizeOneByOne(vertex2elements, element2entity, badElts, par, res); else - Msg::Error("Unknown strategy %d for mesh optimization", par.patch->strategy); + Msg::Error("Unknown strategy %d for mesh optimization", par.patchDef->strategy); if (par.verbose > 0) { if (res.success == 1)