diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index a6a0498312de1560caf1d11b31b2758cd59aa685..096ecef238410b86bfe9ab28ac798f6afa500a82 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -111,7 +111,7 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
if(ges.getSign() == 0){ // oops
Msg::Error("Something wrong in edge loop of size=%d, no sign !", wire.size());
for (std::list<GEdge* >::iterator it = wire.begin(); it != wire.end(); it++){
- printf("GEdge=%d begin=%g end =%d \n", (*it)->tag(), (*it)->getBeginVertex()->tag(), (*it)->getEndVertex()->tag() );
+ printf("GEdge=%d begin=%d end =%d \n", (*it)->tag(), (*it)->getBeginVertex()->tag(), (*it)->getEndVertex()->tag() );
}
break;
}
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index bd895b65bad6d027214d870f8524af710b771eae..c53b3f4ce755602087c7323b887ed1df726489a4 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -854,8 +854,11 @@ GPoint GFace::closestPoint(const SPoint3 & queryPoint, const double initialGuess
bool GFace::containsParam(const SPoint2 &pt) const
{
+
Range<double> uu = parBounds(0);
Range<double> vv = parBounds(1);
+ //printf("p =%g %g uu.low=%g uu.high=%g vv.low=%g vv.high=%g\n", pt.x(), pt.y(), uu.low(), uu.high(), vv.low(), vv.high());
+
if((pt.x() >= uu.low() && pt.x() <= uu.high()) &&
(pt.y() >= vv.low() && pt.y() <= vv.high()))
return true;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 12f4ea05c9f6b4e4237dd7af113ff44cc650e9b1..72a71cba2a2f4e753b7ce2855c07660fe8e9ef3a 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -30,6 +30,8 @@
#include "qualityMeasures.h"
#include "Field.h"
#include "OS.h"
+#include "Octree.h"
+#include "MElementOctree.h"
#include "HighOrder.h"
#include "meshGEdge.h"
#include "meshPartitionOptions.h"
@@ -400,42 +402,9 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
if (isMeshed) return true;
}
- // build face normal
- std::list<GEdge*>::iterator ite = edges.begin();
- std::list<int>::iterator itd = dir.begin();
- // SVector3 tan1, tan2;
- // for (int i=0; i<2; i++){
- // MVertex *v1 = (*ite)->lines[0]->getVertex(0);
- // MVertex *v2 = (*ite)->lines[0]->getVertex(1);
- // SPoint2 p1, p2;
- // if (*itd > 0.0){
- // //printf("v1=%g %g v2=%g %g \n", v1->x(), v1->y(), v2->x(), v2->y());
- // reparamMeshVertexOnFace(v1, gf, p1);
- // reparamMeshVertexOnFace(v2, gf, p2);
- // //printf("orient OK p1=%g %g p2=%g %g \n", p1.x(), p1.y(), p2.x(), p2.y());
- // }
- // else{
- // //printf("v1=%g %g v2=%g %g \n", v1->x(), v1->y(), v2->x(), v2->y());
- // reparamMeshVertexOnFace(v2, gf, p1);
- // reparamMeshVertexOnFace(v1, gf, p2);
- // //printf("orient -1 p1=%g %g p2=%g %g \n", p1.x(), p1.y(), p2.x(), p2.y());
- // }
- // //printf("edge=%d tan=%g %g \n", (*ite)->tag(), p2.x()-p1.x(),p2.y()-p1.y());
- // if (i==0) tan1 = SVector3(p2.x()-p1.x(),p2.y()-p1.y(),0.0);
- // else tan2 = SVector3(p2.x()-p1.x(),p2.y()-p1.y(),0.0);
- // ite++; itd++;
- // }
- // tan1.normalize(); tan2.normalize();
- // SVector3 next = crossprod(tan1,tan2);
- // next.normalize();
- // printf("next =%g %g %g \n", next.x(), next.y(), next.z());
- SVector3 next(0., 0., 1.0);
-
// build a set with all points of the boundaries
std::set<MVertex*> all_vertices;
- std::map<SPoint2, SVector3> pt2Normal;
- ite = edges.begin();
- itd = dir.begin();
+ std::list<GEdge*>::iterator ite = edges.begin();
while(ite != edges.end()){
if((*ite)->isSeam(gf)) return false;
if(!(*ite)->isMeshDegenerated()){
@@ -444,40 +413,11 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
all_vertices.insert(v1);
all_vertices.insert(v2);
-
- SPoint2 p1, p2;
- if (*itd > 0.0){
- reparamMeshVertexOnFace(v1, gf, p1);
- reparamMeshVertexOnFace(v2, gf, p2);
- }
- else{
- reparamMeshVertexOnFace(v2, gf, p1);
- reparamMeshVertexOnFace(v1, gf, p2);
- }
- SVector3 tan(p2.x()-p1.x(),p2.y()-p1.y(),0.0);
- tan.normalize();
- SVector3 ne = crossprod(tan, next);
- ne.normalize();
- std::map<SPoint2, SVector3>::iterator it = pt2Normal.find(p1);
- if (it == pt2Normal.end())
- pt2Normal.insert(std::make_pair(p1,ne));
- else{
- SVector3 n1 = it->second;
- it->second = n1+ne;
- }
- std::map<SPoint2, SVector3>::iterator it2 = pt2Normal.find(p2);
- if (it2 == pt2Normal.end())
- pt2Normal.insert(std::make_pair(p2,ne));
- else{
- SVector3 n1 = it2->second;
- it2->second = n1+ne;
- }
-
}
}
else
printf("edge %d degenerated mesh \n", (*ite)->tag());
- ++ite; ++itd;
+ ++ite;
}
std::list<GEdge*> emb_edges = gf->embeddedEdges();
@@ -550,9 +490,8 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
// use a divide & conquer type algorithm to create a triangulation.
// We add to the triangulation a box with 4 points that encloses the
// domain.
- std::map<SPoint2, SVector3> pt2NormalNEW;
+ DocRecord doc(points.size() + 4);
{
- DocRecord doc(points.size() + 4);
for(unsigned int i = 0; i < points.size(); i++){
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
(double)RAND_MAX;
@@ -563,14 +502,6 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
doc.points[i].data = points[i];
doc.points[i].adjacent = NULL;
- SPoint2 p(points[i]->u, points[i]->v);
- SPoint2 pNEW(points[i]->u+XX, points[i]->v+YY);
- std::map<SPoint2, SVector3>::iterator it = pt2Normal.find(p);
- if (it != pt2Normal.end() ){
- pt2NormalNEW.insert(std::make_pair(pNEW,it->second));
- }
- else
- Msg::Error("no point found \n");
}
// increase the size of the bounding box
@@ -626,32 +557,23 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
ite = edges.begin();
while(ite != edges.end()){
if(!(*ite)->isMeshDegenerated())
- recoverEdge
- (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
++ite;
}
ite = emb_edges.begin();
while(ite != emb_edges.end()){
if(!(*ite)->isMeshDegenerated())
- recoverEdge
- (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
++ite;
}
Msg::Debug("Recovering %d mesh Edges", edgesToRecover.size());
-
- if (Msg::GetVerbosity() == 10){
- doc.Voronoi();
- doc.makePosView("voronoi.pos", gf);
- doc.printMedialAxis(pt2NormalNEW, "skeleton.pos", gf);
- }
-
+
// effectively recover the medge
ite = edges.begin();
while(ite != edges.end()){
if(!(*ite)->isMeshDegenerated()){
- recoverEdge
- (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
}
++ite;
}
@@ -733,8 +655,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
ite = emb_edges.begin();
while(ite != emb_edges.end()){
if(!(*ite)->isMeshDegenerated())
- recoverEdge
- (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
++ite;
}
@@ -802,6 +723,38 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
outputScalarField(m->triangles, name, 1);
}
+ {
+ std::list<BDS_Face*>::iterator itt = m->triangles.begin();
+ while (itt != m->triangles.end()){
+ BDS_Face *t = *itt;
+ if(!t->deleted){
+ BDS_Point *n[4];
+ t->getNodes(n);
+ MVertex *v1 = recoverMap[n[0]];
+ MVertex *v2 = recoverMap[n[1]];
+ MVertex *v3 = recoverMap[n[2]];
+ if(!n[3]){
+ if(v1 != v2 && v1 != v3 && v2 != v3)
+ gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ }
+ else{
+ MVertex *v4 = recoverMap[n[3]];
+ gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
+ }
+ }
+ ++itt;
+ }
+ }
+ Octree * _octree = buildMElementOctree(gf->model());
+ if (Msg::GetVerbosity() == 10){
+ doc.Voronoi();
+ doc.makePosView("voronoi.pos", gf);
+ doc.printMedialAxis(_octree, "skeleton.pos", gf);
+ }
+ Octree_Delete(_octree);
+ gf->triangles.clear();
+ gf->quadrangles.clear();
+
int nb_swap;
//outputScalarField(m->triangles, "beforeswop.pos",1);
Msg::Debug("Delaunizing the initial mesh");
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 7cfb9069bb12c61682a9f45741b3aa43b6618e2c..856cc49d954fb7b68e0d29e40c1bf573546db91a 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -273,7 +273,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
char opts[128];
buildTetgenStructure(gr, in, numberedV);
if (Msg::GetVerbosity() == 10)
- sprintf(opts, "peVv");
+ sprintf(opts, "peVvS0s");
else
sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q': (Msg::GetVerbosity() > 6)? 'V': '\0');
try{
@@ -326,7 +326,8 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
gr->set(faces);
// now do insertion of points
- insertVerticesInRegion(gr);
+ //insertVerticesInRegion(gr);
+
#endif
}
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index b313a24d5afe742ca249aa9f6a4fe9a11977eeba..8fbccdf9ff2d4f31957520d6cfb8f586439a9dcc 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -576,7 +576,7 @@ void DocRecord::makePosView(std::string fileName, GFace *gf)
double pc[2] = {(double)points[i].where.h, (double)points[i].where.v};
if (!onHull(i)){
GPoint p0(pc[0], pc[1], 0.0);
- if (gf) p0 = gf->point(pc[0], pc[1]);
+ //if (gf) p0 = gf->point(pc[0], pc[1]);
fprintf(f,"SP(%g,%g,%g) {%g};\n",p0.x(),p0.y(),p0.z(),(double)i);
voronoiCell (i,pts);
for (unsigned int j = 0; j < pts.size(); j++){
@@ -584,10 +584,10 @@ void DocRecord::makePosView(std::string fileName, GFace *gf)
SPoint2 pp2 = pts[(j+1)%pts.size()];
GPoint p1(pp1.x(), pp1.y(), 0.0);
GPoint p2(pp2.x(), pp2.y(), 0.0);
- if (gf) {
- p1 = gf->point(p1.x(), p1.y());
- p2 = gf->point(p2.x(), p2.y());
- }
+ // if (gf) {
+ // p1 = gf->point(p1.x(), p1.y());
+ // p2 = gf->point(p2.x(), p2.y());
+ // }
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z(),
(double)i,(double)i);
@@ -602,8 +602,7 @@ void DocRecord::makePosView(std::string fileName, GFace *gf)
fclose(f);
}
-void DocRecord::printMedialAxis(std::map<SPoint2, SVector3> &pt2Normal,
- std::string fileName, GFace *gf)
+void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf)
{
FILE *f = fopen(fileName.c_str(),"w");
@@ -613,10 +612,6 @@ void DocRecord::printMedialAxis(std::map<SPoint2, SVector3> &pt2Normal,
std::vector<SPoint2> pts;
SPoint2 pc((double)points[i].where.h, (double)points[i].where.v);
if (!onHull(i)){
- std::map<SPoint2, SVector3>::const_iterator it = pt2Normal.find(pc);
- if (it == pt2Normal.end()) printf("pt not found \n");
- SVector3 n = it->second;
- //fprintf(f,"VP(%g,%g,%g) {%g,%g,%g};\n",pc.x(),pc.y(), 0.0, n.x(), n.y(), n.z());
GPoint p0(pc[0], pc[1], 0.0);
if (gf) p0 = gf->point(pc[0], pc[1]);
fprintf(f,"SP(%g,%g,%g) {%g};\n",p0.x(),p0.y(),p0.z(),(double)i);
@@ -632,12 +627,16 @@ void DocRecord::printMedialAxis(std::map<SPoint2, SVector3> &pt2Normal,
p1 = gf->point(p1.x(), p1.y());
p2 = gf->point(p2.x(), p2.y());
}
- if (dot(v1,n) > 0.0 && dot(v2,n) > 0.0){
+ double P1[3] = {p1.x(), p1.y(), p1.z()};
+ double P2[3] = {p2.x(), p2.y(), p2.z()};
+ MElement *m1 = (MElement*)Octree_Search(P1, _octree);
+ MElement *m2 = (MElement*)Octree_Search(P2, _octree);
+ if (m1 && m2){
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
p1.x(),p1.y(),p1.z(),
p2.x(),p2.y(),p2.z(),
(double)i,(double)i);
- }
+ }
}
}
}
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index 87eeee10fab48a1d5ae1295a3641d6f259fd6dd0..843b2a0e22ef424e2265c0e49bb889121d623f32 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -10,6 +10,7 @@
#include "fullMatrix.h"
#include "SPoint2.h"
#include "simpleFunction.h"
+#include "Octree.h"
class binding;
class GFace;
@@ -95,7 +96,7 @@ class DocRecord{
int hullSize() {return _hullSize;}
int onHull(PointNumero i) { return std::binary_search(_hull, _hull+_hullSize, i); }
void makePosView(std::string, GFace *gf=NULL);
- void printMedialAxis(std::map<SPoint2, SVector3> &pt2Normal, std::string, GFace *gf=NULL);
+ void printMedialAxis(Octree *_octree, std::string, GFace *gf=NULL);
double Lloyd (int);
void voronoiCell (PointNumero pt, std::vector<SPoint2> &pts) const;
static void registerBindings(binding *b);
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index d38bd96e3781fb526266271136958fc9da8445ad..73d61795855e0cf268c6e69adeefb14671f872e1 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -901,46 +901,258 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
}
}
-void signedDistancesPointsLine (std::vector<double>&distances,
- std::vector<SPoint3>&closePts,
- const std::vector<SPoint3> &pts,
- const SPoint3 &p1,
- const SPoint3 &p2){
-
+void signedDistancePointLine(const SPoint3 &p1,const SPoint3 &p2,const SPoint3 &p, double &d, SPoint3 &closePt){
SVector3 t1 = p2 - p1;
const double n2t1 = dot(t1, t1);
+ SVector3 pp1 = p - p1;
+ const double t12 = dot(pp1, t1) / n2t1;
+ d = 1.e10;
+ bool found = false;
+ if (t12 >= 0 && t12 <= 1.){
+ d = std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
+ closePt = p1 + (p2 - p1) * t12;
+ found = true;
+ }
+ if (p.distance(p1) < fabs(d)){
+ closePt = p1;
+ d = std::min(fabs(d), p.distance(p1));
+ }
+ if (p.distance(p2) < fabs(d)){
+ closePt = p2;
+ d = std::min(fabs(d), p.distance(p2));
+ }
+}
+
+void signedDistancesPointsLine (std::vector<double>&distances,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2){
- const unsigned pts_size=pts.size();
distances.clear();
- distances.resize(pts_size);
+ distances.resize(pts.size());
closePts.clear();
- closePts.resize(pts_size);
+ closePts.resize(pts.size());
double d;
- for (unsigned int i = 0; i <pts_size;++i){
- const SPoint3 &p = pts[i];
- SVector3 pp1 = p - p1;
- const double t12 = dot(pp1, t1) / n2t1;
- d = 1.e10;
- bool found = false;
+ for (unsigned int i = 0; i < pts.size();i++){
SPoint3 closePt;
- if (t12 >= 0 && t12 <= 1.){
- d = std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
- closePt = p1 + (p2 - p1) * t12;
- found = true;
- }
+ const SPoint3 &p = pts[i];
+ signedDistancePointLine(p1,p2,p,d,closePt);
+ distances[i] = d;
+ closePts[i] = closePt;
+ }
+}
+
+void changeReferential(const int direction,const SPoint3 &p,const SPoint3 &closePt,const SPoint3 &p1,const SPoint3 &p2,double* xp,double* yp,double* otherp,double* x,double* y,double* other){
+ if (direction==1){
+ const SPoint3 &d1=SPoint3(1.0,0.0,0.0);
+ const SPoint3 &d=SPoint3(p2.x()-p1.x(),p2.y()-p1.y(),p2.z()-p1.z());
+ double norm=sqrt( d.x()*d.x()+d.y()*d.y()+d.z()*d.z() );
+ const SPoint3 &dn=SPoint3(d.x()/norm,d.y()/norm,d.z()/norm);
+ const SPoint3 &d3=SPoint3(d1.y()*dn.z()-d1.z()*dn.y(),d1.z()*dn.x()-d1.x()*dn.z(),d1.x()*dn.y()-d1.y()*dn.x());
+ norm=sqrt( d3.x()*d3.x()+d3.y()*d3.y()+d3.z()*d3.z() );
+ const SPoint3 &d3n=SPoint3(d3.x()/norm,d3.y()/norm,d3.z()/norm);
+ const SPoint3 &d2=SPoint3(d3n.y()*d1.z()-d3n.z()*d1.y(),d3n.z()*d1.x()-d3n.x()*d1.z(),d3n.x()*d1.y()-d3n.y()*d1.x());
+ norm=sqrt( d2.x()*d2.x()+d2.y()*d2.y()+d2.z()*d2.z() );
+ const SPoint3 &d2n=SPoint3(d2.x()/norm,d2.y()/norm,d2.z()/norm);
+ *xp=p.x()*d1.x()+p.y()*d1.y()+p.z()*d1.z();
+ *yp=p.x()*d3n.x()+p.y()*d3n.y()+p.z()*d3n.z();
+ *otherp=p.x()*d2n.x()+p.y()*d2n.y()+p.z()*d2n.z();
+ *x=closePt.x()*d1.x()+closePt.y()*d1.y()+closePt.z()*d1.z();
+ *y=closePt.x()*d3n.x()+closePt.y()*d3n.y()+closePt.z()*d3n.z();
+ *other=closePt.x()*d2n.x()+closePt.y()*d2n.y()+closePt.z()*d2n.z();
+ }else{
+ const SPoint3 &d2=SPoint3(0.0,1.0,0.0);
+ const SPoint3 &d=SPoint3(p2.x()-p1.x(),p2.y()-p1.y(),p2.z()-p1.z());
+ double norm=sqrt( d.x()*d.x()+d.y()*d.y()+d.z()*d.z() );
+ const SPoint3 &dn=SPoint3(d.x()/norm,d.y()/norm,d.z()/norm);
+ const SPoint3 &d3=SPoint3(dn.y()*d2.z()-dn.z()*d2.y(),dn.z()*d2.x()-dn.x()*d2.z(),dn.x()*d2.y()-dn.y()*d2.x());
+ norm=sqrt( d3.x()*d3.x()+d3.y()*d3.y()+d3.z()*d3.z() );
+ const SPoint3 &d3n=SPoint3(d3.x()/norm,d3.y()/norm,d3.z()/norm);
+ const SPoint3 &d1=SPoint3(d2.y()*d3n.z()-d2.z()*d3n.y(),d2.z()*d3n.x()-d2.x()*d3n.z(),d2.x()*d3n.y()-d2.y()*d3n.x());
+ norm=sqrt( d1.x()*d1.x()+d1.y()*d1.y()+d1.z()*d1.z() );
+ const SPoint3 &d1n=SPoint3(d1.x()/norm,d1.y()/norm,d1.z()/norm);
+ *xp=p.x()*d2.x()+p.y()*d2.y()+p.z()*d2.z();
+ *yp=p.x()*d3n.x()+p.y()*d3n.y()+p.z()*d3n.z();
+ *otherp=p.x()*d1n.x()+p.y()*d1n.y()+p.z()*d1n.z();
+ *x=closePt.x()*d2.x()+closePt.y()*d2.y()+closePt.z()*d2.z();
+ *y=closePt.x()*d3n.x()+closePt.y()*d3n.y()+closePt.z()*d3n.z();
+ *other=closePt.x()*d1n.x()+closePt.y()*d1n.y()+closePt.z()*d1n.z();
+ }
+}
- if (p.distance(p1) < fabs(d)){
- closePt = p1;
- d = std::min(fabs(d), p.distance(p1));
+int computeDistanceRatio(const double &y, const double &yp,const double &x,const double &xp, double *distance, const double &r1, const double &r2){
+ double b;
+ double a;
+ if (y==yp){
+ b=-y;
+ a=0.0;
+ }else{
+ if (x==xp){
+ b=-x;
+ a==0.0;
+ }else{
+ b=(xp*y-x*yp)/(yp-y);
+ if (yp==0.0){
+ a=-(b+x)/y;
+ }else{
+ a=-(b+xp)/yp;
+ }
+ }
+ }
+ double ae;
+ double be;
+ double ce;
+ double da=r1*r1;
+ double db=r2*r2;
+ if (y==yp){
+ ae=1.0/da;
+ be=-(2*x)/da;
+ ce=(x*x/da)-1.0;
+ }else{
+ if (x==xp){
+ ae = 1.0/db;
+ be = -(2.0*y)/db;
+ ce = (y*y/db)-1.0;
+ }else{
+ if (fabs(a)<0.00001){
+ ae = 1.0/db;
+ be = -(2.0*y)/db;
+ ce = (y*y/db)-1.0;
+ }else{
+ double a2=a*a;
+ ae=(1.0/da)+(1.0/(db*a2));
+ be=(2.0*y)/(db*a)+(2.0*b)/(a2*db)-((2.0*x)/da);
+ ce=(x*x)/da+(b*b)/(db*a2)+(2.0*b*y)/(a*db)+(y*y/db)-1.0;
+ }
+ }
+ }
+ double rho=be*be-4*ae*ce;
+ double x1,x2,y1,y2,propdist;
+ if (rho<0) {
+ return 1;
+ }else{
+ x1=-(be+sqrt(rho))/(2.0*ae);
+ x2=(-be+sqrt(rho))/(2.0*ae);
+ if (y==yp){
+ y1=-b;
+ y2=-b;
+ }else{
+ if (x==xp){
+ y1=x1;
+ y2=x2;
+ x1=-b;
+ x2=-b;
+ }else{
+ if (fabs(a)<0.00001){
+ y1=x1;
+ y2=x2;
+ x1=-b;
+ x2=-b;
+ }else{
+ y1=-(b+x1)/a;
+ y2=-(b+x2)/a;
+ }
+ }
}
- if (p.distance(p2) < fabs(d)){
- closePt = p2;
- d = std::min(fabs(d), p.distance(p2));
+ if (x1==x2){
+ propdist=(y1-y)/(yp-y);
+ if(propdist<0.0){
+ propdist=(y2-y)/(yp-y);
+ }
+ }else{
+ if (xp!=x){
+ propdist=(x1-x)/(xp-x);
+ if (propdist<0.0){
+ propdist=(x2-x)/(xp-x);
+ }
+ }else{
+ if (yp!=y){
+ propdist=(y1-y)/(yp-y);
+ if(propdist<0.0){
+ propdist=(y2-y)/(yp-y);
+ }
+ }else{
+ propdist=0.01;
+ }
+ }
}
+ *distance=propdist;
+ return 0;
+ }
+}
+
+void signedDistancesPointsEllipseLine (std::vector<double>&distances,
+ std::vector<double> &distancesE,
+ std::vector<int>&isInYarn,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2){
+
+ distances.clear();
+ distances.resize(pts.size());
+ distancesE.clear();
+ distancesE.resize(pts.size());
+ isInYarn.clear();
+ isInYarn.resize(pts.size());
+ closePts.clear();
+ closePts.resize(pts.size());
+ double d;
+ for (unsigned int i = 0; i < pts.size();i++){
+ SPoint3 closePt;
+ const SPoint3 &p = pts[i];
+ signedDistancePointLine(p1,p2,p,d,closePt);
distances[i] = d;
closePts[i] = closePt;
+ int direction;
+ double distancesMin;
+ if (!(p.x()==closePt.x() && p.y()==closePt.y() && p.z()==closePt.z())){
+ double xp,yp,x,y,otherp,other,propdist;
+ if (p1.x()==p2.x()){
+ direction=1;
+ if (fabs(closePt.x()-0.0)<0.00000001) isInYarn[i]=1;
+ if (fabs(closePt.x()-2.2)<0.00000001) isInYarn[i]=4;
+ if (fabs(closePt.x()-4.4)<0.00000001) isInYarn[i]=2;
+ if (fabs(closePt.x()-6.6)<0.00000001) isInYarn[i]=5;
+ if (fabs(closePt.x()-8.8)<0.00000001) isInYarn[i]=3;
+ if (fabs(closePt.x()-11.0)<0.00000001) isInYarn[i]=1;
+ }else{
+ if (p1.y()==p2.y()){
+ direction=2;
+ if (fabs(closePt.y()-0.0)<0.00000001) isInYarn[i]=6;
+ if (fabs(closePt.y()-2.2)<0.00000001) isInYarn[i]=7;
+ if (fabs(closePt.y()-4.4)<0.00000001) isInYarn[i]=8;
+ if (fabs(closePt.y()-6.6)<0.00000001) isInYarn[i]=9;
+ if (fabs(closePt.y()-8.8)<0.00000001) isInYarn[i]=10;
+ if (fabs(closePt.y()-11.0)<0.00000001) isInYarn[i]=6;
+ }else{
+ printf("WTF %lf %lf\n",closePt.x(),closePt.y());
+ }
+ }
+ changeReferential(direction,p,closePt,p1,p2,&xp,&yp,&otherp,&x,&y,&other);
+ int result;
+ if (fabs(other-otherp)>0.01){
+ result=1;
+ }else{
+ result=computeDistanceRatio(y,yp,x,xp,&propdist,1.1,0.0875);
+ }
+ if (result==1){
+ distancesE[i] = 1.e10;
+ isInYarn[i]=0;
+ }else{
+ if (propdist<1.0){
+ isInYarn[i]=0;
+ distancesE[i]=(1.0/propdist)-1.0;
+ }else{
+ distancesE[i]=(1.0-(1.0/propdist))/3.0;
+ }
+ }
+ }else{
+ isInYarn[i]=0;
+ distancesE[i]=1000000.0;
+ }
}
-
}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 7a5d38018eb5eea4bade98efc6f8274f4937c5c4..8748b0da5ef716d715fbaac5751059924d3c058a 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -90,9 +90,21 @@ void signedDistancesPointsTriangle(std::vector<double> &distances,
const SPoint3 &p1,
const SPoint3 &p2,
const SPoint3 &p3);
-void signedDistancesPointsLine(std::vector<double> &distances,
- std::vector<SPoint3> &closePts,
- const std::vector<SPoint3> &pts,
- const SPoint3 &p1,
- const SPoint3 &p2);
+void signedDistancePointLine(const SPoint3 &p1,const SPoint3 &p2,const SPoint3 &p, double &distance, SPoint3 &closePt);
+void signedDistancesPointsLine (std::vector<double>&distances,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2);
+
+void changeReferential(const int direction,const SPoint3 &p,const SPoint3 &closePt,const SPoint3 &p1,const SPoint3 &p2,double* xp,double* yp,double* otherp,double* x,double* y,double* other);
+int computeDistanceRatio(const double &y, const double &yp,const double &x,const double &xp, double *distance, const double &r1, const double &r2);
+
+void signedDistancesPointsEllipseLine (std::vector<double>&distances,
+ std::vector<double>&distancesE,
+ std::vector<int>&isInYarn,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2);
#endif
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 9bc823c5787b9162696bef129a1b6ba46b3e8a43..e7a84a7732c8b3f54382283f66bf67581a924e61 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -24,12 +24,12 @@
StringXNumber DistanceOptions_Number[] = {
- {GMSH_FULLRC, "Point", NULL, 0.},
- {GMSH_FULLRC, "Line", NULL, 0.},
- {GMSH_FULLRC, "Surface", NULL, 0.},
+ {GMSH_FULLRC, "PhysPoint", NULL, 0.},
+ {GMSH_FULLRC, "PhysLine", NULL, 0.},
+ {GMSH_FULLRC, "PhysSurface", NULL, 0.},
{GMSH_FULLRC, "Computation", NULL, -1},
- {GMSH_FULLRC, "Min Scale", NULL, -1},
- {GMSH_FULLRC, "Max Scale", NULL, -1},
+ {GMSH_FULLRC, "MinScale", NULL, -1},
+ {GMSH_FULLRC, "MaxScale", NULL, -1},
{GMSH_FULLRC, "Orthogonal", NULL, -1}
};
@@ -171,10 +171,23 @@ PView *GMSH_DistancePlugin::execute(PView *v)
std::vector<GEntity*> _entities;
GModel::current()->getEntities(_entities);
+ GEntity* ge = _entities[_entities.size()-1];
+ int integrationPointTetra[2];
+ if (type==-100){
+ integrationPointTetra[0]=1;
+ integrationPointTetra[1]=4;
+ }else{
+ integrationPointTetra[0]=0;
+ integrationPointTetra[1]=0;
+ }
- int totNumNodes = 0;
- for (unsigned int i = 0; i < _entities.size() ; i++)
- totNumNodes += _entities[i]->mesh_vertices.size();
+ int numnodes = 0;
+ for(unsigned int i = 0; i < _entities.size()-1; i++)
+ numnodes += _entities[i]->mesh_vertices.size();
+ int totNodes=numnodes + _entities[_entities.size()-1]->mesh_vertices.size();
+ printf("%d\n",totNodes);
+ int order=ge->getMeshElement(0)->getPolynomialOrder();
+ int totNumNodes = totNodes+ge->getNumMeshElements()*integrationPointTetra[order-1];
if (totNumNodes ==0) {
Msg::Error("This plugin needs a mesh !");
@@ -192,7 +205,48 @@ PView *GMSH_DistancePlugin::execute(PView *v)
distances.reserve(totNumNodes);
pt2Vertex.reserve(totNumNodes);
- for (int i = 0; i < totNumNodes; i++) distances.push_back(1.e22);
+ std::map<MVertex*,double> _distanceE_map;
+ std::map<MVertex*,int> _isInYarn_map;
+ std::vector<int> index;
+ std::vector<double> distancesE;
+ std::vector<int> isInYarn;
+ std::vector<SPoint3> closePts;
+ std::vector<double> distances2;
+ std::vector<double> distancesE2;
+ std::vector<int> isInYarn2;
+ std::vector<SPoint3> closePts2;
+
+ if (type==-100){
+ index.clear();
+ distancesE.clear();
+ closePts.clear();
+ isInYarn.clear();
+ isInYarn.reserve(totNumNodes);
+ closePts.reserve(totNumNodes);
+ distancesE.reserve(totNumNodes);
+ index.reserve(totNumNodes);
+ distances2.clear();
+ distancesE2.clear();
+ closePts2.clear();
+ isInYarn2.clear();
+ distances2.reserve(totNumNodes);
+ isInYarn2.reserve(totNumNodes);
+ closePts2.reserve(totNumNodes);
+ distancesE2.reserve(totNumNodes);
+ }
+
+ for (int i = 0; i < totNumNodes; i++) {
+ distances.push_back(1.e22);
+ if (type==-100){
+ distancesE.push_back(1.e22);
+ isInYarn.push_back(0);
+ closePts.push_back(SPoint3(0.,0.,0.));
+ distances2.push_back(1.e22);
+ distancesE2.push_back(1.e22);
+ isInYarn2.push_back(0);
+ closePts2.push_back(SPoint3(0.,0.,0.));
+ }
+ }
int k = 0;
for(unsigned int i = 0; i < _entities.size(); i++){
@@ -202,11 +256,43 @@ PView *GMSH_DistancePlugin::execute(PView *v)
MVertex *v = ge->mesh_vertices[j];
pts.push_back(SPoint3(v->x(), v->y(),v->z()));
_distance_map.insert(std::make_pair(v, 0.0));
+ if (type==-100){
+ index.push_back(v->getIndex());
+ _isInYarn_map.insert(std::make_pair(v, 0));
+ _distanceE_map.insert(std::make_pair(v, 0.0));
+ }
pt2Vertex[k] = v;
k++;
}
}
-
+
+ if (type==-100){
+ double jac[3][3];
+ for (unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement(i);
+ IntPt *ptsi;
+ int nptsi;
+ double uvw[4][3];
+ e->getIntegrationPoints(e->getPolynomialOrder(),&nptsi, &ptsi);
+ for(int j = 0; j < 4; j++) {
+ double xyz[3] = {e->getVertex(j)->x(),
+ e->getVertex(j)->y(),
+ e->getVertex(j)->z()};
+ e->xyz2uvw(xyz, uvw[j]);
+ }
+
+ for(int ip = 0; ip < nptsi; ip++){
+ const double u = ptsi[ip].pt[0];
+ const double v = ptsi[ip].pt[1];
+ const double w = ptsi[ip].pt[2];
+ const double weight = ptsi[ip].weight;
+ const double detJ = e->getJacobian(u, v, w, jac);
+ SPoint3 p;
+ e->pnt(u, v, w, p);
+ pts.push_back(p);
+ }
+ }
+ }
// Compute geometrical distance to mesh boundaries
//------------------------------------------------------
if (type < 0.0 ){
@@ -230,37 +316,155 @@ PView *GMSH_DistancePlugin::execute(PView *v)
for(unsigned int k = 0; k < g2->getNumMeshElements(); k++){
std::vector<double> iDistances;
std::vector<SPoint3> iClosePts;
+ std::vector<double> iDistancesE;
+ std::vector<int> iIsInYarn;
MElement *e = g2->getMeshElement(k);
MVertex *v1 = e->getVertex(0);
MVertex *v2 = e->getVertex(1);
SPoint3 p1(v1->x(), v1->y(), v1->z());
SPoint3 p2(v2->x(), v2->y(), v2->z());
- if(e->getNumVertices() == 2){
- signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+ if((e->getNumVertices() == 2 and order==1) or (e->getNumVertices() == 3 and order==2)){
+ if (type==-100){
+// if ( !((p1.x()==p2.x()) & (p1.y()==p2.y()) & (p1.z()==p2.z())) ){
+ signedDistancesPointsEllipseLine(iDistances, iDistancesE, iIsInYarn, iClosePts, pts, p1,p2);
+// }
+ }else{
+ signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+ }
}
- else if(e->getNumVertices() == 3){
+ else if(e->getNumVertices() == 3 and order==1){
MVertex *v3 = e->getVertex(2);
SPoint3 p3 (v3->x(),v3->y(),v3->z());
signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
}
for (unsigned int kk = 0; kk< pts.size(); kk++) {
- if (std::abs(iDistances[kk]) < distances[kk]){
- distances[kk] = std::abs(iDistances[kk]);
- MVertex *v = pt2Vertex[kk];
- _distance_map[v] = distances[kk];
+ if (type==-100){
+ if( !((p1.x()==p2.x()) & (p1.y()==p2.y()) & (p1.z()==p2.z()))){
+// printf("%d %d %d %lf %lf %lf %lf\n",kk,isInYarn[kk],iIsInYarn[kk],distances[kk],iDistances[kk],distancesE[kk],iDistancesE[kk]);
+ if (iIsInYarn[kk]>0){
+ if (isInYarn[kk]==0){
+ distances[kk] = iDistances[kk];
+ distancesE[kk]= iDistancesE[kk];
+ isInYarn[kk] = iIsInYarn[kk];
+ closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),iClosePts[kk].z());
+ }else{
+ if (isInYarn[kk]!=iIsInYarn[kk]){
+ if (isInYarn2[kk]==0){
+ distances2[kk] = iDistances[kk];
+ distancesE2[kk]= iDistancesE[kk];
+ isInYarn2[kk] = iIsInYarn[kk];
+ closePts2[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),iClosePts[kk].z());
+ }else{
+ if (isInYarn2[kk]==iIsInYarn[kk]){
+ if (iDistancesE[kk] < distancesE2[kk]){
+ distances2[kk] = iDistances[kk];
+ distancesE2[kk]= iDistancesE[kk];
+ closePts2[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),iClosePts[kk].z());
+ }
+ }
+ }
+ }else{
+ if (iDistancesE[kk] < distancesE[kk]){
+ distances[kk] = iDistances[kk];
+ distancesE[kk]= iDistancesE[kk];
+ closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),iClosePts[kk].z());
+ }
+ }
+ }
+ }else{
+ if (isInYarn[kk]==0){
+ if (iDistancesE[kk] < distancesE[kk]){
+ distances[kk] = iDistances[kk];
+ distancesE[kk]= iDistancesE[kk];
+ closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),iClosePts[kk].z());
+ }
+ }
+ }
+ }
+ }else{
+ if (std::abs(iDistances[kk]) < distances[kk]){
+ distances[kk] = std::abs(iDistances[kk]);
+ MVertex *v = pt2Vertex[kk];
+ _distance_map[v] = distances[kk];
+ }
}
}
}
}
}
+ if (type==-100){
+ for (unsigned int kk = 0; kk< pts.size(); kk++) {
+ if (isInYarn2[kk]>0){
+ if (distancesE2[kk]>distancesE[kk]){
+ distances[kk]=distances2[kk];
+ distancesE[kk]=distancesE2[kk];
+ isInYarn[kk]=isInYarn2[kk];
+ closePts[kk]=closePts2[kk];
+ }
+ }
+ if (kk<totNodes){
+ MVertex *v = pt2Vertex[kk];
+ _distance_map[v] = distancesE[kk];
+ _distanceE_map[v] = distances[kk];
+ _isInYarn_map[v] = isInYarn[kk];
+ }
+ }
+ }
if (!existEntity){
if (id_pt != 0) Msg::Error("The Physical Point does not exist !");
if (id_line != 0) Msg::Error("The Physical Line does not exist !");
if (id_face != 0) Msg::Error("The Physical Surface does not exist !");
return view;
}
+ printView(_entities, _distance_map);
+ if (type==-100){
+
+ Msg::Info("Writing integrationPointInYarn.pos");
+ FILE* f5 = fopen("integrationPointInYarn.pos","w");
+ FILE* f6 = fopen("integrationPointInYarn.bin","wb");
+ FILE* f7 = fopen("integrationPointInYarn.txt","w");
+ int j=0;
+ fprintf(f5,"View \"integrationPointInYarn\"{\n");
+ for (std::vector<int>::iterator it = isInYarn.begin(); it !=isInYarn.end(); it++) {
+ if (j>=totNodes){
+ int iPIY=*it;
+ fwrite(&iPIY,sizeof(int),1,f6);
+ fprintf(f7,"%d %lf %lf %lf\n",iPIY,pts[j].x(), pts[j].y(), pts[j].z());
+ fprintf(f5,"SP(%g,%g,%g){%d};\n",
+ pts[j].x(), pts[j].y(), pts[j].z(),
+ *it);
+ }
+ j++;
+ }
+ fclose(f6);
+ fclose(f7);
+ fprintf(f5,"};\n");
+ fclose(f5);
- printView(_entities, _distance_map);
+ Msg::Info("Writing isInYarn.pos");
+ FILE * f4 = fopen("isInYarn.pos","w");
+ fprintf(f4,"View \"isInYarn\"{\n");
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement(i);
+ fprintf(f4,"SS(");
+ std::vector<int> inYarn;
+ for(int j = 0; j < e->getNumVertices(); j++) {
+ MVertex *v = e->getVertex(j);
+ if(j) fprintf(f4,",%g,%g,%g",v->x(),v->y(), v->z());
+ else fprintf(f4,"%g,%g,%g", v->x(),v->y(), v->z());
+ std::map<MVertex*, int>::iterator it = _isInYarn_map.find(v);
+ inYarn.push_back(it->second);
+ }
+ fprintf(f4,"){");
+ for (int i=0; i<inYarn.size(); i++){
+ if (i) fprintf(f4,",%d", inYarn[i]);
+ else fprintf(f4,"%d", inYarn[i]);
+ }
+ fprintf(f4,"};\n");
+ }
+ fprintf(f4,"};\n");
+ fclose(f4);
+ }
}
@@ -359,12 +563,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
//compute also orthogonal vector to distance field
// A Uortho = -C DIST
//------------------------------------------------
- if (ortho > 0 && _maxDim != 2){
- Msg::Error("The orthogonal field is only implemented for a 2D distance field !");
- }
- else if (ortho > 0 && _maxDim == 2){
+ if (ortho > 0){
#if defined(HAVE_SOLVER)
-
+
#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys2 = new linearSystemCSRTaucs<double>;
#else
@@ -477,4 +678,3 @@ PView *GMSH_DistancePlugin::execute(PView *v)
return view;
}
-