Skip to content
Snippets Groups Projects
Commit 15e842b9 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

fix compile

parent 1ebc4baf
No related branches found
No related tags found
No related merge requests found
......@@ -17,9 +17,12 @@
#include "Numeric.h"
#include "GaussLegendre1D.h"
#include "Context.h"
#if defined(HAVE_MESH)
#include "meshGFaceOptimize.h"
#include "meshGFaceLloyd.h"
#include "BackgroundMesh.h"
#endif
#if defined(HAVE_BFGS)
#include "stdafx.h"
......@@ -28,20 +31,6 @@
#define SQU(a) ((a)*(a))
class data_wrapper{
private:
const GFace* gf;
SPoint3 point;
public:
data_wrapper() {gf = NULL; point = SPoint3();}
~data_wrapper() {}
const GFace* get_face() {return gf;}
void set_face(const GFace* face) {gf = face;}
SPoint3 get_point() {return point;}
void set_point(SPoint3 _point) {point = SPoint3(_point);}
};
GFace::GFace(GModel *model, int tag)
: GEntity(model, tag), r1(0), r2(0), compound(0), va_geom_triangles(0)
{
......@@ -589,11 +578,10 @@ void GFace::getMeanPlaneData(double plan[3][3]) const
plan[i][j] = meanPlane.plan[i][j];
}
void GFace::computeMeshSizeFieldAccuracy(double &avg,double &max_e, double &min_e,
int &nE, int &GS)
{
#if defined(HAVE_MESH)
std::set<MEdge, Less_Edge> es;
for(unsigned int i = 0; i < getNumMeshElements(); i++){
MElement *e = getMeshElement(i);
......@@ -621,13 +609,13 @@ void GFace::computeMeshSizeFieldAccuracy(double &avg,double &max_e, double &min_
double correctLC = 0.5*(l1+l2);
double lone = it->length()/correctLC;
if (lone > oneoversqr2 && lone < sqr2) GS++;
double add = lone >1 ? (1. / lone) - 1. : lone - 1.;
//double add = lone > 1 ? (1. / lone) - 1. : lone - 1.;
avg += lone >1 ? (1. / lone) - 1. : lone - 1.;
max_e = std::max(max_e, lone);
min_e = std::min(min_e, lone);
}
avg = 100*exp(1./nE*avg);
#endif
}
double GFace::curvatureDiv(const SPoint2 &param) const
......@@ -893,8 +881,23 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p, bool onSurface) const
}
#if defined(HAVE_BFGS)
class data_wrapper{
private:
const GFace* gf;
SPoint3 point;
public:
data_wrapper() { gf = NULL; point = SPoint3(); }
~data_wrapper() {}
const GFace* get_face() { return gf; }
void set_face(const GFace* face) { gf = face; }
SPoint3 get_point() { return point; }
void set_point(SPoint3 _point) { point = SPoint3(_point); }
};
// Callback function for BFGS
void bfgs_callback(const alglib::real_1d_array& x, double& func, alglib::real_1d_array& grad, void* ptr) {
void bfgs_callback(const alglib::real_1d_array& x, double& func, alglib::real_1d_array& grad, void* ptr)
{
data_wrapper* w = static_cast<data_wrapper*>(ptr);
SPoint3 p = w->get_point();
const GFace* gf = w->get_face();
......@@ -909,19 +912,19 @@ void bfgs_callback(const alglib::real_1d_array& x, double& func, alglib::real_1d
// Value of the gradient
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(x[0], x[1]));
grad[0] = -(p.x() - pnt.x()) * der.left().x() - (p.y() - pnt.y()) * der.left().y() - (p.z() - pnt.z()) * der.left().z();
grad[1] = -(p.x() - pnt.x()) * der.right().x() - (p.y() - pnt.y()) * der.right().y() - (p.z() - pnt.z()) * der.right().z();
// printf("func %22.15E Gradients %22.15E %22.15E der %g %g %g\n",func, grad[0], grad[1],der.left().x(),der.left().y(),der.left().z());
grad[0] = -(p.x() - pnt.x()) * der.left().x() - (p.y() - pnt.y()) * der.left().y()
- (p.z() - pnt.z()) * der.left().z();
grad[1] = -(p.x() - pnt.x()) * der.right().x() - (p.y() - pnt.y()) * der.right().y()
- (p.z() - pnt.z()) * der.right().z();
// printf("func %22.15E Gradients %22.15E %22.15E der %g %g %g\n",
// func, grad[0], grad[1],der.left().x(),der.left().y(),der.left().z());
}
#endif
GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const
{
#if defined(HAVE_BFGS)
//
// Creating the optimisation problem
//
// printf("STARTING OPTIMIZATION\n");
alglib::ae_int_t dim = 2;
......@@ -948,7 +951,9 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
double dist = queryPoint.distance(spnt);
// fprintf(F,"SP(%g,%g,%g) {%g};\n",pnt.x(), pnt.y(), pnt.z(),dist);
// printf("lmocal (%12.5E %12.5E) (point) : %12.5E %12.5E %12.5E (query) : %12.5E %12.5E %12.5E DSIT %12.5E\n",u,v, pnt.x(), pnt.y(), pnt.z(),queryPoint.x(),queryPoint.y(),queryPoint.z(),
// printf("lmocal (%12.5E %12.5E) (point) : %12.5E %12.5E %12.5E (query) : "
// "%12.5E %12.5E %12.5E DSIT %12.5E\n",u,v, pnt.x(), pnt.y(), pnt.z(),
// queryPoint.x(),queryPoint.y(),queryPoint.z(),
// dist);
if (dist < min_dist) {
// printf("min_dist %f\n", dist);
......@@ -976,9 +981,7 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
minlbfgscreate(2, corr, x, state);
//
// Defining the stopping criterion
//
double epsg = 1.e-12;
double epsf = 0.0;
double epsx = 0.0;
......@@ -986,26 +989,22 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
minlbfgssetcond(state,epsg,epsf,epsx,maxits);
//
// Solving the problem
//
data_wrapper w;
w.set_point(queryPoint);
w.set_face(this);
minlbfgsoptimize(state,bfgs_callback,NULL,&w);
//
// Getting the results
//
alglib::minlbfgsreport rep;
minlbfgsresults(state,x,rep);
GPoint pntF = point(x[0], x[1]);
if (rep.terminationtype != 4){
// printf("Initial conditions (point) : %f %f %f local (%g %g) Looking for %f %f %f DIST = %12.5E at the end (%f %f) %f %f %f\n",
// printf("Initial conditions (point) : %f %f %f local (%g %g) "
// "Looking for %f %f %f DIST = %12.5E at the end (%f %f) %f %f %f\n",
// pnt.x(), pnt.y(), pnt.z(),min_u,min_v,
// queryPoint.x(),queryPoint.y(),queryPoint.z(),
// min_dist,x[0],x[1],pntF.x(), pntF.y(), pntF.z());
......@@ -1015,7 +1014,8 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
// ( queryPoint.z() - pntF.z()) * ( queryPoint.z() - pntF.z());
// if (sqrt(DDD) > 1.e-4)
/*
printf("Initial conditions (point) : %f %f %f local (%g %g) Looking for %f %f %f DIST = %12.5E at the end (%f %f) %f %f %f termination %d\n",
printf("Initial conditions (point) : %f %f %f local (%g %g) Looking for %f %f %f "
"DIST = %12.5E at the end (%f %f) %f %f %f termination %d\n",
pnt.x(), pnt.y(), pnt.z(),min_u,min_v,
queryPoint.x(),queryPoint.y(),queryPoint.z(),
min_dist,x[0],x[1],pntF.x(), pntF.y(), pntF.z(),rep.terminationtype);
......
......@@ -625,7 +625,7 @@ void Centerline::createSplitCompounds()
int num_gec = NE+i+1;
Msg::Info("Create Compound Line (%d) = %d discrete edge",
num_gec, pe->tag());
GEdge *gec = current->addCompoundEdge(e_compound,num_gec);
/* GEdge *gec = */ current->addCompoundEdge(e_compound,num_gec);
//gec->meshAttributes.Method = MESH_TRANSFINITE;
//gec->meshAttributes.nbPointsTransfinite = nbPoints;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment