From 105e922749ed7f81c8090c878f896244b856939e Mon Sep 17 00:00:00 2001
From: geay <anthony.geay@cea.fr>
Date: Wed, 12 Feb 2014 12:54:59 +0100
Subject: [PATCH] Addition of new reference coords including degenerated cells.

---
 .../GaussPoints/InterpKernelGaussCoords.cxx   | 443 +++++++++++++++++-
 .../GaussPoints/InterpKernelGaussCoords.hxx   |  23 +
 2 files changed, 442 insertions(+), 24 deletions(-)

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
index bca5fb2af..121a2c255 100644
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
@@ -187,13 +187,51 @@ bool GaussInfo::isSatisfy()
   return anIsSatisfy;
 }
 
+std::vector<double> GaussInfo::NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector<double>& inputArray)
+{
+  std::size_t sz(inputArray.size()),dim((std::size_t)inputDim);
+  if(dim==0)
+    throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid dimension ! Must be !=0 !");
+  if(sz%dim!=0)
+    throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid input array ! Inconsistent with the given dimension !");
+  const CellModel& cm(CellModel::GetCellModel(ct));
+  std::size_t baseDim((std::size_t)cm.getDimension());
+  if(baseDim==dim)
+    return inputArray;
+  std::size_t nbOfItems(sz/dim);
+  std::vector<double> ret(nbOfItems*baseDim);
+  if(baseDim>dim)
+    {
+      for(std::size_t i=0;i<nbOfItems;i++)
+        {
+          std::size_t j=0;
+          for(;j<dim;j++)
+            ret[i*baseDim+j]=inputArray[i*dim+j];
+          for(;j<baseDim;j++)
+            ret[i*baseDim+j]=0.;
+        }
+    }
+  else
+    {
+      for(std::size_t i=0;i<nbOfItems;i++)
+        {
+          std::size_t j=0;
+          for(;j<baseDim;j++)
+            ret[i*baseDim+j]=inputArray[i*dim+j];
+        }
+    }
+  return ret;
+}
+
+typedef void (*MapToShapeFunction)(GaussInfo& obj);
+
 /*!
  * Initialize the internal vectors
  */
 void GaussInfo::initLocalInfo()
 {
   bool aSatify = false;
-  const CellModel& cellModel=CellModel::GetCellModel(_my_geometry);
+  const CellModel& cellModel(CellModel::GetCellModel(_my_geometry));
   switch( _my_geometry ) 
     {
     case NORM_SEG2:
@@ -240,17 +278,19 @@ void GaussInfo::initLocalInfo()
       break;
 
     case NORM_QUAD4:
-      _my_local_ref_dim = 2;
-      _my_local_nb_ref  = 4;
-      quad4aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          quad4bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
+      {
+        _my_local_ref_dim = 2;
+        _my_local_nb_ref  = 4;
+        MapToShapeFunction QUAD4PTR[]={Quad4aInit,Quad4bInit,Quad4cInit,Quad4DegSeg2Init};
+        std::size_t NB_OF_QUAD4PTR(sizeof(QUAD4PTR)/sizeof(MapToShapeFunction));
+        for(std::size_t i=0;i<NB_OF_QUAD4PTR && !aSatify;i++)
+          {
+            (QUAD4PTR[i])(*this);
+            aSatify = isSatisfy();
+          }
+        CHECK_MACRO;
+        break;
+      }
       break;
 
     case NORM_QUAD8:
@@ -332,6 +372,21 @@ void GaussInfo::initLocalInfo()
       break;
 
     case NORM_PENTA6:
+      {
+        _my_local_ref_dim = 3;
+        _my_local_nb_ref  = 6;
+        MapToShapeFunction PENTA6PTR[]={Penta6aInit,Penta6bInit,Penta6DegTria3aInit,Penta6DegTria3bInit};
+        std::size_t NB_OF_PENTA6PTR(sizeof(PENTA6PTR)/sizeof(MapToShapeFunction));
+        for(std::size_t i=0;i<NB_OF_PENTA6PTR && !aSatify;i++)
+          {
+            (PENTA6PTR[i])(*this);
+            aSatify = isSatisfy();
+          }
+        CHECK_MACRO;
+        break;
+      }
+
+
       _my_local_ref_dim = 3;
       _my_local_nb_ref  = 6;
       penta6aInit();
@@ -360,18 +415,19 @@ void GaussInfo::initLocalInfo()
       break;
 
     case NORM_HEXA8:
-      _my_local_ref_dim = 3;
-      _my_local_nb_ref  = 8;
-      hexa8aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          hexa8bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
+      {
+        _my_local_ref_dim = 3;
+        _my_local_nb_ref  = 8;
+        MapToShapeFunction HEXA8PTR[]={Hexa8aInit,Hexa8bInit,Hexa8DegQuad4aInit,Hexa8DegQuad4bInit,Hexa8DegQuad4cInit};
+        std::size_t NB_OF_HEXA8PTR(sizeof(HEXA8PTR)/sizeof(MapToShapeFunction));
+        for(std::size_t i=0;i<NB_OF_HEXA8PTR && !aSatify;i++)
+          {
+            (HEXA8PTR[i])(*this);
+            aSatify = isSatisfy();
+          }
+        CHECK_MACRO;
+        break;
+      }
 
     case NORM_HEXA20:
       _my_local_ref_dim = 3;
@@ -661,6 +717,67 @@ void GaussInfo::quad4bInit()
    SHAPE_FUN_MACRO_END;
 }
 
+void GaussInfo::quad4cInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   break;
+
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
+   funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
+   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;
+}
+
+/*!
+ * This shapefunc map is same as degenerated seg2Init
+ */
+void GaussInfo::quad4DegSeg2Init()
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.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] =  0.0;
+   break;
+ case  3:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*(1.0 - gc[0]);
+   funValue[1] = 0.5*(1.0 + gc[0]);
+   funValue[2] = 0.;
+   funValue[3] = 0.;
+   SHAPE_FUN_MACRO_END;
+}
 
 /*!
  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
@@ -1440,6 +1557,103 @@ void GaussInfo::penta6bInit()
    funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
    SHAPE_FUN_MACRO_END;
 }
+
+/*!
+ * This shapefunc map is same as degenerated tria3aInit
+ */
+void GaussInfo::penta6DegTria3aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*(1.0 + gc[1]);
+   funValue[1] = -0.5*(gc[0] + gc[1]);
+   funValue[2] = 0.5*(1.0 + gc[0]);
+   funValue[3] = 0.;
+   funValue[4] = 0.;
+   funValue[5] = 0.;
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * This shapefunc map is same as degenerated tria3bInit
+ */
+void GaussInfo::penta6DegTria3bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 1.0 - gc[0] - gc[1];
+   funValue[1] = gc[0];
+   funValue[2] = gc[1];
+   funValue[3] = 0.;
+   funValue[4] = 0.;
+   funValue[5] = 0.;
+   SHAPE_FUN_MACRO_END;
+}
+
 /*!
  * Init Pentahedron Reference coordinates and Shape function.
  * Case A.
@@ -1780,6 +1994,187 @@ void GaussInfo::hexa8bInit()
    SHAPE_FUN_MACRO_END;
 }
 
+/*!
+ * This shapefunc map is same as degenerated quad4bInit
+ */
+void GaussInfo::hexa8DegQuad4aInit()
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  6:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  7:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+  LOCAL_COORD_MACRO_END;
+  
+  SHAPE_FUN_MACRO_BEGIN;
+  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[4] = 0.;
+  funValue[5] = 0.;
+  funValue[6] = 0.;
+  funValue[7] = 0.;
+  SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * This shapefunc map is same as degenerated quad4bInit
+ */
+void GaussInfo::hexa8DegQuad4bInit()
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  6:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  7:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
+   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
+   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
+   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
+   funValue[4] = 0.;
+   funValue[5] = 0.;
+   funValue[6] = 0.;
+   funValue[7] = 0.;
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * This shapefunc map is same as degenerated quad4cInit
+ */
+void GaussInfo::hexa8DegQuad4cInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  6:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  7:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
+   funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
+   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
+   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
+   funValue[4] = 0. ;
+   funValue[5] = 0. ;
+   funValue[6] = 0. ;
+   funValue[7] = 0. ;
+   SHAPE_FUN_MACRO_END;
+}
+
 /*!
  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
  * Case A.
diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
index 35194aecc..f05603096 100644
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
@@ -54,6 +54,8 @@ namespace INTERP_KERNEL
     INTERPKERNEL_EXPORT const double* getFunctionValues( const int theGaussId ) const;
 
     INTERPKERNEL_EXPORT void initLocalInfo();
+    
+    INTERPKERNEL_EXPORT static std::vector<double> NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector<double>& inputArray);
 
   protected:
 
@@ -70,7 +72,13 @@ namespace INTERP_KERNEL
     void tria6bInit();
 
     void quad4aInit();
+    static void Quad4aInit(GaussInfo& obj) { obj.quad4aInit(); }
     void quad4bInit();
+    static void Quad4bInit(GaussInfo& obj) { obj.quad4bInit(); }
+    void quad4cInit();
+    static void Quad4cInit(GaussInfo& obj) { obj.quad4cInit(); }
+    void quad4DegSeg2Init();
+    static void Quad4DegSeg2Init(GaussInfo& obj) { obj.quad4DegSeg2Init(); }
     void quad8aInit();
     void quad8bInit();
     void quad9aInit();
@@ -87,12 +95,27 @@ namespace INTERP_KERNEL
     void pyra13bInit();
 
     void penta6aInit();
+    static void Penta6aInit(GaussInfo& obj) { obj.penta6aInit(); }
     void penta6bInit();
+    static void Penta6bInit(GaussInfo& obj) { obj.penta6bInit(); }
+    void penta6DegTria3aInit();
+    static void Penta6DegTria3aInit(GaussInfo& obj) { obj.penta6DegTria3aInit(); }
+    void penta6DegTria3bInit();
+    static void Penta6DegTria3bInit(GaussInfo& obj) { obj.penta6DegTria3bInit(); }
+    
     void penta15aInit();
     void penta15bInit();
 
     void hexa8aInit();
+    static void Hexa8aInit(GaussInfo& obj) { obj.hexa8aInit(); }
     void hexa8bInit();
+    static void Hexa8bInit(GaussInfo& obj) { obj.hexa8bInit(); }
+    void hexa8DegQuad4aInit();
+    static void Hexa8DegQuad4aInit(GaussInfo& obj) { obj.hexa8DegQuad4aInit(); }
+    void hexa8DegQuad4bInit();
+    static void Hexa8DegQuad4bInit(GaussInfo& obj) { obj.hexa8DegQuad4bInit(); }
+    void hexa8DegQuad4cInit();
+    static void Hexa8DegQuad4cInit(GaussInfo& obj) { obj.hexa8DegQuad4cInit(); }
     void hexa20aInit();
     void hexa20bInit();
 
-- 
2.39.2