From 791f862310efa99bc01010d1cfc44c4c89df4970 Mon Sep 17 00:00:00 2001
From: geay <anthony.geay@cea.fr>
Date: Fri, 7 Mar 2014 15:49:54 +0100
Subject: [PATCH] Addition of TRI7 shape functions. And addition of management
 of NORM_POINT1 type.

---
 .../GaussPoints/InterpKernelGaussCoords.cxx   | 92 ++++++++++++++++---
 .../GaussPoints/InterpKernelGaussCoords.hxx   |  7 +-
 .../MEDCouplingFieldDiscretization.cxx        | 12 +++
 .../MEDCouplingFieldDiscretization.hxx        |  1 +
 src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 24 +++++
 5 files changed, 123 insertions(+), 13 deletions(-)

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
index 14f8af006..b049d6475 100644
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
@@ -234,6 +234,14 @@ void GaussInfo::initLocalInfo()
   const CellModel& cellModel(CellModel::GetCellModel(_my_geometry));
   switch( _my_geometry ) 
     {
+    case NORM_POINT1:
+      _my_local_ref_dim = 0;
+      _my_local_nb_ref  = 1;
+      point1Init();
+      aSatify = isSatisfy();
+      CHECK_MACRO;
+      break;
+
     case NORM_SEG2:
       _my_local_ref_dim = 1;
       _my_local_nb_ref  = 2;
@@ -277,6 +285,14 @@ void GaussInfo::initLocalInfo()
         }
       break;
 
+    case NORM_TRI7:
+      _my_local_ref_dim = 2;
+      _my_local_nb_ref  = 7;
+      tria7aInit();
+      aSatify = isSatisfy();
+      CHECK_MACRO;
+      break;
+
     case NORM_QUAD4:
       {
         _my_local_ref_dim = 2;
@@ -401,18 +417,19 @@ void GaussInfo::initLocalInfo()
       break;
 
     case NORM_PENTA15:
-      _my_local_ref_dim = 3;
-      _my_local_nb_ref  = 15;
-      penta15aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          penta15bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
+      {
+        _my_local_ref_dim = 3;
+        _my_local_nb_ref  = 15;
+        MapToShapeFunction PENTA15PTR[]={Penta15aInit,Penta15bInit};
+        std::size_t NB_OF_PENTA15PTR(sizeof(PENTA15PTR)/sizeof(MapToShapeFunction));
+        for(std::size_t i=0;i<NB_OF_PENTA15PTR && !aSatify;i++)
+          {
+            (PENTA15PTR[i])(*this);
+            aSatify = isSatisfy();
+          }
+        CHECK_MACRO;
+        break;
+      }
 
     case NORM_HEXA8:
       {
@@ -465,6 +482,12 @@ const double* GaussInfo::getFunctionValues( const int theGaussId ) const
   return &_my_function_value[ _my_nb_ref*theGaussId ];
 }
 
+void GaussInfo::point1Init()
+{
+  double *funValue(&_my_function_value[0]);
+  funValue[0] = 1. ;
+}
+
 /*!
  * Init Segment 2 Reference coordinates ans Shape function.
  */
@@ -657,6 +680,51 @@ void GaussInfo::tria6bInit()
    SHAPE_FUN_MACRO_END;
 }
 
