X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FInterpolationUtils.hxx;h=2c3cb65012b4e369cda1f5604cbd592889997524;hb=4b4c848900cf2fb2cd8cb585c6ff21c7d284b87b;hp=18d550d32029b38d72a928581787f4418603e52f;hpb=a58b30a1b00a008f265f7d8d455bc444cca06cbc;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index 18d550d32..2c3cb6501 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -1,9 +1,9 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D +// Copyright (C) 2007-2016 CEA/DEN, EDF R&D // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -23,6 +23,7 @@ #include "INTERPKERNELDefines.hxx" #include "InterpKernelException.hxx" +#include "VolSurfUser.hxx" #include "NormalizedUnstructuredMesh.hxx" @@ -413,7 +414,7 @@ namespace INTERP_KERNEL } /*! - * 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. @@ -421,7 +422,7 @@ namespace INTERP_KERNEL */ inline void barycentric_coords(const std::vector& n, const double *p, double *bc) { - enum { _X, _Y, _Z }; + enum { _XX=0, _YY, _ZZ }; switch(n.size()) { case 2: @@ -435,8 +436,8 @@ namespace INTERP_KERNEL { // 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::min()) ) @@ -447,7 +448,7 @@ namespace INTERP_KERNEL // 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; @@ -462,9 +463,9 @@ namespace INTERP_KERNEL // 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; @@ -535,6 +536,136 @@ namespace INTERP_KERNEL } } + /*! + * 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& 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& 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. */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ @@ -710,7 +841,7 @@ namespace INTERP_KERNEL class AngleLess { public: - bool operator()(std::pairtheta1, std::pair theta2) + bool operator()(std::pairtheta1, std::pair theta2) const { double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second); double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second); @@ -734,7 +865,7 @@ namespace INTERP_KERNEL inline std::vector reconstruct_polygon(const std::vector& V) { - std::size_t taille=V.size(); + int taille((int)V.size()); //VB : why 6 ? @@ -749,7 +880,7 @@ namespace INTERP_KERNEL COS[0]=1.0; SIN[0]=0.0; //angle[0]=0.0; - for(std::size_t i=0; i Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]); COS[i+1]=Trigo[0]; @@ -765,7 +896,7 @@ namespace INTERP_KERNEL Pt_ordonne.reserve(taille); // std::multimap Ordre; std::multimap,int, AngleLess> CosSin; - for(std::size_t i=0;i sw(3,false); + double inpVect2[3]; + std::transform(inpVect,inpVect+3,inpVect2,std::ptr_fun(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 */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/