#include "INTERPKERNELDefines.hxx"
#include "InterpKernelException.hxx"
+#include "VolSurfUser.hxx"
#include "NormalizedUnstructuredMesh.hxx"
/* calcul la surface d'un triangle */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
+ inline double Surf_Tri(const double *P_1,const double *P_2,const double *P_3)
{
double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]);
double Surface = 0.5*fabs(A);
//fonction qui calcul le determinant des vecteurs: P3P1 et P3P2
//(cf doc CGAL).
- inline double mon_determinant(const double* P_1,
- const double* P_2,
- const double* P_3)
+ inline double mon_determinant(const double *P_1,
+ const double *P_2,
+ const double *P_3)
{
double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
return mon_det;
{
double X=P_1[0]-P_2[0];
double Y=P_1[1]-P_2[1];
- double norme=sqrt(X*X+Y*Y);
- return norme;
+ return sqrt(X*X+Y*Y);
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
}
/*!
- * Calculate barycentric coordinates of a point p with respect to triangle or tetra verices.
+ * Calculate barycentric coordinates of a point p with respect to triangle or tetra vertices.
* This method makes 2 assumptions :
* - this is a simplex
* - spacedim == meshdim. For TRI3 and TRI6 spaceDim is expected to be equal to 2 and for TETRA4 spaceDim is expected to be equal to 3.
*/
inline void barycentric_coords(const std::vector<const double*>& n, const double *p, double *bc)
{
- enum { _X, _Y, _Z };
+ enum { _XX=0, _YY, _ZZ };
switch(n.size())
{
case 2:
{ // TRIA3
// matrix 2x2
double
- T11 = n[0][_X]-n[2][_X], T12 = n[1][_X]-n[2][_X],
- T21 = n[0][_Y]-n[2][_Y], T22 = n[1][_Y]-n[2][_Y];
+ T11 = n[0][_XX]-n[2][_XX], T12 = n[1][_XX]-n[2][_XX],
+ T21 = n[0][_YY]-n[2][_YY], T22 = n[1][_YY]-n[2][_YY];
// matrix determinant
double Tdet = T11*T22 - T12*T21;
if ( (std::fabs( Tdet) ) < (std::numeric_limits<double>::min()) )
// matrix inverse
double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
// vector
- double r11 = p[_X]-n[2][_X], r12 = p[_Y]-n[2][_Y];
+ double r11 = p[_XX]-n[2][_XX], r12 = p[_YY]-n[2][_YY];
// barycentric coordinates: mutiply matrix by vector
bc[0] = (t11 * r11 + t12 * r12)/Tdet;
bc[1] = (t21 * r11 + t22 * r12)/Tdet;
// bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4
double T[3][4]=
- {{ n[0][_X]-n[3][_X], n[1][_X]-n[3][_X], n[2][_X]-n[3][_X], p[_X]-n[3][_X] },
- { n[0][_Y]-n[3][_Y], n[1][_Y]-n[3][_Y], n[2][_Y]-n[3][_Y], p[_Y]-n[3][_Y] },
- { n[0][_Z]-n[3][_Z], n[1][_Z]-n[3][_Z], n[2][_Z]-n[3][_Z], p[_Z]-n[3][_Z] }};
+ {{ n[0][_XX]-n[3][_XX], n[1][_XX]-n[3][_XX], n[2][_XX]-n[3][_XX], p[_XX]-n[3][_XX] },
+ { n[0][_YY]-n[3][_YY], n[1][_YY]-n[3][_YY], n[2][_YY]-n[3][_YY], p[_YY]-n[3][_YY] },
+ { n[0][_ZZ]-n[3][_ZZ], n[1][_ZZ]-n[3][_ZZ], n[2][_ZZ]-n[3][_ZZ], p[_ZZ]-n[3][_ZZ] }};
if ( !solveSystemOfEquations<3>( T, bc ) )
bc[0]=1., bc[1] = bc[2] = bc[3] = 0;
}
}
+ /*!
+ * Compute barycentric coords of point \a point relative to segment S [segStart,segStop] in 3D space.
+ * If point \a point is further from S than eps false is returned.
+ * If point \a point projection on S is outside S false is also returned.
+ * If point \a point is closer from S than eps and its projection inside S true is returned and \a bc0 and \a bc1 output parameter set.
+ */
+ inline bool IsPointOn3DSeg(const double segStart[3], const double segStop[3], const double point[3], double eps, double& bc0, double& bc1)
+ {
+ double AB[3]={segStop[0]-segStart[0],segStop[1]-segStart[1],segStop[2]-segStart[2]},AP[3]={point[0]-segStart[0],point[1]-segStart[1],point[2]-segStart[2]};
+ double l_AB(sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]));
+ double AP_dot_AB((AP[0]*AB[0]+AP[1]*AB[1]+AP[2]*AB[2])/(l_AB*l_AB));
+ double projOfPOnAB[3]={segStart[0]+AP_dot_AB*AB[0],segStart[1]+AP_dot_AB*AB[1],segStart[2]+AP_dot_AB*AB[2]};
+ double V_dist_P_AB[3]={point[0]-projOfPOnAB[0],point[1]-projOfPOnAB[1],point[2]-projOfPOnAB[2]};
+ double dist_P_AB(sqrt(V_dist_P_AB[0]*V_dist_P_AB[0]+V_dist_P_AB[1]*V_dist_P_AB[1]+V_dist_P_AB[2]*V_dist_P_AB[2]));
+ if(dist_P_AB>=eps)
+ return false;//to far from segment [segStart,segStop]
+ if(AP_dot_AB<-eps || AP_dot_AB>1.+eps)
+ return false;
+ AP_dot_AB=std::max(AP_dot_AB,0.); AP_dot_AB=std::min(AP_dot_AB,1.);
+ bc0=1.-AP_dot_AB; bc1=AP_dot_AB;
+ return true;
+ }
+
+ /*!
+ * Calculate pseudo barycentric coordinates of a point p with respect to the quadrangle vertices.
+ * This method makes the assumption that:
+ * - spacedim == meshdim (2 here).
+ * - the point is within the quad
+ * Quadratic elements are not supported yet.
+ *
+ * A quadrangle can be described as 3 vectors, one point being taken as the origin.
+ * Denoting A, B, C the three other points, any point P within the quad is written as
+ * P = xA+ yC + xy(B-A-C)
+ * This method solve those 2 equations (one per component) for x and y.
+ *
+
+ A------B
+ | |
+ | |
+ 0------C
+ */
+ inline void quad_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
+ {
+ double prec = 1.0e-14;
+ enum { _XX=0, _YY, _ZZ };
+
+ if(n.size() != 4)
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::quad_mapped_coords : unrecognized geometric type! Only QUAD4 supported.");
+
+ double A[2] = {n[1][_XX] - n[0][_XX], n[1][_YY] - n[0][_YY]};
+ double B[2] = {n[2][_XX] - n[0][_XX], n[2][_YY] - n[0][_YY]};
+ double C[2] = {n[3][_XX] - n[0][_XX], n[3][_YY] - n[0][_YY]};
+ double N[2] = {B[_XX] - A[_XX] - C[_XX], B[_YY] - A[_YY] - C[_YY]};
+ double P[2] = {p[_XX] - n[0][_XX], p[_YY] - n[0][_YY]};
+
+ // degenerated case: a rectangle:
+ if (fabs(N[0]) < prec && fabs(N[1]) < prec)
+ {
+ double det = C[0]*A[1] -C[1]*A[0];
+ if (fabs(det) < prec)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords() has a degenerated 2x2 system!");
+ bc[0] = (P[0]*A[1]-P[1]*A[0])/det;
+ bc[1] = (P[1]*C[0]-P[0]*C[1])/det;
+ return;
+ }
+ double b,c ,a = A[1]*N[0]-A[0]*N[1];
+ bool cas1;
+ if (fabs(a) > 1.0e-14)
+ {
+ b = A[1]*C[0]+N[1]*P[0]-N[0]*P[1]-A[0]*C[1];
+ c = P[0]*C[1] - P[1]*C[0];
+ cas1 = true;
+ }
+ else
+ {
+ a = -C[1]*N[0]+C[0]*N[1];
+ b = A[1]*C[0]-N[1]*P[0]+N[0]*P[1]-A[0]*C[1];
+ c = -P[0]*A[1] + P[1]*A[0];
+ cas1 = false;
+ }
+ double delta = b*b - 4.0*a*c;
+ if (delta < 0.0)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): imaginary solutions!");
+ bc[1] = 0.5*(-b+sqrt(delta))/a;
+ if (bc[1] < -prec || bc[1] > (1.0+prec))
+ bc[1] = 0.5*(-b-sqrt(delta))/a;
+ if (bc[1] < -prec || bc[1] > (1.0+prec))
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+ if (cas1)
+ {
+ double denom = C[0]+bc[1]*N[0];
+ if (fabs(denom) < prec)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+ bc[0] = (P[0]-bc[1]*A[0])/denom;
+ if (bc[0] < -prec || bc[0] > (1.0+prec))
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+ }
+ else
+ {
+ bc[0] = bc[1];
+ double denom = A[1]+bc[0]*N[1];
+ if (fabs(denom) < prec)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
+ bc[1] = (P[1]-bc[0]*C[1])/denom;
+ if (bc[1] < -prec || bc[1] > (1.0+prec))
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
+ }
+ }
+
+ /*!
+ * Doing as in quad_mapped_coords() would lead to a 4th order equation ... So go simpler here:
+ * orthogonal distance to each pair of parallel faces is computed. The ratio gives a number in [0,1]
+ *
+ * Conventions:
+ * - for HEXA8, point F (5) is taken to be the origin (see med file ref connec):
+ * 0 ------ 3
+ /| /|
+ / | / |
+ 1 ------ 2 |
+ | | | |
+ | | | |
+ | 4-----|- 7
+ | / | /
+ 5 ------ 6
+
+ *
+ */
+
+ inline void cuboid_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
+ {
+ double prec = 1.0e-14;
+ enum { _XX=0, _YY };
+ if (n.size() != 8)
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: unrecognized geometric type! Only HEXA8 supported.");
+
+ double dx1, dx2, dy1, dy2, dz1, dz2;
+ dx1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[5],n[1]);
+ dx2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[7],n[3],n[2]);
+
+ dy1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[6],n[2]);
+ dy2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[0],n[3]);
+
+ dz1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[4],n[7]);
+ dz2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[1],n[2],n[3]);
+
+ if (dx1 < -prec || dx2 < -prec || dy1 < -prec || dy2 < -prec || dz1 < -prec || dz2 < -prec)
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: point outside HEXA8");
+
+ bc[0] = dx1+dx2 < prec ? 0.5 : dx1/(dx1+dx2);
+ bc[1] = dy1+dy2 < prec ? 0.5 : dy1/(dy1+dy2);
+ bc[2] = dz1+dz2 < prec ? 0.5 : dz1/(dz1+dz2);
+ }
+
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/* calcul la surface d'un polygone. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
}
/*! Function that compares two angles from the values of the pairs (sin,cos)*/
- /*! Angles are considered in [0, 2Pi] bt are not computed explicitely */
+ /*! Angles are considered in [0, 2Pi] bt are not computed explicitly */
class AngleLess
{
public:
}
}
+ /*!
+ * Find a vector orthogonal to the input vector
+ */
+ inline void orthogonalVect3(const double inpVect[3], double outVect[3])
+ {
+ std::vector<bool> sw(3,false);
+ double inpVect2[3];
+ std::transform(inpVect,inpVect+3,inpVect2,std::ptr_fun<double,double>(fabs));
+ std::size_t posMin(std::distance(inpVect2,std::min_element(inpVect2,inpVect2+3)));
+ sw[posMin]=true;
+ std::size_t posMax(std::distance(inpVect2,std::max_element(inpVect2,inpVect2+3)));
+ if(posMax==posMin)
+ { posMax=(posMin+1)%3; }
+ sw[posMax]=true;
+ std::size_t posMid(std::distance(sw.begin(),std::find(sw.begin(),sw.end(),false)));
+ outVect[posMin]=0.; outVect[posMid]=1.; outVect[posMax]=-inpVect[posMid]/inpVect[posMax];
+ }
+
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/* Computes the dot product of a and b */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
// just to be able to compile
}
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- /* Checks wether point A is inside the quadrangle BCDE */
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* Checks whether point A is inside the quadrangle BCDE */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
const double* E,double* ABC, double* ADE)