diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index c9de40f03e48492fd4a8dc1a7953e9158400081b..b398a48ff366aa747e4453685be3be66249d3b7c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -22,6 +22,8 @@ OptHOM::OptHOM(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &
{
};
+
+
// Contribution of the element Jacobians to the objective function value and gradients (2D version)
bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
{
@@ -31,9 +33,7 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
- // mesh.scaledJac(iEl,sJ);
std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl)); // Gradients of scaled Jacobians
- // mesh.gradScaledJac(iEl,gSJ);
mesh.scaledJacAndGradients (iEl,sJ,gSJ);
for (int l = 0; l < mesh.nBezEl(iEl); l++) {
@@ -127,8 +127,9 @@ void OptHOM::recalcJacDist()
minJac = 1.e300;
maxJac = -1.e300;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
- std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
- mesh.scaledJac(iEl,sJ);
+ std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
+ std::vector<double> dumGSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl)); // (Dummy) gradients of scaled Jacobians
+ mesh.scaledJacAndGradients (iEl,sJ,dumGSJ);
for (int l = 0; l < mesh.nBezEl(iEl); l++) {
minJac = std::min(minJac, sJ[l]);
maxJac = std::max(maxJac, sJ[l]);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index ebe132ac949df7e5ccc66ce126b10be21692689a..23c635348d9d9542958af57c10aa9c9f737757ec 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -8,11 +8,11 @@
-std::map<int, std::vector<double> > Mesh::_jacBez;
std::map<int, fullMatrix<double> > Mesh::_gradShapeFunctions;
std::map<int, fullMatrix<double> > Mesh::_lag2Bez;
+
fullMatrix<double> Mesh::computeGSF(const polynomialBasis *lagrange, const bezierBasis *bezier)
{
// bezier points are defined in the [0,1] x [0,1] quad
@@ -30,59 +30,6 @@ fullMatrix<double> Mesh::computeGSF(const polynomialBasis *lagrange, const bezie
}
-std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
-{
- int nbNodes = lagrange->points.size1();
-
- // bezier points are defined in the [0,1] x [0,1] quad
- fullMatrix<double> bezierPoints = bezier->points;
- if (lagrange->parentType == TYPE_QUA) {
- for (int i = 0; i < bezierPoints.size1(); ++i) {
- bezierPoints(i, 0) = -1 + 2 * bezierPoints(i, 0);
- bezierPoints(i, 1) = -1 + 2 * bezierPoints(i, 1);
- }
- }
-
- fullMatrix<double> allDPsi;
- lagrange->df(bezierPoints, allDPsi);
- int size = bezier->points.size1();
- std::vector<double> JB;
- for (int d = 0; d < lagrange->dimension; ++d) {
- size *= nbNodes;
- }
- JB.resize(size, 0.);
- for (int k = 0; k < bezier->points.size1(); ++k) {
- fullMatrix<double> dPsi(allDPsi, k * 3, 3);
- if (lagrange->dimension == 2) {
- for (int i = 0; i < nbNodes; i++) {
- for (int j = 0; j < nbNodes; j++) {
- double Jij = dPsi(i, 0) * dPsi(j, 1) - dPsi(i, 1) * dPsi(j,0);
- for (int l = 0; l < bezier->points.size1(); l++) {
- JB[indJB2DBase(nbNodes,l,i,j)] += bezier->matrixLag2Bez(l, k) * Jij;
- }
- }
- }
- }
- if (lagrange->dimension == 3) {
- for (int i = 0; i < nbNodes; i++) {
- for (int j = 0; j < nbNodes; j++) {
- for (int m = 0; m < nbNodes; m++) {
- double Jijm =
- (dPsi(j, 1) * dPsi(m, 2) - dPsi(j, 2) * dPsi(m, 1)) * dPsi(i, 0)
- + (dPsi(j, 2) * dPsi(m, 0) - dPsi(j, 0) * dPsi(m, 2)) * dPsi(i, 1)
- + (dPsi(j, 0) * dPsi(m, 1) - dPsi(j, 1) * dPsi(m, 0)) * dPsi(i, 2);
- for (int l = 0; l < bezier->points.size1(); l++) {
- JB[indJB3DBase(nbNodes,l,i,j,m)] += bezier->matrixLag2Bez(l, k) * Jijm;
- }
- }
- }
- }
- }
- }
- return JB;
-}
-
-
Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFix, int method) :
_ge(ge)
@@ -125,8 +72,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
_el[iEl] = el;
const polynomialBasis *lagrange = el->getFunctionSpace();
const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
- if (_jacBez.find(lagrange->type) == _jacBez.end()) {
- _jacBez[lagrange->type] = computeJB(lagrange, bezier);
+ if (_lag2Bez.find(lagrange->type) == _lag2Bez.end()) {
_gradShapeFunctions[lagrange->type] = computeGSF(lagrange, bezier);
_lag2Bez[lagrange->type] = bezier->matrixLag2Bez;
}
@@ -440,145 +386,6 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<
}
-void Mesh::scaledJac(int iEl, std::vector<double> &sJ)
-{
- const std::vector<double> &jacBez = _jacBez[_el[iEl]->getTypeForMSH()];
- if (_dim == 2) {
- SVector3 &n = _normEl[iEl];
- if (projJac) {
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- sJ[l] = 0.;
- for (int i = 0; i < _nNodEl[iEl]; i++) {
- int &iVi = _el2V[iEl][i];
- for (int j = 0; j < _nNodEl[iEl]; j++) {
- int &iVj = _el2V[iEl][j];
- sJ[l] += jacBez[indJB2D(iEl,l,i,j)]
- * (_xyz[iVi].x() * _xyz[iVj].y() * n.z() - _xyz[iVi].x() * _xyz[iVj].z() * n.y()
- + _xyz[iVi].y() * _xyz[iVj].z() * n.x());
- }
- }
- }
- }
- else
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- sJ[l] = 0.;
- for (int i = 0; i < _nNodEl[iEl]; i++) {
- int &iVi = _el2V[iEl][i];
- for (int j = 0; j < _nNodEl[iEl]; j++) {
- int &iVj = _el2V[iEl][j];
- sJ[l] += jacBez[indJB2D(iEl,l,i,j)] * _xyz[iVi].x() * _xyz[iVj].y();
- }
- }
- sJ[l] *= n.z();
- }
- }
- else {
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- sJ[l] = 0.;
- for (int i = 0; i < _nNodEl[iEl]; i++) {
- int &iVi = _el2V[iEl][i];
- for (int j = 0; j < _nNodEl[iEl]; j++) {
- int &iVj = _el2V[iEl][j];
- for (int m = 0; m < _nNodEl[iEl]; m++) {
- int &iVm = _el2V[iEl][m];
- sJ[l] += jacBez[indJB3D(iEl,l,i,j,m)] * _xyz[iVi].x() * _xyz[iVj].y() * _xyz[iVm].z();
- }
- }
- }
- sJ[l] *= _invStraightJac[iEl];
- }
- }
-
-}
-
-
-
-void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
-{
- const std::vector<double> &jacBez = _jacBez[_el[iEl]->getTypeForMSH()];
- if (_dim == 2) {
- int iPC = 0;
- SVector3 n = _normEl[iEl];
- if (projJac) {
- for (int i = 0; i < _nNodEl[iEl]; i++) {
- int &iFVi = _el2FV[iEl][i];
- if (iFVi >= 0) {
- std::vector<SPoint3> gXyzV(_nBezEl[iEl],SPoint3(0.,0.,0.));
- std::vector<SPoint3> gUvwV(_nBezEl[iEl]);
- for (int m = 0; m < _nNodEl[iEl]; m++) {
- int &iVm = _el2V[iEl][m];
- const double vpx = _xyz[iVm].y() * n.z() - _xyz[iVm].z() * n.y();
- const double vpy = -_xyz[iVm].x() * n.z() + _xyz[iVm].z() * n.x();
- const double vpz = _xyz[iVm].x() * n.y() - _xyz[iVm].y() * n.x();
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- gXyzV[l][0] += jacBez[indJB2D(iEl,l,i,m)] * vpx;
- gXyzV[l][1] += jacBez[indJB2D(iEl,l,i,m)] * vpy;
- gXyzV[l][2] += jacBez[indJB2D(iEl,l,i,m)] * vpz;
- }
- }
- _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
- if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
- }
- iPC += _nPCFV[iFVi];
- }
- }
- }
- else
- for (int i = 0; i < _nNodEl[iEl]; i++) {
- int &iFVi = _el2FV[iEl][i];
- if (iFVi >= 0) {
- std::vector<SPoint3> gXyzV(_nBezEl[iEl],SPoint3(0.,0.,0.));
- std::vector<SPoint3> gUvwV(_nBezEl[iEl]);
- for (int m = 0; m < _nNodEl[iEl]; m++) {
- int &iVm = _el2V[iEl][m];
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- gXyzV[l][0] += jacBez[indJB2D(iEl,l,i,m)] * _xyz[iVm].y() * n.z();
- gXyzV[l][1] += jacBez[indJB2D(iEl,l,m,i)] * _xyz[iVm].x() * n.z();
- }
- }
- _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
- if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
- }
- iPC += _nPCFV[iFVi];
- }
- }
- }
- else {
- int iPC = 0;
- for (int i = 0; i < _nNodEl[iEl]; i++) {
- int &iFVi = _el2FV[iEl][i];
- if (iFVi >= 0) {
- std::vector<SPoint3> gXyzV(_nBezEl[iEl],SPoint3(0.,0.,0.));
- std::vector<SPoint3> gUvwV(_nBezEl[iEl]);
- for (int a = 0; a < _nNodEl[iEl]; a++) {
- int &iVa = _el2V[iEl][a];
- for (int b = 0; b < _nNodEl[iEl]; b++) {
- int &iVb = _el2V[iEl][b];
- for (int l = 0; l < _nBezEl[iEl]; l++) {
- gXyzV[l][0] += jacBez[indJB3D(iEl,l,i,a,b)] * _xyz[iVa].y() * _xyz[iVb].z() * _invStraightJac[iEl];
- gXyzV[l][1] += jacBez[indJB3D(iEl,l,a,i,b)] * _xyz[iVa].x() * _xyz[iVb].z() * _invStraightJac[iEl];
- gXyzV[l][2] += jacBez[indJB3D(iEl,l,a,b,i)] * _xyz[iVa].x() * _xyz[iVb].y() * _invStraightJac[iEl];
- }
- }
- }
- _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
- for (int l = 0; l < _nBezEl[iEl]; 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 Mesh::pcScale(int iFV, std::vector<double> &scale)
{
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 25f194beac99eaed9ec48752354c5fb0fed3c4d4..bac927e1a3139010324bfe48c3296cd07359f980 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -30,9 +30,6 @@ public:
inline const int &indPCEl(int iEl, int iPC) { return _indPCEl[iEl][iPC]; }
inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
- void scaledJac(int iEl, std::vector<double> &sJ);
- void gradScaledJac(int iEl, std::vector<double> &gSJ);
- // a faster version that computes both jacobians and their gradients
void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; }
@@ -81,14 +78,12 @@ private:
ParamCoord *_pc;
bool projJac; // Using "projected" Jacobians or not
- static std::map<int, std::vector<double> > _jacBez;
static std::map<int, fullMatrix<double> > _gradShapeFunctions; // gradients of shape functions at Bezier points
static std::map<int, fullMatrix<double> > _lag2Bez; // gradients of shape functions at Bezier points
int addVert(MVertex* vert);
int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
SVector3 getNormalEl(int iEl);
- static std::vector<double> computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier);
static fullMatrix<double> computeGSF(const polynomialBasis *lagrange, const bezierBasis *bezier);
static inline int indJB2DBase(int nNod, int l, int i, int j) { return (l*nNod+i)*nNod+j; }
inline int indJB2D(int iEl, int l, int i, int j) { return indJB2DBase(_nNodEl[iEl],l,i,j); }
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
index bfe1fc18fb688c2f951edfbfe2ee2a98c54be18d..cdb3e3c8db07d182b61f63b46123bd7cd430a5d0 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
@@ -4,6 +4,8 @@
#include "MVertex.h"
#include "ParamCoord.h"
+
+
ParamCoordSurf::ParamCoordSurf(GEntity *ge)
{
if ((ge->dim() == 2) && (ge->geomType() != GEntity::DiscreteSurface)) _gf = static_cast<GFace*>(ge);
@@ -171,6 +173,8 @@ void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::v
}
+
+
void ParamCoordParent::exportParamCoord(MVertex *v, const SPoint3 &uvw)
{
for (int d = 0; d < v->onWhat()->dim(); ++d) {