From: Anthony Geay <anthony.geay@edf.fr>
Date: Thu, 8 Sep 2022 06:24:40 +0000 (+0200)
Subject: Management of derivatives of shape functions of NORM_HEXA8,NORM_PENTA6,NORM_PYRA5... 
X-Git-Tag: V9_10_0a1~6
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=083874a866268323aa868f9e8b383bbae8851725;p=tools%2Fmedcoupling.git

Management of derivatives of shape functions of NORM_HEXA8,NORM_PENTA6,NORM_PYRA5,NORM_PENTA15,NORM_HEXA20,NORM_HEXA27
---

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__':