+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
#include <mpi.h>
#include "CommInterface.hxx"
#include "ElementLocator.hxx"
#include <set>
#include <limits>
using namespace std;
+using namespace MED_EN;
namespace ParaMEDMEM
delete [] _domain_bounding_boxes;
}
+void ElementLocator::_getElemsIntesectingBndBox(const int *conn, const int *conn_index,
+ int nbElems, int firstElemId, double* distant_bb,
+ set<int>& elems)
+{
+ int dim=_local_mesh->getSpaceDimension();
+ double* elem_bb=new double[2*dim];
+ double* empty_bb=new double[2*dim];
+
+ for (int i=0; i<dim; i++)
+ {
+ empty_bb[i*2]=std::numeric_limits<double>::max();
+ empty_bb[i*2+1]=-std::numeric_limits<double>::max();
+ }
+ const double* coords = _local_mesh->getCoordinates(MED_FULL_INTERLACE);
+ const size_t size = size_t( sizeof(double)*2*dim );
+ for (int ielem=0; ielem<nbElems; ielem++)
+ {
+ memcpy( elem_bb, empty_bb, size );
+ 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(elems.end(), ielem+firstElemId);
+ }
+ }
+ delete[] elem_bb;
+ delete[] empty_bb;
+}
/*! Procedure for exchanging mesh between a distant proc and a local processor
\param idistantrank proc id on distant group
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]=std::numeric_limits<double>::max();
- elem_bb[i*2+1]=-std::numeric_limits<double>::max();
- }
- 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);
- }
- }
-
- //send_mesh contains null pointer if elems is empty
- 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=0;
- if (elems.size()>0)
- {
- 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++;
- }
- }
+ 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;
+
+ //defining pointers to med
+ const int *conn, *conn_index;
+ // standard elements
+ // ------------------
+ int nbelem = _local_mesh->getNumberOfElements(MED_CELL,
+ MED_ALL_ELEMENTS);
+ if ( nbelem > 0)
+ {
+ conn=_local_mesh->getConnectivity(MED_FULL_INTERLACE,
+ MED_NODAL,
+ MED_CELL,
+ MED_ALL_ELEMENTS);
+ conn_index= _local_mesh->getConnectivityIndex(
+ MED_NODAL,
+ MED_CELL);
+ _getElemsIntesectingBndBox( conn, conn_index, nbelem, 1, distant_bb, elems );
+ }
+ // polygons
+ // ---------
+ int nbPolygons = _local_mesh->getNumberOfPolygons(MED_CELL);
+ if ( nbPolygons > 0 )
+ {
+ conn=_local_mesh->getPolygonsConnectivity(MED_NODAL,
+ MED_CELL);
+ conn_index= _local_mesh->getPolygonsConnectivityIndex(
+ MED_NODAL,
+ MED_CELL);
+ _getElemsIntesectingBndBox( conn, conn_index, nbPolygons, nbelem+1, distant_bb, elems );
+ }
+ // polyhedrons
+ // ------------
+ int nbPolyhedrons = _local_mesh->getNumberOfPolyhedron();
+ if ( nbPolyhedrons > 0 )
+ {
+ conn=_local_mesh->getPolyhedronConnectivity(MED_NODAL);
+ conn_index= _local_mesh->getPolyhedronIndex(MED_NODAL);
+ const int* face_index = _local_mesh->getPolyhedronFacesIndex();
+ // make polyhedrons index
+ vector<int> index( nbPolyhedrons+1 );
+ index[0] = 1;
+ for (int ielem=1; ielem<=nbPolyhedrons; ielem++) {
+ index[ ielem ] = face_index[ conn_index[ ielem ]-1 ];
+ }
+ _getElemsIntesectingBndBox( conn, &index[0], nbPolyhedrons, nbelem+1, distant_bb, elems );
+ }
+
+ //send_mesh contains null pointer if elems is empty
+ 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=0;
+ if (elems.size()>0)
+ {
+ 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;
}
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);
+ const double* coords = _local_mesh->getCoordinates(MED_FULL_INTERLACE);
int nbnodes = _local_mesh->getNumberOfNodes();
double * minmax=new double [2*dim];
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];
-
- //treatment for non-empty mesh
- int nbtypes=0;
- int nbconn=0;
- int nbelems=0;
+ const int nbint = 9;
+ int* send_buffer = new int[nbint];
+ int* recv_buffer = new int[nbint];
+
+ //treatment for non-empty mesh
+ int nbtypes =0;
+ int nbconn =0;
+ int nbelems =0;
+ int nbpolys =0;
+ int nbpolyface=0;// nb of polyhedra faces
+ int nbpolyconn=0;
if (local_mesh !=0)
{
- nbtypes= local_mesh->getNumberOfTypes(MED_EN::MED_CELL);
- nbconn = local_mesh->getConnectivityLength(MED_EN::MED_FULL_INTERLACE, MED_EN::MED_NODAL,MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
- nbelems = local_mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
-
+ nbtypes = local_mesh->getNumberOfTypes(MED_CELL);
+ nbelems = local_mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
+ if ( nbelems )
+ nbconn = local_mesh->getConnectivityLength(MED_FULL_INTERLACE, MED_NODAL,MED_CELL,
+ MED_ALL_ELEMENTS);
+ nbpolys = local_mesh->getNumberOfPolygons(MED_CELL);
+ if ( nbpolys > 0 ) {
+ nbpolyface = 0;
+ nbpolyconn = local_mesh->getPolygonsConnectivityLength(MED_NODAL, MED_CELL);
+ } else {
+ nbpolys = local_mesh->getNumberOfPolyhedron();
+ if ( nbpolys > 0 ) {
+ nbpolyface = local_mesh->getNumberOfPolyhedronFaces();
+ nbpolyconn = local_mesh->getPolyhedronConnectivityLength(MED_NODAL);
+ }
+ }
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;
+ send_buffer[6] = nbpolys;
+ send_buffer[7] = nbpolyface;
+ send_buffer[8] = nbpolyconn;
}
else
{
- for (int i=0; i<6; i++)
+ for (int i=0; i<nbint; 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];
-
+ 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, nbint, MPI_INT, iprocdistant_in_union, 1112,
+ recv_buffer, nbint, 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];
+ int distant_nbpolys = recv_buffer[6];
+ int distant_nbpolyface= recv_buffer[7];
+ int distant_nbpolyconn= recv_buffer[8];
+
delete[] send_buffer;
delete[] recv_buffer;
-
+
//Second stage : exchanging connectivity buffers
- int nb_integers = nbtypes*2+nbconn+1+nbelems;
+ int nb_integers = 2*nbtypes+1+nbconn+nbelems;
+ if ( nbpolys )
+ nb_integers += nbpolys+nbpolyconn+nbpolys+1+nbpolyface+1;
send_buffer = new int[nb_integers];
- const MED_EN::medGeometryElement* types=0;
- const int* conn = 0;
- const int* global_numbering=0;
- int* ptr_buffer=send_buffer;
- if (local_mesh!=0)
- {
- conn=local_mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE, MED_EN::MED_NODAL,
- MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
-
- global_numbering = local_mesh->getGlobalNumberingIndex(MED_EN::MED_CELL);
- types = local_mesh->getTypes(MED_EN::MED_CELL);
-
- //copying the data in the integer buffer
-
- 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];
-
- //exchanging integer buffer
- 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);
-
+ const medGeometryElement* types= 0;
+ const int* global_numbering = 0;
+ const int* conn = 0;
+ const int* poly_conn = 0;
+ const int* poly_index = 0;
+ const int* face_index = 0;
+ int* ptr_buffer = send_buffer;
+ if (local_mesh!=0)
+ {
+ if ( nbelems ) {
+ conn=local_mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL,
+ MED_CELL, MED_ALL_ELEMENTS);
+
+ global_numbering = local_mesh->getGlobalNumberingIndex(MED_CELL);
+ types = local_mesh->getTypes(MED_CELL);
+ }
+ if ( nbpolys )
+ if ( nbpolyface ) {
+ poly_conn = local_mesh->getPolyhedronConnectivity(MED_NODAL);
+ poly_index = local_mesh->getPolyhedronIndex(MED_NODAL);
+ face_index = local_mesh->getPolyhedronFacesIndex();
+ } else {
+ poly_conn = local_mesh->getPolygonsConnectivity(MED_NODAL,MED_CELL);
+ poly_index = local_mesh->getPolygonsConnectivityIndex(MED_NODAL,MED_CELL);
+ }
+ //copying the data in the integer buffer
+
+ if ( nbelems ) {
+ 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+nbpolys)*sizeof(int));
+ if ( poly_conn ) {
+ ptr_buffer+=(nbelems+nbpolys);
+ memcpy(ptr_buffer, poly_conn, nbpolyconn*sizeof(int));
+ ptr_buffer+=nbpolyconn;
+ memcpy(ptr_buffer, poly_index, (nbpolys+1)*sizeof(int));
+ ptr_buffer+=(nbpolys+1);
+ if ( nbpolyface )
+ memcpy(ptr_buffer, face_index, (nbpolyface+1)*sizeof(int));
+ }
+ }
+ //preparing the recv buffers
+ int nb_recv_integers = 2*distant_nb_types+1+distant_nb_conn+distant_nb_elems;
+ if ( distant_nbpolys > 0 )
+ nb_recv_integers += distant_nbpolys+distant_nbpolyconn+distant_nbpolys+1+distant_nbpolyface+1;
+ recv_buffer=new int[nb_recv_integers];
+
+ //exchanging integer buffer
+ 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=0;
- double* coords=0;
- double* recv_coords=0;
-
- if (local_mesh!=0)
- {
- nb_send_floats = local_mesh->getSpaceDimension()* local_mesh->getNumberOfNodes();
- coords = const_cast<double*> (local_mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE));
- }
-
- if (nb_recv_floats>0)
- recv_coords = new double[nb_recv_floats];
+ int nb_send_floats =0;
+ double* coords =0;
+ double* recv_coords=0;
+
+ if (local_mesh!=0)
+ {
+ nb_send_floats = local_mesh->getSpaceDimension()* local_mesh->getNumberOfNodes();
+ coords = const_cast<double*> (local_mesh->getCoordinates(MED_FULL_INTERLACE));
+ }
+
+ if (nb_recv_floats>0)
+ recv_coords = new double[nb_recv_floats];
+
+ 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);
- 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 && distant_space_dim !=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);
+ meshing->setCoordinates(distant_space_dim, distant_nnodes, recv_coords, "CARTESIAN", MED_FULL_INTERLACE);
+ meshing->setNumberOfTypes(distant_nb_types,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++)
+ if ( distant_nb_elems )
{
- meshing->setConnectivity(recv_buffer_ptr, MED_EN::MED_CELL,recv_buffer[i]);
- recv_buffer_ptr+=nbtypes[i]*(recv_buffer[i]%100);
+ medGeometryElement* types_vector = new medGeometryElement[distant_nb_types];
+ for (int i=0; i<distant_nb_types; i++)
+ types_vector[i]=(medGeometryElement)recv_buffer_ptr[i];
+ meshing->setTypes(types_vector, 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_CELL);
+
+ for (int i=0; i<distant_nb_types; i++)
+ {
+ meshing->setConnectivity(recv_buffer_ptr, MED_CELL,recv_buffer[i]);
+ recv_buffer_ptr+=nbtypes[i]*(recv_buffer[i]%100);
+ }
+ delete[] nbtypes;
}
- distant_ids_recv=new int [distant_nb_elems];
- for (int i=0; i<distant_nb_elems; i++)
+ distant_ids_recv=new int [distant_nb_elems+distant_nbpolys];
+ for (int i=0; i<distant_nb_elems+distant_nbpolys; i++)
{
distant_ids_recv[i]=*recv_buffer_ptr++;
}
meshing->setMeshDimension(distant_mesh_dim);
+ if ( distant_nbpolys )
+ {
+ poly_conn = recv_buffer_ptr;
+ recv_buffer_ptr+=distant_nbpolyconn;
+ poly_index = recv_buffer_ptr;
+ recv_buffer_ptr+=distant_nbpolys+1;
+ if ( distant_nbpolyface ) {
+ face_index = recv_buffer_ptr;
+ meshing->setPolyhedraConnectivity( poly_index, face_index, poly_conn,
+ distant_nbpolys, MED_CELL );
+ } else {
+ meshing->setPolygonsConnectivity( poly_index, poly_conn,
+ distant_nbpolys, MED_CELL );
+ }
+ }
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
-
+ delete[] recv_coords; // the coordinates are present if the recv_buffer is not empty
+
+ {
+ MEDMEM::STRING file(getenv("TMP"));
+ file<<"/_" << this <<"_"<<iprocdistant_in_union<<"_local_mesh.med";
+ local_mesh->addDriver(MEDMEM::MED_DRIVER, file, "local_mesh");
+ cout << "WRITE " << file << endl;
+ local_mesh->write();
+ }
+ {
+ MEDMEM::STRING file(getenv("TMP"));
+ file<<"/_" << this <<"_"<<iprocdistant_in_union<<"_distant_mesh.med";
+ distant_mesh->addDriver(MEDMEM::MED_DRIVER, file, "distant_mesh");
+ cout << "WRITE " << file << endl;
+ distant_mesh->write();
+ }
}
MEDMEM::MESH* ElementLocator::_meshFromElems(set<int>& elems)
{
- //returns null pointer if there are no elems in the mesh
- if (elems.size()==0) return 0;
+ //returns null pointer if there are no elems in the mesh
+ if (elems.size()==0) return 0;
+
+ int dim = _local_mesh->getSpaceDimension();
+ int nbelem = _local_mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
//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;
+ const int *conn_mesh, *conn_index;
+ const double* coords = _local_mesh->getCoordinates(MED_FULL_INTERLACE);
+ if ( nbelem ) {
+ conn_mesh=_local_mesh->getConnectivity(MED_FULL_INTERLACE,
+ MED_NODAL,
+ MED_CELL,
+ MED_ALL_ELEMENTS);
+ conn_index=_local_mesh->getConnectivityIndex(
+ MED_NODAL,
+ MED_CELL);
+ }
+ map<int,int> old2newNodes;
int nbconn=0;
- map<MED_EN::medGeometryElement,int> nbelems_per_type;
- for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++)
+ map<medGeometryElement,int> nbelems_per_type;
+ set<int>::const_iterator eIter,eEnd = elems.end(), polyBeg = elems.upper_bound( nbelem );
+ for (eIter=elems.begin(); eIter!=polyBeg; eIter++) // loop on standard elements
{
- 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);
+ int n1 = conn_index[*eIter-1], n2 = conn_index[*eIter];
+ for (int inode = n1; inode < n2; inode++)
+ old2newNodes.insert( make_pair( conn_mesh[inode-1], old2newNodes.size()+1 ));
+ nbconn+=n2-n1;
+ medGeometryElement type = _local_mesh->getElementType(MED_CELL,*eIter);
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++)
+ int dimmax=0;
+ if ( !nbelems_per_type.empty() )
+ dimmax = nbelems_per_type.rbegin()->first / 100;
+
+ MEDMEM::MESHING* meshing = new MEDMEM::MESHING();
+ meshing->setNumberOfTypes(nbelems_per_type.size(),MED_CELL);
+
+ if ( polyBeg != eEnd ) // poly elements
{
- small2big[i]=*iter;
- big2small[*iter]=i;
- i++;
+ int nbPolygons = _local_mesh->getNumberOfPolygons(MED_CELL);
+ int nbPolyhedrons = _local_mesh->getNumberOfPolyhedron();
+ if ( nbPolygons > 0 ) // 2D
+ {
+ if ( dimmax < 2 ) dimmax = 2;
+ meshing->setMeshDimension(dimmax);
+
+ conn_mesh=_local_mesh->getPolygonsConnectivity(MED_NODAL,
+ MED_CELL);
+ conn_index=_local_mesh->getPolygonsConnectivityIndex(
+ MED_NODAL,
+ MED_CELL);
+ vector< int > polyConn, polyIndex;
+ polyConn.reserve( _local_mesh->getPolygonsConnectivityLength(MED_NODAL,MED_CELL));
+ polyIndex.reserve( nbPolygons + 1 );
+ polyIndex.push_back( 1 );
+ for ( eIter = polyBeg; eIter!=eEnd; eIter++)
+ {
+ int n1 = conn_index[*eIter-1], n2 = conn_index[*eIter];
+ for (int inode = n1; inode < n2; inode++)
+ {
+ map<int,int>::iterator old_new =
+ old2newNodes.insert( make_pair( conn_mesh[inode-1], old2newNodes.size()+1 )).first;
+ polyConn.push_back( old_new->second );
+ }
+ polyIndex.push_back( polyConn.size() + 1 );
+ }
+ meshing->setPolygonsConnectivity( & polyIndex[0],
+ & polyConn[0],
+ polyIndex.size()-1,
+ MED_CELL);
+ }
+ else if ( nbPolyhedrons > 0 ) // 3D
+ {
+ if ( dimmax < 3 ) dimmax = 3;
+ meshing->setMeshDimension(dimmax);
+
+ conn_mesh =_local_mesh->getPolyhedronConnectivity(MED_NODAL);
+ conn_index =_local_mesh->getPolyhedronIndex(MED_NODAL);
+ const int* face_index = _local_mesh->getPolyhedronFacesIndex();
+
+ vector< int > polyConn, polyIndex, faceIndex;
+ polyConn.reserve( _local_mesh->getPolyhedronConnectivityLength(MED_NODAL) );
+ polyIndex.reserve( nbPolyhedrons + 1 );
+ polyIndex.push_back( 1 );
+ faceIndex.reserve( _local_mesh->getNumberOfPolyhedronFaces() + 1 );
+ faceIndex.push_back( 1 );
+ for ( eIter = polyBeg; eIter!=eEnd; eIter++)
+ {
+ int f1 = conn_index[*eIter-1], f2 = conn_index[*eIter];
+ int n1 = face_index[f1-1], n2 = face_index[f2-1];
+ for (int inode = n1; inode < n2; inode++)
+ {
+ map<int,int>::iterator old_new =
+ old2newNodes.insert( make_pair( conn_mesh[inode-1], old2newNodes.size()+1 )).first;
+ polyConn.push_back( old_new->second );
+ }
+ for (int iface = f1; iface < f2; ++iface )
+ {
+ faceIndex.push_back( faceIndex.back() + ( face_index[iface] - face_index[iface-1] ));
+ }
+ polyIndex.push_back( faceIndex.size() );
+ }
+ meshing->setPolyhedraConnectivity( & polyIndex[0],
+ & faceIndex[0],
+ & polyConn[0],
+ polyIndex.size()-1,
+ MED_CELL);
+ }
}
+
int* conn= new int[nbconn];
- double* new_coords = new double[nodes.size()*_local_mesh->getSpaceDimension()];
+ double* new_coords = new double[old2newNodes.size()*dim];
int index=0;
- for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++)
+ for (set<int>::const_iterator iter=elems.begin(); iter!=polyBeg; iter++)
{
for (int inode = conn_index[*iter-1]-1; inode < conn_index[*iter]-1; inode++)
{
- conn[index]=big2small[conn_mesh[inode]]+1;
+ conn[index]=old2newNodes[conn_mesh[inode]];
index++;
- }
+ }
}
- index=0;
- for (set<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
+ const size_t size = size_t( sizeof(double)*dim );
+ for (map<int,int>::const_iterator iter=old2newNodes.begin(); iter!=old2newNodes.end(); iter++)
{
- int dim = _local_mesh->getSpaceDimension();
- for (int i=0; i<dim;i++)
- {
- new_coords[index]=coords[(*iter-1)*dim+i];
- index++;
- }
+ memcpy( new_coords+dim*(iter->second-1), coords+dim*(iter->first-1), size );
}
-
+
int* nbtypes=new int[nbelems_per_type.size()];
- MED_EN::medGeometryElement* new_types = new MED_EN::medGeometryElement[nbelems_per_type.size()];
+ medGeometryElement* new_types = new 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);
+ for (map<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++;
+ }
+ meshing->setMeshDimension(dimmax);
+ meshing->setCoordinates(_local_mesh->getSpaceDimension(), old2newNodes.size(), new_coords, string("CARTESIAN"), MED_FULL_INTERLACE);
+ meshing->setTypes(new_types,MED_CELL);
+ meshing->setNumberOfElements(nbtypes,MED_CELL);
- int dimmax=0;
- int* conn_ptr= conn;
+ 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]);
+ meshing->setConnectivity(conn_ptr, MED_CELL,new_types[i]);
conn_ptr+=nbtypes[i]*(new_types[i]%100);
- if (new_types[i]/100>dimmax) dimmax=new_types[i]/100;
}
- meshing->setMeshDimension(dimmax);
- delete [] small2big;
delete [] nbtypes;
delete [] conn;
delete [] new_coords;
