+// Copyright (C) 2007-2008 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.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
#ifndef _INTERPOLATIONUTILS_HXX_
#define _INTERPOLATIONUTILS_HXX_
#include <vector>
#include <algorithm>
#include <iostream>
+#include <float.h>
namespace INTERP_KERNEL
{
return Bary;
}
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Calculate barycentric coordinates of a 2D point p */
+ /* with respect to the triangle verices. */
+ /* triaCoords are in full interlace */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+ inline void barycentric_coodrs(const double* triaCoords, const double* p, double* bc)
+ {
+ // matrix 2x2
+ double
+ T11 = triaCoords[0]-triaCoords[4], T12 = triaCoords[2]-triaCoords[4],
+ T21 = triaCoords[1]-triaCoords[5], T22 = triaCoords[3]-triaCoords[5];
+ // matrix determinant
+ double Tdet = T11*T22 - T12*T21;
+ if ( fabs( Tdet ) < DBL_MIN ) {
+ bc[0]=1; bc[1]=0; bc[2]=0;
+ return;
+ }
+ // matrix inverse
+ double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
+ // vector
+ double r11 = p[0]-triaCoords[4], r12 = p[1]-triaCoords[5];
+ // barycentric coordinates: mutiply matrix by vector
+ bc[0] = (t11 * r11 + t12 * r12)/Tdet;
+ bc[1] = (t21 * r11 + t22 * r12)/Tdet;
+ bc[2] = 1. - bc[0] - bc[1];
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /* Calculate barycentric coordinates of a point p */
+ /* with respect to triangle or tetra verices. */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+ inline void barycentric_coodrs(const std::vector<const double*>& n, const double* p, double* bc)
+ {
+ enum { _X, _Y, _Z };
+ if ( n.size() == 3 ) // 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];
+ // matrix determinant
+ double Tdet = T11*T22 - T12*T21;
+ if ( std::fabs( Tdet ) < DBL_MIN ) {
+ bc[0]=1; bc[1]=bc[2]=0; // no solution
+ return;
+ }
+ // 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];
+ // barycentric coordinates: mutiply matrix by vector
+ bc[0] = (t11 * r11 + t12 * r12)/Tdet;
+ bc[1] = (t21 * r11 + t22 * r12)/Tdet;
+ bc[2] = 1. - bc[0] - bc[1];
+ }
+ else // TETRA4
+ {
+ bc[3]=0; // for no solution
+
+ // Find bc by solving system of 3 equations
+ // Equation for X:
+ // bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4
+ const int nbCol=4, nbRow=3;
+
+ double T[nbRow][nbCol]=
+ {{ 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] }};
+
+ // make upper triangular matrix
+
+ int iR[nbRow] = { 0, 1, 2 };
+
+ for ( int i = 0; i < 2; ++i ) // nullify 2 rows
+ {
+ // swap rows to have max value of i-th column in i-th row
+ double max = std::fabs( T[ iR[i] ][i] );
+ for ( int r = i+1; r < nbRow; ++r ) {
+ double t = std::fabs( T[ iR[r] ][i] );
+ if ( t > max ) {
+ max = t;
+ std::swap( iR[r], iR[i] );
+ }
+ }
+ if ( max < DBL_MIN ) {
+ bc[0]=1; bc[1]=bc[2]=bc[3]=0;
+ return; // no solution
+ }
+ // make 0 below T[i][i] (actually we do not modify i-th column)
+ double* tUpRow = T[ iR[i] ];
+ for ( int r = i+1; r < nbRow; ++r ) {
+ double* tRow = T[ iR[r] ];
+ double coef = tRow[ i ] / tUpRow[ i ];
+ for ( int c = i+1; c < nbCol; ++c )
+ tRow[ c ] -= tUpRow[ c ] * coef;
+ }
+ }
+ double* tRow = T[ iR[2] ];
+ if ( std::fabs( tRow[ 2 ] ) < DBL_MIN ) {
+ bc[0]=1; bc[1]=bc[2]=bc[3]=0;
+ return; // no solution
+ }
+ tRow[ 3 ] /= tRow[ 2 ];
+ // calculate solution: backward substitution
+
+ int r0 = iR[0], r1 = iR[1]/*, r2 = iR[2], r3 = 3*/;
+// bc[ r2 ] = tRow[ 3 ] / tRow[ 2 ];
+
+// tRow = T[ r1 ];
+// bc[ r1 ] = (tRow[ 3 ] - bc[r2]*tRow[ 2 ]) / tRow[ 1 ];
+
+// tRow = T[ r0 ];
+// bc[ r0 ] = (tRow[ 3 ] - bc[r2]*tRow[ 2 ] - bc[r1]*tRow[ 1 ]) / tRow[ 0 ];
+
+// bc[ r3 ] = 1. - bc[0] - bc[1] - bc[2];
+ bc[ 2 ] = tRow[ 3 ];
+
+ tRow = T[ iR[1] ];
+ bc[ 1 ] = (tRow[ 3 ] - bc[2]*tRow[ 2 ]) / tRow[ 1 ];
+
+ tRow = T[ iR[0] ];
+ bc[ 0 ] = (tRow[ 3 ] - bc[2]*tRow[ 2 ] - bc[1]*tRow[ 1 ]) / tRow[ 0 ];
+
+ bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2];
+ }
+ }
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/* calcul la surface d'un polygone. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
