From: Anthony Geay <anthony.geay@edf.fr>
Date: Wed, 24 Aug 2022 07:55:30 +0000 (+0200)
Subject: Jacobian is computed externally + test
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=995c03b6405921fdeedf8756557dddf4a13eaf38;p=tools%2Fmedcoupling.git

Jacobian is computed externally + test
---

diff --git a/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx
index 721ac00c9..68eb57acb 100644
--- a/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx
+++ b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx
@@ -150,7 +150,8 @@ namespace INTERP_KERNEL
    * check is true if the routine has converged to a local minimum.
    */
   template <class T>
-  void SolveWithNewton(std::vector<double> &x, bool &check, T &vecfunc)
+  void SolveWithNewtonWithJacobian(std::vector<double> &x, bool &check, T &vecfunc,
+    std::function<void(const std::vector<double>&, const std::vector<double>&, INTERP_KERNEL::DenseMatrix&)> jacobian)
   {
     const mcIdType MAXITS=200;
     const double TOLF=1.0e-8,TOLMIN=1.0e-12,STPMX=100.0;
@@ -160,7 +161,6 @@ namespace INTERP_KERNEL
     std::vector<double> g(n),p(n),xold(n);
     INTERP_KERNEL::DenseMatrix fjac(n,n);
     FMin<T> fmin(vecfunc);
-    JacobianCalculator<T> fdjac(vecfunc);
     std::vector<double> &fvec=fmin.getVector();
     f=fmin(x);
     test=0.0;
@@ -176,7 +176,7 @@ namespace INTERP_KERNEL
     stpmax=STPMX*std::fmax(std::sqrt(sum),double(n));
     for (its=0;its<MAXITS;its++)
     {
-      fjac=fdjac(x,fvec);
+      jacobian(x,fvec,fjac);//fjac=fdjac(x,fvec);
       for (i=0;i<n;i++)
       {
         sum=0.0;
@@ -217,6 +217,21 @@ namespace INTERP_KERNEL
       if (test < TOLX)
           return;
     }
-    THROW_IK_EXCEPTION("MAXITS exceeded in newt");
+    THROW_IK_EXCEPTION("MAXITS exceeded in SolveWithNewtonWithJacobian");
+  }
+
+  /*!
+   * check is false on normal return.
+   * check is true if the routine has converged to a local minimum.
+   */
+  template <class T>
+  void SolveWithNewton(std::vector<double> &x, bool &check, T &vecfunc)
+  {
+    JacobianCalculator<T> fdjac(vecfunc);
+    auto myJacobian = [&fdjac,vecfunc](const std::vector<double>& x, const std::vector<double>& fvec, INTERP_KERNEL::DenseMatrix& fjac)
+    {
+      fjac = fdjac(x,fvec);
+    };
+    SolveWithNewtonWithJacobian(x,check,vecfunc,myJacobian);
   }
 }
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx
index 4e0138911..fe696ab81 100755
--- a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx
@@ -20,6 +20,7 @@
 
 #include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
+#include "InterpKernelDenseMatrix.hxx"
 #include "InterpKernelRootsMultiDim.hxx"
 #include "MEDCouplingUMesh.hxx"
 #include "InterpolationHelper.txx"
@@ -117,14 +118,30 @@ bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector<double
   const double *refCoo(gl.getRefCoords().data());
   INTERP_KERNEL::NormalizedCellType ct(gl.getType());
   Functor func(gl,nbPtsInCell,ptsInCell.data(),locInReal);
