return false;
}
-void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
+void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *ptsCoo, double *res) const
{
- MCAuto<DataArrayDouble> res2( this->getValueOnMulti(arr,mesh,loc,1) );
+ MCAuto<DataArrayDouble> res2( this->getValueOnMulti(arr,mesh,ptsCoo,1) );
std::copy(res2->begin(),res2->end(),res);
}
-/*void GetValueOnMulti()
+void MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts,
+ std::function<void(const MEDCouplingGaussLocalization&, const std::vector<mcIdType>&)> customFunc)
{
-
-}*/
+ const double *coordsOfMesh( umesh->getCoords()->begin() );
+ MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh);
+ BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) );
+ for(mcIdType iPt = 0 ; iPt < nbOfPts ; ++iPt)
+ {
+ std::vector<mcIdType> elems;
+ tree.getElementsAroundPoint(ptsCoo+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(ptsCoo + iPt*3,ptsCoo + (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); } );
+ std::vector<double> locInRef(3);
+ if( IsInside3D(gl,ptsInCell,gsCoo.data(),locInRef.data()) )
+ {
+ gl.setGaussCoords(locInRef);
+ customFunc(gl,conn);
+ found = true;
+ }
+ }
+ if(!found)
+ THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << ptsCoo[0] << " Y=" << ptsCoo[1] << " Z=" << ptsCoo[2] << " !");
+ }
+}
DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const
{
std::size_t nbCompo( arr->getNumberOfComponents() );
MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
ret->alloc(nbOfTargetPoints,nbCompo);
- 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 )
+
+ auto arrayFeeder = [arr,&res,nbCompo](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
{
- std::vector<mcIdType> elems;
- tree.getElementsAroundPoint(loc+3*iPt,elems);
- bool found(false);
- for(auto cellId = elems.cbegin() ; cellId != elems.cend() && !found ; ++cellId)
+ MCAuto<DataArrayDouble> resVector( gl.getShapeFunctionValues() );
{
- 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 + 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); } );
- 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)
{
- 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);
- }
- }
+ {
+ res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp);
+ }
}
- found = true;
- }
+ res += nbCompo;
}
- if(!found)
- THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << loc[0] << " Y=" << loc[1] << " Z=" << loc[2] << " !");
- }
+ };
+
+ GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfTargetPoints,arrayFeeder);
return ret.retn();
}
#include "MEDCouplingCMesh.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
#include "MEDCouplingNormalizedCartesianMesh.txx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "Interpolation1D.txx"
#include "Interpolation2DCurve.hxx"
}
}
_matrix=m;
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(srcNbElem);
+ synchronizeSizeOfSideMatricesAfterMatrixComputation(srcNbElem);
}
int MEDCouplingRemapper::prepareInterpKernelOnly()
}
else
throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(nbCols);
- declareAsNew();
+ synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols);
return 1;
}
buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
target_mesh->getMesh3DIds()->getConstPointer());
//
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
- declareAsNew();
+ synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols2D*nbCols1D);
return 1;
}
ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
nullifiedTinyCoeffInCrudeMatrixAbs(0.);
//
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
- declareAsNew();
+ synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
return 1;
}
}
nullifiedTinyCoeffInCrudeMatrixAbs(0.);
//
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
- declareAsNew();
+ synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
return 1;
}
}
nullifiedTinyCoeffInCrudeMatrixAbs(0.);
//
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
- declareAsNew();
+ synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
return 1;
}
for(const mcIdType *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
_matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
}
- _deno_multiply.clear();
- _deno_multiply.resize(_matrix.size());
- _deno_reverse_multiply.clear();
- _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
- declareAsNew();
+ synchronizeSizeOfSideMatricesAfterMatrixComputation( srcLoc->getNumberOfTuples() );
return 1;
}
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();
+ const MEDCouplingUMesh *srcUMesh( dynamic_cast<const MEDCouplingUMesh *>(_src_ft->getMesh()) );
+ const MEDCouplingPointSet *trgMesh( dynamic_cast<const MEDCouplingPointSet *>(_target_ft->getMesh()) );
+ if( !srcUMesh )
+ THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D UMesh supported as source !");
+ if( !trgMesh )
+ THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D PointSet mesh supported as target !");
+
+ _matrix.clear();
+ _matrix.resize(trgMesh->getNumberOfNodes());
+ mcIdType rowId(0);
+
+ auto matrixFeeder = [this,&rowId](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
+ {
+ auto& row = this->_matrix[rowId++];
+ MCAuto<DataArrayDouble> resVector( gl.getShapeFunctionValues() );
+ for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt)
+ {
+ row[ conn[iPt] ] = resVector->getIJ(0,iPt);
+ }
+ };
+
+ MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(srcUMesh,trgMesh->getCoords()->begin(),trgMesh->getNumberOfNodes(),matrixFeeder);
+ synchronizeSizeOfSideMatricesAfterMatrixComputation( srcUMesh->getNumberOfNodes() );
return 1;
}
throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
}
+void MEDCouplingRemapper::synchronizeSizeOfSideMatricesAfterMatrixComputation(mcIdType nbOfColsInMatrix)
+{
+ _deno_multiply.clear();
+ _deno_multiply.resize(_matrix.size());
+ _deno_reverse_multiply.clear();
+ _deno_reverse_multiply.resize(nbOfColsInMatrix);
+ declareAsNew();
+}
+
/*!
* This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
* This method returns 3 information (2 in output parameters and 1 in return).
