From 6abca093aeb4cfa629252f9c863395f1680a80f5 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 19 Feb 2016 13:50:41 +0000
Subject: [PATCH] refactorization of reparametrization in progress
---
Geo/MVertex.cpp | 8 ++
Geo/discreteDiskFace.cpp | 130 ++++++++++++++++++++++++++--
Geo/discreteDiskFace.h | 3 +
Mesh/meshGFace.cpp | 2 +-
Mesh/meshGFaceDelaunayInsertion.cpp | 14 ++-
Numeric/Numeric.cpp | 2 +-
6 files changed, 150 insertions(+), 9 deletions(-)
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 874e84f8be..c6bf343a7b 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -362,6 +362,14 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> ¶ms
return;
}
+#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
+ if (gf->geomType() == GEntity::DiscreteDiskSurface ){
+ discreteDiskFace *gfc = (discreteDiskFace*) gf;
+ params.push_back(gfc->parFromVertex(v));
+ return ;
+ }
+#endif
+
if(v->onWhat()->dim() == 0){
GVertex *gv = (GVertex*)v->onWhat();
std::list<GEdge*> ed = gv->edges();
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 810db95119..20df1506e0 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -452,7 +452,13 @@ GPoint discreteDiskFace::point(const double par1, const double par2) const
SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
{
- if (v->onWhat()->dim()==2)Msg::Fatal("discreteDiskFace::parFromVertex should not be called with a vertex that is classified on a surface");
+ if (v->onWhat() == this){
+ double uu,vv;
+ v->getParameter(0,uu);
+ v->getParameter(1,vv);
+ return SPoint2(uu,vv);
+ }
+
std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y());
// The 1D mesh has been re-done
@@ -481,19 +487,20 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
SVector3 discreteDiskFace::normal(const SPoint2 ¶m) const
{
- return SVector3();
+ return GFace::normal(param);
}
double discreteDiskFace::curvatureMax(const SPoint2 ¶m) const
{
-
- return false;
+ throw;
+ return false;
}
double discreteDiskFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const
{
- return false;
+ throw;
+ return false;
}
Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
@@ -505,7 +512,10 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
MTriangle* tri = NULL;
if (ddft) tri = ddft->tri;
- else Msg::Error("discreteDiskFace::firstDer << triangle not found");
+ else {
+ Msg::Warning("discreteDiskFace::firstDer << triangle not found %g %g",param[0],param[1]);
+ return Pair<SVector3, SVector3>(SVector3(1,0,0), SVector3(0,1,0));
+ }
std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it0 =
firstDerivatives.find(tri->getVertex(0));
@@ -583,4 +593,112 @@ void discreteDiskFace::putOnView()
#endif
}
+// useful for mesh generators ----------------------------------------
+// Intersection of a circle and a plane
+GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
+ const SVector3 &p, const double &d,
+ double uv[2]) const
+{
+ SVector3 n = crossprod(n1,n2);
+ n.normalize();
+ for (int i=-1;i<(int)discrete_triangles.size();i++){
+ discreteDiskFaceTriangle *ct;
+ double U,V;
+ if (i == -1) getTriangleUV(uv[0],uv[1], &ct, U,V);
+ else ct = &_ddft[i];
+ if (!ct) continue;
+ // we have two planes, defined with n1,n2 and t1,t2
+ // we have then a direction m = (n1 x n2) x (t1 x t2)
+ // and a point p that defines a straight line
+ // the point is situated in the half plane defined
+ // by direction n1 and point p (exclude the one on the
+ // other side)
+
+ SVector3 t1 = ct->v2 - ct->v1;
+ SVector3 t2 = ct->v3 - ct->v1;
+
+ // let us first compute point q to be the intersection
+ // of the plane of the triangle with the line L = p + t n1
+ // Compute w = t1 x t2 = (a,b,c) and write a point of the plabe as
+ // X(x,y,z) = ax + by + cz - (v1_x a + v1_y b + v1_z * c) = 0
+ // Then
+ // a (p_x + t n1_x) + a (p_y + t n1_y) + c (p_z + t n1_z) -
+ // (v1_x a + v1_y b + v1_z * c) = 0
+ // gives t
+
+ SVector3 w = crossprod(t1,t2);
+ double t = d;
+ // if everything is not coplanar ...
+ if (fabs(dot(n1,w)) > 1.e-12){
+ t = dot(ct->v1-p,w)/dot(n1,w);
+ }
+ SVector3 q = p + n1*t;
+
+ // consider the line that intersects both planes in its
+ // parametric form
+ // X(x,y,z) = q + t m
+ // We have now to intersect that line with the sphere
+ // (x-p_x)^2 + (y-p_y)^2 + (z-p_z)^2 = d^2
+ // Inserting the parametric form of the line into the sphere gives
+ // (t m_x + q_x - p_x)^2 + (t m_y + q_y - p_y)^2 + (t m_z + q_z - p_z)^2 = d^2
+ // t^2 (m_x^2 + m_y^2 + m_z^2) + 2 t (m_x (q_x - p_x) + m_y (q_y - p_z) +
+ // m_z (q_z - p_z)) + ((q_x - p_x)^2 + (q_y - p_y)^2 + (q_z - p_z)^2 - d^2) = 0
+ // t^2 m^2 + 2 t (m . (q-p)) + ((q-p)^2 - d^2) = 0
+ // Let us call ta and tb the two roots
+ // they correspond to two points s_1 = q + ta m and s_2 = q + tb m
+ // we choose the one that corresponds to (s_i - p) . n1 > 0
+ SVector3 m = crossprod(w,n);
+
+
+ const double a = dot(m,m);
+ const double b = 2*dot(m,q-p);
+ const double c = dot(q-p,q-p) - d*d;
+ const double delta = b*b-4*a*c;
+
+ if (delta >= 0){
+ const double ta = (-b + sqrt(delta)) / (2.*a);
+ const double tb = (-b - sqrt(delta)) / (2.*a);
+ SVector3 s1 = q + m * ta;
+ SVector3 s2 = q + m * tb;
+ SVector3 s;
+ if (dot(s1-p,n1) > 0){
+ s = s1;
+ }
+ else if (dot(s2-p,n1) > 0){
+ s = s2;
+ }
+ else continue;
+
+ // we have now to look if the point is inside the triangle
+ // s = v1 + u t1 + v t2
+ // we know the system has a solution because s is in the plane
+ // defined by v1 and v2 we solve then
+ // (s - v1) . t1 = u t1^2 + v t2 . t1
+ // (s - v1) . t2 = u t1 . t2 + v t2^2
+
+ double mat[2][2], b[2],uv[2];
+ mat[0][0] = dot(t1,t1);
+ mat[1][1] = dot(t2,t2);
+ mat[0][1] = mat[1][0] = dot(t1,t2);
+ b[0] = dot(s-ct->v1,t1) ;
+ b[1] = dot(s-ct->v1,t2) ;
+ sys2x2(mat,b,uv);
+ // check now if the point is inside the triangle
+ if (uv[0] >= -1.e-6 && uv[1] >= -1.e-6 &&
+ 1.-uv[0]-uv[1] >= -1.e-6 ) {
+ SVector3 pp =
+ ct->p1 * ( 1.-uv[0]-uv[1] ) +
+ ct->p2 *uv[0] +
+ ct->p3 *uv[1] ;
+ return GPoint(s.x(),s.y(),s.z(),this,pp);
+ }
+ }
+ }
+ GPoint pp(0);
+ pp.setNoSuccess();
+ Msg::Debug("ARGG no success intersection circle");
+ return pp;
+}
+
+
#endif
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 4d553f8d60..6aa25db7a7 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -60,6 +60,9 @@ class discreteDiskFace : public GFace {
virtual void secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
GEntity::GeomType geomType() const { return DiscreteDiskSurface; }
+ GPoint intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
+ const SVector3 &p, const double &d,
+ double uv[2]) const;
protected:
//------------------------------------------------
// a copy of the mesh that should not be destroyed
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 0aceb8e58a..f2b605a949 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2331,7 +2331,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// BOUNDARY LAYER
modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf);
-
+
if(algoDelaunay2D(gf)){
if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL)
bowyerWatsonFrontal(gf, &equivalence, ¶metricCoordinates);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 85ca7063d1..df8d26e558 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -24,6 +24,7 @@
#include "Field.h"
#include "GModel.h"
#include "GFaceCompound.h"
+#include "discreteDiskFace.h"
#include "intersectCurveSurface.h"
#include "HilbertCurve.h"
@@ -1125,7 +1126,6 @@ void bowyerWatson(GFace *gf, int MAXPNT,
}
#endif
transferDataStructure(gf, AllTris, DATA);
- // removeThreeTrianglesNodes(gf);
}
/*
@@ -1311,6 +1311,18 @@ void optimalPointFrontalB (GFace *gf,
// P = d * (n1 cos(t) + n2 sin(t)) that is on the surface
// so we have to find t, starting with t = 0
+
+ if (gf->geomType() == GEntity::DiscreteDiskSurface){
+ discreteDiskFace *ddf = dynamic_cast<discreteDiskFace*> (gf);
+ if (ddf){
+ GPoint gp = ddf->intersectionWithCircle(n1,n2,middle,d,newPoint);
+ if (gp.succeeded()){
+ newPoint[0] = gp.u();
+ newPoint[1] = gp.v();
+ return ;
+ }
+ }
+ }
if (gf->geomType() == GEntity::CompoundSurface){
GFaceCompound *gfc = dynamic_cast<GFaceCompound*> (gf);
if (gfc){
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 51f2b6bcd1..95b79f35d7 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -332,7 +332,7 @@ void circumCenterXY(double *p1, double *p2, double *p3, double *res)
d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
if(d == 0.0) {
- // Msg::Warning("Colinear points in circum circle computation");
+ // Msg::Warning("Colinear points in circum circle computation");
res[0] = res[1] = -99999.;
return ;
}
--
GitLab