Back to index

salome-med  6.5.0
OverlapInterpolationMatrix.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "OverlapInterpolationMatrix.hxx"
00021 #include "ParaMESH.hxx"
00022 #include "ParaFIELD.hxx"
00023 #include "ProcessorGroup.hxx"
00024 #include "TranslationRotationMatrix.hxx"
00025 #include "Interpolation.hxx"
00026 #include "Interpolation1D.txx"
00027 #include "Interpolation2DCurve.hxx"
00028 #include "Interpolation2D.txx"
00029 #include "Interpolation3DSurf.hxx"
00030 #include "Interpolation3D.txx"
00031 #include "Interpolation3D2D.txx"
00032 #include "Interpolation2D1D.txx"
00033 #include "MEDCouplingUMesh.hxx"
00034 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
00035 #include "InterpolationOptions.hxx"
00036 #include "NormalizedUnstructuredMesh.hxx"
00037 #include "ElementLocator.hxx"
00038 #include "InterpKernelAutoPtr.hxx"
00039 
00040 #include <algorithm>
00041 
00042 using namespace std;
00043 
00044 namespace ParaMEDMEM
00045 {
00046   OverlapInterpolationMatrix::OverlapInterpolationMatrix(ParaFIELD *source_field,
00047                                                          ParaFIELD *target_field,
00048                                                          const ProcessorGroup& group,
00049                                                          const DECOptions& dec_options,
00050                                                          const INTERP_KERNEL::InterpolationOptions& i_opt):
00051     INTERP_KERNEL::InterpolationOptions(i_opt),
00052     DECOptions(dec_options),
00053     _source_field(source_field),
00054     _target_field(target_field),
00055     _source_support(source_field->getSupport()->getCellMesh()),
00056     _target_support(target_field->getSupport()->getCellMesh()),
00057     _mapping(group),
00058     _group(group)
00059   {
00060     int nbelems = source_field->getField()->getNumberOfTuples();
00061     _row_offsets.resize(nbelems+1);
00062     _coeffs.resize(nbelems);
00063     _target_volume.resize(nbelems);
00064   }
00065 
00066   void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayInt *ids)
00067   {
00068     _mapping.keepTracksOfSourceIds(procId,ids);
00069   }
00070 
00071   void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayInt *ids)
00072   {
00073     _mapping.keepTracksOfTargetIds(procId,ids);
00074   }
00075 
00076   OverlapInterpolationMatrix::~OverlapInterpolationMatrix()
00077   {
00078   }
00079 
00080   void OverlapInterpolationMatrix::addContribution(const MEDCouplingPointSet *src, const DataArrayInt *srcIds, const std::string& srcMeth, int srcProcId,
00081                                                    const MEDCouplingPointSet *trg, const DataArrayInt *trgIds, const std::string& trgMeth, int trgProcId)
00082   {
00083     std::string interpMethod(srcMeth);
00084     interpMethod+=trgMeth;
00085     //creating the interpolator structure
00086     vector<map<int,double> > surfaces;
00087     int colSize=0;
00088     //computation of the intersection volumes between source and target elements
00089     const MEDCouplingUMesh *trgC=dynamic_cast<const MEDCouplingUMesh *>(trg);
00090     const MEDCouplingUMesh *srcC=dynamic_cast<const MEDCouplingUMesh *>(src);
00091     if ( src->getMeshDimension() == -1 )
00092       {
00093         if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==2)
00094           {
00095             MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trgC);
00096             INTERP_KERNEL::Interpolation2D interpolation(*this);
00097             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
00098           }
00099         else if(trgC->getMeshDimension()==3 && trgC->getSpaceDimension()==3)
00100           {
00101             MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(trgC);
00102             INTERP_KERNEL::Interpolation3D interpolation(*this);
00103             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
00104           }
00105         else if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==3)
00106           {
00107             MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(trgC);
00108             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00109             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
00110           }
00111         else
00112           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
00113       }
00114     else if ( trg->getMeshDimension() == -1 )
00115       {
00116         if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==2)
00117           {
00118             MEDCouplingNormalizedUnstructuredMesh<2,2> local_mesh_wrapper(srcC);
00119             INTERP_KERNEL::Interpolation2D interpolation(*this);
00120             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
00121           }
00122         else if(srcC->getMeshDimension()==3 && srcC->getSpaceDimension()==3)
00123           {
00124             MEDCouplingNormalizedUnstructuredMesh<3,3> local_mesh_wrapper(srcC);
00125             INTERP_KERNEL::Interpolation3D interpolation(*this);
00126             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
00127           }
00128         else if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==3)
00129           {
00130             MEDCouplingNormalizedUnstructuredMesh<3,2> local_mesh_wrapper(srcC);
00131             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00132             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
00133           }
00134         else
00135           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
00136       }
00137     else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 3
00138               && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
00139       {
00140         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
00141         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
00142         
00143         INTERP_KERNEL::Interpolation3D2D interpolator (*this);
00144         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
00145         target_wrapper.releaseTempArrays();
00146         source_wrapper.releaseTempArrays();
00147       }
00148     else if ( src->getMeshDimension() == 3 && trg->getMeshDimension() == 2
00149               && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
00150       {
00151         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
00152         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
00153         
00154         INTERP_KERNEL::Interpolation3D2D interpolator (*this);
00155         vector<map<int,double> > surfacesTranspose;
00156         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());//not a bug target in source.
00157         TransposeMatrix(surfacesTranspose,colSize,surfaces);
00158         colSize=surfacesTranspose.size();
00159         target_wrapper.releaseTempArrays();
00160         source_wrapper.releaseTempArrays();
00161       }
00162     else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2
00163               && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
00164       {
00165         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
00166         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
00167         
00168         INTERP_KERNEL::Interpolation2D1D interpolator (*this);
00169         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
00170         target_wrapper.releaseTempArrays();
00171         source_wrapper.releaseTempArrays();
00172       }
00173     else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 1
00174               && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
00175       {
00176         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
00177         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
00178         
00179         INTERP_KERNEL::Interpolation2D1D interpolator (*this);
00180         vector<map<int,double> > surfacesTranspose;
00181         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfacesTranspose,interpMethod.c_str());//not a bug target in source.
00182         TransposeMatrix(surfacesTranspose,colSize,surfaces);
00183         colSize=surfacesTranspose.size();
00184         target_wrapper.releaseTempArrays();
00185         source_wrapper.releaseTempArrays();
00186       }
00187     else if (trg->getMeshDimension() != _source_support->getMeshDimension())
00188       {
00189         throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
00190       }
00191     else if( src->getMeshDimension() == 1
00192              && src->getSpaceDimension() == 1 )
00193       {
00194         MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(trgC);
00195         MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(srcC);
00196 
00197         INTERP_KERNEL::Interpolation1D interpolation(*this);
00198         colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
00199         target_wrapper.releaseTempArrays();
00200         source_wrapper.releaseTempArrays();
00201       }
00202     else if( trg->getMeshDimension() == 1
00203              && trg->getSpaceDimension() == 2 )
00204       {
00205         MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(trgC);
00206         MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(srcC);
00207 
00208         INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
00209         colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
00210         target_wrapper.releaseTempArrays();
00211         source_wrapper.releaseTempArrays();
00212       }
00213     else if ( trg->getMeshDimension() == 2
00214               && trg->getSpaceDimension() == 3 )
00215       {
00216         MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(trgC);
00217         MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(srcC);
00218 
00219         INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
00220         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
00221         target_wrapper.releaseTempArrays();
00222         source_wrapper.releaseTempArrays();
00223       }
00224     else if ( trg->getMeshDimension() == 2
00225               && trg->getSpaceDimension() == 2)
00226       {
00227         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
00228         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
00229 
00230         INTERP_KERNEL::Interpolation2D interpolator (*this);
00231         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
00232         target_wrapper.releaseTempArrays();
00233         source_wrapper.releaseTempArrays();
00234       }
00235     else if ( trg->getMeshDimension() == 3
00236               && trg->getSpaceDimension() == 3 )
00237       {
00238         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
00239         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
00240 
00241         INTERP_KERNEL::Interpolation3D interpolator (*this);
00242         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
00243         target_wrapper.releaseTempArrays();
00244         source_wrapper.releaseTempArrays();
00245       }
00246     else
00247       {
00248         throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
00249       }
00250     bool needSourceSurf=isSurfaceComputationNeeded(srcMeth);
00251     MEDCouplingFieldDouble *source_triangle_surf=0;
00252     if(needSourceSurf)
00253       source_triangle_surf=src->getMeasureField(getMeasureAbsStatus());
00254     //
00255     fillDistributedMatrix(surfaces,srcIds,srcProcId,trgIds,trgProcId);
00256     //
00257     if(needSourceSurf)
00258       source_triangle_surf->decrRef();
00259   }
00260 
00264   void OverlapInterpolationMatrix::fillDistributedMatrix(const std::vector< std::map<int,double> >& res,
00265                                                          const DataArrayInt *srcIds, int srcProc,
00266                                                          const DataArrayInt *trgIds, int trgProc)
00267   {
00268     _mapping.addContributionST(res,srcIds,srcProc,trgIds,trgProc);
00269   }
00270 
00275   void OverlapInterpolationMatrix::prepare(const std::vector< std::vector<int> >& procsInInteraction)
00276   {
00277     if(_source_support)
00278       _mapping.prepare(procsInInteraction,_target_field->getField()->getNumberOfTuplesExpected());
00279     else
00280       _mapping.prepare(procsInInteraction,0);
00281   }
00282 
00283   void OverlapInterpolationMatrix::computeDeno()
00284   {
00285     if(_target_field->getField()->getNature()==ConservativeVolumic)
00286       _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected());
00287     else
00288       throw INTERP_KERNEL::Exception("Policy Not implemented yet : only ConservativeVolumic defined !");
00289   }
00290 
00291   void OverlapInterpolationMatrix::multiply()
00292   {
00293     _mapping.multiply(_source_field->getField(),_target_field->getField());
00294   }
00295 
00296   void OverlapInterpolationMatrix::transposeMultiply()
00297   {
00298     _mapping.transposeMultiply(_target_field->getField(),_source_field->getField());
00299   }
00300   
00301   bool OverlapInterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const
00302   {
00303     return method=="P0";
00304   }
00305 
00306   void OverlapInterpolationMatrix::TransposeMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
00307   {
00308     matOut.resize(nbColsMatIn);
00309     int id=0;
00310     for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
00311       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00312         matOut[(*iter2).first][id]=(*iter2).second;
00313   }
00314 }