From: Anthony Geay Date: Thu, 18 Aug 2022 09:08:14 +0000 (+0200) Subject: Small refactoring of interpolation and implementation of localization of points regar... X-Git-Tag: V9_9_1b1^0 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=17d4c9c3365f9de973b5db520e2f50ba2e016a0a;p=tools%2Fmedcoupling.git Small refactoring of interpolation and implementation of localization of points regarding associated reference FE element. --- diff --git a/doc/developer/doxygen/doxfiles/reference/fields/discretization.dox b/doc/developer/doxygen/doxfiles/reference/fields/discretization.dox index abb5dd13b..31f4d4b32 100644 --- a/doc/developer/doxygen/doxfiles/reference/fields/discretization.dox +++ b/doc/developer/doxygen/doxfiles/reference/fields/discretization.dox @@ -21,6 +21,7 @@ A field can be supported by: - Gauss points: built with \ref MEDCoupling::TypeOfField "ON_GAUSS_PT" - Gauss points on nodes per element: built with \ref MEDCoupling::TypeOfField "ON_GAUSS_NE" - Kriging points: built with \ref MEDCoupling::TypeOfField "ON_NODES_KR" + - On points using standard FE description : built with \ref MEDCoupling::TypeOfField "ON_NODES_FE" The spatial discretization is at the center of the \ref interpolation "interpolation" mechanisms, since one of the main interpolation parameter is indeed specifying from which source discretization diff --git a/src/INTERP_KERNEL/BBTree.txx b/src/INTERP_KERNEL/BBTree.txx index b9bb357fd..227eab087 100644 --- a/src/INTERP_KERNEL/BBTree.txx +++ b/src/INTERP_KERNEL/BBTree.txx @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -31,7 +32,6 @@ constexpr double BBTREE_DFT_EPSILON = 1e-12; template class BBTree { - private: BBTree* _left; BBTree* _right; @@ -47,7 +47,7 @@ private: static const int MIN_NB_ELEMS=15; static const int MAX_LEVEL=20; public: - + BBTree() = default; /*! Constructor of the bounding box tree \param bbs pointer to the [xmin1 xmax1 ymin1 ymax1 xmin2 xmax2 ...] array containing the bounding boxes that are to be indexed. @@ -67,7 +67,6 @@ public: BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12); \endcode */ - BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON): _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon) { @@ -76,24 +75,26 @@ public: _terminal=true; } - double* nodes=new double [nbelems]; - _elems.resize(nbelems); - for (ConnType i=0; i::max(); + { + std::unique_ptr nodes( new double [nbelems] ); + _elems.resize(nbelems); + for (ConnType i=0; i(nodes, nodes+nbelems/2, nodes+nbelems); - double median = *(nodes+nbelems/2); - delete[] nodes; + std::nth_element(nodes.get(), nodes.get()+nbelems/2, nodes.get()+nbelems); + median = *(nodes.get()+nbelems/2); + } // std:: cout << *median < new_elems_left; @@ -140,8 +141,6 @@ public: _right=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_right.size(),_epsilon); } - - ~BBTree() { @@ -150,13 +149,11 @@ public: } - /*! returns in \a elems the list of elements potentially intersecting the bounding box pointed to by \a bb \param bb pointer to query bounding box \param elems list of elements (given in 0-indexing that is to say in \b C \b mode) intersecting the bounding box */ - void getIntersectingElems(const double* bb, std::vector& elems) const { // terminal node : return list of elements intersecting bb @@ -234,7 +231,6 @@ public: \param xx pointer to query point coords \param elems list of elements (given in 0-indexing) intersecting the bounding box */ - void getElementsAroundPoint(const double* xx, std::vector& elems) const { // terminal node : return list of elements intersecting bb @@ -272,8 +268,6 @@ public: _right->getElementsAroundPoint(xx,elems); } - - ConnType size() { if (_terminal) return _nbelems; diff --git a/src/INTERP_KERNEL/BBTreeStandAlone.txx b/src/INTERP_KERNEL/BBTreeStandAlone.txx new file mode 100644 index 000000000..00ea437df --- /dev/null +++ b/src/INTERP_KERNEL/BBTreeStandAlone.txx @@ -0,0 +1,38 @@ +// Copyright (C) 2022 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, 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 +// 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 +// + +#pragma once + +#include "BBTree.txx" +#include + +/*! + * Wrapper over BBTree to deal with ownership of bbox double array. + */ +template +class BBTreeStandAlone +{ +private: + std::unique_ptr _bbox; + BBTree _effective; +public: + BBTreeStandAlone(std::unique_ptr&& bbs, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON):_bbox(std::move(bbs)),_effective(_bbox.get(),nullptr,0,nbelems,epsilon) { } + void getIntersectingElems(const double* bb, std::vector& elems) const { _effective.getIntersectingElems(bb,elems); } + void getElementsAroundPoint(const double* xx, std::vector& elems) const { _effective.getElementsAroundPoint(xx,elems); } +}; diff --git a/src/INTERP_KERNEL/Bases/InterpKernelException.hxx b/src/INTERP_KERNEL/Bases/InterpKernelException.hxx index 7703b7b46..f4f05bd62 100644 --- a/src/INTERP_KERNEL/Bases/InterpKernelException.hxx +++ b/src/INTERP_KERNEL/Bases/InterpKernelException.hxx @@ -25,6 +25,7 @@ #include #include +#include namespace INTERP_KERNEL { diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index ab192e4f6..fe9e99204 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -36,10 +36,31 @@ namespace INTERP_KERNEL * */ BoundingBox::BoundingBox(const double** pts, const unsigned numPts) - :_coords(new double[6]) { - assert(numPts > 0); + initializeWith(pts,numPts); + } + + void BoundingBox::fillInXMinXmaxYminYmaxZminZmaxFormat(double data[6]) const + { + data[0] = this->getCoordinate(BoundingBox::XMIN); + data[1] = this->getCoordinate(BoundingBox::XMAX); + data[2] = this->getCoordinate(BoundingBox::YMIN); + data[3] = this->getCoordinate(BoundingBox::YMAX); + data[4] = this->getCoordinate(BoundingBox::ZMIN); + data[5] = this->getCoordinate(BoundingBox::ZMAX); + } + /** + * Constructor creating box from an array of the points corresponding + * to the vertices of the element. + * Each point is represented by an array of three doubles. + * + * @param pts array of points + * @param numPts number of vertices + * + */ + void BoundingBox::initializeWith(const double** pts, const unsigned numPts) + { // initialize with first two points const double *pt0(pts[0]); @@ -63,11 +84,8 @@ namespace INTERP_KERNEL * @param box1 the first box * @param box2 the second box */ - BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box2) - : _coords(new double[6]) + BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box2) { - assert(_coords != 0); - for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1)) { _coords[c] = std::min(box1._coords[c], box2._coords[c]); @@ -77,15 +95,6 @@ namespace INTERP_KERNEL assert(isValid()); } - /** - * Destructor - * - */ - BoundingBox::~BoundingBox() - { - delete [] _coords; - } - /** * Determines if the intersection with a given box is empty * diff --git a/src/INTERP_KERNEL/BoundingBox.hxx b/src/INTERP_KERNEL/BoundingBox.hxx index d453e8a62..34a47b17c 100644 --- a/src/INTERP_KERNEL/BoundingBox.hxx +++ b/src/INTERP_KERNEL/BoundingBox.hxx @@ -33,15 +33,20 @@ namespace INTERP_KERNEL class INTERPKERNEL_EXPORT BoundingBox { public: - /// Enumeration representing the six coordinates that define the bounding box enum BoxCoord { XMIN = 0, YMIN = 1, ZMIN = 2, XMAX = 3, YMAX = 4, ZMAX = 5 }; - + + BoundingBox() { } + BoundingBox(const double** pts, const unsigned numPts); BoundingBox(const BoundingBox& box1, const BoundingBox& box2); - ~BoundingBox(); + ~BoundingBox() { } + + void fillInXMinXmaxYminYmaxZminZmaxFormat(double data[6]) const; + + void initializeWith(const double** pts, const unsigned numPts); bool isDisjointWith(const BoundingBox& box) const; @@ -54,6 +59,8 @@ namespace INTERP_KERNEL inline void dumpCoords() const; void toCompactData(double data[6]) const; + + BoundingBox& operator=(const BoundingBox& box) = delete; private: @@ -62,12 +69,9 @@ namespace INTERP_KERNEL /// disallow copying BoundingBox(const BoundingBox& box); - /// disallow assignment - BoundingBox& operator=(const BoundingBox& box); - /// Vector containing the coordinates of the box /// interlaced in the order XMIN, YMIN, ZMIN, XMAX, YMAX, ZMAX - double* _coords; + double _coords[6]; }; diff --git a/src/INTERP_KERNEL/CMakeLists.txt b/src/INTERP_KERNEL/CMakeLists.txt index 5de054ed0..bc07afc1c 100644 --- a/src/INTERP_KERNEL/CMakeLists.txt +++ b/src/INTERP_KERNEL/CMakeLists.txt @@ -62,6 +62,9 @@ SET(interpkernel_SOURCES ExprEval/InterpKernelValue.cxx ExprEval/InterpKernelAsmX86.cxx GaussPoints/InterpKernelGaussCoords.cxx + LinearAlgebra/InterpKernelDenseMatrix.cxx + LinearAlgebra/InterpKernelLUDecomp.cxx + LinearAlgebra/InterpKernelQRDecomp.cxx ) INCLUDE_DIRECTORIES( @@ -70,6 +73,7 @@ INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR}/Geometric2D ${CMAKE_CURRENT_SOURCE_DIR}/ExprEval ${CMAKE_CURRENT_SOURCE_DIR}/GaussPoints + ${CMAKE_CURRENT_SOURCE_DIR}/LinearAlgebra ) IF (NOT DEFINED MSVC) diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx index 00c8ed03a..61feaf166 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx @@ -108,16 +108,25 @@ const double GaussInfo::HEXA27A_REF[81]={-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, } //--------------------------------------------------------------- -#define SHAPE_FUN_MACRO_BEGIN \ +#define SHAPE_FUN_MACRO_BEGIN \ for( int gaussId = 0 ; gaussId < _my_nb_gauss ; gaussId++ ) \ - { \ - double* funValue = &_my_function_value[ gaussId * _my_nb_ref ]; \ + { \ + double* funValue = &_my_function_value[ gaussId * _my_nb_ref ]; \ const double* gc = &_my_gauss_coord[ gaussId * getGaussCoordDim() ]; //--------------------------------------------------------------- -#define SHAPE_FUN_MACRO_END \ +#define SHAPE_FUN_MACRO_END \ } +#define DEV_SHAPE_FUN_MACRO_BEGIN \ + for( int gaussId = 0 ; gaussId < _my_nb_gauss ; gaussId++ ) \ + { \ + double *devFunValue = _my_derivative_func_value.data() + gaussId * getReferenceCoordDim() * _my_nb_ref; \ + const double *gc = _my_gauss_coord.data() + gaussId * getGaussCoordDim(); + +#define DEV_SHAPE_FUN_MACRO_END \ + } + #define CHECK_MACRO \ if( ! aSatify ) \ { \ @@ -157,9 +166,9 @@ GaussInfo::GaussInfo( NormalizedCellType theGeometry, _my_nb_ref(theNbRef), _my_reference_coord(theReferenceCoord) { - //Allocate shape function values _my_function_value.resize( _my_nb_gauss * _my_nb_ref ); + _my_derivative_func_value.resize( _my_nb_gauss * _my_nb_ref * getReferenceCoordDim() ); } /*! @@ -453,6 +462,80 @@ std::vector GaussInfo::NormalizeCoordinatesIfNecessary(NormalizedCellTyp return ret; } +std::vector GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellType ct) +{ + switch(ct) + { + case INTERP_KERNEL::NORM_SEG2: + return std::vector(SEG2A_REF,SEG2A_REF+sizeof(SEG2A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_SEG3: + return std::vector(SEG3_REF,SEG3_REF+sizeof(SEG3_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_TRI3: + return std::vector(TRIA3A_REF,TRIA3A_REF+sizeof(TRIA3A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_TRI6: + return std::vector(TRIA6A_REF,TRIA6A_REF+sizeof(TRIA6A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_TRI7: + return std::vector(TRIA7A_REF,TRIA7A_REF+sizeof(TRIA7A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_QUAD4: + return std::vector(QUAD4A_REF,QUAD4A_REF+sizeof(QUAD4A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_QUAD8: + return std::vector(QUAD8A_REF,QUAD8A_REF+sizeof(QUAD8A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_QUAD9: + return std::vector(QUAD9A_REF,QUAD9A_REF+sizeof(QUAD9A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_TETRA4: + return std::vector(TETRA4A_REF,TETRA4A_REF+sizeof(TETRA4A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_TETRA10: + return std::vector(TETRA10A_REF,TETRA10A_REF+sizeof(TETRA10A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_PYRA5: + return std::vector(PYRA5A_REF,PYRA5A_REF+sizeof(PYRA5A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_PYRA13: + return std::vector(PYRA13A_REF,PYRA13A_REF+sizeof(PYRA13A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_PENTA6: + return std::vector(PENTA6A_REF,PENTA6A_REF+sizeof(PENTA6A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_PENTA15: + return std::vector(PENTA15A_REF,PENTA15A_REF+sizeof(PENTA15A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_PENTA18: + return std::vector(PENTA18A_REF,PENTA18A_REF+sizeof(PENTA18A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_HEXA8: + return std::vector(HEXA8A_REF,HEXA8A_REF+sizeof(HEXA8A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_HEXA20: + return std::vector(HEXA20A_REF,HEXA20A_REF+sizeof(HEXA20A_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_HEXA27: + return std::vector(HEXA27A_REF,HEXA27A_REF+sizeof(HEXA27A_REF)/sizeof(double)); + } + THROW_IK_EXCEPTION("Input type " << ct << "is not managed by GetDefaultReferenceCoordinatesOf") +} + +/*! + * Returns true if \a ptInRefCoo is in reference cell of type \a ct (regarding GetDefaultReferenceCoordinatesOf) + * \sa GetDefaultReferenceCoordinatesOf + */ +bool GaussInfo::IsInOrOutForReference(NormalizedCellType ct, const double *ptInRefCoo, double eps) +{ + switch(ct) + { + case INTERP_KERNEL::NORM_SEG2: + case INTERP_KERNEL::NORM_SEG3: + { + return std::fabs(ptInRefCoo[0]) < 1.0+eps; + } + case INTERP_KERNEL::NORM_QUAD4: + case INTERP_KERNEL::NORM_QUAD8: + case INTERP_KERNEL::NORM_QUAD9: + { + return std::find_if(ptInRefCoo,ptInRefCoo+2,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+2; + } + case INTERP_KERNEL::NORM_HEXA8: + case INTERP_KERNEL::NORM_HEXA20: + case INTERP_KERNEL::NORM_HEXA27: + { + return std::find_if(ptInRefCoo,ptInRefCoo+3,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+3; + } + default: + THROW_IK_EXCEPTION("IsInOrOutForReference : not implemented for this geo type !"); + } +} + typedef void (*MapToShapeFunction)(GaussInfo& obj); /*! @@ -726,9 +809,17 @@ void GaussInfo::initLocalInfo() /** * Return shape function value by node id */ -const double* GaussInfo::getFunctionValues( const int theGaussId ) const +const double *GaussInfo::getFunctionValues( const int theGaussId ) const +{ + return _my_function_value.data() + _my_nb_ref*theGaussId ; +} + +/** + * Return the derivative of shape function value by node id + */ +const double *GaussInfo::getDerivativeOfShapeFunctionAt( const int theGaussId ) const { - return &_my_function_value[ _my_nb_ref*theGaussId ]; + return _my_derivative_func_value.data() + _my_nb_ref*getReferenceCoordDim()*theGaussId; } void GaussInfo::point1Init() @@ -2875,6 +2966,90 @@ void GaussInfo::hexa20aInit() funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]); SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = 0.125*(1. + gc[1] + gc[2] + 2* gc[0]) * (1.0 - gc[1])*(1.0 - gc[2]); + devFunValue[1] = 0.125*(1.0 - gc[0])*(1. + gc[0] + gc[2] + 2 * gc[1])*(1.0 - gc[2]); + devFunValue[2] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[0] + gc[1] + 2 *gc[2]); + + devFunValue[3] = 0.125*(-1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + devFunValue[4] = 0.125*(1.0 + gc[0])* (1 - gc[0] + gc[2] + 2*gc[1]) *(1.0 - gc[2]); + devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])* (1 - gc[0] + gc[1] + 2*gc[2]); + + devFunValue[6] = 0.125*(-1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + devFunValue[7] = 0.125*(1.0 + gc[0])* (-1 + gc[0] - gc[2] + 2*gc[1]) *(1.0 - gc[2]); + devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[0] - gc[1] + 2*gc[2]); + + devFunValue[9] = 0.125*(1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + devFunValue[10] = 0.125*(1.0 - gc[0])*(-1 - gc[0] - gc[2] + 2*gc[1])*(1.0 - gc[2]); + devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[0] - gc[1] + 2*gc[2]); + + devFunValue[12] = 0.125*(1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + devFunValue[13] = 0.125*(1.0 - gc[0])*(1 + gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]); + devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1] + 2*gc[2]); + + devFunValue[15] = 0.125*(-1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + devFunValue[16] = 0.125*(1.0 + gc[0])*(1 - gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]); + devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1] + 2*gc[2]); + + devFunValue[18] = 0.125*(-1.0 + gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + devFunValue[19] = 0.125*(1.0 + gc[0])*(-1 + gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]); + devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1] + 2*gc[2]); + + devFunValue[21] = 0.125*(1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + devFunValue[22] = 0.125*(1.0 - gc[0])*(-1 - gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]); + devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1] + 2*gc[2]); + + devFunValue[24] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + devFunValue[25] = 0.25*(1.0 - gc[0]*gc[0])*(-1)*(1.0 - gc[2]); + devFunValue[26] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(-1.0); + + devFunValue[27] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[2]); + devFunValue[28] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]); + devFunValue[29] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(-1.0); + + devFunValue[30] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + devFunValue[31] = 0.25*(1.0 - gc[0]*gc[0])*(1.0)*(1.0 - gc[2]); + devFunValue[32] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(-1.0); + + devFunValue[33] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 - gc[2]); + devFunValue[34] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]); + devFunValue[35] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(-1.0); + + devFunValue[36] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 - gc[1]); + devFunValue[37] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(-1.0); + devFunValue[38] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]); + + devFunValue[39] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[1]); + devFunValue[40] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(-1.0); + devFunValue[41] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]); + + devFunValue[42] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[1]); + devFunValue[43] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0]); + devFunValue[44] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]); + + devFunValue[45] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 + gc[1]); + devFunValue[46] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0]); + devFunValue[47] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]); + + devFunValue[48] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + devFunValue[49] = 0.25*(1.0 - gc[0]*gc[0])*(-1.0)*(1.0 + gc[2]); + devFunValue[50] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]); + + devFunValue[51] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[2]); + devFunValue[52] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]); + devFunValue[53] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]); + + devFunValue[54] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + devFunValue[55] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[2]); + devFunValue[56] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]); + + devFunValue[57] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 + gc[2]); + devFunValue[58] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]); + devFunValue[59] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]); + + SHAPE_FUN_MACRO_END; } /*! @@ -3191,6 +3366,118 @@ void GaussInfo::hexa27aInit() funValue[26]=(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]); SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.); + devFunValue[1] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.); + devFunValue[2] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.); + + devFunValue[3] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.); + devFunValue[4] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.); + devFunValue[5] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.); + + devFunValue[6] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.); + devFunValue[7] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.); + devFunValue[8] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.); + + devFunValue[9] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.); + devFunValue[10] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.); + devFunValue[11] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.); + + devFunValue[12] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.); + devFunValue[13] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.); + devFunValue[14] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.); + + devFunValue[15] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.); + devFunValue[16] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.); + devFunValue[17] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.); + + devFunValue[18] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.); + devFunValue[19] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.); + devFunValue[20] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.); + + devFunValue[21] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.); + devFunValue[22] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.); + devFunValue[23] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.); + + devFunValue[24] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.); + devFunValue[25] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]-1.); + devFunValue[26] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.); + + devFunValue[27] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.); + devFunValue[28] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]-1.); + devFunValue[29] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]-1.); + + devFunValue[30] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.); + devFunValue[31] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]-1.); + devFunValue[32] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.); + + devFunValue[33] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.); + devFunValue[34] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]-1.); + devFunValue[35] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]-1.); + + devFunValue[36] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.); + devFunValue[37] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]+1.); + devFunValue[38] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.); + + devFunValue[39] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.); + devFunValue[40] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]+1.); + devFunValue[41] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]+1.); + + devFunValue[42] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.); + devFunValue[43] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]+1.); + devFunValue[44] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.); + + devFunValue[45] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.); + devFunValue[46] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]+1.); + devFunValue[47] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]+1.); + + devFunValue[48] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]); + devFunValue[49] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]); + devFunValue[50] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]); + + devFunValue[51] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]); + devFunValue[52] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]); + devFunValue[53] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]); + + devFunValue[54] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]); + devFunValue[55] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]); + devFunValue[56] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]); + + devFunValue[57] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]); + devFunValue[58] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]); + devFunValue[59] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]); + + devFunValue[60] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.); + devFunValue[61] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]-1.); + devFunValue[62] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]-1.); + + devFunValue[63] = 0.5*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]); + devFunValue[64] = 0.5*gc[0]*(gc[0]-1.)*(-2*gc[1])*(1.-gc[2]*gc[2]); + devFunValue[65] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(-2*gc[2]); + + devFunValue[66] = 0.5*(-2*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]); + devFunValue[67] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*(1.-gc[2]*gc[2]); + devFunValue[68] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(-2*gc[2]); + + devFunValue[69] = 0.5*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]); + devFunValue[70] = 0.5*gc[0]*(gc[0]+1.)*(-2*gc[1])*(1.-gc[2]*gc[2]); + devFunValue[71] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(-2*gc[2]); + + devFunValue[72] = 0.5*(-2*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]); + devFunValue[73] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*(1.-gc[2]*gc[2]); + devFunValue[74] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(-2*gc[2]); + + devFunValue[75] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.); + devFunValue[76] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]+1.); + devFunValue[77] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]+1.); + + devFunValue[78] = (-2*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]); + devFunValue[79] = (1.-gc[0]*gc[0])*(-2*gc[1])*(1.-gc[2]*gc[2]); + devFunValue[80] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(-2*gc[2]); + + DEV_SHAPE_FUN_MACRO_END; } //////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx index f5b50b6fa..5afdfa7d2 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx @@ -57,12 +57,14 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT GaussInfo convertToLinear() const; - INTERPKERNEL_EXPORT const double* getFunctionValues( const int theGaussId ) const; + INTERPKERNEL_EXPORT const double *getFunctionValues( const int theGaussId ) const; + INTERPKERNEL_EXPORT const double *getDerivativeOfShapeFunctionAt( const int theGaussId ) const; INTERPKERNEL_EXPORT void initLocalInfo(); INTERPKERNEL_EXPORT static std::vector NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector& inputArray); - + INTERPKERNEL_EXPORT static std::vector GetDefaultReferenceCoordinatesOf(NormalizedCellType ct); + INTERPKERNEL_EXPORT static bool IsInOrOutForReference(NormalizedCellType ct, const double *ptInRefCoo, double eps); public: static const double SEG2A_REF[2]; static const double SEG2B_REF[2]; @@ -193,6 +195,7 @@ namespace INTERP_KERNEL int _my_local_nb_ref; //Nb of the local reference coordinates DataVector _my_function_value; //Shape Function values + DataVector _my_derivative_func_value; //Derivative of the shape function }; diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx new file mode 100644 index 000000000..38e4a6128 --- /dev/null +++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx @@ -0,0 +1,53 @@ +// Copyright (C) 2022 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, 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 +// 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 +// + +#pragma once + +#include "INTERPKERNELDefines.hxx" +#include "MCIdType.hxx" + +namespace INTERP_KERNEL +{ + template + class INTERPKERNEL_EXPORT DenseMatrixT + { + private: + mcIdType nn; + mcIdType mm; + T **v; + public: + DenseMatrixT(); + DenseMatrixT(mcIdType n, mcIdType m); + DenseMatrixT(mcIdType n, mcIdType m, const T &a); + DenseMatrixT(mcIdType n, mcIdType m, const T *a); + DenseMatrixT(const DenseMatrixT &rhs); + DenseMatrixT & operator=(const DenseMatrixT &rhs); + using value_type = T; + //! subscripting: pointer to row i + T* operator[](const mcIdType i) { return v[i]; } + const T* operator[](const mcIdType i) const { return v[i]; } + mcIdType nrows() const { return nn; } + mcIdType ncols() const { return mm; } + void resize(mcIdType newn, mcIdType newm); + void assign(mcIdType newn, mcIdType newm, const T &a); + ~DenseMatrixT(); + }; + + using DenseMatrix = DenseMatrixT; +} diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx new file mode 100644 index 000000000..9705f4fdb --- /dev/null +++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx @@ -0,0 +1,141 @@ +// Copyright (C) 2022 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, 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 +// 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 +// + +#pragma once + +#include "InterpKernelDenseMatrix.hxx" + +namespace INTERP_KERNEL +{ + template + DenseMatrixT::DenseMatrixT() : nn(0), mm(0), v(nullptr) {} + + template + DenseMatrixT::DenseMatrixT(mcIdType n, mcIdType m) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr) + { + mcIdType i,nel=m*n; + if (v) v[0] = nel>0 ? new T[nel] : nullptr; + for (i=1;i + DenseMatrixT::DenseMatrixT(mcIdType n, mcIdType m, const T &a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr) + { + mcIdType i,j,nel=m*n; + if (v) v[0] = nel>0 ? new T[nel] : nullptr; + for (i=1; i< n; i++) v[i] = v[i-1] + m; + for (i=0; i< n; i++) for (j=0; j + DenseMatrixT::DenseMatrixT(mcIdType n, mcIdType m, const T *a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr) + { + mcIdType i,j,nel=m*n; + if (v) v[0] = nel>0 ? new T[nel] : nullptr; + for (i=1; i< n; i++) v[i] = v[i-1] + m; + for (i=0; i< n; i++) for (j=0; j + DenseMatrixT::DenseMatrixT(const DenseMatrixT &rhs) : nn(rhs.nn), mm(rhs.mm), v(nn>0 ? new T*[nn] : nullptr) + { + mcIdType i,j,nel=mm*nn; + if (v) v[0] = nel>0 ? new T[nel] : nullptr; + for (i=1; i< nn; i++) v[i] = v[i-1] + mm; + for (i=0; i< nn; i++) for (j=0; j + DenseMatrixT & DenseMatrixT::operator=(const DenseMatrixT &rhs) + { + if (this != &rhs) { + mcIdType i,j,nel; + if (nn != rhs.nn || mm != rhs.mm) { + if ( v ) { + delete[] (v[0]); + delete[] (v); + } + nn=rhs.nn; + mm=rhs.mm; + v = nn>0 ? new T*[nn] : nullptr; + nel = mm*nn; + if (v) v[0] = nel>0 ? new T[nel] : nullptr; + for (i=1; i< nn; i++) v[i] = v[i-1] + mm; + } + for (i=0; i< nn; i++) for (j=0; j + void DenseMatrixT::resize(mcIdType newn, mcIdType newm) + { + mcIdType i,nel; + if (newn != nn || newm != mm) + { + if ( v ) + { + delete[] (v[0]); + delete[] (v); + } + nn = newn; + mm = newm; + v = nn>0 ? new T*[nn] : nullptr; + nel = mm*nn; + if (v) v[0] = nel>0 ? new T[nel] : nullptr; + for (i=1; i< nn; i++) v[i] = v[i-1] + mm; + } + } + + template + void DenseMatrixT::assign(mcIdType newn, mcIdType newm, const T& a) + { + mcIdType i,j,nel; + if (newn != nn || newm != mm) + { + if ( v ) + { + delete[] (v[0]); + delete[] (v); + } + nn = newn; + mm = newm; + v = nn>0 ? new T*[nn] : nullptr; + nel = mm*nn; + if (v) v[0] = nel>0 ? new T[nel] : nullptr; + for (i=1; i< nn; i++) v[i] = v[i-1] + mm; + } + for (i=0; i< nn; i++) for (j=0; j + DenseMatrixT::~DenseMatrixT() + { + if ( v ) + { + delete[] (v[0]); + delete[] (v); + } + } +} diff --git a/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx new file mode 100644 index 000000000..68eb57acb --- /dev/null +++ b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx @@ -0,0 +1,237 @@ +// Copyright (C) 2022 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, 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 +// 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 +// + +// Implementation coming from Numerical Recipes in C of 1994 (version 2.04) p 384 +// Root finding and non linear sets of equation using line searches and backtracking + +#include "InterpKernelDenseMatrix.hxx" +#include "InterpKernelLUDecomp.hxx" +#include "InterpKernelQRDecomp.hxx" +#include "InterpKernelException.hxx" +#include "MCIdType.hxx" + +#include +#include +#include + +template +inline T sqr(const T a) { return a*a; } + +namespace INTERP_KERNEL +{ + template + void LineSearches(const std::vector& xold, const double fold, const std::vector& g, std::vector &p, + std::vector &x, double &f, const double stpmax, bool &check, T &func) + { + const double ALF=1.0e-4, TOLX=std::numeric_limits::epsilon(); + double a,alam,alam2=0.0,alamin,b,disc,f2=0.0; + double rhs1,rhs2,slope=0.0,sum=0.0,temp,test,tmplam; + mcIdType i,n=xold.size(); + check=false; + for (i=0;i stpmax) + for (i=0;i= 0.0) THROW_IK_EXCEPTION("Roundoff problem in LineSearches."); + test=0.0; + for (i=0;i test) test=temp; + } + alamin=TOLX/test; + alam=1.0; + for (;;) + { + for (i=0;i0.5*alam) + tmplam=0.5*alam; + } + } + alam2=alam; + f2 = f; + alam=std::fmax(tmplam,0.1*alam); + } + } + template + class JacobianCalculator + { + private: + const double EPS; + T &func; + public: + JacobianCalculator(T &funcc) : EPS(1.0e-8),func(funcc) {} + INTERP_KERNEL::DenseMatrix operator() (const std::vector& x, const std::vector& fvec) + { + mcIdType n=x.size(); + INTERP_KERNEL::DenseMatrix df(n,n); + std::vector xh=x; + for (mcIdType j=0;j f=func(xh); + xh[j]=temp; + for (mcIdType i=0;i + class FMin + { + private: + std::vector fvec; + T &func; + mcIdType n; + public: + FMin(T &funcc) : func(funcc){} + double operator() (const std::vector& x) + { + n=x.size(); + double sum=0; + fvec=func(x); + for (mcIdType i=0;i& getVector() { return fvec; } + }; + + /*! + * check is false on normal return. + * check is true if the routine has converged to a local minimum. + */ + template + void SolveWithNewtonWithJacobian(std::vector &x, bool &check, T &vecfunc, + std::function&, const std::vector&, INTERP_KERNEL::DenseMatrix&)> jacobian) + { + const mcIdType MAXITS=200; + const double TOLF=1.0e-8,TOLMIN=1.0e-12,STPMX=100.0; + const double TOLX=std::numeric_limits::epsilon(); + mcIdType i,j,its,n=x.size(); + double den,f,fold,stpmax,sum,temp,test; + std::vector g(n),p(n),xold(n); + INTERP_KERNEL::DenseMatrix fjac(n,n); + FMin fmin(vecfunc); + std::vector &fvec=fmin.getVector(); + f=fmin(x); + test=0.0; + for (i=0;i test) test=std::abs(fvec[i]); + if (test < 0.01*TOLF) + { + check=false; + return; + } + sum=0.0; + for (i=0;i test) test=std::abs(fvec[i]); + if (test < TOLF) + { + check=false; + return; + } + if (check) + { + test=0.0; + den=std::fmax(f,0.5*double(n)); + for (i=0;i test) test=temp; + } + check=(test < TOLMIN); + return; + } + test=0.0; + for (i=0;i test) test=temp; + } + if (test < TOLX) + return; + } + THROW_IK_EXCEPTION("MAXITS exceeded in SolveWithNewtonWithJacobian"); + } + + /*! + * check is false on normal return. + * check is true if the routine has converged to a local minimum. + */ + template + void SolveWithNewton(std::vector &x, bool &check, T &vecfunc) + { + JacobianCalculator fdjac(vecfunc); + auto myJacobian = [&fdjac,vecfunc](const std::vector& x, const std::vector& fvec, INTERP_KERNEL::DenseMatrix& fjac) + { + fjac = fdjac(x,fvec); + }; + SolveWithNewtonWithJacobian(x,check,vecfunc,myJacobian); + } +} diff --git a/src/INTERP_KERNEL/Interpolation1D0D.txx b/src/INTERP_KERNEL/Interpolation1D0D.txx index dcd25a938..fc1a2e439 100755 --- a/src/INTERP_KERNEL/Interpolation1D0D.txx +++ b/src/INTERP_KERNEL/Interpolation1D0D.txx @@ -73,19 +73,17 @@ namespace INTERP_KERNEL // create BBTree structure // - get bounding boxes std::vector bboxes(2*SPACEDIM*numSrcElems); - ConnType* srcElemIdx = new ConnType[numSrcElems]; for(ConnType i = 0; i < numSrcElems ; ++i) { // get source bboxes in right order srcElems[i]->getBoundingBox()->toCompactData(bboxes.data()+6*i); - srcElemIdx[i] = srcElems[i]->getIndex(); } adjustBoundingBoxes(bboxes); const double *bboxPtr(nullptr); if(numSrcElems>0) bboxPtr=&bboxes[0]; - BBTree tree(bboxPtr, srcElemIdx, 0, numSrcElems); + BBTree tree(bboxPtr, nullptr, 0, numSrcElems); const ConnType *trgConnPtr(targetMesh.getConnectivityPtr()),*trgConnIPtr(targetMesh.getConnectivityIndexPtr()); const ConnType *srcConnPtr(srcMesh.getConnectivityPtr()),*srcConnIPtr(srcMesh.getConnectivityIndexPtr()); const double *trgCooPtr(targetMesh.getCoordinatesPtr()),*srcCooPtr(srcMesh.getCoordinatesPtr()); @@ -109,7 +107,6 @@ namespace INTERP_KERNEL } } } - delete [] srcElemIdx; for(ConnType i = 0 ; i < numSrcElems ; ++i) delete srcElems[i]; return srcMesh.getNumberOfNodes(); diff --git a/src/INTERP_KERNEL/Interpolation2D3D.hxx b/src/INTERP_KERNEL/Interpolation2D3D.hxx index b2698b18e..21fcedd4e 100644 --- a/src/INTERP_KERNEL/Interpolation2D3D.hxx +++ b/src/INTERP_KERNEL/Interpolation2D3D.hxx @@ -58,9 +58,9 @@ namespace INTERP_KERNEL protected: template - void performAdjustmentOfBB(Intersector3D* intersector, std::vector& bbox) const + void performAdjustmentOfBB(Intersector3D* intersector, double *bbox, std::size_t sz) const { - intersector->adjustBoundingBoxes(bbox,InterpolationOptions::getBoundingBoxAdjustment(),InterpolationOptions::getBoundingBoxAdjustmentAbs()); + intersector->adjustBoundingBoxes(bbox,sz,InterpolationOptions::getBoundingBoxAdjustment(),InterpolationOptions::getBoundingBoxAdjustmentAbs()); } private: diff --git a/src/INTERP_KERNEL/Interpolation2D3D.txx b/src/INTERP_KERNEL/Interpolation2D3D.txx index 6069c689e..73960e815 100755 --- a/src/INTERP_KERNEL/Interpolation2D3D.txx +++ b/src/INTERP_KERNEL/Interpolation2D3D.txx @@ -32,7 +32,7 @@ #include "PointLocator3DIntersectorP1P0.txx" #include "PolyhedronIntersectorP1P1.txx" #include "PointLocator3DIntersectorP1P1.txx" -#include "Log.hxx" +#include "InterpolationHelper.txx" #include "BBTree.txx" @@ -67,24 +67,13 @@ namespace INTERP_KERNEL { typedef typename MyMeshType::MyConnType ConnType; // create MeshElement objects corresponding to each element of the two meshes - const ConnType numSrcElems = srcMesh.getNumberOfElements(); const ConnType numTargetElems = targetMesh.getNumberOfElements(); - LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); + LOG(2, "Target mesh has " << numTargetElems << " elements "); - std::vector*> srcElems(numSrcElems); - std::vector*> targetElems(numTargetElems); + DuplicateFacesType intersectFaces; - std::map*, ConnType> indices; - DuplicateFacesType intersectFaces; - - for(ConnType i = 0 ; i < numSrcElems ; ++i) - srcElems[i] = new MeshElement(i, srcMesh); - - for(ConnType i = 0 ; i < numTargetElems ; ++i) - targetElems[i] = new MeshElement(i, targetMesh); - - Intersector3D* intersector=0; + std::unique_ptr< Intersector3D > intersector; std::string methC = InterpolationOptions::filterInterpolationMethod(method); const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel()); if(methC=="P0P0") @@ -92,12 +81,12 @@ namespace INTERP_KERNEL switch(InterpolationOptions::getIntersectionType()) { case Triangulation: - intersector=new Polyhedron3D2DIntersectorP0P0(targetMesh, + intersector.reset( new Polyhedron3D2DIntersectorP0P0(targetMesh, srcMesh, dimCaracteristic, getPrecision(), intersectFaces, - getSplittingPolicy()); + getSplittingPolicy()) ); break; default: throw INTERP_KERNEL::Exception("Invalid 2D to 3D intersection type for P0P0 interp specified : must be Triangulation."); @@ -109,56 +98,29 @@ namespace INTERP_KERNEL matrix.resize(intersector->getNumberOfRowsOfResMatrix()); // create BBTree structure - // - get bounding boxes - std::vector bboxes(6 * numSrcElems); - ConnType* srcElemIdx = new ConnType[numSrcElems]; - for(ConnType i = 0; i < numSrcElems ; ++i) - { - // get source bboxes in right order - const BoundingBox* box = srcElems[i]->getBoundingBox(); - bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); - bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); - bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); - bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); - bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); - bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); - - // source indices have to begin with zero for BBox, I think - srcElemIdx[i] = srcElems[i]->getIndex(); - } - // [ABN] Adjust 2D bounding box (those might be flat in the cases where the 2D surf are perfectly aligned with the axis) - performAdjustmentOfBB(intersector, bboxes); - - BBTree<3,ConnType> tree(bboxes.data(), srcElemIdx, 0, numSrcElems, 0.); + BBTreeStandAlone<3,ConnType> tree( BuildBBTreeWithAdjustment(srcMesh,[this,&intersector](double *bbox, typename MyMeshType::MyConnType sz){ this->performAdjustmentOfBB(intersector.get(),bbox,sz); }) ); // for each target element, get source elements with which to calculate intersection // - calculate intersection by calling intersectCells for(ConnType i = 0; i < numTargetElems; ++i) { - const BoundingBox* box = targetElems[i]->getBoundingBox(); - const ConnType targetIdx = targetElems[i]->getIndex(); + MeshElement trgMeshElem(i, targetMesh); + + const BoundingBox *box = trgMeshElem.getBoundingBox(); // get target bbox in right order double targetBox[6]; - targetBox[0] = box->getCoordinate(BoundingBox::XMIN); - targetBox[1] = box->getCoordinate(BoundingBox::XMAX); - targetBox[2] = box->getCoordinate(BoundingBox::YMIN); - targetBox[3] = box->getCoordinate(BoundingBox::YMAX); - targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); - targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + box->fillInXMinXmaxYminYmaxZminZmaxFormat(targetBox); std::vector intersectElems; tree.getIntersectingElems(targetBox, intersectElems); if ( !intersectElems.empty() ) - intersector->intersectCells(targetIdx, intersectElems, matrix); - + intersector->intersectCells(i, intersectElems, matrix); } - delete [] srcElemIdx; - DuplicateFacesType::iterator iter; for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter) { @@ -171,16 +133,6 @@ namespace INTERP_KERNEL // free allocated memory ConnType ret=intersector->getNumberOfColsOfResMatrix(); - delete intersector; - - for(ConnType i = 0 ; i < numSrcElems ; ++i) - { - delete srcElems[i]; - } - for(ConnType i = 0 ; i < numTargetElems ; ++i) - { - delete targetElems[i]; - } return ret; } diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx index 355217c96..c93dffcbb 100755 --- a/src/INTERP_KERNEL/Interpolation3D.txx +++ b/src/INTERP_KERNEL/Interpolation3D.txx @@ -47,10 +47,12 @@ #else // use BBTree class -#include "BBTree.txx" +#include "InterpolationHelper.txx" #endif +#include + namespace INTERP_KERNEL { /** @@ -77,35 +79,23 @@ namespace INTERP_KERNEL template typename MyMeshType::MyConnType Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method) { - typedef typename MyMeshType::MyConnType ConnType; + using ConnType = typename MyMeshType::MyConnType; // create MeshElement objects corresponding to each element of the two meshes - const ConnType numSrcElems = srcMesh.getNumberOfElements(); const ConnType numTargetElems = targetMesh.getNumberOfElements(); - LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); - - std::vector*> srcElems(numSrcElems); - std::vector*> targetElems(numTargetElems); - - std::map*, ConnType> indices; + LOG(2, "Target mesh has " << numTargetElems << " elements "); - for(ConnType i = 0 ; i < numSrcElems ; ++i) - srcElems[i] = new MeshElement(i, srcMesh); - - for(ConnType i = 0 ; i < numTargetElems ; ++i) - targetElems[i] = new MeshElement(i, targetMesh); - - Intersector3D* intersector=0; + std::unique_ptr> intersector; std::string methC = InterpolationOptions::filterInterpolationMethod(method); if(methC=="P0P0") { switch(InterpolationOptions::getIntersectionType()) { case Triangulation: - intersector=new PolyhedronIntersectorP0P0(targetMesh, srcMesh, getSplittingPolicy()); + intersector.reset( new PolyhedronIntersectorP0P0(targetMesh, srcMesh, getSplittingPolicy()) ); break; case PointLocator: - intersector=new PointLocator3DIntersectorP0P0(targetMesh, srcMesh, getPrecision()); + intersector.reset( new PointLocator3DIntersectorP0P0(targetMesh, srcMesh, getPrecision()) ); break; default: throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangle or PointLocator."); @@ -116,10 +106,10 @@ namespace INTERP_KERNEL switch(InterpolationOptions::getIntersectionType()) { case Triangulation: - intersector=new PolyhedronIntersectorP0P1(targetMesh, srcMesh, getSplittingPolicy()); + intersector.reset( new PolyhedronIntersectorP0P1(targetMesh, srcMesh, getSplittingPolicy()) ); break; case PointLocator: - intersector=new PointLocator3DIntersectorP0P1(targetMesh, srcMesh, getPrecision()); + intersector.reset( new PointLocator3DIntersectorP0P1(targetMesh, srcMesh, getPrecision()) ); break; default: throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P1 interp specified : must be Triangle or PointLocator."); @@ -130,13 +120,13 @@ namespace INTERP_KERNEL switch(InterpolationOptions::getIntersectionType()) { case Triangulation: - intersector=new PolyhedronIntersectorP1P0(targetMesh, srcMesh, getSplittingPolicy()); + intersector.reset( new PolyhedronIntersectorP1P0(targetMesh, srcMesh, getSplittingPolicy()) ); break; case PointLocator: - intersector=new PointLocator3DIntersectorP1P0(targetMesh, srcMesh, getPrecision()); + intersector.reset( new PointLocator3DIntersectorP1P0(targetMesh, srcMesh, getPrecision()) ); break; case Barycentric: - intersector=new PolyhedronIntersectorP1P0Bary(targetMesh, srcMesh, getSplittingPolicy()); + intersector.reset( new PolyhedronIntersectorP1P0Bary(targetMesh, srcMesh, getSplittingPolicy()) ); break; default: throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P0 interp specified : must be Triangle, PointLocator or Barycentric."); @@ -147,16 +137,16 @@ namespace INTERP_KERNEL switch(InterpolationOptions::getIntersectionType()) { case Triangulation: - intersector=new PolyhedronIntersectorP1P1(targetMesh, srcMesh, getSplittingPolicy()); + intersector.reset( new PolyhedronIntersectorP1P1(targetMesh, srcMesh, getSplittingPolicy()) ); break; case PointLocator: - intersector=new PointLocator3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()); + intersector.reset( new PointLocator3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()) ); break; case Barycentric: - intersector=new Barycentric3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()); + intersector.reset( new Barycentric3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()) ); break; case MappedBarycentric: - intersector=new MappedBarycentric3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()); + intersector.reset( new MappedBarycentric3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()) ); break; default: throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle, PointLocator, Barycentric or MappedBarycentric."); @@ -303,71 +293,31 @@ namespace INTERP_KERNEL } #else // Use BBTree - - // create BBTree structure - // - get bounding boxes - double* bboxes = new double[6 * numSrcElems]; - ConnType* srcElemIdx = new ConnType[numSrcElems]; - for(ConnType i = 0; i < numSrcElems ; ++i) - { - // get source bboxes in right order - const BoundingBox* box = srcElems[i]->getBoundingBox(); - bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); - bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); - bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); - bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); - bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); - bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); - - // source indices have to begin with zero for BBox, I think - srcElemIdx[i] = srcElems[i]->getIndex(); - } - - BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems); + // create BBTree structure + BBTreeStandAlone<3,ConnType> tree( BuildBBTree(srcMesh) ); // for each target element, get source elements with which to calculate intersection // - calculate intersection by calling intersectCells for(ConnType i = 0; i < numTargetElems; ++i) { - const BoundingBox* box = targetElems[i]->getBoundingBox(); - const ConnType targetIdx = targetElems[i]->getIndex(); + MeshElement trgMeshElem(i, targetMesh); + + const BoundingBox *box = trgMeshElem.getBoundingBox(); // get target bbox in right order double targetBox[6]; - targetBox[0] = box->getCoordinate(BoundingBox::XMIN); - targetBox[1] = box->getCoordinate(BoundingBox::XMAX); - targetBox[2] = box->getCoordinate(BoundingBox::YMIN); - targetBox[3] = box->getCoordinate(BoundingBox::YMAX); - targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); - targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + box->fillInXMinXmaxYminYmaxZminZmaxFormat(targetBox); std::vector intersectElems; tree.getIntersectingElems(targetBox, intersectElems); if ( !intersectElems.empty() ) - intersector->intersectCells(targetIdx,intersectElems,result); + intersector->intersectCells(i,intersectElems,result); } - delete [] bboxes; - delete [] srcElemIdx; - #endif - // free allocated memory - ConnType ret=intersector->getNumberOfColsOfResMatrix(); - - delete intersector; - - for(ConnType i = 0 ; i < numSrcElems ; ++i) - { - delete srcElems[i]; - } - for(ConnType i = 0 ; i < numTargetElems ; ++i) - { - delete targetElems[i]; - } - return ret; - + return intersector->getNumberOfColsOfResMatrix(); } } diff --git a/src/INTERP_KERNEL/Interpolation3D1D.cxx b/src/INTERP_KERNEL/Interpolation3D1D.cxx index 28e9c438a..f64146c40 100644 --- a/src/INTERP_KERNEL/Interpolation3D1D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D1D.cxx @@ -27,16 +27,13 @@ namespace INTERP_KERNEL Interpolation3D1D::Interpolation3D1D(const InterpolationOptions& io):Interpolation(io) { } - /** - * Inspired from PlanarIntersector::adjustBoundingBoxes - */ - void Interpolation3D1D::adjustBoundingBoxes(std::vector& bbox) + void Interpolation3D1D::adjustBoundingBoxes(double *bbox, std::size_t sz) { const int SPACE_DIM = 3; const double adj = getBoundingBoxAdjustmentAbs(); const double adjRel = getBoundingBoxAdjustment(); - std::size_t size = bbox.size()/(2*SPACE_DIM); + std::size_t size = sz/(2*SPACE_DIM); for (std::size_t i=0; i::max(); @@ -52,4 +49,12 @@ namespace INTERP_KERNEL } } } + + /** + * Inspired from PlanarIntersector::adjustBoundingBoxes + */ + void Interpolation3D1D::adjustBoundingBoxes(std::vector& bbox) + { + adjustBoundingBoxes(bbox.data(),bbox.size()); + } } diff --git a/src/INTERP_KERNEL/Interpolation3D1D.hxx b/src/INTERP_KERNEL/Interpolation3D1D.hxx index 2015bb210..187309316 100755 --- a/src/INTERP_KERNEL/Interpolation3D1D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D1D.hxx @@ -16,18 +16,18 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// Author : A Bruneton (CEA/DEN) - -#pragma once - -#include "INTERPKERNELDefines.hxx" -#include "Interpolation.hxx" -#include "NormalizedUnstructuredMesh.hxx" -#include "InterpolationOptions.hxx" - -#include - -namespace INTERP_KERNEL +// Author : A Bruneton (CEA/DEN) + +#pragma once + +#include "INTERPKERNELDefines.hxx" +#include "Interpolation.hxx" +#include "NormalizedUnstructuredMesh.hxx" +#include "InterpolationOptions.hxx" + +#include + +namespace INTERP_KERNEL { /** * \class Interpolation3D1D @@ -35,15 +35,16 @@ namespace INTERP_KERNEL * Can be seen as a specialization of Interpolation3D, and allows notably the adjustment of bounding boxes. */ - class INTERPKERNEL_EXPORT Interpolation3D1D : public Interpolation - { - public: - Interpolation3D1D(); - Interpolation3D1D(const InterpolationOptions& io); - template - typename MyMeshType::MyConnType interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method); - private: - void adjustBoundingBoxes(std::vector& bbox); - }; -} - + class INTERPKERNEL_EXPORT Interpolation3D1D : public Interpolation + { + public: + Interpolation3D1D(); + Interpolation3D1D(const InterpolationOptions& io); + template + typename MyMeshType::MyConnType interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method); + private: + void adjustBoundingBoxes(double *bbox, std::size_t sz); + void adjustBoundingBoxes(std::vector& bbox); + }; +} + diff --git a/src/INTERP_KERNEL/Interpolation3D1D.txx b/src/INTERP_KERNEL/Interpolation3D1D.txx index fdec33c6c..6a619ee14 100755 --- a/src/INTERP_KERNEL/Interpolation3D1D.txx +++ b/src/INTERP_KERNEL/Interpolation3D1D.txx @@ -29,9 +29,7 @@ #include "PointLocator3DIntersectorP1P1.txx" #include "Log.hxx" -#include "BBTree.txx" - -#include +#include "InterpolationHelper.txx" namespace INTERP_KERNEL { @@ -47,21 +45,9 @@ namespace INTERP_KERNEL typedef typename MyMeshType::MyConnType ConnType; // create MeshElement objects corresponding to each element of the two meshes - const ConnType numSrcElems = srcMesh.getNumberOfElements(); const ConnType numTargetElems = targetMesh.getNumberOfElements(); - LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); - - std::vector< std::unique_ptr< MeshElement > > srcElems(numSrcElems); - std::vector< std::unique_ptr< MeshElement > > targetElems(numTargetElems); - - std::map*, ConnType> indices; - - for(ConnType i = 0 ; i < numSrcElems ; ++i) - srcElems[i].reset( new MeshElement(i, srcMesh) ); - - for(ConnType i = 0 ; i < numTargetElems ; ++i) - targetElems[i].reset( new MeshElement(i, targetMesh) ); + LOG(2, "Target mesh has " << numTargetElems << " elements "); std::unique_ptr< Intersector3D > intersector; std::string methC = InterpolationOptions::filterInterpolationMethod(method); @@ -82,52 +68,25 @@ namespace INTERP_KERNEL // create empty maps for all source elements result.resize(intersector->getNumberOfRowsOfResMatrix()); - // create BBTree structure - // - get bounding boxes - std::vector bboxes(6*numSrcElems); - std::unique_ptr srcElemIdx{ new ConnType[numSrcElems] }; - for(ConnType i = 0; i < numSrcElems ; ++i) - { - // get source bboxes in right order - const BoundingBox* box = srcElems[i]->getBoundingBox(); - bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); - bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); - bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); - bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); - bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); - bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); - - srcElemIdx[i] = srcElems[i]->getIndex(); - } - - adjustBoundingBoxes(bboxes); - const double *bboxPtr = nullptr; - if(numSrcElems>0) - bboxPtr=bboxes.data(); - BBTree<3,ConnType> tree(bboxPtr, srcElemIdx.get(), 0, numSrcElems); + BBTreeStandAlone<3,ConnType> tree( BuildBBTreeWithAdjustment(srcMesh,[this](double *bbox, typename MyMeshType::MyConnType sz){ this->adjustBoundingBoxes(bbox,sz); }) ); // for each target element, get source elements with which to calculate intersection // - calculate intersection by calling intersectCells for(ConnType i = 0; i < numTargetElems; ++i) { - const BoundingBox* box = targetElems[i]->getBoundingBox(); - const ConnType targetIdx = targetElems[i]->getIndex(); - + MeshElement trgMeshElem(i, targetMesh); + + const BoundingBox* box = trgMeshElem.getBoundingBox(); // get target bbox in right order double targetBox[6]; - targetBox[0] = box->getCoordinate(BoundingBox::XMIN); - targetBox[1] = box->getCoordinate(BoundingBox::XMAX); - targetBox[2] = box->getCoordinate(BoundingBox::YMIN); - targetBox[3] = box->getCoordinate(BoundingBox::YMAX); - targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); - targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + box->fillInXMinXmaxYminYmaxZminZmaxFormat(targetBox); std::vector intersectElems; tree.getIntersectingElems(targetBox, intersectElems); if ( !intersectElems.empty() ) - intersector->intersectCells(targetIdx,intersectElems,result); + intersector->intersectCells(i,intersectElems,result); } return intersector->getNumberOfColsOfResMatrix(); } diff --git a/src/INTERP_KERNEL/InterpolationHelper.txx b/src/INTERP_KERNEL/InterpolationHelper.txx new file mode 100755 index 000000000..48ca55d1f --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationHelper.txx @@ -0,0 +1,58 @@ +// Copyright (C) 2022 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, 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 +// 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 +// +// Author : Anthony Geay (EDF R&D) + +#pragma once + +#include "BBTreeStandAlone.txx" +#include "MeshElement.txx" +#include "Log.hxx" + +#include +#include + +namespace INTERP_KERNEL +{ + template + BBTreeStandAlone<3,typename MyMeshType::MyConnType> BuildBBTree(const MyMeshType& srcMesh) + { + return BuildBBTreeWithAdjustment(srcMesh,[](double *,typename MyMeshType::MyConnType){}); + } + + template + BBTreeStandAlone<3,typename MyMeshType::MyConnType> BuildBBTreeWithAdjustment(const MyMeshType& srcMesh, std::function bboxAdjuster) + { + using ConnType = typename MyMeshType::MyConnType; + const ConnType numSrcElems = srcMesh.getNumberOfElements(); + LOG(2, "Source mesh has " << numSrcElems << " elements"); + // create BBTree structure + // - get bounding boxes + const ConnType nbElts = 6 * numSrcElems; + std::unique_ptr bboxes( new double[nbElts] ); + for(ConnType i = 0; i < numSrcElems ; ++i) + { + MeshElement srcElem(i,srcMesh); + // get source bboxes in right order + const BoundingBox *box( srcElem.getBoundingBox() ); + box->fillInXMinXmaxYminYmaxZminZmaxFormat(bboxes.get()+6*i); + } + bboxAdjuster(bboxes.get(),nbElts); + return BBTreeStandAlone<3,ConnType>(std::move(bboxes),numSrcElems); + } +} diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelDenseMatrix.cxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelDenseMatrix.cxx new file mode 100644 index 000000000..141b3b935 --- /dev/null +++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelDenseMatrix.cxx @@ -0,0 +1,22 @@ +// Copyright (C) 2022 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, 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 +// 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 +// + +#include "InterpKernelDenseMatrix.txx" + +template class INTERPKERNEL_EXPORT INTERP_KERNEL::DenseMatrixT; diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx new file mode 100644 index 000000000..3416e5ad2 --- /dev/null +++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx @@ -0,0 +1,142 @@ +// Copyright (C) 2022 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, 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 +// 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 +// + +// Implementation coming from Numerical Recipes in C of 1994 (version 2.04) + +#include "InterpKernelLUDecomp.hxx" +#include "InterpKernelException.hxx" + +#include + +using namespace INTERP_KERNEL; + +LUDecomp::LUDecomp(const INTERP_KERNEL::DenseMatrix& a) : n(a.nrows()), lu(a), aref(a), indx(n) +{ + const double TINY=1.0e-40; + mcIdType i,imax,j,k; + double big,temp; + std::vector vv(n); + d=1.0; + for (i=0;i big) big=temp; + if (big == 0.0) THROW_IK_EXCEPTION("Singular matrix in LUDecomp"); + vv[i]=1.0/big; + } + for (k=0;k big) + { + big=temp; + imax=i; + } + } + if (k != imax) { + for (j=0;j& b, std::vector &x) +{ + mcIdType i,ii=0,ip,j; + double sum; + if (b.size() != ToSizeT(n) || x.size() != ToSizeT(n)) + THROW_IK_EXCEPTION("LUDecomp::solve bad sizes"); + for (i=0;i=0;i--) { + sum=x[i]; + for (j=i+1;j xx(n); + for (j=0;j& b, std::vector &x) +{ + mcIdType i,j; + std::vector r(n); + for (i=0;i + +namespace INTERP_KERNEL +{ + class LUDecomp + { + private: + mcIdType n; + INTERP_KERNEL::DenseMatrix lu; + std::vector indx; + double d; + const INTERP_KERNEL::DenseMatrix& aref; + public: + LUDecomp (const INTERP_KERNEL::DenseMatrix& a); + void solve(const std::vector& b, std::vector &x); + void solve(const INTERP_KERNEL::DenseMatrix& b, INTERP_KERNEL::DenseMatrix &x); + void inverse(INTERP_KERNEL::DenseMatrix &ainv); + double det(); + void mprove(const std::vector& b, std::vector &x); + }; +} diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx new file mode 100644 index 000000000..e180a9b8a --- /dev/null +++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx @@ -0,0 +1,170 @@ +// Copyright (C) 2022 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, 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 +// 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 +// + +// Implementation coming from Numerical Recipes in C of 1994 (version 2.04) + +#include "InterpKernelQRDecomp.hxx" +#include "InterpKernelException.hxx" + +#include + +using namespace INTERP_KERNEL; + +template +inline T sqr(const T a) { return a*a; } + +template +inline T sign(const T &a, const T &b) { return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); } + +QRDecomp::QRDecomp(const INTERP_KERNEL::DenseMatrix& a): n(a.nrows()), qt(n,n), r(a), sing(false) +{ + mcIdType i,j,k; + std::vector c(n), d(n); + double scale,sigma,sum,tau; + for (k=0;k& b, std::vector &x) +{ + qtmult(b,x); + rsolve(x,x); +} + +void QRDecomp::qtmult(const std::vector& b, std::vector &x) +{ + mcIdType i,j; + double sum; + for (i=0;i& b, std::vector &x) +{ + mcIdType i,j; + double sum; + if (sing) THROW_IK_EXCEPTION("attempting solve in a singular QR"); + for (i=n-1;i>=0;i--) { + sum=b[i]; + for (j=i+1;j& u, const std::vector& v) +{ + mcIdType i,k; + std::vector w(u); + for (k=n-1;k>=0;k--) + if (w[k] != 0.0) break; + if (k < 0) k=0; + for (i=k-1;i>=0;i--) { + rotate(i,w[i],-w[i+1]); + if (w[i] == 0.0) + w[i]=std::abs(w[i+1]); + else if (std::abs(w[i]) > std::abs(w[i+1])) + w[i]=std::abs(w[i])*std::sqrt(1.0+sqr(w[i+1]/w[i])); + else w[i]=std::abs(w[i+1])*std::sqrt(1.0+sqr(w[i]/w[i+1])); + } + for (i=0;i= 0.0 ? 1.0 : -1.0); + } + else if (std::abs(a) > std::abs(b)) + { + fact=b/a; + c=sign(1.0/std::sqrt(1.0+(fact*fact)),a); + s=fact*c; + } + else + { + fact=a/b; + s=sign(1.0/std::sqrt(1.0+(fact*fact)),b); + c=fact*s; + } + for (j=i;j + +namespace INTERP_KERNEL +{ + struct QRDecomp + { + private: + mcIdType n; + INTERP_KERNEL::DenseMatrix qt; + INTERP_KERNEL::DenseMatrix r; + bool sing; + public: + QRDecomp(const INTERP_KERNEL::DenseMatrix& a); + void solve(const std::vector& b, std::vector &x); + void qtmult(const std::vector& b, std::vector &x); + void rsolve(const std::vector& b, std::vector &x); + void update(const std::vector& u, const std::vector& v); + void rotate(const mcIdType i, const double a, const double b); + }; +} diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx index b0e7b8272..900e5d533 100644 --- a/src/INTERP_KERNEL/MeshElement.hxx +++ b/src/INTERP_KERNEL/MeshElement.hxx @@ -41,28 +41,27 @@ namespace INTERP_KERNEL template MeshElement(const ConnType index, const MyMeshType& mesh); - ~MeshElement(); - - ConnType getIndex() const { return _index; } + MeshElement() = default; + + template + void assign(const ConnType index, const MyMeshType& mesh); + + ~MeshElement() { } nbnodesincelltype getNumberOfNodes() const { return _number; } - const BoundingBox* getBoundingBox() const { return _box; } + const BoundingBox *getBoundingBox() const { return &_box; } + + MeshElement& operator=(const MeshElement& elem) = delete; private: /// disallow copying MeshElement(const MeshElement& elem); - /// disallow assignment - MeshElement& operator=(const MeshElement& elem); - - /// global number of the element - const ConnType _index; - nbnodesincelltype _number; /// bounding box of the element - does not change after having been initialised - BoundingBox* _box; + BoundingBox _box; }; /** diff --git a/src/INTERP_KERNEL/MeshElement.txx b/src/INTERP_KERNEL/MeshElement.txx index 2a21431ec..7a716cff7 100755 --- a/src/INTERP_KERNEL/MeshElement.txx +++ b/src/INTERP_KERNEL/MeshElement.txx @@ -28,6 +28,7 @@ #include #include #include +#include namespace INTERP_KERNEL { @@ -40,36 +41,30 @@ namespace INTERP_KERNEL */ template template - MeshElement::MeshElement(const ConnType index, const MyMeshType& mesh) - : _index(index), _number( 0 ), _box(nullptr) + MeshElement::MeshElement(const ConnType index, const MyMeshType& mesh): _number( 0 ) + { + this->assign(index,mesh); + } + + template + template + void MeshElement::assign(const ConnType index, const MyMeshType& mesh) { auto numberCore = mesh.getNumberOfNodesOfElement(OTT::indFC(index)); if(numberCore < std::numeric_limits::max()) { _number = static_cast< nbnodesincelltype >(numberCore); - const double**vertices = new const double*[_number]; + std::unique_ptr vertices( new const double*[_number] ); for( nbnodesincelltype i = 0 ; i < _number ; ++i) vertices[i] = getCoordsOfNode(i , OTT::indFC(index), mesh); - // create bounding box - _box = new BoundingBox(vertices,_number); - delete [] vertices; + _box.initializeWith(vertices.get(),_number); } else { - std::cerr << "ERROR at index " << index << " : exceeding capacity !" << std::endl; + THROW_IK_EXCEPTION("ERROR at index " << index << " : exceeding capacity !"); } } - - /** - * Destructor - * - */ - template - MeshElement::~MeshElement() - { - delete _box; - } ///////////////////////////////////////////////////////////////////// /// ElementBBoxOrder ///////////// diff --git a/src/INTERP_KERNEL/TargetIntersector.hxx b/src/INTERP_KERNEL/TargetIntersector.hxx index 84f29e87f..c0613c7de 100644 --- a/src/INTERP_KERNEL/TargetIntersector.hxx +++ b/src/INTERP_KERNEL/TargetIntersector.hxx @@ -52,6 +52,7 @@ namespace INTERP_KERNEL virtual ~TargetIntersector() { } void adjustBoundingBoxes(std::vector& bbox, double adjustmentEps, double adjustmentEpsAbs); + void adjustBoundingBoxes(double *bbox, std::size_t sz, double adjustmentEps, double adjustmentEpsAbs); }; } diff --git a/src/INTERP_KERNEL/TargetIntersector.txx b/src/INTERP_KERNEL/TargetIntersector.txx index 752784632..0cd9149fa 100644 --- a/src/INTERP_KERNEL/TargetIntersector.txx +++ b/src/INTERP_KERNEL/TargetIntersector.txx @@ -26,6 +26,12 @@ namespace INTERP_KERNEL { + template + void TargetIntersector::adjustBoundingBoxes(std::vector& bbox, double adjustmentEps, double adjustmentEpsAbs) + { + this->adjustBoundingBoxes(bbox.data(),bbox.size(),adjustmentEps,adjustmentEpsAbs); + } + /*! Readjusts a set of bounding boxes so that they are extended in all dimensions for avoiding missing interesting intersections @@ -34,9 +40,9 @@ namespace INTERP_KERNEL @param adjustmentEpsAbs absolute adjustment value (added on each side of the BBox in each dimension) */ template - void TargetIntersector::adjustBoundingBoxes(std::vector& bbox, double adjustmentEps, double adjustmentEpsAbs) + void TargetIntersector::adjustBoundingBoxes(double *bbox, std::size_t sz, double adjustmentEps, double adjustmentEpsAbs) { - std::size_t size = bbox.size()/(2*SPACEDIM); + std::size_t size = sz/(2*SPACEDIM); for (std::size_t i=0; i::max(); diff --git a/src/MEDCoupling/CMakeLists.txt b/src/MEDCoupling/CMakeLists.txt index 8a18a7a55..5bd1fccdd 100644 --- a/src/MEDCoupling/CMakeLists.txt +++ b/src/MEDCoupling/CMakeLists.txt @@ -37,6 +37,7 @@ INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL/Geometric2D ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL/ExprEval ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL/GaussPoints + ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL/LinearAlgebra ) SET(medcoupling_SOURCES @@ -61,6 +62,7 @@ SET(medcoupling_SOURCES MEDCouplingStructuredMesh.cxx MEDCouplingTimeDiscretization.cxx MEDCouplingFieldDiscretization.cxx + MEDCouplingFieldDiscretizationOnNodesFE.cxx MEDCouplingRefCountObject.cxx MEDCouplingPointSet.cxx MEDCouplingFieldTemplate.cxx diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index 693f24b30..d7cf86b66 100755 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -19,6 +19,7 @@ // Author : Anthony Geay (EDF R&D) #include "MEDCouplingFieldDiscretization.hxx" +#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx" #include "MEDCouplingCMesh.hxx" #include "MEDCouplingUMesh.hxx" #include "MEDCouplingFieldDouble.hxx" @@ -44,26 +45,16 @@ const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12; const char MEDCouplingFieldDiscretizationP0::REPR[]="P0"; -const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS; - const char MEDCouplingFieldDiscretizationP1::REPR[]="P1"; -const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES; - const mcIdType MEDCouplingFieldDiscretizationPerCell::DFT_INVALID_LOCID_VALUE=-1; const char MEDCouplingFieldDiscretizationGauss::REPR[]="GAUSS"; -const TypeOfField MEDCouplingFieldDiscretizationGauss::TYPE=ON_GAUSS_PT; - const char MEDCouplingFieldDiscretizationGaussNE::REPR[]="GSSNE"; -const TypeOfField MEDCouplingFieldDiscretizationGaussNE::TYPE=ON_GAUSS_NE; - const char MEDCouplingFieldDiscretizationKriging::REPR[]="KRIGING"; -const TypeOfField MEDCouplingFieldDiscretizationKriging::TYPE=ON_NODES_KR; - // doc is here http://www.code-aster.org/V2/doc/default/fr/man_r/r3/r3.01.01.pdf const double MEDCouplingFieldDiscretizationGaussNE::FGP_POINT1[1]={0.}; const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG2[2]={1.,1.}; @@ -142,6 +133,8 @@ MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField return new MEDCouplingFieldDiscretizationGaussNE; case MEDCouplingFieldDiscretizationKriging::TYPE: return new MEDCouplingFieldDiscretizationKriging; + case MEDCouplingFieldDiscretizationOnNodesFE::TYPE: + return new MEDCouplingFieldDiscretizationOnNodesFE; default: throw INTERP_KERNEL::Exception("Chosen discretization is not implemented yet."); } @@ -159,22 +152,30 @@ TypeOfField MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(const s return MEDCouplingFieldDiscretizationGaussNE::TYPE; if(repr==MEDCouplingFieldDiscretizationKriging::REPR) return MEDCouplingFieldDiscretizationKriging::TYPE; + if(repr==MEDCouplingFieldDiscretizationOnNodesFE::REPR) + return MEDCouplingFieldDiscretizationOnNodesFE::TYPE; throw INTERP_KERNEL::Exception("Representation does not match with any field discretization !"); } std::string MEDCouplingFieldDiscretization::GetTypeOfFieldRepr(TypeOfField type) { - if(type==MEDCouplingFieldDiscretizationP0::TYPE) - return MEDCouplingFieldDiscretizationP0::REPR; - if(type==MEDCouplingFieldDiscretizationP1::TYPE) - return MEDCouplingFieldDiscretizationP1::REPR; - if(type==MEDCouplingFieldDiscretizationGauss::TYPE) - return MEDCouplingFieldDiscretizationGauss::REPR; - if(type==MEDCouplingFieldDiscretizationGaussNE::TYPE) - return MEDCouplingFieldDiscretizationGaussNE::REPR; - if(type==MEDCouplingFieldDiscretizationKriging::TYPE) - return MEDCouplingFieldDiscretizationKriging::REPR; - throw INTERP_KERNEL::Exception("GetTypeOfFieldRepr : Representation does not match with any field discretization !"); + switch(type) + { + case MEDCouplingFieldDiscretizationP0::TYPE: + return MEDCouplingFieldDiscretizationP0::REPR; + case MEDCouplingFieldDiscretizationP1::TYPE: + return MEDCouplingFieldDiscretizationP1::REPR; + case MEDCouplingFieldDiscretizationGauss::TYPE: + return MEDCouplingFieldDiscretizationGauss::REPR; + case MEDCouplingFieldDiscretizationGaussNE::TYPE: + return MEDCouplingFieldDiscretizationGaussNE::REPR; + case MEDCouplingFieldDiscretizationKriging::TYPE: + return MEDCouplingFieldDiscretizationKriging::REPR; + case MEDCouplingFieldDiscretizationOnNodesFE::TYPE: + return MEDCouplingFieldDiscretizationOnNodesFE::REPR; + default: + throw INTERP_KERNEL::Exception("GetTypeOfFieldRepr : Representation does not match with any field discretization !"); + } } bool MEDCouplingFieldDiscretization::isEqual(const MEDCouplingFieldDiscretization *other, double eps) const @@ -480,20 +481,6 @@ void MEDCouplingFieldDiscretization::RenumberEntitiesFromN2OArr(const mcIdType * } } -template -MCAuto MEDCouplingFieldDiscretization::EasyAggregate(std::vector& fds) -{ - if(fds.empty()) - throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : input array is empty"); - for(const MEDCouplingFieldDiscretization * it : fds) - { - const FIELD_DISC *itc(dynamic_cast(it)); - if(!itc) - throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : same field discretization expected for all input discretizations !"); - } - return fds[0]->clone(); -} - MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization() { } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index e1004c571..2eab0a945 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -124,10 +124,10 @@ namespace MEDCoupling public: MEDCOUPLING_EXPORT TypeOfField getEnum() const; MEDCOUPLING_EXPORT std::string getClassName() const override { return std::string("MEDCouplingFieldDiscretizationP0"); } - MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const; + MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const override; MEDCOUPLING_EXPORT std::string getStringRepr() const; MEDCOUPLING_EXPORT const char *getRepr() const; - MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const; + MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const override; MEDCOUPLING_EXPORT mcIdType getNumberOfTuplesExpectedRegardingCode(const std::vector& code, const std::vector& idsPerType) const; MEDCOUPLING_EXPORT mcIdType getNumberOfTuples(const MEDCouplingMesh *mesh) const; MEDCOUPLING_EXPORT mcIdType getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const; @@ -135,14 +135,14 @@ namespace MEDCoupling MEDCOUPLING_EXPORT void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector& arrays, const mcIdType *old2NewBg, bool check); MEDCOUPLING_EXPORT DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const; - MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const; + MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const override; MEDCOUPLING_EXPORT void computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const mcIdType *tupleIdsBg, const mcIdType *tupleIdsEnd, DataArrayIdType *&cellRestriction, DataArrayIdType *&trueTupleRestriction) const; MEDCOUPLING_EXPORT void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const; - MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const; - MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; + MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override; + MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override; MEDCOUPLING_EXPORT void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, mcIdType i, mcIdType j, mcIdType k, double *res) const; - MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const; + MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override; MEDCOUPLING_EXPORT void renumberValuesOnNodes(double epsOnVals, const mcIdType *old2New, mcIdType newNbOfNodes, DataArrayDouble *arr) const; MEDCOUPLING_EXPORT void renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const mcIdType *old2New, mcIdType newSz, DataArrayDouble *arr) const; MEDCOUPLING_EXPORT void renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const mcIdType *new2old, mcIdType newSz, DataArrayDouble *arr) const; @@ -153,7 +153,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT void reprQuickOverview(std::ostream& stream) const; public: static const char REPR[]; - static const TypeOfField TYPE; + static constexpr TypeOfField TYPE = ON_CELLS; }; class MEDCouplingFieldDiscretizationOnNodes : public MEDCouplingFieldDiscretization @@ -184,19 +184,19 @@ namespace MEDCoupling public: MEDCOUPLING_EXPORT TypeOfField getEnum() const; MEDCOUPLING_EXPORT std::string getClassName() const override { return std::string("MEDCouplingFieldDiscretizationP1"); } - MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const; + MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const override; MEDCOUPLING_EXPORT std::string getStringRepr() const; MEDCOUPLING_EXPORT const char *getRepr() const; - MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const; - MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const; - MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const; - MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; - MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const; + MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const override; + MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const override; + MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override; + MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override; + MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override; MEDCOUPLING_EXPORT void reprQuickOverview(std::ostream& stream) const; MEDCOUPLING_EXPORT MCAuto aggregate(std::vector& fds) const override; public: static const char REPR[]; - static const TypeOfField TYPE; + static constexpr TypeOfField TYPE = ON_NODES; protected: MEDCOUPLING_EXPORT void getValueInCell(const MEDCouplingMesh *mesh, mcIdType cellId, const DataArrayDouble *arr, const double *loc, double *res) const; }; @@ -222,7 +222,7 @@ namespace MEDCoupling std::size_t getHeapMemorySizeWithoutChildren() const; std::vector getDirectChildrenWithNull() const; void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const; - bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const; + bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const override; bool isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const; void renumberCells(const mcIdType *old2NewBg, bool check); protected: @@ -238,9 +238,9 @@ namespace MEDCoupling MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationGauss(); MEDCOUPLING_EXPORT TypeOfField getEnum() const; MEDCOUPLING_EXPORT std::string getClassName() const override { return std::string("MEDCouplingFieldDiscretizationGauss"); } - MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const; + MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const override; MEDCOUPLING_EXPORT bool isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const; - MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const; + MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const override; MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clonePart(const mcIdType *startCellIds, const mcIdType *endCellIds) const; MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clonePartRange(mcIdType beginCellIds, mcIdType endCellIds, mcIdType stepCellIds) const; MEDCOUPLING_EXPORT std::string getStringRepr() const; @@ -255,7 +255,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const; MEDCOUPLING_EXPORT void computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const mcIdType *tupleIdsBg, const mcIdType *tupleIdsEnd, DataArrayIdType *&cellRestriction, DataArrayIdType *&trueTupleRestriction) const; - MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const; + MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const override; MEDCOUPLING_EXPORT void getTinySerializationIntInformation(std::vector& tinyInfo) const; MEDCOUPLING_EXPORT void getTinySerializationDbleInformation(std::vector& tinyInfo) const; MEDCOUPLING_EXPORT void finishUnserialization(const std::vector& tinyInfo); @@ -264,10 +264,10 @@ namespace MEDCoupling MEDCOUPLING_EXPORT void checkForUnserialization(const std::vector& tinyInfo, const DataArrayIdType *arr); MEDCOUPLING_EXPORT double getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da, mcIdType cellId, mcIdType nodeIdInCell, int compoId) const; MEDCOUPLING_EXPORT void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const; - MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const; - MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; + MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override; + MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override; MEDCOUPLING_EXPORT void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, mcIdType i, mcIdType j, mcIdType k, double *res) const; - MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const; + MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override; MEDCOUPLING_EXPORT MEDCouplingMesh *buildSubMeshData(const MEDCouplingMesh *mesh, const mcIdType *start, const mcIdType *end, DataArrayIdType *&di) const; MEDCOUPLING_EXPORT MEDCouplingMesh *buildSubMeshDataRange(const MEDCouplingMesh *mesh, mcIdType beginCellIds, mcIdType endCellIds, mcIdType stepCellIds, mcIdType& beginOut, mcIdType& endOut, mcIdType& stepOut, DataArrayIdType *&di) const; MEDCOUPLING_EXPORT DataArrayIdType *computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const mcIdType *startCellIds, const mcIdType *endCellIds) const; @@ -301,7 +301,7 @@ namespace MEDCoupling void commonUnserialization(const std::vector& tinyInfo); public: static const char REPR[]; - static const TypeOfField TYPE; + static constexpr TypeOfField TYPE = ON_GAUSS_PT; private: std::vector _loc; }; @@ -315,10 +315,10 @@ namespace MEDCoupling MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationGaussNE(); MEDCOUPLING_EXPORT TypeOfField getEnum() const; MEDCOUPLING_EXPORT std::string getClassName() const override { return std::string("MEDCouplingFieldDiscretizationGaussNE"); } - MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const; + MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const override; MEDCOUPLING_EXPORT std::string getStringRepr() const; MEDCOUPLING_EXPORT const char *getRepr() const; - MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const; + MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const override; MEDCOUPLING_EXPORT mcIdType getNumberOfTuplesExpectedRegardingCode(const std::vector& code, const std::vector& idsPerType) const; MEDCOUPLING_EXPORT mcIdType getNumberOfTuples(const MEDCouplingMesh *mesh) const; MEDCOUPLING_EXPORT mcIdType getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const; @@ -329,13 +329,13 @@ namespace MEDCoupling MEDCOUPLING_EXPORT void integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const; MEDCOUPLING_EXPORT void computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const mcIdType *tupleIdsBg, const mcIdType *tupleIdsEnd, DataArrayIdType *&cellRestriction, DataArrayIdType *&trueTupleRestriction) const; - MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const; + MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const override; MEDCOUPLING_EXPORT double getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da, mcIdType cellId, mcIdType nodeIdInCell, int compoId) const; MEDCOUPLING_EXPORT void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const; - MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const; - MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; + MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override; + MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override; MEDCOUPLING_EXPORT void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, mcIdType i, mcIdType j, mcIdType k, double *res) const; - MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const; + MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override; MEDCOUPLING_EXPORT MEDCouplingMesh *buildSubMeshData(const MEDCouplingMesh *mesh, const mcIdType *start, const mcIdType *end, DataArrayIdType *&di) const; MEDCOUPLING_EXPORT MEDCouplingMesh *buildSubMeshDataRange(const MEDCouplingMesh *mesh, mcIdType beginCellIds, mcIdType endCellIds, mcIdType stepCellIds, mcIdType& beginOut, mcIdType& endOut, mcIdType& stepOut, DataArrayIdType *&di) const; MEDCOUPLING_EXPORT DataArrayIdType *computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const mcIdType *startCellIds, const mcIdType *endCellIds) const; @@ -351,7 +351,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other); public: static const char REPR[]; - static const TypeOfField TYPE; + static constexpr TypeOfField TYPE = ON_GAUSS_NE; static const double FGP_POINT1[1]; static const double FGP_SEG2[2]; static const double FGP_SEG3[3]; @@ -418,13 +418,13 @@ namespace MEDCoupling MEDCOUPLING_EXPORT TypeOfField getEnum() const; MEDCOUPLING_EXPORT std::string getClassName() const override { return std::string("MEDCouplingFieldDiscretizationKriging"); } MEDCOUPLING_EXPORT const char *getRepr() const; - MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const; + MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const override; MEDCOUPLING_EXPORT std::string getStringRepr() const; - MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const; - MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const; - MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const; - MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; - MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const; + MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const override; + MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const override; + MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override; + MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override; + MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override; MEDCOUPLING_EXPORT void reprQuickOverview(std::ostream& stream) const; MEDCOUPLING_EXPORT MCAuto aggregate(std::vector& fds) const override; public://specific part @@ -440,8 +440,22 @@ namespace MEDCoupling MEDCOUPLING_EXPORT static DataArrayDouble *PerformDriftOfVec(const DataArrayDouble *arr, mcIdType isDrift); public: static const char REPR[]; - static const TypeOfField TYPE; + static constexpr TypeOfField TYPE = ON_NODES_KR; }; + + template + MCAuto MEDCouplingFieldDiscretization::EasyAggregate(std::vector& fds) + { + if(fds.empty()) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : input array is empty"); + for(const MEDCouplingFieldDiscretization * it : fds) + { + const FIELD_DISC *itc(dynamic_cast(it)); + if(!itc) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : same field discretization expected for all input discretizations !"); + } + return fds[0]->clone(); + } } #endif diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx new file mode 100755 index 000000000..fe696ab81 --- /dev/null +++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx @@ -0,0 +1,274 @@ +// Copyright (C) 2007-2022 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, 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 +// 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 +// +// Author : Anthony Geay (EDF R&D) + +#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx" +#include "MEDCouplingNormalizedUnstructuredMesh.txx" +#include "InterpKernelDenseMatrix.hxx" +#include "InterpKernelRootsMultiDim.hxx" +#include "MEDCouplingUMesh.hxx" +#include "InterpolationHelper.txx" +#include "InterpKernelGaussCoords.hxx" + +#include + +using namespace MEDCoupling; + +const char MEDCouplingFieldDiscretizationOnNodesFE::REPR[]="FE"; + +std::string MEDCouplingFieldDiscretizationOnNodesFE::getStringRepr() const +{ + return std::string(REPR); +} + +void MEDCouplingFieldDiscretizationOnNodesFE::reprQuickOverview(std::ostream& stream) const +{ + stream << "NodeFE spatial discretization."; +} + +MCAuto MEDCouplingFieldDiscretizationOnNodesFE::aggregate(std::vector& fds) const +{ + return EasyAggregate(fds); +} + +bool MEDCouplingFieldDiscretizationOnNodesFE::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const +{ + if(!other) + { + reason="other spatial discretization is NULL, and this spatial discretization (Node FE) is defined."; + return false; + } + const MEDCouplingFieldDiscretizationOnNodesFE *otherC=dynamic_cast(other); + bool ret=otherC!=0; + if(!ret) + reason="Spatial discrtization of this is ON_NODES_FE, which is not the case of other."; + return ret; +} + +/*! + * This method is simply called by MEDCouplingFieldDiscretization::deepCopy. It performs the deep copy of \a this. + * + * \sa MEDCouplingFieldDiscretization::deepCopy. + */ +MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationOnNodesFE::clone() const +{ + return new MEDCouplingFieldDiscretizationOnNodesFE; +} + +void MEDCouplingFieldDiscretizationOnNodesFE::checkCompatibilityWithNature(NatureOfField nat) const +{ + if(nat!=IntensiveMaximum) + throw INTERP_KERNEL::Exception("Invalid nature for NodeFE field : expected IntensiveMaximum !"); +} + +MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const +{ + if(!mesh) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField : mesh instance specified is NULL !"); + throw INTERP_KERNEL::Exception("getMeasureField on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !"); +} + +class Functor +{ +private: + const MEDCouplingGaussLocalization* _gl; + std::size_t _nb_pts_in_cell; + const double *_pts_in_cell; + const double *_point; +public: + Functor(const MEDCouplingGaussLocalization& gl, std::size_t nbPtsInCell, const double *ptsInCell, const double point[3]):_gl(&gl),_nb_pts_in_cell(nbPtsInCell), + _pts_in_cell(ptsInCell),_point(point) { } + std::vector operator()(const std::vector& x) + { + MEDCouplingGaussLocalization gl(_gl->getType(),_gl->getRefCoords(),x,{1.0}); + MCAuto shapeFunc = gl.getShapeFunctionValues(); + const double *shapeFuncPtr( shapeFunc->begin() ); + std::vector ret(3,0); + for(std::size_t iPt = 0; iPt < _nb_pts_in_cell ; ++iPt) + { + for(short iDim = 0 ; iDim < 3 ; ++iDim) + ret[iDim] += shapeFuncPtr[iPt] * _pts_in_cell[3*iPt + iDim]; + } + ret[0] -= _point[0]; ret[1] -= _point[1]; ret[2] -= _point[2]; + return ret; + } +}; + +bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector& ptsInCell, const double locInReal[3], double locInRef[3]) +{ + constexpr double EPS_IN_OUT = 1e-12; + std::size_t nbPtsInCell(ptsInCell.size()/3); + bool ret(false); + const double *refCoo(gl.getRefCoords().data()); + INTERP_KERNEL::NormalizedCellType ct(gl.getType()); + Functor func(gl,nbPtsInCell,ptsInCell.data(),locInReal); + + auto myJacobian = [&gl,nbPtsInCell,ptsInCell](const std::vector& x, const std::vector&, INTERP_KERNEL::DenseMatrix& jacobian) + { + MEDCouplingGaussLocalization mygl(gl.getType(),gl.getRefCoords(),x,{1.0}); + MCAuto shapeFunc = mygl.getDerivativeOfShapeFunctionValues(); + for(std::size_t i = 0 ; i < 3 ; ++i) + for(std::size_t j = 0 ; j < 3 ; ++j) + { + double res = 0.0; + for( std::size_t k = 0 ; k < nbPtsInCell ; ++k ) + res += ptsInCell[k*3+i] * shapeFunc->getIJ(0,3*k+j); + jacobian[ i ][ j ] = res; + } + }; + + // loop on refcoords as initialization point for Newton algo. vini is the initialization vector of Newton. + for(std::size_t attemptId = 0 ; attemptId < nbPtsInCell ; ++attemptId) + { + std::vector vini(refCoo + attemptId*3, refCoo + (attemptId+1)*3); + try + { + bool check(true); + //INTERP_KERNEL::SolveWithNewton(vini,check,func); + INTERP_KERNEL::SolveWithNewtonWithJacobian(vini,check,func,myJacobian); + ret = (check==false);//looks strange but OK regarding newt (SolveWithNewton) at page 387 of numerical recipes for semantic of check parameter + } + catch( INTERP_KERNEL::Exception& ex ) + { ret = false; }// Something get wrong during Newton process + if(ret) + {//Newton has converged. Now check if it converged to a point inside cell + if( ! INTERP_KERNEL::GaussInfo::IsInOrOutForReference(ct,vini.data(),EPS_IN_OUT) ) + {// converged but locInReal has been detected outside of cell + ret = false; + } + } + if(ret) + { + locInRef[0] = vini[0]; locInRef[1] = vini[1]; locInRef[2] = vini[2]; + return ret; + } + } + std::fill(locInRef,locInRef+3,std::numeric_limits::max()); + return false; +} + +void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *ptsCoo, double *res) const +{ + MCAuto res2( this->getValueOnMulti(arr,mesh,ptsCoo,1) ); + std::copy(res2->begin(),res2->end(),res); +} + +void MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts, + std::function&)> customFunc) +{ + const double *coordsOfMesh( umesh->getCoords()->begin() ); + MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh); + BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) ); + for(mcIdType iPt = 0 ; iPt < nbOfPts ; ++iPt) + { + std::vector elems; + tree.getElementsAroundPoint(ptsCoo+3*iPt,elems); + bool found(false); + for(auto cellId = elems.cbegin() ; cellId != elems.cend() && !found ; ++cellId) + { + INTERP_KERNEL::NormalizedCellType gt( umesh->getTypeOfCell(*cellId) ); + std::vector conn; + umesh->getNodeIdsOfCell(*cellId,conn); + MCAuto refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) ); + std::vector refCooCpp(refCoo->begin(),refCoo->end()); + std::vector gsCoo(ptsCoo + iPt*3,ptsCoo + (iPt+1)*3); + MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.}); + std::vector ptsInCell; ptsInCell.reserve(conn.size()*gl.getDimension()); + std::for_each( conn.cbegin(), conn.cend(), [coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*3,coordsOfMesh+(c+1)*3); } ); + std::vector locInRef(3); + if( IsInside3D(gl,ptsInCell,gsCoo.data(),locInRef.data()) ) + { + gl.setGaussCoords(locInRef); + customFunc(gl,conn); + found = true; + } + } + if(!found) + THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << ptsCoo[0] << " Y=" << ptsCoo[1] << " Z=" << ptsCoo[2] << " !"); + } +} + +const MEDCouplingUMesh *MEDCouplingFieldDiscretizationOnNodesFE::checkConfig3D(const MEDCouplingMesh *mesh) const +{ + const MEDCouplingUMesh *umesh( dynamic_cast(mesh) ); + if( !umesh ) + THROW_IK_EXCEPTION("getValueOn : not implemented yet for type != MEDCouplingUMesh !"); + if(umesh->getSpaceDimension() != 3 || umesh->getMeshDimension() != 3) + THROW_IK_EXCEPTION("getValueOn : implemented only for meshes with spacedim == 3 and meshdim == 3 !"); + umesh->checkConsistency(); + return umesh; +} + +DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const +{ + if(!arr || !arr->isAllocated()) + throw INTERP_KERNEL::Exception("getValueOnMulti : input array is null or not allocated !"); + mcIdType nbOfRows=getNumberOfMeshPlaces(mesh); + if(arr->getNumberOfTuples()!=nbOfRows) + { + THROW_IK_EXCEPTION( "getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !") + } + const MEDCouplingUMesh *umesh = checkConfig3D(mesh); + std::size_t nbCompo( arr->getNumberOfComponents() ); + MCAuto ret(DataArrayDouble::New()); + ret->alloc(nbOfTargetPoints,nbCompo); + double *res( ret->getPointer() ); + + auto arrayFeeder = [arr,&res,nbCompo](const MEDCouplingGaussLocalization& gl, const std::vector& conn) + { + MCAuto resVector( gl.getShapeFunctionValues() ); + { + std::for_each(res,res+nbCompo,[](double& v) { v = 0.0; }); + for(std::size_t iComp = 0 ; iComp < nbCompo ; ++iComp) + for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt) + { + { + res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp); + } + } + res += nbCompo; + } + }; + + GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfTargetPoints,arrayFeeder); + return ret.retn(); +} + +/*! + * Returns for each \a nbOfPoints point in \a loc its coordinate in reference element. + */ +MCAuto MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement(const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const +{ + const MEDCouplingUMesh *umesh = checkConfig3D(mesh); + MCAuto ret(DataArrayDouble::New()); + ret->alloc(nbOfPoints,3); + double *retPtr(ret->getPointer() ); + + auto arrayFeeder = [&retPtr](const MEDCouplingGaussLocalization& gl, const std::vector& conn) + { + std::vector resVector( gl.getGaussCoords() ); + { + std::copy(resVector.begin(),resVector.end(),retPtr); + retPtr += 3; + } + }; + + GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfPoints,arrayFeeder); + return ret; +} diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx new file mode 100644 index 000000000..378d0964c --- /dev/null +++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx @@ -0,0 +1,58 @@ +// Copyright (C) 2022 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, 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 +// 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 +// +// Author : Anthony Geay (EDF R&D) + +#pragma once + +#include "MEDCouplingFieldDiscretization.hxx" + +#include + +namespace MEDCoupling +{ + /*! + * Class in charge to implement FE functions with shape functions + */ + class MEDCouplingFieldDiscretizationOnNodesFE : public MEDCouplingFieldDiscretizationOnNodes + { + public: + MEDCOUPLING_EXPORT TypeOfField getEnum() const override { return TYPE; } + MEDCOUPLING_EXPORT std::string getClassName() const override { return std::string("MEDCouplingFieldDiscretizationOnNodesFE"); } + MEDCOUPLING_EXPORT const char *getRepr() const override { return REPR; } + MEDCOUPLING_EXPORT std::string getStringRepr() const override; + MEDCOUPLING_EXPORT void reprQuickOverview(std::ostream& stream) const override; + MEDCOUPLING_EXPORT MCAuto aggregate(std::vector& fds) const override; + MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const override; + MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization *clone() const override; + MEDCOUPLING_EXPORT void checkCompatibilityWithNature(NatureOfField nat) const override; + MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override; + MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override; + MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override; + public: + MEDCOUPLING_EXPORT MCAuto getCooInRefElement(const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const; + public: + MEDCOUPLING_EXPORT static void GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts, + std::function&)> customFunc); + private: + const MEDCouplingUMesh *checkConfig3D(const MEDCouplingMesh *mesh) const; + public: + static const char REPR[]; + static constexpr TypeOfField TYPE = ON_NODES_FE; + }; +} diff --git a/src/MEDCoupling/MEDCouplingFieldT.txx b/src/MEDCoupling/MEDCouplingFieldT.txx index 125c314b4..b4c13257f 100644 --- a/src/MEDCoupling/MEDCouplingFieldT.txx +++ b/src/MEDCoupling/MEDCouplingFieldT.txx @@ -449,7 +449,8 @@ namespace MEDCoupling * \ref MEDCoupling::ON_NODES "ON_NODES", * \ref MEDCoupling::ON_GAUSS_PT "ON_GAUSS_PT", * \ref MEDCoupling::ON_GAUSS_NE "ON_GAUSS_NE", - * \ref MEDCoupling::ON_NODES_KR "ON_NODES_KR"). + * \ref MEDCoupling::ON_NODES_KR "ON_NODES_KR", + * \ref MEDCoupling::ON_NODES_FE "ON_NODES_FE"). * * For example, \a this is a field on cells lying on a mesh that have 10 cells, \a partBg contains the following cell ids [3,7,6]. * Then the returned field will lie on mesh having 3 cells and will contain 3 tuples. @@ -512,7 +513,8 @@ namespace MEDCoupling * \ref MEDCoupling::ON_NODES "ON_NODES", * \ref MEDCoupling::ON_GAUSS_PT "ON_GAUSS_PT", * \ref MEDCoupling::ON_GAUSS_NE "ON_GAUSS_NE", - * \ref MEDCoupling::ON_NODES_KR "ON_NODES_KR"). + * \ref MEDCoupling::ON_NODES_KR "ON_NODES_KR", + * \ref MEDCoupling::ON_NODES_FE "ON_NODES_FE"). * * For example, \a this is a field on cells lying on a mesh that have 10 cells, \a part contains following cell ids [3,7,6]. * Then the returned field will lie on mesh having 3 cells and the returned field will contain 3 tuples.
diff --git a/src/MEDCoupling/MEDCouplingGaussLocalization.cxx b/src/MEDCoupling/MEDCouplingGaussLocalization.cxx index 3026d8ac1..e59f54741 100644 --- a/src/MEDCoupling/MEDCouplingGaussLocalization.cxx +++ b/src/MEDCoupling/MEDCouplingGaussLocalization.cxx @@ -76,15 +76,20 @@ void MEDCouplingGaussLocalization::checkConsistencyLight() const int MEDCouplingGaussLocalization::getDimension() const { if(_weight.empty()) - return -1; + THROW_IK_EXCEPTION("getDimension : weight is empty !"); return (int)_gauss_coord.size()/(int)_weight.size(); } int MEDCouplingGaussLocalization::getNumberOfPtsInRefCell() const { - int dim=getDimension(); - if(dim==0) - return -1; + if(_gauss_coord.empty()) + { + if(!_weight.empty()) + THROW_IK_EXCEPTION("getNumberOfPtsInRefCell : gauss_coords are empty whereas weights are not empty !"); + const INTERP_KERNEL::CellModel& cm = INTERP_KERNEL::CellModel::GetCellModel(_type); + return ((int)_ref_coord.size()) / ((int)cm.getDimension()); + } + int dim( getDimension() ); return (int)_ref_coord.size()/dim; } @@ -241,6 +246,46 @@ MCAuto MEDCouplingGaussLocalization::buildRefCell() const return MCAuto(ret->buildUnstructured()); } +/*! + * This method returns an array containing for each Gauss Points in \a this, function values relative to the points of the + * reference cell. Number of components of returned array is equal to the number of points of the reference cell. + */ +MCAuto MEDCouplingGaussLocalization::getShapeFunctionValues() const +{ + MCAuto ret(DataArrayDouble::New()); + int nbGaussPt(getNumberOfGaussPt()),nbPtsRefCell(getNumberOfPtsInRefCell()),dim(getDimension()); + ret->alloc(nbGaussPt,nbPtsRefCell); + double *retPtr(ret->getPointer()); + for(int iGaussPt = 0 ; iGaussPt < nbGaussPt ; ++iGaussPt) + { + std::vector curGaussPt(_gauss_coord.begin()+iGaussPt*dim,_gauss_coord.begin()+(iGaussPt+1)*dim); + INTERP_KERNEL::GaussInfo gi(_type,curGaussPt,1,_ref_coord,nbPtsRefCell); + gi.initLocalInfo(); + const double *funcVal( gi.getFunctionValues(0) ); + std::copy(funcVal,funcVal+nbPtsRefCell,retPtr); + retPtr += nbPtsRefCell; + } + return ret; +} + +MCAuto MEDCouplingGaussLocalization::getDerivativeOfShapeFunctionValues() const +{ + MCAuto ret(DataArrayDouble::New()); + int nbGaussPt(getNumberOfGaussPt()),nbPtsRefCell(getNumberOfPtsInRefCell()),dim(getDimension()); + ret->alloc(nbGaussPt,nbPtsRefCell*dim); + double *retPtr(ret->getPointer()); + for(int iGaussPt = 0 ; iGaussPt < nbGaussPt ; ++iGaussPt) + { + std::vector curGaussPt(_gauss_coord.begin()+iGaussPt*dim,_gauss_coord.begin()+(iGaussPt+1)*dim); + INTERP_KERNEL::GaussInfo gi(_type,curGaussPt,1,_ref_coord,nbPtsRefCell); + gi.initLocalInfo(); + const double *devOfFuncVal( gi.getDerivativeOfShapeFunctionAt(0) ); + std::copy(devOfFuncVal,devOfFuncVal+nbPtsRefCell*dim,retPtr); + retPtr += nbPtsRefCell*dim; + } + return ret; +} + /*! * This method sets the comp_th component of ptIdInCell_th point coordinate of reference element of type this->_type. * @throw if not 0<=ptIdInCell& v1, std::transform(tmp.begin(),tmp.end(),tmp.begin(),[](double c){return fabs(c);}); return *std::max_element(tmp.begin(),tmp.end()) MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type) +{ + std::vector retCpp(INTERP_KERNEL::GaussInfo::GetDefaultReferenceCoordinatesOf(type)); + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type); + auto nbDim(cm.getDimension()); + std::size_t sz(retCpp.size()); + MCAuto ret(DataArrayDouble::New()); + if( sz%std::size_t(nbDim) != 0 ) + THROW_IK_EXCEPTION("GetDefaultReferenceCoordinatesOf : unexpected size of defaut array : " << sz << " % " << nbDim << " != 0 !"); + ret->alloc(sz/size_t(nbDim),nbDim); + std::copy(retCpp.begin(),retCpp.end(),ret->getPointer()); + return ret; +} diff --git a/src/MEDCoupling/MEDCouplingGaussLocalization.hxx b/src/MEDCoupling/MEDCouplingGaussLocalization.hxx index 56e061e93..49cb19da3 100644 --- a/src/MEDCoupling/MEDCouplingGaussLocalization.hxx +++ b/src/MEDCoupling/MEDCouplingGaussLocalization.hxx @@ -54,6 +54,8 @@ namespace MEDCoupling // MEDCOUPLING_EXPORT MCAuto localizePtsInRefCooForEachCell(const DataArrayDouble *ptsInRefCoo, const MEDCouplingUMesh *mesh) const; MEDCOUPLING_EXPORT MCAuto buildRefCell() const; + MEDCOUPLING_EXPORT MCAuto getShapeFunctionValues() const; + MEDCOUPLING_EXPORT MCAuto getDerivativeOfShapeFunctionValues() const; // MEDCOUPLING_EXPORT const std::vector& getRefCoords() const { return _ref_coord; } MEDCOUPLING_EXPORT double getRefCoord(int ptIdInCell, int comp) const; @@ -70,6 +72,7 @@ namespace MEDCoupling // MEDCOUPLING_EXPORT static MEDCouplingGaussLocalization BuildNewInstanceFromTinyInfo(mcIdType dim, const std::vector& tinyData); MEDCOUPLING_EXPORT static bool AreAlmostEqual(const std::vector& v1, const std::vector& v2, double eps); + MEDCOUPLING_EXPORT static MCAuto GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type); private: int checkCoherencyOfRequest(mcIdType gaussPtIdInCell, int comp) const; private: diff --git a/src/MEDCoupling/MEDCouplingMemArray.txx b/src/MEDCoupling/MEDCouplingMemArray.txx index 9cb18e6fc..16a07683f 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.txx +++ b/src/MEDCoupling/MEDCouplingMemArray.txx @@ -5866,7 +5866,7 @@ struct NotInRange throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : only single component allowed !"); std::size_t nbOfElements=this->getNumberOfTuples(); if(nbOfElements<2) - throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : 1 tuple at least must be present in 'this' !"); + throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : 2 tuples at least must be present in 'this' !"); const T *ptr=this->getConstPointer(); DataArrayType *ret=DataArrayType::New(); ret->alloc(nbOfElements-1,1); diff --git a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx index d3b7ff2a1..7934d0562 100644 --- a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx +++ b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx @@ -35,7 +35,7 @@ class MEDCouplingNormalizedUnstructuredMesh public: static const int MY_SPACEDIM=SPACEDIM; static const int MY_MESHDIM=MESHDIM; - typedef mcIdType MyConnType; + using MyConnType = mcIdType; static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE; public: MEDCouplingNormalizedUnstructuredMesh(const MEDCoupling::MEDCouplingPointSet *mesh); diff --git a/src/MEDCoupling/MEDCouplingRefCountObject.hxx b/src/MEDCoupling/MEDCouplingRefCountObject.hxx index 674eb9bc5..6dd1ee5c2 100644 --- a/src/MEDCoupling/MEDCouplingRefCountObject.hxx +++ b/src/MEDCoupling/MEDCouplingRefCountObject.hxx @@ -45,7 +45,8 @@ namespace MEDCoupling ON_NODES = 1, ON_GAUSS_PT = 2, ON_GAUSS_NE = 3, - ON_NODES_KR = 4 + ON_NODES_KR = 4, + ON_NODES_FE = 5 } TypeOfField; //! The various temporal discretization of a field diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index 8f620ffe6..231f29bba 100755 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -27,6 +27,7 @@ #include "MEDCouplingCMesh.hxx" #include "MEDCouplingNormalizedUnstructuredMesh.txx" #include "MEDCouplingNormalizedCartesianMesh.txx" +#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx" #include "Interpolation1D.txx" #include "Interpolation2DCurve.hxx" @@ -123,10 +124,7 @@ void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, } } _matrix=m; - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(srcNbElem); + synchronizeSizeOfSideMatricesAfterMatrixComputation(srcNbElem); } int MEDCouplingRemapper::prepareInterpKernelOnly() @@ -181,6 +179,8 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnly() { case 0: return prepareNotInterpKernelOnlyGaussGauss(); + case 1: + return prepareNotInterpKernelOnlyFEFE(); default: { std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !"; @@ -689,11 +689,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() } else throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension"); - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(nbCols); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols); return 1; } @@ -727,11 +723,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyEE() buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D, target_mesh->getMesh3DIds()->getConstPointer()); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(nbCols2D*nbCols1D); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols2D*nbCols1D); return 1; } @@ -783,11 +775,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUC() ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix); nullifiedTinyCoeffInCrudeMatrixAbs(0.); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(src_mesh->getNumberOfCells()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells()); return 1; } @@ -837,11 +825,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCU() } nullifiedTinyCoeffInCrudeMatrixAbs(0.); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(src_mesh->getNumberOfCells()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells()); return 1; } @@ -890,11 +874,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCC() } nullifiedTinyCoeffInCrudeMatrixAbs(0.); // - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(src_mesh->getNumberOfCells()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells()); return 1; } @@ -949,11 +929,43 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss() for(const mcIdType *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++) _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.; } - _deno_multiply.clear(); - _deno_multiply.resize(_matrix.size()); - _deno_reverse_multiply.clear(); - _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples()); - declareAsNew(); + synchronizeSizeOfSideMatricesAfterMatrixComputation( srcLoc->getNumberOfTuples() ); + return 1; +} + +int MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE() +{ + if(getIntersectionType()!=INTERP_KERNEL::PointLocator) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE : The intersection type is not supported ! Only PointLocator is supported for FE->FE interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !"); + MCAuto trgLoc=_target_ft->getLocalizationOfDiscr(); + mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents()); + if(trgSpaceDim!=3) + THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for target !") + if(_src_ft->getMesh()->getSpaceDimension() != 3) + THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for source !") + const MEDCouplingUMesh *srcUMesh( dynamic_cast(_src_ft->getMesh()) ); + const MEDCouplingPointSet *trgMesh( dynamic_cast(_target_ft->getMesh()) ); + if( !srcUMesh ) + THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D UMesh supported as source !"); + if( !trgMesh ) + THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D PointSet mesh supported as target !"); + + _matrix.clear(); + _matrix.resize(trgMesh->getNumberOfNodes()); + mcIdType rowId(0); + + auto matrixFeeder = [this,&rowId](const MEDCouplingGaussLocalization& gl, const std::vector& conn) + { + auto& row = this->_matrix[rowId++]; + MCAuto resVector( gl.getShapeFunctionValues() ); + for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt) + { + row[ conn[iPt] ] = resVector->getIJ(0,iPt); + } + }; + + MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(srcUMesh,trgMesh->getCoords()->begin(),trgMesh->getNumberOfNodes(),matrixFeeder); + synchronizeSizeOfSideMatricesAfterMatrixComputation( srcUMesh->getNumberOfNodes() ); return 1; } @@ -965,9 +977,11 @@ int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel { if(method=="GAUSSGAUSS") return 0; + if(method=="FEFE") + return 1; std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : "; oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method."; - oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !"; + oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS FEFE !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } @@ -1029,6 +1043,15 @@ void MEDCouplingRemapper::checkPrepare() const throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !"); } +void MEDCouplingRemapper::synchronizeSizeOfSideMatricesAfterMatrixComputation(mcIdType nbOfColsInMatrix) +{ + _deno_multiply.clear(); + _deno_multiply.resize(_matrix.size()); + _deno_reverse_multiply.clear(); + _deno_reverse_multiply.resize(nbOfColsInMatrix); + declareAsNew(); +} + /*! * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft. * This method returns 3 information (2 in output parameters and 1 in return). diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx index 6966f6575..23bf893b3 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.hxx +++ b/src/MEDCoupling/MEDCouplingRemapper.hxx @@ -89,12 +89,14 @@ namespace MEDCoupling // int prepareNotInterpKernelOnly(); int prepareNotInterpKernelOnlyGaussGauss(); + int prepareNotInterpKernelOnlyFEFE(); // static int CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method); // bool isInterpKernelOnlyOrNotOnly() const; void updateTime() const; void checkPrepare() const; + void synchronizeSizeOfSideMatricesAfterMatrixComputation(mcIdType nbOfColsInMatrix); std::string checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const; void releaseData(bool matrixSuppression); void restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index 58979fbdc..bd3edc096 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1029,6 +1029,57 @@ class MEDCouplingBasicsTest7(unittest.TestCase): self.assertTrue( all([len(elt)==1 for elt in res]) ) self.assertTrue( all([elt[0]>0.99 and elt[0]<1.01 for elt in res]) ) + @unittest.skipUnless(MEDCouplingHasNumPyBindings(),"requires numpy") + def testShapeFuncAndDerivative0(self): + """ + Test values returned by MEDCoupling on HEXA27 element of shape function and its derivatives. + See https://www.code-aster.org/V2/doc/v10/fr/man_r/r3/r3.01.01.pdf + """ + import numpy as np + + ref_coords_hexa27_med = [[-1.0, -1.0, -1.0], [-1.0, 1.0, -1.0], [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [-1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [1.0, 0.0, -1.0], [0.0, -1.0, -1.0], [-1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [1.0, 0.0, 1.0], [0.0, -1.0, 1.0], [-1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]] + + def coor2index(coor): + zeMap = {-1.0 : 0, 0.0 : 2 , 1.0 : 1} + return zeMap[coor] + + vcoor2index = np.vectorize( coor2index ) + node2ijk_hexa27_med = vcoor2index( np.array(ref_coords_hexa27_med) ) + + def N_1d_quad(x): + return np.array([-0.5*x*(1-x), 0.5*x*(x+1), 1.-x*x]) + + def N_3d_hexa27(x, i, j, k): + return N_1d_quad(x[0])[i]*N_1d_quad(x[1])[j]*N_1d_quad(x[2])[k] + + def N_hexa27(x): + return np.array([N_3d_hexa27(x, *node2ijk_hexa27_med[node,:]) for node in range(27)]) + + # Implementing shape function derivatives + def diff_N_1d_quad(x): + return np.array([x-0.5, x+0.5, -2.*x]) + + def diff_N_3d_hexa27(x, i, j, k): + return np.array([diff_N_1d_quad(x[0])[i]*N_1d_quad(x[1])[j] *N_1d_quad(x[2])[k], + N_1d_quad(x[0])[i] *diff_N_1d_quad(x[1])[j]*N_1d_quad(x[2])[k], + N_1d_quad(x[0])[i] *N_1d_quad(x[1])[j] *diff_N_1d_quad(x[2])[k]]) + + def diff_N_hexa27(x): + return np.array([diff_N_3d_hexa27(x, *node2ijk_hexa27_med[node,:]) for node in range(27)]) + # computation of ref values + posInRefCoord = [-0.85685375,-0.90643355,-0.90796825] + ref = N_hexa27( np.array(posInRefCoord) ) + ref2 = diff_N_hexa27( np.array(posInRefCoord) ) + # computation using MEDCoupling + gl = MEDCouplingGaussLocalization(NORM_HEXA27,sum(ref_coords_hexa27_med,[]),posInRefCoord,[1]) + mcShapeFunc = gl.getShapeFunctionValues() + mcShapeFunc.rearrange(1) + self.assertTrue( mcShapeFunc.isEqual(DataArrayDouble(ref),1e-12) ) + + mvDevOfShapeFunc = gl.getDerivativeOfShapeFunctionValues() + mvDevOfShapeFunc.rearrange(1) + ref2_mc = DataArrayDouble(ref2) ; ref2_mc.rearrange(1) + self.assertTrue( mvDevOfShapeFunc.isEqual( ref2_mc, 1e-12) ) pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index c760c94d2..7acf25f32 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -47,6 +47,7 @@ #include "MEDCouplingFieldOverTime.hxx" #include "MEDCouplingDefinitionTime.hxx" #include "MEDCouplingFieldDiscretization.hxx" +#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx" #include "MEDCouplingCartesianAMRMesh.hxx" #include "MEDCouplingAMRAttribute.hxx" #include "MEDCouplingMatrix.hxx" @@ -217,6 +218,7 @@ typedef long mcPyPtrType; %feature("autodoc", "1"); %feature("docstring"); +%newobject MEDCoupling::MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement; %newobject MEDCoupling::MEDCouplingField::buildMeasureField; %newobject MEDCoupling::MEDCouplingField::getLocalizationOfDiscr; %newobject MEDCoupling::MEDCouplingField::computeTupleIdsToSelectFromCellIds; @@ -478,6 +480,9 @@ typedef long mcPyPtrType; %newobject MEDCoupling::DenseMatrix::__mul__; %newobject MEDCoupling::MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell; %newobject MEDCoupling::MEDCouplingGaussLocalization::buildRefCell; +%newobject MEDCoupling::MEDCouplingGaussLocalization::getShapeFunctionValues; +%newobject MEDCoupling::MEDCouplingGaussLocalization::getDerivativeOfShapeFunctionValues; +%newobject MEDCoupling::MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf; %newobject MEDCoupling::MEDCouplingSkyLineArray::BuildFromPolyhedronConn; %newobject MEDCoupling::MEDCouplingSkyLineArray::getSuperIndexArray; %newobject MEDCoupling::MEDCouplingSkyLineArray::getIndexArray; @@ -503,6 +508,7 @@ typedef long mcPyPtrType; %feature("unref") MEDCouplingFieldDiscretizationGauss "$this->decrRef();" %feature("unref") MEDCouplingFieldDiscretizationGaussNE "$this->decrRef();" %feature("unref") MEDCouplingFieldDiscretizationKriging "$this->decrRef();" +%feature("unref") MEDCouplingFieldDiscretizationOnNodesFE "$this->decrRef();" %feature("unref") MEDCouplingFieldDouble "$this->decrRef();" %feature("unref") MEDCouplingFieldFloat "$this->decrRef();" %feature("unref") MEDCouplingFieldInt32 "$this->decrRef();" @@ -636,7 +642,8 @@ namespace MEDCoupling ON_NODES = 1, ON_GAUSS_PT = 2, ON_GAUSS_NE = 3, - ON_NODES_KR = 4 + ON_NODES_KR = 4, + ON_NODES_FE = 5 } TypeOfField; typedef enum @@ -1304,6 +1311,24 @@ namespace MEDCoupling MCAuto ret(self->buildRefCell()); return ret.retn(); } + + DataArrayDouble *getShapeFunctionValues() const + { + MCAuto ret(self->getShapeFunctionValues()); + return ret.retn(); + } + + DataArrayDouble *getDerivativeOfShapeFunctionValues() const + { + MCAuto ret(self->getDerivativeOfShapeFunctionValues()); + return ret.retn(); + } + + static DataArrayDouble *GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type) + { + MCAuto ret(MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(type)); + return ret.retn(); + } } }; diff --git a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i index 1ffd26910..deb4bee70 100644 --- a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i +++ b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i @@ -469,4 +469,21 @@ namespace MEDCoupling } } }; + + class MEDCouplingFieldDiscretizationOnNodesFE : public MEDCouplingFieldDiscretizationOnNodes + { + public: + %extend + { + DataArrayDouble *getCooInRefElement(const MEDCouplingMesh *mesh, PyObject *locs) + { + mcIdType sw,nbPts; + double v0; MEDCoupling::DataArrayDouble *v1(0); MEDCoupling::DataArrayDoubleTuple *v2(0); std::vector v3; + const double *inp=convertObjToPossibleCpp5_Safe2(locs,sw,v0,v1,v2,v3,"wrap of MEDCouplingFieldDouble::getValueOnMulti", + mesh->getSpaceDimension(),true,nbPts); + MCAuto ret(self->getCooInRefElement(mesh,inp,nbPts)); + return ret.retn(); + } + } + }; } diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index cbb2c2f6d..ef5df95cd 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -1586,6 +1586,62 @@ class MEDCouplingBasicsTest(unittest.TestCase): delta = abs(csr_new[0,0]-ref_value)/ref_value self.assertTrue(delta < 1e-3) + def testFEFE_0(self): + """ + Test to stress localisation of target points into a source mesh using standard reference FE elements. + """ + gt = NORM_HEXA27 + nbPtsInCell = MEDCouplingUMesh.GetNumberOfNodesOfGeometricType( gt ) + coo = DataArrayDouble([(9.0, 18.0, 27.0), (9.0, 22.0, 27.0), (11.0, 22.0, 27.0), (11.0, 18.0, 27.0), (9.0, 18.0, 33.0), (9.0, 22.0, 33.0), (11.0, 22.0, 33.0), (11.0, 18.0, 33.0), (8.8, 20.0, 26.4), (10.0, 21.6, 27.6), (11.2, 20.0, 26.4), (10.0, 18.4, 27.6), (8.8, 20.0, 33.6), (10.0, 21.6, 32.4), (11.2, 20.0, 33.6), (10.0, 18.4, 32.4), (8.8, 17.6, 30.0), (9.2, 21.6, 30.0), (11.2, 22.4, 30.0), (10.8, 18.4, 30.0), (10.0, 20.0, 26.4), (9.2, 20.0, 30.0), (10.0, 22.4, 30.0), (10.8, 20.0, 30.0), (10.0, 17.6, 30.0), (10.0, 20.0, 32.4), (10.0, 20.0, 30.0)]) + m = MEDCouplingUMesh("mesh",3) + m.setCoords(coo) + m.allocateCells() + m.insertNextCell(gt,list(range(nbPtsInCell))) + + inPts = DataArrayDouble( [(9.1, 18.2, 27.3), (9.1, 21.8, 27.3), (10.9, 21.8, 27.3), (10.9, 18.2, 27.3), (9.1, 18.2, 32.7), (9.1, 21.8, 32.7), (10.9, 21.8, 32.7), (10.9, 18.2, 32.7), (8.9, 20.0, 26.7), (10.0, 21.4, 27.9), (11.1, 20.0, 26.7), (10.0, 18.6, 27.9), (8.9, 20.0, 33.3), (10.0, 21.4, 32.1), (11.1, 20.0, 33.3), (10.0, 18.6, 32.1), (8.9, 17.8, 30.0), (9.3, 21.4, 30.0), (11.1, 22.2, 30.0), (10.7, 18.6, 30.0), (10.0, 20.0, 26.7), (9.3, 20.0, 30.0), (10.0, 22.2, 30.0), (10.7, 20.0, 30.0), (10.0, 17.8, 30.0), (10.0, 20.0, 32.1), (10.0, 20.0, 30.0)] ) + + outPts = DataArrayDouble( [(8.9, 17.8, 26.7), (8.9, 22.2, 26.7), (11.1, 22.2, 26.7), (11.1, 17.8, 26.7), (8.9, 17.8, 33.3), (8.9, 22.2, 33.3), (11.1, 22.2, 33.3), (11.1, 17.8, 33.3), (8.7, 20.0, 26.1), (10.0, 21.8, 27.3), (11.3, 20.0, 26.1), (10.0, 18.2, 27.3), (8.7, 20.0, 33.9), (10.0, 21.8, 32.7), (11.3, 20.0, 33.9), (10.0, 18.2, 32.7), (8.7, 17.4, 30.0), (9.1, 21.8, 30.0), (11.3, 22.6, 30.0), (10.9, 18.2, 30.0), (10.0, 20.0, 26.1), (9.1, 20.0, 30.0), (10.0, 22.6, 30.0), (10.9, 20.0, 30.0), (10.0, 17.4, 30.0), (10.0, 20.0, 32.7), (14.0, 20.0, 30.0)] ) + + srcField = MEDCouplingFieldDouble(ON_NODES_FE) + srcField.setMesh(m) + arr = DataArrayDouble(nbPtsInCell) + arr.iota() + srcField.setArray(arr) + + ref_val0 = DataArrayDouble( [6.5782430384766535, 5.505760346014328, 7.80996256527073, 7.643290943690746, 9.758765408487054, 9.06408508454036, 11.003779543627997, 11.205026363340515, 10.56416007071563, 18.44359561721225, 12.499588353132655, 19.85351355074463, 14.186041573114885, 19.339743214851023, 16.084629460041207, 20.892007336402663, 17.269258227200577, 19.549962126638338, 19.039562190136063, 21.648928068870756, 20.667409503475078, 22.062499999998867, 22.562500000009678, 23.812499999505995, 24.395833333387696, 25.62468592706991, 26.] ) + + eps = 1e-8 + + for i in range(nbPtsInCell): + self.assertTrue( abs( srcField.getValueOn( inPts[i] )[0]-ref_val0[i] ) < eps ) + + self.assertTrue( srcField.getValueOnMulti( inPts ).isEqual(ref_val0,eps) ) + + srcFt=MEDCouplingFieldTemplate(srcField) + trgFt=MEDCouplingFieldTemplate(ON_NODES_FE) + trgMesh = MEDCouplingUMesh.Build0DMeshFromCoords( inPts ) + trgFt.setMesh(trgMesh) + rem = MEDCouplingRemapper() + rem.setIntersectionType( PointLocator ) + rem.prepareEx(srcFt,trgFt) + # scan content of matrix computed by remapper + mat = rem.getCrudeMatrix() + self.assertEqual( len(mat) , nbPtsInCell ) + for irow,row in enumerate(mat): + self.assertTrue( abs( sum([y for x,y in row.items()]) - 1.0) < eps ) + self.assertEqual( irow , [x for x,y in row.items() if y == max([y for x,y in row.items()])][0] ) # check that max coeff is for irow elt (due to construction of inPts ) + + # ask for MEDCouplingFieldDiscretizationOnNodesFE instance to compute coordination into ref element + sd = srcField.getDiscretization() + coosInRefElemFoundByNewton = sd.getCooInRefElement(srcField.getMesh(),inPts) + + for zePt,cooInRefElemFoundByNewton in zip(inPts,coosInRefElemFoundByNewton): + # now check by performing refCoo -> realCoo + ftest = MEDCouplingFieldDouble(ON_GAUSS_PT) + ftest.setMesh(srcField.getMesh()) + ftest.setGaussLocalizationOnType(gt,sum([list(elt) for elt in MEDCouplingGaussLocalization.GetDefaultReferenceCoordinatesOf(gt).getValuesAsTuple()],[]),list(cooInRefElemFoundByNewton),[1]) + self.assertTrue ( float( (ftest.getLocalizationOfDiscr()-zePt).magnitude() ) < 1e-4 ) + def checkMatrix(self,mat1,mat2,nbCols,eps): self.assertEqual(len(mat1),len(mat2)) for i in range(len(mat1)): diff --git a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i index 775ba4b48..e6029473e 100644 --- a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i +++ b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i @@ -28,6 +28,7 @@ #include "MEDCouplingMappedExtrudedMesh.hxx" #include "MEDCoupling1GTUMesh.hxx" #include "MEDCouplingFieldDiscretization.hxx" +#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx" #include "MEDCouplingMultiFields.hxx" #include "MEDCouplingPartDefinition.hxx" #include "MEDCouplingCartesianAMRMesh.hxx" @@ -121,6 +122,8 @@ static PyObject *convertFieldDiscretization(MEDCoupling::MEDCouplingFieldDiscret ret=SWIG_NewPointerObj(reinterpret_cast(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationGaussNE,owner); if(dynamic_cast(fd)) ret=SWIG_NewPointerObj(reinterpret_cast(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationKriging,owner); + if(dynamic_cast(fd)) + ret=SWIG_NewPointerObj(reinterpret_cast(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationOnNodesFE,owner); if(!ret) throw INTERP_KERNEL::Exception("Not recognized type of field discretization on downcast !"); return ret;