From 4729a552f3f1a9126f02960394cf105147623df3 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Tue, 23 Aug 2022 16:01:44 +0200 Subject: [PATCH] First sucessful use of MEDCouplingFieldDiscretizationOnNodesFE::getValueOn on hexa27 --- .../GaussPoints/InterpKernelGaussCoords.cxx | 30 +++++++++ .../GaussPoints/InterpKernelGaussCoords.hxx | 1 + .../InterpKernelRootsMultiDim.hxx | 16 +++-- .../LinearAlgebra/InterpKernelQRDecomp.cxx | 2 +- src/MEDCoupling/CMakeLists.txt | 1 + ...EDCouplingFieldDiscretizationOnNodesFE.cxx | 65 ++++++++++++++++++- 6 files changed, 105 insertions(+), 10 deletions(-) diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx index 36e61963d..61feaf166 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx @@ -506,6 +506,36 @@ std::vector GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellTy THROW_IK_EXCEPTION("Input type " << ct << "is not managed by GetDefaultReferenceCoordinatesOf") } +/*! + * Returns true if \a ptInRefCoo is in reference cell of type \a ct (regarding GetDefaultReferenceCoordinatesOf) + * \sa GetDefaultReferenceCoordinatesOf + */ +bool GaussInfo::IsInOrOutForReference(NormalizedCellType ct, const double *ptInRefCoo, double eps) +{ + switch(ct) + { + case INTERP_KERNEL::NORM_SEG2: + case INTERP_KERNEL::NORM_SEG3: + { + return std::fabs(ptInRefCoo[0]) < 1.0+eps; + } + case INTERP_KERNEL::NORM_QUAD4: + case INTERP_KERNEL::NORM_QUAD8: + case INTERP_KERNEL::NORM_QUAD9: + { + return std::find_if(ptInRefCoo,ptInRefCoo+2,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+2; + } + case INTERP_KERNEL::NORM_HEXA8: + case INTERP_KERNEL::NORM_HEXA20: + case INTERP_KERNEL::NORM_HEXA27: + { + return std::find_if(ptInRefCoo,ptInRefCoo+3,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+3; + } + default: + THROW_IK_EXCEPTION("IsInOrOutForReference : not implemented for this geo type !"); + } +} + typedef void (*MapToShapeFunction)(GaussInfo& obj); /*! diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx index 4497180e2..5afdfa7d2 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx @@ -64,6 +64,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT static std::vector NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector& inputArray); INTERPKERNEL_EXPORT static std::vector GetDefaultReferenceCoordinatesOf(NormalizedCellType ct); + INTERPKERNEL_EXPORT static bool IsInOrOutForReference(NormalizedCellType ct, const double *ptInRefCoo, double eps); public: static const double SEG2A_REF[2]; static const double SEG2B_REF[2]; diff --git a/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx index 3fba59fdb..721ac00c9 100644 --- a/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx +++ b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx @@ -55,7 +55,7 @@ namespace INTERP_KERNEL test=0.0; for (i=0;i test) test=temp; } alamin=TOLX/test; @@ -93,7 +93,7 @@ namespace INTERP_KERNEL } alam2=alam; f2 = f; - alam=std::max(tmplam,0.1*alam); + alam=std::fmax(tmplam,0.1*alam); } } template @@ -145,6 +145,10 @@ namespace INTERP_KERNEL std::vector& getVector() { return fvec; } }; + /*! + * check is false on normal return. + * check is true if the routine has converged to a local minimum. + */ template void SolveWithNewton(std::vector &x, bool &check, T &vecfunc) { @@ -169,7 +173,7 @@ namespace INTERP_KERNEL } sum=0.0; for (i=0;i test) test=temp; } check=(test < TOLMIN); @@ -207,7 +211,7 @@ namespace INTERP_KERNEL test=0.0; for (i=0;i test) test=temp; } if (test < TOLX) diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx index 78afd9829..e180a9b8a 100644 --- a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx +++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx @@ -39,7 +39,7 @@ QRDecomp::QRDecomp(const INTERP_KERNEL::DenseMatrix& a): n(a.nrows()), qt(n,n), double scale,sigma,sum,tau; for (k=0;k @@ -81,10 +83,67 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField throw INTERP_KERNEL::Exception("getMeasureField on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !"); } +class Functor +{ +private: + const MEDCouplingGaussLocalization* _gl; + std::size_t _nb_pts_in_cell; + const double *_pts_in_cell; + const double *_point; +public: + Functor(const MEDCouplingGaussLocalization& gl, std::size_t nbPtsInCell, const double *ptsInCell, const double point[3]):_gl(&gl),_nb_pts_in_cell(nbPtsInCell), + _pts_in_cell(ptsInCell),_point(point) { } + std::vector operator()(const std::vector& x) + { + MEDCouplingGaussLocalization gl(_gl->getType(),_gl->getRefCoords(),x,{1.0}); + MCAuto shapeFunc = gl.getShapeFunctionValues(); + const double *shapeFuncPtr( shapeFunc->begin() ); + std::vector ret(3,0); + for(std::size_t iPt = 0; iPt < _nb_pts_in_cell ; ++iPt) + { + for(short iDim = 0 ; iDim < 3 ; ++iDim) + ret[iDim] += shapeFuncPtr[iPt] * _pts_in_cell[3*iPt + iDim]; + } + ret[0] -= _point[0]; ret[1] -= _point[1]; ret[2] -= _point[2]; + return ret; + } +}; + bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector& ptsInCell, const double locInReal[3], double locInRef[3]) { - locInRef[0] = 0; locInRef[1] = 0; locInRef[2] = 0; - return true; + constexpr double EPS_IN_OUT = 1e-12; + std::size_t nbPtsInCell(ptsInCell.size()/3); + bool ret(false); + 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. + for(std::size_t attemptId = 0 ; attemptId < nbPtsInCell ; ++attemptId) + { + std::vector vini(refCoo + attemptId*3, refCoo + (attemptId+1)*3); + try + { + bool check(true); + INTERP_KERNEL::SolveWithNewton(vini,check,func); + 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 ) + { ret = false; }// Something get wrong during Newton process + 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 + ret = false; + } + } + if(ret) + { + locInRef[0] = vini[0]; locInRef[1] = vini[1]; locInRef[2] = vini[2]; + return ret; + } + } + std::fill(locInRef,locInRef+3,std::numeric_limits::max()); + return false; } void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const @@ -133,7 +192,7 @@ DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const umesh->getNodeIdsOfCell(*cellId,conn); MCAuto refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) ); std::vector refCooCpp(refCoo->begin(),refCoo->end()); - std::vector gsCoo(loc,loc+3); + std::vector gsCoo(loc + iPt*3,loc + (iPt+1)*3); MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.}); std::vector 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); } ); -- 2.39.2