-  // loop on refcoords as initialization point for Newton algo.
+
+  auto myJacobian = [&gl,nbPtsInCell,ptsInCell](const std::vector<double>& x, const std::vector<double>&, INTERP_KERNEL::DenseMatrix& jacobian)
+  {
+    MEDCouplingGaussLocalization mygl(gl.getType(),gl.getRefCoords(),x,{1.0});
+    MCAuto<DataArrayDouble> shapeFunc = mygl.getDerivativeOfShapeFunctionValues();
+    for(std::size_t i = 0 ; i < 3 ; ++i)
+      for(std::size_t j = 0 ; j < 3 ; ++j)
+      {
+        double res = 0.0;
+        for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
+          res += ptsInCell[k*3+i] * shapeFunc->getIJ(0,3*k+j);
+        jacobian[ i ][ j ] = res;
+      }
+  };
+
+  // loop on refcoords as initialization point for Newton algo. vini is the initialization vector of Newton.
   for(std::size_t attemptId = 0 ; attemptId < nbPtsInCell ; ++attemptId)
   {
     std::vector<double> vini(refCoo + attemptId*3, refCoo + (attemptId+1)*3);
     try
     {
       bool check(true);
-      INTERP_KERNEL::SolveWithNewton(vini,check,func);
+      //INTERP_KERNEL::SolveWithNewton(vini,check,func);
+      INTERP_KERNEL::SolveWithNewtonWithJacobian(vini,check,func,myJacobian);
       ret = (check==false);//looks strange but OK regarding newt (SolveWithNewton) at page 387 of numerical recipes for semantic of check parameter
     }
     catch( INTERP_KERNEL::Exception& ex )
@@ -132,7 +149,7 @@ bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector<double
     if(ret)
     {//Newton has converged. Now check if it converged to a point inside cell
       if( ! INTERP_KERNEL::GaussInfo::IsInOrOutForReference(ct,vini.data(),EPS_IN_OUT) )
-      {// converged but detected as outside
+      {// converged but locInReal has been detected outside of cell
         ret = false;
       }
     }
@@ -187,7 +204,7 @@ void MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(const
   }
 }
 
-DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const
+const MEDCouplingUMesh *MEDCouplingFieldDiscretizationOnNodesFE::checkConfig3D(const MEDCouplingMesh *mesh) const
 {
   const MEDCouplingUMesh *umesh( dynamic_cast<const MEDCouplingUMesh *>(mesh) );
   if( !umesh )
@@ -195,14 +212,19 @@ DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const
   if(umesh->getSpaceDimension() != 3 || umesh->getMeshDimension() != 3)
     THROW_IK_EXCEPTION("getValueOn : implemented only for meshes with spacedim == 3 and meshdim == 3 !");
   umesh->checkConsistency();
+  return umesh;
+}
+
+DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const
+{
   if(!arr || !arr->isAllocated())
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti : input array is null or not allocated !");
+    throw INTERP_KERNEL::Exception("getValueOnMulti : input array is null or not allocated !");
   mcIdType nbOfRows=getNumberOfMeshPlaces(mesh);
   if(arr->getNumberOfTuples()!=nbOfRows)
   {
-    THROW_IK_EXCEPTION( "MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !")
+    THROW_IK_EXCEPTION( "getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !")
   }
-  umesh->checkConsistency();
+  const MEDCouplingUMesh *umesh = checkConfig3D(mesh);
   std::size_t nbCompo( arr->getNumberOfComponents() );
   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
   ret->alloc(nbOfTargetPoints,nbCompo);
@@ -227,3 +249,26 @@ DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const
   GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfTargetPoints,arrayFeeder);
   return ret.retn();
 }
