From 4bb7a01cda6b7dc506a76b7f3da4b1ccf25d67df Mon Sep 17 00:00:00 2001 From: vbd Date: Wed, 1 Aug 2007 06:36:59 +0000 Subject: [PATCH] adding the possibility to send data from the target (idle) side and receiving it from the source (working) side. --- src/ParaMEDMEM/InterpolationMatrix.cxx | 322 +++++++++++++++---------- src/ParaMEDMEM/InterpolationMatrix.hxx | 16 +- src/ParaMEDMEM/IntersectionDEC.cxx | 108 +++++---- src/ParaMEDMEM/IntersectionDEC.hxx | 2 +- src/ParaMEDMEM/MxN_Mapping.cxx | 87 +++++++ src/ParaMEDMEM/MxN_Mapping.hxx | 1 + 6 files changed, 350 insertions(+), 186 deletions(-) diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index 1452d9aee..58f9b406c 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -6,10 +6,12 @@ #include "Interpolation3DSurf.hxx" /*! \class InterpolationMatrix + This class enables the storage of an interpolation matrix Wij mapping -source field Sj to target field Ti via Ti=Wij.Sj. +source field Sj to target field Ti via Ti=Vi^(-1).Wij.Sj. The matrix is built and stored on the processors belonging to the source group. + */ namespace ParaMEDMEM @@ -22,22 +24,26 @@ namespace ParaMEDMEM \param local_group processor group containing the local processors \param distant_group processor group containing the distant processors \param method interpolation method - */ + */ InterpolationMatrix::InterpolationMatrix( const ParaMEDMEM::ParaMESH& source_support, -const ProcessorGroup& local_group, -const ProcessorGroup& distant_group, +const ProcessorGroup& source_group, +const ProcessorGroup& target_group, const string& method): -_source_support(*source_support.getMesh()), -_local_group(local_group), -_distant_group(distant_group), -_mapping(local_group, distant_group), -_method(method) + _source_support(*source_support.getMesh()), + _mapping(source_group, target_group), + _method(method), + _source_group(source_group), + _target_group(target_group), + _source_volume(0) + { int nbelems= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); - _row_offsets.resize(nbelems+1,0); - + _row_offsets.resize(nbelems+1); + for (int i=0; i > surfaces = interpolator->interpol_maillages(distant_support,_source_support); - delete interpolator; + delete interpolator; + if (surfaces.size() != source_size) + throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix"); + + //computing the vectors containing the source and target element volumes + MEDMEM::SUPPORT target_support(&distant_support,"all cells", MED_EN::MED_CELL); + MEDMEM::FIELD* target_triangle_surf = distant_support.getArea(&target_support); + MEDMEM::SUPPORT source_support (const_cast(&_source_support),"all cells", MED_EN::MED_CELL); + MEDMEM::FIELD* source_triangle_surf = _source_support.getArea(&source_support); + + //storing the source volumes + _source_volume.resize(source_size); + for (int i=0; igetValueIJ(i+1,1); + + //loop over the elements to build the interpolation + //matrix structures + for (int ielem=0; ielem < surfaces.size(); ielem++) + { + _row_offsets[ielem+1]+=surfaces[ielem].size(); + // _source_indices.push_back(make_pair(iproc_distant, ielem)); + for (map::const_iterator iter=surfaces[ielem].begin(); + iter!=surfaces[ielem].end(); + iter++) + { + //surface of the target triangle + double surf = target_triangle_surf->getValueIJ(iter->first,1); + //locating the (iproc, itriangle) pair in the list of columns + vector >::iterator iter2 = find (_col_offsets.begin(), _col_offsets.end(),make_pair(iproc_distant,iter->first)); + int col_id; + + + if (iter2==_col_offsets.end()) + { + //(iproc, itriangle) is not registered in the list + //of distant elements - if (surfaces.size() != nbelems) - throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix"); + _col_offsets.push_back(make_pair(iproc_distant,iter->first)); + col_id =_col_offsets.size(); + _mapping.addElementFromSource(iproc_distant,distant_elems[iter->first-1]); + _target_volume.push_back(surf); + } + else + { + col_id=iter2-_col_offsets.begin()+1; + } - MEDMEM::SUPPORT support(&distant_support,"all cells", MED_EN::MED_CELL); - MEDMEM::FIELD* triangle_surf = distant_support.getArea(&support); + //the non zero coefficient is stored + //ielem is the row, + //col_id is the number of the column + //iter->second is the value of the coefficient - for (int ielem=0; ielem < surfaces.size(); ielem++) - { - _row_offsets[ielem+1]+=surfaces[ielem].size(); - // _source_indices.push_back(make_pair(iproc_distant, ielem)); - for (map::const_iterator iter=surfaces[ielem].begin(); - iter!=surfaces[ielem].end(); - iter++) - { - //surface of the source triangle - double surf = triangle_surf->getValueIJ(iter->first,1); - //oriented surfaces are not considered - surf = (surf>0)?surf:-surf; - //locating the (iproc, itriangle) pair in the list of columns - vector >::iterator iter2 = find (_col_offsets.begin(), _col_offsets.end(),make_pair(iproc_distant,iter->first)); - int col_id; - - //(iproc, itriangle) is not registered in the list - //of distant elements - if (iter2==_col_offsets.end()) - { - _col_offsets.push_back(make_pair(iproc_distant,iter->first)); - col_id =_col_offsets.size(); - _mapping.addElementFromSource(iproc_distant,distant_elems[iter->first-1]); - } - else - { - col_id=iter2-_col_offsets.begin()+1; - } - //the column indirection number is registered - //with the interpolation coefficient - //actual column number is obtained with _col_offsets[col_id] - _col_numbers.push_back(col_id); - _coeffs.push_back(iter->second/surf); - } - } - delete triangle_surf; + _coeffs[ielem].push_back(make_pair(col_id,iter->second)); + + } + } + delete source_triangle_surf; + delete target_triangle_surf; } - -// intersect (_local_support,_distant_support); -// const int* conn_source_index = _local_support->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_ALL_ELEMENTS); -// const int* conn_source = distant_support->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_ALL_ELEMENTS); -// const double* coords_source = _local_support->getCoords(MED_EN::MED_FULL_INTERLACE); -// for (int ielem=0; ielem< _local_support->getNumberOfElements(MED_EN::MED_CELL); ielem++) -// { -// MED_EN::medGeometryElementType type = distant_support->getElementType(ielem); -// switch (type) -// { -// -// case MED_EN::TRIA3 : -// double* coords=new int[3*dim]; -// for (ielem_distant=0; ielem_distantgetNumberOfElements(MED_EN::MED_CELL);ielem_distant++) -// { -// int ielem_target=*iter; -// -// for (int i=0; i<3; i++) -// { -// int node = conn[conn_source_index[ielem]+i]; -// for (int idim=0; i int(1.dV) -// // 1 -> int(x.dV) -// // 2 -> int(y.dV) -// // 3 -> int(z.dV) -// if (intersection_moments[0]>0) -// { -// _source_indices.push_back(make_pair(iproc_distant, ielem); -// _col_numbers.push_back(ielem_target); -// _row_offsets[irow]++; -// _coeffs.push_back(intersection_moments[0]); -// for (int i=0; i& field) const { - double* target_value = new double[_col_offsets.size()* field.getNumberOfComponents()]; - for (int i=0; i< _col_offsets.size()* field.getNumberOfComponents();i++) - { - target_value[i]=0.0; - } - int nbrows= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); - for (int irow=0; irow target_value(_col_offsets.size()* field.getNumberOfComponents(),0.0); + + //computing the matrix multiply on source side + if (_source_group.containsMyRank()) + { + int nbcomp = field.getNumberOfComponents(); + int nbrows= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + + // performing W.S + // W is the intersection matrix + // S is the source vector + + for (int irow=0; irow0) - delete[] target_value; } - +/*! + \brief performs s=WTt, where t is the target field, s is the source field, WT is the transpose matrix from W + + The call to this method must be called both on the working side +and on the idle side. On the working side, the target vector T is received and the + vector S=VS^(-1).(WT.T) is computed to update the field. +On the idle side, no computation is done, but the field is sent. + + \param field source field on processors involved on the source side, target field on processors on the target side + */ +void InterpolationMatrix::transposeMultiply(MEDMEM::FIELD& field) const +{ + vector source_value(_col_offsets.size()* field.getNumberOfComponents(),0.0); + _mapping.reverseSendRecv(&source_value[0],field); + + //treatment of the transpose matrix multiply on the source side + if (_source_group.containsMyRank()) + { + int nbcomp = field.getNumberOfComponents(); + int nbrows = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + + vector target_value(nbrows,0.0); + + //performing WT.T + //WT is W transpose + //T is the target vector + for (int irow=0; irow&) const; + void transposeMultiply(MEDMEM::FIELD&)const; void prepare(); int getNbRows() const {return _row_offsets.size();} private: // vector > _source_indices; vector _row_offsets; - vector _col_numbers; + //vector _col_numbers; vector > _col_offsets; - vector _coeffs; + // vector _coeffs; const MEDMEM::MESH& _source_support; MxN_Mapping _mapping; string _method; - const ProcessorGroup& _local_group; - const ProcessorGroup& _distant_group; - - + const ProcessorGroup& _source_group; + const ProcessorGroup& _target_group; + vector _target_volume; + vector _source_volume; + vector > > _coeffs; }; } diff --git a/src/ParaMEDMEM/IntersectionDEC.cxx b/src/ParaMEDMEM/IntersectionDEC.cxx index 78a15ab71..142ea8d40 100644 --- a/src/ParaMEDMEM/IntersectionDEC.cxx +++ b/src/ParaMEDMEM/IntersectionDEC.cxx @@ -29,85 +29,93 @@ IntersectionDEC::~IntersectionDEC() { if (_interpolation_matrix !=0) delete _interpolation_matrix; + } /*! Synchronization process for exchanging topologies */ void IntersectionDEC::synchronize() { - + const ParaMEDMEM::ParaMESH* para_mesh = _local_field->getSupport()->getMesh(); + cout <<"size of Interpolation Matrix"<containsMyRank()) { - const ParaMEDMEM::ParaMESH* para_mesh = _source_field->getSupport()->getMesh(); - _interpolation_matrix = new InterpolationMatrix (*para_mesh, *_source_group,*_target_group,"P0"); - //locate the distant meshes ElementLocator locator(*para_mesh, *_target_group); - + MESH* distant_mesh=0; int* distant_ids=0; for (int i=0; i<_target_group->size(); i++) - { - // int idistant_proc = (i+_source_group->myRank())%_target_group->size(); - int idistant_proc=i; - - //gathers pieces of the target meshes that can intersect the local mesh - locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids); - - if (distant_mesh !=0) - { - //adds the contribution of the distant mesh on the local one - int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc); - cout <<"add contribution from proc "<myRank()<addContribution(*distant_mesh,idistant_proc_in_union,distant_ids); - - delete distant_mesh; - delete[] distant_ids; - distant_mesh=0; - distant_ids=0; - } - } - - } - - if (_target_field!=0) + { + // int idistant_proc = (i+_source_group->myRank())%_target_group->size(); + int idistant_proc=i; + + //gathers pieces of the target meshes that can intersect the local mesh + locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids); + + if (distant_mesh !=0) + { + //adds the contribution of the distant mesh on the local one + int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc); + cout <<"add contribution from proc "<myRank()<addContribution(*distant_mesh,idistant_proc_in_union,distant_ids); + + delete distant_mesh; + delete[] distant_ids; + distant_mesh=0; + distant_ids=0; + } + } + + } + + if (_target_group->containsMyRank()) { - const ParaMEDMEM::ParaMESH* para_mesh = _target_field->getSupport()->getMesh(); - _interpolation_matrix = new InterpolationMatrix (*para_mesh, *_source_group,*_target_group,"P0"); - ElementLocator locator(*para_mesh, *_source_group); MESH* distant_mesh=0; int* distant_ids=0; for (int i=0; i<_source_group->size(); i++) - { - // int idistant_proc = (i+_target_group->myRank())%_source_group->size(); - int idistant_proc=i; - //gathers pieces of the target meshes that can intersect the local mesh - locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids); - cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<myRank())%_source_group->size(); + int idistant_proc=i; + //gathers pieces of the target meshes that can intersect the local mesh + locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids); + cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<prepare(); + // _volume_vector=para_mesh->getVolume(); + } void IntersectionDEC::recvData() { - _interpolation_matrix->multiply(*_target_field->getField()); + if (_source_group->containsMyRank()) + _interpolation_matrix->transposeMultiply(*_local_field->getField()); + else if (_target_group->containsMyRank()) + _interpolation_matrix->multiply(*_local_field->getField()); } void IntersectionDEC::sendData() { - _interpolation_matrix->multiply(*_source_field->getField()); + if (_source_group->containsMyRank()) + _interpolation_matrix->multiply(*_local_field->getField()); + else if (_target_group->containsMyRank()) + _interpolation_matrix->transposeMultiply(*_local_field->getField()); } + + } diff --git a/src/ParaMEDMEM/IntersectionDEC.hxx b/src/ParaMEDMEM/IntersectionDEC.hxx index ae9f6bc37..d4ddf3acd 100644 --- a/src/ParaMEDMEM/IntersectionDEC.hxx +++ b/src/ParaMEDMEM/IntersectionDEC.hxx @@ -1,7 +1,6 @@ #ifndef INTERSECTIONDEC_HXX_ #define INTERSECTIONDEC_HXX_ - namespace ParaMEDMEM { class DEC; @@ -39,6 +38,7 @@ namespace ParaMEDMEM string _method; InterpolationMatrix* _interpolation_matrix; + }; } diff --git a/src/ParaMEDMEM/MxN_Mapping.cxx b/src/ParaMEDMEM/MxN_Mapping.cxx index 8dc6fa749..729f1bd8c 100644 --- a/src/ParaMEDMEM/MxN_Mapping.cxx +++ b/src/ParaMEDMEM/MxN_Mapping.cxx @@ -203,6 +203,7 @@ void MxN_Mapping::sendRecv(double* sendfield, MEDMEM::FIELD& field) cons } //building the buffer of the elements to be sent vector offsets = _send_proc_offsets; + for (int i=0; i<_sending_ids.size();i++) { int iproc = _sending_ids[i].first; @@ -244,5 +245,91 @@ void MxN_Mapping::sendRecv(double* sendfield, MEDMEM::FIELD& field) cons } +/*! Exchanging field data between two groups of processes + * + * \param field MEDMEM field containing the values to be sent + * + * The ids that were defined by addElementFromSource method + * are sent. + */ +void MxN_Mapping::reverseSendRecv(double* recvfield, MEDMEM::FIELD& field) const +{ + CommInterface comm_interface=_union_group->getCommInterface(); + const MPIProcessorGroup* group = static_cast(_union_group); + + int nbcomp=field.getNumberOfComponents(); + double* sendbuf=0; + double* recvbuf=0; + if (_recv_ids.size() >0) + sendbuf = new double[_recv_ids.size()*nbcomp]; + if (_sending_ids.size()>0) + recvbuf = new double[_sending_ids.size()*nbcomp]; + + int* sendcounts = new int[_union_group->size()]; + int* senddispls=new int[_union_group->size()]; + int* recvcounts=new int[_union_group->size()]; + int* recvdispls=new int[_union_group->size()]; + + for (int i=0; i< _union_group->size(); i++) + { + sendcounts[i]=nbcomp*(_recv_proc_offsets[i+1]-_recv_proc_offsets[i]); + senddispls[i]=nbcomp*(_recv_proc_offsets[i]); + recvcounts[i]=nbcomp*(_send_proc_offsets[i+1]-_send_proc_offsets[i]); + recvdispls[i]=nbcomp*(_send_proc_offsets[i]); + } + //building the buffer of the elements to be sent + vector offsets = _recv_proc_offsets; + + for (int iproc=0; iproc<_union_group->size();iproc++) + for (int i=_recv_proc_offsets[iproc]; i<_recv_proc_offsets[iproc+1]; i++) + { + for (int icomp=0; icompgetComm(); + comm_interface.allToAllV(sendbuf, sendcounts, senddispls, MPI_DOUBLE, + recvbuf, recvcounts, recvdispls, MPI_DOUBLE, + *comm); + + + //setting zero values in the field + // for (int i=0; i< field.getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);i++) + // for (int j=0; jsize()]; i++) + { + for (int icomp=0; icomp& field); void sendRecv(double* field, MEDMEM::FIELD& field) const ; + void reverseSendRecv(double* field, MEDMEM::FIELD& field) const ; private : // ProcessorGroup& _local_group; -- 2.39.2