WIP

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx
index 4494bec34082f7067567498047f3a0402812b4ac..db791fadb2ff1e67632330a6f6573d8d0f0ca51b 100755 (executable)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx
@@ -81,25 +81,31 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField
   throw INTERP_KERNEL::Exception("getMeasureField on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !");
 }
 
+bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector<double>& ptsInCell, const double locInReal[3], double locInRef[3])
+{
+  locInRef[0] = 0; locInRef[1] = 0; locInRef[2] = 0;
+  return true;
+}
+
 void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
+{
+  MCAuto<DataArrayDouble> res2( this->getValueOnMulti(arr,mesh,loc,1) );
+  std::copy(res2->begin(),res2->end(),res);
+}
+
+/*void GetValueOnMulti()
+{
+
+}*/
+
+DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const
 {
   const MEDCouplingUMesh *umesh( dynamic_cast<const MEDCouplingUMesh *>(mesh) );
   if( !umesh )
     THROW_IK_EXCEPTION("getValueOn : not implemented yet for type != MEDCouplingUMesh !");
   if(umesh->getSpaceDimension() != 3 || umesh->getMeshDimension() != 3)
     THROW_IK_EXCEPTION("getValueOn : implemented only for meshes with spacedim == 3 and meshdim == 3 !");
-  MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh);
-  BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) );
-  std::vector<mcIdType> elems;
-  tree.getElementsAroundPoint(loc,elems);
-  
-  /*MCAuto<DataArrayDouble> res2(this->getValueOnMulti(arr,mesh,loc,1));
-  std::copy(res2->begin(),res2->end(),res);*/
-  throw INTERP_KERNEL::Exception("getValueOn on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !");
-}
-
-DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const
-{
+  umesh->checkConsistency();
   if(!arr || !arr->isAllocated())
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti : input array is null or not allocated !");
   mcIdType nbOfRows=getNumberOfMeshPlaces(mesh);
@@ -107,12 +113,50 @@ DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const
   {
     THROW_IK_EXCEPTION( "MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !")
   }
-  throw INTERP_KERNEL::Exception("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !");
-  /*mcIdType nbCols(-1);
-  std::size_t nbCompo=arr->getNumberOfComponents();
-  MCAuto<DataArrayDouble> m(computeEvaluationMatrixOnGivenPts(mesh,loc,nbOfTargetPoints,nbCols));
+  umesh->checkConsistency();
+  std::size_t nbCompo( arr->getNumberOfComponents() );
   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
   ret->alloc(nbOfTargetPoints,nbCompo);
-  INTERP_KERNEL::matrixProduct(m->begin(),nbOfTargetPoints,nbCols,arr->begin(),nbOfRows,ToIdType(nbCompo),ret->getPointer());
-  return ret.retn();*/
+  const double *coordsOfMesh( umesh->getCoords()->begin() );
+  MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh);
+  BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) );
+  double *res( ret->getPointer() );
+  for(mcIdType iPt = 0 ; iPt < nbOfTargetPoints ; ++iPt, res+=nbCompo )
+  {
+    std::vector<mcIdType> elems;
+    tree.getElementsAroundPoint(loc+3*iPt,elems);
+    bool found(false);
+    for(auto cellId = elems.cbegin() ; cellId != elems.cend() && !found ; ++cellId)
+    {
+      INTERP_KERNEL::NormalizedCellType gt( umesh->getTypeOfCell(*cellId) );
+      std::vector<mcIdType> conn;
+      umesh->getNodeIdsOfCell(*cellId,conn);
+      MCAuto<DataArrayDouble> refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) );
+      std::vector<double> refCooCpp(refCoo->begin(),refCoo->end());
+      std::vector<double> gsCoo(loc,loc+3); 
+      MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.});
+      std::vector<double> ptsInCell; ptsInCell.reserve(conn.size()*gl.getDimension());
+      std::for_each( conn.cbegin(), conn.cend(), [coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*3,coordsOfMesh+(c+1)*3); } );
+      std::vector<double> locInRef(3);
+      if( IsInside3D(gl,ptsInCell,gsCoo.data(),locInRef.data()) )
+      {
+        gl.setGaussCoords(locInRef);
+        MCAuto<DataArrayDouble> resVector( gl.getShapeFunctionValues() );
+        {
+          std::for_each(res,res+nbCompo,[](double& v) { v = 0.0; });
+          for(std::size_t iComp = 0 ; iComp < nbCompo ; ++iComp)
+            for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt)
+            {
+              {
+                res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp);
+              }
+            }
+        }
+        found = true;
+      }
+    }
+    if(!found)
+      THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << loc[0] << " Y=" << loc[1] << " Z=" << loc[2] << " !");
+  }
+  return ret.retn();
 }

diff --git a/src/MEDCoupling/MEDCouplingMemArray.txx b/src/MEDCoupling/MEDCouplingMemArray.txx
index 9cb18e6fcce2424666f56e0248f59c2b129ba02f..16a07683f314eae43e8d051a0f90d985e9a0a33b 100755 (executable)
--- a/src/MEDCoupling/MEDCouplingMemArray.txx
+++ b/src/MEDCoupling/MEDCouplingMemArray.txx
@@ -5866,7 +5866,7 @@ struct NotInRange
       throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : only single component allowed !");
     std::size_t nbOfElements=this->getNumberOfTuples();
     if(nbOfElements<2)
-      throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : 1 tuple at least must be present in 'this' !");
+      throw INTERP_KERNEL::Exception("DataArrayInt::deltaShiftIndex : 2 tuples at least must be present in 'this' !");
     const T *ptr=this->getConstPointer();
     DataArrayType *ret=DataArrayType::New();
     ret->alloc(nbOfElements-1,1);

diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx
index 8f620ffe644cc7b2fc441c37e9c09ec2d66d6b19..c85f44ccceaea60876eb798983a895d85abf40c8 100755 (executable)
--- a/src/MEDCoupling/MEDCouplingRemapper.cxx
+++ b/src/MEDCoupling/MEDCouplingRemapper.cxx
@@ -181,6 +181,8 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnly()
   {
     case 0:
       return prepareNotInterpKernelOnlyGaussGauss();
+    case 1:
+      return prepareNotInterpKernelOnlyFEFE();
     default:
       {
         std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
@@ -957,6 +959,21 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
   return 1;
 }
 
+int MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE()
+{
+  if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
+    throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE : The intersection type is not supported ! Only PointLocator is supported for FE->FE interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !");
+  MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
+  mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents());
+  if(trgSpaceDim!=3)
+    THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for target !")
+  if(_src_ft->getMesh()->getSpaceDimension() != 3)
+    THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for source !")
+  
+  declareAsNew();
+  return 1;
+}
+
 /*!
  * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
  * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
@@ -965,9 +982,11 @@ int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel
 {
   if(method=="GAUSSGAUSS")
     return 0;
+  if(method=="FEFE")
+    return 1;
   std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
   oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
-  oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
+  oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS FEFE !";
   throw INTERP_KERNEL::Exception(oss.str().c_str());
 }
 

diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx
index 6966f6575348cc9023b4683878d8f194f45fba82..2acdafea1b183c0d7cb47c5917013543da6fa26e 100644 (file)
--- a/src/MEDCoupling/MEDCouplingRemapper.hxx
+++ b/src/MEDCoupling/MEDCouplingRemapper.hxx
@@ -89,6 +89,7 @@ namespace MEDCoupling
     //
     int prepareNotInterpKernelOnly();
     int prepareNotInterpKernelOnlyGaussGauss();
+    int prepareNotInterpKernelOnlyFEFE();
     //
     static int CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method);
     //