+
+/*!
+ * Returns for each \a nbOfPoints point in \a loc its coordinate in reference element.
+ */
+MCAuto<DataArrayDouble> MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement(const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const
+{
+  const MEDCouplingUMesh *umesh = checkConfig3D(mesh);
+  MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+  ret->alloc(nbOfPoints,3);
+  double *retPtr(ret->getPointer() );
+  
+  auto arrayFeeder = [&retPtr](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
+  {
+    std::vector<double> resVector( gl.getGaussCoords() );
+    {
+      std::copy(resVector.begin(),resVector.end(),retPtr);
+      retPtr += 3;
+    }
+  };
+
+  GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfPoints,arrayFeeder);
+  return ret;
+}
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx
index 3aa28072b..378d0964c 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.hxx
@@ -44,9 +44,13 @@ namespace MEDCoupling
       MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const override;
       MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const override;
       MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const override;
+    public:
+      MEDCOUPLING_EXPORT MCAuto<DataArrayDouble> getCooInRefElement(const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const;
     public:
       MEDCOUPLING_EXPORT static void GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts,
   std::function<void(const MEDCouplingGaussLocalization&, const std::vector<mcIdType>&)> customFunc);
+    private:
+      const MEDCouplingUMesh *checkConfig3D(const MEDCouplingMesh *mesh) const;
     public:
       static const char REPR[];
       static constexpr TypeOfField TYPE = ON_NODES_FE;
diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i
index 1e25eb35a..7acf25f32 100644
--- a/src/MEDCoupling_Swig/MEDCouplingCommon.i
+++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i
@@ -47,6 +47,7 @@
 #include "MEDCouplingFieldOverTime.hxx"
 #include "MEDCouplingDefinitionTime.hxx"
 #include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
 #include "MEDCouplingCartesianAMRMesh.hxx"
 #include "MEDCouplingAMRAttribute.hxx"
 #include "MEDCouplingMatrix.hxx"
@@ -217,6 +218,7 @@ typedef long mcPyPtrType;
 %feature("autodoc", "1");
 %feature("docstring");
 
+%newobject MEDCoupling::MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement;
 %newobject MEDCoupling::MEDCouplingField::buildMeasureField;
 %newobject MEDCoupling::MEDCouplingField::getLocalizationOfDiscr;
 %newobject MEDCoupling::MEDCouplingField::computeTupleIdsToSelectFromCellIds;
@@ -506,6 +508,7 @@ typedef long mcPyPtrType;
 %feature("unref") MEDCouplingFieldDiscretizationGauss "$this->decrRef();"
 %feature("unref") MEDCouplingFieldDiscretizationGaussNE "$this->decrRef();"
 %feature("unref") MEDCouplingFieldDiscretizationKriging "$this->decrRef();"
+%feature("unref") MEDCouplingFieldDiscretizationOnNodesFE "$this->decrRef();"
 %feature("unref") MEDCouplingFieldDouble "$this->decrRef();"
 %feature("unref") MEDCouplingFieldFloat "$this->decrRef();"
 %feature("unref") MEDCouplingFieldInt32 "$this->decrRef();"
diff --git a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i
index 1ffd26910..deb4bee70 100644
--- a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i
+++ b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i
@@ -469,4 +469,21 @@ namespace MEDCoupling
       }
     }
   };
+
+  class MEDCouplingFieldDiscretizationOnNodesFE : public MEDCouplingFieldDiscretizationOnNodes
+  {
+    public:
+    %extend
+    {
+      DataArrayDouble *getCooInRefElement(const MEDCouplingMesh *mesh, PyObject *locs)
+      {
+        mcIdType sw,nbPts;
+        double v0; MEDCoupling::DataArrayDouble *v1(0); MEDCoupling::DataArrayDoubleTuple *v2(0); std::vector<double> v3;
+        const double *inp=convertObjToPossibleCpp5_Safe2(locs,sw,v0,v1,v2,v3,"wrap of MEDCouplingFieldDouble::getValueOnMulti",
+                                                         mesh->getSpaceDimension(),true,nbPts);
+        MCAuto<DataArrayDouble> ret(self->getCooInRefElement(mesh,inp,nbPts));
+        return ret.retn();
+      }
+    }
+  };
 }
diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
index cbb2c2f6d..ef5df95cd 100644
--- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
@@ -1586,6 +1586,62 @@ class MEDCouplingBasicsTest(unittest.TestCase):
             delta = abs(csr_new[0,0]-ref_value)/ref_value
             self.assertTrue(delta < 1e-3)
 
