From eb89ffa7833f70341c6daf6f5b9b848e10a57711 Mon Sep 17 00:00:00 2001 From: abn Date: Tue, 24 Nov 2015 16:09:40 +0100 Subject: [PATCH] OverlapDEC: bug fix. Bounding box adjustment was negative. Also: documentation + typedefs to make code easier to read. To come: more tests + new algo for job distribution. --- src/MEDCoupling/MEDCouplingField.cxx | 2 +- src/MEDCoupling/MEDCouplingMemArray.cxx | 5 + src/MEDCoupling/MEDCouplingMemArray.hxx | 1 + src/MEDCoupling/MEDCouplingRemapper.cxx | 50 ++- src/ParaMEDMEM/ElementLocator.cxx | 2 +- src/ParaMEDMEM/OverlapDEC.cxx | 14 +- src/ParaMEDMEM/OverlapDEC.hxx | 2 +- src/ParaMEDMEM/OverlapElementLocator.cxx | 160 ++++---- src/ParaMEDMEM/OverlapElementLocator.hxx | 51 ++- src/ParaMEDMEM/OverlapInterpolationMatrix.cxx | 123 +++--- src/ParaMEDMEM/OverlapInterpolationMatrix.hxx | 65 +-- src/ParaMEDMEM/OverlapMapping.cxx | 378 ++++++++++-------- src/ParaMEDMEM/OverlapMapping.hxx | 83 ++-- src/ParaMEDMEMTest/CMakeLists.txt | 2 +- src/ParaMEDMEMTest/ParaMEDMEMTest.hxx | 19 +- .../ParaMEDMEMTest_Gauthier1.cxx | 1 + .../ParaMEDMEMTest_OverlapDEC.cxx | 349 +++++++++++++--- 17 files changed, 806 insertions(+), 501 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingField.cxx b/src/MEDCoupling/MEDCouplingField.cxx index 463a44bbb..cb69a6de7 100644 --- a/src/MEDCoupling/MEDCouplingField.cxx +++ b/src/MEDCoupling/MEDCouplingField.cxx @@ -503,7 +503,7 @@ MEDCouplingField::MEDCouplingField(const MEDCouplingField& other, bool deepCopy) /*! * Returns a new MEDCouplingMesh constituted by some cells of the underlying mesh of \a - * this filed, and returns ids of entities (nodes, cells, Gauss points) lying on the + * this field, and returns ids of entities (nodes, cells, Gauss points) lying on the * specified cells. The cells to include to the result mesh are specified by an array of * cell ids. The new mesh shares the coordinates array with the underlying mesh. * \param [in] start - an array of cell ids to include to the result mesh. diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 751741283..998ac654e 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -1628,6 +1628,11 @@ DataArrayDouble *DataArrayDouble::selectByTupleId(const int *new2OldBg, const in return ret.retn(); } +DataArrayDouble *DataArrayDouble::selectByTupleId(const DataArrayInt & di) const +{ + return selectByTupleId(di.getConstPointer(), di.getConstPointer()+di.getNumberOfTuples()); +} + /*! * Returns a shorten and permuted copy of \a this array. The new DataArrayDouble is * of size \a new2OldEnd - \a new2OldBg and it's values are permuted as required by diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index 1fc39204c..201b7601e 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -253,6 +253,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayDouble *renumberR(const int *new2Old) const; MEDCOUPLING_EXPORT DataArrayDouble *renumberAndReduce(const int *old2New, int newNbOfTuple) const; MEDCOUPLING_EXPORT DataArrayDouble *selectByTupleId(const int *new2OldBg, const int *new2OldEnd) const; + MEDCOUPLING_EXPORT DataArrayDouble *selectByTupleId(const DataArrayInt & di) const; MEDCOUPLING_EXPORT DataArrayDouble *selectByTupleIdSafe(const int *new2OldBg, const int *new2OldEnd) const; MEDCOUPLING_EXPORT DataArrayDouble *selectByTupleId2(int bg, int end2, int step) const; MEDCOUPLING_EXPORT DataArray *selectByTupleRanges(const std::vector >& ranges) const; diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index 18741778f..8dad02890 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -107,29 +107,41 @@ int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const ME int MEDCouplingRemapper::prepareInterpKernelOnly() { int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType(); + // *** Remember: +// typedef enum +// { +// UNSTRUCTURED = 5, +// CARTESIAN = 7, +// EXTRUDED = 8, +// CURVE_LINEAR = 9, +// SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10, +// SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11, +// IMAGE_GRID = 12 +// } MEDCouplingMeshType; + switch(meshInterpType) { - case 90: - case 91: - case 165: - case 181: - case 170: - case 171: - case 186: - case 187: - case 85://Unstructured-Unstructured + case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED + case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED + case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED + case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED + case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED + case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED + case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED + case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED + case 85: // UNSTRUCTURED - UNSTRUCTURED return prepareInterpKernelOnlyUU(); - case 167: - case 183: - case 87://Unstructured-Cartesian + case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN + case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN + case 87: // UNSTRUCTURED - CARTESIAN return prepareInterpKernelOnlyUC(); - case 122: - case 123: - case 117://Cartesian-Unstructured + case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED + case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED + case 117: // CARTESIAN - UNSTRUCTURED return prepareInterpKernelOnlyCU(); - case 119://Cartesian-Cartesian + case 119: // CARTESIAN - CARTESIAN return prepareInterpKernelOnlyCC(); - case 136://Extruded-Extruded + case 136: // EXTRUDED - EXTRUDED return prepareInterpKernelOnlyEE(); default: throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !"); @@ -462,7 +474,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces(); if(!duplicateFaces.empty()) { - std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n"; + std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n"; for(std::map >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++) { oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : "; @@ -866,7 +878,7 @@ int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel /*! * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal - * to IK_ONLY_PREFERED = 0 ) , which method will be applied. If \c true is returned the INTERP_KERNEL only method should be applied to \c false the \b not + * to IK_ONLY_PREFERED, which method will be applied. If \c true is returned the INTERP_KERNEL only method should be applied to \c false the \b not * only INTERP_KERNEL method should be applied. */ bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index 74e27545a..1374387ee 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -209,9 +209,9 @@ namespace ParaMEDMEM double* local_bb = _domain_bounding_boxes+_union_group->myRank()*2*_local_cell_mesh_space_dim; double* distant_bb = _domain_bounding_boxes+irank*2*_local_cell_mesh_space_dim; + const double eps = 1e-12; for (int idim=0; idim < _local_cell_mesh_space_dim; idim++) { - const double eps = 1e-12; bool intersects = (distant_bb[idim*2]getField()) + throw INTERP_KERNEL::Exception("OverlapDEC::synchronize(): currently, having a void source field on a proc is not allowed!"); + if (!_target_field || !_target_field->getField()) + throw INTERP_KERNEL::Exception("OverlapDEC::synchronize(): currently, having a void target field on a proc is not allowed!"); + if (_target_field->getField()->getNumberOfComponents() != _source_field->getField()->getNumberOfComponents()) + throw INTERP_KERNEL::Exception("OverlapDEC::synchronize(): source and target field have different number of components!"); delete _interpolation_matrix; _interpolation_matrix=new OverlapInterpolationMatrix(_source_field,_target_field,*_group,*this,*this); - OverlapElementLocator locator(_source_field,_target_field,*_group); + OverlapElementLocator locator(_source_field,_target_field,*_group, getBoundingBoxAdjustmentAbs()); locator.copyOptions(*this); locator.exchangeMeshes(*_interpolation_matrix); std::vector< std::pair > jobs=locator.getToDoList(); @@ -290,7 +296,7 @@ namespace ParaMEDMEM const DataArrayInt *trgIds=locator.getTargetIds((*it).second); _interpolation_matrix->addContribution(src,srcIds,srcMeth,(*it).first,trg,trgIds,trgMeth,(*it).second); } - _interpolation_matrix->prepare(locator.getProcsInInteraction()); + _interpolation_matrix->prepare(locator.getProcsToSendFieldData()); _interpolation_matrix->computeDeno(); } diff --git a/src/ParaMEDMEM/OverlapDEC.hxx b/src/ParaMEDMEM/OverlapDEC.hxx index b7b9b8c2f..0bf57b2ea 100644 --- a/src/ParaMEDMEM/OverlapDEC.hxx +++ b/src/ParaMEDMEM/OverlapDEC.hxx @@ -43,7 +43,7 @@ namespace ParaMEDMEM void synchronize(); void attachSourceLocalField(ParaFIELD *field, bool ownPt=false); void attachTargetLocalField(ParaFIELD *field, bool ownPt=false); - ProcessorGroup *getGrp() { return _group; } + ProcessorGroup *getGroup() { return _group; } bool isInGroup() const; private: bool _own_group; diff --git a/src/ParaMEDMEM/OverlapElementLocator.cxx b/src/ParaMEDMEM/OverlapElementLocator.cxx index 51560e145..134fc4008 100644 --- a/src/ParaMEDMEM/OverlapElementLocator.cxx +++ b/src/ParaMEDMEM/OverlapElementLocator.cxx @@ -39,20 +39,24 @@ using namespace std; namespace ParaMEDMEM { - OverlapElementLocator::OverlapElementLocator(const ParaFIELD *sourceField, const ParaFIELD *targetField, const ProcessorGroup& group) + const int OverlapElementLocator::START_TAG_MESH_XCH = 1140; + + OverlapElementLocator::OverlapElementLocator(const ParaFIELD *sourceField, const ParaFIELD *targetField, + const ProcessorGroup& group, double epsAbs) : _local_source_field(sourceField), _local_target_field(targetField), _local_source_mesh(0), _local_target_mesh(0), _domain_bounding_boxes(0), - _group(group) + _group(group), + _epsAbs(epsAbs) { if(_local_source_field) _local_source_mesh=_local_source_field->getSupport()->getCellMesh(); if(_local_target_field) _local_target_mesh=_local_target_field->getSupport()->getCellMesh(); _comm=getCommunicator(); - computeBoundingBoxes(); + computeBoundingBoxesAndTodoList(); } OverlapElementLocator::~OverlapElementLocator() @@ -66,7 +70,7 @@ namespace ParaMEDMEM return group->getComm(); } - void OverlapElementLocator::computeBoundingBoxes() + void OverlapElementLocator::computeBoundingBoxesAndTodoList() { CommInterface comm_interface=_group.getCommInterface(); const MPIProcessorGroup* group=static_cast (&_group); @@ -113,19 +117,21 @@ namespace ParaMEDMEM for(int j=0;j<_group.size();j++) { if(intersectsBoundingBox(i,j)) + { _proc_pairs[i].push_back(j); +// std::cout << "(" << _group.myRank() << ") PROC pair: " << i << "," << j << std::endl; + } } - // OK now let's assigning as balanced as possible, job to each proc of group - std::vector< std::vector< std::pair > > pairsToBeDonePerProc(_group.size()); + std::vector< std::vector< ProcCouple > > pairsToBeDonePerProc(_group.size()); int i=0; for(std::vector< std::vector< int > >::const_iterator it1=_proc_pairs.begin();it1!=_proc_pairs.end();it1++,i++) for(std::vector< int >::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) { if(pairsToBeDonePerProc[i].size()<=pairsToBeDonePerProc[*it2].size())//it includes the fact that i==*it2 - pairsToBeDonePerProc[i].push_back(std::pair(i,*it2)); + pairsToBeDonePerProc[i].push_back(ProcCouple(i,*it2)); else - pairsToBeDonePerProc[*it2].push_back(std::pair(i,*it2)); + pairsToBeDonePerProc[*it2].push_back(ProcCouple(i,*it2)); } //Keeping todo list of current proc. _to_do_list contains a set of pair where at least _group.myRank() appears once. //This proc will be in charge to perform interpolation of any of element of '_to_do_list' @@ -135,20 +141,32 @@ namespace ParaMEDMEM int myProcId=_group.myRank(); _to_do_list=pairsToBeDonePerProc[myProcId]; - //Feeding now '_procs_to_send'. A same id can appears twice. The second parameter in pair means what to send true=source, false=target - _procs_to_send.clear(); +// std::stringstream scout; +// scout << "(" << myProcId << ") my TODO list is: "; +// for (std::vector< ProcCouple >::const_iterator itdbg=_to_do_list.begin(); itdbg!=_to_do_list.end(); itdbg++) +// scout << "(" << (*itdbg).first << "," << (*itdbg).second << ")"; +// std::cout << scout.str() << "\n"; + + // Feeding now '_procs_to_send*'. A same id can appears twice. The second parameter in pair means what + // to send true=source, false=target + _procs_to_send_mesh.clear(); + _procs_to_send_field.clear(); for(int i=_group.size()-1;i>=0;i--) - if(i!=myProcId) - { - const std::vector< std::pair >& anRemoteProcToDoList=pairsToBeDonePerProc[i]; - for(std::vector< std::pair >::const_iterator it=anRemoteProcToDoList.begin();it!=anRemoteProcToDoList.end();it++) - { - if((*it).first==myProcId) - _procs_to_send.push_back(std::pair(i,true)); - if((*it).second==myProcId) - _procs_to_send.push_back(std::pair(i,false)); - } - } + { + const std::vector< ProcCouple >& anRemoteProcToDoList=pairsToBeDonePerProc[i]; + for(std::vector< ProcCouple >::const_iterator it=anRemoteProcToDoList.begin();it!=anRemoteProcToDoList.end();it++) + { + if((*it).first==myProcId) + { + if(i!=myProcId) + _procs_to_send_mesh.push_back(Proc_SrcOrTgt(i,true)); + _procs_to_send_field.push_back((*it).second); + } + if((*it).second==myProcId) + if(i!=myProcId) + _procs_to_send_mesh.push_back(Proc_SrcOrTgt(i,false)); + } + } } /*! @@ -159,39 +177,40 @@ namespace ParaMEDMEM { int myProcId=_group.myRank(); //starting to receive every procs whose id is lower than myProcId. - std::vector< std::pair > toDoListForFetchRemaining; - for(std::vector< std::pair >::const_iterator it=_to_do_list.begin();it!=_to_do_list.end();it++) + std::vector< ProcCouple > toDoListForFetchRemaining; + for(std::vector< ProcCouple >::const_iterator it=_to_do_list.begin();it!=_to_do_list.end();it++) { - if((*it).first!=(*it).second) + int first = (*it).first, second = (*it).second; + if(first!=second) { - if((*it).first==myProcId) + if(first==myProcId) { - if((*it).second((*it).first,(*it).second)); + toDoListForFetchRemaining.push_back(ProcCouple(first,second)); } else {//(*it).second==myProcId - if((*it).first((*it).first,(*it).second)); + toDoListForFetchRemaining.push_back(ProcCouple(first,second)); } } } //sending source or target mesh to remote procs - for(std::vector< std::pair >::const_iterator it2=_procs_to_send.begin();it2!=_procs_to_send.end();it2++) + for(std::vector< Proc_SrcOrTgt >::const_iterator it2=_procs_to_send_mesh.begin();it2!=_procs_to_send_mesh.end();it2++) sendLocalMeshTo((*it2).first,(*it2).second,matrix); //fetching remaining meshes - for(std::vector< std::pair >::const_iterator it=toDoListForFetchRemaining.begin();it!=toDoListForFetchRemaining.end();it++) + for(std::vector< ProcCouple >::const_iterator it=toDoListForFetchRemaining.begin();it!=toDoListForFetchRemaining.end();it++) { if((*it).first!=(*it).second) { if((*it).first==myProcId) - receiveRemoteMesh((*it).second,false); + receiveRemoteMeshFrom((*it).second,false); else//(*it).second==myProcId - receiveRemoteMesh((*it).first,true); + receiveRemoteMeshFrom((*it).first,true); } } } @@ -211,8 +230,8 @@ namespace ParaMEDMEM int myProcId=_group.myRank(); if(myProcId==procId) return _local_source_mesh; - std::pair p(procId,true); - std::map, MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > >::const_iterator it=_remote_meshes.find(p); + Proc_SrcOrTgt p(procId,true); + std::map::const_iterator it=_remote_meshes.find(p); return (*it).second; } @@ -221,8 +240,8 @@ namespace ParaMEDMEM int myProcId=_group.myRank(); if(myProcId==procId) return 0; - std::pair p(procId,true); - std::map, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p); + Proc_SrcOrTgt p(procId,true); + std::map::const_iterator it=_remote_elems.find(p); return (*it).second; } @@ -231,8 +250,8 @@ namespace ParaMEDMEM int myProcId=_group.myRank(); if(myProcId==procId) return _local_target_mesh; - std::pair p(procId,false); - std::map, MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > >::const_iterator it=_remote_meshes.find(p); + Proc_SrcOrTgt p(procId,false); + std::map::const_iterator it=_remote_meshes.find(p); return (*it).second; } @@ -241,11 +260,12 @@ namespace ParaMEDMEM int myProcId=_group.myRank(); if(myProcId==procId) return 0; - std::pair p(procId,false); - std::map, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > >::const_iterator it=_remote_elems.find(p); + Proc_SrcOrTgt p(procId,false); + std::map::const_iterator it=_remote_elems.find(p); return (*it).second; } + bool OverlapElementLocator::intersectsBoundingBox(int isource, int itarget) const { const double *source_bb=_domain_bounding_boxes+isource*2*2*_local_space_dim; @@ -253,9 +273,8 @@ namespace ParaMEDMEM for (int idim=0; idim < _local_space_dim; idim++) { - const double eps = -1e-12;//tony to change - bool intersects = (target_bb[idim*2] elems=local_mesh->getCellsInBoundingBox(distant_bb,getBoundingBoxAdjustment()); - DataArrayInt *idsToSend; - MEDCouplingPointSet *send_mesh=static_cast(field->getField()->buildSubMeshData(elems->begin(),elems->end(),idsToSend)); + AutoDAInt elems=local_mesh->getCellsInBoundingBox(distant_bb,getBoundingBoxAdjustment()); + DataArrayInt *old2new_map; + MEDCouplingPointSet *send_mesh=static_cast(field->getField()->buildSubMeshData(elems->begin(),elems->end(),old2new_map)); if(sourceOrTarget) - matrix.keepTracksOfSourceIds(procId,idsToSend);//Case#1 in Step2 of main algorithm. + matrix.keepTracksOfSourceIds(procId,old2new_map); else - matrix.keepTracksOfTargetIds(procId,idsToSend);//Case#0 in Step2 of main algorithm. - sendMesh(procId,send_mesh,idsToSend); + matrix.keepTracksOfTargetIds(procId,old2new_map); + sendMesh(procId,send_mesh,old2new_map); send_mesh->decrRef(); - idsToSend->decrRef(); + old2new_map->decrRef(); } /*! * This method recieves source remote mesh on proc 'procId' if sourceOrTarget==True * This method recieves target remote mesh on proc 'procId' if sourceOrTarget==False */ - void OverlapElementLocator::receiveRemoteMesh(int procId, bool sourceOrTarget) + void OverlapElementLocator::receiveRemoteMeshFrom(int procId, bool sourceOrTarget) { - DataArrayInt *da=0; + DataArrayInt *old2new_map=0; MEDCouplingPointSet *m=0; - receiveMesh(procId,m,da); - std::pair p(procId,sourceOrTarget); + receiveMesh(procId,m,old2new_map); + Proc_SrcOrTgt p(procId,sourceOrTarget); _remote_meshes[p]=m; - _remote_elems[p]=da; + _remote_elems[p]=old2new_map; } void OverlapElementLocator::sendMesh(int procId, const MEDCouplingPointSet *mesh, const DataArrayInt *idsToSend) const { CommInterface comInterface=_group.getCommInterface(); + // First stage : exchanging sizes vector tinyInfoLocalD;//tinyInfoLocalD not used for the moment vector tinyInfoLocal; @@ -325,16 +345,16 @@ namespace ParaMEDMEM int lgth[2]; lgth[0]=tinyInfoLocal.size(); lgth[1]=idsToSend->getNbOfElems(); - comInterface.send(&lgth,2,MPI_INT,procId,1140,*_comm); - comInterface.send(&tinyInfoLocal[0],tinyInfoLocal.size(),MPI_INT,procId,1141,*comm); + comInterface.send(&lgth,2,MPI_INT,procId,START_TAG_MESH_XCH,*_comm); + comInterface.send(&tinyInfoLocal[0],tinyInfoLocal.size(),MPI_INT,procId,START_TAG_MESH_XCH+1,*comm); // DataArrayInt *v1Local=0; DataArrayDouble *v2Local=0; mesh->serialize(v1Local,v2Local); - comInterface.send(v1Local->getPointer(),v1Local->getNbOfElems(),MPI_INT,procId,1142,*comm); - comInterface.send(v2Local->getPointer(),v2Local->getNbOfElems(),MPI_DOUBLE,procId,1143,*comm); + comInterface.send(v1Local->getPointer(),v1Local->getNbOfElems(),MPI_INT,procId,START_TAG_MESH_XCH+2,*comm); + comInterface.send(v2Local->getPointer(),v2Local->getNbOfElems(),MPI_DOUBLE,procId,START_TAG_MESH_XCH+3,*comm); //finished for mesh, ids now - comInterface.send(const_cast(idsToSend->getConstPointer()),lgth[1],MPI_INT,procId,1144,*comm); + comInterface.send(const_cast(idsToSend->getConstPointer()),lgth[1],MPI_INT,procId,START_TAG_MESH_XCH+4,*comm); // v1Local->decrRef(); v2Local->decrRef(); @@ -346,19 +366,19 @@ namespace ParaMEDMEM MPI_Status status; const MPI_Comm *comm=getCommunicator(); CommInterface comInterface=_group.getCommInterface(); - comInterface.recv(lgth,2,MPI_INT,procId,1140,*_comm,&status); + comInterface.recv(lgth,2,MPI_INT,procId,START_TAG_MESH_XCH,*_comm,&status); std::vector tinyInfoDistant(lgth[0]); ids=DataArrayInt::New(); ids->alloc(lgth[1],1); - comInterface.recv(&tinyInfoDistant[0],lgth[0],MPI_INT,procId,1141,*comm,&status); + comInterface.recv(&tinyInfoDistant[0],lgth[0],MPI_INT,procId,START_TAG_MESH_XCH+1,*comm,&status); mesh=MEDCouplingPointSet::BuildInstanceFromMeshType((MEDCouplingMeshType)tinyInfoDistant[0]); std::vector unusedTinyDistantSts; vector tinyInfoDistantD(1);//tinyInfoDistantD not used for the moment DataArrayInt *v1Distant=DataArrayInt::New(); DataArrayDouble *v2Distant=DataArrayDouble::New(); mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts); - comInterface.recv(v1Distant->getPointer(),v1Distant->getNbOfElems(),MPI_INT,procId,1142,*comm,&status); - comInterface.recv(v2Distant->getPointer(),v2Distant->getNbOfElems(),MPI_DOUBLE,procId,1143,*comm,&status); + comInterface.recv(v1Distant->getPointer(),v1Distant->getNbOfElems(),MPI_INT,procId,START_TAG_MESH_XCH+2,*comm,&status); + comInterface.recv(v2Distant->getPointer(),v2Distant->getNbOfElems(),MPI_DOUBLE,procId,START_TAG_MESH_XCH+3,*comm,&status); mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts); //finished for mesh, ids now comInterface.recv(ids->getPointer(),lgth[1],MPI_INT,procId,1144,*comm,&status); diff --git a/src/ParaMEDMEM/OverlapElementLocator.hxx b/src/ParaMEDMEM/OverlapElementLocator.hxx index 6ce2677f2..566395f9a 100644 --- a/src/ParaMEDMEM/OverlapElementLocator.hxx +++ b/src/ParaMEDMEM/OverlapElementLocator.hxx @@ -38,15 +38,17 @@ namespace ParaMEDMEM class ProcessorGroup; class OverlapInterpolationMatrix; + typedef std::pair ProcCouple; // a couple of proc IDs, typically used to define a exchange betw 2 procs + class OverlapElementLocator : public INTERP_KERNEL::InterpolationOptions { public: - OverlapElementLocator(const ParaFIELD *sourceField, const ParaFIELD *targetField, const ProcessorGroup& group); + OverlapElementLocator(const ParaFIELD *sourceField, const ParaFIELD *targetField, const ProcessorGroup& group,double epsAbs); virtual ~OverlapElementLocator(); const MPI_Comm *getCommunicator() const; void exchangeMeshes(OverlapInterpolationMatrix& matrix); std::vector< std::pair > getToDoList() const { return _to_do_list; } - std::vector< std::vector< int > > getProcsInInteraction() const { return _proc_pairs; } + std::vector< int > getProcsToSendFieldData() const { return _procs_to_send_field; } std::string getSourceMethod() const; std::string getTargetMethod() const; const MEDCouplingPointSet *getSourceMesh(int procId) const; @@ -54,36 +56,49 @@ namespace ParaMEDMEM const MEDCouplingPointSet *getTargetMesh(int procId) const; const DataArrayInt *getTargetIds(int procId) const; private: - void computeBoundingBoxes(); + void computeBoundingBoxesAndTodoList(); bool intersectsBoundingBox(int i, int j) const; void sendLocalMeshTo(int procId, bool sourceOrTarget, OverlapInterpolationMatrix& matrix) const; - void receiveRemoteMesh(int procId, bool sourceOrTarget); + void receiveRemoteMeshFrom(int procId, bool sourceOrTarget); void sendMesh(int procId, const MEDCouplingPointSet *mesh, const DataArrayInt *idsToSend) const; void receiveMesh(int procId, MEDCouplingPointSet* &mesh, DataArrayInt *&ids) const; private: + typedef MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > AutoMCPointSet; + typedef MEDCouplingAutoRefCountObjectPtr< DataArrayInt > AutoDAInt; + typedef std::pair Proc_SrcOrTgt; // a key indicating a proc ID and whether the data is for source mesh/field or target mesh/field + + static const int START_TAG_MESH_XCH; + const ParaFIELD *_local_source_field; const ParaFIELD *_local_target_field; + int _local_space_dim; MEDCouplingPointSet *_local_source_mesh; MEDCouplingPointSet *_local_target_mesh; - std::vector _distant_cell_meshes; - std::vector _distant_face_meshes; - //! of size _group.size(). Contains for each source proc i, the ids of proc j the targets interact with. This vector is common for all procs in _group. + + /*! of size _group.size(). Contains for each source proc i, the ids of proc j the targets interact with. + This vector is common for all procs in _group. */ std::vector< std::vector< int > > _proc_pairs; - //! list of interpolations couple to be done - std::vector< std::pair > _to_do_list; - std::vector< std::pair > _procs_to_send; - std::map, MEDCouplingAutoRefCountObjectPtr< MEDCouplingPointSet > > _remote_meshes; - std::map, MEDCouplingAutoRefCountObjectPtr< DataArrayInt > > _remote_elems; + //! list of interpolation couples to be done by this proc only. This is a simple extraction of the member _pairsToBeDonePerProc + std::vector< ProcCouple > _to_do_list; + //! list of procs the local proc will have to send mesh data to: + std::vector< Proc_SrcOrTgt > _procs_to_send_mesh; + /*! list of procs the local proc will have to send field data to for the final matrix-vector computation: + * This can be different from _procs_to_send_mesh because interpolation matrix bits are computed on a potentially + * different proc than the target one. */ + std::vector< int > _procs_to_send_field; + //! Set of distant meshes + std::map< Proc_SrcOrTgt, AutoMCPointSet > _remote_meshes; + //! Set of cell ID mappings for the above distant meshes (because only part of the meshes are exchanged) + std::map< Proc_SrcOrTgt, AutoDAInt > _remote_elems; double* _domain_bounding_boxes; - const ProcessorGroup& _group; + //! bounding box absolute adjustment + double _epsAbs; + std::vector _distant_proc_ids; + + const ProcessorGroup& _group; const MPI_Comm *_comm; - //Attributes only used by lazy side - //std::vector _values_added; - //std::vector< std::vector > _ids_per_working_proc; - //std::vector< std::vector > _ids_per_working_proc3; - //std::vector< std::vector > _values_per_working_proc; }; } diff --git a/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx b/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx index b57541bb8..95b7e1b94 100644 --- a/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx +++ b/src/ParaMEDMEM/OverlapInterpolationMatrix.cxx @@ -58,10 +58,6 @@ namespace ParaMEDMEM _mapping(group), _group(group) { - 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) @@ -78,13 +74,20 @@ namespace ParaMEDMEM { } + // TODO? Merge with MEDCouplingRemapper::prepareInterpKernelOnlyUU() ? + /**! + * @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::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) { std::string interpMethod(srcMeth); interpMethod+=trgMeth; //creating the interpolator structure - vector > surfaces; + vector sparse_matrix_part; int colSize=0; //computation of the intersection volumes between source and target elements const MEDCouplingUMesh *trgC=dynamic_cast(trg); @@ -95,19 +98,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 +121,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"); @@ -142,9 +145,7 @@ namespace ParaMEDMEM 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(); + 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 ) @@ -153,12 +154,10 @@ namespace ParaMEDMEM 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(); + 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=matrixTranspose.size(); } else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 ) @@ -167,9 +166,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 +175,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=matrixTranspose.size(); } else if (trg->getMeshDimension() != _source_support->getMeshDimension()) { @@ -196,9 +191,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 +200,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 +209,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 +218,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,45 +227,30 @@ 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. - */ - 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] + * 'procsToSendField' gives the list of procs field data has to be sent to. + * See OverlapElementLocator::computeBoundingBoxesAndTodoList() */ - 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() @@ -299,17 +271,18 @@ namespace ParaMEDMEM _mapping.transposeMultiply(_target_field->getField(),_source_field->getField()); } - bool OverlapInterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const - { - return method=="P0"; - } +// 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, + int 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; } } diff --git a/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx b/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx index 2190e9add..2447079ee 100644 --- a/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx +++ b/src/ParaMEDMEM/OverlapInterpolationMatrix.hxx @@ -54,7 +54,7 @@ namespace ParaMEDMEM void 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); - void prepare(const std::vector< std::vector >& procsInInteraction); + void prepare(const std::vector< int > & procsToSendField); void computeDeno(); @@ -63,68 +63,19 @@ namespace ParaMEDMEM void transposeMultiply(); virtual ~OverlapInterpolationMatrix(); -#if 0 - void addContribution(MEDCouplingPointSet& distant_support, int iproc_distant, - const int* distant_elems, const std::string& srcMeth, const std::string& targetMeth); - void finishContributionW(ElementLocator& elementLocator); - void finishContributionL(ElementLocator& elementLocator); - void multiply(MEDCouplingFieldDouble& field) const; - void transposeMultiply(MEDCouplingFieldDouble& field)const; - void prepare(); - int getNbRows() const { return _row_offsets.size(); } - MPIAccessDEC* getAccessDEC() { return _mapping.getAccessDEC(); } private: - void computeConservVolDenoW(ElementLocator& elementLocator); - void computeIntegralDenoW(ElementLocator& elementLocator); - void computeRevIntegralDenoW(ElementLocator& elementLocator); - void computeGlobConstraintDenoW(ElementLocator& elementLocator); - void computeConservVolDenoL(ElementLocator& elementLocator); - void computeIntegralDenoL(ElementLocator& elementLocator); - void computeRevIntegralDenoL(ElementLocator& elementLocator); - - void computeLocalColSum(std::vector& res) const; - void computeLocalRowSum(const std::vector& distantProcs, std::vector >& resPerProcI, - std::vector >& resPerProcD) const; - void computeGlobalRowSum(ElementLocator& elementLocator, std::vector >& denoStrorage, std::vector >& denoStrorageInv); - void computeGlobalColSum(std::vector >& denoStrorage); - void resizeGlobalColSum(std::vector >& denoStrorage); - void fillDSFromVM(int iproc_distant, const int* distant_elems, const std::vector< std::map >& values, MEDCouplingFieldDouble *surf); - void serializeMe(std::vector< std::vector< std::map > >& data1, std::vector& data2) const; - void initialize(); - void findAdditionnalElements(ElementLocator& elementLocator, std::vector >& elementsToAdd, - const std::vector >& resPerProcI, const std::vector >& globalIdsPartial); - void addGhostElements(const std::vector& distantProcs, const std::vector >& elementsToAdd); - int mergePolicies(const std::vector& policyPartial); - void mergeRowSum(const std::vector< std::vector >& rowsPartialSumD, const std::vector< std::vector >& globalIdsPartial, - std::vector& globalIdsLazySideInteraction, std::vector& sumCorresponding); - void mergeRowSum2(const std::vector< std::vector >& globalIdsPartial, std::vector< std::vector >& rowsPartialSumD, - const std::vector& globalIdsLazySideInteraction, const std::vector& sumCorresponding); - void mergeRowSum3(const std::vector< std::vector >& globalIdsPartial, std::vector< std::vector >& rowsPartialSumD); - void mergeCoeffs(const std::vector& procsInInteraction, const std::vector< std::vector >& rowsPartialSumI, - const std::vector >& globalIdsPartial, std::vector >& denoStrorageInv); - void divideByGlobalRowSum(const std::vector& distantProcs, const std::vector >& resPerProcI, - const std::vector >& resPerProcD, std::vector >& deno); -#endif - private: - bool isSurfaceComputationNeeded(const std::string& method) const; - void fillDistributedMatrix(const std::vector< std::map >& res, - const DataArrayInt *srcIds, int srcProc, - const DataArrayInt *trgIds, int trgProc); - static void TransposeMatrix(const std::vector >& matIn, int nbColsMatIn, std::vector >& matOut); + + static void TransposeMatrix(const std::vector& matIn, int nbColsMatIn, + std::vector& matOut); +// bool isSurfaceComputationNeeded(const std::string& method) const; private: - ParaMEDMEM::ParaFIELD *_source_field; - ParaMEDMEM::ParaFIELD *_target_field; - std::vector _row_offsets; - std::map, int > _col_offsets; + ParaFIELD *_source_field; + ParaFIELD *_target_field; MEDCouplingPointSet *_source_support; MEDCouplingPointSet *_target_support; - OverlapMapping _mapping; + OverlapMapping _mapping; const ProcessorGroup& _group; - std::vector< std::vector > _target_volume; - std::vector > > _coeffs; - std::vector > _deno_multiply; - std::vector > _deno_reverse_multiply; }; } diff --git a/src/ParaMEDMEM/OverlapMapping.cxx b/src/ParaMEDMEM/OverlapMapping.cxx index abb7a1d1b..6871f1f35 100644 --- a/src/ParaMEDMEM/OverlapMapping.cxx +++ b/src/ParaMEDMEM/OverlapMapping.cxx @@ -36,46 +36,51 @@ OverlapMapping::OverlapMapping(const ProcessorGroup& group):_group(group) } /*! - * This method keeps tracks of source ids to know in step 6 of main algorithm, which tuple ids to send away. - * This method incarnates item#1 of step2 algorithm. + * Keeps the link between a given a proc holding source mesh data, and the corresponding cell IDs. */ void OverlapMapping::keepTracksOfSourceIds(int procId, DataArrayInt *ids) { ids->incrRef(); - _src_ids_st2.push_back(ids); - _src_proc_st2.push_back(procId); + _sent_src_ids_st2.push_back(ids); + _sent_src_proc_st2.push_back(procId); } /*! - * This method keeps tracks of target ids to know in step 6 of main algorithm. - * This method incarnates item#0 of step2 algorithm. - */ + * Same as keepTracksOfSourceIds() but for target mesh data. +*/ void OverlapMapping::keepTracksOfTargetIds(int procId, DataArrayInt *ids) { ids->incrRef(); - _trg_ids_st2.push_back(ids); - _trg_proc_st2.push_back(procId); + _sent_trg_ids_st2.push_back(ids); + _sent_trg_proc_st2.push_back(procId); } /*! - * This method stores from a matrix in format Target(rows)/Source(cols) for a source procId 'srcProcId' and for a target procId 'trgProcId'. - * All ids (source and target) are in format of local ids. + * This method stores in the local members the contribution coming from a matrix in format + * Target(rows)/Source(cols) for a source procId 'srcProcId' and for a target procId 'trgProcId'. + * All IDs received here (source and target) are in the format of local IDs. + * + * @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 OverlapMapping::addContributionST(const std::vector< std::map >& matrixST, const DataArrayInt *srcIds, int srcProcId, const DataArrayInt *trgIds, int trgProcId) +void OverlapMapping::addContributionST(const std::vector< SparseDoubleVec >& matrixST, const DataArrayInt *srcIds, int srcProcId, const DataArrayInt *trgIds, int trgProcId) { _matrixes_st.push_back(matrixST); _source_proc_id_st.push_back(srcProcId); _target_proc_id_st.push_back(trgProcId); - if(srcIds) + if(srcIds) // source mesh part is remote {//item#1 of step2 algorithm in proc m. Only to know in advanced nb of recv ids [ (0,1) computed on proc1 and Matrix-Vector on proc1 ] _nb_of_src_ids_proc_st2.push_back(srcIds->getNumberOfTuples()); _src_ids_proc_st2.push_back(srcProcId); } - else + else // source mesh part is local {//item#0 of step2 algorithm in proc k std::set s; - for(std::vector< std::map >::const_iterator it1=matrixST.begin();it1!=matrixST.end();it1++) - for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) + // For all source IDs (=col indices) in the sparse matrix: + for(std::vector< SparseDoubleVec >::const_iterator it1=matrixST.begin();it1!=matrixST.end();it1++) + for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) s.insert((*it2).first); _src_ids_zip_st2.resize(_src_ids_zip_st2.size()+1); _src_ids_zip_st2.back().insert(_src_ids_zip_st2.back().end(),s.begin(),s.end()); @@ -84,14 +89,17 @@ void OverlapMapping::addContributionST(const std::vector< std::map > } /*! - * '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]. + * This method is in charge to send matrices in AlltoAll mode. + * + * 'procsToSendField' gives the list of procs field data has to be sent to. + * See OverlapElementLocator::computeBoundingBoxesAndTodoList() * - * This method is in charge to send matrixes in AlltoAll mode. - * After the call of this method 'this' contains the matrixST for all source elements of the current proc + * After the call of this method, 'this' contains the matrixST for all source cells of the current proc */ -void OverlapMapping::prepare(const std::vector< std::vector >& procsInInteraction, int nbOfTrgElems) +void OverlapMapping::prepare(const std::vector< int >& procsToSendField, int nbOfTrgElems) { +// printMatrixesST(); + CommInterface commInterface=_group.getCommInterface(); const MPIProcessorGroup *group=static_cast(&_group); const MPI_Comm *comm=group->getComm(); @@ -101,12 +109,6 @@ void OverlapMapping::prepare(const std::vector< std::vector >& procsInInter INTERP_KERNEL::AutoPtr nbsend3=new int[grpSize]; std::fill(nbsend,nbsend+grpSize,0); int myProcId=_group.myRank(); - _proc_ids_to_recv_vector_st.clear(); - int curProc=0; - for(std::vector< std::vector >::const_iterator it1=procsInInteraction.begin();it1!=procsInInteraction.end();it1++,curProc++) - if(std::find((*it1).begin(),(*it1).end(),myProcId)!=(*it1).end()) - _proc_ids_to_recv_vector_st.push_back(curProc); - _proc_ids_to_send_vector_st=procsInInteraction[myProcId]; for(std::size_t i=0;i<_matrixes_st.size();i++) if(_source_proc_id_st[i]==myProcId) nbsend[_target_proc_id_st[i]]=_matrixes_st[i].size(); @@ -149,9 +151,46 @@ void OverlapMapping::prepare(const std::vector< std::vector >& procsInInter bigArrRecv2,bigArrDRecv2,nbrecv3,nbrecv4); //updating _src_ids_zip_st2 and _src_ids_zip_st2 with received matrix. updateZipSourceIdsForFuture(); - //finish to fill _the_matrix_st with already in place matrix in _matrixes_st + //finish to fill _the_matrix_st with already in place matrix in _matrixes_st (local computation) finishToFillFinalMatrixST(); - //printTheMatrix(); + // Prepare proc lists for future field data exchange (mutliply()): + fillProcToSendRcvForMultiply(procsToSendField); + // Make some space on local proc: + _matrixes_st.clear(); + +// std::stringstream scout; +// scout << "\n(" << myProcId << ") to send:"; +// for (std::vector< int >::const_iterator itdbg=_proc_ids_to_send_vector_st.begin(); itdbg!=_proc_ids_to_send_vector_st.end(); itdbg++) +// scout << "," << *itdbg; +// scout << "\n(" << myProcId << ") to recv:"; +// for (std::vector< int >::const_iterator itdbg=_proc_ids_to_recv_vector_st.begin(); itdbg!=_proc_ids_to_recv_vector_st.end(); itdbg++) +// scout << "," << *itdbg; +// std::cout << scout.str() << "\n"; +// +// printTheMatrix(); +} + +/** + * Fill the members _proc_ids_to_send_vector_st and _proc_ids_to_recv_vector_st + * indicating for each proc to which proc it should send its source field data + * and from which proc it should receive source field data. + * + * Should be called once finishToFillFinalMatrixST() has been executed, i.e. when the + * local matrices are completly ready. + */ +void OverlapMapping::fillProcToSendRcvForMultiply(const std::vector< int >& procsToSendField) +{ +// _proc_ids_to_send_vector_st.clear(); + _proc_ids_to_recv_vector_st.clear(); + // Receiving side is easy - just inspect non-void terms in the final matrix + int i=0; + std::vector< std::vector< SparseDoubleVec > >::const_iterator it; + std::vector< SparseDoubleVec >::const_iterator it2; + for(it=_the_matrix_st.begin(); it!=_the_matrix_st.end(); it++, i++) + _proc_ids_to_recv_vector_st.push_back(_the_matrix_st_source_proc_id[i]); + + // Sending side was computed in OverlapElementLocator::computeBoundingBoxesAndTodoList() + _proc_ids_to_send_vector_st = procsToSendField; } /*! @@ -169,11 +208,11 @@ void OverlapMapping::computeDenoGlobConstraint() for(std::size_t j=0;j& mToFill=_the_deno_st[i][j]; - const std::map& m=_the_matrix_st[i][j]; - for(std::map::const_iterator it=m.begin();it!=m.end();it++) + SparseDoubleVec& mToFill=_the_deno_st[i][j]; + const SparseDoubleVec& m=_the_matrix_st[i][j]; + for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++) sum+=(*it).second; - for(std::map::const_iterator it=m.begin();it!=m.end();it++) + for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++) mToFill[(*it).first]=sum; } } @@ -184,7 +223,6 @@ void OverlapMapping::computeDenoGlobConstraint() */ void OverlapMapping::computeDenoConservativeVolumic(int nbOfTuplesTrg) { - CommInterface commInterface=_group.getCommInterface(); int myProcId=_group.myRank(); // _the_deno_st.clear(); @@ -193,24 +231,24 @@ void OverlapMapping::computeDenoConservativeVolumic(int nbOfTuplesTrg) std::vector deno(nbOfTuplesTrg); for(std::size_t i=0;i >& mat=_the_matrix_st[i]; + const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i]; int curSrcId=_the_matrix_st_source_proc_id[i]; - std::vector::iterator isItem1=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),curSrcId); + std::vector::iterator isItem1=std::find(_sent_trg_proc_st2.begin(),_sent_trg_proc_st2.end(),curSrcId); int rowId=0; - if(isItem1==_trg_proc_st2.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids. + if(isItem1==_sent_trg_proc_st2.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids. { - for(std::vector< std::map >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) - for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) + for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) + for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) deno[rowId]+=(*it2).second; } else {//item0 of step2 main algo. More complicated. std::vector::iterator fnd=isItem1;//std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),curSrcId); - int locId=std::distance(_trg_proc_st2.begin(),fnd); - const DataArrayInt *trgIds=_trg_ids_st2[locId]; + int locId=std::distance(_sent_trg_proc_st2.begin(),fnd); + const DataArrayInt *trgIds=_sent_trg_ids_st2[locId]; const int *trgIds2=trgIds->getConstPointer(); - for(std::vector< std::map >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) - for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) + for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) + for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) deno[trgIds2[rowId]]+=(*it2).second; } } @@ -218,26 +256,26 @@ void OverlapMapping::computeDenoConservativeVolumic(int nbOfTuplesTrg) for(std::size_t i=0;i >& mat=_the_matrix_st[i]; + const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i]; int curSrcId=_the_matrix_st_source_proc_id[i]; - std::vector::iterator isItem1=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),curSrcId); - std::vector< std::map >& denoM=_the_deno_st[i]; + std::vector::iterator isItem1=std::find(_sent_trg_proc_st2.begin(),_sent_trg_proc_st2.end(),curSrcId); + std::vector< SparseDoubleVec >& denoM=_the_deno_st[i]; denoM.resize(mat.size()); - if(isItem1==_trg_proc_st2.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids. + if(isItem1==_sent_trg_proc_st2.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids. { int rowId=0; - for(std::vector< std::map >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) - for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) + for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) + for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) denoM[rowId][(*it2).first]=deno[rowId]; } else { std::vector::iterator fnd=isItem1; - int locId=std::distance(_trg_proc_st2.begin(),fnd); - const DataArrayInt *trgIds=_trg_ids_st2[locId]; + int locId=std::distance(_sent_trg_proc_st2.begin(),fnd); + const DataArrayInt *trgIds=_sent_trg_ids_st2[locId]; const int *trgIds2=trgIds->getConstPointer(); - for(std::vector< std::map >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) - for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) + for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++) + for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) denoM[rowId][(*it2).first]=deno[trgIds2[rowId]]; } } @@ -275,8 +313,8 @@ void OverlapMapping::serializeMatrixStep0ST(const int *nbOfElemsSrc, int *&bigAr int start=offsets[_target_proc_id_st[i]]; int *work=bigArr+start; *work=0; - const std::vector< std::map >& mat=_matrixes_st[i]; - for(std::vector< std::map >::const_iterator it=mat.begin();it!=mat.end();it++,work++) + const std::vector< SparseDoubleVec >& mat=_matrixes_st[i]; + for(std::vector< SparseDoubleVec >::const_iterator it=mat.begin();it!=mat.end();it++,work++) work[1]=work[0]+(*it).size(); } } @@ -295,6 +333,7 @@ void OverlapMapping::serializeMatrixStep0ST(const int *nbOfElemsSrc, int *&bigAr /*! * This method performs step#1 and step#2/3. It returns the size of expected array to get allToAllV. + * It is where the locally computed matrices are serialized to be sent to adequate final proc. */ int OverlapMapping::serializeMatrixStep1ST(const int *nbOfElemsSrc, const int *recvStep0, const int *countStep0, const int *offsStep0, int *&bigArrI, double *&bigArrD, int *count, int *offsets, @@ -322,9 +361,9 @@ int OverlapMapping::serializeMatrixStep1ST(const int *nbOfElemsSrc, const int *r { if(_source_proc_id_st[i]==myProcId) { - const std::vector< std::map >& mat=_matrixes_st[i]; + const std::vector< SparseDoubleVec >& mat=_matrixes_st[i]; int lgthToSend=0; - for(std::vector< std::map >::const_iterator it=mat.begin();it!=mat.end();it++) + for(std::vector< SparseDoubleVec >::const_iterator it=mat.begin();it!=mat.end();it++) lgthToSend+=(*it).size(); count[_target_proc_id_st[i]]=lgthToSend; fullLgth+=lgthToSend; @@ -341,11 +380,11 @@ int OverlapMapping::serializeMatrixStep1ST(const int *nbOfElemsSrc, const int *r { if(_source_proc_id_st[i]==myProcId) { - const std::vector< std::map >& mat=_matrixes_st[i]; - for(std::vector< std::map >::const_iterator it1=mat.begin();it1!=mat.end();it1++) + const std::vector< SparseDoubleVec >& mat=_matrixes_st[i]; + for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++) { int j=0; - for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++,j++) + for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++,j++) { bigArrI[fullLgth+j]=(*it2).first; bigArrD[fullLgth+j]=(*it2).second; @@ -396,9 +435,10 @@ void OverlapMapping::unserializationST(int nbOfTrgElems, } /*! - * This method should be called when all remote matrix with sourceProcId==thisProcId have been retrieved and are in 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id' - * and 'this->_the_matrix_st_target_ids'. - * This method finish the job of filling 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id' by putting candidates in 'this->_matrixes_st' into them. + * This method should be called when all remote matrix with sourceProcId==thisProcId have been retrieved and are + * in 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id' and 'this->_the_matrix_st_target_ids'. + * This method finish the job of filling 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id' + * by putting candidates in 'this->_matrixes_st' into them (i.e. local computation result). */ void OverlapMapping::finishToFillFinalMatrixST() { @@ -417,64 +457,11 @@ void OverlapMapping::finishToFillFinalMatrixST() for(int i=0;i >& mat=_matrixes_st[i]; + const std::vector& mat=_matrixes_st[i]; _the_matrix_st[j]=mat; _the_matrix_st_source_proc_id.push_back(_source_proc_id_st[i]); j++; } - _matrixes_st.clear(); -} - -/*! - * This method performs the operation of target ids broadcasting. - */ -void OverlapMapping::prepareIdsToSendST() -{ - CommInterface commInterface=_group.getCommInterface(); - const MPIProcessorGroup *group=static_cast(&_group); - const MPI_Comm *comm=group->getComm(); - int grpSize=_group.size(); - _source_ids_to_send_st.clear(); - _source_ids_to_send_st.resize(grpSize); - INTERP_KERNEL::AutoPtr nbsend=new int[grpSize]; - std::fill(nbsend,nbsend+grpSize,0); - for(std::size_t i=0;i<_the_matrix_st_source_proc_id.size();i++) - nbsend[_the_matrix_st_source_proc_id[i]]=_the_matrix_st_source_ids[i].size(); - INTERP_KERNEL::AutoPtr nbrecv=new int[grpSize]; - commInterface.allToAll(nbsend,1,MPI_INT,nbrecv,1,MPI_INT,*comm); - // - INTERP_KERNEL::AutoPtr nbsend2=new int[grpSize]; - std::copy((int *)nbsend,((int *)nbsend)+grpSize,(int *)nbsend2); - INTERP_KERNEL::AutoPtr nbsend3=new int[grpSize]; - nbsend3[0]=0; - for(int i=1;i bigDataSend=new int[sendSz]; - for(std::size_t i=0;i<_the_matrix_st_source_proc_id.size();i++) - { - int offset=nbsend3[_the_matrix_st_source_proc_id[i]]; - std::copy(_the_matrix_st_source_ids[i].begin(),_the_matrix_st_source_ids[i].end(),((int *)nbsend3)+offset); - } - INTERP_KERNEL::AutoPtr nbrecv2=new int[grpSize]; - INTERP_KERNEL::AutoPtr nbrecv3=new int[grpSize]; - std::copy((int *)nbrecv,((int *)nbrecv)+grpSize,(int *)nbrecv2); - nbrecv3[0]=0; - for(int i=1;i bigDataRecv=new int[recvSz]; - // - commInterface.allToAllV(bigDataSend,nbsend2,nbsend3,MPI_INT, - bigDataRecv,nbrecv2,nbrecv3,MPI_INT, - *comm); - for(int i=0;i0) - { - _source_ids_to_send_st[i].insert(_source_ids_to_send_st[i].end(),((int *)bigDataRecv)+nbrecv3[i],((int *)bigDataRecv)+nbrecv3[i]+nbrecv2[i]); - } - } } /*! @@ -499,20 +486,29 @@ void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCoupl nbsend2[0]=0; nbrecv2[0]=0; std::vector valsToSend; + + /* + * FIELD VALUE XCHGE + */ for(int i=0;i::const_iterator isItem1=std::find(_src_proc_st2.begin(),_src_proc_st2.end(),i); + std::vector::const_iterator isItem1=std::find(_sent_src_proc_st2.begin(),_sent_src_proc_st2.end(),i); MEDCouplingAutoRefCountObjectPtr vals; - if(isItem1!=_src_proc_st2.end())//item1 of step2 main algo + if(isItem1!=_sent_src_proc_st2.end()) // source data was sent to proc 'i' { - int id=std::distance(_src_proc_st2.begin(),isItem1); - vals=fieldInput->getArray()->selectByTupleId(_src_ids_st2[id]->getConstPointer(),_src_ids_st2[id]->getConstPointer()+_src_ids_st2[id]->getNumberOfTuples()); + // Prepare local field data to send to proc i + int id=std::distance(_sent_src_proc_st2.begin(),isItem1); + vals=fieldInput->getArray()->selectByTupleId(*_sent_src_ids_st2[id]); } - else + else // no source data was sent to proc 'i' {//item0 of step2 main algo - int id=std::distance(_src_ids_zip_proc_st2.begin(),std::find(_src_ids_zip_proc_st2.begin(),_src_ids_zip_proc_st2.end(),i)); + std::vector::const_iterator isItem22 = std::find(_src_ids_zip_proc_st2.begin(),_src_ids_zip_proc_st2.end(),i ); + int id=std::distance(_src_ids_zip_proc_st2.begin(),isItem22); + if (isItem22 == _src_ids_zip_proc_st2.end()) + std::cout << "(" << myProcId << ") has end iterator!!!!!\n"; vals=fieldInput->getArray()->selectByTupleId(&(_src_ids_zip_st2[id])[0],&(_src_ids_zip_st2[id])[0]+_src_ids_zip_st2[id].size()); } nbsend[i]=vals->getNbOfElems(); @@ -520,8 +516,8 @@ void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCoupl } if(std::find(_proc_ids_to_recv_vector_st.begin(),_proc_ids_to_recv_vector_st.end(),i)!=_proc_ids_to_recv_vector_st.end()) { - std::vector::const_iterator isItem0=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),i); - if(isItem0==_trg_proc_st2.end())//item1 of step2 main algo [ (0,1) computed on proc1 and Matrix-Vector on proc1 ] + std::vector::const_iterator isItem0=std::find(_sent_trg_proc_st2.begin(),_sent_trg_proc_st2.end(),i); + if(isItem0==_sent_trg_proc_st2.end())//item1 of step2 main algo [ (0,1) computed on proc1 and Matrix-Vector on proc1 ] { std::vector::const_iterator it1=std::find(_src_ids_proc_st2.begin(),_src_ids_proc_st2.end(),i); if(it1!=_src_ids_proc_st2.end()) @@ -529,12 +525,12 @@ void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCoupl int id=std::distance(_src_ids_proc_st2.begin(),it1); nbrecv[i]=_nb_of_src_ids_proc_st2[id]*nbOfCompo; } - else if(i==myProcId) + else if(i==myProcId) // diagonal element (i.e. proc #i talking to same proc #i) { - nbrecv[i]=fieldInput->getNumberOfTuplesExpected()*nbOfCompo; + nbrecv[i] = nbsend[i]; } else - throw INTERP_KERNEL::Exception("Plouff ! send email to anthony.geay@cea.fr ! "); + throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): Internal error at field values exchange!"); } else {//item0 of step2 main algo [ (2,1) computed on proc2 but Matrix-Vector on proc1 ] [(1,0) computed on proc1 but Matrix-Vector on proc0] @@ -549,8 +545,19 @@ void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCoupl nbrecv2[i]=nbrecv2[i-1]+nbrecv[i-1]; } INTERP_KERNEL::AutoPtr bigArr=new double[nbrecv2[grpSize-1]+nbrecv[grpSize-1]]; + +// scout << "(" << myProcId << ") nbsend :" << nbsend[0] << "," << nbsend[1] << "," << nbsend[2] << "\n"; +// scout << "(" << myProcId << ") nbrecv :" << nbrecv[0] << "," << nbrecv[1] << "," << nbrecv[2] << "\n"; +// std::cout << scout.str() << "\n"; + /* + * *********************** ALL-TO-ALL + */ commInterface.allToAllV(&valsToSend[0],nbsend,nbsend2,MPI_DOUBLE, bigArr,nbrecv,nbrecv2,MPI_DOUBLE,*comm); + /* + * + * TARGET FIELD COMPUTATION (matrix-vec computation) + */ fieldOutput->getArray()->fillWithZero(); INTERP_KERNEL::AutoPtr tmp=new double[nbOfCompo]; for(int i=0;i >& mat=_the_matrix_st[id]; - const std::vector< std::map >& deno=_the_deno_st[id]; - std::vector::const_iterator isItem0=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),i); - if(isItem0==_trg_proc_st2.end())//item1 of step2 main algo [ (0,1) computed on proc1 and Matrix-Vector on proc1 ] + const std::vector< SparseDoubleVec >& mat=_the_matrix_st[id]; + const std::vector< SparseDoubleVec >& deno=_the_deno_st[id]; + std::vector::const_iterator isItem0=std::find(_sent_trg_proc_st2.begin(),_sent_trg_proc_st2.end(),i); + if(isItem0==_sent_trg_proc_st2.end())//item1 of step2 main algo [ (0,1) computed on proc1 and Matrix-Vector on proc1 ] { int nbOfTrgTuples=mat.size(); for(int j=0;j& mat1=mat[j]; - const std::map& deno1=deno[j]; - std::map::const_iterator it4=deno1.begin(); - for(std::map::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it4++) + const SparseDoubleVec& mat1=mat[j]; + const SparseDoubleVec& deno1=deno[j]; + SparseDoubleVec::const_iterator it4=deno1.begin(); + for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it4++) { - std::transform(bigArr+nbrecv2[i]+((*it3).first)*nbOfCompo,bigArr+nbrecv2[i]+((*it3).first+1)*(nbOfCompo),(double *)tmp,std::bind2nd(std::multiplies(),(*it3).second/(*it4).second)); + std::transform(bigArr+nbrecv2[i]+((*it3).first)*nbOfCompo, + bigArr+nbrecv2[i]+((*it3).first+1)*(nbOfCompo), + (double *)tmp, + std::bind2nd(std::multiplies(),(*it3).second/(*it4).second)); std::transform((double *)tmp,(double *)tmp+nbOfCompo,pt,pt,std::plus()); } } @@ -589,21 +599,24 @@ void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCoupl int newId=0; for(std::vector::const_iterator it=zipIds.begin();it!=zipIds.end();it++,newId++) zipCor[*it]=newId; - int id2=std::distance(_trg_proc_st2.begin(),std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),i)); - const DataArrayInt *tgrIds=_trg_ids_st2[id2]; + int id2=std::distance(_sent_trg_proc_st2.begin(),std::find(_sent_trg_proc_st2.begin(),_sent_trg_proc_st2.end(),i)); + const DataArrayInt *tgrIds=_sent_trg_ids_st2[id2]; const int *tgrIds2=tgrIds->getConstPointer(); int nbOfTrgTuples=mat.size(); for(int j=0;j& mat1=mat[j]; - const std::map& deno1=deno[j]; - std::map::const_iterator it5=deno1.begin(); - for(std::map::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++) + const SparseDoubleVec& mat1=mat[j]; + const SparseDoubleVec& deno1=deno[j]; + SparseDoubleVec::const_iterator it5=deno1.begin(); + for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++) { std::map::const_iterator it4=zipCor.find((*it3).first); if(it4==zipCor.end()) - throw INTERP_KERNEL::Exception("Hmmmmm send e mail to anthony.geay@cea.fr !"); - std::transform(bigArr+nbrecv2[i]+((*it4).second)*nbOfCompo,bigArr+nbrecv2[i]+((*it4).second+1)*(nbOfCompo),(double *)tmp,std::bind2nd(std::multiplies(),(*it3).second/(*it5).second)); + throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): Internal error in final value computation!"); + std::transform(bigArr+nbrecv2[i]+((*it4).second)*nbOfCompo, + bigArr+nbrecv2[i]+((*it4).second+1)*(nbOfCompo), + (double *)tmp, + std::bind2nd(std::multiplies(),(*it3).second/(*it5).second)); std::transform((double *)tmp,(double *)tmp+nbOfCompo,pt+tgrIds2[j]*nbOfCompo,pt+tgrIds2[j]*nbOfCompo,std::plus()); } } @@ -621,11 +634,16 @@ void OverlapMapping::transposeMultiply(const MEDCouplingFieldDouble *fieldInput, } /*! - * This method should be called immediately after _the_matrix_st has been filled with remote computed matrix put in this proc for Matrix-Vector. - * This method computes for these matrix the minimal set of source ids corresponding to the source proc id. + * This method should be called immediately after _the_matrix_st has been filled with remote computed matrix + * put in this proc for Matrix-Vector. + * This method computes for these matrix the minimal set of source ids corresponding to the source proc id. + * */ void OverlapMapping::updateZipSourceIdsForFuture() { + /* When it is called, only the bits received from other processors (i.e. the remotely executed jobs) are in the + big matrix _the_matrix_st. */ + CommInterface commInterface=_group.getCommInterface(); int myProcId=_group.myRank(); int nbOfMatrixRecveived=_the_matrix_st_source_proc_id.size(); @@ -634,20 +652,22 @@ void OverlapMapping::updateZipSourceIdsForFuture() int curSrcProcId=_the_matrix_st_source_proc_id[i]; if(curSrcProcId!=myProcId) { - const std::vector< std::map >& mat=_the_matrix_st[i]; + const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i]; _src_ids_zip_proc_st2.push_back(curSrcProcId); _src_ids_zip_st2.resize(_src_ids_zip_st2.size()+1); std::set s; - for(std::vector< std::map >::const_iterator it1=mat.begin();it1!=mat.end();it1++) - for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) + for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++) + for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) s.insert((*it2).first); _src_ids_zip_st2.back().insert(_src_ids_zip_st2.back().end(),s.begin(),s.end()); + +// std::stringstream scout; +// scout << "(" << myProcId << ") pushed " << _src_ids_zip_st2.back().size() << " values for tgt proc " << curSrcProcId; +// std::cout << scout.str() << "\n"; } } } -// #include - // void OverlapMapping::printTheMatrix() const // { // CommInterface commInterface=_group.getCommInterface(); @@ -655,19 +675,55 @@ void OverlapMapping::updateZipSourceIdsForFuture() // const MPI_Comm *comm=group->getComm(); // int grpSize=_group.size(); // int myProcId=_group.myRank(); -// std::cerr << "I am proc #" << myProcId << std::endl; +// std::stringstream oscerr; // int nbOfMat=_the_matrix_st.size(); -// std::cerr << "I do manage " << nbOfMat << "matrix : "<< std::endl; +// oscerr << "(" << myProcId << ") I hold " << nbOfMat << " matrix(ces) : "<< std::endl; // for(int i=0;i >& locMat=_the_matrix_st[i]; -// for(std::vector< std::map >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++) +// oscerr << " - Matrix #" << i << " coming from source proc #" << _the_matrix_st_source_proc_id[i] << ": "; +// const std::vector< SparseDoubleVec >& locMat=_the_matrix_st[i]; +// int j = 0; +// for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++) // { -// for(std::map::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) -// std::cerr << "(" << (*it2).first << "," << (*it2).second << "), "; -// std::cerr << std::endl; +// oscerr << " Target Cell #" << j; +// for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) +// oscerr << " (" << (*it2).first << "," << (*it2).second << "), "; +// oscerr << std::endl; // } // } -// std::cerr << "*********" << std::endl; +// oscerr << "*********" << std::endl; +// +// // Hope this will be flushed in one go: +// std::cerr << oscerr.str() << std::endl; +//// if(myProcId != 0) +//// MPI_Barrier(MPI_COMM_WORLD); // } +// +// void OverlapMapping::printMatrixesST() const +// { +// CommInterface commInterface=_group.getCommInterface(); +// const MPIProcessorGroup *group=static_cast(&_group); +// const MPI_Comm *comm=group->getComm(); +// int grpSize=_group.size(); +// int myProcId=_group.myRank(); +// std::stringstream oscerr; +// int nbOfMat=_matrixes_st.size(); +// oscerr << "(" << myProcId << ") I hold " << nbOfMat << " LOCAL matrix(ces) : "<< std::endl; +// for(int i=0;i& locMat=_matrixes_st[i]; +// int j = 0; +// for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++) +// { +// oscerr << " Target Cell #" << j; +// for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++) +// oscerr << " (" << (*it2).first << "," << (*it2).second << "), "; +// oscerr << std::endl; +// } +// } +// oscerr << "*********" << std::endl; +// +// // Hope this will be flushed in one go: +// std::cerr << oscerr.str() << std::endl; +// } diff --git a/src/ParaMEDMEM/OverlapMapping.hxx b/src/ParaMEDMEM/OverlapMapping.hxx index cfb06b1bb..2a5b32e4b 100644 --- a/src/ParaMEDMEM/OverlapMapping.hxx +++ b/src/ParaMEDMEM/OverlapMapping.hxx @@ -32,7 +32,9 @@ namespace ParaMEDMEM class DataArrayInt; class MEDCouplingFieldDouble; - /* + typedef std::map SparseDoubleVec; + + /*! * Internal class, not part of the public API. * * Used by the impl of OverlapInterpolationMatrix, plays an equivalent role than what the NxM_Mapping @@ -42,17 +44,19 @@ namespace ParaMEDMEM class OverlapMapping { public: + OverlapMapping(const ProcessorGroup& group); void keepTracksOfSourceIds(int procId, DataArrayInt *ids); void keepTracksOfTargetIds(int procId, DataArrayInt *ids); - void addContributionST(const std::vector< std::map >& matrixST, const DataArrayInt *srcIds, int srcProcId, const DataArrayInt *trgIds, int trgProcId); - void prepare(const std::vector< std::vector >& procsInInteraction, int nbOfTrgElems); + void addContributionST(const std::vector< SparseDoubleVec >& matrixST, const DataArrayInt *srcIds, int srcProcId, const DataArrayInt *trgIds, int trgProcId); + void prepare(const std::vector< int >& procsToSendField, int nbOfTrgElems); void computeDenoConservativeVolumic(int nbOfTuplesTrg); void computeDenoGlobConstraint(); // void multiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput) const; void transposeMultiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput); private: + void fillProcToSendRcvForMultiply(const std::vector< int >& procsToSendField); void serializeMatrixStep0ST(const int *nbOfElemsSrc, int *&bigArr, int *count, int *offsets, int *countForRecv, int *offsetsForRecv) const; int serializeMatrixStep1ST(const int *nbOfElemsSrc, const int *recvStep0, const int *countStep0, const int *offsStep0, @@ -61,36 +65,65 @@ namespace ParaMEDMEM void unserializationST(int nbOfTrgElems, const int *nbOfElemsSrcPerProc, const int *bigArrRecv, const int *bigArrRecvCounts, const int *bigArrRecvOffs, const int *bigArrRecv2, const double *bigArrDRecv2, const int *bigArrRecv2Count, const int *bigArrRecv2Offs); void finishToFillFinalMatrixST(); - void prepareIdsToSendST(); void updateZipSourceIdsForFuture(); - //void printTheMatrix() const; + + // Debug +// void printMatrixesST() const; +// void printTheMatrix() const; private: const ProcessorGroup &_group; - //! vector of ids - std::vector< MEDCouplingAutoRefCountObjectPtr > _src_ids_st2;//item #1 - std::vector< int > _src_proc_st2;//item #1 - std::vector< MEDCouplingAutoRefCountObjectPtr > _trg_ids_st2;//item #0 - std::vector< int > _trg_proc_st2;//item #0 - std::vector< int > _nb_of_src_ids_proc_st2;//item #1 - std::vector< int > _src_ids_proc_st2;//item #1 - std::vector< std::vector > _src_ids_zip_st2;//same size as _src_ids_zip_proc_st2. Sorted. specifies for each id the corresponding ids to send. This is for item0 of Step2 of main algorithm - std::vector< int > _src_ids_zip_proc_st2; - //! vector of matrixes the first entry correspond to source proc id in _source_ids_st - std::vector< std::vector< std::map > > _matrixes_st; - std::vector< std::vector > _source_ids_st; + /**! Vector of DAInt of cell identifiers. The 2 following class members work in pair. For a proc ID i, + * first member gives an old2new map for the local part of the source mesh that has been sent. + * Second member gives proc ID. */ + std::vector< MEDCouplingAutoRefCountObjectPtr > _sent_src_ids_st2; + //! see above _sent_src_ids_st2 + std::vector< int > _sent_src_proc_st2; + + //! See _src_ids_st2 and _sent_src_proc_st2. Same for target mesh. + std::vector< MEDCouplingAutoRefCountObjectPtr > _sent_trg_ids_st2; + //! See _src_ids_st2 and _sent_src_proc_st2. Same for target mesh. + std::vector< int > _sent_trg_proc_st2; + + + /**! Vector of matrixes (partial interpolation ratios), result of the local interpolator run. + * Indexing shared with _source_proc_id_st, and _target_proc_id_st. */ + std::vector< std::vector< SparseDoubleVec > > _matrixes_st; + //! See _matrixes_st - vec of source proc IDs std::vector< int > _source_proc_id_st; - std::vector< std::vector > _target_ids_st; + //! See _matrixes_st - vec of target proc IDs std::vector< int > _target_proc_id_st; - //! the matrix for matrix-vector product. The first dimension the set of target procs that interacts with local source mesh. The second dimension correspond to nb of local source ids. - std::vector< std::vector< std::map > > _the_matrix_st; + + //! Vector of remote remote proc IDs for source mesh. Indexing shared with _nb_of_src_ids_proc_st2 + std::vector< int > _src_ids_proc_st2; + //! Number of cells in the mesh/mapping received from the remote proc i for source mesh. See _src_ids_proc_st2 above + std::vector< int > _nb_of_src_ids_proc_st2; + + /**! Specifies for each remote proc ID (given in _src_ids_zip_proc_st2 below) the corresponding local + * source cell IDs to use/send. Same indexing as _src_ids_zip_proc_st2. Sorted. + * On a given proc, those two members contain exactly the same set of cell identifiers as what is given + * in the locally held interpolation matrices. */ + std::vector< std::vector > _src_ids_zip_st2; + //! Vector of remote proc ID to which the local source mapping above corresponds. See _src_ids_zip_st2 above. + std::vector< int > _src_ids_zip_proc_st2; + + /**! THE matrix for matrix-vector product. The first dimension is indexed in the set of target procs + * that interacts with local source mesh. The second dim is the pseudo id of source proc. + * Same indexing as _the_matrix_st_source_proc_id */ + std::vector< std::vector< SparseDoubleVec > > _the_matrix_st; + //! See _the_matrix_st above. List of source proc IDs contributing to _the_matrix_st std::vector< int > _the_matrix_st_source_proc_id; - std::vector< std::vector > _the_matrix_st_source_ids; - std::vector< std::vector< std::map > > _the_deno_st; - //! this attribute stores the proc ids that wait for data from this proc ids for matrix-vector computation + + //! Proc IDs to which data will be sent (originating this current proc) for matrix-vector computation std::vector< int > _proc_ids_to_send_vector_st; + //! Proc IDs from which data will be received (on this current proc) for matrix-vector computation std::vector< int > _proc_ids_to_recv_vector_st; - //! this attribute is of size _group.size(); for each procId in _group _source_ids_to_send_st[procId] contains tupleId to send abroad - std::vector< std::vector > _source_ids_to_send_st; + + // Denominators (computed from the numerator matrix) + std::vector< std::vector< SparseDoubleVec > > _the_deno_st; + +// std::vector< std::vector > _the_matrix_st_source_ids; +// //! this attribute is of size _group.size(); for each procId in _group _source_ids_to_send_st[procId] contains tupleId to send abroad +// std::vector< std::vector > _source_ids_to_send_st; }; } diff --git a/src/ParaMEDMEMTest/CMakeLists.txt b/src/ParaMEDMEMTest/CMakeLists.txt index 513031360..4d021f565 100644 --- a/src/ParaMEDMEMTest/CMakeLists.txt +++ b/src/ParaMEDMEMTest/CMakeLists.txt @@ -68,7 +68,7 @@ SET(ParaMEDMEMTest_SOURCES ADD_LIBRARY(ParaMEDMEMTest SHARED ${ParaMEDMEMTest_SOURCES}) SET_TARGET_PROPERTIES(ParaMEDMEMTest PROPERTIES COMPILE_FLAGS "") -TARGET_LINK_LIBRARIES(ParaMEDMEMTest paramedmem paramedloader ${CPPUNIT_LIBRARIES}) +TARGET_LINK_LIBRARIES(ParaMEDMEMTest paramedmem paramedloader medcouplingremapper ${CPPUNIT_LIBRARIES}) INSTALL(TARGETS ParaMEDMEMTest DESTINATION ${SALOME_INSTALL_LIBS}) SET(TESTSParaMEDMEM) diff --git a/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx b/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx index a8bf2b474..54db5e2e0 100644 --- a/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx +++ b/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx @@ -48,7 +48,10 @@ class ParaMEDMEMTest : public CppUnit::TestFixture CPPUNIT_TEST(testInterpKernelDECPartialProcs); CPPUNIT_TEST(testInterpKernelDEC3DSurfEmptyBBox); CPPUNIT_TEST(testOverlapDEC1); - + CPPUNIT_TEST(testOverlapDEC2); +// CPPUNIT_TEST(testOverlapDEC_LMEC_seq); +// CPPUNIT_TEST(testOverlapDEC_LMEC_para); +// CPPUNIT_TEST(testSynchronousEqualInterpKernelWithoutInterpNativeDEC_2D); CPPUNIT_TEST(testSynchronousEqualInterpKernelWithoutInterpDEC_2D); CPPUNIT_TEST(testSynchronousEqualInterpKernelDEC_2D); @@ -61,11 +64,11 @@ class ParaMEDMEMTest : public CppUnit::TestFixture CPPUNIT_TEST(testAsynchronousSlowerSourceInterpKernelDEC_2D); CPPUNIT_TEST(testAsynchronousSlowSourceInterpKernelDEC_2D); CPPUNIT_TEST(testAsynchronousFastSourceInterpKernelDEC_2D); -#ifdef MED_ENABLE_FVM - //can be added again after FVM correction for 2D - // CPPUNIT_TEST(testNonCoincidentDEC_2D); - CPPUNIT_TEST(testNonCoincidentDEC_3D); -#endif +//#ifdef MED_ENABLE_FVM +// //can be added again after FVM correction for 2D +// // CPPUNIT_TEST(testNonCoincidentDEC_2D); +// CPPUNIT_TEST(testNonCoincidentDEC_3D); +//#endif CPPUNIT_TEST(testStructuredCoincidentDEC); CPPUNIT_TEST(testStructuredCoincidentDEC); CPPUNIT_TEST(testICoco1); @@ -104,6 +107,9 @@ public: void testInterpKernelDECPartialProcs(); void testInterpKernelDEC3DSurfEmptyBBox(); void testOverlapDEC1(); + void testOverlapDEC2(); + void testOverlapDEC_LMEC_seq(); + void testOverlapDEC_LMEC_para(); #ifdef MED_ENABLE_FVM void testNonCoincidentDEC_2D(); void testNonCoincidentDEC_3D(); @@ -155,6 +161,7 @@ private: void testInterpKernelDEC_2D_(const char *srcMeth, const char *targetMeth); void testInterpKernelDEC2_2D_(const char *srcMeth, const char *targetMeth); void testInterpKernelDEC_3D_(const char *srcMeth, const char *targetMeth); + }; // to automatically remove temporary files from disk diff --git a/src/ParaMEDMEMTest/ParaMEDMEMTest_Gauthier1.cxx b/src/ParaMEDMEMTest/ParaMEDMEMTest_Gauthier1.cxx index cc97ede18..e11b5a944 100644 --- a/src/ParaMEDMEMTest/ParaMEDMEMTest_Gauthier1.cxx +++ b/src/ParaMEDMEMTest/ParaMEDMEMTest_Gauthier1.cxx @@ -249,6 +249,7 @@ void ParaMEDMEMTest::testGauthier1() void ParaMEDMEMTest::testGauthier2() { + std::cout << "testGauthier2\n"; double valuesExpected1[2]={0.,0.}; double valuesExpected2[2]={0.95,0.970625}; diff --git a/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx b/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx index c7a610461..c1f5bb92f 100644 --- a/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx +++ b/src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx @@ -33,42 +33,166 @@ #include -void ParaMEDMEMTest::testOverlapDEC1() -{ - std::string srcM("P0"); - std::string targetM("P0"); - int size; - int rank; - MPI_Comm_size(MPI_COMM_WORLD,&size); - MPI_Comm_rank(MPI_COMM_WORLD,&rank); +using namespace std; - if (size != 3) return ; - - int nproc = 3; - std::set procs; - - for (int i=0; i MUMesh; +typedef MEDCouplingAutoRefCountObjectPtr MFDouble; + +//void ParaMEDMEMTest::testOverlapDEC_LMEC_seq() +//{ +// // T_SC_Trio_src.med -- "SupportOf_" +// // T_SC_Trio_dst.med -- "SupportOf_T_SC_Trio" +// // h_TH_Trio_src.med -- "SupportOf_" +// // h_TH_Trio_dst.med -- "SupportOf_h_TH_Trio" +// string rep("/export/home/adrien/support/antoine_LMEC/"); +// string src_mesh_nam(rep + string("T_SC_Trio_src.med")); +// string tgt_mesh_nam(rep + string("T_SC_Trio_dst.med")); +//// string src_mesh_nam(rep + string("h_TH_Trio_src.med")); +//// string tgt_mesh_nam(rep + string("h_TH_Trio_dst.med")); +// MUMesh src_mesh=MEDLoader::ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0); +// MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0); +//// MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_h_TH_Trio",0); +// +// MFDouble srcField = MEDCouplingFieldDouble::New(ON_CELLS, ONE_TIME); +// srcField->setMesh(src_mesh); +// DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1); +// dad->fillWithValue(1.0); +// srcField->setArray(dad); +// srcField->setNature(ConservativeVolumic); +// +// MEDCouplingRemapper remap; +// remap.setOrientation(2); // always consider surface intersections as absolute areas. +// remap.prepare(src_mesh, tgt_mesh, "P0P0"); +// MFDouble tgtField = remap.transferField(srcField, 1.0e+300); +// tgtField->setName("result"); +// string out_nam(rep + string("adrien.med")); +// MEDLoader::WriteField(out_nam,tgtField, true); +// cout << "wrote: " << out_nam << "\n"; +// double integ1 = 0.0, integ2 = 0.0; +// srcField->integral(true, &integ1); +// tgtField->integral(true, &integ2); +//// tgtField->reprQuickOverview(cout); +// CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8); +// +// dad->decrRef(); +//} +// +//void ParaMEDMEMTest::testOverlapDEC_LMEC_para() +//{ +// using namespace ParaMEDMEM; +// +// int size; +// int rank; +// MPI_Comm_size(MPI_COMM_WORLD,&size); +// MPI_Comm_rank(MPI_COMM_WORLD,&rank); +// +// if (size != 1) return ; +// +// int nproc = 1; +// std::set procs; +// +// for (int i=0; isetMesh(src_mesh); +// DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1); +// dad->fillWithValue(1.0); +// srcField->setArray(dad); +// srcField->setNature(ConservativeVolumic); +// +// ComponentTopology comptopo; +// parameshS = new ParaMESH(src_mesh,*dec.getGroup(),"source mesh"); +// parafieldS = new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo); +// parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint +// parafieldS->getField()->setArray(dad); +// +// // **** TARGET +// parameshT=new ParaMESH(tgt_mesh,*dec.getGroup(),"target mesh"); +// parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo); +// parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint +// parafieldT->getField()->getArray()->fillWithValue(1.0e300); +//// valsT[0]=7.; +// } +// dec.setOrientation(2); +// dec.attachSourceLocalField(parafieldS); +// dec.attachTargetLocalField(parafieldT); +// dec.synchronize(); +// dec.sendRecvData(true); +// // +// if(rank==0) +// { +// double integ1 = 0.0, integ2 = 0.0; +// MEDCouplingFieldDouble * tgtField; +// +// srcField->integral(true, &integ1); +// tgtField = parafieldT->getField(); +//// tgtField->reprQuickOverview(cout); +// tgtField->integral(true, &integ2); +// tgtField->setName("result"); +// string out_nam(rep + string("adrien_para.med")); +// MEDLoader::WriteField(out_nam,tgtField, true); +// cout << "wrote: " << out_nam << "\n"; +// CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8); +// } +// delete parafieldS; +// delete parafieldT; +// delete parameshS; +// delete parameshT; +// +// MPI_Barrier(MPI_COMM_WORLD); +//} + +void prepareData1(int rank, ProcessorGroup * grp, + MEDCouplingUMesh *& meshS, MEDCouplingUMesh *& meshT, + ParaMESH*& parameshS, ParaMESH*& parameshT, + ParaFIELD*& parafieldS, ParaFIELD*& parafieldT) +{ if(rank==0) { const double coordsS[10]={0.,0.,0.5,0.,1.,0.,0.,0.5,0.5,0.5}; const double coordsT[6]={0.,0.,1.,0.,1.,1.}; - meshS=ParaMEDMEM::MEDCouplingUMesh::New(); + meshS=MEDCouplingUMesh::New(); meshS->setMeshDimension(2); - ParaMEDMEM::DataArrayDouble *myCoords=ParaMEDMEM::DataArrayDouble::New(); + DataArrayDouble *myCoords=DataArrayDouble::New(); myCoords->alloc(5,2); std::copy(coordsS,coordsS+10,myCoords->getPointer()); meshS->setCoords(myCoords); @@ -78,16 +202,16 @@ void ParaMEDMEMTest::testOverlapDEC1() meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS); meshS->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connS+4); meshS->finishInsertingCells(); - ParaMEDMEM::ComponentTopology comptopo; - parameshS=new ParaMEDMEM::ParaMESH(meshS,*dec.getGrp(),"source mesh"); - parafieldS=new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::NO_TIME,parameshS,comptopo); - parafieldS->getField()->setNature(ParaMEDMEM::ConservativeVolumic);//IntegralGlobConstraint + ComponentTopology comptopo; + parameshS=new ParaMESH(meshS, *grp,"source mesh"); + parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo); + parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint double *valsS=parafieldS->getField()->getArray()->getPointer(); valsS[0]=7.; valsS[1]=8.; // - meshT=ParaMEDMEM::MEDCouplingUMesh::New(); + meshT=MEDCouplingUMesh::New(); meshT->setMeshDimension(2); - myCoords=ParaMEDMEM::DataArrayDouble::New(); + myCoords=DataArrayDouble::New(); myCoords->alloc(3,2); std::copy(coordsT,coordsT+6,myCoords->getPointer()); meshT->setCoords(myCoords); @@ -96,9 +220,9 @@ void ParaMEDMEMTest::testOverlapDEC1() meshT->allocateCells(1); meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT); meshT->finishInsertingCells(); - parameshT=new ParaMEDMEM::ParaMESH(meshT,*dec.getGrp(),"target mesh"); - parafieldT=new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::NO_TIME,parameshT,comptopo); - parafieldT->getField()->setNature(ParaMEDMEM::ConservativeVolumic);//IntegralGlobConstraint + parameshT=new ParaMESH(meshT,*grp,"target mesh"); + parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo); + parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint double *valsT=parafieldT->getField()->getArray()->getPointer(); valsT[0]=7.; } @@ -107,9 +231,9 @@ void ParaMEDMEMTest::testOverlapDEC1() { const double coordsS[10]={1.,0.,0.5,0.5,1.,0.5,0.5,1.,1.,1.}; const double coordsT[6]={0.,0.,0.5,0.5,0.,1.}; - meshS=ParaMEDMEM::MEDCouplingUMesh::New(); + meshS=MEDCouplingUMesh::New(); meshS->setMeshDimension(2); - ParaMEDMEM::DataArrayDouble *myCoords=ParaMEDMEM::DataArrayDouble::New(); + DataArrayDouble *myCoords=DataArrayDouble::New(); myCoords->alloc(5,2); std::copy(coordsS,coordsS+10,myCoords->getPointer()); meshS->setCoords(myCoords); @@ -119,16 +243,17 @@ void ParaMEDMEMTest::testOverlapDEC1() meshS->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connS); meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS+3); meshS->finishInsertingCells(); - ParaMEDMEM::ComponentTopology comptopo; - parameshS=new ParaMEDMEM::ParaMESH(meshS,*dec.getGrp(),"source mesh"); - parafieldS=new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::NO_TIME,parameshS,comptopo); - parafieldS->getField()->setNature(ParaMEDMEM::ConservativeVolumic);//IntegralGlobConstraint + ComponentTopology comptopo; + parameshS=new ParaMESH(meshS,*grp,"source mesh"); + parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo); + parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint double *valsS=parafieldS->getField()->getArray()->getPointer(); - valsS[0]=9.; valsS[1]=11.; + valsS[0]=9.; + valsS[1]=11.; // - meshT=ParaMEDMEM::MEDCouplingUMesh::New(); + meshT=MEDCouplingUMesh::New(); meshT->setMeshDimension(2); - myCoords=ParaMEDMEM::DataArrayDouble::New(); + myCoords=DataArrayDouble::New(); myCoords->alloc(3,2); std::copy(coordsT,coordsT+6,myCoords->getPointer()); meshT->setCoords(myCoords); @@ -137,9 +262,9 @@ void ParaMEDMEMTest::testOverlapDEC1() meshT->allocateCells(1); meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT); meshT->finishInsertingCells(); - parameshT=new ParaMEDMEM::ParaMESH(meshT,*dec.getGrp(),"target mesh"); - parafieldT=new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::NO_TIME,parameshT,comptopo); - parafieldT->getField()->setNature(ParaMEDMEM::ConservativeVolumic);//IntegralGlobConstraint + parameshT=new ParaMESH(meshT,*grp,"target mesh"); + parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo); + parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint double *valsT=parafieldT->getField()->getArray()->getPointer(); valsT[0]=8.; } @@ -148,9 +273,9 @@ void ParaMEDMEMTest::testOverlapDEC1() { const double coordsS[8]={0.,0.5, 0.5,0.5, 0.,1., 0.5,1.}; const double coordsT[6]={0.5,0.5,0.,1.,1.,1.}; - meshS=ParaMEDMEM::MEDCouplingUMesh::New(); + meshS=MEDCouplingUMesh::New(); meshS->setMeshDimension(2); - ParaMEDMEM::DataArrayDouble *myCoords=ParaMEDMEM::DataArrayDouble::New(); + DataArrayDouble *myCoords=DataArrayDouble::New(); myCoords->alloc(4,2); std::copy(coordsS,coordsS+8,myCoords->getPointer()); meshS->setCoords(myCoords); @@ -159,16 +284,16 @@ void ParaMEDMEMTest::testOverlapDEC1() meshS->allocateCells(1); meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS); meshS->finishInsertingCells(); - ParaMEDMEM::ComponentTopology comptopo; - parameshS=new ParaMEDMEM::ParaMESH(meshS,*dec.getGrp(),"source mesh"); - parafieldS=new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::NO_TIME,parameshS,comptopo); - parafieldS->getField()->setNature(ParaMEDMEM::ConservativeVolumic);//IntegralGlobConstraint + ComponentTopology comptopo; + parameshS=new ParaMESH(meshS,*grp,"source mesh"); + parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo); + parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint double *valsS=parafieldS->getField()->getArray()->getPointer(); valsS[0]=10.; // - meshT=ParaMEDMEM::MEDCouplingUMesh::New(); + meshT=MEDCouplingUMesh::New(); meshT->setMeshDimension(2); - myCoords=ParaMEDMEM::DataArrayDouble::New(); + myCoords=DataArrayDouble::New(); myCoords->alloc(3,2); std::copy(coordsT,coordsT+6,myCoords->getPointer()); meshT->setCoords(myCoords); @@ -177,12 +302,112 @@ void ParaMEDMEMTest::testOverlapDEC1() meshT->allocateCells(1); meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT); meshT->finishInsertingCells(); - parameshT=new ParaMEDMEM::ParaMESH(meshT,*dec.getGrp(),"target mesh"); - parafieldT=new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::NO_TIME,parameshT,comptopo); - parafieldT->getField()->setNature(ParaMEDMEM::ConservativeVolumic);//IntegralGlobConstraint + parameshT=new ParaMESH(meshT,*grp,"target mesh"); + parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo); + parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint double *valsT=parafieldT->getField()->getArray()->getPointer(); valsT[0]=9.; } + +} + +/*! Test case from the official doc of the OverlapDEC. + * WARNING: bounding boxes are tweaked here to make the case more interesting (i.e. to avoid an all to all exchange + * between all procs). + */ +void ParaMEDMEMTest::testOverlapDEC1() +{ + int size, rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + // char hostname[256]; + // printf("(%d) PID %d on localhost ready for attach\n", rank, getpid()); + // fflush(stdout); + + if (size != 3) return ; + + int nproc = 3; + std::set procs; + + for (int i=0; igetField()->getArray()->getIJ(0,0),1e-12); + } + if(rank==1) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12); + } + if(rank==2) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12); + } + delete parafieldS; + delete parafieldT; + delete parameshS; + delete parameshT; + meshS->decrRef(); + meshT->decrRef(); + + MPI_Barrier(MPI_COMM_WORLD); +} + +/*! + * Same as testOverlapDEC1() but with regular bounding boxes. If you're looking for a nice debug case, + * testOverlapDEC1() is more appropriate. + */ +void ParaMEDMEMTest::testOverlapDEC2() +{ + int size, rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + + if (size != 3) return ; + + int nproc = 3; + std::set procs; + + for (int i=0; i