src/MEDCoupling/MEDCouplingFieldDiscretizationOnNodesFE.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (EDF R&D)
  20
  21 #include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
  22 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
  23 #include "InterpKernelDenseMatrix.hxx"
  24 #include "InterpKernelRootsMultiDim.hxx"
  25 #include "MEDCouplingUMesh.hxx"
  26 #include "InterpolationHelper.txx"
  27 #include "InterpKernelGaussCoords.hxx"
  28
  29 #include <sstream>
  30
  31 using namespace MEDCoupling;
  32
  33 const char MEDCouplingFieldDiscretizationOnNodesFE::REPR[]="FE";
  34
  35 std::string MEDCouplingFieldDiscretizationOnNodesFE::getStringRepr() const
  36 {
  37   return std::string(REPR);
  38 }
  39
  40 void MEDCouplingFieldDiscretizationOnNodesFE::reprQuickOverview(std::ostream& stream) const
  41 {
  42   stream << "NodeFE spatial discretization.";
  43 }
  44
  45 MCAuto<MEDCouplingFieldDiscretization> MEDCouplingFieldDiscretizationOnNodesFE::aggregate(std::vector<const MEDCouplingFieldDiscretization *>& fds) const
  46 {
  47   return EasyAggregate<MEDCouplingFieldDiscretizationOnNodesFE>(fds);
  48 }
  49
  50 bool MEDCouplingFieldDiscretizationOnNodesFE::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
  51 {
  52   if(!other)
  53     {
  54       reason="other spatial discretization is NULL, and this spatial discretization (Node FE) is defined.";
  55       return false;
  56     }
  57   const MEDCouplingFieldDiscretizationOnNodesFE *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationOnNodesFE *>(other);
  58   bool ret=otherC!=0;
  59   if(!ret)
  60     reason="Spatial discrtization of this is ON_NODES_FE, which is not the case of other.";
  61   return ret;
  62 }
  63
  64 /*!
  65  * This method is simply called by MEDCouplingFieldDiscretization::deepCopy. It performs the deep copy of \a this.
  66  *
  67  * \sa MEDCouplingFieldDiscretization::deepCopy.
  68  */
  69 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationOnNodesFE::clone() const
  70 {
  71   return new MEDCouplingFieldDiscretizationOnNodesFE;
  72 }
  73
  74 void MEDCouplingFieldDiscretizationOnNodesFE::checkCompatibilityWithNature(NatureOfField nat) const
  75 {
  76   if(nat!=IntensiveMaximum)
  77     throw INTERP_KERNEL::Exception("Invalid nature for NodeFE field : expected IntensiveMaximum !");
  78 }
  79
  80 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
  81 {
  82   if(!mesh)
  83     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField : mesh instance specified is NULL !");
  84   throw INTERP_KERNEL::Exception("getMeasureField on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !");
  85 }
  86
  87 class Functor
  88 {
  89 private:
  90   const MEDCouplingGaussLocalization* _gl;
  91   std::size_t _nb_pts_in_cell;
  92   const double *_pts_in_cell;
  93   const double *_point;
  94 public:
  95   Functor(const MEDCouplingGaussLocalization& gl, std::size_t nbPtsInCell, const double *ptsInCell, const double point[3]):_gl(&gl),_nb_pts_in_cell(nbPtsInCell),
  96   _pts_in_cell(ptsInCell),_point(point) { }
  97   std::vector<double> operator()(const std::vector<double>& x)
  98   {
  99     MEDCouplingGaussLocalization gl(_gl->getType(),_gl->getRefCoords(),x,{1.0});
 100     MCAuto<DataArrayDouble> shapeFunc = gl.getShapeFunctionValues();
 101     const double *shapeFuncPtr( shapeFunc->begin() );
 102     std::vector<double> ret(3,0);
 103     for(std::size_t iPt = 0; iPt < _nb_pts_in_cell ; ++iPt)
 104     {
 105       for(short iDim = 0 ; iDim < 3 ; ++iDim)
 106         ret[iDim] += shapeFuncPtr[iPt] * _pts_in_cell[3*iPt + iDim];
 107     }
 108     ret[0] -= _point[0]; ret[1] -= _point[1]; ret[2] -= _point[2];
 109     return ret;
 110   }
 111 };
 112
 113 bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector<double>& ptsInCell, const double locInReal[3], double locInRef[3])
 114 {
 115   constexpr double EPS_IN_OUT = 1e-12;
 116   std::size_t nbPtsInCell(ptsInCell.size()/3);
 117   bool ret(false);
 118   const double *refCoo(gl.getRefCoords().data());
 119   INTERP_KERNEL::NormalizedCellType ct(gl.getType());
 120   Functor func(gl,nbPtsInCell,ptsInCell.data(),locInReal);
 121
 122   auto myJacobian = [&gl,nbPtsInCell,ptsInCell](const std::vector<double>& x, const std::vector<double>&, INTERP_KERNEL::DenseMatrix& jacobian)
 123   {
 124     MEDCouplingGaussLocalization mygl(gl.getType(),gl.getRefCoords(),x,{1.0});
 125     MCAuto<DataArrayDouble> shapeFunc = mygl.getDerivativeOfShapeFunctionValues();
 126     for(std::size_t i = 0 ; i < 3 ; ++i)
 127       for(std::size_t j = 0 ; j < 3 ; ++j)
 128       {
 129         double res = 0.0;
 130         for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
 131           res += ptsInCell[k*3+i] * shapeFunc->getIJ(0,3*k+j);
 132         jacobian[ i ][ j ] = res;
 133       }
 134   };
 135
 136   // loop on refcoords as initialization point for Newton algo. vini is the initialization vector of Newton.
 137   for(std::size_t attemptId = 0 ; attemptId < nbPtsInCell ; ++attemptId)
 138   {
 139     std::vector<double> vini(refCoo + attemptId*3, refCoo + (attemptId+1)*3);
 140     try
 141     {
 142       bool check(true);
 143       //INTERP_KERNEL::SolveWithNewton(vini,check,func);
 144       INTERP_KERNEL::SolveWithNewtonWithJacobian(vini,check,func,myJacobian);
 145       ret = (check==false);//looks strange but OK regarding newt (SolveWithNewton) at page 387 of numerical recipes for semantic of check parameter
 146     }
 147     catch( INTERP_KERNEL::Exception& ex )
 148       { ret = false; }// Something get wrong during Newton process
 149     if(ret)
 150     {//Newton has converged. Now check if it converged to a point inside cell
 151       if( ! INTERP_KERNEL::GaussInfo::IsInOrOutForReference(ct,vini.data(),EPS_IN_OUT) )
 152       {// converged but locInReal has been detected outside of cell
 153         ret = false;
 154       }
 155     }
 156     if(ret)
 157     {
 158       locInRef[0] = vini[0]; locInRef[1] = vini[1]; locInRef[2] = vini[2];
 159       return ret;
 160     }
 161   }
 162   std::fill(locInRef,locInRef+3,std::numeric_limits<double>::max());
 163   return false;
 164 }
 165
 166 void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *ptsCoo, double *res) const
 167 {
 168   MCAuto<DataArrayDouble> res2( this->getValueOnMulti(arr,mesh,ptsCoo,1) );
 169   std::copy(res2->begin(),res2->end(),res);
 170 }
 171
 172 void MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts,
 173   std::function<void(const MEDCouplingGaussLocalization&, const std::vector<mcIdType>&)> customFunc)
 174 {
 175   const double *coordsOfMesh( umesh->getCoords()->begin() );
 176   MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh);
 177   BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) );
 178   for(mcIdType iPt = 0 ; iPt < nbOfPts ; ++iPt)
 179   {
 180     std::vector<mcIdType> elems;
 181     tree.getElementsAroundPoint(ptsCoo+3*iPt,elems);
 182     bool found(false);
 183     for(auto cellId = elems.cbegin() ; cellId != elems.cend() && !found ; ++cellId)
 184     {
 185       INTERP_KERNEL::NormalizedCellType gt( umesh->getTypeOfCell(*cellId) );
 186       std::vector<mcIdType> conn;
 187       umesh->getNodeIdsOfCell(*cellId,conn);
 188       MCAuto<DataArrayDouble> refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) );
 189       std::vector<double> refCooCpp(refCoo->begin(),refCoo->end());
 190       std::vector<double> gsCoo(ptsCoo + iPt*3,ptsCoo + (iPt+1)*3);
 191       MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.});
 192       std::vector<double> ptsInCell; ptsInCell.reserve(conn.size()*gl.getDimension());
 193       std::for_each( conn.cbegin(), conn.cend(), [coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*3,coordsOfMesh+(c+1)*3); } );
 194       std::vector<double> locInRef(3);
 195       if( IsInside3D(gl,ptsInCell,gsCoo.data(),locInRef.data()) )
 196       {
 197         gl.setGaussCoords(locInRef);
 198         customFunc(gl,conn);
 199         found = true;
 200       }
 201     }
 202     if(!found)
 203       THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << ptsCoo[0] << " Y=" << ptsCoo[1] << " Z=" << ptsCoo[2] << " !");
 204   }
 205 }
 206
 207 const MEDCouplingUMesh *MEDCouplingFieldDiscretizationOnNodesFE::checkConfig3D(const MEDCouplingMesh *mesh) const
 208 {
 209   const MEDCouplingUMesh *umesh( dynamic_cast<const MEDCouplingUMesh *>(mesh) );
 210   if( !umesh )
 211     THROW_IK_EXCEPTION("getValueOn : not implemented yet for type != MEDCouplingUMesh !");
 212   if(umesh->getSpaceDimension() != 3 || umesh->getMeshDimension() != 3)
 213     THROW_IK_EXCEPTION("getValueOn : implemented only for meshes with spacedim == 3 and meshdim == 3 !");
 214   umesh->checkConsistency();
 215   return umesh;
 216 }
 217
 218 DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const
 219 {
 220   if(!arr || !arr->isAllocated())
 221     throw INTERP_KERNEL::Exception("getValueOnMulti : input array is null or not allocated !");
 222   mcIdType nbOfRows=getNumberOfMeshPlaces(mesh);
 223   if(arr->getNumberOfTuples()!=nbOfRows)
 224   {
 225     THROW_IK_EXCEPTION( "getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !")
 226   }
 227   const MEDCouplingUMesh *umesh = checkConfig3D(mesh);
 228   std::size_t nbCompo( arr->getNumberOfComponents() );
 229   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 230   ret->alloc(nbOfTargetPoints,nbCompo);
 231   double *res( ret->getPointer() );
 232
 233   auto arrayFeeder = [arr,&res,nbCompo](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
 234   {
 235     MCAuto<DataArrayDouble> resVector( gl.getShapeFunctionValues() );
 236     {
 237       std::for_each(res,res+nbCompo,[](double& v) { v = 0.0; });
 238       for(std::size_t iComp = 0 ; iComp < nbCompo ; ++iComp)
 239         for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt)
 240         {
 241           {
 242             res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp);
 243           }
 244         }
 245       res += nbCompo;
 246     }
 247   };
 248
 249   GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfTargetPoints,arrayFeeder);
 250   return ret.retn();
 251 }
 252
 253 /*!
 254  * Returns for each \a nbOfPoints point in \a loc its coordinate in reference element.
 255  */
 256 MCAuto<DataArrayDouble> MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement(const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const
 257 {
 258   const MEDCouplingUMesh *umesh = checkConfig3D(mesh);
 259   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 260   ret->alloc(nbOfPoints,3);
 261   double *retPtr(ret->getPointer() );
 262
 263   auto arrayFeeder = [&retPtr](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
 264   {
 265     std::vector<double> resVector( gl.getGaussCoords() );
 266     {
 267       std::copy(resVector.begin(),resVector.end(),retPtr);
 268       retPtr += 3;
 269     }
 270   };
 271
 272   GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfPoints,arrayFeeder);
 273   return ret;
 274 }