From: Anthony Geay <anthony.geay@edf.fr>
Date: Tue, 24 Jan 2023 07:45:50 +0000 (+0100)
Subject: WIP
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=2123f76a368328bb8b1413d2f5eef78b41d93924;p=tools%2Fmedcoupling.git

WIP
---

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
index 657e4a898..fbb224eeb 100644
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
@@ -1103,8 +1103,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;
 }
 
 /*!
@@ -1253,6 +1269,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;
 }
 
 /*!
diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
index 38e4a6128..63080077f 100644
--- a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
+++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
@@ -47,6 +47,7 @@ namespace INTERP_KERNEL
     void resize(mcIdType newn, mcIdType newm);
     void assign(mcIdType newn, mcIdType newm, const T &a);
     ~DenseMatrixT();
+    T determinant() const;
   };
 
   using DenseMatrix = DenseMatrixT<double>;
diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
index 9705f4fdb..2a262c31b 100644
--- a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
+++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
@@ -20,6 +20,7 @@
 #pragma once
 
 #include "InterpKernelDenseMatrix.hxx"
+#include "InterpKernelException.hxx"
 
 namespace INTERP_KERNEL
 {
@@ -138,4 +139,22 @@ namespace INTERP_KERNEL
       delete[] (v);
     }
   }
+
+  template<class T>
+  T Determinant22(const T *m)
+  { return m[0]*m[3]-m[1]*m[2]; }
+
+  template<class T>
+  T Determinant33(const T *m)
+  { return m[0]*(m[4]*m[8]-m[7]*m[5])-m[1]*(m[3]*m[8]-m[6]*m[5])+m[2]*(m[3]*m[7]-m[6]*m[4]);}
+  
+  template <class T>
+  T DenseMatrixT<T>::determinant() const
+  {
+    if(nn==2 && mm==2)
+      return Determinant22(v[0]);
+    if(nn==3 && mm==3)
+      return Determinant33(v[0]);
+    THROW_IK_EXCEPTION("DenseMatrixT::determinant : only 2x2 and 3x3 implemented !");
+  }
 }
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index bf1397a4c..9fdf4ece7 100755
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -1745,8 +1745,8 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
           const MEDCouplingGaussLocalization& loc=_loc[locId];
           mcIdType nbOfGaussPt=loc.getNumberOfGaussPt();
           INTERP_KERNEL::AutoPtr<double> weights=new double[nbOfGaussPt];
-          double sum=std::accumulate(loc.getWeights().begin(),loc.getWeights().end(),0.);
-          std::transform(loc.getWeights().begin(),loc.getWeights().end(),(double *)weights,std::bind(std::multiplies<double>(),std::placeholders::_1,1./sum));
+          //double sum=std::accumulate(loc.getWeights().begin(),loc.getWeights().end(),0.);
+          //std::transform(loc.getWeights().begin(),loc.getWeights().end(),(double *)weights,std::bind(std::multiplies<double>(),std::placeholders::_1,1./sum));
           for(const mcIdType *cellId=curIds->begin();cellId!=curIds->end();cellId++)
           {
             std::vector<mcIdType> conn;
@@ -1756,18 +1756,18 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
             std::size_t nbPtsInCell(ptsInCell.size()/3);
             INTERP_KERNEL::DenseMatrix jacobian(2,2);
             MCAuto<DataArrayDouble> shapeFunc = loc.getDerivativeOfShapeFunctionValues();
-            for(std::size_t i = 0 ; i < 2 ; ++i)
-              for(std::size_t j = 0 ; j < 2 ; ++j)
-              {
-                double res = 0.0;
-                for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
-                  res += ptsInCell[k*2+i] * shapeFunc->getIJ(0,2*k+j);
-                jacobian[ i ][ j ] = res;
-              }
-
-
-            for(mcIdType j=0;j<nbOfGaussPt;j++)
-              arrPtr[offsetPtr[*cellId]+j]=weights[j];//*volPtr[*cellId]
+            for(mcIdType iGPt = 0 ; iGPt < nbOfGaussPt ; ++iGPt)
+            {
+              for(std::size_t i = 0 ; i < 2 ; ++i)
+                for(std::size_t j = 0 ; j < 2 ; ++j)
+                {
+                  double res = 0.0;
+                  for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
+                    res += ptsInCell[k*2+i] * shapeFunc->getIJ(iGPt,2*k+j);
+                  jacobian[ i ][ j ] = res;
+                }
+              arrPtr[offsetPtr[*cellId]+iGPt]=jacobian.determinant();
+            }
           }
         }
       else