[EDF26877] : Improvement of precision of integral computation of fields on Gauss...

[tools/medcoupling.git] / src / INTERP_KERNEL / GaussPoints / InterpKernelGaussCoords.cxx

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
index 657e4a898befd9ab820f4e3b9d3d0abcdf827b49..05435efe39dfecdc8d835ce90c47784cb504d109 100644 (file)
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
@@ -25,6 +25,7 @@
 #include <math.h>
 #include <algorithm>
 #include <sstream>
+#include <cmath>
 
 using namespace INTERP_KERNEL;
 
@@ -502,8 +503,9 @@ std::vector<double> GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellTy
       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));
+    default:
+      THROW_IK_EXCEPTION("Input type " << ct << "is not managed by GetDefaultReferenceCoordinatesOf")
   }
-  THROW_IK_EXCEPTION("Input type " << ct << "is not managed by GetDefaultReferenceCoordinatesOf")
 }
 
 /*!
@@ -914,6 +916,19 @@ void GaussInfo::tria3aInit()
   funValue[1] = -0.5*(gc[0] + gc[1]);
   funValue[2] = 0.5*(1.0 + gc[0]);
   SHAPE_FUN_MACRO_END;
+  
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = 0.0 ;
+  devFunValue[1] = 0.5 ;
+
+  devFunValue[2] = -0.5;
+  devFunValue[3] = -0.5;
+
+  devFunValue[4] = 0.5;
+  devFunValue[5] = 0.0;
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 /*!
@@ -942,6 +957,19 @@ void GaussInfo::tria3bInit()
   funValue[1] = gc[0];
   funValue[2] = gc[1];
   SHAPE_FUN_MACRO_END;
+
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = -1.0 ;
+  devFunValue[1] = -1.0 ;
+
+  devFunValue[2] = 1.0 ;
+  devFunValue[3] = 0.0 ;
+
+  devFunValue[4] = 0.0 ;
+  devFunValue[5] = 1.0 ;
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 /*!
@@ -985,6 +1013,28 @@ void GaussInfo::tria6aInit()
   funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
   funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
   SHAPE_FUN_MACRO_END;
+  
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = 0.0;
+  devFunValue[1] = 0.5*( 2*gc[1] + 1.0 );
+
+  devFunValue[2] = 0.5*( 2*gc[0] + 2.0*gc[1] + 1.0);
+  devFunValue[3] = 0.5*( 2*gc[1] + 2.0*gc[0] + 1.0);
+
+  devFunValue[4] = gc[0] + 0.5;
+  devFunValue[5] = 0.0;
+
+  devFunValue[6] = -1.0*(1.0 + gc[1]);
+  devFunValue[7] = -1.0*(2*gc[1]+gc[0]+1.0);
+
+  devFunValue[8] = -1.0*(2*gc[0]+gc[1]+1.0);
+  devFunValue[9] = -1.0*(1.0 + gc[0]);
+
+  devFunValue[10] = 0.0;
+  devFunValue[11] = (2*gc[1]+2.0);
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 /*!
@@ -1028,6 +1078,28 @@ void GaussInfo::tria6bInit()
   funValue[4] = 4.0*gc[0]*gc[1];
   funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
   SHAPE_FUN_MACRO_END;
+
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = 4*gc[0] + 4*gc[1] - 3.0 ;
+  devFunValue[1] = 4*gc[1] + 4*gc[0] - 3.0 ;
+
+  devFunValue[2] = 4*gc[0] - 1.0 ;
+  devFunValue[3] = 0.0 ;
+
+  devFunValue[4] = 0.0 ;
+  devFunValue[5] = 4*gc[1] - 1.0 ;
+
+  devFunValue[6] = -8.0*gc[0] - 4.0 * gc[1] + 4.0 ;
+  devFunValue[7] = -4.0*gc[0] ; 
+
+  devFunValue[8] = 4.0*gc[1];
+  devFunValue[9] = 4.0*gc[0];
+
+  devFunValue[10] = -4.0*gc[1] ;
+  devFunValue[11] = -8.0*gc[1] - 4.0*gc[0] + 4.0 ; 
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 void GaussInfo::tria7aInit()
@@ -1103,8 +1175,24 @@ void GaussInfo::quad4aInit()
   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
-  funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
+  funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0]);
   SHAPE_FUN_MACRO_END;
+  
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = -0.25*(1.0 + gc[1]);
+  devFunValue[1] = 0.25*(1.0 - gc[0]);
+
+  devFunValue[2] = -0.25*(1.0 - gc[1]);
+  devFunValue[3] = -0.25*(1.0 - gc[0]);
+
+  devFunValue[4] = 0.25*(1.0 - gc[1]);
+  devFunValue[5] = -0.25*(1.0 + gc[0]);
+
+  devFunValue[6] = 0.25*(1.0 + gc[1]);
+  devFunValue[7] = 0.25*(1.0 + gc[0]);
+  
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 /*!
@@ -1138,6 +1226,22 @@ void GaussInfo::quad4bInit()
   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
   SHAPE_FUN_MACRO_END;
+  
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = -0.25*(1.0 - gc[1]);
+  devFunValue[1] = -0.25*(1.0 - gc[0]);
+
+  devFunValue[2] = 0.25*(1.0 - gc[1]);
+  devFunValue[3] = -0.25*(1.0 + gc[0]);
+
+  devFunValue[4] = 0.25*(1.0 + gc[1]);
+  devFunValue[5] = 0.25*(1.0 + gc[0]);
+
+  devFunValue[6] = -0.25*(1.0 + gc[1]);
+  devFunValue[7] = 0.25*(1.0 - gc[0]);
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 void GaussInfo::quad4cInit() 
@@ -1168,6 +1272,22 @@ void GaussInfo::quad4cInit()
    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
    funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
    SHAPE_FUN_MACRO_END;
+   
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = -0.25*(1.0 - gc[1]);
+  devFunValue[1] = -0.25*(1.0 - gc[0]);
+
+  devFunValue[2] = -0.25*(1.0 + gc[1]);
+  devFunValue[3] = 0.25*(1.0 - gc[0]);
+
+  devFunValue[4] = 0.25*(1.0 + gc[0]);
+  devFunValue[5] = 0.25*(1.0 + gc[1]);
+
+  devFunValue[6] = 0.25*(1.0 - gc[1]);
+  devFunValue[7] = -0.25*(1.0 + gc[0]);
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 /*!
@@ -1253,6 +1373,34 @@ void GaussInfo::quad8aInit()
   funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
   funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
   SHAPE_FUN_MACRO_END;
+
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = 0.25*(1.0 + gc[1])*(2*gc[0]-gc[1]);
+  devFunValue[1] = 0.25*(1.0 - gc[0])*(2*gc[1]-gc[0]);
+
+  devFunValue[2] = 0.25*(1.0 - gc[1])*(2*gc[0]+gc[1]);
+  devFunValue[3] = 0.25*(1.0 - gc[0])*(2*gc[1]+gc[0]);
+
+  devFunValue[4] = 0.25*(1.0 - gc[1])*(2*gc[0]-gc[1]);
+  devFunValue[5] = 0.25*(1.0 + gc[0])*(2*gc[1]-gc[0]);
+
+  devFunValue[6] = 0.25*(1.0 + gc[1])*(2*gc[0]+gc[1]);
+  devFunValue[7] = 0.25*(1.0 + gc[0])*(2*gc[1]+gc[0]);
+
+  devFunValue[8] = -0.5*(1.0 - gc[1])*(1.0 + gc[1]);
+  devFunValue[9] = 0.5*(1.0 - gc[0])*(-2*gc[1]);
+
+  devFunValue[10] = 0.5*(1.0 - gc[1])*(-2*gc[0]);
+  devFunValue[11] = -0.5*(1.0 - gc[0])*(1.0 + gc[0]);
+
+  devFunValue[12] = 0.5*(1.0 - gc[1])*(1.0 + gc[1]);
+  devFunValue[13] = 0.5*(1.0 + gc[0])*(-2*gc[1]);
+
+  devFunValue[14] = 0.5*(1.0 + gc[1])*(-2*gc[0]);
+  devFunValue[15] = 0.5*(1.0 - gc[0])*(1.0 + gc[0]);
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 /*!
@@ -1306,6 +1454,34 @@ void GaussInfo::quad8bInit()
   funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
   funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
   SHAPE_FUN_MACRO_END;
+  
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = 0.25*(1.0 - gc[1])*(2*gc[0] + gc[1] );
+  devFunValue[1] = 0.25*(1.0 - gc[0])*(2*gc[1] + gc[0] );
+
+  devFunValue[2] = 0.25*(1.0 - gc[1])*(2*gc[0] - gc[1] );
+  devFunValue[3] = 0.25*(1.0 + gc[0])*(2*gc[1] - gc[0] );
+
+  devFunValue[4] = 0.25*(1.0 + gc[1])*(2*gc[0] + gc[1]);
+  devFunValue[5] = 0.25*(1.0 + gc[0])*(2*gc[1] + gc[0]);
+
+  devFunValue[6] = 0.25*( 1.0 + gc[1])*(2*gc[0] - gc[1] );
+  devFunValue[7] = 0.25*(1.0 - gc[0])*(2*gc[1] - gc[0]);
+
+  devFunValue[8] = -(1.0 - gc[1])*gc[0];
+  devFunValue[9] = -0.5*(1.0 - gc[0]*gc[0]);
+
+  devFunValue[10] = 0.5*(1.0 - gc[1]*gc[1]);
+  devFunValue[11] = -(1.0 + gc[0])*gc[1];
+
+  devFunValue[12] = -(1.0 + gc[1])*gc[0];
+  devFunValue[13] = 0.5*(1.0 - gc[0]*gc[0]);
+
+  devFunValue[14] = -0.5*(1.0 - gc[1]*gc[1]);
+  devFunValue[15] = -(1.0 - gc[0])*gc[1];
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 void GaussInfo::quad9aInit()