}
//---------------------------------------------------------------
-#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 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[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();
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
};
* 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::getFunctionValues() const
+MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::getShapeFunctionValues() const
{
MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
int nbGaussPt(getNumberOfGaussPt()),nbPtsRefCell(getNumberOfPtsInRefCell()),dim(getDimension());
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
//
MEDCOUPLING_EXPORT MCAuto<DataArrayDouble> localizePtsInRefCooForEachCell(const DataArrayDouble *ptsInRefCoo, const MEDCouplingUMesh *mesh) const;
MEDCOUPLING_EXPORT MCAuto<MEDCouplingUMesh> buildRefCell() const;
- MEDCOUPLING_EXPORT MCAuto<DataArrayDouble> getFunctionValues() 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;
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
%newobject MEDCoupling::DenseMatrix::__mul__;
%newobject MEDCoupling::MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell;
%newobject MEDCoupling::MEDCouplingGaussLocalization::buildRefCell;
-%newobject MEDCoupling::MEDCouplingGaussLocalization::getFunctionValues;
+%newobject MEDCoupling::MEDCouplingGaussLocalization::getShapeFunctionValues;
+%newobject MEDCoupling::MEDCouplingGaussLocalization::getDerivativeOfShapeFunctionValues;
%newobject MEDCoupling::MEDCouplingSkyLineArray::BuildFromPolyhedronConn;
%newobject MEDCoupling::MEDCouplingSkyLineArray::getSuperIndexArray;
%newobject MEDCoupling::MEDCouplingSkyLineArray::getIndexArray;
return ret.retn();
}
- DataArrayDouble *getFunctionValues() const
+ DataArrayDouble *getShapeFunctionValues() const
{
- MCAuto<DataArrayDouble> ret(self->getFunctionValues());
+ MCAuto<DataArrayDouble> ret(self->getShapeFunctionValues());
+ return ret.retn();
+ }
+
+ DataArrayDouble *getDerivativeOfShapeFunctionValues() const
+ {
+ MCAuto<DataArrayDouble> ret(self->getDerivativeOfShapeFunctionValues());
return ret.retn();
}
}
