Skip to content
Snippets Groups Projects
Commit 6abca093 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

refactorization of reparametrization in progress

parent a1f20a24
Branches
Tags
No related merge requests found
...@@ -362,6 +362,14 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params ...@@ -362,6 +362,14 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params
return; 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){ if(v->onWhat()->dim() == 0){
GVertex *gv = (GVertex*)v->onWhat(); GVertex *gv = (GVertex*)v->onWhat();
std::list<GEdge*> ed = gv->edges(); std::list<GEdge*> ed = gv->edges();
......
...@@ -452,7 +452,13 @@ GPoint discreteDiskFace::point(const double par1, const double par2) const ...@@ -452,7 +452,13 @@ GPoint discreteDiskFace::point(const double par1, const double par2) const
SPoint2 discreteDiskFace::parFromVertex(MVertex *v) 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); std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y()); if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y());
// The 1D mesh has been re-done // The 1D mesh has been re-done
...@@ -481,18 +487,19 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const ...@@ -481,18 +487,19 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
SVector3 discreteDiskFace::normal(const SPoint2 &param) const SVector3 discreteDiskFace::normal(const SPoint2 &param) const
{ {
return SVector3(); return GFace::normal(param);
} }
double discreteDiskFace::curvatureMax(const SPoint2 &param) const double discreteDiskFace::curvatureMax(const SPoint2 &param) const
{ {
throw;
return false; return false;
} }
double discreteDiskFace::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin, double discreteDiskFace::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const double *curvMax, double *curvMin) const
{ {
throw;
return false; return false;
} }
...@@ -505,7 +512,10 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const ...@@ -505,7 +512,10 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 &param) const
MTriangle* tri = NULL; MTriangle* tri = NULL;
if (ddft) tri = ddft->tri; 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 = std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it0 =
firstDerivatives.find(tri->getVertex(0)); firstDerivatives.find(tri->getVertex(0));
...@@ -583,4 +593,112 @@ void discreteDiskFace::putOnView() ...@@ -583,4 +593,112 @@ void discreteDiskFace::putOnView()
#endif #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 #endif
...@@ -60,6 +60,9 @@ class discreteDiskFace : public GFace { ...@@ -60,6 +60,9 @@ class discreteDiskFace : public GFace {
virtual void secondDer(const SPoint2 &param, virtual void secondDer(const SPoint2 &param,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const; SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
GEntity::GeomType geomType() const { return DiscreteDiskSurface; } 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: protected:
//------------------------------------------------ //------------------------------------------------
// a copy of the mesh that should not be destroyed // a copy of the mesh that should not be destroyed
......
...@@ -24,6 +24,7 @@ ...@@ -24,6 +24,7 @@
#include "Field.h" #include "Field.h"
#include "GModel.h" #include "GModel.h"
#include "GFaceCompound.h" #include "GFaceCompound.h"
#include "discreteDiskFace.h"
#include "intersectCurveSurface.h" #include "intersectCurveSurface.h"
#include "HilbertCurve.h" #include "HilbertCurve.h"
...@@ -1125,7 +1126,6 @@ void bowyerWatson(GFace *gf, int MAXPNT, ...@@ -1125,7 +1126,6 @@ void bowyerWatson(GFace *gf, int MAXPNT,
} }
#endif #endif
transferDataStructure(gf, AllTris, DATA); transferDataStructure(gf, AllTris, DATA);
// removeThreeTrianglesNodes(gf);
} }
/* /*
...@@ -1311,6 +1311,18 @@ void optimalPointFrontalB (GFace *gf, ...@@ -1311,6 +1311,18 @@ void optimalPointFrontalB (GFace *gf,
// P = d * (n1 cos(t) + n2 sin(t)) that is on the surface // P = d * (n1 cos(t) + n2 sin(t)) that is on the surface
// so we have to find t, starting with t = 0 // 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){ if (gf->geomType() == GEntity::CompoundSurface){
GFaceCompound *gfc = dynamic_cast<GFaceCompound*> (gf); GFaceCompound *gfc = dynamic_cast<GFaceCompound*> (gf);
if (gfc){ if (gfc){
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment