X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FParaMEDMEM%2FOverlapInterpolationMatrix.cxx;h=00b64df0c7b5262ccae49e84d749516ef917e16f;hb=b832b15337be013a56e0976170e5e235b89fcb03;hp=77404f179eee1f4d2421d85f7aee7c0f1d3970c3;hpb=fb512e2b77325290aaa2b4c9fd8f22d5949b6369;p=tools%2Fmedcoupling.git diff --git a/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx b/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx index 77404f179..00b64df0c 100644 --- a/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx +++ b/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2014 CEA/DEN, EDF R&D +// Copyright (C) 2007-2024 CEA, EDF // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -29,9 +29,10 @@ #include "Interpolation2D.txx" #include "Interpolation3DSurf.hxx" #include "Interpolation3D.txx" -#include "Interpolation3D2D.txx" +#include "Interpolation2D3D.txx" #include "Interpolation2D1D.txx" #include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" #include "MEDCouplingNormalizedUnstructuredMesh.txx" #include "InterpolationOptions.hxx" #include "NormalizedUnstructuredMesh.hxx" @@ -42,34 +43,30 @@ using namespace std; -namespace ParaMEDMEM +namespace MEDCoupling { OverlapInterpolationMatrix::OverlapInterpolationMatrix(ParaFIELD *source_field, ParaFIELD *target_field, const ProcessorGroup& group, const DECOptions& dec_options, - const INTERP_KERNEL::InterpolationOptions& i_opt): + const INTERP_KERNEL::InterpolationOptions& i_opt, + const OverlapElementLocator & locator): INTERP_KERNEL::InterpolationOptions(i_opt), DECOptions(dec_options), _source_field(source_field), _target_field(target_field), _source_support(source_field->getSupport()->getCellMesh()), _target_support(target_field->getSupport()->getCellMesh()), - _mapping(group), - _group(group) + _mapping(group, locator) { - int nbelems = source_field->getField()->getNumberOfTuples(); - _row_offsets.resize(nbelems+1); - _coeffs.resize(nbelems); - _target_volume.resize(nbelems); } - void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayInt *ids) + void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayIdType *ids) { _mapping.keepTracksOfSourceIds(procId,ids); } - void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayInt *ids) + void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayIdType *ids) { _mapping.keepTracksOfTargetIds(procId,ids); } @@ -78,14 +75,23 @@ namespace ParaMEDMEM { } - void OverlapInterpolationMatrix::addContribution(const MEDCouplingPointSet *src, const DataArrayInt *srcIds, const std::string& srcMeth, int srcProcId, - const MEDCouplingPointSet *trg, const DataArrayInt *trgIds, const std::string& trgMeth, int trgProcId) + // TODO? Merge with MEDCouplingRemapper::prepareInterpKernelOnlyUU() ? + /**! + * Local run (on this proc) of the sequential interpolation algorithm. + * + * @param srcIds is null if the source mesh is on the local proc + * @param trgIds is null if the source mesh is on the local proc + * + * One of the 2 is necessarily null (the two can be null together) + */ + void OverlapInterpolationMatrix::computeLocalIntersection(const MEDCouplingPointSet *src, const DataArrayIdType *srcIds, const std::string& srcMeth, int srcProcId, + const MEDCouplingPointSet *trg, const DataArrayIdType *trgIds, const std::string& trgMeth, int trgProcId) { std::string interpMethod(srcMeth); interpMethod+=trgMeth; //creating the interpolator structure - vector > surfaces; - int colSize=0; + vector sparse_matrix_part; + mcIdType colSize=0; //computation of the intersection volumes between source and target elements const MEDCouplingUMesh *trgC=dynamic_cast(trg); const MEDCouplingUMesh *srcC=dynamic_cast(src); @@ -95,19 +101,19 @@ namespace ParaMEDMEM { MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trgC); INTERP_KERNEL::Interpolation2D interpolation(*this); - colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth); + colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth); } else if(trgC->getMeshDimension()==3 && trgC->getSpaceDimension()==3) { MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(trgC); INTERP_KERNEL::Interpolation3D interpolation(*this); - colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth); + colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth); } else if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==3) { MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(trgC); INTERP_KERNEL::Interpolation3DSurf interpolation(*this); - colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth); + colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth); } else throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh"); @@ -118,19 +124,19 @@ namespace ParaMEDMEM { MEDCouplingNormalizedUnstructuredMesh<2,2> local_mesh_wrapper(srcC); INTERP_KERNEL::Interpolation2D interpolation(*this); - colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth); + colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth); } else if(srcC->getMeshDimension()==3 && srcC->getSpaceDimension()==3) { MEDCouplingNormalizedUnstructuredMesh<3,3> local_mesh_wrapper(srcC); INTERP_KERNEL::Interpolation3D interpolation(*this); - colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth); + colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth); } else if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==3) { MEDCouplingNormalizedUnstructuredMesh<3,2> local_mesh_wrapper(srcC); INTERP_KERNEL::Interpolation3DSurf interpolation(*this); - colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth); + colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth); } else throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh"); @@ -141,10 +147,8 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC); MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC); - INTERP_KERNEL::Interpolation3D2D interpolator (*this); - colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + INTERP_KERNEL::Interpolation2D3D interpolator (*this); + colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod); } else if ( src->getMeshDimension() == 3 && trg->getMeshDimension() == 2 && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 ) @@ -152,13 +156,11 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC); MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC); - INTERP_KERNEL::Interpolation3D2D interpolator (*this); - vector > surfacesTranspose; - colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);//not a bug target in source. - TransposeMatrix(surfacesTranspose,colSize,surfaces); - colSize=surfacesTranspose.size(); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + INTERP_KERNEL::Interpolation2D3D interpolator (*this); + vector matrixTranspose; + colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,sparse_matrix_part,interpMethod);//not a bug target in source. + TransposeMatrix(matrixTranspose,colSize,sparse_matrix_part); + colSize=ToIdType(matrixTranspose.size()); } else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 ) @@ -167,9 +169,7 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC); INTERP_KERNEL::Interpolation2D1D interpolator (*this); - colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod); } else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 1 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 ) @@ -178,12 +178,10 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC); INTERP_KERNEL::Interpolation2D1D interpolator (*this); - vector > surfacesTranspose; - colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfacesTranspose,interpMethod);//not a bug target in source. - TransposeMatrix(surfacesTranspose,colSize,surfaces); - colSize=surfacesTranspose.size(); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + vector matrixTranspose; + colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,matrixTranspose,interpMethod);//not a bug target in source. + TransposeMatrix(matrixTranspose,colSize,sparse_matrix_part); + colSize=ToIdType(matrixTranspose.size()); } else if (trg->getMeshDimension() != _source_support->getMeshDimension()) { @@ -196,9 +194,7 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(srcC); INTERP_KERNEL::Interpolation1D interpolation(*this); - colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod); } else if( trg->getMeshDimension() == 1 && trg->getSpaceDimension() == 2 ) @@ -207,9 +203,7 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(srcC); INTERP_KERNEL::Interpolation2DCurve interpolation(*this); - colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod); } else if ( trg->getMeshDimension() == 2 && trg->getSpaceDimension() == 3 ) @@ -218,9 +212,7 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(srcC); INTERP_KERNEL::Interpolation3DSurf interpolator (*this); - colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod); } else if ( trg->getMeshDimension() == 2 && trg->getSpaceDimension() == 2) @@ -229,9 +221,7 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC); INTERP_KERNEL::Interpolation2D interpolator (*this); - colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod); } else if ( trg->getMeshDimension() == 3 && trg->getSpaceDimension() == 3 ) @@ -240,58 +230,59 @@ namespace ParaMEDMEM MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC); INTERP_KERNEL::Interpolation3D interpolator (*this); - colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod); - target_wrapper.releaseTempArrays(); - source_wrapper.releaseTempArrays(); + colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod); } else { - throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions "); + throw INTERP_KERNEL::Exception("No interpolator exists for these mesh and space dimensions!"); } - bool needSourceSurf=isSurfaceComputationNeeded(srcMeth); - MEDCouplingFieldDouble *source_triangle_surf=0; - if(needSourceSurf) - source_triangle_surf=src->getMeasureField(getMeasureAbsStatus()); - // - fillDistributedMatrix(surfaces,srcIds,srcProcId,trgIds,trgProcId); - // - if(needSourceSurf) - source_triangle_surf->decrRef(); + /* Fill distributed matrix: + In sparse_matrix_part rows refer to target, and columns (=first param of map in SparseDoubleVec) + refer to source. + */ + _mapping.addContributionST(sparse_matrix_part,srcIds,srcProcId,trgIds,trgProcId); } /*! - * \b res rows refers to target and column (first param of map) to source. + * 'procsToSendField' gives the list of procs field data has to be sent to. + * See OverlapElementLocator::computeBoundingBoxesAndTodoList() */ - void OverlapInterpolationMatrix::fillDistributedMatrix(const std::vector< std::map >& res, - const DataArrayInt *srcIds, int srcProc, - const DataArrayInt *trgIds, int trgProc) - { - _mapping.addContributionST(res,srcIds,srcProc,trgIds,trgProc); - } - - /*! - * 'procsInInteraction' gives the global view of interaction between procs. - * In 'procsInInteraction' for a proc with id i, is in interaction with procs listed in procsInInteraction[i] - */ - void OverlapInterpolationMatrix::prepare(const std::vector< std::vector >& procsInInteraction) + void OverlapInterpolationMatrix::prepare(const std::vector< int >& procsToSendField) { if(_source_support) - _mapping.prepare(procsInInteraction,_target_field->getField()->getNumberOfTuplesExpected()); + _mapping.prepare(procsToSendField,_target_field->getField()->getNumberOfTuplesExpected()); else - _mapping.prepare(procsInInteraction,0); + _mapping.prepare(procsToSendField,0); } - void OverlapInterpolationMatrix::computeDeno() + void OverlapInterpolationMatrix::computeSurfacesAndDeno() { - if(_target_field->getField()->getNature()==ConservativeVolumic) + if(_target_field->getField()->getNature()==IntensiveMaximum) _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected()); else - throw INTERP_KERNEL::Exception("Policy Not implemented yet : only ConservativeVolumic defined !"); + throw INTERP_KERNEL::Exception("OverlapDEC: Policy not set (did you call setNature()?) or not implemented yet: only IntensiveMaximum!"); +// { +// if(_target_field->getField()->getNature()==IntensiveConservation) +// { +// MCAuto f; +// int orient = getOrientation(); // From InterpolationOptions inheritance +// if(orient == 2) // absolute areas +// f = _target_support->getMeasureField(true); +// else +// if(orient == 0) // relative areas +// f = _target_support->getMeasureField(false); +// else +// throw INTERP_KERNEL::Exception("OverlapDEC: orientation policy not impl. yet!"); +// _mapping.computeDenoRevIntegral(*f->getArray()); +// } +// else +// throw INTERP_KERNEL::Exception("OverlapDEC: Policy not implemented yet: only IntensiveMaximum and IntensiveConservation defined!"); +// } } - void OverlapInterpolationMatrix::multiply() + void OverlapInterpolationMatrix::multiply(double default_val) { - _mapping.multiply(_source_field->getField(),_target_field->getField()); + _mapping.multiply(_source_field->getField(),_target_field->getField(), default_val); } void OverlapInterpolationMatrix::transposeMultiply() @@ -299,17 +290,13 @@ namespace ParaMEDMEM _mapping.transposeMultiply(_target_field->getField(),_source_field->getField()); } - bool OverlapInterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const - { - return method=="P0"; - } - - void OverlapInterpolationMatrix::TransposeMatrix(const std::vector >& matIn, int nbColsMatIn, std::vector >& matOut) + void OverlapInterpolationMatrix::TransposeMatrix(const std::vector& matIn, + mcIdType nbColsMatIn, std::vector& matOut) { matOut.resize(nbColsMatIn); int id=0; - for(std::vector >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++) - for(std::map::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++) + for(std::vector::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++) + for(SparseDoubleVec::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++) matOut[(*iter2).first][id]=(*iter2).second; } }