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

High order optimizer: moved functions for computation of distance to CAD to separate files

parent d91001df
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 407 additions and 513 deletions
contrib/HighOrderMeshOptimizer/CMakeLists.txt
+ 1
− 0
View file @ 759b7618
......@@ -8,6 +8,7 @@ set(SRC
OptHOM.cpp
OptHomRun.cpp
OptHomIntegralBoundaryDist.cpp
OptHomCADDist.cpp
ParamCoord.cpp
SuperEl.cpp
OptHomElastic.cpp
......
contrib/HighOrderMeshOptimizer/OptHOM.cpp
+ 119
− 502
View file @ 759b7618
......@@ -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
contrib/HighOrderMeshOptimizer/OptHOM.h
+ 0
− 11
View file @ 759b7618
......@@ -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)
......
contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp 0 → 100644
+ 168
− 0
View file @ 759b7618
// 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;
}
contrib/HighOrderMeshOptimizer/OptHomCADDist.h 0 → 100644
+ 23
− 0
View file @ 759b7618
// 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_ */
contrib/HighOrderMeshOptimizer/OptHomObjContribCADDistOld.h 0 → 100644
+ 96
− 0
View file @ 759b7618
// 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_ */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment