From: vbd Date: Mon, 16 Apr 2007 15:04:34 +0000 (+0000) Subject: adding the IntersectionDEC to the suite of tools X-Git-Tag: trio_trio_coupling~106 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=54d3dc86da1e70afb7f32120ac4158d3b753ce59;p=tools%2Fmedcoupling.git adding the IntersectionDEC to the suite of tools --- diff --git a/src/ParaMEDMEM/BlockTopology.cxx b/src/ParaMEDMEM/BlockTopology.cxx index 0469ed468..f95c12841 100644 --- a/src/ParaMEDMEM/BlockTopology.cxx +++ b/src/ParaMEDMEM/BlockTopology.cxx @@ -139,7 +139,7 @@ BlockTopology::BlockTopology(const ProcessorGroup& group, int nb_elem):_proc_gro } _cycle_type.resize(1); _cycle_type[0]=ParaMEDMEM::Block; - + delete[] nbelems_per_proc; } BlockTopology::~BlockTopology() diff --git a/src/ParaMEDMEM/CommInterface.hxx b/src/ParaMEDMEM/CommInterface.hxx index 1c4c3c822..5c9b799b9 100644 --- a/src/ParaMEDMEM/CommInterface.hxx +++ b/src/ParaMEDMEM/CommInterface.hxx @@ -43,7 +43,15 @@ public: MPI_Comm comm) const {return MPI_Allgather(sendbuf,sendcount, sendtype, recvbuf, recvcount, recvtype, comm); } - + int sendRecv(void* sendbuf, int sendcount, MPI_Datatype sendtype, + int dest, int sendtag, void* recvbuf, int recvcount, + MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm, + MPI_Status* status) + { + return + MPI_Sendrecv(sendbuf, sendcount, sendtype, dest, sendtag, recvbuf, + recvcount, recvtype, source, recvtag, comm,status); + } int worldSize() const {int size; MPI_Comm_size(MPI_COMM_WORLD, &size); return size;} }; diff --git a/src/ParaMEDMEM/DEC.cxx b/src/ParaMEDMEM/DEC.cxx index 27c93d0b0..b702ec21b 100644 --- a/src/ParaMEDMEM/DEC.cxx +++ b/src/ParaMEDMEM/DEC.cxx @@ -12,21 +12,35 @@ namespace ParaMEDMEM { + + DEC::DEC(ProcessorGroup& local_group, ProcessorGroup& distant_group):_source_field(0),_target_field(0), + _source_group(&local_group), _target_group(&distant_group) + { + + _union_group = local_group.fuse(distant_group); + } + + DEC::~DEC() + { + delete _union_group; + } + void DEC::attachTargetField(const ParaFIELD* field) { - _target_field=field; - //if (field!=0) - //{ - //BlockTopology* topo=dynamic_cast(field->getTopology()); - _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface()); - //} + _target_field=field; + //if (field!=0) + //{ + //BlockTopology* topo=dynamic_cast(field->getTopology()); + _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface()); + //} } void DEC::attachSourceField(const ParaFIELD* field) -{_source_field=field; -//if (field!=0) -//{ -// BlockTopology* topo=dynamic_cast(field->getTopology()); - _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface()); - //} +{ + _source_field=field; + //if (field!=0) + //{ + // BlockTopology* topo=dynamic_cast(field->getTopology()); + _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface()); + //} } } diff --git a/src/ParaMEDMEM/DEC.hxx b/src/ParaMEDMEM/DEC.hxx index eb7be98a9..5f6a8b8fe 100644 --- a/src/ParaMEDMEM/DEC.hxx +++ b/src/ParaMEDMEM/DEC.hxx @@ -10,6 +10,7 @@ class DEC { public: DEC():_source_field(0),_target_field(0){} + DEC(ProcessorGroup& local_group, ProcessorGroup& distant_group); void attachTargetField(const ParaFIELD* field); void attachSourceField(const ParaFIELD* field) ; virtual void prepareSourceDE()=0; @@ -17,13 +18,16 @@ public: virtual void recvData()=0; virtual void sendData()=0; virtual void synchronize()=0; - virtual ~DEC(){} + virtual ~DEC(); virtual void computeProcGroup(){}; protected: const ParaFIELD* _source_field; const ParaFIELD* _target_field; //! Processor group representing the union of target and source processors - ProcessorGroup* _group; + ProcessorGroup* _union_group; + ProcessorGroup* _source_group; + ProcessorGroup* _target_group; + const CommInterface* _comm_interface; }; diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx new file mode 100644 index 000000000..74e344b79 --- /dev/null +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -0,0 +1,376 @@ +#include +#include "CommInterface.hxx" +#include "ElementLocator.hxx" +#include "Topology.hxx" +#include "BlockTopology.hxx" +#include "ParaMESH.hxx" +#include "ProcessorGroup.hxx" +#include "MPIProcessorGroup.hxx" +#include "ParaSUPPORT.hxx" +#include "MEDMEM_Mesh.hxx" +#include "MEDMEM_Meshing.hxx" + +#include +namespace ParaMEDMEM +{ + +const double HUGE = 1e300; + +ElementLocator::ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group) +:_local_mesh(mesh.getMesh()), +_local_group(*mesh.getBlockTopology()->getProcGroup()), +_distant_group(distant_group) +{ + _union_group = _local_group.fuse(distant_group); + _computeBoundingBoxes(); +} + +ElementLocator::ElementLocator(const ParaSUPPORT& support, const ProcessorGroup& distant_group) +:_local_group(*support.getTopology()->getProcGroup()), +_distant_group(distant_group), +_union_group(_local_group.fuse(distant_group)) +{ + throw ("Element Locator SUPPORT constructor not implemented yet"); +} + +ElementLocator::~ElementLocator() +{ + delete _union_group; + delete [] _domain_bounding_boxes; +} + + +/*! Procedure for exchanging mesh between a distant proc and a local processor +\param idistantrank proc id on distant group +\param distant_mesh on return , points to a local reconstruction of the distant mesh +\param distant_ids on return, contains a vector defining a correspondence between the distant ids and the ids of the local reconstruction +*/ + +void ElementLocator::exchangeMesh(int idistantrank, MEDMEM::MESH*& distant_mesh, int*& distant_ids) +{ + int dim=_local_mesh->getSpaceDimension(); + int rank= _union_group->translateRank(&_distant_group,idistantrank); + if (find(_distant_proc_ids.begin(), _distant_proc_ids.end(),rank)==_distant_proc_ids.end()) + return; + + set elems; + double* distant_bb = _domain_bounding_boxes+rank*2*dim; + double* elem_bb=new double[2*dim]; + + //defining pointers to med + const int* conn=_local_mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE, + MED_EN::MED_NODAL, + MED_EN::MED_CELL, + MED_EN::MED_ALL_ELEMENTS); + const int* conn_index= _local_mesh->getConnectivityIndex( + MED_EN::MED_NODAL, + MED_EN::MED_CELL); + const double* coords = _local_mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE); + + for (int ielem=0; ielem<_local_mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);ielem++) + { + for (int i=0; ielem_bb[idim*2+1]? coords[node*dim+idim]: elem_bb[idim*2+1]; + } + } + if (_intersectsBoundingBox(elem_bb, distant_bb, dim)) + { + elems.insert(ielem+1); + } + } + + MEDMEM::MESH* send_mesh= _meshFromElems(elems); + + // Constituting an array containing the ids of the elements that are + // going to be sent to the distant subdomain. + // This array enables the correct redistribution of the data when the + // interpolated field is transmitted to the target array + + int* distant_ids_send = new int[elems.size()]; + + int index=0; + for (std::set::const_iterator iter = elems.begin(); iter!= elems.end(); iter++) + { + distant_ids_send[index]=*iter; + index++; + } + + _exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids); + delete[] distant_ids_send; + delete[] elem_bb; + delete send_mesh; +} + +void ElementLocator::_computeBoundingBoxes() +{ + CommInterface comm_interface=_union_group->getCommInterface(); + int dim = _local_mesh->getSpaceDimension(); + _domain_bounding_boxes = new double[2*dim*_union_group->size()]; + const double* coords = _local_mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE); + + int nbnodes = _local_mesh->getNumberOfNodes(); + double * minmax=new double [2*dim]; + for (int idim=0; idimcoords[i*dim+idim]?minmax[idim*2+1]:coords[i*dim+idim]); + } + MPIProcessorGroup* group=static_cast (_union_group); + const MPI_Comm* comm = group->getComm(); + comm_interface.allGather(minmax, 2*dim, MPI_DOUBLE, + _domain_bounding_boxes,2*dim, MPI_DOUBLE, + *comm); + + for (int i=0; i< _distant_group.size(); i++) + { + int rank= _union_group->translateRank(&_distant_group,i); + if (_intersectsBoundingBox(rank)) + _distant_proc_ids.push_back(rank); + } + delete[] minmax; +} + +bool ElementLocator::_intersectsBoundingBox(int irank) +{ + int dim=_local_mesh->getSpaceDimension(); + double* local_bb = _domain_bounding_boxes+_union_group->myRank()*2*dim; + double* distant_bb = _domain_bounding_boxes+irank*2*dim; + for (int idim=0; idim < _local_mesh->getSpaceDimension(); idim++) + { + bool intersects = distant_bb[idim*2]getCommInterface(); + + // First stage : exchanging sizes + int* send_buffer = new int[6]; + int* recv_buffer = new int[6]; + int nbtypes= local_mesh->getNumberOfTypes(MED_EN::MED_CELL); + int nbconn = local_mesh->getConnectivityLength(MED_EN::MED_FULL_INTERLACE, MED_EN::MED_NODAL,MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS); + int nbelems = local_mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + if (local_mesh !=0) + { + send_buffer[0] = local_mesh->getSpaceDimension(); + send_buffer[1] = local_mesh->getMeshDimension(); + send_buffer[2] = local_mesh->getNumberOfNodes(); + send_buffer[3] = nbelems; + send_buffer[4] = nbtypes; + send_buffer[5] = nbconn; + } + else + { + for (int i=0; i<6; i++) + send_buffer[i]=0; + } + MPIProcessorGroup* group=static_cast (_union_group); + const MPI_Comm* comm=(group->getComm()); + MPI_Status status; + // iproc_distant is the number of proc in distant group + // it must be converted to union numbering before communication + int iprocdistant_in_union = group->translateRank(&_distant_group, iproc_distant); + comm_interface.sendRecv(send_buffer, 6, MPI_INT, iprocdistant_in_union, 1112, recv_buffer, 6, MPI_INT,iprocdistant_in_union,1112, *comm, &status); + int distant_space_dim = recv_buffer[0]; + int distant_mesh_dim = recv_buffer[1]; + int distant_nnodes = recv_buffer[2]; + int distant_nb_elems = recv_buffer[3]; + int distant_nb_types = recv_buffer[4]; + int distant_nb_conn = recv_buffer[5]; + + delete[] send_buffer; + delete[] recv_buffer; + + //Second stage : exchanging connectivity buffers + int nb_integers = nbtypes*2+nbconn+1+nbelems; + send_buffer = new int[nb_integers]; + const MED_EN::medGeometryElement* types = local_mesh->getTypes(MED_EN::MED_CELL); + int* ptr_buffer=send_buffer; + const int* conn = local_mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE, MED_EN::MED_NODAL, + MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS); + const int* global_numbering = local_mesh->getGlobalNumberingIndex(MED_EN::MED_CELL); + + memcpy(ptr_buffer, types, nbtypes*sizeof(int)); + ptr_buffer+=nbtypes; + memcpy(ptr_buffer, global_numbering, (nbtypes+1)*sizeof(int)); + ptr_buffer+=nbtypes+1; + memcpy(ptr_buffer,conn, nbconn*sizeof(int)); + ptr_buffer+=nbconn; + memcpy(ptr_buffer, distant_ids_send, nbelems*sizeof(int)); + + //preparing the recv buffers + int nb_recv_integers = 2*distant_nb_types+1+distant_nb_conn+distant_nb_elems; + recv_buffer=new int[nb_recv_integers]; + + comm_interface.sendRecv(send_buffer, nb_integers, MPI_INT, iprocdistant_in_union, 1111, + recv_buffer, nb_recv_integers, MPI_INT, iprocdistant_in_union,1111, + *comm, &status); + + if (nb_integers>0) delete[] send_buffer; + //Third stage : exchanging coordinates + + int nb_recv_floats = distant_space_dim*distant_nnodes; + int nb_send_floats = local_mesh->getSpaceDimension()* local_mesh->getNumberOfNodes(); + double* recv_coords = new double[nb_recv_floats]; + double* coords = const_cast (local_mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE)); + comm_interface.sendRecv(coords, nb_send_floats, MPI_DOUBLE, iprocdistant_in_union, 1112, + recv_coords, nb_recv_floats, MPI_DOUBLE, iprocdistant_in_union, 1112, + *group->getComm(), &status); + + //Reconstructing an image of the distant mesh locally + + if (nb_recv_integers>0) + { + MEDMEM::MESHING* meshing = new MEDMEM::MESHING (); + int* recv_buffer_ptr = recv_buffer; + meshing->setCoordinates(distant_space_dim, distant_nnodes, recv_coords, "CARTESIAN", MED_EN::MED_FULL_INTERLACE); + meshing->setNumberOfTypes(distant_nb_types,MED_EN::MED_CELL); + + // converting the types from int to medGeometryElement + MED_EN::medGeometryElement* types_vector = new MED_EN::medGeometryElement[distant_nb_types]; + for (int i=0; isetTypes(types_vector, MED_EN::MED_CELL); + delete[] types_vector; + + recv_buffer_ptr+=distant_nb_types; + int* nbtypes = new int[distant_nb_types]; + for (int i=0; isetNumberOfElements( nbtypes, MED_EN::MED_CELL); + + for (int i=0; isetConnectivity(recv_buffer_ptr, MED_EN::MED_CELL,recv_buffer[i]); + recv_buffer_ptr+=nbtypes[i]*(recv_buffer[i]%100); + } + distant_ids_recv=new int [distant_nb_elems]; + for (int i=0; i0) + delete[] recv_coords; // the coordinates are present if the recv_buffer is not empty + +} + +MEDMEM::MESH* ElementLocator::_meshFromElems(set& elems) +{ + //defining pointers to med + const int* conn_mesh=_local_mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE, + MED_EN::MED_NODAL, + MED_EN::MED_CELL, + MED_EN::MED_ALL_ELEMENTS); + const int* conn_index= _local_mesh->getConnectivityIndex( + MED_EN::MED_NODAL, + MED_EN::MED_CELL); + const double* coords = _local_mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE); + set nodes; + int nbconn=0; + map nbelems_per_type; + for (set::const_iterator iter=elems.begin(); iter!=elems.end(); iter++) + { + for (int inode = conn_index[*iter-1]-1; inode < conn_index[*iter]-1; inode++) + nodes.insert(conn_mesh[inode]); + nbconn+=conn_index[*iter]-conn_index[*iter-1]; + MED_EN::medGeometryElement type = _local_mesh->getElementType(MED_EN::MED_CELL,*iter); + nbelems_per_type[type]++; + } + int* small2big=new int[nodes.size()]; + map big2small; + int i=0; + for (set::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++) + { + small2big[i]=*iter; + big2small[*iter]=i; + i++; + } + int* conn= new int[nbconn]; + double* new_coords = new double[nodes.size()*_local_mesh->getSpaceDimension()]; + int index=0; + for (set::const_iterator iter=elems.begin(); iter!=elems.end(); iter++) + { + for (int inode = conn_index[*iter-1]-1; inode < conn_index[*iter]-1; inode++) + { + conn[index]=big2small[conn_mesh[inode]]+1; + index++; + } + } + index=0; + for (set::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++) + { + int dim = _local_mesh->getSpaceDimension(); + for (int i=0; i::const_iterator iter= nbelems_per_type.begin(); iter!=nbelems_per_type.end(); iter++) + { + nbtypes[index]=iter->second; + new_types[index]=iter->first; + index++; + } + MEDMEM::MESHING* meshing = new MEDMEM::MESHING(); + meshing->setCoordinates(_local_mesh->getSpaceDimension(), nodes.size(), new_coords, string("CARTESIAN"), MED_EN::MED_FULL_INTERLACE); + meshing->setNumberOfTypes(nbelems_per_type.size(),MED_EN::MED_CELL); + meshing->setTypes(new_types,MED_EN::MED_CELL); + meshing->setNumberOfElements(nbtypes,MED_EN::MED_CELL); + + int* conn_ptr= conn; + for (int i=0; isetConnectivity(conn_ptr, MED_EN::MED_CELL,new_types[i]); + conn_ptr+=nbtypes[i]*(new_types[i]%100); + } + delete [] small2big; + delete [] nbtypes; + delete [] conn; + delete [] new_coords; + delete [] new_types; + return meshing; +} + +} diff --git a/src/ParaMEDMEM/ElementLocator.hxx b/src/ParaMEDMEM/ElementLocator.hxx new file mode 100644 index 000000000..af3243c9e --- /dev/null +++ b/src/ParaMEDMEM/ElementLocator.hxx @@ -0,0 +1,47 @@ +#ifndef ELEMENTLOCATOR_HXX_ +#define ELEMENTLOCATOR_HXX_ + +#include +#include + +namespace MEDMEM +{ + class MESH; +} +namespace ParaMEDMEM +{ +class ParaMESH; +class ProcessorGroup; +class ParaSUPPORT; +class InterpolationMatrix; + + +class ElementLocator +{ +public: + ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group); + ElementLocator(const ParaSUPPORT& support, const ProcessorGroup& distant_group); + virtual ~ElementLocator(); + void exchangeMesh(int idistantrank, MEDMEM::MESH*& distant_mesh, int*& distant_ids); + +private: + MEDMEM::MESH* _local_mesh; + std::vector _distant_meshes; + double* _domain_bounding_boxes; + const ProcessorGroup& _distant_group; + const ProcessorGroup& _local_group; + ProcessorGroup* _union_group; + std::vector _distant_proc_ids; + //InterpolationMatrix _matrix; + //MxN_Mapping _mapping; + + void _computeBoundingBoxes(); + bool _intersectsBoundingBox(int irank); + bool _intersectsBoundingBox(double* bb1, double* bb2, int dim); + void _exchangeMesh(MEDMEM::MESH* local_mesh, MEDMEM::MESH*& distant_mesh, int iproc_distant, const int* distant_ids_send, int*& distant_ids_recv); + MEDMEM::MESH* _meshFromElems(std::set& elems); +}; + +} + +#endif /*ELEMENTLOCATOR_HXX_*/ diff --git a/src/ParaMEDMEM/ExplicitMapping.hxx b/src/ParaMEDMEM/ExplicitMapping.hxx new file mode 100644 index 000000000..a4fae11d1 --- /dev/null +++ b/src/ParaMEDMEM/ExplicitMapping.hxx @@ -0,0 +1,155 @@ +#ifndef EXPLICIT_MAPPING_HXX_ +#define EXPLICIT_MAPPING_HXX_ + +#include +#include +#include + +namespace ParaMEDMEM +{ + class ExplicitMapping + { + public: + + ExplicitMapping():_numbers(0), _domains(0), _commbuffer(0) {}; + + ~ExplicitMapping() + { + if (_domains!=0) delete[] _domains; + if (_numbers!=0) delete[] _numbers; + if (_commbuffer!=0) delete[] _commbuffer; + } + void pushBackElem(pair idistant) { + _mapping.push_back(idistant); + } + void setDistantElem(int ilocal, pair idistant) + { + _mapping[ilocal]=idistant; + } + + int nbDistantDomains() + { + if (_distant_domains.empty()) + { + for (vector >::const_iterator iter= _mapping.begin(); + iter!=_mapping.end(); + iter++) + _distant_domains.insert(iter->first); + } + return _distant_domains.size(); + } + + pair getDistantNumbering(int ielem)const + { + return _mapping[ielem]; + } + + + int getDistantDomain(int i) + { + if (_domains==0) + computeNumbers(); + + return _domains[i]; + } + + int getNbDistantElems(int i) + { + if (_numbers==0) + computeNumbers(); + return _numbers[i]; + } + + int* serialize(int idproc) + { + _commbuffer=new int[_mapping.size()*2]; + vector offsets(_distant_domains.size()); + offsets[0]=0; + for (int i=1; i<_distant_domains.size();i++) + offsets[i]=offsets[i-1]+_numbers[i-1]; + + for (int i=0; i< _mapping.size(); i++) + { + int offset= offsets[_mapping[i].first]; + _commbuffer[offset*2]=idproc; + _commbuffer[offset*2+1]=_mapping[i].second; + offsets[_mapping[i].first]++; + } + return _commbuffer; + } + + void unserialize(int nbprocs, int* sizes,int nbtarget, int* targetrank, int* commbuffer) + { + int total_size=0; + for (int i=0; i< nbprocs; i++) + total_size+=sizes[i]; + + _mapping.resize(total_size); + _buffer_index=new int[total_size]; + int indmap=0; + for (int i=0; i0) + { + _numbers[index]=sizes[targetrank[i]]; + _domains[index]=i; + index++; + } + } + _send_counts=new int[nbprocs]; + for (int i=0; i > _mapping; + set _distant_domains; + int* _numbers; + int* _domains; + int* _commbuffer; + int* _buffer_index; + int* _send_counts; + + void computeNumbers() + { + map counts; + if (_numbers==0) + { + _numbers=new int[nbDistantDomains()]; + _domains=new int[nbDistantDomains()]; + for (int i=0; i< _mapping.size(); i++) + { + if ( counts.find(_mapping[i].first) == counts.end()) + counts.insert(make_pair(_mapping[i].first,1)); + else + (counts[_mapping[i].first])++; + } + int counter=0; + for (map::const_iterator iter=counts.begin(); + iter!=counts.end(); + iter++) + { + _numbers[counter]=iter->second; + _domains[counter]=iter->first; + counter++; + } + } + } + }; + +} +#endif diff --git a/src/ParaMEDMEM/ExplicitTopology.cxx b/src/ParaMEDMEM/ExplicitTopology.cxx index 3ba799f73..c5be133df 100644 --- a/src/ParaMEDMEM/ExplicitTopology.cxx +++ b/src/ParaMEDMEM/ExplicitTopology.cxx @@ -52,8 +52,8 @@ ExplicitTopology::ExplicitTopology(const ExplicitTopology& topo, int nb_componen _proc_group = topo._proc_group; _nb_elems = topo._nb_elems; _nb_components = nb_components; - _loc2glob=new int[2*_nb_elems]; - for (int i=0; i<2*_nb_elems; i++) + _loc2glob=new int[_nb_elems]; + for (int i=0; i<_nb_elems; i++) { _loc2glob[i]=topo._loc2glob[i]; } @@ -63,6 +63,7 @@ ExplicitTopology::ExplicitTopology(const ExplicitTopology& topo, int nb_componen ExplicitTopology::~ExplicitTopology() { + if (_loc2glob != 0) delete[] _loc2glob; } diff --git a/src/ParaMEDMEM/INTERPOLATION_2D.cxx b/src/ParaMEDMEM/INTERPOLATION_2D.cxx new file mode 100755 index 000000000..27eadc667 --- /dev/null +++ b/src/ParaMEDMEM/INTERPOLATION_2D.cxx @@ -0,0 +1,879 @@ +#include "INTERPOLATION_2D.hxx" +using namespace std; +using namespace MED_EN; +using namespace MEDMEM; + + + +struct monMaillageS +{ + int _NbMaille; + int _NbNoeud; + const int* _Connect; + const double* _Coord; + const int* _ReversNConnect; + const int* _ReversNConnectIndex; + double _DimCaracteristic; +}; + +struct monMaillageP +{ + int _NbMaille; + int _NbNoeud; + const int* _Connect; + const double* _Coord; + double _DimCaracteristic; +}; + + + +INTERPOLATION_2D::INTERPOLATION_2D() +{ + _Precision=1.0E-12; + _DimCaracteristic=1; + _SecondMethodOfLocalisation=1; + _PrintLevel=0; +} + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* Options : */ +/* Precision : for geometric computation */ +/* SecondMethodOfLocalisation : 0/1 */ +/* PrintLevel : between 0 and 3 */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +void INTERPOLATION_2D::setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel) +{ + _Precision=Precision; + _SecondMethodOfLocalisation=SecondMethodOfLocalisation; + _PrintLevel=PrintLevel; +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* calcul la surface d'un triangle */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + +double INTERPOLATION_2D::Surf_Tri(const double* P_1,const double* P_2,const double* P_3) +{ + double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]); + double Surface = 1/2.0*abs(A); + return Surface; +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* fonction qui calcul le déterminant */ +/* de deux vecteur(cf doc CGAL). */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ + +//fonction qui calcul le déterminant des vecteurs: P3P1 et P3P2 +//(cf doc CGAL). + +double INTERPOLATION_2D::mon_determinant(const double* P_1,const double* P_2, + const double* P_3) +{ + double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]); + return mon_det; +} + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ +//calcul la norme du vecteur P1P2 + +double INTERPOLATION_2D::norme_vecteur(const double* P_1,const double* P_2) +{ + double X=P_1[0]-P_2[0]; + double Y=P_1[1]-P_2[1]; + double norme=sqrt(X*X+Y*Y); + return norme; +} + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* calcul le cos et le sin de l'angle P1P2,P1P3 */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + +vector INTERPOLATION_2D::calcul_cos_et_sin(const double* P_1,const double* P_2, + const double* P_3) +{ + + vector Vect; + double P1_P2=norme_vecteur(P_1,P_2); + double P2_P3=norme_vecteur(P_2,P_3); + double P3_P1=norme_vecteur(P_3,P_1); + + double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3; + double D=2.0*P1_P2*P3_P1; + double COS=N/D; + Vect.push_back(COS); + double V=mon_determinant(P_2,P_3,P_1); + double D_1=P1_P2*P3_P1; + double SIN=V/D_1; + Vect.push_back(SIN); + + return Vect; + +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/*fonction pour vérifier qu'un point n'a pas déja été considérer dans */ +/* le vecteur et le rajouter au vecteur sinon. */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + +void INTERPOLATION_2D::verif_point_dans_vect(const double* P, vector& V) +{ + + int taille=V.size(); + // bool isPresent=false; + for(int i=0;i0.0?dx:-dx; + if (dx>_Precision) + continue; + double dy=P[1]-V[2*i+1]; + dy=dy>0.0?dy:-dy; + if (dy<_Precision) + + return; + // if( sqrt((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))/_DimCaracteristic<_Precision) + // isPresent=true; + + } + // if(!isPresent) + //{ + + V.push_back(P[0]); + V.push_back(P[1]); + + //} +} + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* calcul les coordonnées du barycentre d'un polygone */ +/* le vecteur en entrée est constitué des coordonnées */ +/* des sommets du polygone */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + +vector INTERPOLATION_2D::bary_poly(const vector& V) +{ + vector Bary; + int taille=V.size(); + double x=0; + double y=0; + + for (int i=0;i& Poly) +{ + + double Surface=0; + for (int i=0; i<(Poly.size())/2-2; i++) + { + double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] ); + Surface=Surface + Surf ; + } + return Surface ; +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ +/* calcul de l'intersection de deux segments: segments P1P2 avec P3P4 */ +/* . Si l'intersection est non nulle et si celle-ci n'est */ +/* n'est pas déjà contenue dans Vect on la rajoute à Vect */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ + +void INTERPOLATION_2D::inters_de_segment(const double * P_1,const double * P_2, + const double * P_3,const double * P_4,vector& Vect) +{ + + + // calcul du déterminant de P1_1P_2 et P_3P_4. + double det=(P_2[0]-P_1[0])*(P_4[1]-P_3[1])-(P_4[0]-P_3[0])*(P_2[1]-P_1[1]); + + + if(abs(det)/_DimCaracteristic>_Precision) + { + + double k_1=1/((P_2[0]-P_1[0])*(P_3[1]-P_4[1])+(P_4[0]-P_3[0])*(P_2[1]-P_1[1])) + *((P_3[1]-P_4[1])*(P_3[0]-P_1[0])+(P_4[0]-P_3[0])*(P_3[1]-P_1[1])); + + if( k_1>=0 && k_1<=1) + { + double k_2=1/((P_4[0]-P_3[0])*(P_1[1]-P_2[1])+(P_2[0]-P_1[0])*(P_4[1]-P_3[1])) + *((P_1[1]-P_2[1])*(P_1[0]-P_3[0])+(P_2[0]-P_1[0])*(P_1[1]-P_3[1])); + + + if( k_2>=0 && k_2<=1) + { + double P_0[2]; + P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]); + P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]); + verif_point_dans_vect(P_0,Vect); + + } + } + } +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* fonction qui teste si un point est dans une maille */ +/* point: P_0 */ +/* P_1, P_2, P_3 sommet des mailles */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + +int INTERPOLATION_2D::point_dans_triangle(const double* P_0,const double* P_1, + const double* P_2,const double* P_3) +{ + + int A=0; + double det_1=mon_determinant(P_1,P_3,P_0); + double det_2=mon_determinant(P_3,P_2,P_0); + double det_3=mon_determinant(P_2,P_1,P_0); + if( (det_1>=0 && det_2>=0 && det_3>=0) || (det_1<=0 && det_2<=0 && det_3<=0) ) + { + A=1; + } + + return A; +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* fonction qui rajoute les sommet du triangle P dans le vecteur V */ +/* si ceux-ci sont compris dans le triangle S et ne sont pas déjà dans */ +/* V. */ +/*sommets de P: P_1, P_2, P_3 */ +/*sommets de S: P_4, P_5, P_6 */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + +void INTERPOLATION_2D::rajou_sommet_triangl(const double* P_1,const double* P_2,const double* P_3, + const double* P_4,const double* P_5,const double* P_6,vector& V) +{ + + int A_1=point_dans_triangle(P_1,P_4,P_5,P_6); + if(A_1==1) + {verif_point_dans_vect(P_1,V); + // if (V.size()==1) +// { +// // all nodes are necessarily in the triangle +// verif_point_dans_vect(P_2,V); +// verif_point_dans_vect(P_3,V); +// return; +// } + } + int A_2=point_dans_triangle(P_2,P_4,P_5,P_6); + if(A_2==1) + {verif_point_dans_vect(P_2,V);} + + int A_3=point_dans_triangle(P_3,P_4,P_5,P_6); + if(A_3==1) + {verif_point_dans_vect(P_3,V);} + + + +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ +/* calcul l'intersection de deux triangles */ +/* P_1, P_2, P_3: sommets du premier triangle */ +/* P_4, P_5, P_6: sommets du deuxième triangle */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ + +vector INTERPOLATION_2D::intersec_de_triangle(const double* P_1,const double* P_2, + const double* P_3,const double* P_4,const double* P_5,const double* P_6) +{ + + //debug cout << " T1 : " << *P_1 << " , " << *(P_1+1) << " , " << *P_2 << " , " << *(P_2+1) << " , " << *P_3 << " , " << *(P_3+1)<< endl; + //debug cout << " T2 : " << *P_4 << " , " << *(P_4+1) << " , " << *P_5 << " , " << *(P_5+1) << " , " << *P_6 << " , " << *(P_6+1)<< endl; + vector Vect; + + inters_de_segment(P_1,P_2,P_4,P_5,Vect); + inters_de_segment(P_1,P_2,P_5,P_6,Vect); + inters_de_segment(P_1,P_2,P_6,P_4,Vect); + inters_de_segment(P_2,P_3,P_4,P_5,Vect); + inters_de_segment(P_2,P_3,P_5,P_6,Vect); + inters_de_segment(P_2,P_3,P_6,P_4,Vect); + inters_de_segment(P_3,P_1,P_4,P_5,Vect); + inters_de_segment(P_3,P_1,P_5,P_6,Vect); + inters_de_segment(P_3,P_1,P_6,P_4,Vect); + rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect); + rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect); + //debug cout << " Inter : "; + //debug for (int i=0; i INTERPOLATION_2D::reconstruct_polygon(const vector& V) +{ + + int taille=V.size(); + if(taille<=6) + {return V;} + else + { + double COS[taille/2]; + double SIN[taille/2]; + double angle[taille/2]; + vector Bary=bary_poly(V); + COS[0]=1.0; + SIN[0]=0.0; + angle[0]=0.0; + for(int i=0; i Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]); + COS[i+1]=Trigo[0]; + SIN[i+1]=Trigo[1]; + if(SIN[i+1]>=0) + {angle[i+1]=acos(COS[i+1]);} + else + {angle[i+1]=-acos(COS[i+1]);} + } + + //ensuite on ordonne les angles. + vector Pt_ordonne; + Pt_ordonne.reserve(taille); + map Ordre; + for(int i=0;i::iterator mi; + for(mi=Ordre.begin();mi!=Ordre.end();mi++) + { + int j=(*mi).second; + Pt_ordonne.push_back(V[2*j]); + Pt_ordonne.push_back(V[2*j+1]); + } + return Pt_ordonne; + } +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */ +/*fonction pour trouver les mailles voisines d'une maille triangle.*/ +/* (mailles voisines par les faces) */ +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */ + +vector INTERPOLATION_2D::touv_maill_vois(int i_S,const monMaillageS& MaS) +{ + + int N_1=MaS._Connect[3*(i_S-1)]; + int N_2=MaS._Connect[3*(i_S-1)+1]; + int N_3=MaS._Connect[3*(i_S-1)+2]; + vector tab; + int I_1=MaS._ReversNConnectIndex[N_1-1]-1;// Dans MED le n°des noeuds commencent à 1. + int I_2=MaS._ReversNConnectIndex[N_2-1]-1; + int I_3=MaS._ReversNConnectIndex[N_3-1]-1; + int F_1=MaS._ReversNConnectIndex[N_1]-1; + int F_2=MaS._ReversNConnectIndex[N_2]-1; + int F_3=MaS._ReversNConnectIndex[N_3]-1; + + + vector V_12; + V_12.reserve(2); + insert_iterator > ib_1(V_12,V_12.begin()); + vector V_23; + V_23.reserve(2); + insert_iterator > ib_2(V_23,V_23.begin()); + vector V_13; + V_13.reserve(2); + insert_iterator > ib_3(V_13,V_13.begin()); + + + set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_2, + MaS._ReversNConnect+F_2,ib_1); + set_intersection(MaS._ReversNConnect+I_2,MaS._ReversNConnect+F_2,MaS._ReversNConnect+I_3, + MaS._ReversNConnect+F_3,ib_2); + set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_3, + MaS._ReversNConnect+F_3,ib_3); + + + for(int i=0;i<(V_12).size();i++) + { + if(V_12[i]!=i_S) + { + tab.push_back(V_12[i]); + break; + } + } + + for(int i=0;i& V) +{ + int taille=V.size(); + int A=0; + for(int i=0;i& localis ) +{ + + //calcul des coordonnées des barycentre des mailles de P + + + //calcule des coordonnees du barycentre d'une cellule + std::auto_ptr Support ( new SUPPORT(&myMesh_P,"monSupport",MED_CELL) ); + std::auto_ptr > Barycentre (myMesh_P.getBarycenter(Support.get() )); + const double* Bary=Barycentre->getValue(); + + + //localisation des barycentres des mailles de P dans les mailles de S. + Meta_Wrapper<2> fromWrapper(MaS._NbNoeud,(double *)MaS._Coord,(MEDMEM::CONNECTIVITY*)myMesh_S.getConnectivityptr()); + Meta_Wrapper<2> toWrapper(MaP._NbMaille,(double *)Bary); + Meta_Mapping<2> mapping(&fromWrapper,&toWrapper); + mapping.Cree_Mapping(0); + for(int i=0;i monInterpol_S(myMesh_S); + for(int i_P=0;i_P Inter=intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1), + + MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1)); + + if(Inter.size()>0) + { + //debug cout << "Localise : maille proche barycentre trouvée : " << N_Maille << endl; + Test=1; + localis[i_P]=N_Maille; + break; + } + } + if(Test==0) + { + int NoeuMaillP[3]={NP_1,NP_2,NP_3}; + for(int Num=0;Num<3;Num++) + { + int Num_noeud=NoeuMaillP[Num]; + double coord_i_P_0[2]; + //VB + coord_i_P_0[0]= MaP._Coord[2*(Num_noeud-1)]; + coord_i_P_0[1]=MaP._Coord[2*(Num_noeud-1)+1]; + //VB + int Vois_B=monInterpol_S.getNearestNode(coord_i_P_0)+1; + int Ni=MaS._ReversNConnectIndex[Vois_B-1]; + int Nf=MaS._ReversNConnectIndex[Vois_B]; + int diff=Nf-Ni; + + for(int j=0;j Inter=intersec_de_triangle(MaS._Coord+2*(N_1-1),MaS._Coord+2*(N_2-1), + MaS._Coord+2*(N_3-1),MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1) ); + + if(Inter.size()>0) + { + //debug cout << "Localise : maille proche sommet " << Num+1 << " trouvée : " << N_Maille << endl; + Test=1; + localis[i_P]=N_Maille; + break;//on sort de la boucle sur les maille commune à un noeud de i_S + } + } + if(Test==1) + {break;}//on sort de la boucle sur les noeuds de i_P + } + } + } + } + } +} + + + + + + + + + +/* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +/* fonction qui permet de remplir le vecteur donnant la surface d'intersection */ +/* de la maille i_P du maillage P (projeté) avec la maille i_S du maillage S (support) */ +/* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + + +inline void INTERPOLATION_2D::rempli_vect_d_intersect(int i_P,int i_S1,const monMaillageP& MaP,const monMaillageS& MaS, + vector >& Surface_d_intersect,FIELD& myField_P) +{ + + + //on récupere le n° des noeuds. + + //debug cout << "\nintersect maille " << i_P << endl; + int NP_1=MaP._Connect[3*(i_P-1)]; + int NP_2=MaP._Connect[3*(i_P-1)+1]; + int NP_3=MaP._Connect[3*(i_P-1)+2]; + + + //on calcule la surface de la maille i_P + + double Surf_i_P =Surf_Tri(MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1), + MaP._Coord+2*(NP_3-1)); + + double Surface = 0 ; + vector maill_deja_vu_S ; + + maill_deja_vu_S.push_back(i_S1); + + for (int M_S=0; M_S_Precision && M_S!=maill_deja_vu_S.size() ) + { + int i_S=maill_deja_vu_S[M_S]; + + int NS_1=MaS._Connect[3*(i_S-1)]; + int NS_2=MaS._Connect[3*(i_S-1)+1]; + int NS_3=MaS._Connect[3*(i_S-1)+2]; + + + vector Inter=intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1), + MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1)); + + + + + //on teste l'intersection. + int taille=Inter.size()/2; + //debug cout << " -> maille source : " << i_S << " , nb intersection=" << taille << endl; + + /* CAS1 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + if(taille==0) + {int Rien=0;} + + /* CAS 2 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + else if(taille==1 || taille==2) + { + //on ne remplit pas le vecteur des correspondances n° maille-surfaces + //d'intersections mais on récupère le numéro des mailles voisines de la maille i_S. + vector M_Vois=touv_maill_vois(i_S,MaS); + for(int i=0;i< M_Vois.size();i++) + {verif_maill_dans_vect(M_Vois[i],maill_deja_vu_S);} + } + + /* CAS 3 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ + else if(taille==3) + { + double Surf_inter = Surf_Poly(Inter) ; + + //on remplit le vecteur des correspondances n° maille-surfaces d'intersections. + Surface_d_intersect[i_P-1][i_S]=Surf_inter; + + // on récupère le numéro des mailles voisines de la maille i_S. + vector M_Vois=touv_maill_vois(i_S,MaS); + for(int i_M=0;i_M<(M_Vois).size();i_M++) + {verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);} + + Surface = Surface + Surf_inter ; + } + + /* CAS 4 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + + else //taille>3 + { + vector Poly_Inter=reconstruct_polygon(Inter); + + double Surf_inter =Surf_Poly(Poly_Inter) ; + + //on remplit le vecteur des correspondances n° maille-surfaces d'intersection + (Surface_d_intersect[i_P-1])[i_S]=Surf_inter; + + // on récupère le numéro des mailles voisines de la maille i_S. + vector M_Vois=touv_maill_vois(i_S,MaS); + for(int i_M=0;i_M<(M_Vois).size();i_M++) + {verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);} + + Surface = Surface + Surf_inter ; + } + } + //debug cout << " surface = " << Surface << endl << flush; + } + + /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + //on rentre la fraction totale de la maille i_P qui a été considérer lors de l'interpolation. + double Fract=Surface/Surf_i_P; + myField_P.setValueIJ(i_P,1,Fract); + +} + + +/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ +//fonction principale pour interpoler deux maillages triangulaires. +vector > INTERPOLATION_2D::interpol_maillages(const MEDMEM::MESH& myMesh_S, + const MEDMEM::MESH& myMesh_P) +{ + + + /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + // on vérifie que les deux maillages sont formés de triangles. + int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_CELL); + int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_CELL); + if ( NumberOfCellsTypes_S != 1 || NumberOfCellsTypes_P != 1) + { cout<<"l'un au moins des deux maillages n'est pas triangulaires"< > BoxS=myMesh_S.getBoundingBox(); + vector > BoxP=myMesh_P.getBoundingBox(); + double AreaS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0])+(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1])); + MaS._DimCaracteristic=sqrt(4./sqrt(3.)*AreaS/MaS._NbMaille); + double AreaP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0])+(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1])); + MaP._DimCaracteristic=sqrt(4./sqrt(3.)*AreaP/MaP._NbMaille); + _DimCaracteristic=min(MaS._DimCaracteristic,MaP._DimCaracteristic); + if (_PrintLevel) + { + cout<<"recherche de la dimension caractéristique des maillages 2D :"< localis; + localis.reserve(MaP._NbMaille); + MEDMEM::MESH* ptr_S = const_cast(&myMesh_S); + MEDMEM::MESH* ptr_P = const_cast(&myMesh_P); + + cout << "INTERPOLATION_2D::local_iP_dansS"< MAP_init; + vector > Surface_d_intersect(MaP._NbMaille,MAP_init);//on initialise. + + /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + //création d'un champ pour rentrer les fractions de mailles considérées. + + std::auto_ptr mySupport_P (new SUPPORT(const_cast(&myMesh_P),"Support on all CELLS",MED_CELL) ); + MEDMEM::FIELD myField_P(mySupport_P.get(),1); + string NameF=" M2 dans M1" ; + myField_P.setName(NameF); + string NameC="fracSurf"; + string ComponentsNames[1]={NameC}; + myField_P.setComponentsNames(ComponentsNames); + myField_P.setComponentsDescriptions(ComponentsNames); + string ComponentsUnits[1]={"%"}; + myField_P.setMEDComponentsUnits(ComponentsUnits); + + /*à remplacer par: + WriteFractSurfInFile(vector< > >& FractVol,MEDMEM::MESH* myMesh1, + MEDMEM::MESH* myMesh2,string FileName,string NameMeshP,string NameMeshS) + */ + + /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ + //boucle sur les mailles de P. + //Coeur de l'algorithme + + for(int i_P=0; i_P0) + rempli_vect_d_intersect(i_P+1,i_S_init,MaP,MaS,Surface_d_intersect,myField_P); + + } + + if (_PrintLevel >= 2) + { + cout<::iterator surface; + for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++) + cout<<" ("<=3) +// { +//// string NameFile="fractSurf.med"; +//// int id1=myMesh_P.addDriver(MED_DRIVER,NameFile,"M1"); +//// int id2=myMesh_S.addDriver(MED_DRIVER,NameFile,"M2"); +//// int id3=myField_P.addDriver(MED_DRIVER,NameFile,NameF); +//// myMesh_P.write(id1); +//// myMesh_S.write(id2); +//// myField_P.write(id3); +//// +//// cout<& FractSurf,MEDMEM::MESH* myMesh_P, + MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S) +{ + //On sauvegarde le champ crée sur la maillage P. + string FileName="FracSurf"+NameMesh_S+"_to_"+NameMesh_P; + int id1=(*myMesh_P).addDriver(MED_DRIVER,FileName,"MP"); + int id2=(*myMesh_S).addDriver(MED_DRIVER,FileName,"MS"); + int id3=FractSurf.addDriver(MED_DRIVER,FileName,"% MP"); + (*myMesh_P).write(id1); + (*myMesh_S).write(id2); + FractSurf.write(id3); +} diff --git a/src/ParaMEDMEM/INTERPOLATION_2D.hxx b/src/ParaMEDMEM/INTERPOLATION_2D.hxx new file mode 100755 index 000000000..58026ff48 --- /dev/null +++ b/src/ParaMEDMEM/INTERPOLATION_2D.hxx @@ -0,0 +1,150 @@ +#ifndef _INTERPOLATION_2D_HXX_ +#define _INTERPOLATION_2D_HXX_ + + +#include "MEDMEM_InterpolationHighLevelObjects.hxx" +#include "MEDMEM_Mesh.hxx" +#include "MEDMEM_Field.hxx" +#include "MEDMEM_Support.hxx" +#include "MEDMEM_Family.hxx" +#include "MEDMEM_Group.hxx" +#include "MEDMEM_Coordinate.hxx" +#include "MEDMEM_Connectivity.hxx" +#include "MEDMEM_CellModel.hxx" +#include "MEDMEM_DriverFactory.hxx" +#include "MEDMEM_Meshing.hxx" +#include "MEDMEM_define.hxx" +#include "MEDMEM_Interpolation.hxx" + + +#include +#include +#include +#include +#include +#include +#include +#include + + +struct monMaillageS; + + +struct monMaillageP; + + +class INTERPOLATION_2D +{ + + + private: + double _Precision; + double _DimCaracteristic; + int _SecondMethodOfLocalisation; + int _PrintLevel; + // Méthodes publiques + public: + + INTERPOLATION_2D(); + + + // precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions + void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel); + + //pour calculer la surface d'un triangle. + double Surf_Tri(const double* P_1,const double* P_2, + const double* P_3); + + //calcul déterminant de trois points (cf doc CGAL) + double mon_determinant(const double* P_1,const double* P_2, + const double* P_3); + + //calcul la norme d'un vecteur. + double norme_vecteur(const double* P_1,const double* P_2); + + //calcul le cos et le sin. + vector calcul_cos_et_sin(const double* P_1, + const double* P_2,const double* P_3); + + //fonction pour vérifier qu'un point n'a pas déja été considérer dans + //le vecteur et le rajouter au vecteur sinon. + void verif_point_dans_vect(const double* P, vector& V); + + //calcul les coordonnées du barycentre + //d'un polygone + vector bary_poly(const vector& V); + + //fonction qui calcul la surface d'un polygone + double Surf_Poly(const vector& Poly); + + //calcul de l'intersection de deux segments. + void inters_de_segment(const double* P_1,const double* P_2, + const double* P_3,const double* P_4,vector& Vect); + + //fonction qui teste si un point est dans une maille + int point_dans_triangle(const double* P_0,const double* P_1, + const double* P_2,const double* P_3); + + //fonction qui rajouter les sommets du triangle P2P2P3 si ceux-ci sont compris + // dans le triangle P4P5P6 et n'ont pas encore été considérer. + //Rq:point du triangle donné dans ordre quelconque. + void rajou_sommet_triangl(const double* P_1,const double* P_2, + const double* P_3,const double* P_4,const double* P_5, + const double* P_6,vector& V); + + //calcul l'intersection de deux triangles (sommets du triangle donné dans + //ordre quelconque) + vector intersec_de_triangle(const double* P_1,const double* P_2, + const double* P_3,const double * P_4,const double * P_5, + const double * P_6); + + //fonction pour reconstituer un polygon convexe à partir + //d'un nuage de point. + vector reconstruct_polygon(const vector& V); + + //fonction pour trouver les mailles voisines d'une maille triangle. + vector touv_maill_vois(int i_S,const monMaillageS& MaS); + + //fonction pour vérifier qu'un n°de maille n'a pas déja été considérer + // dans le vecteur et le rajouter au vecteur sinon. + void INTERPOLATION_2D::verif_maill_dans_vect(int Num, vector& V); + + + //localise le barycentre des mailles de P dans le maillage S + void local_iP_dans_S(MEDMEM::MESH& myMesh_S,MEDMEM::MESH& myMesh_P, + monMaillageP& MaP,monMaillageS& MaS,vector& localis ); + + + //fonction qui permet de remplir le vecteur donnant la surface d'intersection + //de la mailles i_P du maillage projetté avec la maille i_S du maillage support. + inline void rempli_vect_d_intersect(int i_P,int i_S,const monMaillageP& MaP,const monMaillageS& MaS, + vector >& Surface_d_intersect, + FIELD& myField_P); + + + //fonction principale pour interpoler deux maillages triangulaires. + std::vector > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2); + + + // fonction qui écrit les résultats d'annelise dans un fichier: + // pour chaque maille du maillage 1 on stoque la fraction volumique + // considéré lors de l'algorithme + + + + + + void WriteFractSurfInFile(MEDMEM::FIELD& FractSurf,MEDMEM::MESH* myMesh_P, + MEDMEM::MESH* myMesh_S,string FileName,string NameMesh_P,string NameMesh_S); + + + + // MEDMEM::FIELD* createField(); + // ... + + private: + + +}; + +#endif diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx new file mode 100644 index 000000000..508760181 --- /dev/null +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -0,0 +1,177 @@ +#include "ParaMESH.hxx" +#include "ProcessorGroup.hxx" +#include "MxN_Mapping.hxx" +#include "InterpolationMatrix.hxx" +#include "INTERPOLATION_2D.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. +The matrix is built and stored on the processors belonging to the source +group. +*/ + +namespace ParaMEDMEM +{ + + /*! + Creates an empty matrix structure linking two distributed supports. + The method must be called by all processors belonging to source and target groups. + \param source_support local support + \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 string& method): +_source_support(*source_support.getMesh()), +_local_group(local_group), +_distant_group(distant_group), +_mapping(local_group, distant_group), +_method(method) +{ + int nbelems= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + _row_offsets.resize(nbelems+1,0); + + +} + +InterpolationMatrix::~InterpolationMatrix() +{ +} + +/*! +adds the contribution of a distant subdomain to the interpolation matrix +\param distant_support local representation of the distant subdomain +\param iproc_distant id of the distant subdomain (in the distant group) +\param distant_elems mapping between the local representation of the subdomain and the actual elem ids on the distant subdomain + */ +void InterpolationMatrix::addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems) +{ + INTERPOLATION_2D interpolator; + int nbelems= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + //_row_offsets.resize(nbelems); + vector > surfaces = interpolator.interpol_maillages(distant_support,_source_support); + if (surfaces.size() != nbelems) + throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix"); + + MEDMEM::SUPPORT support(&distant_support,"all cells", MED_EN::MED_CELL); + MEDMEM::FIELD* triangle_surf = distant_support.getArea(&support); + + 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++) + { + double surf = triangle_surf->getValueIJ(iter->first,1); + surf = (surf>0)?surf:-surf; + vector >::iterator iter2 = find (_col_offsets.begin(), _col_offsets.end(),make_pair(iproc_distant,iter->first)); + int col_id; + if (iter2==_col_offsets.end()) + { + _col_offsets.push_back(make_pair(iproc_distant,iter->first)); + col_id =_col_offsets.size(); + _mapping.addElementFromSource(col_id,iproc_distant,distant_elems[iter->first-1]); + } + else + { + col_id=iter2-_col_offsets.begin()+1; + } + _col_numbers.push_back(col_id); + _coeffs.push_back(iter->second/surf); + + } + } + delete 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_indices.size(); + int nbrows= _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + for (int irow=0; irow0) + delete[] target_value; +} + + +} diff --git a/src/ParaMEDMEM/InterpolationMatrix.hxx b/src/ParaMEDMEM/InterpolationMatrix.hxx new file mode 100644 index 000000000..ac686308c --- /dev/null +++ b/src/ParaMEDMEM/InterpolationMatrix.hxx @@ -0,0 +1,43 @@ +#ifndef INTERPOLATIONMATRIX_HXX_ +#define INTERPOLATIONMATRIX_HXX_ + +#include "MEDMEM_Field.hxx" +#include "MxN_Mapping.hxx" + +namespace ParaMEDMEM +{ + class InterpolationMatrix + { + public: + //InterpolationMatrix(); + InterpolationMatrix(const ParaMEDMEM::ParaMESH& source_support, + const ProcessorGroup& local_group, + const ProcessorGroup& distant_group, + const string& method); + + //InterpolationMatrix(const MEDMEM::MESH& source_support, const string& method); + //InterpolationMatrix(const MEDMEM::MESH& target_support, const MEDMEM::MESH& source_support); + virtual ~InterpolationMatrix(); + void addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems); + void multiply(MEDMEM::FIELD&) const; + void prepare(); + int getNbRows() const {return _row_offsets.size();} + + private: + // vector > _source_indices; + vector _row_offsets; + vector _col_numbers; + vector > _col_offsets; + vector _coeffs; + const MEDMEM::MESH& _source_support; + MxN_Mapping _mapping; + string _method; + const ProcessorGroup& _local_group; + const ProcessorGroup& _distant_group; + + + }; + +} + +#endif /*INTERPOLATIONMATRIX_HXX_*/ diff --git a/src/ParaMEDMEM/IntersectionDEC.cxx b/src/ParaMEDMEM/IntersectionDEC.cxx new file mode 100644 index 000000000..346aede57 --- /dev/null +++ b/src/ParaMEDMEM/IntersectionDEC.cxx @@ -0,0 +1,115 @@ +#include +#include "CommInterface.hxx" +#include "Topology.hxx" +#include "BlockTopology.hxx" +#include "ComponentTopology.hxx" +#include "ParaFIELD.hxx" +#include "MPIProcessorGroup.hxx" +#include "ParaMESH.hxx" +#include "ParaSUPPORT.hxx" +#include "DEC.hxx" +#include "InterpolationMatrix.hxx" +#include "IntersectionDEC.hxx" +#include "ElementLocator.hxx" + + +namespace ParaMEDMEM +{ + +IntersectionDEC::IntersectionDEC() +{ +} + +IntersectionDEC::IntersectionDEC(ProcessorGroup& local_group, ProcessorGroup& distant_group): + DEC(local_group, distant_group),_interpolation_matrix(0) +{ +} + +IntersectionDEC::~IntersectionDEC() +{ + if (_interpolation_matrix !=0) + delete _interpolation_matrix; +} + +/*! Synchronization process for exchanging topologies + */ +void IntersectionDEC::synchronize() +{ + + //setting up the communication DEC on both sides + if (_source_field!=0) + { + 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(); + // 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) + { + 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(); + + //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(); +} + + +void IntersectionDEC::recvData() +{ + _interpolation_matrix->multiply(*_target_field->getField()); +} + +void IntersectionDEC::sendData() +{ + _interpolation_matrix->multiply(*_source_field->getField()); +} + +} + + + diff --git a/src/ParaMEDMEM/IntersectionDEC.hxx b/src/ParaMEDMEM/IntersectionDEC.hxx new file mode 100644 index 000000000..ae9f6bc37 --- /dev/null +++ b/src/ParaMEDMEM/IntersectionDEC.hxx @@ -0,0 +1,46 @@ +#ifndef INTERSECTIONDEC_HXX_ +#define INTERSECTIONDEC_HXX_ + + +namespace ParaMEDMEM +{ + class DEC; + class InterpolationMatrix; + + class IntersectionDEC:public DEC + { + public: + IntersectionDEC(); + IntersectionDEC(ProcessorGroup& local_group, ProcessorGroup& distant_group); + + virtual ~IntersectionDEC(); + + void synchronize(); + + void recvData(); + + void sendData(); + + void prepareSourceDE(){}; + void prepareTargetDE(){}; + + private : + + //Number of distant points to be located locally + int _nb_distant_points; + + //coordinates of distant points + const double* _distant_coords; + + //local element number containing the distant points + const int* _distant_locations; + + //inerpolation method + string _method; + + InterpolationMatrix* _interpolation_matrix; + }; +} + + +#endif diff --git a/src/ParaMEDMEM/MPIProcessorGroup.cxx b/src/ParaMEDMEM/MPIProcessorGroup.cxx index 3f6330753..f433bcc81 100644 --- a/src/ParaMEDMEM/MPIProcessorGroup.cxx +++ b/src/ParaMEDMEM/MPIProcessorGroup.cxx @@ -38,8 +38,6 @@ ProcessorGroup(interface, proc_ids) // copying proc_ids in ranks copy::const_iterator,int*> (proc_ids.begin(), proc_ids.end(), ranks); - for (int i=0; i procs = _proc_ids; + const set& distant_proc_ids = group.getProcIDs(); + for (set::const_iterator iter=distant_proc_ids.begin(); iter!=distant_proc_ids.end(); iter++) + { + procs.insert(*iter); + } + return new MPIProcessorGroup(_comm_interface,procs); +} + } diff --git a/src/ParaMEDMEM/MPIProcessorGroup.hxx b/src/ParaMEDMEM/MPIProcessorGroup.hxx index aa0d49ea4..8ebb98908 100644 --- a/src/ParaMEDMEM/MPIProcessorGroup.hxx +++ b/src/ParaMEDMEM/MPIProcessorGroup.hxx @@ -17,7 +17,7 @@ public: MPIProcessorGroup(const CommInterface& interface, set proc_ids); MPIProcessorGroup (const ProcessorGroup& proc_group, set proc_ids); virtual ~MPIProcessorGroup(); - void fuse (ProcessorGroup&){}; + virtual ProcessorGroup* fuse (const ProcessorGroup&) const; void intersect (ProcessorGroup&){}; int myRank() const {int rank; MPI_Comm_rank(_comm,&rank); return rank;} bool containsMyRank() const { int rank; MPI_Group_rank(_group, &rank); return (rank!=MPI_UNDEFINED);} diff --git a/src/ParaMEDMEM/Makefile.in b/src/ParaMEDMEM/Makefile.in index 3cd41cd22..2778fa0e8 100644 --- a/src/ParaMEDMEM/Makefile.in +++ b/src/ParaMEDMEM/Makefile.in @@ -52,9 +52,14 @@ ComponentTopology.hxx\ ExplicitTopology.hxx\ ParaFIELD.hxx\ DEC.hxx\ +MxN_Mapping.hxx\ +INTERPOLATION_2D.hxx\ StructuredCoincidentDEC.hxx\ +InterpolationMatrix.hxx\ +IntersectionDEC.hxx\ UnstructuredParaSUPPORT.hxx\ ExplicitCoincidentDEC.hxx\ +ElementLocator.hxx\ ExplicitMapping.hxx # Libraries targets @@ -71,9 +76,14 @@ StructuredParaSUPPORT.cxx\ ComponentTopology.cxx\ ParaFIELD.cxx\ DEC.cxx\ +INTERPOLATION_2D.cxx\ +MxN_Mapping.cxx\ +InterpolationMatrix.cxx\ StructuredCoincidentDEC.cxx\ ExplicitCoincidentDEC.cxx\ +IntersectionDEC.cxx\ UnstructuredParaSUPPORT.cxx\ +ElementLocator.cxx\ ExplicitTopology.cxx @@ -85,12 +95,12 @@ BIN_SERVER_IDL = BIN_CLIENT_IDL = TEST_PROGS = test_ProcessorGroup test_BlockTopology test_ParaStructuredSupport \ -test_ParaField test_DEC test_UnstructuredDEC test_ExplicitDEC +test_ParaField test_DEC test_UnstructuredDEC test_ExplicitDEC test_IntersectionDEC LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome -CPPFLAGS+=$(MED2_INCLUDES) $(MPI_INCLUDES) $(LAM_INCLUDES) -I/data/tmpawa/vb144235/lam_install/include +CPPFLAGS+=$(MED2_INCLUDES) $(MPI_INCLUDES) $(LAM_INCLUDES) -I/data/tmpawa/vb144235/lam_install/include -I/data/tmpawa/vb144235/fvm/src -DFVM_HAVE_MPI CXXFLAGS+=@CXXTMPDPTHFLAGS@ CPPFLAGS+=$(BOOST_CPPFLAGS) diff --git a/src/ParaMEDMEM/MxN_Mapping.cxx b/src/ParaMEDMEM/MxN_Mapping.cxx new file mode 100644 index 000000000..ea93fa0c4 --- /dev/null +++ b/src/ParaMEDMEM/MxN_Mapping.cxx @@ -0,0 +1,243 @@ +#include "CommInterface.hxx" +#include "ProcessorGroup.hxx" +#include "MPIProcessorGroup.hxx" +#include "MxN_Mapping.hxx" + +namespace ParaMEDMEM +{ + +MxN_Mapping::MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group) +: _union_group(local_group.fuse(distant_group)) +{ + _send_proc_offsets.resize(_union_group->size()+1,0); + _recv_proc_offsets.resize(_union_group->size()+1,0); + +} + +MxN_Mapping::~MxN_Mapping() +{ + delete _union_group; +} + +void MxN_Mapping::addElementFromSource(int local_element, int distant_proc, int distant_element) +{ + _sending_ids.push_back(make_pair(distant_proc,distant_element)); + // _local_ids.push_back(local_element); + for (int i=distant_proc; i<_union_group->size(); i++) + _send_proc_offsets[i+1]++; +} + +void MxN_Mapping::prepareSendRecv() +{ + CommInterface comm_interface=_union_group->getCommInterface(); + // sending count pattern + int* nbsend=new int[_union_group->size()]; + int* nbrecv=new int[_union_group->size()]; + for (int i=0; i<_union_group->size(); i++) + { + nbsend[i]=_send_proc_offsets[i+1]-_send_proc_offsets[i]; + } + + MPIProcessorGroup* group = static_cast(_union_group); + const MPI_Comm* comm=group->getComm(); + comm_interface.allToAll(nbsend, 1, MPI_INT, + nbrecv, 1, MPI_INT, + *comm); + + for (int i=0; i<_union_group->size(); i++) + { + for (int j=i+1;j<_union_group->size()+1; j++) + _recv_proc_offsets[j]+=nbrecv[i]; + + } + + delete[] nbsend; + delete[] nbrecv; + + _recv_ids.resize(_recv_proc_offsets[_union_group->size()]); + int* isendbuf=0; + int* irecvbuf=0; + if (_sending_ids.size()>0) + isendbuf = new int[_sending_ids.size()]; + if (_recv_ids.size()>0) + irecvbuf = new int[_recv_ids.size()]; + 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]=_send_proc_offsets[i+1]-_send_proc_offsets[i]; + senddispls[i]=_send_proc_offsets[i]; + recvcounts[i]=_recv_proc_offsets[i+1]-_recv_proc_offsets[i]; + recvdispls[i]=_recv_proc_offsets[i]; + } + vector offsets = _send_proc_offsets; + for (int i=0; i<_sending_ids.size();i++) + { + int iproc = _sending_ids[i].first; + isendbuf[offsets[iproc]]=_sending_ids[i].second; + offsets[iproc]++; + } + comm_interface.allToAllV(isendbuf, sendcounts, senddispls, MPI_INT, + irecvbuf, recvcounts, recvdispls, MPI_INT, + *comm); + + for (int i=0; i< _recv_proc_offsets[_union_group->size()]; i++) + _recv_ids[i]=irecvbuf[i]; + + if (_sending_ids.size()>0) + delete[] isendbuf; + if (_recv_ids.size()>0) + delete[] irecvbuf; + delete[] sendcounts; + delete[]recvcounts; + delete[]senddispls; + delete[] recvdispls; +} + + +/*! 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::sendRecv(MEDMEM::FIELD& field) +{ + CommInterface comm_interface=_union_group->getCommInterface(); + const MPIProcessorGroup* group = static_cast(_union_group); + + int nbcomp=field.getNumberOfComponents(); + double* sendbuf=0; + double* recvbuf=0; + if (_sending_ids.size() >0) + sendbuf = new double[_sending_ids.size()*nbcomp]; + if (_recv_ids.size()>0) + recvbuf = new double[_recv_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*(_send_proc_offsets[i+1]-_send_proc_offsets[i]); + senddispls[i]=nbcomp*(_send_proc_offsets[i]); + recvcounts[i]=nbcomp*(_recv_proc_offsets[i+1]-_recv_proc_offsets[i]); + recvdispls[i]=nbcomp*(_recv_proc_offsets[i]); + } + //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; + for (int icomp=0; icompgetComm(); + comm_interface.allToAllV(sendbuf, sendcounts, senddispls, MPI_INT, + recvbuf, recvcounts, recvdispls, MPI_INT, + *comm); + + + //setting the received values in the field + double* recvptr; + for (int i=0; i< _recv_proc_offsets[_union_group->size()]; i++) + { + for (int icomp=0; icomp& 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 (_sending_ids.size() >0) + sendbuf = new double[_sending_ids.size()*nbcomp]; + if (_recv_ids.size()>0) + recvbuf = new double[_recv_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*(_send_proc_offsets[i+1]-_send_proc_offsets[i]); + senddispls[i]=nbcomp*(_send_proc_offsets[i]); + recvcounts[i]=nbcomp*(_recv_proc_offsets[i+1]-_recv_proc_offsets[i]); + recvdispls[i]=nbcomp*(_recv_proc_offsets[i]); + } + //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; + 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 + +#include "MEDMEM_Field.hxx" + + +namespace ParaMEDMEM +{ +class ProcessorGroup; + +class MxN_Mapping +{ +public: + MxN_Mapping(); + MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group); + virtual ~MxN_Mapping(); + void addElementFromSource(int local_elem, int distant_proc, int distant_elem); + void prepareSendRecv(); + void sendRecv(MEDMEM::FIELD& field); + void sendRecv(double* field, MEDMEM::FIELD& field) const ; + +private : +// ProcessorGroup& _local_group; +// ProcessorGroup& _distant_group; + ProcessorGroup* _union_group; + int _nb_comps; + std::vector > _sending_ids; + std::vector _recv_ids; + std::vector _send_proc_offsets; + std::vector _recv_proc_offsets; + +}; + +} + +#endif /*MXN_MAPPING_HXX_*/ diff --git a/src/ParaMEDMEM/NonCoincidentDEC.cxx b/src/ParaMEDMEM/NonCoincidentDEC.cxx new file mode 100644 index 000000000..c4617d1ec --- /dev/null +++ b/src/ParaMEDMEM/NonCoincidentDEC.cxx @@ -0,0 +1,238 @@ +#include +#include "CommInterface.hxx" +#include "Topology.hxx" +#include "BlockTopology.hxx" +#include "ComponentTopology.hxx" +#include "ParaFIELD.hxx" +#include "MPIProcessorGroup.hxx" +#include "DEC.hxx" +#include "NonCoincidentDEC.hxx" + +extern "C" { +#include +#include +#include +#include +} + +namespace ParaMEDMEM +{ + fvm_nodal_t* medmemMeshToFVMMesh(const MEDMEM::MESH* mesh) + { + // create an FVM structure from the paramesh structure + fvm_nodal_t * fvm_nodal = fvm_nodal_create(mesh->getName().c_str()); + + //loop on cell types + int nbtypes = mesh->getNumberOfTypes(MED_EN::MED_CELL); + const MED_EN::medGeometryElement* types = mesh->getTypes(MED_EN::MED_CELL); + for (int itype=0; itypegetNumberOfElements(MED_EN::MED_CELL, types[itype]); + fvm_lnum_t* conn = const_cast (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, types[itype])); + + fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,0); + + } + return fvm_nodal; + } + + fvm_nodal_t* medmemSupportToFVMMesh(const MEDMEM::SUPPORT* support) + { + + // create an FVM structure from the paramesh structure + fvm_nodal_t * fvm_nodal = fvm_nodal_create(support->getName().c_str()); + + const MEDMEM::MESH* mesh= support->getMesh(); + + //loop on cell types + MED_EN::medEntityMesh entity = support->getEntity(); + + int nbtypes = support->getNumberOfTypes(); + const MED_EN::medGeometryElement* types = support->getTypes(); + int ioffset=0; + const int* type_offset = support->getNumberIndex(); + + //browsing through all types + for (int itype=0; itypegetNumberOfElements(types[itype]); + + //for a partial support, defining the element numbers that are taken into + //account in the support + fvm_lnum_t* elem_numbers=0; + if (!support->isOnAllElements()) + { + elem_numbers = const_cast (support->getNumber(types[itype])); + + //creating work arrays to store list of elems for partial suports + if (itype>0) + { + fvm_lnum_t* temp = new int[nbelems]; + for (int i=0; i< nbelems; i++) + temp[i] = elem_numbers [i]-ioffset; + ioffset+=type_offset[itype]; + elem_numbers = temp; + } + } + //retrieving original mesh connectivity + fvm_lnum_t* conn = const_cast (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,entity, types[itype])); + + // adding the elements to the FVM structure + fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,elem_numbers); + + //cleaning work arrays (for partial supports) + if (!support->isOnAllElements() && itype>0) + delete[] elem_numbers; + + } + return fvm_nodal; + } + +NonCoincidentDEC::NonCoincidentDEC() +{ +} + +NonCoincidentDEC::~NonCoincidentDEC() +{ +} + +/*! Synchronization process for exchanging topologies + */ +void NonCoincidentDEC::synchronize() +{ + + //initializing FVM parallel environment + fvm_parall_set_mpi_comm(MPI_COMM_WORLD); + + + //setting up the communication DEC on both sides + + if (_source_field!=0) + { + MEDMEM::MESH* mesh = _source_field->getField()->getSupport()->getMesh(); + fvm_nodal_t* source_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh); + + const ProcessorGroup* proc_group = _source_field->getTopology()->getProcGroup(); + int target_size = proc_group->getCommInterface().worldSize() - + proc_group->size() ; + int start_rank= proc_group->size(); + const MPI_Comm* comm = (dynamic_cast (proc_group))->getComm(); + fvm_locator_t* locator = + fvm_locator_create(1e-6, + *comm, + target_size, + start_rank); + + fvm_locator_set_nodal(locator, + source_nodal, + 2, + 0, + 0, + 0); + + + _nb_distant_points = fvm_locator_get_n_dist_points(locator); + _distant_coords = fvm_locator_get_dist_coords(locator); + _distant_locations = fvm_locator_get_dist_locations(locator); + + } + if (_target_field!=0) + { + MEDMEM::MESH* mesh = _target_field->getField()->getSupport()->getMesh(); + + fvm_nodal_t* target_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh); + const ProcessorGroup* proc_group = _target_field->getTopology()->getProcGroup(); + int source_size = proc_group->getCommInterface().worldSize() -proc_group->size() ; + int start_rank= 0 ; + const MPI_Comm* comm = (dynamic_cast (proc_group))->getComm(); + fvm_locator_t* locator = + fvm_locator_create(1e-6, + *comm, + source_size, + start_rank); + int nbcells = mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS); + SUPPORT support(mesh, MED_EN::MED_CELL); + const double* coords = (mesh->getBarycenter(&support))->getValue(); + fvm_locator_set_nodal(locator, + target_nodal, + 2, + nbcells, + 0, + coords); + } +} + + +void NonCoincidentDEC::recvData() +{ + int nbelems = _target_field->getField()->getSupport()->getMesh()->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS); + int nbcomp = _target_field->getField()->getNumberOfComponents(); + double* values = new double [nbelems*nbcomp]; + fvm_locator_exchange_point_var(_locator, + 0, + values, + sizeof(double), + nbcomp); + _target_field->getField()->setValue(values); +} + +void NonCoincidentDEC::sendData() +{ + const double* values=_source_field->getField()->getValue(); + int nbcomp = _source_field->getField()->getNumberOfComponents(); + double* distant_values = new double [_nb_distant_points*nbcomp]; + for (int i=0; i<_nb_distant_points; i++) + for (int j=0; j + +typedef enum {NN} InterpolationMethod; + +namespace ParaMEDMEM +{ + class DEC; + + class NonCoincidentDEC:public DEC + { + public: + NonCoincidentDEC(); + + virtual ~NonCoincidentDEC(); + + void synchronize(); + + void recvData(); + + void sendData(); + + void prepareSourceDE(){}; + void prepareTargetDE(){}; + + void setInterpolationMethod(InterpolationMethod method) + {_method=method;} + + private : + // Structure for computing the localization + // of remote nodes on local mesh + fvm_locator_t* _locator; + + //Number of distant points to be located locally + int _nb_distant_points; + + //coordinates of distant points + const double* _distant_coords; + + //local element number containing the distant points + const int* _distant_locations; + + //inerpolation method + InterpolationMethod _method; + }; +} + + +#endif diff --git a/src/ParaMEDMEM/ParaFIELD.cxx b/src/ParaMEDMEM/ParaFIELD.cxx index 3251a61db..07aada33d 100644 --- a/src/ParaMEDMEM/ParaFIELD.cxx +++ b/src/ParaMEDMEM/ParaFIELD.cxx @@ -4,8 +4,10 @@ #include "ComponentTopology.hxx" #include "ParaSUPPORT.hxx" #include "StructuredParaSUPPORT.hxx" +#include "UnstructuredParaSUPPORT.hxx" #include "ExplicitCoincidentDEC.hxx" #include "StructuredCoincidentDEC.hxx" + #include "ParaFIELD.hxx" #include "ParaMESH.hxx" @@ -16,7 +18,8 @@ namespace ParaMEDMEM ParaFIELD::ParaFIELD(const ParaSUPPORT* para_support, const ComponentTopology& component_topology) :_support(para_support), -_component_topology(component_topology) + _component_topology(component_topology), + _field(0) { if (dynamic_cast(para_support)!=0) {const BlockTopology* source_topo = dynamic_cast(para_support->getTopology()); @@ -61,12 +64,34 @@ _component_topology(component_topology) _field->setIterationNumber(0); _field->setOrderNumber(0); _field->setTime(0.0); + delete[] compnames; + delete[] compdesc; } +/*! Constructor creating the ParaFIELD + * from a given FIELD + */ +ParaFIELD::ParaFIELD(MEDMEM::FIELD* subdomain_field, const ProcessorGroup& proc_group): +_field(subdomain_field), +_support(new UnstructuredParaSUPPORT(subdomain_field->getSupport(), proc_group)), +_component_topology(ComponentTopology(_field->getNumberOfComponents())) +{ + const ExplicitTopology* source_topo= + dynamic_cast (_support->getTopology()); + _topology=new ExplicitTopology(*source_topo,_component_topology.nbLocalComponents()); +} + ParaFIELD::ParaFIELD(MEDMEM::driverTypes driver_type, const string& file_name, const string& driver_name, const ComponentTopology& component_topology) - throw (MEDEXCEPTION):_component_topology(component_topology){} -ParaFIELD::~ParaFIELD(){} + throw (MEDEXCEPTION):_component_topology(component_topology),_field(0){} + +ParaFIELD::~ParaFIELD() +{ + if (_field!=0) + delete _field; + if (_topology!=0) + delete _topology; +} void ParaFIELD::write(MEDMEM::driverTypes driverType, const string& fileName, const string& meshName){ // Topology* topo = dynamic_cast (_topology); diff --git a/src/ParaMEDMEM/ParaFIELD.hxx b/src/ParaMEDMEM/ParaFIELD.hxx index 7e5fac29d..f560066fe 100644 --- a/src/ParaMEDMEM/ParaFIELD.hxx +++ b/src/ParaMEDMEM/ParaFIELD.hxx @@ -7,6 +7,7 @@ namespace MEDMEM{ class MEDEXCEPTION; + template class FIELD; } @@ -14,6 +15,7 @@ namespace ParaMEDMEM { class ComponentTopology; class ParaSUPPORT; +class ProcessorGroup; class ParaFIELD { @@ -24,11 +26,14 @@ public: ParaFIELD(MEDMEM::driverTypes driver_type, const string& file_name, const string& driver_name, const ComponentTopology& component_topology) throw (MEDMEM::MEDEXCEPTION); + ParaFIELD(MEDMEM::FIELD* field, const ProcessorGroup& group); + virtual ~ParaFIELD(); void write(MEDMEM::driverTypes driverType, const string& fileName="", const string& meshName=""); void synchronizeTarget(ParaFIELD* source_field); void synchronizeSource(ParaFIELD* target_field); MEDMEM::FIELD* getField() const {return _field;} + const ParaSUPPORT* getSupport() const {return _support;} Topology* getTopology() const {return _topology;} int nbComponents() const {return _component_topology.nbComponents();} private: diff --git a/src/ParaMEDMEM/ParaMESH.cxx b/src/ParaMEDMEM/ParaMESH.cxx index 45d1579a4..6b9a5f064 100644 --- a/src/ParaMEDMEM/ParaMESH.cxx +++ b/src/ParaMEDMEM/ParaMESH.cxx @@ -48,13 +48,14 @@ throw (MEDMEM::MEDEXCEPTION){ string mesh; int idomain; string host; - + string medfilename; + for (int i=0; i<=domain_id;i++) { //reading information about the domain - asciiinput >> mesh >> idomain >> meshstring >> host >> _medfilename; + asciiinput >> mesh >> idomain >> meshstring >> host >> medfilename; if (idomain!=i+1) { @@ -62,7 +63,7 @@ throw (MEDMEM::MEDEXCEPTION){ throw (MEDEXCEPTION("Error : domain must be written from 1 to N in asciifile descriptor")); } strcpy(meshname,meshstring.c_str()); - strcpy(file,_medfilename.c_str()); + strcpy(file,medfilename.c_str()); } _name=meshstring; /////////////////////////////////////////// @@ -80,7 +81,7 @@ throw (MEDMEM::MEDEXCEPTION){ char joint_description[MED_TAILLE_DESC]; char name[MED_TAILLE_NOM]; char name_distant[MED_TAILLE_NOM]; - cout << "arguments"<< fid<<" "<localToGlobal(make_pair(_my_domain_id,0)); + for (int i=0; i proc_ids): _comm_interface(proc_group.getCommInterface()){} virtual ~ProcessorGroup(){} - virtual void fuse (ProcessorGroup&)=0; + virtual ProcessorGroup* fuse (const ProcessorGroup&) const=0; virtual void intersect (ProcessorGroup&)=0; bool contains(int rank) const {return _proc_ids.find(rank)!=_proc_ids.end();}; virtual bool containsMyRank() const=0; @@ -28,6 +28,7 @@ public: virtual int translateRank(const ProcessorGroup*, int) const =0; virtual ProcessorGroup* createComplementProcGroup() const =0; virtual ProcessorGroup* createProcGroup() const=0; + virtual const std::set& getProcIDs()const {return _proc_ids;} protected: const CommInterface& _comm_interface; std::set _proc_ids; diff --git a/src/ParaMEDMEM/UnstructuredParaSUPPORT.cxx b/src/ParaMEDMEM/UnstructuredParaSUPPORT.cxx index f793bc72f..300fb9283 100644 --- a/src/ParaMEDMEM/UnstructuredParaSUPPORT.cxx +++ b/src/ParaMEDMEM/UnstructuredParaSUPPORT.cxx @@ -1,5 +1,6 @@ #include "Topology.hxx" #include "ExplicitTopology.hxx" +#include "BlockTopology.hxx" #include "ParaMESH.hxx" #include "UnstructuredParaSUPPORT.hxx" #include "MEDMEM_Support.hxx" @@ -7,13 +8,25 @@ namespace ParaMEDMEM { -/*! Constructor on all elements from a MESH */ +/*! Constructor on from a ParaMESH and a local support*/ UnstructuredParaSUPPORT::UnstructuredParaSUPPORT(const ParaMESH* const mesh, const MEDMEM::SUPPORT* support): -_entity(support->getEntity()), -_explicit_topology(new ExplicitTopology(*support)) +ParaSUPPORT(mesh,support), +_entity(support->getEntity()) { _mesh=mesh; _support=support; + _explicit_topology=new ExplicitTopology(*(ParaSUPPORT*)this); +} + +/*! Constructor on from a ProcessorGroup and a local support*/ +UnstructuredParaSUPPORT::UnstructuredParaSUPPORT(const MEDMEM::SUPPORT* support, const ProcessorGroup& proc_group): +_entity(support->getEntity()) +{ + ostringstream name ("mesh associated with support "); + name << support->getName(); + _mesh=new ParaMESH(*support->getMesh(), proc_group, name.str() ); + _support=support; + _explicit_topology=new ExplicitTopology(*(ParaSUPPORT*)this); } UnstructuredParaSUPPORT::UnstructuredParaSUPPORT(const ParaMESH* const mesh, const MED_EN::medEntityMesh entity): @@ -21,14 +34,14 @@ UnstructuredParaSUPPORT::UnstructuredParaSUPPORT(const ParaMESH* const mesh, con _entity(entity), _explicit_topology(new ExplicitTopology(*this)) { - //_mesh=mesh; - // _support=new SUPPORT(_mesh->getMesh(), "support on all entities", entity); - //_explicit_topology(new ExplicitTopology(*support)); } UnstructuredParaSUPPORT::~UnstructuredParaSUPPORT() { delete _support; + delete _explicit_topology; } + + }//end of namespace ParaMEDMEM diff --git a/src/ParaMEDMEM/UnstructuredParaSUPPORT.hxx b/src/ParaMEDMEM/UnstructuredParaSUPPORT.hxx index cc2ad70b8..5a82fb7c1 100644 --- a/src/ParaMEDMEM/UnstructuredParaSUPPORT.hxx +++ b/src/ParaMEDMEM/UnstructuredParaSUPPORT.hxx @@ -23,12 +23,14 @@ namespace ParaMEDMEM public: UnstructuredParaSUPPORT(const ParaMESH* const mesh, const MEDMEM::SUPPORT* support ); - UnstructuredParaSUPPORT(const ParaMESH* const mesh, const MED_EN::medEntityMesh entity); - virtual ~UnstructuredParaSUPPORT(); - const Topology* getTopology() const {return _explicit_topology;} + UnstructuredParaSUPPORT(const MEDMEM::SUPPORT* support , const ProcessorGroup& group); + UnstructuredParaSUPPORT(const ParaMESH* const mesh, const MED_EN::medEntityMesh entity); + virtual ~UnstructuredParaSUPPORT(); + const Topology* UnstructuredParaSUPPORT::getTopology() const + {return _explicit_topology;} private: - const ExplicitTopology* const _explicit_topology; + const ExplicitTopology* _explicit_topology; const MED_EN::medEntityMesh _entity; }; diff --git a/src/ParaMEDMEM/test_IntersectionDEC.cxx b/src/ParaMEDMEM/test_IntersectionDEC.cxx new file mode 100644 index 000000000..bf849b42b --- /dev/null +++ b/src/ParaMEDMEM/test_IntersectionDEC.cxx @@ -0,0 +1,176 @@ +#include +#include +#include +#include +#include + +#include "ProcessorGroup.hxx" +#include "Topology.hxx" +#include "BlockTopology.hxx" +#include "CommInterface.hxx" + + +#include "MPIProcessorGroup.hxx" +#include "MEDMEM_Mesh.hxx" +#include "MEDMEM_Group.hxx" +#include "MEDMEM_Support.hxx" +#include "ParaMESH.hxx" +#include "UnstructuredParaSUPPORT.hxx" +#include "ComponentTopology.hxx" +#include "ParaFIELD.hxx" +#include "IntersectionDEC.cxx" + +using namespace std; +using namespace ParaMEDMEM; +using namespace MEDMEM; +int main(int argc, char** argv) +{ + string testname="ParaMEDMEM - test #1 -"; + MPI_Init(&argc, &argv); + int size; + int rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + if (argc !=2 || size < 2) + { + cout << "usage :"< test_NonCoincidentDEC "<=2)"< self_procs; + set procs_source; + set procs_target; + + for (int i=0; icontainsMyRank()) + { + string master = "/home/vb144235/resources/square128000_split"; + cout <<"loading source"<getBlockTopology(); + ostringstream strstream; + strstream < mesh_groups; + // + // for (int i=0; igetName() == "boundary") + // boundary_group = mesh_groups[i]; + // } + // + boundary_group = new SUPPORT(mesh_bis, "support on all cells", MED_EN::MED_CELL); + } + + //loading the geometry for the target group + if (target_group->containsMyRank()) + { + string master = "/home/vb144235/resources/square30000_split"; + target_mesh=new ParaMESH(MED_DRIVER,master,*target_group); + topo_target=target_mesh->getBlockTopology(); + ostringstream strstream; + strstream << master<<(rank-nproc_source+1)<<".med"; + //strstream<<"/home/vb144235/resources/square30000.med"; + ostringstream meshname ; + meshname<< "Mesh_1_"< mesh_groups; + // + // for (int i=0; igetName() == "boundary") + // boundary_group = mesh_groups[i]; + // } + boundary_group = new SUPPORT(mesh_bis, "support on all cells", MED_EN::MED_CELL); + + } + + + //attaching a DEC to the source group + if (source_group->containsMyRank()) + { + + //UnstructuredParaSUPPORT source_support(boundary_group, source_group); + UnstructuredParaSUPPORT source_support(source_mesh_bis, boundary_group); + ComponentTopology source_comp(1); + cout << "setting up field"<getNbLocalElements(); + cout << "Source field nb elems on rank : "<setValue(value); + cout <<"creating intersectionDEC"<write(MED_DRIVER,"/home/vb144235/tmp/sourcesquare"); + source_field.write(MED_DRIVER,"/home/vb144235/tmp/sourcesquare","boundary"); + delete[] value; + } + + //attaching a DEC to the target group + if (target_group->containsMyRank()) + { + + UnstructuredParaSUPPORT target_support(target_mesh_bis, boundary_group); + //UnstructuredParaSUPPORT target_support(boundary_group, target_group); + ComponentTopology target_comp(1); + ParaFIELD target_field(&target_support,target_comp); + int nb_local = target_field.getTopology()->getNbLocalElements(); + double * value= new double[nb_local]; + for(int ielem=0; ielemsetValue(value); + IntersectionDEC dec(*source_group,*target_group); + dec.attachTargetField(&target_field); + dec.synchronize(); + dec.recvData(); + target_mesh->write(MED_DRIVER, "/home/vb144235/tmp/targetsquare"); + target_field.write(MED_DRIVER, "/home/vb144235/tmp/targetsquare", "boundary"); + delete[] value; + } + + MPI_Barrier(MPI_COMM_WORLD); + MPI_Finalize(); + return 0; +} + + diff --git a/src/ParaMEDMEM/test_NonCoincidentDEC.cxx b/src/ParaMEDMEM/test_NonCoincidentDEC.cxx new file mode 100644 index 000000000..6246c3417 --- /dev/null +++ b/src/ParaMEDMEM/test_NonCoincidentDEC.cxx @@ -0,0 +1,155 @@ +#include +#include +#include +#include +#include + +#include "ProcessorGroup.hxx" +#include "Topology.hxx" +#include "BlockTopology.hxx" +#include "CommInterface.hxx" + + +#include "MPIProcessorGroup.hxx" +#include "MEDMEM_Mesh.hxx" +#include "ParaMESH.hxx" +#include "UnstructuredParaSUPPORT.hxx" +#include "ComponentTopology.hxx" +#include "ParaFIELD.hxx" + + +using namespace std; +using namespace ParaMEDMEM; +using namespace MEDMEM; +int main(int argc, char** argv) +{ + string testname="ParaMEDMEM - test #1 -"; + MPI_Init(&argc, &argv); + int size; + int rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + if (argc !=2 || size < 2) + { + cout << "usage :"< test_NonCoincidentDEC "<=2)"< self_procs; + set procs_source; + set procs_target; + + for (int i=0; icontainsMyRank()) + { + string master = "/home/vb144235/resources/fluid_split"; + source_mesh=new ParaMESH(MED_DRIVER,master,*source_group); + topo_source=source_mesh->getBlockTopology(); + ostringstream strstream (master); + strstream << "_"< mesh_groups; + const MEDMEM::GROUP* boundary_group; + for (int i=0; icontainsMyRank()) + { + string master = "/home/vb144235/resources/blade_split"; + target_mesh=new ParaMESH(MED_DRIVER,master,*target_group); + topo_target=target_mesh->getBlockTopology(); + ostringstream strstream (master); + strstream << "_"<<(rank-nproc_source)< mesh_groups; + const MEDMEM::GROUP* boundary_group; + for (int i=0; icontainsMyRank()) + { + + //UnstructuredParaSUPPORT source_support(boundary_group, source_group); + UnstructuredParaSUPPORT source_support(source_mesh_bis, boundary_group); + ComponentTopology source_comp(1); + ParaFIELD source_field (source_support,*source_comp); + + int nb_local = source_field.getTopology()->getNbLocalElements(); + cout << "Source field nb elems on rank : "<setValue(value); + NonCoincidentDEC dec; + dec.attachSourceField(&source_field); + dec.synchronize(); + dec.sendData(); + + source_mesh.write(MED_DRIVER,"/home/vb144235/tmp/sourceexp"); + source_field.write(MED_DRIVER,"/home/vb144235/tmp/sourceexp","boundary"); + } + + //attaching a DEC to the target group + if (target_group->containsMyRank()) + { + UnstructuredParaSUPPORT target_support(target_mesh_bis, boundary_group); + //UnstructuredParaSUPPORT target_support(boundary_group, target_group); + ComponentTopology target_comp(1); + ParaFIELD target_field(target_support,*target_comp); + NonCoincidentDEC dec; + dec.attachTargetField(&target_field); + dec.synchronize(); + dec.recvData(); + target_mesh.write(MED_DRIVER, "/home/vb144235/tmp/targetexp"); + target_field.write(MED_DRIVER, "/home/vb144235/tmp/targetexp", "boundary"); + } + + MPI_Barrier(MPI_COMM_WORLD); + MPI_Finalize(); + return 0; +} + + + +