+    def testFEFE_0(self):
+        """
+        Test to stress localisation of target points into a source mesh using standard reference FE elements.
+        """
+        gt = NORM_HEXA27
+        nbPtsInCell = MEDCouplingUMesh.GetNumberOfNodesOfGeometricType( gt )
+        coo = DataArrayDouble([(9.0, 18.0, 27.0), (9.0, 22.0, 27.0), (11.0, 22.0, 27.0), (11.0, 18.0, 27.0), (9.0, 18.0, 33.0), (9.0, 22.0, 33.0), (11.0, 22.0, 33.0), (11.0, 18.0, 33.0), (8.8, 20.0, 26.4), (10.0, 21.6, 27.6), (11.2, 20.0, 26.4), (10.0, 18.4, 27.6), (8.8, 20.0, 33.6), (10.0, 21.6, 32.4), (11.2, 20.0, 33.6), (10.0, 18.4, 32.4), (8.8, 17.6, 30.0), (9.2, 21.6, 30.0), (11.2, 22.4, 30.0), (10.8, 18.4, 30.0), (10.0, 20.0, 26.4), (9.2, 20.0, 30.0), (10.0, 22.4, 30.0), (10.8, 20.0, 30.0), (10.0, 17.6, 30.0), (10.0, 20.0, 32.4), (10.0, 20.0, 30.0)])
+        m = MEDCouplingUMesh("mesh",3)
+        m.setCoords(coo)
+        m.allocateCells()
+        m.insertNextCell(gt,list(range(nbPtsInCell)))
+
+        inPts = DataArrayDouble( [(9.1, 18.2, 27.3), (9.1, 21.8, 27.3), (10.9, 21.8, 27.3), (10.9, 18.2, 27.3), (9.1, 18.2, 32.7), (9.1, 21.8, 32.7), (10.9, 21.8, 32.7), (10.9, 18.2, 32.7), (8.9, 20.0, 26.7), (10.0, 21.4, 27.9), (11.1, 20.0, 26.7), (10.0, 18.6, 27.9), (8.9, 20.0, 33.3), (10.0, 21.4, 32.1), (11.1, 20.0, 33.3), (10.0, 18.6, 32.1), (8.9, 17.8, 30.0), (9.3, 21.4, 30.0), (11.1, 22.2, 30.0), (10.7, 18.6, 30.0), (10.0, 20.0, 26.7), (9.3, 20.0, 30.0), (10.0, 22.2, 30.0), (10.7, 20.0, 30.0), (10.0, 17.8, 30.0), (10.0, 20.0, 32.1), (10.0, 20.0, 30.0)] )
+
+        outPts = DataArrayDouble( [(8.9, 17.8, 26.7), (8.9, 22.2, 26.7), (11.1, 22.2, 26.7), (11.1, 17.8, 26.7), (8.9, 17.8, 33.3), (8.9, 22.2, 33.3), (11.1, 22.2, 33.3), (11.1, 17.8, 33.3), (8.7, 20.0, 26.1), (10.0, 21.8, 27.3), (11.3, 20.0, 26.1), (10.0, 18.2, 27.3), (8.7, 20.0, 33.9), (10.0, 21.8, 32.7), (11.3, 20.0, 33.9), (10.0, 18.2, 32.7), (8.7, 17.4, 30.0), (9.1, 21.8, 30.0), (11.3, 22.6, 30.0), (10.9, 18.2, 30.0), (10.0, 20.0, 26.1), (9.1, 20.0, 30.0), (10.0, 22.6, 30.0), (10.9, 20.0, 30.0), (10.0, 17.4, 30.0), (10.0, 20.0, 32.7), (14.0, 20.0, 30.0)] )
+
+        srcField = MEDCouplingFieldDouble(ON_NODES_FE)
+        srcField.setMesh(m)
+        arr = DataArrayDouble(nbPtsInCell)
+        arr.iota()
+        srcField.setArray(arr)
+
+        ref_val0 = DataArrayDouble( [6.5782430384766535, 5.505760346014328, 7.80996256527073, 7.643290943690746, 9.758765408487054, 9.06408508454036, 11.003779543627997, 11.205026363340515, 10.56416007071563, 18.44359561721225, 12.499588353132655, 19.85351355074463, 14.186041573114885, 19.339743214851023, 16.084629460041207, 20.892007336402663, 17.269258227200577, 19.549962126638338, 19.039562190136063, 21.648928068870756, 20.667409503475078, 22.062499999998867, 22.562500000009678, 23.812499999505995, 24.395833333387696, 25.62468592706991, 26.] )
+
+        eps = 1e-8
+
+        for i in range(nbPtsInCell):
+            self.assertTrue( abs( srcField.getValueOn( inPts[i] )[0]-ref_val0[i] ) < eps )
+
+        self.assertTrue( srcField.getValueOnMulti( inPts ).isEqual(ref_val0,eps) )
+
+        srcFt=MEDCouplingFieldTemplate(srcField)
+        trgFt=MEDCouplingFieldTemplate(ON_NODES_FE)
+        trgMesh = MEDCouplingUMesh.Build0DMeshFromCoords( inPts )
+        trgFt.setMesh(trgMesh)
+        rem = MEDCouplingRemapper()
+        rem.setIntersectionType( PointLocator )
+        rem.prepareEx(srcFt,trgFt)
+        # scan content of matrix computed by remapper
+        mat = rem.getCrudeMatrix()
+        self.assertEqual( len(mat) , nbPtsInCell )
+        for irow,row in enumerate(mat):
+            self.assertTrue( abs( sum([y for x,y in row.items()]) - 1.0) < eps )
+            self.assertEqual( irow , [x for x,y in row.items() if y == max([y for x,y in row.items()])][0] ) # check that max coeff is for irow elt (due to construction of inPts )
+
+        # ask for MEDCouplingFieldDiscretizationOnNodesFE instance to compute coordination into ref element
+        sd = srcField.getDiscretization()
+        coosInRefElemFoundByNewton = sd.getCooInRefElement(srcField.getMesh(),inPts)
+
+        for zePt,cooInRefElemFoundByNewton in zip(inPts,coosInRefElemFoundByNewton):
+            # now check by performing refCoo -> realCoo
+            ftest = MEDCouplingFieldDouble(ON_GAUSS_PT)
+            ftest.setMesh(srcField.getMesh())
+            ftest.setGaussLocalizationOnType(gt,sum([list(elt) for elt in MEDCouplingGaussLocalization.GetDefaultReferenceCoordinatesOf(gt).getValuesAsTuple()],[]),list(cooInRefElemFoundByNewton),[1])
+            self.assertTrue ( float( (ftest.getLocalizationOfDiscr()-zePt).magnitude() ) < 1e-4 )
+
     def checkMatrix(self,mat1,mat2,nbCols,eps):
         self.assertEqual(len(mat1),len(mat2))
         for i in range(len(mat1)):
