[EDF26877] : management of computation of measure field on field on Gauss Point in...

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
index 05435efe39dfecdc8d835ce90c47784cb504d109..72c49b896bf6bd79fcf8edfaddbc893ee7c58d7b 100644 (file)
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
@@ -848,6 +848,14 @@ void GaussInfo::seg2aInit()
   funValue[0] = 0.5*(1.0 - gc[0]);
   funValue[1] = 0.5*(1.0 + gc[0]);
   SHAPE_FUN_MACRO_END;
+  
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = -0.5;
+
+  devFunValue[1] = 0.5;
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 void GaussInfo::seg2bInit() 
@@ -865,6 +873,14 @@ void GaussInfo::seg2bInit()
   funValue[0] = 1.0 - gc[0];
   funValue[1] = gc[0];
   SHAPE_FUN_MACRO_END;
+  
+  DEV_SHAPE_FUN_MACRO_BEGIN;
+
+  devFunValue[0] = -1.0 ;
+
+  devFunValue[1] = 1.0 ;
+
+  DEV_SHAPE_FUN_MACRO_END;
 }
 
 /*!

diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
index 63080077f3c0d39c472f6f4e8d3255987177f3bc..f363e21222ed33c3f111bdf074dd712bbb04a71c 100644 (file)
--- a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
+++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
@@ -48,6 +48,7 @@ namespace INTERP_KERNEL
     void assign(mcIdType newn, mcIdType newm, const T &a);
     ~DenseMatrixT();
     T determinant() const;
+    T toJacobian() const;
   };
 
   using DenseMatrix = DenseMatrixT<double>;

diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
index 2a262c31b14648b9a9ba5aecb3f736040bb2f609..3ae1aa2c30b706714cde61f480b540049a595a99 100644 (file)
--- a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
+++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
@@ -21,6 +21,9 @@
 
 #include "InterpKernelDenseMatrix.hxx"
 #include "InterpKernelException.hxx"
+#include "VectorUtils.hxx"
+
+#include <cmath>
 
 namespace INTERP_KERNEL
 {
@@ -151,10 +154,33 @@ namespace INTERP_KERNEL
   template <class T>
   T DenseMatrixT<T>::determinant() const
   {
+    if(nn==1 && mm==1)
+      return v[0][0];
     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 !");
+    THROW_IK_EXCEPTION("DenseMatrixT::determinant : only 1x1, 2x2 and 3x3 implemented !");
+  }
+  
+  template <class T>
+  T DenseMatrixT<T>::toJacobian() const
+  {
+    if(nn == mm)
+      return determinant();
+    const T *vPtr(this->v[0]);
+    if(nn==3 && mm==1)
+      return norm(vPtr);
+    if(nn==2 && mm==1)
+      return std::sqrt(vPtr[0]*vPtr[0] + vPtr[1]*vPtr[1]);
+    if(nn==3 && mm==2)
+    {
+      T tmp[3];
+      T VA[3] = {v[0][0],v[1][0],v[2][0]};
+      T VB[3] = {v[0][1],v[1][1],v[2][1]};
+      cross(VA,VB,tmp);
+      return norm(tmp);
+    }
+    THROW_IK_EXCEPTION("DenseMatrixT::toJacobian : only 1x1, 2x1, 3x1, 3x2, 2x2 and 3x3 implemented !");
   }
 }

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index ebf44b23da8a165f829a53e1c80ec9a23096c851..d35dccf516c2aef0eaaee754a414ba97fb2c23a7 100755 (executable)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -1715,8 +1715,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
   MCAuto<MEDCouplingUMesh> umesh(mesh->buildUnstructured());
   const double *coordsOfMesh( umesh->getCoords()->begin() );
   auto spaceDim(mesh->getSpaceDimension());
-  if( spaceDim != mesh->getMeshDimension())
-    THROW_IK_EXCEPTION("MEDCouplingFieldDiscretizationGauss::getMeasureField : only implemented when space dimension and mesh dimension are equal !");
+  auto meshDim(mesh->getMeshDimension());
   MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT);
   ret->setMesh(mesh);
   ret->setDiscretization(const_cast<MEDCouplingFieldDiscretizationGauss *>(this));
@@ -1752,19 +1751,19 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
             std::vector<double> ptsInCell; ptsInCell.reserve(conn.size()*loc.getDimension());
             std::for_each( conn.cbegin(), conn.cend(), [spaceDim,coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*spaceDim,coordsOfMesh+(c+1)*spaceDim); } );
             std::size_t nbPtsInCell(ptsInCell.size()/spaceDim);
-            INTERP_KERNEL::DenseMatrix jacobian(spaceDim,spaceDim);
+            INTERP_KERNEL::DenseMatrix jacobian(spaceDim,meshDim);
             MCAuto<DataArrayDouble> shapeFunc = loc.getDerivativeOfShapeFunctionValues();
             for(mcIdType iGPt = 0 ; iGPt < nbOfGaussPt ; ++iGPt)
             {
               for(auto i = 0 ; i < spaceDim ; ++i)
-                for(auto j = 0 ; j < spaceDim ; ++j)
+                for(auto j = 0 ; j < meshDim ; ++j)
                 {
                   double res = 0.0;
                   for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
-                    res += ptsInCell[spaceDim*k+i] * shapeFunc->getIJ(iGPt,spaceDim*k+j);
+                    res += ptsInCell[spaceDim*k+i] * shapeFunc->getIJ(iGPt,meshDim*k+j);
                   jacobian[ i ][ j ] = res;
                 }
-              arrPtr[offsetPtr[*cellId]+iGPt]=std::abs( jacobian.determinant() )*loc.getWeight(FromIdType<int>(iGPt));
+              arrPtr[offsetPtr[*cellId]+iGPt]=std::abs( jacobian.toJacobian() )*loc.getWeight(FromIdType<int>(iGPt));
             }
           }
         }

diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py
index d8815278c38319232d41d32e93bd11229c6e0feb..89c249da6585bf60d1204039d9b3ff1e37b71523 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py
@@ -1175,6 +1175,36 @@ class MEDCouplingBasicsTest7(unittest.TestCase):
         self.assertTrue( ret.isEqual(DataArrayInt([1, 1, 0, 1, 2, 1, 0, 1, 2]) ) )
         pass
 
+    def testMeasureOnGaussPtMeshDimNotEqualSpaceDim0(self):
+        """
+        [EDF26877] : This test focuses on computation of measure field on field on Gauss Point in the special case where SpaceDim
+        are not eqaul to the meshDim.
+        """
+        seg2 = MEDCouplingUMesh("mesh",1)
+        seg2.setCoords(DataArrayDouble([(3,3),(4,4)]))
+        seg2.allocateCells()
+        seg2.insertNextCell(NORM_SEG2,[0,1])
+        fff=MEDCouplingFieldDouble.New(ON_GAUSS_PT) ; fff.setName("CH1RB") ; fff.setNature(IntensiveMaximum)
+        fff.setMesh(seg2)
+        fff.setGaussLocalizationOnCells([0], [0.,1.], [0.333333333333333], [1.0])
+        disc = fff.getDiscretization()
+        # spaceDim = 2 meshDim = 1
+        self.assertTrue( disc.getMeasureField(seg2,True).getArray().isEqual(DataArrayDouble([sqrt(2.0)]),1e-10) )
+        # spaceDim = 3 meshDim = 1
+        seg2.setCoords(DataArrayDouble([(3,3,3),(4,4,4)]))
+        self.assertTrue( disc.getMeasureField(seg2,True).getArray().isEqual(DataArrayDouble([sqrt(3.0)]),1e-10) )
+        # spaceDim = 3 meshDim = 2
+        tri = MEDCouplingUMesh("mesh",2)
+        tri.setCoords( DataArrayDouble([(0,0,0),(1,1,0),(2,2,2)]) )
+        tri.allocateCells()
+        tri.insertNextCell(NORM_TRI3,[0,1,2])
+        fff=MEDCouplingFieldDouble.New(ON_GAUSS_PT) ; fff.setName("CH1RB") ; fff.setNature(IntensiveMaximum)
+        fff.setMesh(tri)
+        fff.setGaussLocalizationOnCells(list(range(0, 1)), [0., 0., 1., 0., 0., 1.], [0.3333333333333333, 0.3333333333333333], [0.5])
+        disc = fff.getDiscretization()
+        self.assertTrue( disc.getMeasureField(tri,True).getArray().isEqual( tri.getMeasureField(True).getArray(), 1e-10) )
+        pass
+
     pass
 
 if __name__ == '__main__':