+void GaussInfo::tria7aInit()
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case 0:
+   coords[0] = 0.0;
+   coords[1] = 0.0;
+   break;
+ case 1:
+   coords[0] = 1.0;
+   coords[1] = 0.0;
+   break;
+ case 2:
+   coords[0] = 0.0;
+   coords[1] = 1.0;
+   break;
+ case 3:
+   coords[0] = 0.5;
+   coords[1] = 0.0;
+   break;
+ case 4:
+   coords[0] = 0.5;
+   coords[1] = 0.5;
+   break;
+ case 5:
+   coords[0] = 0.0;
+   coords[1] = 0.5;
+   break;
+ case 6:
+   coords[0] = 0.3333333333333333;
+   coords[1] = 0.3333333333333333;
+   break;
+
+  LOCAL_COORD_MACRO_END;
+
+  SHAPE_FUN_MACRO_BEGIN;
+  funValue[0]=1-3*(gc[0]+gc[1])+2*(gc[0]*gc[0]+gc[1]*gc[1])+7*gc[0]*gc[1]-3*gc[0]*gc[1]*(gc[0]+gc[1]);
+  funValue[1]=gc[0]*(-1+2*gc[0]+3*gc[1]-3*gc[1]*(gc[0]+gc[1]));
+  funValue[2]=gc[1]*(-1.+3.*gc[0]+2.*gc[1]-3.*gc[0]*(gc[0]+gc[1]));
+  funValue[3]=4*gc[0]*(1-gc[0]-4*gc[1]+3*gc[1]*(gc[0]+gc[1]));
+  funValue[4]=4*gc[0]*gc[1]*(-2+3*(gc[0]+gc[1]));
+  funValue[5]=4*gc[1]*(1-4*gc[0]-gc[1]+3*gc[0]*(gc[0]+gc[1]));
+  funValue[6]=27*gc[0]*gc[1]*(1-gc[0]-gc[1]);
+  SHAPE_FUN_MACRO_END;
+}
+
 /*!
  * Init Quadrangle Reference coordinates ans Shape function.
  * Case A.
diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
index 465f191fd..07265c65a 100644
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
@@ -60,7 +60,9 @@ namespace INTERP_KERNEL
   protected:
 
     bool isSatisfy();
-
+    
+    void point1Init();
+    
     //1D
     void seg2Init();
     void seg3Init();
@@ -70,6 +72,7 @@ namespace INTERP_KERNEL
     void tria3bInit();
     void tria6aInit();
     void tria6bInit();
+    void tria7aInit();
 
     void quad4aInit();
     static void Quad4aInit(GaussInfo& obj) { obj.quad4aInit(); }
@@ -104,7 +107,9 @@ namespace INTERP_KERNEL
     static void Penta6DegTria3bInit(GaussInfo& obj) { obj.penta6DegTria3bInit(); }
     
     void penta15aInit();
+    static void Penta15aInit(GaussInfo& obj) { obj.penta15aInit(); }
     void penta15bInit();
+    static void Penta15bInit(GaussInfo& obj) { obj.penta15bInit(); }
 
     void hexa8aInit();
     static void Hexa8aInit(GaussInfo& obj) { obj.hexa8aInit(); }
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 65a48342c..4a92d2c8f 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -65,6 +65,7 @@ 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.};
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG3[3]={0.5555555555555556,0.8888888888888888,0.5555555555555556};
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG4[4]={0.347854845137454,0.347854845137454,0.652145154862546,0.652145154862546};
@@ -2269,6 +2270,9 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometric
 {
   switch(geoType)
     {
+    case INTERP_KERNEL::NORM_POINT1:
+      lgth=(int)sizeof(FGP_POINT1)/sizeof(double);
+      return FGP_POINT1;
     case INTERP_KERNEL::NORM_SEG2:
       lgth=(int)sizeof(FGP_SEG2)/sizeof(double);
       return FGP_SEG2;
@@ -2320,6 +2324,9 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetRefCoordsFromGeometricTy
 {
   switch(geoType)
     {
+    case INTERP_KERNEL::NORM_POINT1:
+      lgth=0;
+      return 0;
     case INTERP_KERNEL::NORM_SEG2:
       lgth=(int)sizeof(REF_SEG2)/sizeof(double);
       return REF_SEG2;
@@ -2383,6 +2390,11 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetLocsFromGeometricType(IN
 {
   switch(geoType)
     {
+    case INTERP_KERNEL::NORM_POINT1:
+      {
+	lgth=0;
+	return 0;
+      }
     case INTERP_KERNEL::NORM_SEG2:
       {
         lgth=(int)sizeof(LOC_SEG2)/sizeof(double);
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
index 26f31d952..9c73aa5e6 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
@@ -335,6 +335,7 @@ namespace ParaMEDMEM
   public:
     static const char REPR[];
     static const TypeOfField TYPE;
+    static const double FGP_POINT1[1];
     static const double FGP_SEG2[2];
     static const double FGP_SEG3[3];
     static const double FGP_SEG4[4];
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index de8113d3f..fd9391764 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -14491,6 +14491,30 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertTrue(arrOfDisc2.isEqual(coo,1e-10)) # be less exigent 1e-10 instead of 1e-12 due to shape function sensitivity arount 0.,0.,1. !
         pass
 
+    def testSwig2Tri7GP1(self):
+        coo=DataArrayDouble([[0,0],[0,2],[2,0],[0,1],[1,1],[1,0],[0.6666666666666667,0.6666666666666667]])
+        m=MEDCouplingUMesh("mesh",2) ; m.setCoords(coo)
+        m.allocateCells()
+        # the cell description is exactly those described in the description of TRI7 in MED file 3.0.7 documentation
+        m.insertNextCell(NORM_TRI7,range(7))
+        refCoords=[0.,0.,1.,0.,0.,1.,0.5,0.,0.5,0.5,0.,0.5,0.3333333333333333,0.3333333333333333]
+        gaussCoords=[0.3333333333333333,0.3333333333333333,0.470142064105115,0.470142064105115,0.05971587178977,0.470142064105115,0.470142064105115,0.05971587178977,0.101286507323456,0.101286507323456,0.797426985353088,0.101286507323456,0.101286507323456,0.797426985353088]
+        weights=[0.062969590272413,0.062969590272413,0.062969590272413,0.066197076394253,0.066197076394253,0.066197076394253,0.1125]
+        fGauss=MEDCouplingFieldDouble(ON_GAUSS_PT) ; fGauss.setName("fGauss")
+        fGauss.setMesh(m)
+        fGauss.setGaussLocalizationOnType(NORM_TRI7,refCoords,gaussCoords,weights)
+        arr=DataArrayDouble(fGauss.getNumberOfTuplesExpected()) ; arr.iota()
+        fGauss.setArray(arr)
+        arrOfDisc=fGauss.getLocalizationOfDiscr()
+        self.assertTrue(arrOfDisc.isEqual(DataArrayDouble([0.666666666666667,0.666666666666667,0.9402841282102293,0.9402841282102293,0.9402841282102299,0.11943174357954002,0.11943174357953992,0.9402841282102299,0.20257301464691194,0.20257301464691196,0.20257301464691205,1.5948539707061757,1.5948539707061757,0.20257301464691202],7,2),1e-12))
+        #
+        weights=7*[1]
+        gaussCoords=refCoords
+        fGauss.setGaussLocalizationOnType(NORM_TRI7,refCoords,gaussCoords,weights)
+        arrOfDisc2=fGauss.getLocalizationOfDiscr()
+        self.assertTrue(arrOfDisc2.isEqual(coo,1e-12))
+        pass
+
     def setUp(self):
         pass
     pass
-- 
2.39.2