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);
/*!
test=0.0;
for (i=0;i<n;i++)
{
- temp=std::abs(p[i])/std::max(std::abs(xold[i]),1.0);
+ temp=std::abs(p[i])/std::fmax(std::abs(xold[i]),1.0);
if (temp > test) test=temp;
}
alamin=TOLX/test;
}
alam2=alam;
f2 = f;
- alam=std::max(tmplam,0.1*alam);
+ alam=std::fmax(tmplam,0.1*alam);
}
}
template <class T>
std::vector<double>& getVector() { return fvec; }
};
+ /*!
+ * 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)
{
}
sum=0.0;
for (i=0;i<n;i++) sum += sqr(x[i]);
- stpmax=STPMX*std::max(std::sqrt(sum),double(n));
+ stpmax=STPMX*std::fmax(std::sqrt(sum),double(n));
for (its=0;its<MAXITS;its++)
{
fjac=fdjac(x,fvec);
if (check)
{
test=0.0;
- den=std::max(f,0.5*n);
+ den=std::fmax(f,0.5*double(n));
for (i=0;i<n;i++) {
- temp=std::abs(g[i])*std::max(std::abs(x[i]),1.0)/den;
+ temp=std::abs(g[i])*std::fmax(std::abs(x[i]),1.0)/den;
if (temp > test) test=temp;
}
check=(test < TOLMIN);
test=0.0;
for (i=0;i<n;i++)
{
- temp=(std::abs(x[i]-xold[i]))/std::max(std::abs(x[i]),1.0);
+ temp=(std::abs(x[i]-xold[i]))/std::fmax(std::abs(x[i]),1.0);
if (temp > test) test=temp;
}
if (test < TOLX)
#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
+#include "InterpKernelRootsMultiDim.hxx"
#include "MEDCouplingUMesh.hxx"
#include "InterpolationHelper.txx"
+#include "InterpKernelGaussCoords.hxx"
#include <sstream>
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<double> operator()(const std::vector<double>& x)
+ {
+ MEDCouplingGaussLocalization gl(_gl->getType(),_gl->getRefCoords(),x,{1.0});
+ MCAuto<DataArrayDouble> shapeFunc = gl.getShapeFunctionValues();
+ const double *shapeFuncPtr( shapeFunc->begin() );
+ std::vector<double> 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<double>& 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<double> 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<double>::max());
+ return false;
}
void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
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);
+ std::vector<double> gsCoo(loc + iPt*3,loc + (iPt+1)*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); } );
