From eb63459392b92c43230b6a7e63b28b8faf117bcf Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 17 Dec 2008 15:51:18 +0000 Subject: [PATCH] MEDMEM Industrialization 2008 + /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ + /* Calculate barycentric coordinates of a point p */ + /* with respect to triangle or tetra verices. */ + /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ + + inline void barycentric_coodrs(const std::vector& n, const double* p, double* bc) --- src/INTERP_KERNEL/InterpolationUtils.hxx | 147 +++++++++++++++++++++++ 1 file changed, 147 insertions(+) diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index 63d991e7a..87f5d4ec8 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -1,3 +1,21 @@ +// 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_ @@ -12,6 +30,7 @@ #include #include #include +#include namespace INTERP_KERNEL { @@ -147,6 +166,134 @@ 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& 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. */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ -- 2.39.2