* 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;
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;
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;
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);
}
}
#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
+#include "InterpKernelDenseMatrix.hxx"
#include "InterpKernelRootsMultiDim.hxx"
#include "MEDCouplingUMesh.hxx"
#include "InterpolationHelper.txx"
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 )
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;
}
}
}
}
-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 )
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);
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;
+}
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;
#include "MEDCouplingFieldOverTime.hxx"
#include "MEDCouplingDefinitionTime.hxx"
#include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingCartesianAMRMesh.hxx"
#include "MEDCouplingAMRAttribute.hxx"
#include "MEDCouplingMatrix.hxx"
%feature("autodoc", "1");
%feature("docstring");
+%newobject MEDCoupling::MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement;
%newobject MEDCoupling::MEDCouplingField::buildMeasureField;
%newobject MEDCoupling::MEDCouplingField::getLocalizationOfDiscr;
%newobject MEDCoupling::MEDCouplingField::computeTupleIdsToSelectFromCellIds;
%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();"
}
}
};
+
+ 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();
+ }
+ }
+ };
}
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)):
#include "MEDCouplingMappedExtrudedMesh.hxx"
#include "MEDCoupling1GTUMesh.hxx"
#include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingMultiFields.hxx"
#include "MEDCouplingPartDefinition.hxx"
#include "MEDCouplingCartesianAMRMesh.hxx"
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;
