- 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
#include <algorithm>
#include <iostream>
+#include <memory>
#include <limits>
#include <cmath>
template <int dim, class ConnType = int>
class BBTree
{
-
private:
BBTree* _left;
BBTree* _right;
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.
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)
{
_terminal=true;
}
- double* nodes=new double [nbelems];
- _elems.resize(nbelems);
- for (ConnType i=0; i<nbelems; i++)
- {
- ConnType elem;
- if (elems!=0)
- elem= elems[i];
- else
- elem=i;
+ double median = std::numeric_limits<double>::max();
+ {
+ std::unique_ptr<double[]> nodes( new double [nbelems] );
+ _elems.resize(nbelems);
+ for (ConnType i=0; i<nbelems; i++)
+ {
+ ConnType elem;
+ if (elems)
+ elem= elems[i];
+ else
+ elem=i;
- _elems[i]=elem;
- nodes[i]=bbs[elem*dim*2+(level%dim)*2];
- }
- if (_terminal) { delete[] nodes; return;}
+ _elems[i]=elem;
+ nodes[i]=bbs[elem*dim*2+(level%dim)*2];
+ }
+ if (_terminal) { return; }
- std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
- double median = *(nodes+nbelems/2);
- delete[] nodes;
+ std::nth_element<double*>(nodes.get(), nodes.get()+nbelems/2, nodes.get()+nbelems);
+ median = *(nodes.get()+nbelems/2);
+ }
// std:: cout << *median <<std::endl;
std::vector<ConnType> new_elems_left;
_right=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_right.size(),_epsilon);
}
-
-
~BBTree()
{
}
-
/*! 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<ConnType>& elems) const
{
// terminal node : return list of elements intersecting bb
\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<ConnType>& elems) const
{
// terminal node : return list of elements intersecting bb
_right->getElementsAroundPoint(xx,elems);
}
-
-
ConnType size()
{
if (_terminal) return _nbelems;
--- /dev/null
+// 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 <memory>
+
+/*!
+ * Wrapper over BBTree to deal with ownership of bbox double array.
+ */
+template <int dim, class ConnType>
+class BBTreeStandAlone
+{
+private:
+ std::unique_ptr<double[]> _bbox;
+ BBTree<dim,ConnType> _effective;
+public:
+ BBTreeStandAlone(std::unique_ptr<double[]>&& 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<ConnType>& elems) const { _effective.getIntersectingElems(bb,elems); }
+ void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const { _effective.getElementsAroundPoint(xx,elems); }
+};
#include <string>
#include <exception>
+#include <sstream>
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]);
* @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]);
assert(isValid());
}
- /**
- * Destructor
- *
- */
- BoundingBox::~BoundingBox()
- {
- delete [] _coords;
- }
-
/**
* Determines if the intersection with a given box is empty
*
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;
inline void dumpCoords() const;
void toCompactData(double data[6]) const;
+
+ BoundingBox& operator=(const BoundingBox& box) = delete;
private:
/// 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];
};
ExprEval/InterpKernelValue.cxx
ExprEval/InterpKernelAsmX86.cxx
GaussPoints/InterpKernelGaussCoords.cxx
+ LinearAlgebra/InterpKernelDenseMatrix.cxx
+ LinearAlgebra/InterpKernelLUDecomp.cxx
+ LinearAlgebra/InterpKernelQRDecomp.cxx
)
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)
}
//---------------------------------------------------------------
-#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 ) \
{ \
_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() );
}
/*!
return ret;
}
+std::vector<double> GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellType ct)
+{
+ switch(ct)
+ {
+ case INTERP_KERNEL::NORM_SEG2:
+ return std::vector<double>(SEG2A_REF,SEG2A_REF+sizeof(SEG2A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_SEG3:
+ return std::vector<double>(SEG3_REF,SEG3_REF+sizeof(SEG3_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_TRI3:
+ return std::vector<double>(TRIA3A_REF,TRIA3A_REF+sizeof(TRIA3A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_TRI6:
+ return std::vector<double>(TRIA6A_REF,TRIA6A_REF+sizeof(TRIA6A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_TRI7:
+ return std::vector<double>(TRIA7A_REF,TRIA7A_REF+sizeof(TRIA7A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_QUAD4:
+ return std::vector<double>(QUAD4A_REF,QUAD4A_REF+sizeof(QUAD4A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_QUAD8:
+ return std::vector<double>(QUAD8A_REF,QUAD8A_REF+sizeof(QUAD8A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_QUAD9:
+ return std::vector<double>(QUAD9A_REF,QUAD9A_REF+sizeof(QUAD9A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_TETRA4:
+ return std::vector<double>(TETRA4A_REF,TETRA4A_REF+sizeof(TETRA4A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_TETRA10:
+ return std::vector<double>(TETRA10A_REF,TETRA10A_REF+sizeof(TETRA10A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_PYRA5:
+ return std::vector<double>(PYRA5A_REF,PYRA5A_REF+sizeof(PYRA5A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_PYRA13:
+ return std::vector<double>(PYRA13A_REF,PYRA13A_REF+sizeof(PYRA13A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_PENTA6:
+ return std::vector<double>(PENTA6A_REF,PENTA6A_REF+sizeof(PENTA6A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_PENTA15:
+ return std::vector<double>(PENTA15A_REF,PENTA15A_REF+sizeof(PENTA15A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_PENTA18:
+ return std::vector<double>(PENTA18A_REF,PENTA18A_REF+sizeof(PENTA18A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_HEXA8:
+ return std::vector<double>(HEXA8A_REF,HEXA8A_REF+sizeof(HEXA8A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_HEXA20:
+ return std::vector<double>(HEXA20A_REF,HEXA20A_REF+sizeof(HEXA20A_REF)/sizeof(double));
+ case INTERP_KERNEL::NORM_HEXA27:
+ return std::vector<double>(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);
/*!
/**
* 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()
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;
}
/*!
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;
}
////////////////////////////////////////////////////////////////////////////////////////////////
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<double> NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector<double>& inputArray);
-
+ INTERPKERNEL_EXPORT static std::vector<double> 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];
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
};
--- /dev/null
+// 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 T>
+ 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<double>;
+}
--- /dev/null
+// 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 <class T>
+ DenseMatrixT<T>::DenseMatrixT() : nn(0), mm(0), v(nullptr) {}
+
+ template <class T>
+ DenseMatrixT<T>::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<n;i++) v[i] = v[i-1] + m;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::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<m; j++) v[i][j] = a;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::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<m; j++) v[i][j] = *a++;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::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<mm; j++) v[i][j] = rhs[i][j];
+ }
+
+ /*!
+ * postcondition: normal assignment via copying has been performed
+ * if matrix and rhs were different sizes, matrix
+ * has been resized to match the size of rhs
+ */
+ template <class T>
+ DenseMatrixT<T> & DenseMatrixT<T>::operator=(const DenseMatrixT<T> &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<mm; j++) v[i][j] = rhs[i][j];
+ }
+ return *this;
+ }
+
+ template <class T>
+ void DenseMatrixT<T>::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 <class T>
+ void DenseMatrixT<T>::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<mm; j++) v[i][j] = a;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::~DenseMatrixT()
+ {
+ if ( v )
+ {
+ delete[] (v[0]);
+ delete[] (v);
+ }
+ }
+}
--- /dev/null
+// 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 <vector>
+#include <limits>
+#include <cmath>
+
+template<class T>
+inline T sqr(const T a) { return a*a; }
+
+namespace INTERP_KERNEL
+{
+ template <class T>
+ void LineSearches(const std::vector<double>& xold, const double fold, const std::vector<double>& g, std::vector<double> &p,
+ std::vector<double> &x, double &f, const double stpmax, bool &check, T &func)
+ {
+ const double ALF=1.0e-4, TOLX=std::numeric_limits<double>::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<n;i++) sum += p[i]*p[i];
+ sum=std::sqrt(sum);
+ if (sum > stpmax)
+ for (i=0;i<n;i++)
+ p[i] *= stpmax/sum;
+ for (i=0;i<n;i++)
+ slope += g[i]*p[i];
+ if (slope >= 0.0) THROW_IK_EXCEPTION("Roundoff problem in LineSearches.");
+ test=0.0;
+ for (i=0;i<n;i++)
+ {
+ temp=std::abs(p[i])/std::fmax(std::abs(xold[i]),1.0);
+ if (temp > test) test=temp;
+ }
+ alamin=TOLX/test;
+ alam=1.0;
+ for (;;)
+ {
+ for (i=0;i<n;i++) x[i]=xold[i]+alam*p[i];
+ f=func(x);
+ if (alam < alamin)
+ {
+ for (i=0;i<n;i++) x[i]=xold[i];
+ check=true;
+ return;
+ } else if (f <= fold+ALF*alam*slope) return;
+ else
+ {
+ if (alam == 1.0)
+ tmplam = -slope/(2.0*(f-fold-slope));
+ else
+ {
+ rhs1=f-fold-alam*slope;
+ rhs2=f2-fold-alam2*slope;
+ a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
+ b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
+ if (a == 0.0) tmplam = -slope/(2.0*b);
+ else {
+ disc=b*b-3.0*a*slope;
+ if (disc < 0.0) tmplam=0.5*alam;
+ else if (b <= 0.0) tmplam=(-b+std::sqrt(disc))/(3.0*a);
+ else tmplam=-slope/(b+std::sqrt(disc));
+ }
+ if (tmplam>0.5*alam)
+ tmplam=0.5*alam;
+ }
+ }
+ alam2=alam;
+ f2 = f;
+ alam=std::fmax(tmplam,0.1*alam);
+ }
+ }
+ template <class T>
+ 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<double>& x, const std::vector<double>& fvec)
+ {
+ mcIdType n=x.size();
+ INTERP_KERNEL::DenseMatrix df(n,n);
+ std::vector<double> xh=x;
+ for (mcIdType j=0;j<n;j++)
+ {
+ double temp=xh[j];
+ double h=EPS*std::abs(temp);
+ if (h == 0.0) h=EPS;
+ xh[j]=temp+h;
+ h=xh[j]-temp;
+ std::vector<double> f=func(xh);
+ xh[j]=temp;
+ for (mcIdType i=0;i<n;i++)
+ df[i][j]=(f[i]-fvec[i])/h;
+ }
+ return df;
+ }
+ };
+
+ template <class T>
+ class FMin
+ {
+ private:
+ std::vector<double> fvec;
+ T &func;
+ mcIdType n;
+ public:
+ FMin(T &funcc) : func(funcc){}
+ double operator() (const std::vector<double>& x)
+ {
+ n=x.size();
+ double sum=0;
+ fvec=func(x);
+ for (mcIdType i=0;i<n;i++) sum += sqr(fvec[i]);
+ return 0.5*sum;
+ }
+ std::vector<double>& getVector() { return fvec; }
+ };
+
+ /*!
+ * check is false on normal return.
+ * check is true if the routine has converged to a local minimum.
+ */
+ template <class T>
+ void SolveWithNewtonWithJacobian(std::vector<double> &x, bool &check, T &vecfunc,
+ std::function<void(const std::vector<double>&, const std::vector<double>&, 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<double>::epsilon();
+ mcIdType i,j,its,n=x.size();
+ double den,f,fold,stpmax,sum,temp,test;
+ std::vector<double> g(n),p(n),xold(n);
+ INTERP_KERNEL::DenseMatrix fjac(n,n);
+ FMin<T> fmin(vecfunc);
+ std::vector<double> &fvec=fmin.getVector();
+ f=fmin(x);
+ test=0.0;
+ for (i=0;i<n;i++)
+ if (std::abs(fvec[i]) > test) test=std::abs(fvec[i]);
+ if (test < 0.01*TOLF)
+ {
+ check=false;
+ return;
+ }
+ sum=0.0;
+ for (i=0;i<n;i++) sum += sqr(x[i]);
+ stpmax=STPMX*std::fmax(std::sqrt(sum),double(n));
+ for (its=0;its<MAXITS;its++)
+ {
+ jacobian(x,fvec,fjac);//fjac=fdjac(x,fvec);
+ for (i=0;i<n;i++)
+ {
+ sum=0.0;
+ for (j=0;j<n;j++) sum += fjac[j][i]*fvec[j];
+ g[i]=sum;
+ }
+ for (i=0;i<n;i++) xold[i]=x[i];
+ fold=f;
+ for (i=0;i<n;i++) p[i] = -fvec[i];
+ INTERP_KERNEL::LUDecomp alu(fjac);
+ alu.solve(p,p);
+ LineSearches(xold,fold,g,p,x,f,stpmax,check,fmin);
+ test=0.0;
+ for (i=0;i<n;i++)
+ if (std::abs(fvec[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<n;i++) {
+ temp=std::abs(g[i])*std::fmax(std::abs(x[i]),1.0)/den;
+ if (temp > test) test=temp;
+ }
+ check=(test < TOLMIN);
+ return;
+ }
+ test=0.0;
+ for (i=0;i<n;i++)
+ {
+ temp=(std::abs(x[i]-xold[i]))/std::fmax(std::abs(x[i]),1.0);
+ if (temp > 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 <class T>
+ void SolveWithNewton(std::vector<double> &x, bool &check, T &vecfunc)
+ {
+ JacobianCalculator<T> fdjac(vecfunc);
+ auto myJacobian = [&fdjac,vecfunc](const std::vector<double>& x, const std::vector<double>& fvec, INTERP_KERNEL::DenseMatrix& fjac)
+ {
+ fjac = fdjac(x,fvec);
+ };
+ SolveWithNewtonWithJacobian(x,check,vecfunc,myJacobian);
+ }
+}
// create BBTree structure
// - get bounding boxes
std::vector<double> 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<SPACEDIM,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems);
+ BBTree<SPACEDIM,ConnType> 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());
}
}
}
- delete [] srcElemIdx;
for(ConnType i = 0 ; i < numSrcElems ; ++i)
delete srcElems[i];
return srcMesh.getNumberOfNodes();
protected:
template<class MyMeshType, class MyMatrixType>
- void performAdjustmentOfBB(Intersector3D<MyMeshType,MyMatrixType>* intersector, std::vector<double>& bbox) const
+ void performAdjustmentOfBB(Intersector3D<MyMeshType,MyMatrixType>* intersector, double *bbox, std::size_t sz) const
{
- intersector->adjustBoundingBoxes(bbox,InterpolationOptions::getBoundingBoxAdjustment(),InterpolationOptions::getBoundingBoxAdjustmentAbs());
+ intersector->adjustBoundingBoxes(bbox,sz,InterpolationOptions::getBoundingBoxAdjustment(),InterpolationOptions::getBoundingBoxAdjustmentAbs());
}
private:
#include "PointLocator3DIntersectorP1P0.txx"
#include "PolyhedronIntersectorP1P1.txx"
#include "PointLocator3DIntersectorP1P1.txx"
-#include "Log.hxx"
+#include "InterpolationHelper.txx"
#include "BBTree.txx"
{
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<MeshElement<ConnType>*> srcElems(numSrcElems);
- std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
+ DuplicateFacesType intersectFaces;
- std::map<MeshElement<ConnType>*, ConnType> indices;
- DuplicateFacesType intersectFaces;
-
- for(ConnType i = 0 ; i < numSrcElems ; ++i)
- srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
-
- for(ConnType i = 0 ; i < numTargetElems ; ++i)
- targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
-
- Intersector3D<MyMeshType,MyMatrixType>* intersector=0;
+ std::unique_ptr< Intersector3D<MyMeshType,MyMatrixType> > intersector;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel());
if(methC=="P0P0")
switch(InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,
+ intersector.reset( new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(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.");
matrix.resize(intersector->getNumberOfRowsOfResMatrix());
// create BBTree structure
- // - get bounding boxes
- std::vector<double> 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<ConnType> 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<ConnType> 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)
{
// 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;
}
#else // use BBTree class
-#include "BBTree.txx"
+#include "InterpolationHelper.txx"
#endif
+#include <memory>
+
namespace INTERP_KERNEL
{
/**
template<class MyMeshType, class MatrixType>
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<MeshElement<ConnType>*> srcElems(numSrcElems);
- std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
-
- std::map<MeshElement<ConnType>*, ConnType> indices;
+ LOG(2, "Target mesh has " << numTargetElems << " elements ");
- for(ConnType i = 0 ; i < numSrcElems ; ++i)
- srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
-
- for(ConnType i = 0 ; i < numTargetElems ; ++i)
- targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
-
- Intersector3D<MyMeshType,MatrixType>* intersector=0;
+ std::unique_ptr<Intersector3D<MyMeshType,MatrixType>> intersector;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
if(methC=="P0P0")
{
switch(InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new PolyhedronIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
+ intersector.reset( new PolyhedronIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
break;
case PointLocator:
- intersector=new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ intersector.reset( new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
break;
default:
throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangle or PointLocator.");
switch(InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new PolyhedronIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
+ intersector.reset( new PolyhedronIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
break;
case PointLocator:
- intersector=new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ intersector.reset( new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
break;
default:
throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P1 interp specified : must be Triangle or PointLocator.");
switch(InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new PolyhedronIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
+ intersector.reset( new PolyhedronIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
break;
case PointLocator:
- intersector=new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ intersector.reset( new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
break;
case Barycentric:
- intersector=new PolyhedronIntersectorP1P0Bary<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
+ intersector.reset( new PolyhedronIntersectorP1P0Bary<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
break;
default:
throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P0 interp specified : must be Triangle, PointLocator or Barycentric.");
switch(InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new PolyhedronIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
+ intersector.reset( new PolyhedronIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
break;
case PointLocator:
- intersector=new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ intersector.reset( new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
break;
case Barycentric:
- intersector=new Barycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ intersector.reset( new Barycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
break;
case MappedBarycentric:
- intersector=new MappedBarycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ intersector.reset( new MappedBarycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
break;
default:
throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle, PointLocator, Barycentric or MappedBarycentric.");
}
#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<ConnType> 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<ConnType> 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();
}
}
Interpolation3D1D::Interpolation3D1D(const InterpolationOptions& io):Interpolation<Interpolation3D1D>(io) { }
- /**
- * Inspired from PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes
- */
- void Interpolation3D1D::adjustBoundingBoxes(std::vector<double>& 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<size; i++)
{
double max=- std::numeric_limits<double>::max();
}
}
}
+
+ /**
+ * Inspired from PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes
+ */
+ void Interpolation3D1D::adjustBoundingBoxes(std::vector<double>& bbox)
+ {
+ adjustBoundingBoxes(bbox.data(),bbox.size());
+ }
}
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-// Author : A Bruneton (CEA/DEN)\r
-\r
-#pragma once\r
-\r
-#include "INTERPKERNELDefines.hxx"\r
-#include "Interpolation.hxx"\r
-#include "NormalizedUnstructuredMesh.hxx"\r
-#include "InterpolationOptions.hxx"\r
-\r
-#include <vector>\r
-\r
-namespace INTERP_KERNEL\r
+// Author : A Bruneton (CEA/DEN)
+
+#pragma once
+
+#include "INTERPKERNELDefines.hxx"
+#include "Interpolation.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "InterpolationOptions.hxx"
+
+#include <vector>
+
+namespace INTERP_KERNEL
{
/**
* \class Interpolation3D1D
* Can be seen as a specialization of Interpolation3D, and allows notably the adjustment of bounding boxes.
*/
- class INTERPKERNEL_EXPORT Interpolation3D1D : public Interpolation<Interpolation3D1D>\r
- {\r
- public:\r
- Interpolation3D1D();\r
- Interpolation3D1D(const InterpolationOptions& io);\r
- template<class MyMeshType, class MatrixType>\r
- typename MyMeshType::MyConnType interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method);\r
- private:\r
- void adjustBoundingBoxes(std::vector<double>& bbox);\r
- };\r
-}\r
-\r
+ class INTERPKERNEL_EXPORT Interpolation3D1D : public Interpolation<Interpolation3D1D>
+ {
+ public:
+ Interpolation3D1D();
+ Interpolation3D1D(const InterpolationOptions& io);
+ template<class MyMeshType, class MatrixType>
+ 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<double>& bbox);
+ };
+}
+
#include "PointLocator3DIntersectorP1P1.txx"
#include "Log.hxx"
-#include "BBTree.txx"
-
-#include <memory>
+#include "InterpolationHelper.txx"
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<ConnType> > > srcElems(numSrcElems);
- std::vector< std::unique_ptr< MeshElement<ConnType> > > targetElems(numTargetElems);
-
- std::map<MeshElement<ConnType>*, ConnType> indices;
-
- for(ConnType i = 0 ; i < numSrcElems ; ++i)
- srcElems[i].reset( new MeshElement<ConnType>(i, srcMesh) );
-
- for(ConnType i = 0 ; i < numTargetElems ; ++i)
- targetElems[i].reset( new MeshElement<ConnType>(i, targetMesh) );
+ LOG(2, "Target mesh has " << numTargetElems << " elements ");
std::unique_ptr< Intersector3D<MyMeshType,MatrixType> > intersector;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
// create empty maps for all source elements
result.resize(intersector->getNumberOfRowsOfResMatrix());
- // create BBTree structure
- // - get bounding boxes
- std::vector<double> bboxes(6*numSrcElems);
- std::unique_ptr<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);
-
- 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<ConnType> 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<ConnType> intersectElems;
tree.getIntersectingElems(targetBox, intersectElems);
if ( !intersectElems.empty() )
- intersector->intersectCells(targetIdx,intersectElems,result);
+ intersector->intersectCells(i,intersectElems,result);
}
return intersector->getNumberOfColsOfResMatrix();
}
--- /dev/null
+// 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 <memory>
+#include <functional>
+
+namespace INTERP_KERNEL
+{
+ template<class MyMeshType>
+ BBTreeStandAlone<3,typename MyMeshType::MyConnType> BuildBBTree(const MyMeshType& srcMesh)
+ {
+ return BuildBBTreeWithAdjustment(srcMesh,[](double *,typename MyMeshType::MyConnType){});
+ }
+
+ template<class MyMeshType>
+ BBTreeStandAlone<3,typename MyMeshType::MyConnType> BuildBBTreeWithAdjustment(const MyMeshType& srcMesh, std::function<void(double *,typename MyMeshType::MyConnType)> 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<double[]> bboxes( new double[nbElts] );
+ for(ConnType i = 0; i < numSrcElems ; ++i)
+ {
+ MeshElement<ConnType> 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);
+ }
+}
--- /dev/null
+// 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<double>;
--- /dev/null
+// 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 <sstream>
+
+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<double> vv(n);
+ d=1.0;
+ for (i=0;i<n;i++)
+ {
+ big=0.0;
+ for (j=0;j<n;j++)
+ if ((temp=std::abs(lu[i][j])) > big) big=temp;
+ if (big == 0.0) THROW_IK_EXCEPTION("Singular matrix in LUDecomp");
+ vv[i]=1.0/big;
+ }
+ for (k=0;k<n;k++)
+ {
+ big=0.0;
+ imax=k;
+ for (i=k;i<n;i++)
+ {
+ temp=vv[i]*std::abs(lu[i][k]);
+ if (temp > big)
+ {
+ big=temp;
+ imax=i;
+ }
+ }
+ if (k != imax) {
+ for (j=0;j<n;j++) {
+ temp=lu[imax][j];
+ lu[imax][j]=lu[k][j];
+ lu[k][j]=temp;
+ }
+ d = -d;
+ vv[imax]=vv[k];
+ }
+ indx[k]=imax;
+ if (lu[k][k] == 0.0) lu[k][k]=TINY;
+ for (i=k+1;i<n;i++) {
+ temp=lu[i][k] /= lu[k][k];
+ for (j=k+1;j<n;j++)
+ lu[i][j] -= temp*lu[k][j];
+ }
+ }
+}
+void LUDecomp::solve(const std::vector<double>& b, std::vector<double> &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<n;i++) x[i] = b[i];
+ for (i=0;i<n;i++) {
+ ip=indx[i];
+ sum=x[ip];
+ x[ip]=x[i];
+ if (ii != 0)
+ for (j=ii-1;j<i;j++) sum -= lu[i][j]*x[j];
+ else if (sum != 0.0)
+ ii=i+1;
+ x[i]=sum;
+ }
+ for (i=n-1;i>=0;i--) {
+ sum=x[i];
+ for (j=i+1;j<n;j++) sum -= lu[i][j]*x[j];
+ x[i]=sum/lu[i][i];
+ }
+}
+
+void LUDecomp::solve(const INTERP_KERNEL::DenseMatrix& b, INTERP_KERNEL::DenseMatrix &x)
+{
+ mcIdType i,j,m(b.ncols());
+ if (b.nrows() != n || x.nrows() != n || b.ncols() != x.ncols())
+ THROW_IK_EXCEPTION("LUDecomp::solve bad sizes");
+ std::vector<double> xx(n);
+ for (j=0;j<m;j++) {
+ for (i=0;i<n;i++) xx[i] = b[i][j];
+ solve(xx,xx);
+ for (i=0;i<n;i++) x[i][j] = xx[i];
+ }
+}
+
+void LUDecomp::inverse(INTERP_KERNEL::DenseMatrix &ainv)
+{
+ mcIdType i,j;
+ ainv.resize(n,n);
+ for (i=0;i<n;i++) {
+ for (j=0;j<n;j++) ainv[i][j] = 0.;
+ ainv[i][i] = 1.;
+ }
+ solve(ainv,ainv);
+}
+
+double LUDecomp::det()
+{
+ double dd = d;
+ for (mcIdType i=0;i<n;i++) dd *= lu[i][i];
+ return dd;
+}
+
+void LUDecomp::mprove(const std::vector<double>& b, std::vector<double> &x)
+{
+ mcIdType i,j;
+ std::vector<double> r(n);
+ for (i=0;i<n;i++) {
+ long double sdp = -b[i];
+ for (j=0;j<n;j++)
+ sdp += (long double)aref[i][j] * (long double)x[j];
+ r[i]=double(sdp);
+ }
+ solve(r,r);
+ for (i=0;i<n;i++) x[i] -= r[i];
+}
--- /dev/null
+// 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"
+#include <vector>
+
+namespace INTERP_KERNEL
+{
+ class LUDecomp
+ {
+ private:
+ mcIdType n;
+ INTERP_KERNEL::DenseMatrix lu;
+ std::vector<mcIdType> indx;
+ double d;
+ const INTERP_KERNEL::DenseMatrix& aref;
+ public:
+ LUDecomp (const INTERP_KERNEL::DenseMatrix& a);
+ void solve(const std::vector<double>& b, std::vector<double> &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<double>& b, std::vector<double> &x);
+ };
+}
--- /dev/null
+// 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 <cmath>
+
+using namespace INTERP_KERNEL;
+
+template<class T>
+inline T sqr(const T a) { return a*a; }
+
+template<class T>
+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<double> c(n), d(n);
+ double scale,sigma,sum,tau;
+ for (k=0;k<n-1;k++) {
+ scale=0.0;
+ for (i=k;i<n;i++) scale=std::fmax(scale,std::abs(r[i][k]));
+ if (scale == 0.0) {
+ sing=true;
+ c[k]=d[k]=0.0;
+ } else {
+ for (i=k;i<n;i++) r[i][k] /= scale;
+ for (sum=0.0,i=k;i<n;i++) sum += sqr(r[i][k]);
+ sigma=sign(std::sqrt(sum),r[k][k]);
+ r[k][k] += sigma;
+ c[k]=sigma*r[k][k];
+ d[k] = -scale*sigma;
+ for (j=k+1;j<n;j++) {
+ for (sum=0.0,i=k;i<n;i++) sum += r[i][k]*r[i][j];
+ tau=sum/c[k];
+ for (i=k;i<n;i++) r[i][j] -= tau*r[i][k];
+ }
+ }
+ }
+ d[n-1]=r[n-1][n-1];
+ if (d[n-1] == 0.0) sing=true;
+ for (i=0;i<n;i++) {
+ for (j=0;j<n;j++) qt[i][j]=0.0;
+ qt[i][i]=1.0;
+ }
+ for (k=0;k<n-1;k++) {
+ if (c[k] != 0.0) {
+ for (j=0;j<n;j++) {
+ sum=0.0;
+ for (i=k;i<n;i++)
+ sum += r[i][k]*qt[i][j];
+ sum /= c[k];
+ for (i=k;i<n;i++)
+ qt[i][j] -= sum*r[i][k];
+ }
+ }
+ }
+ for (i=0;i<n;i++) {
+ r[i][i]=d[i];
+ for (j=0;j<i;j++) r[i][j]=0.0;
+ }
+}
+
+void QRDecomp::solve(const std::vector<double>& b, std::vector<double> &x)
+{
+ qtmult(b,x);
+ rsolve(x,x);
+}
+
+void QRDecomp::qtmult(const std::vector<double>& b, std::vector<double> &x)
+{
+ mcIdType i,j;
+ double sum;
+ for (i=0;i<n;i++)
+ {
+ sum = 0.;
+ for (j=0;j<n;j++) sum += qt[i][j]*b[j];
+ x[i] = sum;
+ }
+}
+
+void QRDecomp::rsolve(const std::vector<double>& b, std::vector<double> &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<n;j++) sum -= r[i][j]*x[j];
+ x[i]=sum/r[i][i];
+ }
+}
+void QRDecomp::update(const std::vector<double>& u, const std::vector<double>& v)
+{
+ mcIdType i,k;
+ std::vector<double> 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<n;i++) r[0][i] += w[0]*v[i];
+ for (i=0;i<k;i++)
+ rotate(i,r[i][i],-r[i+1][i]);
+ for (i=0;i<n;i++)
+ if (r[i][i] == 0.0) sing=true;
+}
+
+void QRDecomp::rotate(const mcIdType i, const double a, const double b)
+{
+ mcIdType j;
+ double c,fact,s,w,y;
+ if (a == 0.0)
+ {
+ c=0.0;
+ s=(b >= 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<n;j++)
+ {
+ y=r[i][j];
+ w=r[i+1][j];
+ r[i][j]=c*y-s*w;
+ r[i+1][j]=s*y+c*w;
+ }
+ for (j=0;j<n;j++)
+ {
+ y=qt[i][j];
+ w=qt[i+1][j];
+ qt[i][j]=c*y-s*w;
+ qt[i+1][j]=s*y+c*w;
+ }
+}
--- /dev/null
+// 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)
+
+#pragma once
+
+#include "InterpKernelDenseMatrix.hxx"
+#include <vector>
+
+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<double>& b, std::vector<double> &x);
+ void qtmult(const std::vector<double>& b, std::vector<double> &x);
+ void rsolve(const std::vector<double>& b, std::vector<double> &x);
+ void update(const std::vector<double>& u, const std::vector<double>& v);
+ void rotate(const mcIdType i, const double a, const double b);
+ };
+}
template<class MyMeshType>
MeshElement(const ConnType index, const MyMeshType& mesh);
- ~MeshElement();
-
- ConnType getIndex() const { return _index; }
+ MeshElement() = default;
+
+ template<class MyMeshType>
+ 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;
};
/**
#include <assert.h>
#include <type_traits>
#include <limits>
+#include <memory>
namespace INTERP_KERNEL
{
*/
template<class ConnType>
template<class MyMeshType>
- MeshElement<ConnType>::MeshElement(const ConnType index, const MyMeshType& mesh)
- : _index(index), _number( 0 ), _box(nullptr)
+ MeshElement<ConnType>::MeshElement(const ConnType index, const MyMeshType& mesh): _number( 0 )
+ {
+ this->assign(index,mesh);
+ }
+
+ template<class ConnType>
+ template<class MyMeshType>
+ void MeshElement<ConnType>::assign(const ConnType index, const MyMeshType& mesh)
{
auto numberCore = mesh.getNumberOfNodesOfElement(OTT<typename MyMeshType::MyConnType,MyMeshType::My_numPol>::indFC(index));
if(numberCore < std::numeric_limits<nbnodesincelltype>::max())
{
_number = static_cast< nbnodesincelltype >(numberCore);
- const double**vertices = new const double*[_number];
+ std::unique_ptr<const double*[]> vertices( new const double*[_number] );
for( nbnodesincelltype i = 0 ; i < _number ; ++i)
vertices[i] = getCoordsOfNode(i , OTT<typename MyMeshType::MyConnType,MyMeshType::My_numPol>::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<class ConnType>
- MeshElement<ConnType>::~MeshElement()
- {
- delete _box;
- }
/////////////////////////////////////////////////////////////////////
/// ElementBBoxOrder /////////////
virtual ~TargetIntersector() { }
void adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEps, double adjustmentEpsAbs);
+ void adjustBoundingBoxes(double *bbox, std::size_t sz, double adjustmentEps, double adjustmentEpsAbs);
};
}
namespace INTERP_KERNEL
{
+ template<class MyMeshType, class MyMatrix>
+ void TargetIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& 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
@param adjustmentEpsAbs absolute adjustment value (added on each side of the BBox in each dimension)
*/
template<class MyMeshType, class MyMatrix>
- void TargetIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEps, double adjustmentEpsAbs)
+ void TargetIntersector<MyMeshType,MyMatrix>::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<size; i++)
{
double max=- std::numeric_limits<double>::max();
${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
MEDCouplingStructuredMesh.cxx
MEDCouplingTimeDiscretization.cxx
MEDCouplingFieldDiscretization.cxx
+ MEDCouplingFieldDiscretizationOnNodesFE.cxx
MEDCouplingRefCountObject.cxx
MEDCouplingPointSet.cxx
MEDCouplingFieldTemplate.cxx
// Author : Anthony Geay (EDF R&D)
#include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingCMesh.hxx"
#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingFieldDouble.hxx"
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.};
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.");
}
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
}
}
-template<class FIELD_DISC>
-MCAuto<MEDCouplingFieldDiscretization> MEDCouplingFieldDiscretization::EasyAggregate(std::vector<const MEDCouplingFieldDiscretization *>& 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<const FIELD_DISC *>(it));
- if(!itc)
- throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : same field discretization expected for all input discretizations !");
- }
- return fds[0]->clone();
-}
-
MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization()
{
}
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<mcIdType>& code, const std::vector<const DataArrayIdType *>& idsPerType) const;
MEDCOUPLING_EXPORT mcIdType getNumberOfTuples(const MEDCouplingMesh *mesh) const;
MEDCOUPLING_EXPORT mcIdType getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const;
MEDCOUPLING_EXPORT void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& 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;
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
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<MEDCouplingFieldDiscretization> aggregate(std::vector<const MEDCouplingFieldDiscretization *>& 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;
};
std::size_t getHeapMemorySizeWithoutChildren() const;
std::vector<const BigMemoryObject *> 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:
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;
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<mcIdType>& tinyInfo) const;
MEDCOUPLING_EXPORT void getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const;
MEDCOUPLING_EXPORT void finishUnserialization(const std::vector<double>& tinyInfo);
MEDCOUPLING_EXPORT void checkForUnserialization(const std::vector<mcIdType>& 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;
void commonUnserialization(const std::vector<mcIdType>& tinyInfo);
public:
static const char REPR[];
- static const TypeOfField TYPE;
+ static constexpr TypeOfField TYPE = ON_GAUSS_PT;
private:
std::vector<MEDCouplingGaussLocalization> _loc;
};
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<mcIdType>& code, const std::vector<const DataArrayIdType *>& idsPerType) const;
MEDCOUPLING_EXPORT mcIdType getNumberOfTuples(const MEDCouplingMesh *mesh) const;
MEDCOUPLING_EXPORT mcIdType getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const;
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;
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];
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<MEDCouplingFieldDiscretization> aggregate(std::vector<const MEDCouplingFieldDiscretization *>& fds) const override;
public://specific part
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<class FIELD_DISC>
+ MCAuto<MEDCouplingFieldDiscretization> MEDCouplingFieldDiscretization::EasyAggregate(std::vector<const MEDCouplingFieldDiscretization *>& 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<const FIELD_DISC *>(it));
+ if(!itc)
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : same field discretization expected for all input discretizations !");
+ }
+ return fds[0]->clone();
+ }
}
#endif
--- /dev/null
+// 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 <sstream>
+
+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<MEDCouplingFieldDiscretization> MEDCouplingFieldDiscretizationOnNodesFE::aggregate(std::vector<const MEDCouplingFieldDiscretization *>& fds) const
+{
+ return EasyAggregate<MEDCouplingFieldDiscretizationOnNodesFE>(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<const MEDCouplingFieldDiscretizationOnNodesFE *>(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<double> operator()(const std::vector<double>& x)
+ {
+ MEDCouplingGaussLocalization gl(_gl->getType(),_gl->getRefCoords(),x,{1.0});
+ MCAuto<DataArrayDouble> shapeFunc = gl.getShapeFunctionValues();
+ const double *shapeFuncPtr( shapeFunc->begin() );
+ std::vector<double> 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<double>& 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<double>& x, const std::vector<double>&, INTERP_KERNEL::DenseMatrix& jacobian)
+ {
+ MEDCouplingGaussLocalization mygl(gl.getType(),gl.getRefCoords(),x,{1.0});
+ MCAuto<DataArrayDouble> 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<double> 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<double>::max());
+ return false;
+}
+
+void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *ptsCoo, double *res) const
+{
+ MCAuto<DataArrayDouble> 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<void(const MEDCouplingGaussLocalization&, const std::vector<mcIdType>&)> 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<mcIdType> 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<mcIdType> conn;
+ umesh->getNodeIdsOfCell(*cellId,conn);
+ MCAuto<DataArrayDouble> refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) );
+ std::vector<double> refCooCpp(refCoo->begin(),refCoo->end());
+ std::vector<double> gsCoo(ptsCoo + iPt*3,ptsCoo + (iPt+1)*3);
+ MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.});
+ std::vector<double> 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<double> 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<const MEDCouplingUMesh *>(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<DataArrayDouble> ret(DataArrayDouble::New());
+ ret->alloc(nbOfTargetPoints,nbCompo);
+ double *res( ret->getPointer() );
+
+ auto arrayFeeder = [arr,&res,nbCompo](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
+ {
+ MCAuto<DataArrayDouble> 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<DataArrayDouble> MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement(const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const
+{
+ const MEDCouplingUMesh *umesh = checkConfig3D(mesh);
+ MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+ ret->alloc(nbOfPoints,3);
+ double *retPtr(ret->getPointer() );
+
+ auto arrayFeeder = [&retPtr](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
+ {
+ std::vector<double> resVector( gl.getGaussCoords() );
+ {
+ std::copy(resVector.begin(),resVector.end(),retPtr);
+ retPtr += 3;
+ }
+ };
+
+ GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfPoints,arrayFeeder);
+ return ret;
+}
--- /dev/null
+// 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 <functional>
+
+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<MEDCouplingFieldDiscretization> aggregate(std::vector<const MEDCouplingFieldDiscretization *>& 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<DataArrayDouble> 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<void(const MEDCouplingGaussLocalization&, const std::vector<mcIdType>&)> customFunc);
+ private:
+ const MEDCouplingUMesh *checkConfig3D(const MEDCouplingMesh *mesh) const;
+ public:
+ static const char REPR[];
+ static constexpr TypeOfField TYPE = ON_NODES_FE;
+ };
+}
* \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.
* \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.<br>
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;
}
return MCAuto<MEDCouplingUMesh>(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<DataArrayDouble> MEDCouplingGaussLocalization::getShapeFunctionValues() const
+{
+ MCAuto<DataArrayDouble> 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<double> 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<DataArrayDouble> MEDCouplingGaussLocalization::getDerivativeOfShapeFunctionValues() const
+{
+ MCAuto<DataArrayDouble> 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<double> 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<nbOfNodePerCell or if not 0<=comp<dim
std::transform(tmp.begin(),tmp.end(),tmp.begin(),[](double c){return fabs(c);});
return *std::max_element(tmp.begin(),tmp.end())<eps;
}
+
+MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type)
+{
+ std::vector<double> 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<DataArrayDouble> 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;
+}
//
MEDCOUPLING_EXPORT MCAuto<DataArrayDouble> localizePtsInRefCooForEachCell(const DataArrayDouble *ptsInRefCoo, const MEDCouplingUMesh *mesh) const;
MEDCOUPLING_EXPORT MCAuto<MEDCouplingUMesh> buildRefCell() const;
+ MEDCOUPLING_EXPORT MCAuto<DataArrayDouble> getShapeFunctionValues() const;
+ MEDCOUPLING_EXPORT MCAuto<DataArrayDouble> getDerivativeOfShapeFunctionValues() const;
//
MEDCOUPLING_EXPORT const std::vector<double>& getRefCoords() const { return _ref_coord; }
MEDCOUPLING_EXPORT double getRefCoord(int ptIdInCell, int comp) const;
//
MEDCOUPLING_EXPORT static MEDCouplingGaussLocalization BuildNewInstanceFromTinyInfo(mcIdType dim, const std::vector<mcIdType>& tinyData);
MEDCOUPLING_EXPORT static bool AreAlmostEqual(const std::vector<double>& v1, const std::vector<double>& v2, double eps);
+ MEDCOUPLING_EXPORT static MCAuto<DataArrayDouble> GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type);
private:
int checkCoherencyOfRequest(mcIdType gaussPtIdInCell, int comp) const;
private:
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);
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);
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
#include "MEDCouplingCMesh.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
#include "MEDCouplingNormalizedCartesianMesh.txx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "Interpolation1D.txx"
#include "Interpolation2DCurve.hxx"
}
}
_matrix=m;
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(srcNbElem);
+ synchronizeSizeOfSideMatricesAfterMatrixComputation(srcNbElem);
}
int MEDCouplingRemapper::prepareInterpKernelOnly()
{
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 !";
}
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;
}
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;
}
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;
}
}
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;
}
}
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;
}
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<DataArrayDouble> 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<const MEDCouplingUMesh *>(_src_ft->getMesh()) );
+ const MEDCouplingPointSet *trgMesh( dynamic_cast<const MEDCouplingPointSet *>(_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<mcIdType>& conn)
+ {
+ auto& row = this->_matrix[rowId++];
+ MCAuto<DataArrayDouble> 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;
}
{
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());
}
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).
//
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);
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
#include "MEDCouplingFieldOverTime.hxx"
#include "MEDCouplingDefinitionTime.hxx"
#include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingCartesianAMRMesh.hxx"
#include "MEDCouplingAMRAttribute.hxx"
#include "MEDCouplingMatrix.hxx"
%feature("autodoc", "1");
%feature("docstring");
+%newobject MEDCoupling::MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement;
%newobject MEDCoupling::MEDCouplingField::buildMeasureField;
%newobject MEDCoupling::MEDCouplingField::getLocalizationOfDiscr;
%newobject MEDCoupling::MEDCouplingField::computeTupleIdsToSelectFromCellIds;
%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;
%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();"
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
MCAuto<MEDCouplingUMesh> ret(self->buildRefCell());
return ret.retn();
}
+
+ DataArrayDouble *getShapeFunctionValues() const
+ {
+ MCAuto<DataArrayDouble> ret(self->getShapeFunctionValues());
+ return ret.retn();
+ }
+
+ DataArrayDouble *getDerivativeOfShapeFunctionValues() const
+ {
+ MCAuto<DataArrayDouble> ret(self->getDerivativeOfShapeFunctionValues());
+ return ret.retn();
+ }
+
+ static DataArrayDouble *GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type)
+ {
+ MCAuto<DataArrayDouble> ret(MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(type));
+ return ret.retn();
+ }
}
};
}
}
};
+
+ 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<double> v3;
+ const double *inp=convertObjToPossibleCpp5_Safe2(locs,sw,v0,v1,v2,v3,"wrap of MEDCouplingFieldDouble::getValueOnMulti",
+ mesh->getSpaceDimension(),true,nbPts);
+ MCAuto<DataArrayDouble> ret(self->getCooInRefElement(mesh,inp,nbPts));
+ return ret.retn();
+ }
+ }
+ };
}
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)):
#include "MEDCouplingMappedExtrudedMesh.hxx"
#include "MEDCoupling1GTUMesh.hxx"
#include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingMultiFields.hxx"
#include "MEDCouplingPartDefinition.hxx"
#include "MEDCouplingCartesianAMRMesh.hxx"
ret=SWIG_NewPointerObj(reinterpret_cast<void*>(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationGaussNE,owner);
if(dynamic_cast<MEDCoupling::MEDCouplingFieldDiscretizationKriging *>(fd))
ret=SWIG_NewPointerObj(reinterpret_cast<void*>(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationKriging,owner);
+ if(dynamic_cast<MEDCoupling::MEDCouplingFieldDiscretizationOnNodesFE *>(fd))
+ ret=SWIG_NewPointerObj(reinterpret_cast<void*>(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationOnNodesFE,owner);
if(!ret)
throw INTERP_KERNEL::Exception("Not recognized type of field discretization on downcast !");
return ret;
