refactor!: remove adm_local/ directory

[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingFieldDiscretization.cxx

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index d7cf86b66e9e3ed794092e7c3af5237daba5f18c..1e50212834c1af5cc6f15e83b60dd0cb5b25c227 100755 (executable)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2022  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024  CEA, EDF
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -30,6 +30,7 @@
 #include "InterpKernelAutoPtr.hxx"
 #include "InterpKernelGaussCoords.hxx"
 #include "InterpKernelMatrixTools.hxx"
+#include "InterpKernelDenseMatrix.hxx"
 
 #include <set>
 #include <list>
@@ -1711,8 +1712,10 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
 {
   if(!mesh)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : mesh instance specified is NULL !");
-  MCAuto<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isAbs);
-  const double *volPtr=vol->getArray()->begin();
+  MCAuto<MEDCouplingUMesh> umesh(mesh->buildUnstructured());
+  const double *coordsOfMesh( umesh->getCoords()->begin() );
+  auto spaceDim(mesh->getSpaceDimension());
+  auto meshDim(mesh->getMeshDimension());
   MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT);
   ret->setMesh(mesh);
   ret->setDiscretization(const_cast<MEDCouplingFieldDiscretizationGauss *>(this));
@@ -1721,7 +1724,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
   _discr_per_cell->checkAllocated();
   if(_discr_per_cell->getNumberOfComponents()!=1)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array defined but with nb of components different from 1 !");
-  if(_discr_per_cell->getNumberOfTuples()!=vol->getNumberOfTuples())
+  if(_discr_per_cell->getNumberOfTuples()!=mesh->getNumberOfCells())
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array defined but mismatch between nb of cells of mesh and size of spatial disr array !");
   MCAuto<DataArrayIdType> offset=getOffsetArr(mesh);
   MCAuto<DataArrayDouble> arr=DataArrayDouble::New(); arr->alloc(getNumberOfTuples(mesh),1);
@@ -1741,11 +1744,28 @@ 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));
           for(const mcIdType *cellId=curIds->begin();cellId!=curIds->end();cellId++)
-            for(mcIdType j=0;j<nbOfGaussPt;j++)
-              arrPtr[offsetPtr[*cellId]+j]=weights[j]*volPtr[*cellId];
+          {
+            std::vector<mcIdType> conn;
+            umesh->getNodeIdsOfCell(*cellId,conn);
+            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,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 < meshDim ; ++j)
+                {
+                  double res = 0.0;
+                  for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
+                    res += ptsInCell[spaceDim*k+i] * shapeFunc->getIJ(iGPt,meshDim*k+j);
+                  jacobian[ i ][ j ] = res;
+                }
+              arrPtr[offsetPtr[*cellId]+iGPt]=std::abs( jacobian.toJacobian() )*loc.getWeight(FromIdType<int>(iGPt));
+            }
+          }
         }
       else
         {