diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index f241e458cd25c51fe188f5ab89e1a8b8684af4e5..baaa04f621e6cd8ea7cc9eb18cdadff2a3fc10e8 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -5,7 +5,6 @@
set(SRC
boundaryLayersData.cpp
- distanceToMesh.cpp
intersectCurveSurface.cpp
GEntity.cpp STensor3.cpp
GVertex.cpp GEdge.cpp GFace.cpp GRegion.cpp
diff --git a/Geo/distanceToMesh.cpp b/Geo/distanceToMesh.cpp
deleted file mode 100644
index 195cdf1c7709bd23898283e1132f448cf11cd80e..0000000000000000000000000000000000000000
--- a/Geo/distanceToMesh.cpp
+++ /dev/null
@@ -1,96 +0,0 @@
-#include "distanceToMesh.h"
-#include "MElement.h"
-#include "Numeric.h"
-
-distanceToMesh::distanceToMesh(GModel *gm, std::string physical, int nbClose)
- : _computed (false) , _gm(gm), _nbClose(nbClose)
-{
- std::map<int, std::vector<GEntity*> > groups [4];
- gm->getPhysicalGroups(groups);
- for (GModel::piter it = gm->firstPhysicalName() ;
- it != gm->lastPhysicalName() ; ++it){
- if (it->second == physical){
- _entities = groups[it->first.first][it->first.second];
- }
- }
-}
-
-
-#if defined(HAVE_ANN)
-distanceToMesh::~distanceToMesh()
-{
- if (_computed){
- delete [] _index;
- delete [] _dist;
- delete _kdtree;
- annDeallocPts (_nodes);
- }
-}
-
-void distanceToMesh::setup()
-{
- _computed = true;
- std::set<MVertex *> _all;
- for (unsigned int i=0;i<_entities.size();i++){
- for (unsigned int k = 0; k < _entities[i]->getNumMeshElements(); k++) {
- for (int j = 0; j< _entities[i]->getMeshElement(k)->getNumVertices();j++){
- _all.insert(_entities[i]->getMeshElement(k)->getVertex(j));
- _v2e.insert (std::make_pair(_entities[i]->getMeshElement(k)->getVertex(j),
- _entities[i]->getMeshElement(k)));
- }
- }
- }
- _nodes = annAllocPts(_all.size(), 3);
- std::set<MVertex*>::iterator itp = _all.begin();
- int ind = 0;
- while (itp != _all.end()){
- MVertex* pt = *itp;
- _nodes[ind][0] = pt->x();
- _nodes[ind][1] = pt->y();
- _nodes[ind][2] = pt->z();
- _vertices.push_back(pt);
- itp++; ind++;
- }
- _kdtree = new ANNkd_tree(_nodes, _all.size(), 3);
- _index = new ANNidx[_nbClose];
- _dist = new ANNdist[_nbClose];
-}
-
-double distanceToMesh::operator () (double x, double y, double z)
-{
- if (!_computed) setup();
- double point[3] = {x,y,z};
- _kdtree->annkSearch(point, _nbClose, _index, _dist);
- std::set<MElement*> elements;
- for (int i=0;i<_nbClose;i++){
- int iVertex = _index [i];
- MVertex *v = _vertices[iVertex];
- for (std::multimap<MVertex*,MElement*>::iterator itm =
- _v2e.lower_bound(v); itm != _v2e.upper_bound(v); ++itm)
- elements.insert (itm->second);
- }
- double minDistance = 1.e22;
- SPoint3 closest;
- for (std::set<MElement*>::iterator it = elements.begin();
- it != elements.end();++it){
- double distance;
- MVertex *v1 = (*it)->getVertex(0);
- MVertex *v2 = (*it)->getVertex(1);
- SPoint3 p1(v1->x(), v1->y(), v1->z());
- SPoint3 p2(v2->x(), v2->y(), v2->z());
- if ((*it)->getDim() == 1){
- SPoint3 pt;
- signedDistancePointLine(p1, p2,SPoint3(x,y,z),distance,pt);
- if (fabs(distance) < fabs(minDistance)){
- minDistance=distance;
- }
- }
- }
- return minDistance;
-}
-#else
-double distanceToMesh::operator () (double x, double y, double z)
-{
- Msg::Fatal ("impossible to compute distance to a mesh without ANN");
-}
-#endif
diff --git a/Geo/distanceToMesh.h b/Geo/distanceToMesh.h
deleted file mode 100644
index e89be2b3bc3dae6b4343575d9774d017a88e2614..0000000000000000000000000000000000000000
--- a/Geo/distanceToMesh.h
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef _DISTANCE_TO_MESH_
-#define _DISTANCE_TO_MESH_
-#include <ANN/ANN.h>
-#include "GModel.h"
-#include "simpleFunction.h"
-class distanceToMesh : public simpleFunction <double>
-{
- bool _computed;
- GModel * _gm;
- const int _nbClose;
-#if defined(HAVE_ANN)
- std::vector<GEntity*> _entities;
- std::vector<MVertex*> _vertices;
- std::multimap<MVertex*,MElement*> _v2e;
- ANNkd_tree *_kdtree;
- ANNpointArray _nodes;
- ANNidxArray _index;
- ANNdistArray _dist;
- void setup ();
-#endif
-public :
- distanceToMesh(GModel *gm, std::string physical, int nbClose = 5);
- double operator () (double x, double y, double z);
-#if defined(HAVE_ANN)
- ~distanceToMesh();
-#endif
-};
-
-#endif
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index c8a24255f827c6488cabf579fb98fdf57d275e1d..33521a7305b1a25d697b1b8299b0607656c4e3f9 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -16,6 +16,7 @@
#include "MElement.h"
#include "Numeric.h"
#include "cartesian.h"
+#include <ANN/ANN.h>
void insertActiveCells(double x, double y, double z, double rmax,
cartesianBox<double> &box)
@@ -848,16 +849,6 @@ double gLevelsetMathEvalAll::operator() (const double x, const double y, const d
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
{
@@ -897,123 +888,108 @@ void gLevelsetMathEvalAll::hessian (double x, double y, double z,
}
}
-gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
- : gLevelsetPrimitive(tag), _box(NULL)
+gLevelsetDistMesh::gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose)
+ : _gm(gm), _nbClose(nbClose)
{
- modelBox = new GModel();
- modelBox->load(box+std::string(".msh"));
- modelGeom = new GModel();
- modelGeom->load(geom);
-
- //EMI FIXME THIS
- int levels = 3;
- // double rmax = 0.1;
- // double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
-
- //FILLING POINTS FROM GEOMBOX
- std::vector<SPoint3> points;
- std::vector<SPoint3> refinePoints;
- for(GModel::viter vit = modelBox->firstVertex(); vit != modelBox->lastVertex(); vit++){
- for(unsigned int k = 0; k < (*vit)->getNumMeshVertices(); k++){
- MVertex *v = (*vit)->getMeshVertex(k);
- SPoint3 p(v->x(), v->y(), v->z());
- points.push_back(p);
+ std::map<int, std::vector<GEntity*> > groups [4];
+ gm->getPhysicalGroups(groups);
+ for (GModel::piter it = gm->firstPhysicalName() ;
+ it != gm->lastPhysicalName() ; ++it){
+ if (it->second == physical){
+ _entities = groups[it->first.first][it->first.second];
}
}
- for (GModel::eiter eit = modelBox->firstEdge(); eit != modelBox->lastEdge(); eit++){
- for(unsigned int k = 0; k < (*eit)->getNumMeshVertices(); k++){
- MVertex *ve = (*eit)->getMeshVertex(k);
- SPoint3 pe(ve->x(), ve->y(), ve->z());
- points.push_back(pe);
- }
+ if (_entities.size() == 0){
+ Msg::Error("distanceToMesh: the physical name '%s' does not exist in the GModel", physical.c_str());
}
- //FILLING POINTS FROM STL
- for (GModel::fiter fit = modelGeom->firstFace(); fit != modelGeom->lastFace(); fit++){
- for(unsigned int k = 0; k < (*fit)->getNumMeshVertices(); k++){
- MVertex *vf = (*fit)->getMeshVertex(k);
- SPoint3 pf(vf->x(), vf->y(), vf->z());
- points.push_back(pf);
- }
- for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
- MElement *e = (*fit)->getMeshElement(k);
- if (e->getType() == TYPE_TRI){
- MVertex *v1 = e->getVertex(0);
- MVertex *v2 = e->getVertex(1);
- MVertex *v3 = e->getVertex(2);
- SPoint3 cg( (v1->x()+v2->x()+v3->x())/3.,
- (v1->y()+v2->y()+v3->y())/3.,
- (v1->z()+v2->z()+v3->z())/3.);
- refinePoints.push_back(cg);
+ //setup
+ std::set<MVertex *> _all;
+ for (unsigned int i=0;i<_entities.size();i++){
+ for (unsigned int k = 0; k < _entities[i]->getNumMeshElements(); k++) {
+ MElement *e = _entities[i]->getMeshElement(k);
+ for (int j = 0; j< e->getNumVertices();j++){
+ MVertex *v = _entities[i]->getMeshElement(k)->getVertex(j);
+ _all.insert(v);
+ _v2e.insert (std::make_pair(v,e));
}
}
}
- //FOR CAD
- //for (GModel::fiter fit = modelGeom->firstFace(); fit != modelGeom->lastFace(); fit++)
- // (*fit)->fillPointCloud(sampling, &points);
-
- if (points.size() == 0) {Msg::Fatal("No points on surfaces \n"); };
-
- SBoundingBox3d bb;
- for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
- for(unsigned int i = 0; i < refinePoints.size(); i++) bb += refinePoints[i];
- //bb.scale(1.01, 1.01, 1.01);
- SVector3 range = bb.max() - bb.min();
- double minLength = std::min( range.x(), std::min(range.y(), range.z()));
- double hmin = minLength / 5;
- int NX = range.x() / hmin;
- int NY = range.y() / hmin;
- int NZ = range.z() / hmin;
- if(NX < 2) NX = 2;
- if(NY < 2) NY = 2;
- if(NZ < 2) NZ = 2;
- double rtube = 2.*hmin; //std::max(lx, std::max(ly, lz))*2.;
-
- Msg::Info(" bounding box min: %g %g %g -- max: %g %g %g",
- bb.min().x(), bb.min().y(), bb.min().z(),
- bb.max().x(), bb.max().y(), bb.max().z());
- Msg::Info(" Nx=%d Ny=%d Nz=%d", NX, NY, NZ);
-
- _box = new cartesianBox<double>(bb.min().x(), bb.min().y(), bb.min().z(),
- SVector3(range.x(), 0, 0),
- SVector3(0, range.y(), 0),
- SVector3(0, 0, range.z()),
- NX, NY, NZ, levels);
- for (int i = 0; i < NX; i++)
- for (int j = 0; j < NY; j++)
- for (int k = 0; k < NZ; k++)
- _box->insertActiveCell(_box->getCellIndex(i, j, k));
-
- cartesianBox<double> *parent = _box, *child;
- while((child = parent->getChildBox())){
- //Msg::Info(" level %d ", child->getLevel());
- for(unsigned int i = 0; i < refinePoints.size(); i++)
- insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
- rtube / pow(2., (levels - child->getLevel())), *child);
- parent = child;
+ _nodes = annAllocPts(_all.size(), 3);
+ std::set<MVertex*>::iterator itp = _all.begin();
+ int ind = 0;
+ while (itp != _all.end()){
+ MVertex* pt = *itp;
+ _nodes[ind][0] = pt->x();
+ _nodes[ind][1] = pt->y();
+ _nodes[ind][2] = pt->z();
+ _vertices.push_back(pt);
+ itp++; ind++;
}
+ _kdtree = new ANNkd_tree(_nodes, _all.size(), 3);
+ _index = new ANNidx[_nbClose];
+ _dist = new ANNdist[_nbClose];
- //Msg::Info("Removing cells to match mesh topology constraints");
- removeBadChildCells(_box);
- removeParentCellsWithChildren(_box);
-
- //Msg::Info("Initializing nodal values in the cartesian grid");
- _box->createNodalValues();
- //Msg::Info("Computing levelset on the cartesian grid");
- computeLevelset(modelGeom, *_box);
+}
- //Msg::Info("Renumbering mesh vertices across levels");
- _box->renumberNodes();
+#if defined(HAVE_ANN)
+gLevelsetDistMesh::~gLevelsetDistMesh()
+{
+ delete [] _index;
+ delete [] _dist;
+ if (_kdtree) delete _kdtree;
+ if (_nodes) annDeallocPts (_nodes);
- _box->writeMSH("yeah.msh", false);
}
-double gLevelsetDistGeom::operator() (const double x, const double y, const double z) const
+double gLevelsetDistMesh::operator () (const double x, const double y, const double z) const
+{
+ double point[3] = {x,y,z};
+ _kdtree->annkSearch(point, _nbClose, _index, _dist);
+ std::set<MElement*> elements;
+ for (int i=0;i<_nbClose;i++){
+ int iVertex = _index [i];
+ MVertex *v = _vertices[iVertex];
+ for (std::multimap<MVertex*,MElement*>::const_iterator itm =
+ _v2e.lower_bound(v); itm != _v2e.upper_bound(v); ++itm)
+ elements.insert (itm->second);
+ }
+ double minDistance = 1.e22;
+ SPoint3 closest;
+ for (std::set<MElement*>::iterator it = elements.begin();
+ it != elements.end();++it){
+ double distance;
+ MVertex *v1 = (*it)->getVertex(0);
+ MVertex *v2 = (*it)->getVertex(1);
+ SPoint3 p1(v1->x(), v1->y(), v1->z());
+ SPoint3 p2(v2->x(), v2->y(), v2->z());
+ SPoint3 pt;
+ if ((*it)->getDim() == 1){
+ signedDistancePointLine(p1, p2,SPoint3(x,y,z),distance,pt);
+ }
+ else if ((*it)->getDim() == 2){
+ MVertex *v3 = (*it)->getVertex(2);
+ SPoint3 p3(v3->x(), v3->y(), v3->z());
+ signedDistancePointTriangle(p1, p2, p3,SPoint3(x,y,z),distance,pt);
+ }
+ else if ((*it)->getDim() == 2){
+ Msg::Error("Cannot compute a dsitance to an entity of dim \n");
+ }
+ if (fabs(distance) < fabs(minDistance)){
+ minDistance=distance;
+ }
+ }
+ return -1.0*minDistance;
+}
+#else
+double gLevelsetDistMesh::operator () (double x, double y, double z)
{
- double dist = _box->getValueContainingPoint(x,y,z);
- return dist;
+ Msg::Fatal ("impossible to compute distance to a mesh without ANN");
}
+#endif
+
+
#if defined (HAVE_POST)
gLevelsetPostView::gLevelsetPostView(int index, int tag)
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index d6716c50b151afbb80c3813578312754f3d09b1b..d62a5e849735e134d1e7e74928d89881c2dc82a7 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -22,7 +22,9 @@
#include "cartesian.h"
#include "simpleFunction.h"
-
+#if defined(HAVE_ANN)
+#include "ANN/ANN.h"
+#endif
#if defined(HAVE_POST)
#include "PView.h"
#include "OctreePost.h"
@@ -306,7 +308,6 @@ 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,
@@ -331,21 +332,29 @@ public:
int type() const { return UNKNOWN; }
};
-class gLevelsetDistGeom: public gLevelsetPrimitive
+class gLevelsetDistMesh: public gLevelsetPrimitive
{
- cartesianBox<double> *_box;
- GModel *modelBox, *modelGeom;
-public:
- gLevelsetDistGeom(std::string g, std::string f, int tag=1);
- ~gLevelsetDistGeom(){
- delete modelBox;
- delete modelGeom;
- if (_box) delete _box;
- }
+ GModel * _gm;
+ const int _nbClose;
+#if defined(HAVE_ANN)
+ std::vector<GEntity*> _entities;
+ std::vector<MVertex*> _vertices;
+ std::multimap<MVertex*,MElement*> _v2e;
+ ANNkd_tree *_kdtree;
+ ANNpointArray _nodes;
+ ANNidxArray _index;
+ ANNdistArray _dist;
+#endif
+public :
+ gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose = 5);
double operator () (const double x, const double y, const double z) const;
+#if defined(HAVE_ANN)
+ ~gLevelsetDistMesh();
+#endif
int type() const { return UNKNOWN; }
};
+
#if defined(HAVE_POST)
class gLevelsetPostView : public gLevelsetPrimitive
{
@@ -630,7 +639,6 @@ public:
gLevelsetNACA00(double x0, double y0, double c, double t);
~gLevelsetNACA00() {}
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,
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 4a41e97ec28dbc0ab702672f5458b978086df13e..ad9c500eba1c80eead023bfbd991daa0c91d9801 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -707,12 +707,12 @@ void meshMetric::scaleMetric( int nbElementsTarget,
SMetric3 m3 = nmt[e->getVertex(2)];
if (_dim == 2){
SMetric3 m = interpolation(m1,m2,m3,0.3333,0.3333);
- N += sqrt(m.determinant()) * e->getVolume() * 3.0;
+ N += sqrt(m.determinant()) * e->getVolume() * 4./sqrt(3.0); //3.0
}
else{
SMetric3 m4 = nmt[e->getVertex(3)];
SMetric3 m = interpolation(m1,m2,m3,m4,0.25,0.25,0.25);
- N += sqrt(m.determinant()) * e->getVolume() * 4.0;
+ N += sqrt(m.determinant()) * e->getVolume() * 12./sqrt(2.0); //4.0;
}
}
double scale = pow ((double)nbElementsTarget/N,2.0/_dim);
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 0dd9bc3f03d206fffe03c58cc50458a07deabe66..d89a17ee96ca43fc2e0e98a762b825600b4b7d86 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -667,15 +667,15 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
fullVector<double> f(N), feps(N), dx(N);
for (int iter = 0; iter < MAXIT; iter++){
- func(x, f, data);
+ func(x, f, data);
- bool isZero = false;
- for (int k=0; k<N; k++) {
- if (f(k) == 0. ) isZero = true;
- else isZero = false;
- if (isZero == false) break;
- }
- if (isZero) break;
+ bool isZero = false;
+ for (int k=0; k<N; k++) {
+ if (f(k) == 0. ) isZero = true;
+ else isZero = false;
+ if (isZero == false) break;
+ }
+ if (isZero) break;
for (int j = 0; j < N; j++){
double h = EPS * fabs(x(j));
@@ -707,9 +707,9 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
}
/*
-min_a f(x+a*d);
+ min_a f(x+a*d);
-f(x+a*d) = f(x) + f'(x) (
+ f(x+a*d) = f(x) + f'(x) (
*/
@@ -739,11 +739,11 @@ void gmshLineSearch(double (*func)(fullVector<double> &, void *), void* data,
if (temp > test) test = temp;
}
/*
- for (int j=0;j<100;j++){
+ for (int j=0;j<100;j++){
double sx = (double)j/99;
for (i=0;i<n;i++) x(i)=xold(i)+10*sx*p(i);
double jzede = (*func)(x,data);
- }
+ }
*/
alamin = TOLX / test;
@@ -752,7 +752,7 @@ void gmshLineSearch(double (*func)(fullVector<double> &, void *), void* data,
for (i = 0; i < n; i++) x(i) = xold(i) + alam*p(i);
f = (*func)(x, data);
// printf("f = %g x = %g %g alam = %g p = %g %g\n",f,x(0),x(1),alam,p(0),p(1));
- if (alam < alamin) {
+ if (alam < alamin) {
for (i = 0; i <n; i++) x(i) = xold(i);
// printf("ALERT : alam %g alamin %g\n",alam,alamin);
check = 1;
@@ -825,48 +825,101 @@ double minimize_grad_fd(double (*func)(fullVector<double> &, void *),
}
/*
-P(p) = p1 + t1 xi + t2 eta
+ P(p) = p1 + t1 xi + t2 eta
-t1 = (p2-p1) ; t2 = (p3-p1) ;
+ t1 = (p2-p1) ; t2 = (p3-p1) ;
-(P(p) - p) = d n
+ (P(p) - p) = d n
-(p1 + t1 xi + t2 eta - p) = d n
-t1 xi + t2 eta + d n = p - p1
+ (p1 + t1 xi + t2 eta - p) = d n
+ t1 xi + t2 eta + d n = p - p1
-| t1x t2x -nx | |xi | |px-p1x|
-| t1y t2y -ny | |eta | = |py-p1y|
-| t1z t2z -nz | |d | |pz-p1z|
+ | t1x t2x -nx | |xi | |px-p1x|
+ | t1y t2y -ny | |eta | = |py-p1y|
+ | t1z t2z -nz | |d | |pz-p1z|
-distance to segment
+ distance to segment
- P(p) = p1 + t (p2-p1)
+ P(p) = p1 + t (p2-p1)
- (p - P(p)) * (p2-p1) = 0
- (p - p1 - t (p2-p1) ) * (p2-p1) = 0
- - t ||p2-p1||^2 + (p-p1)(p2-p1) = 0
+ (p - P(p)) * (p2-p1) = 0
+ (p - p1 - t (p2-p1) ) * (p2-p1) = 0
+ - t ||p2-p1||^2 + (p-p1)(p2-p1) = 0
- t = (p-p1)*(p2-p1)/||p2-p1||^2
+ t = (p-p1)*(p2-p1)/||p2-p1||^2
*/
-void signedDistancesPointsTriangle(std::vector<double> &distances,
- std::vector<SPoint3> &closePts,
- const std::vector<SPoint3> &pts,
- const SPoint3 &p1,
- const SPoint3 &p2,
- const SPoint3 &p3)
+void signedDistancePointTriangle(const SPoint3 &p1,const SPoint3 &p2, const SPoint3 &p3,
+ const SPoint3 &p, double &d, SPoint3 &closePt)
+
{
+
SVector3 t1 = p2 - p1;
SVector3 t2 = p3 - p1;
SVector3 t3 = p3 - p2;
SVector3 n = crossprod(t1, t2);
n.normalize();
-
+ const double n2t1 = dot(t1, t1);
+ const double n2t2 = dot(t2, t2);
+ const double n2t3 = dot(t3, t3);
double mat[3][3] = {{t1.x(), t2.x(), -n.x()},
{t1.y(), t2.y(), -n.y()},
{t1.z(), t2.z(), -n.z()}};
double inv[3][3];
double det = inv3x3(mat, inv);
+ if(det == 0.0) return;
+
+ double u, v;
+ SVector3 pp1 = p - p1;
+ u = (inv[0][0] * pp1.x() + inv[0][1] * pp1.y() + inv[0][2] * pp1.z());
+ v = (inv[1][0] * pp1.x() + inv[1][1] * pp1.y() + inv[1][2] * pp1.z());
+ d = (inv[2][0] * pp1.x() + inv[2][1] * pp1.y() + inv[2][2] * pp1.z());
+ double sign = (d > 0) ? 1. : -1.;
+ if (d == 0) sign = 1.e10;
+ if (u >= 0 && v >= 0 && 1.-u-v >= 0.0){
+ closePt = SPoint3(0.,0.,0.);//TO DO
+ }
+ else {
+ const double t12 = dot(pp1, t1) / n2t1;
+ const double t13 = dot(pp1, t2) / n2t2;
+ SVector3 pp2 = p - p2;
+ const double t23 = dot(pp2, t3) / n2t3;
+ d = 1.e10;
+ SPoint3 closePt;
+ if (t12 >= 0 && t12 <= 1.){
+ d = sign * std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
+ closePt = p1 + (p2 - p1) * t12;
+ }
+ if (t13 >= 0 && t13 <= 1.){
+ if (p.distance(p1 + (p3 - p1) * t13) < fabs(d)) closePt = p1 + (p3 - p1) * t13;
+ d = sign * std::min(fabs(d), p.distance(p1 + (p3 - p1) * t13));
+ }
+ if (t23 >= 0 && t23 <= 1.){
+ if (p.distance(p2 + (p3 - p2) * t23) < fabs(d)) closePt = p2 + (p3 - p2) * t23;
+ d = sign * std::min(fabs(d), p.distance(p2 + (p3 - p2) * t23));
+ }
+ if (p.distance(p1) < fabs(d)){
+ closePt = p1;
+ d = sign * std::min(fabs(d), p.distance(p1));
+ }
+ if (p.distance(p2) < fabs(d)){
+ closePt = p2;
+ d = sign * std::min(fabs(d), p.distance(p2));
+ }
+ if (p.distance(p3) < fabs(d)){
+ closePt = p3;
+ d = sign * std::min(fabs(d), p.distance(p3));
+ }
+ }
+
+}
+void signedDistancesPointsTriangle(std::vector<double> &distances,
+ std::vector<SPoint3> &closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2,
+ const SPoint3 &p3)
+{
const unsigned pts_size = pts.size();
distances.clear();
distances.resize(pts_size);
@@ -876,15 +929,28 @@ void signedDistancesPointsTriangle(std::vector<double> &distances,
for (unsigned int i = 0; i < pts_size; ++i)
distances[i] = 1.e22;
- if(det == 0.0) return;
-
+ SVector3 t1 = p2 - p1;
+ SVector3 t2 = p3 - p1;
+ SVector3 t3 = p3 - p2;
+ SVector3 n = crossprod(t1, t2);
+ n.normalize();
const double n2t1 = dot(t1, t1);
const double n2t2 = dot(t2, t2);
const double n2t3 = dot(t3, t3);
+ double mat[3][3] = {{t1.x(), t2.x(), -n.x()},
+ {t1.y(), t2.y(), -n.y()},
+ {t1.z(), t2.z(), -n.z()}};
+ double inv[3][3];
+ double det = inv3x3(mat, inv);
+ if(det == 0.0) return;
- double u, v, d;
- for (unsigned int i = 0; i < pts_size; ++i){
+ for (unsigned int i = 0; i < pts.size(); i++) {
+ double d;
+ SPoint3 closePt;
const SPoint3 &p = pts[i];
+
+ //signedDistancePointTrianglePrecomputed(p1, p2, p3, p, d, closePt);
+ double u, v;
SVector3 pp1 = p - p1;
u = (inv[0][0] * pp1.x() + inv[0][1] * pp1.y() + inv[0][2] * pp1.z());
v = (inv[1][0] * pp1.x() + inv[1][1] * pp1.y() + inv[1][2] * pp1.z());
@@ -892,8 +958,7 @@ void signedDistancesPointsTriangle(std::vector<double> &distances,
double sign = (d > 0) ? 1. : -1.;
if (d == 0) sign = 1.e10;
if (u >= 0 && v >= 0 && 1.-u-v >= 0.0){
- distances[i] = d;
- closePts[i] = SPoint3(0.,0.,0.);//TO DO
+ closePt = SPoint3(0.,0.,0.);//TO DO
}
else {
const double t12 = dot(pp1, t1) / n2t1;
@@ -901,43 +966,41 @@ void signedDistancesPointsTriangle(std::vector<double> &distances,
SVector3 pp2 = p - p2;
const double t23 = dot(pp2, t3) / n2t3;
d = 1.e10;
- //bool found = false;
SPoint3 closePt;
if (t12 >= 0 && t12 <= 1.){
- d = sign * std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
- closePt = p1 + (p2 - p1) * t12;
- //found = true;
+ d = sign * std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
+ closePt = p1 + (p2 - p1) * t12;
}
if (t13 >= 0 && t13 <= 1.){
- if (p.distance(p1 + (p3 - p1) * t13) < fabs(d)) closePt = p1 + (p3 - p1) * t13;
- d = sign * std::min(fabs(d), p.distance(p1 + (p3 - p1) * t13));
- //found = true;
+ if (p.distance(p1 + (p3 - p1) * t13) < fabs(d)) closePt = p1 + (p3 - p1) * t13;
+ d = sign * std::min(fabs(d), p.distance(p1 + (p3 - p1) * t13));
}
if (t23 >= 0 && t23 <= 1.){
- if (p.distance(p2 + (p3 - p2) * t23) < fabs(d)) closePt = p2 + (p3 - p2) * t23;
- d = sign * std::min(fabs(d), p.distance(p2 + (p3 - p2) * t23));
- //found = true;
+ if (p.distance(p2 + (p3 - p2) * t23) < fabs(d)) closePt = p2 + (p3 - p2) * t23;
+ d = sign * std::min(fabs(d), p.distance(p2 + (p3 - p2) * t23));
}
if (p.distance(p1) < fabs(d)){
- closePt = p1;
- d = sign * std::min(fabs(d), p.distance(p1));
- }
+ closePt = p1;
+ d = sign * std::min(fabs(d), p.distance(p1));
+ }
if (p.distance(p2) < fabs(d)){
- closePt = p2;
- d = sign * std::min(fabs(d), p.distance(p2));
- }
+ closePt = p2;
+ d = sign * std::min(fabs(d), p.distance(p2));
+ }
if (p.distance(p3) < fabs(d)){
- closePt = p3;
- d = sign * std::min(fabs(d), p.distance(p3));
+ closePt = p3;
+ d = sign * std::min(fabs(d), p.distance(p3));
}
- //d = sign * std::min(fabs(d), std::min(std::min(p.distance(p1),
- // p.distance(p2)),p.distance(p3)));
- distances[i] = d;
- closePts[i] = closePt;
}
+ //end signedDistance
+
+ distances[i] = d;
+ closePts[i] = closePt;
}
}
+
+
void signedDistancePointLine(const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p,
double &d, SPoint3 &closePt)
{
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 678bd2aa60404803e944f0a031d9586b50bb113b..91e5cfc2f6343d25b2d946e4a280e1772c5a840f 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -110,12 +110,17 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
fullVector<double> &x, void *data, double relax=1., double tolx=1.e-6);
double minimize_grad_fd(double (*func)(fullVector<double> &, void *),
fullVector<double> &x, void *data);
+
+void signedDistancePointTriangle(const SPoint3 &p1,const SPoint3 &p2, const SPoint3 &p3,
+ const SPoint3 &p, double &d, SPoint3 &closePt);
+
void signedDistancesPointsTriangle(std::vector<double> &distances,
std::vector<SPoint3> &closePts,
const std::vector<SPoint3> &pts,
const SPoint3 &p1,
const SPoint3 &p2,
const SPoint3 &p3);
+
void signedDistancePointLine(const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p,
double &distance, SPoint3 &closePt);
void signedDistancesPointsLine (std::vector<double>&distances,
diff --git a/wrappers/gmshpy/gmshGeo.i b/wrappers/gmshpy/gmshGeo.i
index 512247447dedd8dd6f63bbe83ef94c8986727447..6390afe4355942a9e2431f7538357c49ceddb10e 100644
--- a/wrappers/gmshpy/gmshGeo.i
+++ b/wrappers/gmshpy/gmshGeo.i
@@ -37,7 +37,6 @@
#include "SBoundingBox3d.h"
#include "Curvature.h"
#include "simpleFunction.h"
- #include "distanceToMesh.h"
#include "GeomMeshMatcher.h"
%}
@@ -96,7 +95,6 @@ namespace std {
%include "SBoundingBox3d.h"
%include "Curvature.h"
%include "gmshLevelset.h"
-%include "distanceToMesh.h"
%include "GeomMeshMatcher.h"