diff --git a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i
index 775ba4b48..e6029473e 100644
--- a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i
+++ b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i
@@ -28,6 +28,7 @@
 #include "MEDCouplingMappedExtrudedMesh.hxx"
 #include "MEDCoupling1GTUMesh.hxx"
 #include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
 #include "MEDCouplingMultiFields.hxx"
 #include "MEDCouplingPartDefinition.hxx"
 #include "MEDCouplingCartesianAMRMesh.hxx"
@@ -121,6 +122,8 @@ static PyObject *convertFieldDiscretization(MEDCoupling::MEDCouplingFieldDiscret
     ret=SWIG_NewPointerObj(reinterpret_cast<void*>(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationGaussNE,owner);
   if(dynamic_cast<MEDCoupling::MEDCouplingFieldDiscretizationKriging *>(fd))
     ret=SWIG_NewPointerObj(reinterpret_cast<void*>(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationKriging,owner);
+  if(dynamic_cast<MEDCoupling::MEDCouplingFieldDiscretizationOnNodesFE *>(fd))
+    ret=SWIG_NewPointerObj(reinterpret_cast<void*>(fd),SWIGTYPE_p_MEDCoupling__MEDCouplingFieldDiscretizationOnNodesFE,owner);
   if(!ret)
     throw INTERP_KERNEL::Exception("Not recognized type of field discretization on downcast !");
   return ret;