From 083874a866268323aa868f9e8b383bbae8851725 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Thu, 8 Sep 2022 08:24:40 +0200 Subject: [PATCH] Management of derivatives of shape functions of NORM_HEXA8,NORM_PENTA6,NORM_PYRA5,NORM_PENTA15,NORM_HEXA20,NORM_HEXA27 --- .../GaussPoints/InterpKernelGaussCoords.cxx | 152 ++++++++++++++++++ .../MEDCouplingBasicsTest7.py | 38 +++++ 2 files changed, 190 insertions(+) diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx index 61feaf166..657e4a898 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx @@ -1622,6 +1622,30 @@ void GaussInfo::pyra5aInit() funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]); funValue[4] = gc[2]; SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = 0.5*(gc[0]+1)*(1.0 - gc[2]); + devFunValue[1] = -0.5*gc[1]*(1.0 - gc[2]); + devFunValue[2] = -0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0); + + devFunValue[3] = -0.5*gc[0]*(1.0 - gc[2]); + devFunValue[4] = 0.5*(gc[1]+1)*(1.0 - gc[2]); + devFunValue[5] = -0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0); + + devFunValue[6] = 0.5*(gc[0]-1)*(1.0 - gc[2]); + devFunValue[7] = -0.5*gc[1]*(1.0 - gc[2]); + devFunValue[8] = -0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0); + + devFunValue[9] = -0.5*gc[0]*(1.0 - gc[2]); + devFunValue[10] = 0.5*(1*gc[1]-1)*(1.0 - gc[2]); + devFunValue[11] = -0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0); + + devFunValue[12] = 0.0; + devFunValue[13] = 0.0; + devFunValue[14] = 1.0; + + DEV_SHAPE_FUN_MACRO_END; } /*! * Init Pyramid Reference coordinates ans Shape function. @@ -1916,6 +1940,34 @@ void GaussInfo::penta6aInit() funValue[4] = 0.5*gc[2]*(gc[0] + 1.0); funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = 0.5*gc[1]*(-1.0); + devFunValue[1] = 0.5*(1.0 - gc[0]); + devFunValue[2] = 0.0; + + devFunValue[3] = 0.5*gc[2]*(-1.0); + devFunValue[4] = 0.0; + devFunValue[5] = 0.5*(1.0 - gc[0]); + + devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(-1.0); + devFunValue[7] = 0.5*(-1.0)*(1.0 - gc[0]); + devFunValue[8] = 0.5*(-1.0)*(1.0 - gc[0]); + + devFunValue[9] = 0.5*gc[1]; + devFunValue[10] = 0.5*(gc[0] + 1.0); + devFunValue[11] = 0.0; + + devFunValue[12] = 0.5*gc[2]; + devFunValue[13] = 0.0; + devFunValue[14] = 0.5*(gc[0] + 1.0); + + devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2]); + devFunValue[16] = 0.5*(-1.0)*(1.0 + gc[0]); + devFunValue[17] = 0.5*(-1.0)*(1.0 + gc[0]); + + DEV_SHAPE_FUN_MACRO_END; } /*! @@ -2168,6 +2220,70 @@ void GaussInfo::penta15aInit() funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = 0.5*gc[1]*(2 * gc[0] - 2 * gc[1] + 1.0 ); + devFunValue[1] = (1.0 - gc[0])*(2.0*gc[1] -1 - 0.5*gc[0]); + devFunValue[2] = 0.0; + + devFunValue[3] = 0.5*gc[2]*(2 * gc[0] - 2 * gc[2] + 1.0 ); + devFunValue[4] = 0.0; + devFunValue[5] = (1.0 - gc[0])*(2.0*gc[2] -1 - 0.5*gc[0]); + + devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(2*gc[0] -1.0 + 2*gc[1] + 2*gc[2]); + devFunValue[7] = 0.5*(gc[0] - 1.0)* ( -4*gc[1] -gc[0] -4*gc[2] + 2.0); + devFunValue[8] = 0.5*(gc[0] - 1.0)* ( -4*gc[2] -4*gc[1] -gc[0] + 2.0); + + devFunValue[9] = 0.5*gc[1]*( 2*gc[0] + 2*gc[1] -1.0 ); + devFunValue[10] = (1.0 + gc[0])*(2.0*gc[1] - 1.0 + 0.5*gc[0]); + devFunValue[11] = 0.0; + + devFunValue[12] = 0.5*gc[2]*( 2*gc[0] + 2*gc[2] -1.0 ); + devFunValue[13] = 0.0; + devFunValue[14] = (1.0 + gc[0])*(2.0*gc[2] - 1.0 + 0.5*gc[0]); + + devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2])*(2 * gc[0] - 2.0 * gc[1] - 2.0 * gc[2] + 1.0 ); + devFunValue[16] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[1] + gc[0] - 4.0 * gc[2] + 2.0) ; + devFunValue[17] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[2] - 4*gc[1] + gc[0] + 2.0 ); + + devFunValue[18] = - 2.0*gc[1]*gc[2]; + devFunValue[19] = 2.0*gc[2]*(1.0 - gc[0]); + devFunValue[20] = 2.0*gc[1]*(1.0 - gc[0]); + + devFunValue[21] = -2.0*gc[2]*(1.0 - gc[1] - gc[2]); + devFunValue[22] = -2.0*gc[2]*(1.0 - gc[0]); + devFunValue[23] = (-4*gc[2]-2*gc[1]+2.0)*(1.0 - gc[0]); + + devFunValue[24] = -2.0*gc[1]*(1.0 - gc[1] - gc[2]); + devFunValue[25] = (-4*gc[1]-2*gc[2]+2)*(1.0 - gc[0]); + devFunValue[26] = -2.0*gc[1]*(1.0 - gc[0]); + + devFunValue[27] = gc[1]*(- 2.0 * gc[0]); + devFunValue[28] = (1.0 - gc[0]*gc[0]); + devFunValue[29] = 0.0; + + devFunValue[30] = -2.0 * gc[2] *gc[0]; + devFunValue[31] = 0.0; + devFunValue[32] = (1.0 - gc[0]*gc[0]); + + devFunValue[33] = -2.0 * (1.0 - gc[1] - gc[2]) * gc[0]; + devFunValue[34] = -1.0*(1.0 - gc[0]*gc[0]); + devFunValue[35] = -1.0*(1.0 - gc[0]*gc[0]); + + devFunValue[36] = 2.0*gc[1]*gc[2]; + devFunValue[37] = 2.0*gc[2]*(1.0 + gc[0]); + devFunValue[38] = 2.0*gc[1]*(1.0 + gc[0]); + + devFunValue[39] = 2.0*gc[2]*(1.0 - gc[1] - gc[2]); + devFunValue[40] = -2.0*gc[2]*(1.0 + gc[0]); + devFunValue[41] = (2.0 - 2.0 * gc[1] - 4.0 * gc[2])*(1.0 + gc[0]); + + devFunValue[42] = 2.0*gc[1]*(1.0 - gc[1] - gc[2]); + devFunValue[43] = (2.0 - 4*gc[1] - 2*gc[2])*(1.0 + gc[0]); + devFunValue[44] = -2.0*gc[1]*(1.0 + gc[0]); + + DEV_SHAPE_FUN_MACRO_END; } /*! @@ -2581,6 +2697,42 @@ void GaussInfo::hexa8aInit() funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = 0.125 * (-1.0) * (1.0 - gc[1])*(1.0 - gc[2]); + devFunValue[1] = 0.125 * (1.0 - gc[0]) * (-1.0) * (1.0 - gc[2]); + devFunValue[2] = 0.125 * (1.0 - gc[0]) * (1.0 - gc[1]) *(-1.0); + + devFunValue[3] = 0.125*(1.0 - gc[1])*(1.0 - gc[2]); + devFunValue[4] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 - gc[2]); + devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0); + + devFunValue[6] = 0.125*(1.0 + gc[1])*(1.0 - gc[2]); + devFunValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[2]); + devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0); + + devFunValue[9 ] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 - gc[2]); + devFunValue[10] = 0.125*(1.0 - gc[0])*(1.0 - gc[2]); + devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0); + + devFunValue[12] = 0.125*(-1.0)*(1.0 - gc[1])*(1.0 + gc[2]); + devFunValue[13] = 0.125*(1.0 - gc[0])*(-1.0)*(1.0 + gc[2]); + devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1]); + + devFunValue[15] = 0.125*(1.0 - gc[1])*(1.0 + gc[2]); + devFunValue[16] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 + gc[2]); + devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1]); + + devFunValue[18] = 0.125*(1.0 + gc[1])*(1.0 + gc[2]); + devFunValue[19] = 0.125*(1.0 + gc[0])*(1.0 + gc[2]); + devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1]); + + devFunValue[21] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 + gc[2]); + devFunValue[22] = 0.125*(1.0 - gc[0])*(1.0 + gc[2]); + devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1]); + + DEV_SHAPE_FUN_MACRO_END; } /*! diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index bd3edc096..36d144b4f 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1081,6 +1081,44 @@ class MEDCouplingBasicsTest7(unittest.TestCase): ref2_mc = DataArrayDouble(ref2) ; ref2_mc.rearrange(1) self.assertTrue( mvDevOfShapeFunc.isEqual( ref2_mc, 1e-12) ) + def testShapeFuncAndDerivative1(self): + """ + This test focus + """ + def GetShapeFunc(ref_coord,vec): + gl3 = MEDCouplingGaussLocalization(gt,sum(ref_coord,[]), vec, [1]) + funVal = gl3.getShapeFunctionValues() + funVal.rearrange(1) + return funVal + + def GetDerivative(ref_coord,vec): + gl3 = MEDCouplingGaussLocalization(gt,sum(ref_coord,[]), vec, [1]) + funVal = gl3.getDerivativeOfShapeFunctionValues() + return funVal + vec = [-0.85685375,-0.90643355,-0.90796825] + eps = 1e-6 + + for gt in [NORM_HEXA8,NORM_PENTA6,NORM_PYRA5,NORM_PENTA15,NORM_HEXA20,NORM_HEXA27]: # type of cell for which derivatives are implemented + ref_coord = [list(elt) for elt in MEDCouplingGaussLocalization.GetDefaultReferenceCoordinatesOf(gt).getValuesAsTuple()] + + der_computed = GetDerivative(ref_coord,vec) + der_computed.rearrange(3) + + der_deduced = ( GetShapeFunc(ref_coord,[vec[0]+eps,vec[1],vec[2]])-GetShapeFunc(ref_coord,vec) ) / eps + delta_X = der_computed[:,0]-der_deduced + delta_X.abs() + self.assertTrue(delta_X.findIdsNotInRange(-1e-5,+1e-5).empty()) + + der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1]+eps,vec[2]])-GetShapeFunc(ref_coord,vec) ) / eps + delta_Y = der_computed[:,1]-der_deduced + delta_Y.abs() + self.assertTrue(delta_Y.findIdsNotInRange(-1e-5,+1e-5).empty()) + + der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1],vec[2]+eps])-GetShapeFunc(ref_coord,vec) ) / eps + delta_Z = der_computed[:,2]-der_deduced + delta_Z.abs() + self.assertTrue(delta_Z.findIdsNotInRange(-1e-5,+1e-5).empty()) + pass if __name__ == '__main__': -- 2.39.2