From: eap Date: Wed, 17 Dec 2008 17:45:50 +0000 (+0000) Subject: MEDMEM Industrialization 2008 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=353caff97672f1c0900070ddf86775de86289bd8;p=tools%2Fmedcoupling.git MEDMEM Industrialization 2008 exchange poly meshes --- diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index 8f80feab3..351194a36 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -1,3 +1,21 @@ +// 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 #include "CommInterface.hxx" #include "ElementLocator.hxx" @@ -13,6 +31,7 @@ #include #include using namespace std; +using namespace MED_EN; namespace ParaMEDMEM @@ -41,6 +60,42 @@ ElementLocator::~ElementLocator() delete [] _domain_bounding_boxes; } +void ElementLocator::_getElemsIntesectingBndBox(const int *conn, const int *conn_index, + int nbElems, int firstElemId, double* distant_bb, + set& 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::max(); + empty_bb[i*2+1]=-std::numeric_limits::max(); + } + const double* coords = _local_mesh->getCoordinates(MED_FULL_INTERLACE); + const size_t size = size_t( sizeof(double)*2*dim ); + for (int ielem=0; ielemelem_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 @@ -51,69 +106,80 @@ ElementLocator::~ElementLocator() void ElementLocator::exchangeMesh(int idistantrank, MEDMEM::MESH*& distant_mesh, int*& distant_ids) { int dim=_local_mesh->getSpaceDimension(); - int rank= _union_group->translateRank(&_distant_group,idistantrank); - if (find(_distant_proc_ids.begin(), _distant_proc_ids.end(),rank)==_distant_proc_ids.end()) - return; - - set elems; - double* distant_bb = _domain_bounding_boxes+rank*2*dim; - double* elem_bb=new double[2*dim]; - - //defining pointers to med - const int* conn=_local_mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE, - MED_EN::MED_NODAL, - MED_EN::MED_CELL, - MED_EN::MED_ALL_ELEMENTS); - const int* conn_index= _local_mesh->getConnectivityIndex( - MED_EN::MED_NODAL, - MED_EN::MED_CELL); - const double* coords = _local_mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE); - - for (int ielem=0; ielem<_local_mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);ielem++) - { - for (int i=0; i::max(); - elem_bb[i*2+1]=-std::numeric_limits::max(); - } - for (int inode=conn_index[ielem]; inodeelem_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::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 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 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::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; } @@ -122,7 +188,7 @@ 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); + const double* coords = _local_mesh->getCoordinates(MED_FULL_INTERLACE); int nbnodes = _local_mesh->getNumberOfNodes(); double * minmax=new double [2*dim]; @@ -192,234 +258,387 @@ bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int 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 (_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 (_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 (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 (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; isetTypes(types_vector, MED_EN::MED_CELL); - delete[] types_vector; - - recv_buffer_ptr+=distant_nb_types; - int* nbtypes = new int[distant_nb_types]; - for (int i=0; isetNumberOfElements( nbtypes, MED_EN::MED_CELL); - - for (int i=0; isetConnectivity(recv_buffer_ptr, MED_EN::MED_CELL,recv_buffer[i]); - recv_buffer_ptr+=nbtypes[i]*(recv_buffer[i]%100); + medGeometryElement* types_vector = new medGeometryElement[distant_nb_types]; + for (int i=0; isetTypes(types_vector, MED_CELL); + delete[] types_vector; + + recv_buffer_ptr+=distant_nb_types; + int* nbtypes = new int[distant_nb_types]; + for (int i=0; isetNumberOfElements( nbtypes, MED_CELL); + + for (int i=0; isetConnectivity(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; isetMeshDimension(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 <<"_"<addDriver(MEDMEM::MED_DRIVER, file, "local_mesh"); + cout << "WRITE " << file << endl; + local_mesh->write(); + } + { + MEDMEM::STRING file(getenv("TMP")); + file<<"/_" << this <<"_"<addDriver(MEDMEM::MED_DRIVER, file, "distant_mesh"); + cout << "WRITE " << file << endl; + distant_mesh->write(); + } } MEDMEM::MESH* ElementLocator::_meshFromElems(set& 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 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 old2newNodes; int nbconn=0; - map nbelems_per_type; - for (set::const_iterator iter=elems.begin(); iter!=elems.end(); iter++) + map nbelems_per_type; + set::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 big2small; - int i=0; - for (set::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::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::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::const_iterator iter=elems.begin(); iter!=elems.end(); iter++) + for (set::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::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++) + const size_t size = size_t( sizeof(double)*dim ); + for (map::const_iterator iter=old2newNodes.begin(); iter!=old2newNodes.end(); iter++) { - int dim = _local_mesh->getSpaceDimension(); - for (int i=0; isecond-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::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::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; isetConnectivity(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; diff --git a/src/ParaMEDMEM/ElementLocator.hxx b/src/ParaMEDMEM/ElementLocator.hxx index a76024003..6149117ad 100644 --- a/src/ParaMEDMEM/ElementLocator.hxx +++ b/src/ParaMEDMEM/ElementLocator.hxx @@ -1,3 +1,21 @@ +// 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 +// #ifndef ELEMENTLOCATOR_HXX_ #define ELEMENTLOCATOR_HXX_ @@ -17,12 +35,12 @@ class ParaSUPPORT; class InterpolationMatrix; - class ElementLocator: public INTERP_KERNEL::InterpolationOptions +class ElementLocator: public INTERP_KERNEL::InterpolationOptions { public: - ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group); + ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group); ElementLocator(const ParaSUPPORT& support, const ProcessorGroup& distant_group); - virtual ~ElementLocator(); + virtual ~ElementLocator(); void exchangeMesh(int idistantrank, MEDMEM::MESH*& distant_mesh, int*& distant_ids); private: @@ -35,6 +53,7 @@ private: std::vector _distant_proc_ids; void _computeBoundingBoxes(); + void _getElemsIntesectingBndBox(const int *conn, const int *conn_index, int nbElems, int firstElemId, double* bb, std::set& elems); 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);