diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 1c5d75080094bae84e4b8c7cc7d78dce2d418fb3..2d54c47d4f4af6e5b5f2beb36576fa21de212417 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -322,9 +322,9 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
}
case 3:
{
- dJ = fabs(jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
- jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
- jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
+ dJ = (jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
+ jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
+ jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
break;
}
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 84daa9e986e1a3668235ea92b7ee634a0482a1f2..2ba8c4753815fc4712bd539a81149c8aea63db9d 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -31,9 +31,11 @@ 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);
+ // mesh.scaledJac(iEl,sJ);
std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl)); // Gradients of scaled Jacobians
- mesh.gradScaledJac(iEl,gSJ);
+ // mesh.gradScaledJac(iEl,gSJ);
+ mesh.scaledJacAndGradients (iEl,sJ,gSJ);
+
for (int l = 0; l < mesh.nBezEl(iEl); l++) {
Obj += compute_f(sJ[l]);
const double f1 = compute_f1(sJ[l]);
@@ -198,9 +200,9 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
alglib::mincgreport rep;
mincgcreate(x, state);
alglib::real_1d_array scale;
- calcScale(scale);
- mincgsetscale(state,scale);
- mincgsetprecscale(state);
+ calcScale(scale);
+ mincgsetscale(state,scale);
+ mincgsetprecscale(state);
mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
mincgsetxrep(state, true);
alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
@@ -299,8 +301,8 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
OptHomMessage("Optimization done Range (%g,%g)",minJac,maxJac);
- if (minJac < barrier_min && maxJac > barrier_max) return 0;
- if (minJac > 0.0) return 1;
+ if (minJac > barrier_min && maxJac < barrier_max) return 1;
+ if (minJac > 0.0) return 0;
return -1;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 1525a0a9510b91d04a0878182781b2cd831b6fe8..15b1283ea5236292b5e86696d06c5d1e8108994d 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -9,8 +9,26 @@
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
+ 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);
+ return allDPsi;
+}
+
std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
{
@@ -109,6 +127,8 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
if (_jacBez.find(lagrange->type) == _jacBez.end()) {
_jacBez[lagrange->type] = computeJB(lagrange, bezier);
+ _gradShapeFunctions[lagrange->type] = computeGSF(lagrange, bezier);
+ _lag2Bez[lagrange->type] = bezier->matrixLag2Bez;
}
_nBezEl[iEl] = bezier->points.size1();
_nNodEl[iEl] = lagrange->points.size1();
@@ -301,6 +321,120 @@ void Mesh::updateGEntityPositions()
}
+/*
+ A faster version that computes jacobians and their gradients
+
+ Terms of jacobian are of the form J_{11} = dX / dU
+
+ X = \sum \phi_j X_{j}
+
+ d J_{11} / dX_k = d \phi_k / du
+
+ */
+
+
+
+void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<double> &gSJ) {
+
+ if (_dim == 2 && !projJac){
+ // printf("coucou\n");
+ gradScaledJac(iEl,gSJ);
+ scaledJac(iEl,sJ);
+ return;
+ }
+
+ SVector3 n ;
+ if (_dim == 2)
+ n = _normEl[iEl];
+
+ // std::vector<double> OLD = sJ;
+ // std::vector<double> OLD = gSJ;
+ // gradScaledJac(iEl,OLD);
+ // scaledJac(iEl,OLD);
+
+ fullMatrix<double> &gsf = _gradShapeFunctions[_el[iEl]->getTypeForMSH()];
+ const fullMatrix<double> &l2b = _lag2Bez[_el[iEl]->getTypeForMSH()];
+ const int nbBez = _nBezEl[iEl];
+ const int nbNod = _nNodEl[iEl];
+ fullMatrix<double> JDJ (nbBez,3*nbNod+1);
+ fullMatrix<double> BDB (nbBez,3*nbNod+1);
+
+ double jac[3][3];
+ for (int l = 0; l < nbBez; l++) {
+ fullMatrix<double> dPsi(gsf, l * 3, 3);
+ jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+ jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+ jac[2][0] = jac[2][1] = jac[2][2] = (_dim == 2) ? n.z() : 0.0;
+ const double INVJ = (_dim == 3) ? _invStraightJac[iEl] : 1.0;
+ for (int i = 0; i < nbNod; i++) {
+ int &iVi = _el2V[iEl][i];
+ const double x = _xyz[iVi].x();
+ const double y = _xyz[iVi].y();
+ const double z = _xyz[iVi].z();
+ for (int k = 0; k < 3; k++) {
+ const double gg = dPsi(i, k);
+ jac[k][0] += x * gg;
+ jac[k][1] += y * gg;
+ jac[k][2] += z * gg;
+ }
+ }
+ for (int i = 0; i < nbNod; i++) {
+ JDJ (l,i+0*nbNod) =
+ ( dPsi(i, 0) * jac[1][1] * jac[2][2] + jac[0][2] * dPsi(i, 1) * jac[2][1] +
+ jac[0][1] * jac[1][2] * dPsi(i, 2) - jac[0][2] * jac[1][1] * dPsi(i, 2) -
+ dPsi(i, 0) * jac[1][2] * jac[2][1] - jac[0][1] * dPsi(i, 1) * jac[2][2])
+ * INVJ;
+ JDJ (l,i+1*nbNod) =
+ (jac[0][0] * dPsi(i, 1) * jac[2][2] + jac[0][2] * jac[1][0] * dPsi(i, 2) +
+ dPsi(i, 0) * jac[1][2] * jac[2][0] - jac[0][2] * dPsi(i, 1) * jac[2][0] -
+ jac[0][0] * jac[1][2] * dPsi(i, 2) - dPsi(i, 0) * jac[1][0] * jac[2][2])
+ * INVJ;
+ JDJ (l,i+2*nbNod) =
+ (jac[0][0] * jac[1][1] * dPsi(i, 2) + dPsi(i, 0) * jac[1][0] * jac[2][1] +
+ jac[0][1] * dPsi(i, 1) * jac[2][0] - dPsi(i, 0) * jac[1][1] * jac[2][0] -
+ jac[0][0] * dPsi(i, 1) * jac[2][1] - jac[0][1] * jac[1][0] * dPsi(i, 2))
+ * INVJ;
+ }
+ const double dJ =
+ jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
+ jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
+ jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2];
+ JDJ(l,3*nbNod) = dJ * INVJ;
+ }
+
+ // (N_b x N_b) x (N_b x 3*N_n + 1)
+ l2b.mult(JDJ,BDB);
+
+ // the scaled jacobian
+ // printf("ELEMENT %d\n",iEl);
+ for (int l = 0; l < nbBez; l++) {
+ sJ [l] = BDB (l,3*nbNod);
+ // printf("OLD %12.5E NEW %12.5E\n",OLD[l],sJ[l]);
+ }
+
+ // gradients of the scaled jacobian
+ int iPC = 0;
+ std::vector<SPoint3> gXyzV(nbBez);
+ std::vector<SPoint3> gUvwV(nbBez);
+ for (int i = 0; i < nbNod; i++) {
+ int &iFVi = _el2FV[iEl][i];
+ if (iFVi >= 0) {
+ for (int l = 0; l < nbBez; l++) {
+ gXyzV [l] = SPoint3(BDB(l,i+0*nbNod),BDB(l,i+1*nbNod),BDB(l,i+2*nbNod));
+ }
+ _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
+ for (int l = 0; l < nbBez; 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];
+ // if (fabs(OLD[indGSJ(iEl,l,iPC+2)]-gSJ[indGSJ(iEl,l,iPC+2)]) > 1.e-6)
+ // printf("OLD = %12.5E new = %12.5E\n",OLD[indGSJ(iEl,l,iPC+2)],gSJ[indGSJ(iEl,l,iPC+2)]);
+ }
+ iPC += _nPCFV[iFVi];
+ }
+ }
+}
+
void Mesh::scaledJac(int iEl, std::vector<double> &sJ)
{
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index dcdfd6fc4814c7c8ad4b82057c993a2d1e665e93..25f194beac99eaed9ec48752354c5fb0fed3c4d4 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -32,6 +32,8 @@ public:
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; }
inline double distSq(int iFV);
@@ -80,11 +82,14 @@ private:
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); }
static inline int indJB3DBase(int nNod, int l, int i, int j, int m) { return ((l*nNod+i)*nNod+j)*nNod+m; }
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 8e8ee194785bef4bfac0ddc5a93c3a2b2eb1bf4e..61fb0f0965bf2e4b89ece5425afb63a01a997ab5 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -374,7 +374,7 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
clock_t t1 = clock();
- int samples = 100;
+ int samples = 30;
// int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_FIXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
@@ -425,19 +425,20 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
OptHOM temp(*itf, toOptimizeSplit[i], toFix, method);
temp.recalcJacDist();
temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
- OptHomMessage("Optimizing a blob %i/%i composed of %4d elements minJ %12.5E -- maxJ %12.5E", i, toOptimizeSplit.size(), toOptimizeSplit[i].size(), minJac, maxJac);
+ OptHomMessage("Optimizing a blob %i/%i composed of %4d elements minJ %12.5E -- maxJ %12.5E", i+1, toOptimizeSplit.size(), toOptimizeSplit[i].size(), minJac, maxJac);
p.SUCCESS = std::min(p.SUCCESS,temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax));
+ temp.mesh.updateGEntityPositions();
}
}
else while (1){
std::set<MVertex*> toFix;
- std::set<MElement*> toOptimize;
+ std::set<MElement*> toOptimize;
toFix = filter2D_boundaryLayer(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor, badasses, toOptimize);
OptHOM temp(*itf, toOptimize, toFix, method);
temp.recalcJacDist();
- temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+ temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
// if (minJac < 1.e2)OptHomMessage("Optimizing a blob of %4d elements minJ %12.5E -- maxJ %12.5E",(*itf)->getNumMeshElements(), minJac, maxJac);
std::ostringstream ossI;
ossI << "initial_" << (*itf)->tag() << "ITER_" << ITER << ".msh";
@@ -476,7 +477,21 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
int ITER = 0;
Msg::Info("Model region %4d is considered",(*itr)->tag());
p.SUCCESS = true;
- while (1){
+ if (p.filter == 0) {
+ std::set<MVertex*> toFix;
+ std::set<MElement*> toOptimize;
+ toFix = filterSimple(*itr, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, toOptimize);
+ std::vector<std::set<MElement*> > toOptimizeSplit = splitConnex(toOptimize);
+ for (int i = 0; i < toOptimizeSplit.size(); ++i) {
+ OptHOM temp(*itr, toOptimizeSplit[i], toFix, method);
+ temp.recalcJacDist();
+ temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+ OptHomMessage("Optimizing a blob %i/%i composed of %4d elements minJ %12.5E -- maxJ %12.5E", i+1, toOptimizeSplit.size(), toOptimizeSplit[i].size(), minJac, maxJac);
+ p.SUCCESS = std::min(p.SUCCESS,temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax));
+ temp.mesh.updateGEntityPositions();
+ }
+ }
+ else while (1){
std::set<MVertex*> toFix;
std::set<MElement*> toOptimize;