diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a7b96fbf4f43eadb7ed43ed59d7d046308ec9d5d..40aa04d35e8dc94cc8432bbb8b6dfb9fb88772be 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -2088,7 +2088,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
getTriangle(par1, par2, <, U,V);
if(!lt && _mapping != RBF){
printf("ARRG POINT NOT FOUND--> should improve octree search \n");
- exit(1);
+ //exit(1);
//printf("POINT no success %d tris %d quad \n", triangles.size(), quadrangles.size());
GPoint gp = pointInRemeshedOctree(par1,par2);
gp.setNoSuccess();
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 5f4b4622a8a01e82acfac5f6db09dc92ec00d2da..6d82c1cd2c500ef56ed8520e72bca5f9796f5e9b 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -768,7 +768,7 @@ double gLevelsetPopcorn::operator() (const double x, const double y, const doubl
return val;
}
-gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag)
+gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag)
: gLevelsetPrimitive(tag)
{
std::vector<std::string> expressions(1, f);
@@ -789,6 +789,74 @@ double gLevelsetMathEval::operator() (const double x, const double y, const doub
return 1.;
}
+gLevelsetMathEvalAll::gLevelsetMathEvalAll(std::vector<std::string> expressions, int tag)
+ : gLevelsetPrimitive(tag)
+{
+ _hasDerivatives = true;
+ std::vector<std::string> variables(3);
+ variables[0] = "x";
+ variables[1] = "y";
+ variables[2] = "z";
+ _expr = new mathEvaluator(expressions, variables);
+}
+double gLevelsetMathEvalAll::operator() (const double x, const double y, const double z) const
+{
+ std::vector<double> values(3), res(13);
+ values[0] = x;
+ values[1] = y;
+ values[2] = z;
+ if(_expr->eval(values, res)) return res[0];
+ return 1.;
+}
+std::vector<double> gLevelsetMathEvalAll::getValueAndGradients (const double x, const double y, const double z) const
+{
+ std::vector<double> values(3), res(13);
+ values[0] = x;
+ values[1] = y;
+ values[2] = z;
+ if(_expr->eval(values, res)) return res;
+ return res;
+}
+
+void gLevelsetMathEvalAll::gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const
+{
+
+ std::vector<double> values(3), res(13);
+ values[0] = x;
+ values[1] = y;
+ values[2] = z;
+ if(_expr->eval(values, res)){
+ dfdx = res[1];
+ dfdy = res[2];
+ dfdz = res[3];
+ }
+
+}
+
+void gLevelsetMathEvalAll::hessian (double x, double y, double z,
+ double & dfdxx, double & dfdxy, double & dfdxz,
+ double & dfdyx, double & dfdyy, double & dfdyz,
+ double & dfdzx, double & dfdzy, double & dfdzz) const
+{
+
+ std::vector<double> values(3), res(13);
+ values[0] = x;
+ values[1] = y;
+ values[2] = z;
+ if(_expr->eval(values, res)){
+ dfdxx = res[4];
+ dfdxy = res[5];
+ dfdxz = res[6];
+ dfdyx = res[7];
+ dfdyy = res[8];
+ dfdyz = res[9];
+ dfdzx = res[10];
+ dfdzy = res[11];
+ dfdzz = res[12];
+ }
+}
+
gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
: gLevelsetPrimitive(tag), _box(NULL)
{
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index d4cb0f2d33640eaa2b533a3895db0f0b2877256e..609af08581e5ad09f7c26eb952be673bb3e29ad6 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -293,6 +293,23 @@ public:
int type() const { return UNKNOWN; }
};
+class gLevelsetMathEvalAll: public gLevelsetPrimitive
+{
+ mathEvaluator *_expr;
+public:
+ gLevelsetMathEvalAll(std::vector<std::string> f, int tag=1);
+ ~gLevelsetMathEvalAll(){ if (_expr) delete _expr; }
+ double operator() (const double x, const double y, const double z) const;
+ std::vector<double> getValueAndGradients(const double x, const double y, const double z) const;
+ void gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const;
+ void hessian (double x, double y, double z,
+ double & dfdxx, double & dfdxy, double & dfdxz,
+ double & dfdyx, double & dfdyy, double & dfdyz,
+ double & dfdzx, double & dfdzy, double & dfdzz ) const;
+ int type() const { return UNKNOWN; }
+};
+
class gLevelsetSimpleFunction: public gLevelsetPrimitive
{
simpleFunction<double> *_f;
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index e3cd21babcb9a30f19c1d6af5be9d092f1ae87b5..6e4b3fbd115c4875804459750ad4d4b668835bbd 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -216,7 +216,7 @@ void meshGFaceBamg(GFace *gf){
}
Mesh2 *refinedBamgMesh = 0;
- int iterMax = 11;
+ int iterMax = 41;
for (int k= 0; k < iterMax; k++){
int nbVert = bamgMesh->nv;
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 77944afc148ad3e028d317523f84353c5a5783e5..02a59f6e5ef21726b89b69cfd05cfd2fc47ae215 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -222,6 +222,14 @@ void meshMetric::computeValues()
while (it != _adj.end()) {
std::vector<MElement*> lt = it->second;
MVertex *ver = it->first;
+ //check if function is analytical
+ if (_fct->hasDerivatives()) {
+ printf("YOUPIE WE HAVE DERIVATIVES \n");
+ //int metricNumber = setOfMetrics.size();
+ //setOfMetrics[metricNumber].RECOMPUTE_METRIC = true;
+ exit(1);
+ //COMPUTE_ON_THE_FLY = true;
+ }
vals[ver]= (*_fct)(ver->x(),ver->y(),ver->z());
it++;
}
@@ -808,6 +816,99 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
throw;
}
metr = SMetric3(1.e-22);
+
+ //test linh
+ //if (RECOMPUTE_METRIC)‘{
+ // do computation
+ // call the fgradient function and the hessioan function
+ // _fct->gradient(x,y,z,dfdx,dfdy,dfdz);
+ // _fct->hessian(x,y,z,...);
+ //}
+
+ // //test emi
+ // double signed_dist = sqrt(x*x+y*y)-0.5;
+ // double dist = fabs(signed_dist);
+ // SMetric3 hessian;
+ // hessian(0,0) = y*y/pow(x*x+y*y,1.5); hessian(0,1) = -x*y/pow(x*x+y*y,1.5); hessian(0,2) = 0.0;
+ // hessian(1,0) = -x*y/pow(x*x+y*y,1.5); hessian(1,1) = x*x/pow(x*x+y*y,1.5); hessian(1,2) = 0.0;
+ // hessian(2,0) = 0.0; hessian(2,1) = 0.0; hessian(2,2) = 0.0;
+
+ // SVector3 gVec(x/sqrt(x*x+y*y), y/sqrt(x*x+y*y), 0.0);
+ // const double metric_value_hmax = 1./(hmax*hmax); // Gradient vector
+ // const double gMag = gVec.norm(), invGMag = 1./gMag;
+ // SMetric3 metric;
+
+ // if (signed_dist < _E && signed_dist > _E_moins && gMag != 0.0){
+ // const double metric_value_hmin = 1./(hmin*hmin);
+ // const SVector3 nVec = invGMag*gVec; // Unit normal vector
+ // double lambda_n = 0.; // Eigenvalues of metric for normal & tangential directions
+ // if (_technique==meshMetric::EIGENDIRECTIONS_LINEARINTERP_H){
+ // const double h_dist = hmin + ((hmax-hmin)/_E)*dist;
+ // double beta = CTX::instance()->mesh.smoothRatio;
+ // double h_dista = std::min( (hmin+(dist*log(beta))), hmax);
+ // lambda_n = 1./(h_dista*h_dista);
+ // }
+ // else if(_technique==meshMetric::EIGENDIRECTIONS){
+ // const double maximum_distance = (signed_dist>0.) ? _E : fabs(_E_moins); // ... or linear interpolation between 1/h_min^2 and 1/h_max^2
+ // lambda_n = metric_value_hmin + ((metric_value_hmax-metric_value_hmin)/maximum_distance)*dist;
+ // }
+ // std::vector<SVector3> tVec; // Unit tangential vectors
+ // std::vector<double> kappa; // Curvatures
+ // if (_dim == 2) { // 2D curvature formula: cf. R. Goldman, "Curvature formulas for implicit curves and surfaces", Computer Aided Geometric Design 22 (2005), pp. 632–658
+ // kappa.resize(2);
+ // kappa[0] = fabs(-gVec(1)*(-gVec(1)*hessian(0,0)+gVec(0)*hessian(0,1))+
+ // gVec(0)*(-gVec(1)*hessian(1,0)+gVec(0)*hessian(1,1)))*pow(invGMag,3);
+ // kappa[1] = 1.;
+ // tVec.resize(2);
+ // tVec[0] = SVector3(-nVec(1),nVec(0),0.);
+ // tVec[1] = SVector3(0.,0.,1.);
+ // }
+ // else { // 3D curvature formula: cf. A.G. Belyaev, A.A. Pasko and T.L. Kunii, "Ridges and Ravines on Implicit Surfaces," CGI, pp.530-535, Computer Graphics International 1998 (CGI'98), 1998
+ // fullMatrix<double> ImGG(3,3);
+ // ImGG(0,0) = 1.-gVec(0)*gVec(0); ImGG(0,1) = -gVec(0)*gVec(1); ImGG(0,2) = -gVec(0)*gVec(2);
+ // ImGG(1,0) = -gVec(1)*gVec(0); ImGG(1,1) = 1.-gVec(1)*gVec(1); ImGG(1,2) = -gVec(1)*gVec(2);
+ // ImGG(2,0) = -gVec(2)*gVec(0); ImGG(2,1) = -gVec(2)*gVec(1); ImGG(2,2) = 1.-gVec(2)*gVec(2);
+ // fullMatrix<double> hess(3,3);
+ // hessian.getMat(hess);
+ // fullMatrix<double> gN(3,3); // Gradient of unit normal
+ // gN.gemm(ImGG,hess,1.,0.);
+ // gN.scale(invGMag);
+ // fullMatrix<double> eigVecL(3,3), eigVecR(3,3);
+ // fullVector<double> eigValRe(3), eigValIm(3);
+ // gN.eig(eigValRe,eigValIm,eigVecL,eigVecR,false); // Eigendecomp. of gradient of unit normal
+ // kappa.resize(3); // Store abs. val. of eigenvalues (= principal curvatures only in non-normal directions)
+ // kappa[0] = fabs(eigValRe(0));
+ // kappa[1] = fabs(eigValRe(1));
+ // kappa[2] = fabs(eigValRe(2));
+ // tVec.resize(3); // Store normalized eigenvectors (= principal directions only in non-normal directions)
+ // tVec[0] = SVector3(eigVecR(0,0),eigVecR(1,0),eigVecR(2,0));
+ // tVec[0].normalize();
+ // tVec[1] = SVector3(eigVecR(0,1),eigVecR(1,1),eigVecR(2,1));
+ // tVec[1].normalize();
+ // tVec[2] = SVector3(eigVecR(0,2),eigVecR(1,2),eigVecR(2,2));
+ // tVec[2].normalize();
+ // std::vector<double> tVecDotNVec(3); // Store dot products with normal vector to look for normal direction
+ // tVecDotNVec[0] = fabs(dot(tVec[0],nVec));
+ // tVecDotNVec[1] = fabs(dot(tVec[1],nVec));
+ // tVecDotNVec[2] = fabs(dot(tVec[2],nVec));
+ // const int i_N = max_element(tVecDotNVec.begin(),tVecDotNVec.end())-tVecDotNVec.begin(); // Index of normal dir. = max. dot products (close to 0. in tangential dir.)
+ // kappa.erase(kappa.begin()+i_N); // Remove normal dir.
+ // tVec.erase(tVec.begin()+i_N);
+ // }
+ // const double invh_t0 = (_Np*kappa[0])/6.283185307, invh_t1 = (_Np*kappa[1])/6.283185307; // Inverse of tangential size 0
+ // const double lambda_t0 = invh_t0*invh_t0, lambda_t1 = invh_t1*invh_t1;
+ // const double lambdaP_n = std::min(std::max(lambda_n,metric_value_hmax),metric_value_hmin); // Truncate eigenvalues
+ // const double lambdaP_t0 = std::min(std::max(lambda_t0,metric_value_hmax),metric_value_hmin);
+ // const double lambdaP_t1 = (_dim == 2) ? 1. : std::min(std::max(lambda_t1,metric_value_hmax),metric_value_hmin);
+ // metric = SMetric3(lambdaP_n,lambdaP_t0,lambdaP_t1,nVec,tVec[0],tVec[1]);
+ // }
+ // else{// isotropic metric !
+ // SMetric3 mymetric(metric_value_hmax);
+ // metric = mymetric;
+ // }
+ // metr = metric;
+ // //end test emi
+
SPoint3 xyz(x,y,z), uvw;
double initialTol = MElement::getTolerance();
MElement::setTolerance(1.e-4);
@@ -825,8 +926,6 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
SMetric3 m4 = _nodalMetrics[e->getVertex(3)];
metr = interpolation(m1,m2,m3,m4,uvw[0],uvw[1],uvw[2]);
}
- //if recomputeAnalyticalValues
- //then compute exact metric
}
else{
@@ -840,6 +939,8 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
}
}
}
+
+
}
double meshMetric::getLaplacian (MVertex *v) {
diff --git a/Numeric/mathEvaluator.cpp b/Numeric/mathEvaluator.cpp
index 56c767333a73571b1650b9b9069f695fda3858eb..500e93fe0c470e6670b6a044efb47f8fa2045547 100644
--- a/Numeric/mathEvaluator.cpp
+++ b/Numeric/mathEvaluator.cpp
@@ -60,7 +60,16 @@ bool mathEvaluator::eval(std::vector<double> &values, std::vector<double> &res)
}
catch(smlib::mathex::error e) {
Msg::Error(e.what());
- return false;
+ double eps = 1.e-20;
+ for(unsigned int j = 0; j < values.size(); j++)
+ _variables[j] = values[j] + eps;
+ try {
+ res[i] = _expressions[i]->eval();
+ }
+ catch(smlib::mathex::error e2) {
+ Msg::Error(e2.what());
+ return false;
+ }
}
}
return true;
diff --git a/Numeric/simpleFunction.h b/Numeric/simpleFunction.h
index 931fa46fd0632cb14bc5e779ed332025ba7d8ce6..61754fc94a41166c05335ea8bace28ffc6d09697 100644
--- a/Numeric/simpleFunction.h
+++ b/Numeric/simpleFunction.h
@@ -13,13 +13,22 @@ template <class scalar>
class simpleFunction {
protected:
scalar _val;
+ bool _hasDerivatives;
public :
- simpleFunction(scalar val=0.0) : _val(val) {}
+ simpleFunction(scalar val=0.0) : _val(val), _hasDerivatives(false){}
virtual ~simpleFunction(){}
+ virtual bool hasDerivatives() {return _hasDerivatives;};
virtual scalar operator () (double x, double y, double z) const { return _val; }
virtual void gradient (double x, double y, double z,
scalar & dfdx, scalar & dfdy, scalar & dfdz) const
{ dfdx = dfdy = dfdz = 0.0; }
+ virtual void hessian (double x, double y, double z,
+ scalar & dfdxx, scalar & dfdxy, scalar & dfdxz,
+ scalar & dfdyx, scalar & dfdyy, scalar & dfdyz,
+ scalar & dfdzx, scalar & dfdzy, scalar & dfdzz ) const
+ { dfdxx = dfdxy = dfdxz = 0.0;
+ dfdyx = dfdyy = dfdyz = 0.0;
+ dfdzx = dfdzy = dfdzz = 0.0; }
};
template <class scalar>