X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FParaMEDMEM%2FInterpolationMatrix.cxx;h=1ea54be9f965b7dcf5ab6aa78f687fc12cca6500;hb=659f8c67d0348350e12fde38fe8c4de1ff95dffe;hp=e4605258cc76e05e5f41f7af9e9574de3f7a3da6;hpb=48782c06022ca2caa36f849cb5a29ea4fe2aaa83;p=tools%2Fmedcoupling.git diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index e4605258c..1ea54be9f 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -1,34 +1,43 @@ -// Copyright (C) 2007-2008 CEA/DEN, EDF R&D +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D // -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License. +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // + #include "ParaMESH.hxx" +#include "ParaFIELD.hxx" #include "ProcessorGroup.hxx" #include "MxN_Mapping.hxx" #include "InterpolationMatrix.hxx" #include "TranslationRotationMatrix.hxx" #include "Interpolation.hxx" +#include "Interpolation1D.txx" +#include "Interpolation2DCurve.hxx" #include "Interpolation2D.txx" -#include "Interpolation3DSurf.txx" +#include "Interpolation3DSurf.hxx" #include "Interpolation3D.txx" +#include "Interpolation3D2D.txx" +#include "Interpolation2D1D.txx" +#include "MEDCouplingUMesh.hxx" #include "MEDCouplingNormalizedUnstructuredMesh.txx" #include "InterpolationOptions.hxx" -#include "VolSurfFormulae.hxx" #include "NormalizedUnstructuredMesh.hxx" +#include "ElementLocator.hxx" + +#include // class InterpolationMatrix // This class enables the storage of an interpolation matrix Wij mapping @@ -51,28 +60,23 @@ namespace ParaMEDMEM // param method interpolation method // ==================================================================== - InterpolationMatrix::InterpolationMatrix(ParaMEDMEM::ParaMESH *source_support, + InterpolationMatrix::InterpolationMatrix(const ParaMEDMEM::ParaFIELD *source_field, const ProcessorGroup& source_group, const ProcessorGroup& target_group, const DECOptions& dec_options, const INTERP_KERNEL::InterpolationOptions& interp_options): - _source_support(source_support->getCellMesh()), + INTERP_KERNEL::InterpolationOptions(interp_options), + DECOptions(dec_options), + _source_field(source_field), + _source_support(source_field->getSupport()->getCellMesh()), _mapping(source_group, target_group, dec_options), _source_group(source_group), - _target_group(target_group), - _source_volume(0), - DECOptions(dec_options), - INTERP_KERNEL::InterpolationOptions(interp_options) + _target_group(target_group) { - int nbelems = _source_support->getNumberOfCells(); - + int nbelems = source_field->getField()->getNumberOfTuples(); _row_offsets.resize(nbelems+1); - for (int i=0; igetMeshDimension()) - { - throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions"); - } - std::string interpMethod(srcMeth); - interpMethod+=targetMeth; + std::string interpMethod(targetMeth); + interpMethod+=srcMeth; //creating the interpolator structure vector > surfaces; - int colSize=0; //computation of the intersection volumes between source and target elements + MEDCouplingUMesh *distant_supportC=dynamic_cast(&distant_support); + MEDCouplingUMesh *source_supportC=dynamic_cast(_source_support); + if ( distant_support.getMeshDimension() == -1 ) + { + if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==2) + { + MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(source_supportC); + INTERP_KERNEL::Interpolation2D interpolation(*this); + interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth); + } + else if(source_supportC->getMeshDimension()==3 && source_supportC->getSpaceDimension()==3) + { + MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(source_supportC); + INTERP_KERNEL::Interpolation3D interpolation(*this); + interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth); + } + else if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==3) + { + MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(source_supportC); + INTERP_KERNEL::Interpolation3DSurf interpolation(*this); + interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth); + } + else + throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh"); + } + else if ( source_supportC->getMeshDimension() == -1 ) + { + if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==2) + { + MEDCouplingNormalizedUnstructuredMesh<2,2> distant_mesh_wrapper(distant_supportC); + INTERP_KERNEL::Interpolation2D interpolation(*this); + interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth); + } + else if(distant_supportC->getMeshDimension()==3 && distant_supportC->getSpaceDimension()==3) + { + MEDCouplingNormalizedUnstructuredMesh<3,3> distant_mesh_wrapper(distant_supportC); + INTERP_KERNEL::Interpolation3D interpolation(*this); + interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth); + } + else if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==3) + { + MEDCouplingNormalizedUnstructuredMesh<3,2> distant_mesh_wrapper(distant_supportC); + INTERP_KERNEL::Interpolation3DSurf interpolation(*this); + interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth); + } + else + throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh"); + } + else if ( distant_support.getMeshDimension() == 2 + && _source_support->getMeshDimension() == 3 + && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3) + { + MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC); + INTERP_KERNEL::Interpolation3D2D interpolator (*this); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); + } + else if ( distant_support.getMeshDimension() == 1 + && _source_support->getMeshDimension() == 2 + && distant_support.getSpaceDimension() == 2 && _source_support->getSpaceDimension() == 2) + { + MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC); + INTERP_KERNEL::Interpolation2D1D interpolator (*this); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); + } + else if (distant_support.getMeshDimension() != _source_support->getMeshDimension()) + { + throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions"); + } + else if( distant_support.getMeshDimension() == 1 + && distant_support.getSpaceDimension() == 1 ) + { + MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(source_supportC); - if ( distant_support.getMeshDimension() == 2 - && distant_support.getSpaceDimension() == 3 ) + INTERP_KERNEL::Interpolation1D interpolation(*this); + interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); + } + else if( distant_support.getMeshDimension() == 1 + && distant_support.getSpaceDimension() == 2 ) { - MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(&distant_support); - MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(_source_support); + MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(source_supportC); + + INTERP_KERNEL::Interpolation2DCurve interpolation(*this); + interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); + } + else if ( distant_support.getMeshDimension() == 2 + && distant_support.getSpaceDimension() == 3 ) + { + MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC); INTERP_KERNEL::Interpolation3DSurf interpolator (*this); - colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); - target_wrapper.ReleaseTempArrays(); - source_wrapper.ReleaseTempArrays(); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); } else if ( distant_support.getMeshDimension() == 2 && distant_support.getSpaceDimension() == 2) { - MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(&distant_support); - MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(_source_support); + MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC); INTERP_KERNEL::Interpolation2D interpolator (*this); - colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); - target_wrapper.ReleaseTempArrays(); - source_wrapper.ReleaseTempArrays(); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); } else if ( distant_support.getMeshDimension() == 3 - && distant_support.getSpaceDimension() ==3 ) + && distant_support.getSpaceDimension() == 3 ) { - MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(&distant_support); - MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(_source_support); + MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC); INTERP_KERNEL::Interpolation3D interpolator (*this); - colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); - target_wrapper.ReleaseTempArrays(); - source_wrapper.ReleaseTempArrays(); + interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); } else { throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions "); } - - int source_size=surfaces.size(); + bool needTargetSurf=isSurfaceComputationNeeded(targetMeth); - MEDCouplingFieldDouble *target_triangle_surf = - getSupportVolumes(&distant_support); - MEDCouplingFieldDouble *source_triangle_surf = - getSupportVolumes(_source_support) ; - - // Storing the source volumes - _source_volume.resize(source_size); - for (int i=0; igetIJ(i,0); - } + MEDCouplingFieldDouble *target_triangle_surf=0; + if(needTargetSurf) + target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus()); + fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf); - source_triangle_surf->decrRef(); + if(needTargetSurf) + target_triangle_surf->decrRef(); + } + void InterpolationMatrix::fillDSFromVM(int iproc_distant, const int* distant_elems, const std::vector< std::map >& values, MEDCouplingFieldDouble *surf) + { //loop over the elements to build the interpolation //matrix structures - for (int ielem=0; ielem < surfaces.size(); ielem++) + int source_size=values.size(); + for (int ielem=0; ielem < source_size; ielem++) { - _row_offsets[ielem+1] += surfaces[ielem].size(); - //_source_indices.push_back(make_pair(iproc_distant, ielem)); - - for (map::const_iterator iter = surfaces[ielem].begin(); - iter != surfaces[ielem].end(); - iter++) + _row_offsets[ielem+1] += values[ielem].size(); + for(map::const_iterator iter=values[ielem].begin();iter!=values[ielem].end();iter++) { - //surface of the target triangle - double surf = target_triangle_surf->getIJ(iter->first,0); - + int localId; + if(distant_elems) + localId=distant_elems[iter->first]; + else + localId=iter->first; //locating the (iproc, itriangle) pair in the list of columns - vector >::iterator iter2 = - find(_col_offsets.begin(), _col_offsets.end(), - make_pair(iproc_distant,iter->first)); + map,int >::iterator iter2 = _col_offsets.find(make_pair(iproc_distant,localId)); int col_id; if (iter2 == _col_offsets.end()) { //(iproc, itriangle) is not registered in the list //of distant elements - - _col_offsets.push_back(make_pair(iproc_distant,iter->first)); col_id =_col_offsets.size(); - _mapping.addElementFromSource(iproc_distant, - distant_elems[iter->first]); - _target_volume.push_back(surf); + _col_offsets.insert(make_pair(make_pair(iproc_distant,localId),col_id)); + _mapping.addElementFromSource(iproc_distant,localId); } else { - col_id = iter2 - _col_offsets.begin() + 1; + col_id = iter2->second; } - //the non zero coefficient is stored //ielem is the row, //col_id is the number of the column //iter->second is the value of the coefficient - + if(surf) + { + double surface = surf->getIJ(iter->first,0); + _target_volume[ielem].push_back(surface); + } _coeffs[ielem].push_back(make_pair(col_id,iter->second)); } } - target_triangle_surf->decrRef(); + } + + void InterpolationMatrix::serializeMe(std::vector< std::vector< std::map > >& data1, std::vector& data2) const + { + data1.clear(); + data2.clear(); + const std::vector >& sendingIds=_mapping.getSendingIds(); + std::set procsS; + for(std::vector >::const_iterator iter1=sendingIds.begin();iter1!=sendingIds.end();iter1++) + procsS.insert((*iter1).first); + data1.resize(procsS.size()); + data2.resize(procsS.size()); + std::copy(procsS.begin(),procsS.end(),data2.begin()); + std::map fastProcAcc; + int id=0; + for(std::set::const_iterator iter2=procsS.begin();iter2!=procsS.end();iter2++,id++) + fastProcAcc[*iter2]=id; + int nbOfSrcElt=_coeffs.size(); + for(std::vector< std::vector< std::map > >::iterator iter3=data1.begin();iter3!=data1.end();iter3++) + (*iter3).resize(nbOfSrcElt); + id=0; + for(std::vector< std::vector< std::pair > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,id++) + { + for(std::vector< std::pair >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++) + { + const std::pair& elt=sendingIds[(*iter5).first]; + data1[fastProcAcc[elt.first]][id][elt.second]=(*iter5).second; + } + } + } + + void InterpolationMatrix::initialize() + { + int lgth=_coeffs.size(); + _row_offsets.clear(); _row_offsets.resize(lgth+1); + _coeffs.clear(); _coeffs.resize(lgth); + _target_volume.clear(); _target_volume.resize(lgth); + _col_offsets.clear(); + _mapping.initialize(); + } + + void InterpolationMatrix::finishContributionW(ElementLocator& elementLocator) + { + NatureOfField nature=elementLocator.getLocalNature(); + switch(nature) + { + case ConservativeVolumic: + computeConservVolDenoW(elementLocator); + break; + case Integral: + { + if(!elementLocator.isM1DCorr()) + computeIntegralDenoW(elementLocator); + else + computeGlobConstraintDenoW(elementLocator); + break; + } + case IntegralGlobConstraint: + computeGlobConstraintDenoW(elementLocator); + break; + case RevIntegral: + { + if(!elementLocator.isM1DCorr()) + computeRevIntegralDenoW(elementLocator); + else + computeConservVolDenoW(elementLocator); + break; + } + default: + throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field."); + break; + } + } + + void InterpolationMatrix::finishContributionL(ElementLocator& elementLocator) + { + NatureOfField nature=elementLocator.getLocalNature(); + switch(nature) + { + case ConservativeVolumic: + computeConservVolDenoL(elementLocator); + break; + case Integral: + { + if(!elementLocator.isM1DCorr()) + computeIntegralDenoL(elementLocator); + else + computeConservVolDenoL(elementLocator); + break; + } + case IntegralGlobConstraint: + //this is not a bug doing like ConservativeVolumic + computeConservVolDenoL(elementLocator); + break; + case RevIntegral: + { + if(!elementLocator.isM1DCorr()) + computeRevIntegralDenoL(elementLocator); + else + computeConservVolDenoL(elementLocator); + break; + } + default: + throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field."); + break; + } + } + + void InterpolationMatrix::computeConservVolDenoW(ElementLocator& elementLocator) + { + computeGlobalColSum(_deno_reverse_multiply); + computeGlobalRowSum(elementLocator,_deno_multiply,_deno_reverse_multiply); + } + + void InterpolationMatrix::computeConservVolDenoL(ElementLocator& elementLocator) + { + int pol1=elementLocator.sendPolicyToWorkingSideL(); + if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY) + { + elementLocator.recvFromWorkingSideL(); + elementLocator.sendToWorkingSideL(); + } + else if(ElementLocator::CUMULATIVE_POLICY) + { + //ask for lazy side to deduce ids eventually missing on working side and to send it back. + elementLocator.recvLocalIdsFromWorkingSideL(); + elementLocator.sendCandidatesGlobalIdsToWorkingSideL(); + elementLocator.recvCandidatesForAddElementsL(); + elementLocator.sendAddElementsToWorkingSideL(); + //Working side has updated its eventually missing ids updates its global ids with lazy side procs contribution + elementLocator.recvLocalIdsFromWorkingSideL(); + elementLocator.sendGlobalIdsToWorkingSideL(); + //like no post treatment + elementLocator.recvFromWorkingSideL(); + elementLocator.sendToWorkingSideL(); + } + else + throw INTERP_KERNEL::Exception("Not managed policy detected on lazy side : not implemented !"); + } + + void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator) + { + MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus()); + _deno_multiply.resize(_coeffs.size()); + vector >::iterator iter6=_deno_multiply.begin(); + const double *values=source_triangle_surf->getArray()->getConstPointer(); + for(vector > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++) + { + (*iter6).resize((*iter4).size()); + std::fill((*iter6).begin(),(*iter6).end(),*values); + } + source_triangle_surf->decrRef(); + _deno_reverse_multiply=_target_volume; + } + + void InterpolationMatrix::computeRevIntegralDenoW(ElementLocator& elementLocator) + { + _deno_multiply=_target_volume; + MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus()); + _deno_reverse_multiply.resize(_coeffs.size()); + vector >::iterator iter6=_deno_reverse_multiply.begin(); + const double *values=source_triangle_surf->getArray()->getConstPointer(); + for(vector > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++) + { + (*iter6).resize((*iter4).size()); + std::fill((*iter6).begin(),(*iter6).end(),*values); + } + source_triangle_surf->decrRef(); + } + + /*! + * Nothing to do because surface computation is on working side. + */ + void InterpolationMatrix::computeIntegralDenoL(ElementLocator& elementLocator) + { + } + + /*! + * Nothing to do because surface computation is on working side. + */ + void InterpolationMatrix::computeRevIntegralDenoL(ElementLocator& elementLocator) + { + } + + + void InterpolationMatrix::computeGlobConstraintDenoW(ElementLocator& elementLocator) + { + computeGlobalColSum(_deno_multiply); + computeGlobalRowSum(elementLocator,_deno_reverse_multiply,_deno_multiply); + } + + void InterpolationMatrix::computeGlobalRowSum(ElementLocator& elementLocator, std::vector >& denoStrorage, std::vector >& denoStrorageInv) + { + //stores id in distant procs sorted by lazy procs connected with + vector< vector > rowsPartialSumI; + //stores for each lazy procs connected with, if global info is available and if it's the case the policy + vector policyPartial; + //stores the corresponding values. + vector< vector > rowsPartialSumD; + elementLocator.recvPolicyFromLazySideW(policyPartial); + int pol1=mergePolicies(policyPartial); + if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY) + { + computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD); + elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD); + elementLocator.recvSumFromLazySideW(rowsPartialSumD); + } + else if(pol1==ElementLocator::CUMULATIVE_POLICY) + { + //updateWithNewAdditionnalElements(addingElements); + //stores for each lazy procs connected with, the ids in global mode if it exists (regarding policyPartial). This array has exactly the size of rowsPartialSumI, + //if policyPartial has CUMALATIVE_POLICY in each. + vector< vector > globalIdsPartial; + computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD); + elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI); + elementLocator.recvCandidatesGlobalIdsFromLazyProcsW(globalIdsPartial); + std::vector< std::vector > addingElements; + findAdditionnalElements(elementLocator,addingElements,rowsPartialSumI,globalIdsPartial); + addGhostElements(elementLocator.getDistantProcIds(),addingElements); + rowsPartialSumI.clear(); + globalIdsPartial.clear(); + computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD); + elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI); + elementLocator.recvGlobalIdsFromLazyProcsW(rowsPartialSumI,globalIdsPartial); + // + elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD); + elementLocator.recvSumFromLazySideW(rowsPartialSumD); + mergeRowSum3(globalIdsPartial,rowsPartialSumD); + mergeCoeffs(elementLocator.getDistantProcIds(),rowsPartialSumI,globalIdsPartial,denoStrorageInv); + } + else + throw INTERP_KERNEL::Exception("Not managed policy detected : not implemented !"); + divideByGlobalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD,denoStrorage); + } + + /*! + * @param distantProcs in parameter that indicates which lazy procs are concerned. + * @param resPerProcI out parameter that must be cleared before calling this method. The size of 1st dimension is equal to the size of 'distantProcs'. + * It contains the element ids (2nd dimension) of the corresponding lazy proc. + * @param resPerProcD out parameter with the same format than 'resPerProcI'. It contains corresponding sum values. + */ + void InterpolationMatrix::computeLocalRowSum(const std::vector& distantProcs, std::vector >& resPerProcI, + std::vector >& resPerProcD) const + { + resPerProcI.resize(distantProcs.size()); + resPerProcD.resize(distantProcs.size()); + std::vector res(_col_offsets.size()); + for(vector > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++) + for(vector >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++) + res[(*iter3).first]+=(*iter3).second; + set procsSet; + int id=-1; + const vector >& mapping=_mapping.getSendingIds(); + for(vector >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++) + { + std::pair::iterator,bool> isIns=procsSet.insert((*iter2).first); + if(isIns.second) + id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin(); + resPerProcI[id].push_back((*iter2).second); + resPerProcD[id].push_back(res[iter2-mapping.begin()]); + } + } + + /*! + * This method is only usable when CUMULATIVE_POLICY detected. This method finds elements ids (typically nodes) lazy side that + * are not present in columns of 'this' and that should regarding cumulative merge of elements regarding their global ids. + */ + void InterpolationMatrix::findAdditionnalElements(ElementLocator& elementLocator, std::vector >& elementsToAdd, + const std::vector >& resPerProcI, const std::vector >& globalIdsPartial) + { + std::set globalIds; + int nbLazyProcs=globalIdsPartial.size(); + for(int i=0;i tmp(globalIds.size()); + std::copy(globalIds.begin(),globalIds.end(),tmp.begin()); + globalIds.clear(); + elementLocator.sendCandidatesForAddElementsW(tmp); + elementLocator.recvAddElementsFromLazyProcsW(elementsToAdd); + } + + void InterpolationMatrix::addGhostElements(const std::vector& distantProcs, const std::vector >& elementsToAdd) + { + std::vector< std::vector< std::map > > data1; + std::vector data2; + serializeMe(data1,data2); + initialize(); + int nbOfDistProcs=distantProcs.size(); + for(int i=0;i& eltsForThisProc=elementsToAdd[i]; + if(!eltsForThisProc.empty()) + { + std::vector::iterator iter1=std::find(data2.begin(),data2.end(),procId); + std::map *toFeed=0; + if(iter1!=data2.end()) + {//to test + int rank=iter1-data2.begin(); + toFeed=&(data1[rank].back()); + } + else + { + iter1=std::lower_bound(data2.begin(),data2.end(),procId); + int rank=iter1-data2.begin(); + data2.insert(iter1,procId); + std::vector< std::map > tmp(data1.front().size()); + data1.insert(data1.begin()+rank,tmp); + toFeed=&(data1[rank].back()); + } + for(std::vector::const_iterator iter2=eltsForThisProc.begin();iter2!=eltsForThisProc.end();iter2++) + (*toFeed)[*iter2]=0.; + } + } + // + nbOfDistProcs=data2.size(); + for(int j=0;j& policyPartial) + { + if(policyPartial.empty()) + return ElementLocator::NO_POST_TREATMENT_POLICY; + int ref=policyPartial[0]; + std::vector::const_iterator iter1=std::find_if(policyPartial.begin(),policyPartial.end(),std::bind2nd(std::not_equal_to(),ref)); + if(iter1!=policyPartial.end()) + { + std::ostringstream msg; msg << "Incompatible policies between lazy procs each other : proc # " << iter1-policyPartial.begin(); + throw INTERP_KERNEL::Exception(msg.str().c_str()); + } + return ref; + } + + /*! + * This method introduce global ids aspects in computed 'rowsPartialSumD'. + * As precondition rowsPartialSumD.size()==policyPartial.size()==globalIdsPartial.size(). Foreach i in [0;rowsPartialSumD.size() ) rowsPartialSumD[i].size()==globalIdsPartial[i].size() + * @param rowsPartialSumD : in parameter, Partial row sum computed for each lazy procs connected with. + * @param rowsPartialSumI : in parameter, Corresponding local ids for each lazy procs connected with. + * @param globalIdsPartial : in parameter, the global numbering, of elements connected with. + * @param globalIdsLazySideInteraction : out parameter, constituted from all global ids of lazy procs connected with. + * @para sumCorresponding : out parameter, relative to 'globalIdsLazySideInteraction' + */ + void InterpolationMatrix::mergeRowSum(const std::vector< std::vector >& rowsPartialSumD, const std::vector< std::vector >& globalIdsPartial, + std::vector& globalIdsLazySideInteraction, std::vector& sumCorresponding) + { + std::map sumToReturn; + int nbLazyProcs=rowsPartialSumD.size(); + for(int i=0;i& rowSumOfP=rowsPartialSumD[i]; + const std::vector& globalIdsOfP=globalIdsPartial[i]; + std::vector::const_iterator iter1=rowSumOfP.begin(); + std::vector::const_iterator iter2=globalIdsOfP.begin(); + for(;iter1!=rowSumOfP.end();iter1++,iter2++) + sumToReturn[*iter2]+=*iter1; + } + // + int lgth=sumToReturn.size(); + globalIdsLazySideInteraction.resize(lgth); + sumCorresponding.resize(lgth); + std::vector::iterator iter3=globalIdsLazySideInteraction.begin(); + std::vector::iterator iter4=sumCorresponding.begin(); + for(std::map::const_iterator iter5=sumToReturn.begin();iter5!=sumToReturn.end();iter5++,iter3++,iter4++) + { + *iter3=(*iter5).first; + *iter4=(*iter5).second; + } + } + + /*! + * This method simply reorganize the result contained in 'sumCorresponding' computed by lazy side into 'rowsPartialSumD' with help of 'globalIdsPartial' and 'globalIdsLazySideInteraction' + * + * @param globalIdsPartial : in parameter, global ids sorted by lazy procs + * @param rowsPartialSumD : in/out parameter, with exactly the same size as 'globalIdsPartial' + * @param globalIdsLazySideInteraction : in parameter that represents ALL the global ids of every lazy procs in interaction + * @param sumCorresponding : in parameter with same size as 'globalIdsLazySideInteraction' that stores the corresponding sum of 'globalIdsLazySideInteraction' + */ + void InterpolationMatrix::mergeRowSum2(const std::vector< std::vector >& globalIdsPartial, std::vector< std::vector >& rowsPartialSumD, + const std::vector& globalIdsLazySideInteraction, const std::vector& sumCorresponding) + { + std::map acc; + std::vector::const_iterator iter1=globalIdsLazySideInteraction.begin(); + std::vector::const_iterator iter2=sumCorresponding.begin(); + for(;iter1!=globalIdsLazySideInteraction.end();iter1++,iter2++) + acc[*iter1]=*iter2; + // + int nbLazyProcs=globalIdsPartial.size(); + for(int i=0;i& tmp1=globalIdsPartial[i]; + std::vector& tmp2=rowsPartialSumD[i]; + std::vector::const_iterator iter3=tmp1.begin(); + std::vector::iterator iter4=tmp2.begin(); + for(;iter3!=tmp1.end();iter3++,iter4++) + *iter4=acc[*iter3]; + } } + void InterpolationMatrix::mergeRowSum3(const std::vector< std::vector >& globalIdsPartial, std::vector< std::vector >& rowsPartialSumD) + { + std::map sum; + std::vector< std::vector >::const_iterator iter1=globalIdsPartial.begin(); + std::vector< std::vector >::iterator iter2=rowsPartialSumD.begin(); + for(;iter1!=globalIdsPartial.end();iter1++,iter2++) + { + std::vector::const_iterator iter3=(*iter1).begin(); + std::vector::const_iterator iter4=(*iter2).begin(); + for(;iter3!=(*iter1).end();iter3++,iter4++) + sum[*iter3]+=*iter4; + } + iter2=rowsPartialSumD.begin(); + for(iter1=globalIdsPartial.begin();iter1!=globalIdsPartial.end();iter1++,iter2++) + { + std::vector::const_iterator iter3=(*iter1).begin(); + std::vector::iterator iter4=(*iter2).begin(); + for(;iter3!=(*iter1).end();iter3++,iter4++) + *iter4=sum[*iter3]; + } + } + + /*! + * This method updates this->_coeffs attribute in order to take into account hidden (because having the same global number) similar nodes in _coeffs array. + * If in this->_coeffs two distant element id have the same global id their values will be replaced for each by the sum of the two. + * @param procsInInteraction input parameter : specifies the procId in absolute of distant lazy procs in interaction with + * @param rowsPartialSumI input parameter : local ids of distant lazy procs elements in interaction with + * @param globalIdsPartial input parameter : global ids of distant lazy procs elements in interaction with + */ + void InterpolationMatrix::mergeCoeffs(const std::vector& procsInInteraction, const std::vector< std::vector >& rowsPartialSumI, + const std::vector >& globalIdsPartial, std::vector >& denoStrorageInv) + { + //preparing fast access structures + std::map procT; + int localProcId=0; + for(std::vector::const_iterator iter1=procsInInteraction.begin();iter1!=procsInInteraction.end();iter1++,localProcId++) + procT[*iter1]=localProcId; + int size=procsInInteraction.size(); + std::vector > localToGlobal(size); + for(int i=0;i& myLocalToGlobal=localToGlobal[i]; + const std::vector& locals=rowsPartialSumI[i]; + const std::vector& globals=globalIdsPartial[i]; + std::vector::const_iterator iter3=locals.begin(); + std::vector::const_iterator iter4=globals.begin(); + for(;iter3!=locals.end();iter3++,iter4++) + myLocalToGlobal[*iter3]=*iter4; + } + // + const vector >& mapping=_mapping.getSendingIds(); + std::map globalIdVal; + //accumulate for same global id on lazy part. + for(vector > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++) + for(vector >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++) + { + const std::pair& distantLocalLazyId=mapping[(*iter2).first]; + int localLazyProcId=procT[distantLocalLazyId.first]; + int globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second]; + globalIdVal[globalDistantLazyId]+=(*iter2).second; + } + //perform merge + std::vector >::iterator iter3=denoStrorageInv.begin(); + for(vector > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter3++) + { + double val=(*iter3).back(); + (*iter3).resize((*iter1).size()); + std::vector::iterator iter4=(*iter3).begin(); + for(vector >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter4++) + { + const std::pair& distantLocalLazyId=mapping[(*iter2).first]; + int localLazyProcId=procT[distantLocalLazyId.first]; + int globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second]; + double newVal=globalIdVal[globalDistantLazyId]; + if((*iter2).second!=0.) + (*iter4)=val*newVal/(*iter2).second; + else + (*iter4)=std::numeric_limits::max(); + (*iter2).second=newVal; + } + } + } + + void InterpolationMatrix::divideByGlobalRowSum(const std::vector& distantProcs, const std::vector >& resPerProcI, + const std::vector >& resPerProcD, std::vector >& deno) + { + map fastSums; + int procId=0; + for(vector::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++) + { + const std::vector& currentProcI=resPerProcI[procId]; + const std::vector& currentProcD=resPerProcD[procId]; + vector::const_iterator iter3=currentProcD.begin(); + for(vector::const_iterator iter2=currentProcI.begin();iter2!=currentProcI.end();iter2++,iter3++) + fastSums[_col_offsets[std::make_pair(*iter1,*iter2)]]=*iter3; + } + deno.resize(_coeffs.size()); + vector >::iterator iter6=deno.begin(); + for(vector > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++) + { + (*iter6).resize((*iter4).size()); + vector::iterator iter7=(*iter6).begin(); + for(vector >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++,iter7++) + *iter7=fastSums[(*iter5).first]; + } + } + + void InterpolationMatrix::computeGlobalColSum(std::vector >& denoStrorage) + { + denoStrorage.resize(_coeffs.size()); + vector >::iterator iter2=denoStrorage.begin(); + for(vector > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++) + { + (*iter2).resize((*iter1).size()); + double sumOfCurrentRow=0.; + for(vector >::const_iterator iter3=(*iter1).begin();iter3!=(*iter1).end();iter3++) + sumOfCurrentRow+=(*iter3).second; + std::fill((*iter2).begin(),(*iter2).end(),sumOfCurrentRow); + } + } + + void InterpolationMatrix::resizeGlobalColSum(std::vector >& denoStrorage) + { + vector >::iterator iter2=denoStrorage.begin(); + for(vector > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++) + { + double val=(*iter2).back(); + (*iter2).resize((*iter1).size()); + std::fill((*iter2).begin(),(*iter2).end(),val); + } + } // ================================================================== // The call to this method updates the arrays on the target side @@ -225,11 +839,11 @@ namespace ParaMEDMEM void InterpolationMatrix::prepare() { - int nbelems = _source_support->getNumberOfCells(); + int nbelems = _source_field->getField()->getNumberOfTuples(); for (int ielem=0; ielem < nbelems; ielem++) { _row_offsets[ielem+1]+=_row_offsets[ielem]; - } + } _mapping.prepareSendRecv(); } @@ -254,7 +868,7 @@ namespace ParaMEDMEM //computing the matrix multiply on source side if (_source_group.containsMyRank()) { - int nbrows = _source_support->getNumberOfCells(); + int nbrows = _coeffs.size(); // performing W.S // W is the intersection matrix @@ -269,23 +883,11 @@ namespace ParaMEDMEM { int colid= _coeffs[irow][icol-_row_offsets[irow]].first; double value = _coeffs[irow][icol-_row_offsets[irow]].second; - target_value[(colid-1)*nbcomp+icomp]+=value*coeff_row; + double deno = _deno_multiply[irow][icol-_row_offsets[irow]]; + target_value[colid*nbcomp+icomp]+=value*coeff_row/deno; } } } - - // performing VT^(-1).(W.S) - // where VT^(-1) is the inverse of the diagonal matrix containing - // the volumes of target cells - - for (int i=0; i<_col_offsets.size();i++) - { - for (int icomp=0; icompgetNumberOfComponents(); vector source_value(_col_offsets.size()* nbcomp,0.0); _mapping.reverseSendRecv(&source_value[0],field); @@ -328,7 +929,7 @@ namespace ParaMEDMEM //treatment of the transpose matrix multiply on the source side if (_source_group.containsMyRank()) { - int nbrows = _source_support->getNumberOfCells(); + int nbrows = _coeffs.size(); double *array = field.getArray()->getPointer() ; // Initialization @@ -343,232 +944,19 @@ namespace ParaMEDMEM { int colid = _coeffs[irow][icol-_row_offsets[irow]].first; double value = _coeffs[irow][icol-_row_offsets[irow]].second; - + double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]]; for (int icomp=0; icompisStructured()) - return getSupportUnstructuredVolumes((MEDCouplingUMesh *)mesh); - else - throw INTERP_KERNEL::Exception("Not implemented yet !!!"); - } - - // ==================================================================== - // brief returns the volumes of the cells underlying the field \a field - - // For 2D geometries, the returned field contains the areas. - // For 3D geometries, the returned field contains the volumes. - - // param field field on which cells the volumes are required - // return field containing the volumes - // ==================================================================== - - MEDCouplingFieldDouble* InterpolationMatrix::getSupportUnstructuredVolumes(MEDCouplingUMesh * mesh) - { - int ipt, type ; - int nbelem = mesh->getNumberOfCells() ; - int dim_mesh = mesh->getMeshDimension(); - int dim_space = mesh->getSpaceDimension() ; - double *coords = mesh->getCoords()->getPointer() ; - int *connec = mesh->getNodalConnectivity()->getPointer() ; - int *connec_index = mesh->getNodalConnectivityIndex()->getPointer() ; - - - MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS); - DataArrayDouble* array = DataArrayDouble::New() ; - array->alloc(nbelem, 1) ; - double *area_vol = array->getPointer() ; - - switch (dim_mesh) - { - case 2: // getting the areas - for ( int iel=0 ; ielsetArray(array) ; - array->decrRef(); - field->setMesh(mesh) ; - - return field ; + return method=="P0"; } }