diff --git a/contrib/HighOrderMeshOptimizer/CMakeLists.txt b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
index f7311c8442876ce9bbe0db1d63d3f2279f75b934..29223ca3cd54a09fcdeafb1d3fa363c22560a162 100644
--- a/contrib/HighOrderMeshOptimizer/CMakeLists.txt
+++ b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
@@ -8,6 +8,7 @@ set(SRC
OptHOM.cpp
OptHomRun.cpp
OptHomIntegralBoundaryDist.cpp
+ OptHomCADDist.cpp
ParamCoord.cpp
SuperEl.cpp
OptHomElastic.cpp
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index b255bcc8a6346e71451d565ec9be09590931b7dd..1b2f0307c019a5d2fbdb3a1e2ccf2fe840e5d2a3 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -36,7 +36,8 @@ static int NEVAL = 0;
#include "OptHomRun.h"
#include "GmshMessage.h"
#include "GmshConfig.h"
-#include "ConjugateGradients.h"
+#include "OptHomCADDist.h"
+#include "MElement.h"
#include "MLine.h"
#include "MTriangle.h"
#include "GModel.h"
@@ -112,35 +113,6 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
return true;
}
-bool OptHOM::addJacObjGrad(double &Obj, std::vector<double> &gradObj)
-{
- minJac = 1.e300;
- maxJac = -1.e300;
-
- for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
- // Scaled Jacobians
- std::vector<double> sJ(mesh.nBezEl(iEl));
- // Gradients of scaled Jacobians
- std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));
- mesh.scaledJacAndGradients (iEl,sJ,gSJ);
-
- for (int l = 0; l < mesh.nBezEl(iEl); l++) {
- double f1 = compute_f1(sJ[l], jacBar);
- Obj += compute_f(sJ[l], jacBar);
- if (_optimizeBarrierMax) {
- Obj += compute_f(sJ[l], barrier_min);
- f1 += compute_f1(sJ[l], barrier_min);
- }
- for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++)
- gradObj[mesh.indPCEl(iEl,iPC)] += f1*gSJ[mesh.indGSJ(iEl,l,iPC)];
- minJac = std::min(minJac, sJ[l]);
- maxJac = std::max(maxJac, sJ[l]);
- }
- }
-
- return true;
-}
-
bool OptHOM::addApproximationErrorObjGrad(double factor, double &Obj, alglib::real_1d_array &gradObj, simpleFunction<double>& fct)
{
@@ -173,162 +145,6 @@ static void computeGradSFAtNodes (MElement *e, std::vector<std::vector<SVector3>
}
}
-static double MFaceGFaceDistance (MTriangle *t, GFace *gf, std::vector<std::vector<SVector3> > *gsfT=0, std::map<MVertex*,SVector3> *normalsToCAD=0) {
- const double h = t->maxEdge();
- double jac[3][3];
- double distFace = 0.0;
- // for (int j=0;j<3;j++){
- for (int j=0;j<t->getNumVertices();j++){
- // get parametric coordinates of jth vertex
- // the last line of the jacobian is the normal
- // to the element @ (u_mesh,v_mesh)
-
- if (gsfT){
- double detJ = t->getJacobian((*gsfT)[j],jac);
- }
- else{
- const nodalBasis &elbasis = *t->getFunctionSpace();
- double u_mesh = elbasis.points(j,0);
- double v_mesh = elbasis.points(j,1);
- double detJ = t->getJacobian(u_mesh,v_mesh,0,jac);
- }
-
- SVector3 tg_mesh (jac[2][0],jac[2][1],jac[2][2]);
- tg_mesh.normalize();
-
- SVector3 tg_cad ;
- if (normalsToCAD)tg_cad = (*normalsToCAD)[t->getVertex(j)];
- else {
- SPoint2 p_cad;
- reparamMeshVertexOnFace(t->getVertex (j), gf, p_cad);
- tg_cad = gf->normal(p_cad);
- tg_cad.normalize();
- }
- SVector3 diff1 = (dot(tg_cad, tg_mesh) > 0) ?
- tg_cad - tg_mesh : tg_cad + tg_mesh;
- // printf("%g %g %g vs %g %g %g\n",tg_cad.x(),tg_cad.y(),tg_cad.z(),tg_mesh.x(),tg_mesh.y(),tg_mesh.z());
- distFace += diff1.norm();
- }
- return distFace;
-}
-
-static double MLineGEdgeDistance (MLine *l, GEdge *ge, FILE *f = 0) {
- const nodalBasis &elbasis = *l->getFunctionSpace();
- const double h = .25*0.5*distance (l->getVertex (0), l->getVertex (1) ) / (l->getNumVertices()-1);
- double jac[3][3];
- double distEdge = 0.0;
-
- // if(f)printf("%d\n",l->getNumVertices());
-
- for (int j=0;j<l->getNumVertices();j++){
- double t_mesh = elbasis.points(j,0);
- // if (f) printf("%g ",t_mesh);
- double detJ = l->getJacobian(t_mesh,0,0,jac);
- SVector3 tg_mesh (jac[0][0],jac[0][1],jac[0][2]);
- tg_mesh.normalize();
- double t_cad;
- reparamMeshVertexOnEdge(l->getVertex (j), ge, t_cad);
- SVector3 tg_cad = ge->firstDer(t_cad);
- tg_cad.normalize();
-
- SVector3 diff1 = (dot(tg_cad, tg_mesh) > 0) ?
- tg_cad - tg_mesh : tg_cad + tg_mesh;
-
- if (f){
- fprintf (f,"SP(%g,%g,%g){%g};\n",l->getVertex (j)->x(),
- l->getVertex (j)->y(),l->getVertex (j)->z(),h*diff1.norm());
- }
-
- // SVector3 n = crossprod(tg_cad,tg_mesh);
- // printf("%g %g vs %g %g\n",tg_cad.x(),tg_cad.y(),tg_mesh.x(),tg_mesh.y());
- distEdge += diff1.norm();
- }
- // if(f)printf("\n");
- return h*distEdge;
-}
-
-
-void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,double> &distances){
-
- std::map<MEdge,double,Less_Edge> dist2Edge;
- for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
- if ((*it)->geomType() == GEntity::Line)continue;
- for (unsigned int i=0;i<(*it)->lines.size(); i++){
- double d = MLineGEdgeDistance ( (*it)->lines[i] , *it );
- MEdge e = (*it)->lines[i]->getEdge(0);
- dist2Edge[e] = d;
- }
- }
-
- // printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
-
- std::map<MFace,double,Less_Face> dist2Face;
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
- if ((*it)->geomType() == GEntity::Plane)continue;
- for (unsigned int i=0;i<(*it)->triangles.size(); i++){
- double d = MFaceGFaceDistance ( (*it)->triangles[i] , *it );
- MFace f = (*it)->triangles[i]->getFace(0);
- dist2Face[f] = d;
- }
- }
-
- std::vector<GEntity*> entities;
- gm->getEntities(entities);
- for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
- GEntity* &entity = entities[iEnt];
- if (entity->dim() != dim) continue;
- for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
- MElement *element = entity->getMeshElement(iEl);
- double d = 0.0;
- for (int iEdge = 0; iEdge < element->getNumEdges(); ++iEdge) {
- MEdge e = element->getEdge(iEdge);
- std::map<MEdge,double,Less_Edge>::iterator it = dist2Edge.find(e);
- if(it != dist2Edge.end())d+=it->second;
- }
- for (int iFace = 0; iFace < element->getNumFaces(); ++iFace) {
- MFace f = element->getFace(iFace);
- std::map<MFace,double,Less_Face>::iterator it = dist2Face.find(f);
- if(it != dist2Face.end())d+=it->second;
- }
- distances[element] = d;
- }
- }
-}
-
-
-double distanceToGeometry(GModel *gm)
-{
-
- FILE *f = fopen("toto.pos","w");
- fprintf(f,"View \"\"{\n");
-
- double Obj = 0.0;
-
- for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
- if ((*it)->geomType() == GEntity::Line)continue;
- for (unsigned int i=0;i<(*it)->lines.size(); i++){
- //Obj += MLineGEdgeDistance ( (*it)->lines[i] , *it,f );
- Obj = std::max(MLineGEdgeDistance ( (*it)->lines[i] , *it, f ),Obj);
- }
- }
-
- printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
-
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
- if ((*it)->geomType() == GEntity::Plane)continue;
- // printf("FACE %d with %d triangles\n",(*it)->tag(),(*it)->triangles.size());
- for (unsigned int i=0;i<(*it)->triangles.size(); i++){
- //Obj += MFaceGFaceDistance ( (*it)->triangles[i] , *it );
- Obj = std::max(Obj,MFaceGFaceDistance ( (*it)->triangles[i] , *it ));
- }
- }
-
- printf("DISTANCE TO GEOMETRY : 1D AND 2D PART %22.15E\n",Obj);
- fprintf(f,"};\n");
- fclose(f);
- return Obj;
-}
-
bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gradObj)
{
@@ -355,21 +171,21 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
for (unsigned int i=0;i<(*it)->lines.size(); i++){
doWeCompute[i] = false;
for (unsigned int j=0;j<(*it)->lines[i]->getNumVertices(); j++){
- int index = mesh.getFreeVertexStartIndex((*it)->lines[i]->getVertex(j));
- if (index >=0){
- doWeCompute[i] = true;
- continue;
- }
+ int index = mesh.getFreeVertexStartIndex((*it)->lines[i]->getVertex(j));
+ if (index >=0){
+ doWeCompute[i] = true;
+ continue;
+ }
}
}
std::vector<double> dist((*it)->lines.size());
for (unsigned int i=0;i<(*it)->lines.size(); i++){
if (doWeCompute[i]){
- // compute the distance from the geometry to the mesh
- dist[i] = MLineGEdgeDistance ( (*it)->lines[i] , *it );
- maxDistCAD = std::max(maxDistCAD,dist[i]);
- distCAD += dist [i] * factor;
+ // compute the distance from the geometry to the mesh
+ dist[i] = MLineGEdgeDistance ( (*it)->lines[i] , *it );
+ maxDistCAD = std::max(maxDistCAD,dist[i]);
+ distCAD += dist [i] * factor;
}
}
// be clever to compute the derivative : iterate on all
@@ -379,50 +195,50 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
const double eps = 1.e-6;
for (unsigned int i=0;i<(*it)->lines.size(); i++){
if (doWeCompute[i]){
- for (int j=2 ; j<(*it)->lines[i]->getNumVertices() ; j++){
- MVertex *v = (*it)->lines[i]->getVertex(j);
- int index = mesh.getFreeVertexStartIndex(v);
- // printf("%d %d (%d %d)\n",v->getNum(),index,v->onWhat()->tag(),v->onWhat()->dim());
- if (index >= 0){
- double t;
- v->getParameter(0,t);
- SPoint3 pp (v->x(),v->y(),v->z());
- GPoint gp = (*it)->point(t+eps);
- v->setParameter(0,t+eps);
- v->setXYZ(gp.x(),gp.y(),gp.z());
- double dist2 = MLineGEdgeDistance ( (*it)->lines[i] , *it );
- double deriv = (dist2 - dist[i])/eps;
- v->setXYZ(pp.x(),pp.y(),pp.z());
- v->setParameter(0,t);
- // printf("%g %g %g\n",dist[i],dist2, MLineGEdgeDistance ( (*it)->lines[i] , *it ));
- // get the index of the vertex
- gradObj[index] += deriv * factor;
- }
- }
+ for (int j=2 ; j<(*it)->lines[i]->getNumVertices() ; j++){
+ MVertex *v = (*it)->lines[i]->getVertex(j);
+ int index = mesh.getFreeVertexStartIndex(v);
+ // printf("%d %d (%d %d)\n",v->getNum(),index,v->onWhat()->tag(),v->onWhat()->dim());
+ if (index >= 0){
+ double t;
+ v->getParameter(0,t);
+ SPoint3 pp (v->x(),v->y(),v->z());
+ GPoint gp = (*it)->point(t+eps);
+ v->setParameter(0,t+eps);
+ v->setXYZ(gp.x(),gp.y(),gp.z());
+ double dist2 = MLineGEdgeDistance ( (*it)->lines[i] , *it );
+ double deriv = (dist2 - dist[i])/eps;
+ v->setXYZ(pp.x(),pp.y(),pp.z());
+ v->setParameter(0,t);
+ // printf("%g %g %g\n",dist[i],dist2, MLineGEdgeDistance ( (*it)->lines[i] , *it ));
+ // get the index of the vertex
+ gradObj[index] += deriv * factor;
+ }
+ }
}
// printf("done\n");
// For a low order vertex classified on the GEdge, we recompute
// two distances for the two MLines connected to the vertex
for (unsigned int i=0;i<(*it)->lines.size()-1; i++){
- MVertex *v = (*it)->lines[i]->getVertex(1);
- int index = mesh.getFreeVertexStartIndex(v);
- if (index >= 0){
- double t;
- v->getParameter(0,t);
- SPoint3 pp (v->x(),v->y(),v->z());
- GPoint gp = (*it)->point(t+eps);
- v->setParameter(0,t+eps);
- v->setXYZ(gp.x(),gp.y(),gp.z());
- MLine *l1 = (*it)->lines[i];
- MLine *l2 = (*it)->lines[i+1];
- // printf("%d %d -- %d %d\n",l1->getVertex(0)->getNum(),l1->getVertex(1)->getNum(),l2->getVertex(0)->getNum(),l2->getVertex(1)->getNum());
- double deriv =
- (MLineGEdgeDistance ( l1 , *it ) - dist[i]) /eps +
- (MLineGEdgeDistance ( l2 , *it ) - dist[i+1])/eps;
- v->setXYZ(pp.x(),pp.y(),pp.z());
- v->setParameter(0,t);
- gradObj[index] += deriv * factor;
- }
+ MVertex *v = (*it)->lines[i]->getVertex(1);
+ int index = mesh.getFreeVertexStartIndex(v);
+ if (index >= 0){
+ double t;
+ v->getParameter(0,t);
+ SPoint3 pp (v->x(),v->y(),v->z());
+ GPoint gp = (*it)->point(t+eps);
+ v->setParameter(0,t+eps);
+ v->setXYZ(gp.x(),gp.y(),gp.z());
+ MLine *l1 = (*it)->lines[i];
+ MLine *l2 = (*it)->lines[i+1];
+ // printf("%d %d -- %d %d\n",l1->getVertex(0)->getNum(),l1->getVertex(1)->getNum(),l2->getVertex(0)->getNum(),l2->getVertex(1)->getNum());
+ double deriv =
+ (MLineGEdgeDistance ( l1 , *it ) - dist[i]) /eps +
+ (MLineGEdgeDistance ( l2 , *it ) - dist[i+1])/eps;
+ v->setXYZ(pp.x(),pp.y(),pp.z());
+ v->setParameter(0,t);
+ gradObj[index] += deriv * factor;
+ }
}
}
}
@@ -445,32 +261,32 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
for (unsigned int i=0;i<(*it)->triangles.size(); i++){
doWeCompute[i] = false;
for (unsigned int j=0;j<(*it)->triangles[i]->getNumVertices(); j++){
- int index = mesh.getFreeVertexStartIndex((*it)->triangles[i]->getVertex(j));
- if (index >=0){
- doWeCompute[i] = true;
- }
+ int index = mesh.getFreeVertexStartIndex((*it)->triangles[i]->getVertex(j));
+ if (index >=0){
+ doWeCompute[i] = true;
+ }
}
if (doWeCompute[i]){
- for (unsigned int j=0;j<(*it)->triangles[i]->getNumVertices(); j++){
- MVertex *v = (*it)->triangles[i]->getVertex(j);
- if (normalsToCAD.find(v) == normalsToCAD.end()){
- SPoint2 p_cad;
- reparamMeshVertexOnFace(v, *it, p_cad);
- SVector3 tg_cad = (*it)->normal(p_cad);
- tg_cad.normalize();
- normalsToCAD[v] = tg_cad;
- }
- }
+ for (unsigned int j=0;j<(*it)->triangles[i]->getNumVertices(); j++){
+ MVertex *v = (*it)->triangles[i]->getVertex(j);
+ if (normalsToCAD.find(v) == normalsToCAD.end()){
+ SPoint2 p_cad;
+ reparamMeshVertexOnFace(v, *it, p_cad);
+ SVector3 tg_cad = (*it)->normal(p_cad);
+ tg_cad.normalize();
+ normalsToCAD[v] = tg_cad;
+ }
+ }
}
}
for (unsigned int i=0;i<(*it)->triangles.size(); i++){
// compute the distance from the geometry to the mesh
if(doWeCompute[i]){
- const double d = MFaceGFaceDistance ( (*it)->triangles[i] , *it , &gsfT, &normalsToCAD);
- dist[(*it)->triangles[i]] = d;
- maxDistCAD = std::max(maxDistCAD,d);
- distCAD += d * factor;
+ const double d = MFaceGFaceDistance ( (*it)->triangles[i] , *it , &gsfT, &normalsToCAD);
+ dist[(*it)->triangles[i]] = d;
+ maxDistCAD = std::max(maxDistCAD,d);
+ distCAD += d * factor;
}
}
@@ -478,57 +294,57 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
const double eps = 1.e-6;
for (unsigned int i=0;i<(*it)->triangles.size(); i++){
if(doWeCompute[i]){
- for (unsigned int j=0;j<(*it)->triangles[i]->getNumVertices(); j++){
- // for (; itm !=v2t.end(); ++itm){
- MVertex *v = (*it)->triangles[i]->getVertex(j);
- if(v->onWhat()->dim() == 1){
- int index = mesh.getFreeVertexStartIndex(v);
- if (index >= 0){
- MTriangle *t = (*it)->triangles[i];
- GEdge *ge = v->onWhat()->cast2Edge();
- double t_;
- v->getParameter(0,t_);
- SPoint3 pp (v->x(),v->y(),v->z());
- GPoint gp = ge->point(t_+eps);
- v->setParameter(0,t_+eps);
- v->setXYZ(gp.x(),gp.y(),gp.z());
- const double distT = dist[t];
- double deriv = (MFaceGFaceDistance ( t , *it , &gsfT, &normalsToCAD) - distT) /eps;
- v->setXYZ(pp.x(),pp.y(),pp.z());
- v->setParameter(0,t_);
- gradObj[index] += deriv * factor;
- }
- }
-
- if(v->onWhat() == *it){
- int index = mesh.getFreeVertexStartIndex(v);
- if (index >= 0){
- MTriangle *t = (*it)->triangles[i];
- double uu,vv;
- v->getParameter(0,uu);
- v->getParameter(1,vv);
- SPoint3 pp (v->x(),v->y(),v->z());
-
- const double distT = dist[t];
-
- GPoint gp = (*it)->point(uu+eps,vv);
- v->setParameter(0,uu+eps);
- v->setXYZ(gp.x(),gp.y(),gp.z());
- double deriv = (MFaceGFaceDistance ( t , *it, &gsfT, &normalsToCAD ) - distT) /eps;
- v->setXYZ(pp.x(),pp.y(),pp.z());
- v->setParameter(0,uu);
- gradObj[index] += deriv * factor;
-
- gp = (*it)->point(uu,vv+eps);
- v->setParameter(1,vv+eps);
- v->setXYZ(gp.x(),gp.y(),gp.z());
- deriv = (MFaceGFaceDistance ( t , *it, &gsfT, &normalsToCAD ) - distT) /eps;
- v->setXYZ(pp.x(),pp.y(),pp.z());
- v->setParameter(1,vv);
- gradObj[index+1] += deriv * factor;
- }
- }
- }
+ for (unsigned int j=0;j<(*it)->triangles[i]->getNumVertices(); j++){
+ // for (; itm !=v2t.end(); ++itm){
+ MVertex *v = (*it)->triangles[i]->getVertex(j);
+ if(v->onWhat()->dim() == 1){
+ int index = mesh.getFreeVertexStartIndex(v);
+ if (index >= 0){
+ MTriangle *t = (*it)->triangles[i];
+ GEdge *ge = v->onWhat()->cast2Edge();
+ double t_;
+ v->getParameter(0,t_);
+ SPoint3 pp (v->x(),v->y(),v->z());
+ GPoint gp = ge->point(t_+eps);
+ v->setParameter(0,t_+eps);
+ v->setXYZ(gp.x(),gp.y(),gp.z());
+ const double distT = dist[t];
+ double deriv = (MFaceGFaceDistance ( t , *it , &gsfT, &normalsToCAD) - distT) /eps;
+ v->setXYZ(pp.x(),pp.y(),pp.z());
+ v->setParameter(0,t_);
+ gradObj[index] += deriv * factor;
+ }
+ }
+
+ if(v->onWhat() == *it){
+ int index = mesh.getFreeVertexStartIndex(v);
+ if (index >= 0){
+ MTriangle *t = (*it)->triangles[i];
+ double uu,vv;
+ v->getParameter(0,uu);
+ v->getParameter(1,vv);
+ SPoint3 pp (v->x(),v->y(),v->z());
+
+ const double distT = dist[t];
+
+ GPoint gp = (*it)->point(uu+eps,vv);
+ v->setParameter(0,uu+eps);
+ v->setXYZ(gp.x(),gp.y(),gp.z());
+ double deriv = (MFaceGFaceDistance ( t , *it, &gsfT, &normalsToCAD ) - distT) /eps;
+ v->setXYZ(pp.x(),pp.y(),pp.z());
+ v->setParameter(0,uu);
+ gradObj[index] += deriv * factor;
+
+ gp = (*it)->point(uu,vv+eps);
+ v->setParameter(1,vv+eps);
+ v->setXYZ(gp.x(),gp.y(),gp.z());
+ deriv = (MFaceGFaceDistance ( t , *it, &gsfT, &normalsToCAD ) - distT) /eps;
+ v->setXYZ(pp.x(),pp.y(),pp.z());
+ v->setParameter(1,vv);
+ gradObj[index+1] += deriv * factor;
+ }
+ }
+ }
}
}
}
@@ -539,53 +355,6 @@ bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gr
return true;
}
-bool OptHOM::addBndObjGrad2(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;
-
-}
-
-
-bool OptHOM::addBndObjGrad(double factor, double &Obj, std::vector<double> &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;
-
-}
-
bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
@@ -640,30 +409,6 @@ bool OptHOM::addDistObjGrad(double Fact, double &Obj,
return true;
}
-bool OptHOM::addDistObjGrad(double Fact, double &Obj,
- std::vector<double> &gradObj)
-{
- maxDist = 0;
- avgDist = 0;
- int nbBnd = 0;
-
- for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
- const double dSq = mesh.distSq(iFV), dist = sqrt(dSq);
- Obj += Fact * dSq;
- std::vector<double> gDSq(mesh.nPCFV(iFV));
- mesh.gradDistSq(iFV,gDSq);
- for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++)
- gradObj[mesh.indPCFV(iFV,iPC)] += Fact*gDSq[iPC];
- maxDist = std::max(maxDist, dist);
- avgDist += dist;
- nbBnd++;
- }
- if (nbBnd != 0) avgDist /= nbBnd;
-
- return true;
-}
-
-
// FIXME TEST
// Assume a unit square centered on 0,0
// fct is
@@ -678,25 +423,6 @@ bool OptHOM::addDistObjGrad(double Fact, double &Obj,
// }
//};
-// Gmsh's (cheaper) version of the optimizer
-void OptHOM::evalObjGrad(std::vector<double> &x,
- double &Obj,
- bool gradsNeeded,
- std::vector<double> &gradObj)
-{
- mesh.updateMesh(x.data());
- Obj = 0.;
- for (unsigned int i = 0; i < gradObj.size(); i++) gradObj[i] = 0.;
-
- /// control Jacobians
- addJacObjGrad(Obj, gradObj);
- /// Control distance to the straight sided mesh
- addDistObjGrad(lambda, Obj, gradObj);
- if(_optimizeCAD)
- addBndObjGrad(lambda3, Obj, gradObj);
-
-}
-
void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj,
alglib::real_1d_array &gradObj)
@@ -734,12 +460,6 @@ void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj,
(static_cast<OptHOM*>(HOInst))->evalObjGrad(x, Obj, gradObj);
}
-void evalObjGradFunc(std::vector<double> &x, double &Obj, bool needGrad,
- std::vector<double> &gradObj, void *HOInst)
-{
- (static_cast<OptHOM*>(HOInst))->evalObjGrad(x, Obj, true, gradObj);
-}
-
void OptHOM::recalcJacDist()
{
@@ -803,13 +523,6 @@ void OptHOM::calcScale(alglib::real_1d_array &scale)
}
}
-void OptHOM::OptimPass(std::vector<double> &x, int itMax)
-{
- Msg::Info("--- In-house Optimization pass with initial jac. range (%g, %g), jacBar = %g",
- minJac, maxJac, jacBar);
- GradientDescent (evalObjGradFunc, x, this);
-}
-
void OptHOM::OptimPass(alglib::real_1d_array &x, int itMax)
{
@@ -978,100 +691,4 @@ int OptHOM::optimize(double weight, double weightCAD, double b_min,
return -1;
}
-
-int OptHOM::optimize_inhouse(double weight, double weightCAD, double b_min,
- double b_max, bool optimizeMetricMin, int pInt,
- int itMax, int optPassMax, int optCAD, double distanceMax, double tolerance)
-{
- barrier_min = b_min;
- barrier_max = b_max;
- distance_max = distanceMax;
- progressInterv = pInt;
-// powM = 4;
-// powP = 3;
-
- _optimizeMetricMin = optimizeMetricMin;
- _optimizeCAD = optCAD;
- // Set weights & length scale for non-dimensionalization
- lambda = weight;
- lambda3 = weightCAD;
- geomTol = tolerance;
- std::vector<double> dSq(mesh.nEl());
- mesh.distSqToStraight(dSq);
- const double maxDSq = *max_element(dSq.begin(),dSq.end());
- if (maxDSq < 1.e-10) { // Length scale for non-dim. distance
- std::vector<double> sSq(mesh.nEl());
- mesh.elSizeSq(sSq);
- const double maxSSq = *max_element(sSq.begin(),sSq.end());
- invLengthScaleSq = 1./maxSSq;
- }
- else invLengthScaleSq = 1./maxDSq;
-
- // Set initial guess
- std::vector<double> x(mesh.nPC());
- mesh.getUvw(x.data());
-
- // Calculate initial performance
- recalcJacDist();
- initMaxDist = maxDist;
- initAvgDist = avgDist;
-
- const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
- jacBar = jacBarStart;
-
- _optimizeBarrierMax = false;
- // Calculate initial objective function value and gradient
- initObj = 0.;
- std::vector<double>gradObj(mesh.nPC());
- for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
- evalObjGrad(x, initObj, true, gradObj);
-
- Msg::Info("Start optimizing %d elements (%d vertices, %d free vertices, %d variables) "
- "with min barrier %g and max. barrier %g", mesh.nEl(), mesh.nVert(),
- mesh.nFV(), mesh.nPC(), barrier_min, barrier_max);
-
- int ITER = 0;
- bool minJacOK = true;
-
- while (minJac < barrier_min || (maxDistCAD > distance_max && _optimizeCAD)) {
- const double startMinJac = minJac;
- OptimPass(x, itMax);
- recalcJacDist();
- jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
- if (_optimizeCAD) jacBar = std::min(jacBar,barrier_min);
- if (ITER ++ > optPassMax) {
- minJacOK = (minJac > barrier_min && (maxDistCAD < distance_max || !_optimizeCAD));
- break;
- }
- if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
- Msg::Info("Stagnation in minJac detected, stopping optimization");
- minJacOK = false;
- break;
- }
- }
-
- ITER = 0;
- if (minJacOK && (!_optimizeMetricMin)) {
- _optimizeBarrierMax = true;
- jacBar = 1.1 * maxJac;
- while (maxJac > barrier_max ) {
- const double startMaxJac = maxJac;
- OptimPass(x, itMax);
- recalcJacDist();
- jacBar = 1.1 * maxJac;
- if (ITER ++ > optPassMax) break;
- if (fabs((maxJac-startMaxJac)/startMaxJac) < 0.01) {
- Msg::Info("Stagnation in maxJac detected, stopping optimization");
- break;
- }
- }
- }
-
- Msg::Info("Optimization done Range (%g,%g)", minJac, maxJac);
-
- if (minJac > barrier_min && maxJac < barrier_max) return 1;
- if (minJac > 0.0) return 0;
- return -1;
-}
-
#endif
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index a4927fd7a6596fd01ac085f47e7b9ba0557b3262..7e1c61448989aaf39752ab60fb92538f4f8b8a0b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -54,15 +54,9 @@ public:
int optimize(double lambda, double lambda3, double barrier_min, double barrier_max,
bool optimizeMetricMin, int pInt, int itMax, int optPassMax,
int optimizeCAD, double optCADDistMax, double tolerance);
- int optimize_inhouse(double weight, double weightCAD, double b_min, double b_max,
- bool optimizeMetricMin, int pInt, int itMax, int optPassMax,
- int optCAD, double distanceMax, double tolerance);
void recalcJacDist();
inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
void updateMesh(const alglib::real_1d_array &x);
- void evalObjGrad(std::vector<double> &x, double &Obj, bool gradsNeeded,
- std::vector<double> &gradObj);
-
void evalObjGrad(const alglib::real_1d_array &x, double &Obj,
alglib::real_1d_array &gradObj);
void printProgress(const alglib::real_1d_array &x, double Obj);
@@ -80,16 +74,11 @@ public:
// true : minimize the distance between mesh and CAD
bool addApproximationErrorObjGrad(double Fact, double &Obj, alglib::real_1d_array &gradObj, simpleFunction<double>& fct);
bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
- bool addJacObjGrad(double &Obj, std::vector<double> &);
bool addBndObjGrad (double Fact, double &Obj, alglib::real_1d_array &gradObj);
- bool addBndObjGrad2(double Fact, double &Obj, alglib::real_1d_array &gradObj);
- bool addBndObjGrad(double Fact, double &Obj, std::vector<double> &gradObj);
bool addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj);
- bool addDistObjGrad(double Fact, double &Obj, std::vector<double> &gradObj);
bool addDistObjGrad(double Fact, double &Obj, alglib::real_1d_array &gradObj);
void calcScale(alglib::real_1d_array &scale);
void OptimPass(alglib::real_1d_array &x, int itMax);
- void OptimPass(std::vector<double> &x, int itMax);
};
void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp b/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4644d77435aa66665fc54eef96b489d76a0fc7ad
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
@@ -0,0 +1,168 @@
+// TODO: Header
+
+#include "MLine.h"
+#include "MTriangle.h"
+#include "GModel.h"
+#include "OptHomCADDist.h"
+
+
+double MFaceGFaceDistance(MTriangle *t, GFace *gf,
+ std::vector<std::vector<SVector3> > *gsfT,
+ std::map<MVertex*,SVector3> *normalsToCAD) {
+ const double h = t->maxEdge();
+ double jac[3][3];
+ double distFace = 0.0;
+ // for (int j=0;j<3;j++){
+ for (int j=0;j<t->getNumVertices();j++){
+ // get parametric coordinates of jth vertex
+ // the last line of the jacobian is the normal
+ // to the element @ (u_mesh,v_mesh)
+
+ if (gsfT){
+ double detJ = t->getJacobian((*gsfT)[j],jac);
+ }
+ else{
+ const nodalBasis &elbasis = *t->getFunctionSpace();
+ double u_mesh = elbasis.points(j,0);
+ double v_mesh = elbasis.points(j,1);
+ double detJ = t->getJacobian(u_mesh,v_mesh,0,jac);
+ }
+
+ SVector3 tg_mesh (jac[2][0],jac[2][1],jac[2][2]);
+ tg_mesh.normalize();
+
+ SVector3 tg_cad ;
+ if (normalsToCAD)tg_cad = (*normalsToCAD)[t->getVertex(j)];
+ else {
+ SPoint2 p_cad;
+ reparamMeshVertexOnFace(t->getVertex (j), gf, p_cad);
+ tg_cad = gf->normal(p_cad);
+ tg_cad.normalize();
+ }
+ SVector3 diff1 = (dot(tg_cad, tg_mesh) > 0) ?
+ tg_cad - tg_mesh : tg_cad + tg_mesh;
+ // printf("%g %g %g vs %g %g %g\n",tg_cad.x(),tg_cad.y(),tg_cad.z(),tg_mesh.x(),tg_mesh.y(),tg_mesh.z());
+ distFace += diff1.norm();
+ }
+ return distFace;
+}
+
+
+double MLineGEdgeDistance (MLine *l, GEdge *ge, FILE *f) {
+ const nodalBasis &elbasis = *l->getFunctionSpace();
+ const double h = .25*0.5*distance (l->getVertex (0), l->getVertex (1) ) / (l->getNumVertices()-1);
+ double jac[3][3];
+ double distEdge = 0.0;
+
+ // if(f)printf("%d\n",l->getNumVertices());
+
+ for (int j=0;j<l->getNumVertices();j++){
+ double t_mesh = elbasis.points(j,0);
+ // if (f) printf("%g ",t_mesh);
+ double detJ = l->getJacobian(t_mesh,0,0,jac);
+ SVector3 tg_mesh (jac[0][0],jac[0][1],jac[0][2]);
+ tg_mesh.normalize();
+ double t_cad;
+ reparamMeshVertexOnEdge(l->getVertex (j), ge, t_cad);
+ SVector3 tg_cad = ge->firstDer(t_cad);
+ tg_cad.normalize();
+
+ SVector3 diff1 = (dot(tg_cad, tg_mesh) > 0) ?
+ tg_cad - tg_mesh : tg_cad + tg_mesh;
+
+ if (f){
+ fprintf (f,"SP(%g,%g,%g){%g};\n",l->getVertex (j)->x(),
+ l->getVertex (j)->y(),l->getVertex (j)->z(),h*diff1.norm());
+ }
+
+ // SVector3 n = crossprod(tg_cad,tg_mesh);
+ // printf("%g %g vs %g %g\n",tg_cad.x(),tg_cad.y(),tg_mesh.x(),tg_mesh.y());
+ distEdge += diff1.norm();
+ }
+ // if(f)printf("\n");
+ return h*distEdge;
+}
+
+
+void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,double> &distances){
+
+ std::map<MEdge,double,Less_Edge> dist2Edge;
+ for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
+ if ((*it)->geomType() == GEntity::Line)continue;
+ for (unsigned int i=0;i<(*it)->lines.size(); i++){
+ double d = MLineGEdgeDistance ( (*it)->lines[i] , *it );
+ MEdge e = (*it)->lines[i]->getEdge(0);
+ dist2Edge[e] = d;
+ }
+ }
+
+ // printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
+
+ std::map<MFace,double,Less_Face> dist2Face;
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
+ if ((*it)->geomType() == GEntity::Plane)continue;
+ for (unsigned int i=0;i<(*it)->triangles.size(); i++){
+ double d = MFaceGFaceDistance ( (*it)->triangles[i] , *it );
+ MFace f = (*it)->triangles[i]->getFace(0);
+ dist2Face[f] = d;
+ }
+ }
+
+ std::vector<GEntity*> entities;
+ gm->getEntities(entities);
+ for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
+ GEntity* &entity = entities[iEnt];
+ if (entity->dim() != dim) continue;
+ for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
+ MElement *element = entity->getMeshElement(iEl);
+ double d = 0.0;
+ for (int iEdge = 0; iEdge < element->getNumEdges(); ++iEdge) {
+ MEdge e = element->getEdge(iEdge);
+ std::map<MEdge,double,Less_Edge>::iterator it = dist2Edge.find(e);
+ if(it != dist2Edge.end())d+=it->second;
+ }
+ for (int iFace = 0; iFace < element->getNumFaces(); ++iFace) {
+ MFace f = element->getFace(iFace);
+ std::map<MFace,double,Less_Face>::iterator it = dist2Face.find(f);
+ if(it != dist2Face.end())d+=it->second;
+ }
+ distances[element] = d;
+ }
+ }
+}
+
+
+double distanceToGeometry(GModel *gm)
+{
+
+ FILE *f = fopen("toto.pos","w");
+ fprintf(f,"View \"\"{\n");
+
+ double Obj = 0.0;
+
+ for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
+ if ((*it)->geomType() == GEntity::Line)continue;
+ for (unsigned int i=0;i<(*it)->lines.size(); i++){
+ //Obj += MLineGEdgeDistance ( (*it)->lines[i] , *it,f );
+ Obj = std::max(MLineGEdgeDistance ( (*it)->lines[i] , *it, f ),Obj);
+ }
+ }
+
+ printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
+
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
+ if ((*it)->geomType() == GEntity::Plane)continue;
+ // printf("FACE %d with %d triangles\n",(*it)->tag(),(*it)->triangles.size());
+ for (unsigned int i=0;i<(*it)->triangles.size(); i++){
+ //Obj += MFaceGFaceDistance ( (*it)->triangles[i] , *it );
+ Obj = std::max(Obj,MFaceGFaceDistance ( (*it)->triangles[i] , *it ));
+ }
+ }
+
+ printf("DISTANCE TO GEOMETRY : 1D AND 2D PART %22.15E\n",Obj);
+ fprintf(f,"};\n");
+ fclose(f);
+ return Obj;
+}
+
+
diff --git a/contrib/HighOrderMeshOptimizer/OptHomCADDist.h b/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
new file mode 100644
index 0000000000000000000000000000000000000000..465335331580b57c2b921c46cdcb6f30d2689031
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
@@ -0,0 +1,23 @@
+// TODO: Header
+
+#ifndef OPTHOMCADDIST_H_
+#define OPTHOMCADDIST_H_
+
+
+class GModel;
+class MElement;
+class MLine;
+class MTriangle;
+class GEdge;
+class GFace;
+
+double MFaceGFaceDistance(MTriangle *t, GFace *gf,
+ std::vector<std::vector<SVector3> > *gsfT=0,
+ std::map<MVertex*,SVector3> *normalsToCAD=0);
+double MLineGEdgeDistance (MLine *l, GEdge *ge, FILE *f = 0);
+
+void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*, double> &distances);
+double distanceToGeometry(GModel *gm);
+
+
+#endif /* OPTHOMCADDIST_H_ */
diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h
new file mode 100644
index 0000000000000000000000000000000000000000..944cac46bbb253c870d49a365e58fcacbb9f4e47
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h
@@ -0,0 +1,96 @@
+// TODO: Copyright
+
+#ifndef _OPTHOMOBJCONTRIBCADDISTOLD_H_
+#define _OPTHOMOBJCONTRIBCADDISTOLD_H_
+
+#include "MeshOptObjContrib.h"
+
+
+template<class FuncType>
+class ObjContribCADDistOld : public ObjContrib, public FuncType
+{
+public:
+ ObjContribCADDistOld(double weight, double geomTol);
+ virtual ~ObjContribCADDistOld() {}
+ virtual ObjContrib *copy() const;
+ virtual void initialize(Patch *mesh);
+ virtual bool fail() { return false; }
+ virtual bool addContrib(double &Obj, alglib::real_1d_array &gradObj);
+ virtual void updateParameters() { FuncType::updateParameters(_min, _max); }
+ virtual bool targetReached() { return FuncType::targetReached(_min, _max); }
+ virtual bool stagnated() { return FuncType::stagnated(_min, _max); }
+ virtual void updateMinMax();
+
+protected:
+ Patch *_mesh;
+ double _weight;
+ double _geomTol;
+};
+
+
+template<class FuncType>
+ObjContribCADDistOld<FuncType>::ObjContribCADDistOld(double weight, double geomTol) :
+ ObjContrib("CADDist", FuncType::getNamePrefix()+"CADDist"),
+ _mesh(0), _weight(weight), _geomTol(geomTol)
+{
+}
+
+
+template<class FuncType>
+ObjContrib *ObjContribCADDistOld<FuncType>::copy() const
+{
+ return new ObjContribCADDistOld<FuncType>(*this);
+}
+
+
+template<class FuncType>
+void ObjContribCADDistOld<FuncType>::initialize(Patch *mesh)
+{
+ _mesh = mesh;
+
+ updateMinMax();
+ FuncType::initialize(_min, _max);
+}
+
+
+template<class FuncType>
+bool ObjContribCADDistOld<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
+{
+ _min = BIGVAL;
+ _max = -BIGVAL;
+
+ std::vector<double> gradF;
+ for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ 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;
+ }
+ }
+
+ return true;
+}
+
+
+template<class FuncType>
+void ObjContribCADDistOld<FuncType>::updateMinMax()
+{
+ _min = BIGVAL;
+ _max = -BIGVAL;
+
+ std::vector<double> dumGradF;
+ for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
+ double f;
+ if (_mesh->bndDistAndGradients(iEl, f, dumGradF, _geomTol)) {
+ _min = std::min(_min, f);
+ _max = std::max(_max, f);
+ }
+ }
+}
+
+
+#endif /* _OPTHOMOBJCONTRIBCADDISTOLD_H_ */