}
_cycle_type.resize(1);
_cycle_type[0]=ParaMEDMEM::Block;
-
+ delete[] nbelems_per_proc;
}
BlockTopology::~BlockTopology()
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;}
};
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<BlockTopology*>(field->getTopology());
- _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
- //}
+ _target_field=field;
+ //if (field!=0)
+ //{
+ //BlockTopology* topo=dynamic_cast<BlockTopology*>(field->getTopology());
+ _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
+ //}
}
void DEC::attachSourceField(const ParaFIELD* field)
-{_source_field=field;
-//if (field!=0)
-//{
-// BlockTopology* topo=dynamic_cast<BlockTopology*>(field->getTopology());
- _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
- //}
+{
+ _source_field=field;
+ //if (field!=0)
+ //{
+ // BlockTopology* topo=dynamic_cast<BlockTopology*>(field->getTopology());
+ _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
+ //}
}
}
{
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;
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;
};
--- /dev/null
+#include <mpi.h>
+#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 <set>
+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 <int> 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; i<dim; i++)
+ {
+ elem_bb[i*2]=HUGE;
+ elem_bb[i*2+1]=-HUGE;
+ }
+ for (int inode=conn_index[ielem]; inode<conn_index[ielem+1]; inode++)
+ {
+ int node= conn [inode-1]-1; // - 1 because of MED numbering starting at 1
+
+ for (int idim=0; idim<dim; idim++)
+ {
+ elem_bb[idim*2]=coords[node*dim+idim]<elem_bb[idim*2]? coords[node*dim+idim]: elem_bb[idim*2];
+ elem_bb[idim*2+1]=coords[node*dim+idim]>elem_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<int>::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; idim<dim; idim++)
+ {
+ minmax[idim*2]=HUGE;
+ minmax[idim*2+1]=-HUGE;
+ }
+ for (int i=0; i<nbnodes; i++)
+ for (int idim=0; idim<dim;idim++)
+ {
+ minmax[idim*2]=(minmax[idim*2]<coords[i*dim+idim]?minmax[idim*2]:coords[i*dim+idim]);
+ minmax[idim*2+1]=(minmax[idim*2+1]>coords[i*dim+idim]?minmax[idim*2+1]:coords[i*dim+idim]);
+ }
+ MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_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]<local_bb[idim*2+1] && local_bb[idim*2]<distant_bb[idim*2+1];
+ if (!intersects) return false;
+ }
+ return true;
+}
+
+bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim)
+{
+ for (int idim=0; idim < dim; idim++)
+ {
+ bool intersects = bb1[idim*2]<bb2[idim*2+1] && bb2[idim*2]<bb1[idim*2+1];
+ if (!intersects) return false;
+ }
+ return true;
+}
+
+void ElementLocator::_exchangeMesh(MEDMEM::MESH* local_mesh, MEDMEM::MESH*& distant_mesh, int iproc_distant, const int* distant_ids_send, int*& distant_ids_recv)
+{
+
+ CommInterface comm_interface=_union_group->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<MPIProcessorGroup*> (_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<double*> (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; i<distant_nb_types; i++)
+ types_vector[i]=(MED_EN::medGeometryElement)recv_buffer_ptr[i];
+ meshing->setTypes(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; i<distant_nb_types; i++)
+ nbtypes[i]=recv_buffer_ptr[i+1]-recv_buffer_ptr[i];
+ recv_buffer_ptr+=distant_nb_types+1;
+ meshing->setNumberOfElements( nbtypes, MED_EN::MED_CELL);
+
+ for (int i=0; i<distant_nb_types; i++)
+ {
+ meshing->setConnectivity(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; i<distant_nb_elems; i++)
+ {
+ distant_ids_recv[i]=*recv_buffer_ptr++;
+ }
+
+ distant_mesh=meshing;
+ delete[] recv_buffer;
+ delete[] nbtypes;
+ }
+ if (nb_recv_floats >0)
+ delete[] recv_coords; // the coordinates are present if the recv_buffer is not empty
+
+}
+
+MEDMEM::MESH* ElementLocator::_meshFromElems(set<int>& 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<int> nodes;
+ int nbconn=0;
+ map<MED_EN::medGeometryElement,int> nbelems_per_type;
+ for (set<int>::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<int,int> big2small;
+ int i=0;
+ for (set<int>::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<int>::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<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
+ {
+ int dim = _local_mesh->getSpaceDimension();
+ for (int i=0; i<dim;i++)
+ {
+ new_coords[index]=coords[(*iter-1)*dim+i];
+ index++;
+ }
+ }
+
+ int* nbtypes=new int[nbelems_per_type.size()];
+ MED_EN::medGeometryElement* new_types = new MED_EN::medGeometryElement[nbelems_per_type.size()];
+ index=0;
+ for (map<MED_EN::medGeometryElement,int>::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; i<nbelems_per_type.size(); i++)
+ {
+ meshing->setConnectivity(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;
+}
+
+}
--- /dev/null
+#ifndef ELEMENTLOCATOR_HXX_
+#define ELEMENTLOCATOR_HXX_
+
+#include <vector>
+#include <set>
+
+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<MEDMEM::MESH*> _distant_meshes;
+ double* _domain_bounding_boxes;
+ const ProcessorGroup& _distant_group;
+ const ProcessorGroup& _local_group;
+ ProcessorGroup* _union_group;
+ std::vector<int> _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<int>& elems);
+};
+
+}
+
+#endif /*ELEMENTLOCATOR_HXX_*/
--- /dev/null
+#ifndef EXPLICIT_MAPPING_HXX_
+#define EXPLICIT_MAPPING_HXX_
+
+#include <vector>
+#include <map>
+#include <set>
+
+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<int,int> idistant) {
+ _mapping.push_back(idistant);
+ }
+ void setDistantElem(int ilocal, pair<int,int> idistant)
+ {
+ _mapping[ilocal]=idistant;
+ }
+
+ int nbDistantDomains()
+ {
+ if (_distant_domains.empty())
+ {
+ for (vector <pair<int,int> >::const_iterator iter= _mapping.begin();
+ iter!=_mapping.end();
+ iter++)
+ _distant_domains.insert(iter->first);
+ }
+ return _distant_domains.size();
+ }
+
+ pair <int,int> 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<int> 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; i<nbprocs; i++)
+ for (int ielem=0; ielem<sizes[i]; ielem++)
+ {
+ _mapping[indmap].first=i;
+ _mapping[indmap].second=commbuffer[indmap*2+1];
+ _buffer_index[indmap]=commbuffer[indmap*2+1];
+ indmap++;
+ }
+ _numbers=new int [nbtarget];
+ _domains=new int [nbtarget];
+
+ int index=0;
+ for (int i=0; i<nbtarget; i++)
+ {
+ if (sizes[targetrank[i]]>0)
+ {
+ _numbers[index]=sizes[targetrank[i]];
+ _domains[index]=i;
+ index++;
+ }
+ }
+ _send_counts=new int[nbprocs];
+ for (int i=0; i<nbprocs; i++)
+ _send_counts[i]=sizes[i];
+ }
+
+ int* getBufferIndex() const { return _buffer_index;}
+ int* getCounts() const{return _send_counts;}
+ private:
+ vector <pair<int,int> > _mapping;
+ set<int> _distant_domains;
+ int* _numbers;
+ int* _domains;
+ int* _commbuffer;
+ int* _buffer_index;
+ int* _send_counts;
+
+ void computeNumbers()
+ {
+ map <int,int> 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<int,int>::const_iterator iter=counts.begin();
+ iter!=counts.end();
+ iter++)
+ {
+ _numbers[counter]=iter->second;
+ _domains[counter]=iter->first;
+ counter++;
+ }
+ }
+ }
+ };
+
+}
+#endif
_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];
}
ExplicitTopology::~ExplicitTopology()
{
+ if (_loc2glob != 0) delete[] _loc2glob;
}
--- /dev/null
+#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<double> INTERPOLATION_2D::calcul_cos_et_sin(const double* P_1,const double* P_2,
+ const double* P_3)
+{
+
+ vector<double> 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<double>& V)
+{
+
+ int taille=V.size();
+ // bool isPresent=false;
+ for(int i=0;i<taille/2;i++)
+ {
+ double dx=P[0]-V[2*i];
+ dx=dx>0.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<double> INTERPOLATION_2D::bary_poly(const vector<double>& V)
+{
+ vector<double> Bary;
+ int taille=V.size();
+ double x=0;
+ double y=0;
+
+ for (int i=0;i<taille/2;i++)
+ {
+ x=x+V[2*i];
+ y=y+V[2*i+1];
+ }
+ double A=2*x/taille;
+ double B=2*y/taille;
+ Bary.push_back(A);//taille vecteur=2*nb de points.
+ Bary.push_back(B);
+
+
+ return Bary;
+}
+
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+/* calcul la surface d'un polygone. */
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+double INTERPOLATION_2D::Surf_Poly(const vector<double>& 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<double>& 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<double>& 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<double> 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<double> 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<Vect.size(); ++i)
+ //debug cout << Vect[i] << " ";
+ //debug cout << endl << endl;
+
+ return Vect;
+}
+
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+/* fonction pour reconstituer un polygone convexe à partir */
+/* d'un nuage de point. */
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+vector<double> INTERPOLATION_2D::reconstruct_polygon(const vector<double>& V)
+{
+
+ int taille=V.size();
+ if(taille<=6)
+ {return V;}
+ else
+ {
+ double COS[taille/2];
+ double SIN[taille/2];
+ double angle[taille/2];
+ vector<double> Bary=bary_poly(V);
+ COS[0]=1.0;
+ SIN[0]=0.0;
+ angle[0]=0.0;
+ for(int i=0; i<taille/2-1;i++)
+ {
+ vector<double> 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<double> Pt_ordonne;
+ Pt_ordonne.reserve(taille);
+ map<double,int> Ordre;
+ for(int i=0;i<taille/2;i++)
+ {Ordre[angle[i]]=i;}
+ map <double,int>::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<int> 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<int> 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<int> V_12;
+ V_12.reserve(2);
+ insert_iterator<vector<int> > ib_1(V_12,V_12.begin());
+ vector<int> V_23;
+ V_23.reserve(2);
+ insert_iterator<vector<int> > ib_2(V_23,V_23.begin());
+ vector<int> V_13;
+ V_13.reserve(2);
+ insert_iterator<vector<int> > 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_23.size();i++)
+ {
+ if(V_23[i]!=i_S)
+ {
+ tab.push_back(V_23[i]);
+ break;
+ }
+ }
+
+ for(int i=0;i<V_13.size();i++)
+ {
+ if(V_13[i]!=i_S)
+ {
+ tab.push_back(V_13[i]);
+ break;
+ }
+ }
+ return tab;
+}
+
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+/* 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<int>& V)
+{
+ int taille=V.size();
+ int A=0;
+ for(int i=0;i<taille;i++)
+ {
+ if(Num==V[i])
+ {
+ A=1;
+ break;
+ }
+ }
+ if(A==0)
+ {V.push_back(Num); }
+}
+
+
+
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+/*localise le barycentre des mailles de P dans le maillage S*/
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,monMaillageP& MaP,
+ monMaillageS& MaS,vector<int>& localis )
+{
+
+ //calcul des coordonnées des barycentre des mailles de P
+
+
+ //calcule des coordonnees du barycentre d'une cellule
+ std::auto_ptr<SUPPORT::SUPPORT> Support ( new SUPPORT(&myMesh_P,"monSupport",MED_CELL) );
+ std::auto_ptr<FIELD<double, FullInterlace> > 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<MaP._NbMaille;i++)
+ localis.push_back(mapping.Get_Mapping()[i]+1);
+
+
+ if (_SecondMethodOfLocalisation)
+ {
+ if(_PrintLevel)
+ cout << endl << "Option de localisations complémentaires" << endl;
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ //on vérifie que tout les barycentres ont obtenu un n° de maille
+ //non nul.
+ MEDMEM::INTERPOLATION<2> monInterpol_S(myMesh_S);
+ for(int i_P=0;i_P<MaP._NbMaille;i_P++)
+ {
+ if(localis[i_P]<=0)
+ {
+ //debug cout << endl << "------------------------------------------------" << endl << "Localise : barycentre maille " << i_P << " en dehors" << endl;
+ int Test=0;
+
+ int NP_1=MaP._Connect[3*i_P];
+ int NP_2=MaP._Connect[3*i_P+1];
+ int NP_3=MaP._Connect[3*i_P+2];
+
+ double Coord_bary_i_P[2]={Bary[2*i_P],Bary[2*i_P+1]};
+ int Vois_N=monInterpol_S.getNearestNode(Coord_bary_i_P);
+ int Ni=MaS._ReversNConnectIndex[Vois_N-1];
+ int Nf=MaS._ReversNConnectIndex[Vois_N];
+ int diff=Nf-Ni;
+
+ for(int j=0;j<diff;j++)
+ {
+ int N_Maille=MaS._ReversNConnect[Ni-1+j];
+ int NS_1=MaS._Connect[3*(N_Maille-1)];
+ int NS_2=MaS._Connect[3*(N_Maille-1)+1];
+ int NS_3=MaS._Connect[3*(N_Maille-1)+2];
+
+
+ //debug cout << "mailles sources testées : " << N_Maille << endl;
+ vector<double> 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<diff;j++)
+ {
+ int N_Maille=MaS._ReversNConnect[Ni-1+j];
+ int N_1=MaS._Connect[3*(N_Maille-1)];
+ int N_2=MaS._Connect[3*(N_Maille-1)+1];
+ int N_3=MaS._Connect[3*(N_Maille-1)+2];
+
+ //debug cout << "mailles sources testées : " << N_Maille << endl;
+ vector<double> 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<map<int,double> >& Surface_d_intersect,FIELD<double>& 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<int> maill_deja_vu_S ;
+
+ maill_deja_vu_S.push_back(i_S1);
+
+ for (int M_S=0; M_S<maill_deja_vu_S.size(); M_S++)
+ {
+ if( abs(Surf_i_P-Surface)/Surf_i_P>_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<double> 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<int> 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<int> 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<double> 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<int> 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<map<int,double> > 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"<<endl;}
+ string* Type_S = myMesh_S.getCellTypeNames(MED_CELL);
+ string* Type_P = myMesh_P.getCellTypeNames(MED_CELL);
+ if ( Type_S[0] != "MED_TRIA3" || Type_P[0] != "MED_TRIA3")
+ { cout<<"l'un au moins des deux maillages n'est pas triangulaires"<<endl;}
+ //VB
+ delete[] Type_S;
+ delete[] Type_P;
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ monMaillageS MaS;
+ monMaillageP MaP;
+
+ MaS._NbNoeud = myMesh_S.getNumberOfNodes() ;
+ MaP._NbNoeud = myMesh_P.getNumberOfNodes() ;
+ MaS._NbMaille =myMesh_S.getNumberOfElements(MED_CELL,MED_TRIA3);
+ MaP._NbMaille =myMesh_P.getNumberOfElements(MED_CELL,MED_TRIA3);
+
+ //on charge les connectivités des deux maillages.
+ MaS._Connect =myMesh_S.getConnectivity(MED_FULL_INTERLACE,
+ MED_NODAL, MED_CELL, MED_TRIA3) ;
+ MaP._Connect =myMesh_P.getConnectivity(MED_FULL_INTERLACE,
+ MED_NODAL, MED_CELL, MED_TRIA3) ;
+
+ //on charge les coordonnées des noeuds.
+ MaS._Coord = myMesh_S.getCoordinates(MED_FULL_INTERLACE);
+ MaP._Coord = myMesh_P.getCoordinates(MED_FULL_INTERLACE);
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ //on a besoin de recupérer connectivité inverse de S.
+ MaS._ReversNConnect=
+ myMesh_S.getReverseConnectivity(MED_NODAL);
+ MaS._ReversNConnectIndex=
+ myMesh_S.getReverseConnectivityIndex(MED_NODAL);
+
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+
+ int SpaceDimension_S = myMesh_S.getNumberOfNodes() ;
+ int SpaceDimension_P = myMesh_P.getNumberOfNodes() ;
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //on cherche la dimension caractéristique des maillages
+ vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+ vector<vector<double> > 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 :"<<endl;
+ cout<<" - dimension caractéristique du maillage source : "<<MaS._DimCaracteristic<<endl;
+ cout<<" - dimension caractéristique du maillage projeté: "<<MaS._DimCaracteristic<<endl;
+ cout<<" - d'où la dimension caractéristique: "<<_DimCaracteristic<<endl;
+ }
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //on cherche pour chaque maille i_P du maillage projetté
+ // une maille i_S du maillage S qui chevauche la maille i_P.
+
+ vector<int> localis;
+ localis.reserve(MaP._NbMaille);
+ MEDMEM::MESH* ptr_S = const_cast<MEDMEM::MESH*>(&myMesh_S);
+ MEDMEM::MESH* ptr_P = const_cast<MEDMEM::MESH*>(&myMesh_P);
+
+ cout << "INTERPOLATION_2D::local_iP_dansS"<<endl;
+ local_iP_dans_S(*ptr_S,*ptr_P,MaP,MaS,localis);
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //on crée un tableau où l'on rentrera la valeurs des intersections.
+ cout << "INTERPOLATION_2D::calcul intersec"<<endl;
+ map<int,double> MAP_init;
+ vector<map<int,double> > 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<SUPPORT> mySupport_P (new SUPPORT(const_cast<MEDMEM::MESH*>(&myMesh_P),"Support on all CELLS",MED_CELL) );
+ MEDMEM::FIELD<double> 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<<vector<pair<int,double> > >& 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_P<MaP._NbMaille ; i_P++)
+ {
+ int i_S_init=localis[i_P];
+ if(i_S_init>0)
+ rempli_vect_d_intersect(i_P+1,i_S_init,MaP,MaS,Surface_d_intersect,myField_P);
+
+ }
+
+ if (_PrintLevel >= 2)
+ {
+ cout<<endl<<"Impression des surfaces d'intersection:"<<endl<<endl;
+ cout<<"(maille source, maille cible):surface d'intersection"<<endl;
+ for(int i=0; i< Surface_d_intersect.size();i++)
+ {
+ if(Surface_d_intersect[i].size()!=0)
+ {
+ map<int,double>::iterator surface;
+ for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++)
+ cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+ }
+ }
+ }
+
+// /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+// //On sauvegarde le champ crée sur la maillage P.
+// if (_PrintLevel>=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<<endl<<"Pour chaque mailles de " <<myMesh_S.getName();
+//// cout<<"Impression de la fraction surfaciques totale dans le maillage : "<<myMesh_P.getName()<<endl;
+//// for(int i=0;i<MaP._NbMaille;i++)
+//// cout<<"maille n°"<<i+1<<" de "<<myMesh_P.getName()<<": "<<myField_P.getValueIJ(i+1,1)<<endl;
+// }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ return Surface_d_intersect;
+
+}
+
+
+
+
+
+/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+/* 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<double>& 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);
+}
--- /dev/null
+#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 <cmath>
+#include <map>
+#include <sstream>
+#include <string>
+#include <set>
+#include <algorithm>
+#include <vector>
+#include <complex>
+
+
+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<double> 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<double>& V);
+
+ //calcul les coordonnées du barycentre
+ //d'un polygone
+ vector<double> bary_poly(const vector<double>& V);
+
+ //fonction qui calcul la surface d'un polygone
+ double Surf_Poly(const vector<double>& 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<double>& 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<double>& V);
+
+ //calcul l'intersection de deux triangles (sommets du triangle donné dans
+ //ordre quelconque)
+ vector<double> 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<double> reconstruct_polygon(const vector<double>& V);
+
+ //fonction pour trouver les mailles voisines d'une maille triangle.
+ vector<int> 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<int>& 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<int>& 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<map<int,double> >& Surface_d_intersect,
+ FIELD<double>& myField_P);
+
+
+ //fonction principale pour interpoler deux maillages triangulaires.
+ std::vector<map<int, double> > 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<double>& FractSurf,MEDMEM::MESH* myMesh_P,
+ MEDMEM::MESH* myMesh_S,string FileName,string NameMesh_P,string NameMesh_S);
+
+
+
+ // MEDMEM::FIELD<double>* createField();
+ // ...
+
+ private:
+
+
+};
+
+#endif
--- /dev/null
+#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<map<int,double> > 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<double>* 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<int,double>::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<pair<int,int> >::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_distant<distant_support->getNumberOfElements(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<dim; idim++)
+// coords[idim+dim*i]=coords_source[node*dim+idim];
+// }
+// double* moments = new double[dim+1];
+// intersection (target_support, ielem_target, dim, coords,intersection_moments);
+// }
+//
+// //intersection_moments :
+// // 0 -> 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<dim; i++)
+// _coeffs.push_back(intersection_moments[i+1]);
+// _mapping.addElementFromSource(ielem,iproc_distant,ielem_target);
+// }
+// }
+// irow++;
+// }
+//}
+
+void InterpolationMatrix::prepare()
+{
+ int nbelems = _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ for (int ielem=0; ielem < nbelems; ielem++)
+ {
+ _row_offsets[ielem+1]+=_row_offsets[ielem];
+ }
+ _mapping.prepareSendRecv();
+}
+
+void InterpolationMatrix::multiply(MEDMEM::FIELD<double>& 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; irow<nbrows; irow++)
+ {
+ for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
+ {
+ for (int icomp=0; icomp<field.getNumberOfComponents(); icomp++)
+ {
+ double coeff_row = field.getValueIJ(irow+1,icomp+1);
+ target_value[_col_numbers[icol]-1]+=_coeffs[icol]*coeff_row;
+ }
+ }
+ }
+ _mapping.sendRecv(target_value,field);
+
+ if (_col_offsets.size()>0)
+ delete[] target_value;
+}
+
+
+}
--- /dev/null
+#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<double>&) const;
+ void prepare();
+ int getNbRows() const {return _row_offsets.size();}
+
+ private:
+ // vector<pair <int,int> > _source_indices;
+ vector<int> _row_offsets;
+ vector<int> _col_numbers;
+ vector<pair<int,int> > _col_offsets;
+ vector<double> _coeffs;
+ const MEDMEM::MESH& _source_support;
+ MxN_Mapping _mapping;
+ string _method;
+ const ProcessorGroup& _local_group;
+ const ProcessorGroup& _distant_group;
+
+
+ };
+
+}
+
+#endif /*INTERPOLATIONMATRIX_HXX_*/
--- /dev/null
+#include <mpi.h>
+#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 "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<endl;
+ _interpolation_matrix->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<<endl;
+ if (distant_mesh!=0)
+ {
+ delete distant_mesh;
+ delete[] distant_ids;
+ distant_mesh=0;
+ distant_ids=0;
+ }
+ }
+ }
+ _interpolation_matrix->prepare();
+}
+
+
+void IntersectionDEC::recvData()
+{
+ _interpolation_matrix->multiply(*_target_field->getField());
+}
+
+void IntersectionDEC::sendData()
+{
+ _interpolation_matrix->multiply(*_source_field->getField());
+}
+
+}
+
+
+
--- /dev/null
+#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
// copying proc_ids in ranks
copy<set<int>::const_iterator,int*> (proc_ids.begin(), proc_ids.end(), ranks);
- for (int i=0; i<proc_ids.size();i++)
- cout << "proc id " << ranks[i]<<endl;
_comm_interface.groupIncl(group_world, proc_ids.size(), ranks, &_group);
return new MPIProcessorGroup(_comm_interface, procs);
}
+ProcessorGroup* MPIProcessorGroup::fuse (const ProcessorGroup& group) const
+{
+ set <int> procs = _proc_ids;
+ const set<int>& distant_proc_ids = group.getProcIDs();
+ for (set<int>::const_iterator iter=distant_proc_ids.begin(); iter!=distant_proc_ids.end(); iter++)
+ {
+ procs.insert(*iter);
+ }
+ return new MPIProcessorGroup(_comm_interface,procs);
+}
+
}
MPIProcessorGroup(const CommInterface& interface, set<int> proc_ids);
MPIProcessorGroup (const ProcessorGroup& proc_group, set<int> 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);}
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
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
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)
--- /dev/null
+#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<MPIProcessorGroup*>(_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<int> 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<double>& field)
+{
+ CommInterface comm_interface=_union_group->getCommInterface();
+ const MPIProcessorGroup* group = static_cast<const MPIProcessorGroup*>(_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<int> offsets = _send_proc_offsets;
+ for (int i=0; i<_sending_ids.size();i++)
+ {
+ int iproc = _sending_ids[i].first;
+ for (int icomp=0; icomp<nbcomp; icomp++)
+ sendbuf[offsets[iproc]*nbcomp+icomp]=field.getValueIJ(i+1, icomp+1);
+ offsets[iproc]++;
+ }
+
+ //communication phase
+ const MPI_Comm* comm = group->getComm();
+ 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<nbcomp; icomp++){
+ field.setValueIJ(_recv_ids[i],icomp+1,*recvptr);
+ recvptr++;
+ }
+ }
+
+ if (sendbuf!=0) delete[] sendbuf;
+ if (recvbuf!=0) delete[] recvbuf;
+ 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(double* sendfield, MEDMEM::FIELD<double>& field) const
+{
+ CommInterface comm_interface=_union_group->getCommInterface();
+ const MPIProcessorGroup* group = static_cast<const MPIProcessorGroup*>(_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<int> offsets = _send_proc_offsets;
+ for (int i=0; i<_sending_ids.size();i++)
+ {
+ int iproc = _sending_ids[i].first;
+ for (int icomp=0; icomp<nbcomp; icomp++)
+ // sendbuf[offsets[iproc]*nbcomp+icomp]=sendfield[i+1*nbcomp+icomp];
+ sendbuf[offsets[iproc]*nbcomp+icomp]=sendfield[i*nbcomp+icomp];
+ offsets[iproc]++;
+ }
+
+ //communication phase
+ const MPI_Comm* comm = group->getComm();
+ 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; j<field.getNumberOfComponents(); j++)
+// field.setValueIJ(i+1,j+1,0.0);
+
+ //setting the received values in the field
+ double* recvptr=recvbuf;
+ for (int i=0; i< _recv_proc_offsets[_union_group->size()]; i++)
+ {
+ for (int icomp=0; icomp<nbcomp; icomp++){
+ double temp = field.getValueIJ(_recv_ids[i],icomp+1);
+ field.setValueIJ(_recv_ids[i],icomp+1,temp+*recvptr);
+ recvptr++;
+ }
+ }
+
+ if (sendbuf!=0) delete[] sendbuf;
+ if (recvbuf!=0) delete[] recvbuf;
+ delete[] sendcounts;
+ delete[] recvcounts;
+ delete[] senddispls;
+ delete[] recvdispls;
+
+}
+
+
+}
--- /dev/null
+#ifndef MXN_MAPPING_HXX_
+#define MXN_MAPPING_HXX_
+
+#include <vector>
+
+#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<double>& field);
+ void sendRecv(double* field, MEDMEM::FIELD<double>& field) const ;
+
+private :
+// ProcessorGroup& _local_group;
+// ProcessorGroup& _distant_group;
+ ProcessorGroup* _union_group;
+ int _nb_comps;
+ std::vector<pair<int,int> > _sending_ids;
+ std::vector<int> _recv_ids;
+ std::vector<int> _send_proc_offsets;
+ std::vector<int> _recv_proc_offsets;
+
+};
+
+}
+
+#endif /*MXN_MAPPING_HXX_*/
--- /dev/null
+#include <mpi.h>
+#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 <fvm_parall.h>
+#include <fvm_nodal.h>
+#include <fvm_nodal_append.h>
+#include <fvm_locator.h>
+}
+
+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; itype<nbtypes; itype++)
+ {
+ fvm_element_t fvm_type;
+ switch (types[itype])
+ {
+ case MED_EN::MED_TRIA3 :
+ fvm_type=FVM_FACE_TRIA;
+ break;
+ case MED_EN::MED_QUAD4 :
+ fvm_type=FVM_FACE_QUAD;
+ break;
+ case MED_EN::MED_TETRA4 :
+ fvm_type=FVM_CELL_TETRA;
+ break;
+ case MED_EN::MED_HEXA8 :
+ fvm_type=FVM_CELL_HEXA;
+ break;
+ default:
+ throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
+ break;
+
+ }
+ fvm_lnum_t nbelems = mesh->getNumberOfElements(MED_EN::MED_CELL, types[itype]);
+ fvm_lnum_t* conn = const_cast<fvm_lnum_t*> (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; itype<nbtypes; itype++)
+ {
+ fvm_element_t fvm_type;
+ switch (types[itype])
+ {
+ case MED_EN::MED_TRIA3 :
+ fvm_type=FVM_FACE_TRIA;
+ break;
+ case MED_EN::MED_QUAD4 :
+ fvm_type=FVM_FACE_QUAD;
+ break;
+ case MED_EN::MED_TETRA4 :
+ fvm_type=FVM_CELL_TETRA;
+ break;
+ case MED_EN::MED_HEXA8 :
+ fvm_type=FVM_CELL_HEXA;
+ break;
+ default:
+ throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
+ break;
+
+ }
+ fvm_lnum_t nbelems = support->getNumberOfElements(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<fvm_lnum_t*> (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<fvm_lnum_t*> (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<const MPIProcessorGroup*> (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<const MPIProcessorGroup*> (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 <nbcomp; j++)
+ distant_values[i*nbcomp+j]=values[_distant_locations[i]*nbcomp+j];
+
+ fvm_locator_exchange_point_var(_locator,
+ distant_values,
+ 0,
+ sizeof(double),
+ nbcomp);
+}
+
+}
+
+
+
--- /dev/null
+#ifndef NONCOINCIDENTDEC_HXX_
+#define NONCOINCIDENTDEC_HXX_
+
+#include <fvm_locator.h>
+
+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
#include "ComponentTopology.hxx"
#include "ParaSUPPORT.hxx"
#include "StructuredParaSUPPORT.hxx"
+#include "UnstructuredParaSUPPORT.hxx"
#include "ExplicitCoincidentDEC.hxx"
#include "StructuredCoincidentDEC.hxx"
+
#include "ParaFIELD.hxx"
#include "ParaMESH.hxx"
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<const StructuredParaSUPPORT*>(para_support)!=0)
{const BlockTopology* source_topo = dynamic_cast<const BlockTopology*>(para_support->getTopology());
_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<double>* 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<const ExplicitTopology*> (_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<BlockTopology*> (_topology);
namespace MEDMEM{
class MEDEXCEPTION;
+ template <class T> class FIELD<T>;
}
{
class ComponentTopology;
class ParaSUPPORT;
+class ProcessorGroup;
class ParaFIELD
{
ParaFIELD(MEDMEM::driverTypes driver_type, const string& file_name,
const string& driver_name, const ComponentTopology& component_topology)
throw (MEDMEM::MEDEXCEPTION);
+ ParaFIELD(MEDMEM::FIELD<double>* 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<double>* getField() const {return _field;}
+ const ParaSUPPORT* getSupport() const {return _support;}
Topology* getTopology() const {return _topology;}
int nbComponents() const {return _component_topology.nbComponents();}
private:
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)
{
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;
///////////////////////////////////////////
char joint_description[MED_TAILLE_DESC];
char name[MED_TAILLE_NOM];
char name_distant[MED_TAILLE_NOM];
- cout << "arguments"<< fid<<" "<<file<<" "<<ijoint<<" "<<name<<" "<<joint_description<<" "<<distant<<" "<<name_distant<<endl;
+ // cout << "arguments"<< fid<<" "<<file<<" "<<ijoint<<" "<<name<<" "<<joint_description<<" "<<distant<<" "<<name_distant<<endl;
int ncorr = med_2_2::MEDjointInfo(fid,meshname, ijoint, name,
joint_description,
&distant, name_distant);
END_OF("MEDSPLITTER::MESHCollectionDriver::read()")
};
+/*! Constructor for creating a ParaMESH from a local mesh and
+ * a processor group. Constructor must be called by all the processors
+ * in the group. */
+
+ParaMESH::ParaMESH(MEDMEM::MESH& subdomain_mesh, const ProcessorGroup& proc_group, const string& name):
+_mesh(&subdomain_mesh),
+_name (name),
+_my_domain_id(proc_group.myRank()),
+_block_topology (new BlockTopology(proc_group, subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS)))
+{
+ _cellglobal = new int[subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS)];
+ int offset = _block_topology->localToGlobal(make_pair(_my_domain_id,0));
+ for (int i=0; i<subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS); i++)
+ {
+ _cellglobal[i]=offset+i;
+ }
+}
-ParaMESH::~ParaMESH(){if (_mesh !=0) delete _mesh;};
+ParaMESH::~ParaMESH(){if (_mesh !=0) delete _mesh;};
+
/*! method for writing a distributed MESH
*
* \param driverType type of driver used (MED_DRIVER,VTK_DRIVER)
return _edgeglobal;
case MED_NODE:
return _nodeglobal;
+ default :
+ return 0;
}
}
ParaMESH(MEDMEM::driverTypes driver_type, const std::string& file_name,
const ProcessorGroup& group)
throw (MEDMEM::MEDEXCEPTION);
+ ParaMESH(MEDMEM::MESH& subdomain_mesh, const ProcessorGroup& proc_group, const string& name);
void write(MEDMEM::driverTypes driverType, const std::string& fileName="")
throw (MEDMEM::MEDEXCEPTION);
virtual ~ParaMESH();
MEDMEM::MESH* getMesh() const {return _mesh;}
ParaMEDMEM::BlockTopology* getBlockTopology()const {return _block_topology;}
- const string& getFilename() const {return _medfilename;}
+// const string& getFilename() const {return _medfilename;}
const int* getGlobalNumbering(MED_EN::medEntityMesh)const;
private:
//mesh object underlying the ParaMESH object
int _my_domain_id;
//global topology of the cells
ParaMEDMEM::BlockTopology* _block_topology;
- //name of the local filename (!= masterfilename)
- string _medfilename;
// pointers to global numberings
int* _nodeglobal;
int* _edgeglobal;
{
}
- ParaSUPPORT::ParaSUPPORT(const MEDMEM::SUPPORT& support):_support(&support) {}
+ ParaSUPPORT::ParaSUPPORT(const MEDMEM::SUPPORT& support, const ProcessorGroup& proc_group):_support(&support) {
+ _mesh = new ParaMESH(*(support.getMesh()), proc_group, "mesh from support");
+ }
ParaSUPPORT::~ParaSUPPORT()
{
{
class Topology;
class ParaMESH;
+ class ProcessorGroup;
class ParaSUPPORT
{
public:
ParaSUPPORT();
ParaSUPPORT(const ParaMESH* mesh, const MEDMEM::SUPPORT* support):
_support(support), _mesh(mesh){};
- ParaSUPPORT(const MEDMEM::SUPPORT&);
+ ParaSUPPORT(const MEDMEM::SUPPORT&, const ProcessorGroup& proc_group);
virtual ~ParaSUPPORT();
- virtual const Topology* getTopology() const {};
+ virtual const Topology* getTopology() const{};
virtual const MEDMEM::SUPPORT* getSupport() const {return _support;}
virtual const ParaMESH* getMesh() const {return _mesh;}
protected :
ProcessorGroup (const ProcessorGroup& proc_group, std::set<int> 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;
virtual int translateRank(const ProcessorGroup*, int) const =0;
virtual ProcessorGroup* createComplementProcGroup() const =0;
virtual ProcessorGroup* createProcGroup() const=0;
+ virtual const std::set<int>& getProcIDs()const {return _proc_ids;}
protected:
const CommInterface& _comm_interface;
std::set<int> _proc_ids;
#include "Topology.hxx"
#include "ExplicitTopology.hxx"
+#include "BlockTopology.hxx"
#include "ParaMESH.hxx"
#include "UnstructuredParaSUPPORT.hxx"
#include "MEDMEM_Support.hxx"
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):
_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
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;
};
--- /dev/null
+#include <string>
+#include <vector>
+#include <map>
+#include <iostream>
+#include <mpi.h>
+
+#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 :"<<endl;
+ cout << "mpirun -n <nbprocs> test_NonCoincidentDEC <nbprocsource>"<<endl;
+ cout << " (nbprocs >=2)"<<endl;
+ return 1;
+ }
+ int nproc_source = atoi(argv[1]);
+ set<int> self_procs;
+ set<int> procs_source;
+ set<int> procs_target;
+
+ for (int i=0; i<nproc_source; i++)
+ procs_source.insert(i);
+ for (int i=nproc_source; i<size; i++)
+ procs_target.insert(i);
+ self_procs.insert(rank);
+
+ ParaMEDMEM::CommInterface interface;
+
+ ParaMEDMEM::ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
+ ParaMEDMEM::ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
+ ParaMEDMEM::ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
+
+ ParaMEDMEM::ParaMESH* source_mesh=0;
+ ParaMEDMEM::ParaMESH* target_mesh=0;
+ ParaMEDMEM::ParaMESH* source_mesh_bis = 0;
+ ParaMEDMEM::ParaMESH* target_mesh_bis = 0;
+
+ Topology* topo_source;
+ Topology* topo_target;
+
+ MEDMEM::SUPPORT* boundary_group;
+ //loading the geometry for the source group
+ if (source_group->containsMyRank())
+ {
+ string master = "/home/vb144235/resources/square128000_split";
+ cout <<"loading source"<<endl;
+ source_mesh=new ParaMESH(MED_DRIVER,master,*source_group);
+ cout <<"end of load"<<endl;
+ topo_source=source_mesh->getBlockTopology();
+ ostringstream strstream;
+ strstream <<master<<rank+1<<".med";
+ //strstream<<"/home/vb144235/resources/square128000.med";
+ ostringstream meshname ;
+ meshname<< "Mesh_1_"<< rank+1;
+ // meshname<<"Mesh_1";
+ MEDMEM::MESH* mesh_bis = new MESH(MED_DRIVER,strstream.str(),meshname.str());
+ source_mesh_bis = new ParaMESH(*mesh_bis,*source_group,"mesh_from_mem");
+ // const vector <MEDMEM::GROUP*> mesh_groups;
+ //
+ // for (int i=0; i<mesh_groups.size(); i++)
+ // {
+ // if (mesh_groups[i]->getName() == "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_"<<rank-nproc_source+1;
+ //meshname<<"Mesh_1";
+
+ MEDMEM::MESH* mesh_bis = new MESH(MED_DRIVER,strstream.str(),meshname.str());
+ target_mesh_bis = new ParaMESH(*mesh_bis,*target_group, "mesh_from_mem");
+ // const vector <MEDMEM::GROUP*> mesh_groups;
+ //
+ // for (int i=0; i<mesh_groups.size(); i++)
+ // {
+ // if (mesh_groups[i]->getName() == "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"<<endl;
+ ParaFIELD source_field (&source_support,source_comp);
+
+ int nb_local = source_field.getTopology()->getNbLocalElements();
+ cout << "Source field nb elems on rank : "<<rank<<" : "<<nb_local<<endl;
+ double * value= new double[nb_local];
+ for(int ielem=0; ielem<nb_local;ielem++)
+ value[ielem]=1.0;
+ source_field.getField()->setValue(value);
+ cout <<"creating intersectionDEC"<<endl;
+ IntersectionDEC dec(*source_group,*target_group);
+ dec.attachSourceField(&source_field);
+ dec.synchronize();
+ cout<<"DEC usage"<<endl;
+ dec.sendData();
+ cout <<"writing"<<endl;
+ source_mesh->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; ielem<nb_local;ielem++)
+ value[ielem]=0.0;
+ target_field.getField()->setValue(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;
+}
+
+
--- /dev/null
+#include <string>
+#include <vector>
+#include <map>
+#include <iostream>
+#include <mpi.h>
+
+#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 :"<<endl;
+ cout << "mpirun -n <nbprocs> test_NonCoincidentDEC <nbprocsource>"<<endl;
+ cout << " (nbprocs >=2)"<<endl;
+ return 1;
+ }
+ int nproc_source = atoi(argv[1]);
+ set<int> self_procs;
+ set<int> procs_source;
+ set<int> procs_target;
+
+ for (int i=0; i<nproc_source; i++)
+ procs_source.insert(i);
+ for (int i=nproc_source; i<rank; i++)
+ procs_target.insert(i);
+ self_procs.insert(rank);
+
+ ParaMEDMEM::CommInterface interface;
+
+ ParaMEDMEM::ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
+ ParaMEDMEM::ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
+ ParaMEDMEM::ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
+
+ ParaMEDMEM::ParaMESH* source_mesh=0;
+ ParaMEDMEM::ParaMESH* target_mesh=0;
+ ParaMEDMEM::ParaMESH* source_mesh_bis = 0;
+ ParaMEDMEM::ParaMESH* target_mesh_bis = 0;
+
+ Topology* topo_source;
+ Topology* topo_target;
+
+ //loading the geometry for the source group
+ if (source_group->containsMyRank())
+ {
+ 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 << "_"<<rank<<endl;
+ MEDMEM::MESH mesh_bis = new MESH(MED_DRIVER,strstream.str(),'fluid');
+ source_mesh_bis = new ParaMESH(*mesh_bis,*source_group);
+ const vector <MEDMEM::GROUP*> mesh_groups;
+ const MEDMEM::GROUP* boundary_group;
+ for (int i=0; i<mesh_groups.size(); i++)
+ {
+ if (mesh_groups[i].getName() == "boundary")
+ boundary_group = mesh_groups[i];
+ }
+ }
+
+ //loading the geometry for the target group
+ if (target_group->containsMyRank())
+ {
+ 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)<<endl;
+ MEDMEM::MESH mesh_bis = new MESH(MED_DRIVER,strstream.str(),'Fuse_1');
+ target_mesh_bis = new ParaMESH(*mesh_bis,*target_group);
+ const vector <MEDMEM::GROUP*> mesh_groups;
+ const MEDMEM::GROUP* boundary_group;
+ for (int i=0; i<mesh_groups.size(); i++)
+ {
+ if (mesh_groups[i].getName() == "boundary")
+ boundary_group = mesh_groups[i];
+ }
+ }
+
+ UnstructuredParaSUPPORT* target_support;
+ UnstructuredParaSUPPORT* source_support;
+ ComponentTopology* target_comp;
+ ComponentTopology* source_comp;
+ ParaFIELD* target_field=0;
+ ParaFIELD* source_field=0;
+
+ //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);
+ ParaFIELD source_field (source_support,*source_comp);
+
+ int nb_local = source_field.getTopology()->getNbLocalElements();
+ cout << "Source field nb elems on rank : "<<rank<<" : "<<nb_local<<endl;
+ double * value= new double[nb_local];
+ for(int ielem=0; ielem<nb_local;ielem++)
+ value[ielem]=(double)ielem;
+ source_field.getField()->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;
+}
+
+
+
+
