-// Copyright (C) 2007-2014 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2023 CEA, EDF
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
#include "MEDPARTITIONER_MeshCollection.hxx"
+
+#include "MEDPARTITIONER_ConnectZone.hxx"
+#include "MEDPARTITIONER_Graph.hxx"
#include "MEDPARTITIONER_MeshCollectionDriver.hxx"
-#include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx"
#include "MEDPARTITIONER_MeshCollectionMedAsciiDriver.hxx"
+#include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx"
#include "MEDPARTITIONER_ParaDomainSelector.hxx"
-#include "MEDPARTITIONER_Topology.hxx"
#include "MEDPARTITIONER_ParallelTopology.hxx"
+#include "MEDPARTITIONER_Topology.hxx"
+#include "MEDPARTITIONER_UserGraph.hxx"
+#include "MEDPARTITIONER_Utils.hxx"
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
#include "MEDPARTITIONER_JointFinder.hxx"
#endif
-#include "MEDPARTITIONER_Graph.hxx"
-#include "MEDPARTITIONER_UserGraph.hxx"
-#include "MEDPARTITIONER_Utils.hxx"
-
-#include "MEDLoaderBase.hxx"
-#include "MEDLoader.hxx"
+#include "MCAuto.hxx"
+#include "MEDCouplingFieldDouble.hxx"
#include "MEDCouplingMemArray.hxx"
-#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.hxx"
-#include "MEDCouplingFieldDouble.hxx"
-#include "PointLocator3DIntersectorP0P0.hxx"
-
-#include "MEDCouplingAutoRefCountObjectPtr.hxx"
+#include "MEDCouplingSkyLineArray.hxx"
+#include "MEDCouplingUMesh.hxx"
+#include "MEDLoader.hxx"
+#include "MEDLoaderBase.hxx"
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
#include <mpi.h>
#endif
_domain_selector( 0 ),
_i_non_empty_mesh(-1),
_driver_type(MEDPARTITIONER::MedXml),
- _subdomain_boundary_creates( MyGlobals::_Creates_Boundary_Faces ),
+ _subdomain_boundary_creates( MyGlobals::_Create_Boundary_Faces ),
_family_splitting(false),
_create_empty_groups(false),
_joint_finder(0)
_i_non_empty_mesh(-1),
_name(initialCollection._name),
_driver_type(MEDPARTITIONER::MedXml),
- _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
+ _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
_family_splitting(family_splitting),
_create_empty_groups(create_empty_groups),
_joint_finder(0)
{
- std::vector<std::vector<std::vector<int> > > new2oldIds(initialCollection.getTopology()->nbDomain());
- castCellMeshes(initialCollection, new2oldIds);
+ std::vector<std::vector<std::vector<mcIdType> > > new2oldIds(initialCollection.getTopology()->nbDomain());
+ std::vector<MEDCoupling::DataArrayIdType*> o2nRenumber;
+
+ castCellMeshes(initialCollection, new2oldIds, o2nRenumber );
//defining the name for the collection and the underlying meshes
setName(initialCollection.getName());
//treating faces
/////////////////
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
#endif
//nodeMapping contains the mapping between old nodes and new nodes
// (iolddomain,ioldnode)->(inewdomain,inewnode)
createNodeMapping(initialCollection, nodeMapping);
- std::vector<std::vector<std::vector<int> > > new2oldFaceIds;
+ std::vector<std::vector<std::vector<mcIdType> > > new2oldFaceIds;
castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds);
////////////////////
//treating families
////////////////////
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
#endif
this->getMesh(),
initialCollection.getCellFamilyIds(),
"cellFamily");
- castIntField(initialCollection.getFaceMesh(),
+ castIntField(initialCollection.getFaceMesh(),
this->getFaceMesh(),
initialCollection.getFaceFamilyIds(),
"faceFamily");
//treating groups
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
#endif
std::cout << "treating groups" << std::endl;
_family_info=initialCollection.getFamilyInfo();
_group_info=initialCollection.getGroupInfo();
-
-#ifdef HAVE_MPI2
+
+#ifdef HAVE_MPI
if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
#endif
castAllFields(initialCollection,"cellFieldDouble");
if (_i_non_empty_mesh<0)
{
- for (int i=0; i<_mesh.size(); i++)
+ for (int i=0; i<getNbOfGlobalMeshes(); i++)
{
if (_mesh[i])
{
}
}
+ // find faces common with neighbor domains and put them in groups
+ buildBoundaryFaces();
+
+ //building the connect zones necessary for writing joints
+ buildConnectZones( nodeMapping, o2nRenumber, initialCollection.getTopology()->nbDomain() );
+
+ // delete o2nRenumber
+ for ( size_t i = 0; i < o2nRenumber.size(); ++i )
+ if ( o2nRenumber[i] )
+ o2nRenumber[i]->decrRef();
}
/*!
- Creates the meshes using the topology underlying he mesh collection and the mesh data
+ Creates the meshes using the topology underlying he mesh collection and the mesh data
coming from the ancient collection
\param initialCollection collection from which the data is extracted to create the new meshes
+ \param [out] o2nRenumber returns for each new domain a permutation array returned by sortCellsInMEDFileFrmt()
*/
void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection,
- std::vector<std::vector<std::vector<int> > >& new2oldIds)
+ std::vector<std::vector<std::vector<mcIdType> > >& new2oldIds,
+ std::vector<MEDCoupling::DataArrayIdType*> & o2nRenumber)
{
if (MyGlobals::_Verbose>10)
std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl;
int nbOldDomain=initialCollection.getTopology()->nbDomain();
_mesh.resize(nbNewDomain);
+ o2nRenumber.resize(nbNewDomain,0);
int rank=MyGlobals::_Rank;
//splitting the initial domains into smaller bits
- std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
+ std::vector<std::vector<MEDCoupling::MEDCouplingUMesh*> > splitMeshes;
splitMeshes.resize(nbNewDomain);
for (int inew=0; inew<nbNewDomain; inew++)
{
- splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
+ splitMeshes[inew].resize(nbOldDomain, (MEDCoupling::MEDCouplingUMesh*)0);
}
for (int iold=0; iold<nbOldDomain; iold++)
{
if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
{
- int size=(initialCollection._mesh)[iold]->getNumberOfCells();
- std::vector<int> globalids(size);
+ mcIdType size=(initialCollection._mesh)[iold]->getNumberOfCells();
+ std::vector<mcIdType> globalids(size);
initialCollection.getTopology()->getCellList(iold, &globalids[0]);
- std::vector<int> ilocalnew(size); //local
+ std::vector<mcIdType> ilocalnew(size); //local
std::vector<int> ipnew(size); //idomain old
_topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
}
for (int inew=0; inew<nbNewDomain; inew++)
{
- splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)
+ splitMeshes[inew][iold]=(MEDCoupling::MEDCouplingUMesh*)
(initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
&new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
true);
}
}
}
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if (isParallelMode())
{
//if (MyGlobals::_Verbose>300) std::cout<<"proc "<<rank<<" : castCellMeshes send/receive"<<std::endl;
std::cout << "proc " << rank << " : castCellMeshes fusing" << std::endl;
for (int inew=0; inew<nbNewDomain ;inew++)
{
- std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
-
+ std::vector<const MEDCoupling::MEDCouplingUMesh*> meshes;
+
for (int i=0; i<(int)splitMeshes[inew].size(); i++)
- if (splitMeshes[inew][i]!=0)
+ if (splitMeshes[inew][i]!=0)
if (splitMeshes[inew][i]->getNumberOfCells()>0)
meshes.push_back(splitMeshes[inew][i]);
- if (!isParallelMode()||_domain_selector->isMyDomain(inew))
+ if (!isParallelMode()||_domain_selector->isMyDomain(inew))
+ {
+ if (meshes.size()==0)
{
- if (meshes.size()==0)
- {
_mesh[inew]=CreateEmptyMEDCouplingUMesh();
std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl;
}
else
{
- _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes);
+ _mesh[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshes(meshes);
+ o2nRenumber[inew]=_mesh[inew]->sortCellsInMEDFileFrmt();
bool areNodesMerged;
- int nbNodesMerged;
+ mcIdType nbNodesMerged;
if (meshes.size()>1)
{
- ParaMEDMEM::DataArrayInt* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
+ MEDCoupling::DataArrayIdType* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
array->decrRef(); // array is not used in this case
}
_mesh[inew]->zipCoords();
-
}
}
for (int i=0;i<(int)splitMeshes[inew].size();i++)
}
/*!
- \param initialCollection source mesh collection
+ \param initialCollection source mesh collection
\param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection
*/
void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping)
if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
{
// std::map<pair<double,pair<double, double> >, int > nodeClassifier;
- ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
+ MEDCoupling::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
double* coordsPtr=coords->getPointer();
- dim = coords->getNumberOfComponents();
- int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
+ dim = (int)coords->getNumberOfComponents();
+ mcIdType nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
bbox=new double[nvertices*2*dim];
for (int i=0; i<nvertices*dim;i++)
for (int inew=0; inew<_topology->nbDomain(); inew++)
{
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
//sending meshes for parallel computation
if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold))
_domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold));
else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold))
{
- ParaMEDMEM::MEDCouplingUMesh* mesh;
+ MEDCoupling::MEDCouplingUMesh* mesh;
_domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
- ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords();
+ MEDCoupling::DataArrayDouble* coords = mesh->getCoords();
for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
{
double* coordsPtr=coords->getPointer()+inode*dim;
- vector<int> elems;
+ vector<mcIdType> elems;
tree->getElementsAroundPoint(coordsPtr,elems);
if (elems.size()==0) continue;
nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
#endif
{
- ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords();
+ MEDCoupling::DataArrayDouble* coords = getMesh(inew)->getCoords();
for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
{
double* coordsPtr=coords->getPointer()+inode*dim;
- vector<int> elems;
+ vector<mcIdType> elems;
tree->getElementsAroundPoint(coordsPtr,elems);
if (elems.size()==0) continue;
nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
}
-void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
+void getNodeIds(MEDCoupling::MEDCouplingUMesh& meshOne, MEDCoupling::MEDCouplingUMesh& meshTwo, std::vector<mcIdType>& nodeIds)
{
using std::vector;
using MEDPARTITIONER::BBTreeOfDim;
- if (!&meshOne || !&meshTwo) return; //empty or not existing
+ //if (!&meshOne || !&meshTwo) return; //empty or not existing
double* bbox;
BBTreeOfDim* tree = 0;
- int nv1=meshOne.getNumberOfNodes();
- ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords();
- int dim = coords->getNumberOfComponents();
+ mcIdType nv1=meshOne.getNumberOfNodes();
+ MEDCoupling::DataArrayDouble* coords=meshOne.getCoords();
+ int dim = (int)coords->getNumberOfComponents();
bbox=new double[nv1*2*dim];
double* coordsPtr=coords->getPointer();
}
tree=new BBTreeOfDim( dim, bbox,0,0,nv1,1e-9);
- int nv2=meshTwo.getNumberOfNodes();
+ mcIdType nv2=meshTwo.getNumberOfNodes();
nodeIds.resize(nv2,-1);
coords=meshTwo.getCoords();
for (int inode=0; inode<nv2; inode++)
{
double* coordsPtr2=coords->getPointer()+inode*dim;
- vector<int> elems;
+ vector<mcIdType> elems;
tree->getElementsAroundPoint(coordsPtr2,elems);
if (elems.size()==0) continue;
nodeIds[inode]=elems[0];
faces at the interface are duplicated
*/
void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialCollection,
- const std::multimap<std::pair<int,int>, std::pair<int,int> >& nodeMapping,
- std::vector<std::vector<std::vector<int> > >& new2oldIds)
+ const std::multimap<std::pair<int,mcIdType>, std::pair<int,mcIdType> >& nodeMapping,
+ std::vector<std::vector<std::vector<mcIdType> > >& new2oldIds)
{
//splitMeshes structure will contain the partition of
//the old faces on the new ones
int nbNewDomain=_topology->nbDomain();
int nbOldDomain=initialCollection.getTopology()->nbDomain();
- vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
- vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
+ vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
+ vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
- vector< vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
+ vector< vector<MEDCoupling::MEDCouplingUMesh*> > splitMeshes;
splitMeshes.resize(nbNewDomain);
for (int inew=0; inew<nbNewDomain; inew++)
{
- splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
+ splitMeshes[inew].resize(nbOldDomain, (MEDCoupling::MEDCouplingUMesh*)0);
}
new2oldIds.resize(nbOldDomain);
for (int iold=0; iold<nbOldDomain; iold++) new2oldIds[iold].resize(nbNewDomain);
//initial face mesh known : in my domain
if (meshesCastFrom[iold] != 0)
{
- for (int ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
+ for (mcIdType ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
{
- vector<int> nodes;
+ vector<mcIdType> nodes;
meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
- map <int,int> faces;
+ map <int,mcIdType> faces;
//analysis of element ielem
//counters are set for the element
//are incremented for each target node
//the face is considered as going to target domains if the counter of the domain
//is equal to the number of nodes
- for (int inode=0; inode<(int)nodes.size(); inode++)
+ for (std::size_t inode=0; inode<nodes.size(); inode++)
{
- typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
- int mynode=nodes[inode];
+ typedef multimap<pair<int,mcIdType>,pair<int,mcIdType> >::const_iterator MI;
+ mcIdType mynode=nodes[inode];
pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
for (MI iter=myRange.first; iter!=myRange.second; iter++)
}
}
- for (map<int,int>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
+ for (map<int,mcIdType>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
{
- if (iter->second==(int)nodes.size())
+ if (iter->second==ToIdType(nodes.size()))
//cvw eligible but may be have to be face of a cell of this->getMesh()[inew]?
//it is not sure here...
//done before writeMedfile on option?... see filterFaceOnCell()
if (meshesCastFrom[iold]->getNumberOfCells() > 0)
{
splitMeshes[inew][iold]=
- (ParaMEDMEM::MEDCouplingUMesh*)
+ (MEDCoupling::MEDCouplingUMesh*)
( meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
&new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
true)
}
}
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
//send/receive stuff
if (isParallelMode())
{
- ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
+ MEDCoupling::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
for (int iold=0; iold<nbOldDomain; iold++)
for (int inew=0; inew<nbNewDomain; inew++)
{
meshesCastTo.resize(nbNewDomain);
for (int inew=0; inew<nbNewDomain; inew++)
{
- vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
+ vector<MEDCoupling::MEDCouplingUMesh*> myMeshes;
for (int iold=0; iold<nbOldDomain; iold++)
{
- ParaMEDMEM::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
+ MEDCoupling::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
if (umesh!=0)
if (umesh->getNumberOfCells()>0)
myMeshes.push_back(umesh);
+
}
- ParaMEDMEM::MEDCouplingUMesh *bndMesh = 0;
+ MEDCoupling::MEDCouplingUMesh *bndMesh = 0;
if ( _subdomain_boundary_creates &&
_mesh[inew] &&
_mesh[inew]->getNumberOfCells()>0 )
{
bndMesh =
- ((ParaMEDMEM::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
+ ((MEDCoupling::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
if (bndMesh->getNumberOfCells()>0)
- myMeshes.push_back( bndMesh );
+ myMeshes.push_back( bndMesh );
}
if (myMeshes.size()>0)
{
- meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
+ // Need to merge faces that are both in the initial set of faces and in the boundary mesh
+ // which is the last element in myMeshes:
+ if(bndMesh)
+ {
+ for (int iold2 = 0; iold2 < (int)myMeshes.size()-1; iold2++)
+ myMeshes[iold2]->tryToShareSameCoordsPermute(*bndMesh, 1.e-10);
+ vector<const MEDCoupling::MEDCouplingUMesh*> myMeshes_c;
+ for (auto & mesh: myMeshes) myMeshes_c.push_back(mesh);
+ meshesCastTo[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshesOnSameCoords(myMeshes_c);
+ MCAuto<DataArrayIdType> tmp = meshesCastTo[inew]->zipConnectivityTraducer(2);
+ }
+ else
+ {
+ vector<const MEDCoupling::MEDCouplingUMesh*> myMeshes_c;
+ for (auto & mesh: myMeshes) myMeshes_c.push_back(mesh);
+ meshesCastTo[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshes(myMeshes_c);
+ }
meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef();
}
else
{
- ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
+ MEDCoupling::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
meshesCastTo[inew]=empty;
}
for (int iold=0; iold<nbOldDomain; iold++)
-void MEDPARTITIONER::MeshCollection::castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
- std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
- std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
+void MEDPARTITIONER::MeshCollection::castIntField(std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastFrom,
+ std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastTo,
+ std::vector<MEDCoupling::DataArrayIdType*>& arrayFrom,
std::string nameArrayTo)
{
using std::vector;
- int ioldMax=meshesCastFrom.size();
- int inewMax=meshesCastTo.size();
+ std::size_t ioldMax=meshesCastFrom.size();
+ std::size_t inewMax=meshesCastTo.size();
//preparing bounding box trees for accelerating source-target node identifications
if (MyGlobals::_Verbose>99)
std::cout<<"making accelerating structures"<<std::endl;
std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
- std::vector<ParaMEDMEM::DataArrayDouble*>bbox(ioldMax);
- for (int iold =0; iold< ioldMax; iold++)
+ std::vector<MEDCoupling::DataArrayDouble*>bbox(ioldMax);
+ for (unsigned int iold =0; iold< ioldMax; iold++)
if (isParallelMode() && _domain_selector->isMyDomain(iold))
{
- ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
+ MEDCoupling::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->computeCellCenterOfMass();
bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
sourceCoords->decrRef();
}
// send-recv operations
-#ifdef HAVE_MPI2
- for (int inew=0; inew<inewMax; inew++)
+#ifdef HAVE_MPI
+ for (unsigned int inew=0; inew<inewMax; inew++)
{
- for (int iold=0; iold<ioldMax; iold++)
+ for (unsigned int iold=0; iold<ioldMax; iold++)
{
//sending arrays for distant domains
if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
//send mesh
_domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
//send vector
- int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
- vector<int>sendIds;
+ mcIdType size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
+ vector<mcIdType> sendIds;
if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
if (size>0) //no empty
{
if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
{
//receive mesh
- vector<int> recvIds;
- ParaMEDMEM::MEDCouplingUMesh* recvMesh;
+ vector<mcIdType> recvIds;
+ MEDCoupling::MEDCouplingUMesh* recvMesh;
_domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
//receive vector
if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
}
#endif
- //local contributions and aggregation
- for (int inew=0; inew<inewMax; inew++)
+ //localc ontributions and aggregation
+ for (unsigned int inew=0; inew<inewMax; inew++)
{
- for (int iold=0; iold<ioldMax; iold++)
+ for (unsigned int iold=0; iold<ioldMax; iold++)
if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
{
remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
}
}
- for (int iold=0; iold<ioldMax;iold++)
+ for (unsigned int iold=0; iold<ioldMax;iold++)
if (isParallelMode() && _domain_selector->isMyDomain(iold))
{
bbox[iold]->decrRef();
}
void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
- const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
- const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
- const int* fromArray,
+ const MEDCoupling::MEDCouplingUMesh& sourceMesh,
+ const MEDCoupling::MEDCouplingUMesh& targetMesh,
+ const mcIdType* fromArray,
std::string nameArrayTo,
const BBTreeOfDim* myTree)
{
if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
- ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
+ MEDCoupling::DataArrayDouble* targetCoords=targetMesh.computeCellCenterOfMass();
const double* tc=targetCoords->getConstPointer();
- int targetSize=targetMesh.getNumberOfCells();
- int sourceSize=sourceMesh.getNumberOfCells();
+ mcIdType targetSize=targetMesh.getNumberOfCells();
+ mcIdType sourceSize=sourceMesh.getNumberOfCells();
if (MyGlobals::_Verbose>200)
std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
- std::vector<int> ccI;
+ std::vector<mcIdType> ccI;
std::string str,cle;
str=nameArrayTo+"_toArray";
cle=Cle1ToStr(str,inew);
- int* toArray;
+ mcIdType* toArray;
const BBTreeOfDim* tree;
bool cleantree=false;
- ParaMEDMEM::DataArrayDouble* sourceBBox=0;
- int dim = targetCoords->getNumberOfComponents();
+ MEDCoupling::DataArrayDouble* sourceBBox=0;
+ int dim = (int)targetCoords->getNumberOfComponents();
if (myTree==0)
{
- sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8);
+ sourceBBox=sourceMesh.computeCellCenterOfMass()->computeBBoxPerTuple(1e-8);
tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
cleantree=true;
}
{
if (MyGlobals::_Is0verbose>100)
std::cout << "create " << cle << " size " << targetSize << std::endl;
- ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
+ MEDCoupling::DataArrayIdType* p=MEDCoupling::DataArrayIdType::New();
p->alloc(targetSize,1);
p->fillWithZero();
toArray=p->getPointer();
toArray=_map_dataarray_int.find(cle)->second->getPointer();
}
- std::map< int, int > isource2nb; // count coincident elements
- std::map<int,int>::iterator i2nb;
+ std::map< mcIdType, int > isource2nb; // count coincident elements
+ std::map< mcIdType, int>::iterator i2nb;
- for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
+ for (mcIdType itargetnode=0; itargetnode<targetSize; itargetnode++)
{
- std::vector<int> intersectingElems;
+ std::vector<mcIdType> intersectingElems;
tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
if (intersectingElems.size()!=0)
{
- int isourcenode=intersectingElems[0];
+ mcIdType isourcenode=intersectingElems[0];
if ( intersectingElems.size() > 1 )
{
i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
std::string nameTo="typeData=6"; //resume the type of field casted
// send-recv operations
- int ioldMax=initialCollection.getMesh().size();
- int inewMax=this->getMesh().size();
- int iFieldMax=initialCollection.getFieldDescriptions().size();
+ std::size_t ioldMax=initialCollection.getMesh().size();
+ std::size_t inewMax=this->getMesh().size();
+ std::size_t iFieldMax=initialCollection.getFieldDescriptions().size();
if (MyGlobals::_Verbose>10)
std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
//see collection.prepareFieldDescriptions()
- for (int ifield=0; ifield<iFieldMax; ifield++)
+ for (std::size_t ifield=0; ifield<iFieldMax; ifield++)
{
std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
if (descriptionField.find(nameTo)==std::string::npos)
continue; //only nameTo accepted in Fields name description
-#ifdef HAVE_MPI2
- for (int inew=0; inew<inewMax; inew++)
+#ifdef HAVE_MPI
+ for (unsigned int inew=0; inew<inewMax; inew++)
{
- for (int iold=0; iold<ioldMax; iold++)
+ for (unsigned int iold=0; iold<ioldMax; iold++)
{
//sending arrays for distant domains
if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
{
int target=_domain_selector->getProcessorID(inew);
- ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
+ MEDCoupling::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
if (MyGlobals::_Verbose>10)
std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
SendDataArrayDouble(field, target);
//receive vector
if (MyGlobals::_Verbose>10)
std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
- ParaMEDMEM::DataArrayDouble* field=RecvDataArrayDouble(source);
+ MEDCoupling::DataArrayDouble* field=RecvDataArrayDouble(source);
remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
}
}
}
#endif
//local contributions and aggregation
- for (int inew=0; inew<inewMax; inew++)
+ for (unsigned int inew=0; inew<inewMax; inew++)
{
- for (int iold=0; iold<ioldMax; iold++)
+ for (unsigned int iold=0; iold<ioldMax; iold++)
if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
{
- ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
+ MEDCoupling::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
}
}
}
void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
- ParaMEDMEM::DataArrayDouble* fromArray,
+ MEDCoupling::DataArrayDouble* fromArray,
std::string nameArrayTo,
std::string descriptionField)
//here we use 'cellFamily_ccI inew iold' set in remapIntField
throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
- std::map<std::string,ParaMEDMEM::DataArrayInt*>::iterator it1;
+ std::map<std::string,MEDCoupling::DataArrayIdType*>::iterator it1;
it1=_map_dataarray_int.find(key);
if (it1==_map_dataarray_int.end())
{
return;
}
//create ccI in remapIntField
- ParaMEDMEM::DataArrayInt *ccI=it1->second;
+ MEDCoupling::DataArrayIdType *ccI=it1->second;
if (MyGlobals::_Verbose>300)
std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
- int nbcell=this->getMesh()[inew]->getNumberOfCells();
- int nbcomp=fromArray->getNumberOfComponents();
+ mcIdType nbcell=this->getMesh()[inew]->getNumberOfCells();
+ std::size_t nbcomp=fromArray->getNumberOfComponents();
int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
std::string tag="inewFieldDouble="+IntToStr(inew);
key=descriptionField+SerializeFromString(tag);
- int fromArrayNbOfElem=fromArray->getNbOfElems();
- int fromArrayNbOfComp=fromArray->getNumberOfComponents();
- int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
+ mcIdType fromArrayNbOfElem=fromArray->getNbOfElems();
+ mcIdType fromArrayNbOfComp=ToIdType(fromArray->getNumberOfComponents());
+ mcIdType fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
if (MyGlobals::_Verbose>1000)
{
" nbComponents " << fromArray->getNumberOfComponents() << std::endl;
}
- ParaMEDMEM::DataArrayDouble* field=0;
- std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it2;
+ MEDCoupling::DataArrayDouble* field=0;
+ std::map<std::string,MEDCoupling::DataArrayDouble*>::iterator it2;
it2=_map_dataarray_double.find(key);
if (it2==_map_dataarray_double.end())
{
if (MyGlobals::_Verbose>300)
std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
- field=ParaMEDMEM::DataArrayDouble::New();
+ field=MEDCoupling::DataArrayDouble::New();
_map_dataarray_double[key]=field;
field->alloc(nbcell*nbPtGauss,nbcomp);
field->fillWithZero();
else
{
//replaced by setPartOfValuesAdv if nbPtGauss==1
- int iMax=ccI->getNbOfElems();
- int* pccI=ccI->getPointer();
+ mcIdType iMax=ccI->getNbOfElems();
+ mcIdType* pccI=ccI->getPointer();
double* pField=field->getPointer();
double* pFrom=fromArray->getPointer();
- int itarget, isource, delta=nbPtGauss*nbcomp;
- for (int i=0; i<iMax; i=i+2) //cell
+ mcIdType itarget, isource, delta=ToIdType(nbPtGauss*nbcomp);
+ for (mcIdType i=0; i<iMax; i=i+2) //cell
{
itarget=pccI[i];
isource=pccI[i+1];
if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
throw INTERP_KERNEL::Exception("Error field override");
- int ita=itarget*delta;
- int iso=isource*delta;
- for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
+ mcIdType ita=itarget*delta;
+ mcIdType iso=isource*delta;
+ for (mcIdType k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
+ }
+ }
+}
+
+namespace
+{
+ using namespace MEDCoupling;
+ //================================================================================
+ /*!
+ * \brief Sort correspondence ids of one domain and permute ids of the other accordingly
+ * \param [in,out] ids1 - ids of one domain
+ * \param [in,out] ids2 - ids of another domain
+ * \param [in] delta - a delta to change all ids
+ * \param [in] removeEqual - to remove equal ids
+ * \return DataArrayInt* - array of ids joined back
+ */
+ //================================================================================
+
+ DataArrayIdType* sortCorrespondences( DataArrayIdType* ids1,
+ DataArrayIdType* ids2,
+ int delta,
+ bool removeEqual = false)
+ {
+ // sort
+ MCAuto< DataArrayIdType > renumN2O = ids1->buildPermArrPerLevel();
+ ids1->renumberInPlaceR( renumN2O->begin() );
+ ids2->renumberInPlaceR( renumN2O->begin() );
+
+ if ( removeEqual )
+ {
+ ids1 = ids1->buildUnique();
+ ids2 = ids2->buildUnique();
+ }
+ if ( delta != 0 )
+ {
+ mcIdType * id = ids1->getPointer();
+ for ( ; id < ids1->end(); ++id )
+ ++(*id);
+ id = ids2->getPointer();
+ for ( ; id < ids2->end(); ++id )
+ ++(*id);
+ }
+
+ // join
+ DataArrayIdType* ids12 = DataArrayIdType::Meld( ids1, ids2 ); // two components
+ ids12->rearrange( 1 ); // make one component
+ return ids12;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Renumber ids according to mesh->sortCellsInMEDFileFrmt()
+ * \param [in,out] ids - cell ids to renumber
+ * \param [in] o2nRenumber - renumbering array in "Old to New" mode
+ */
+ //================================================================================
+
+ void renumber( DataArrayIdType* ids, const DataArrayIdType* o2nRenumber )
+ {
+ if ( !ids || !o2nRenumber )
+ return;
+ mcIdType * id = ids->getPointer();
+ const mcIdType * o2n = o2nRenumber->getConstPointer();
+ for ( ; id < ids->end(); ++id )
+ {
+ *id = o2n[ *id ];
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Fill up ConnectZone's stored in _topology with nodal correspondences
+ * \param [in] nodeMapping - mapping between old nodes and new nodes
+ * (iolddomain,ioldnode)->(inewdomain,inewnode)
+ * \param [in] o2nRenumber - renumbering array returned by mesh->sortCellsInMEDFileFrmt()
+ * per a new domain
+ * \param [in] nbInitialDomains - nb of old domains
+ */
+//================================================================================
+
+void MEDPARTITIONER::MeshCollection::buildConnectZones( const NodeMapping& nodeMapping,
+ const std::vector<MEDCoupling::DataArrayIdType*> & o2nRenumber,
+ int nbInitialDomains)
+{
+ if ( !MyGlobals::_Create_Joints || _topology->nbDomain() < 2 )
+ return;
+
+ if ( MyGlobals::_World_Size > 1 )
+ {
+ _topology->getCZ().clear();
+ return; // not implemented for parallel mode
+ }
+
+ // at construction, _topology creates cell correspondences basing on Graph information,
+ // and here we
+ // 1) add node correspondences,
+ // 2) split cell correspondences by cell geometry types
+ // 3) sort ids to be in ascending order
+
+ const int nb_domains = _topology->nbDomain();
+
+ // ==================================
+ // 1) add node correspondences
+ // ==================================
+
+ std::vector< std::vector< std::vector< mcIdType > > > nodeCorresp( nb_domains );
+ for ( int idomain = 0; idomain < nb_domains; ++idomain )
+ {
+ nodeCorresp[ idomain ].resize( nb_domains );
+ }
+
+ NodeMapping::const_iterator nmIt1, nmIt2 = nodeMapping.begin();
+ for ( nmIt1 = nmIt2; nmIt1 != nodeMapping.end(); nmIt1 = nmIt2 )
+ {
+ // look for an "old" node mapped into several "new" nodes in different domains
+ int nbSameOld = 0;
+ while ( ++nmIt2 != nodeMapping.end() && nmIt2->first == nmIt1->first )
+ nbSameOld += ( nmIt2->second != nmIt1->second );
+
+ if ( nbSameOld > 0 )
+ {
+ NodeMapping::const_iterator nmEnd = nmIt2;
+ for ( ; true; ++nmIt1 )
+ {
+ nmIt2 = nmIt1;
+ if ( ++nmIt2 == nmEnd )
+ break;
+ int dom1 = nmIt1->second.first;
+ mcIdType node1 = nmIt1->second.second;
+ for ( ; nmIt2 != nmEnd; ++nmIt2 )
+ {
+ int dom2 = nmIt2->second.first;
+ mcIdType node2 = nmIt2->second.second;
+ if ( dom1 != dom2 )
+ {
+ nodeCorresp[ dom1 ][ dom2 ].push_back( node1 );
+ nodeCorresp[ dom1 ][ dom2 ].push_back( node2 );
+ nodeCorresp[ dom2 ][ dom1 ].push_back( node2 );
+ nodeCorresp[ dom2 ][ dom1 ].push_back( node1 );
+ }
+ }
+ }
+ }
+ }
+
+ // add nodeCorresp to czVec
+
+ std::vector<MEDPARTITIONER::ConnectZone*>& czVec = _topology->getCZ();
+
+ for ( int idomain = 0; idomain < nb_domains; ++idomain )
+ {
+ for ( int idomainNear = 0; idomainNear < nb_domains; ++idomainNear )
+ {
+ std::vector< mcIdType > & corresp = nodeCorresp[ idomain ][ idomainNear ];
+ if ( corresp.empty() )
+ continue;
+
+ MEDPARTITIONER::ConnectZone* cz = 0;
+ for ( size_t i = 0; i < czVec.size() && !cz; ++i )
+ if ( czVec[i] &&
+ czVec[i]->getLocalDomainNumber () == idomain &&
+ czVec[i]->getDistantDomainNumber() == idomainNear )
+ cz = czVec[i];
+
+ if ( !cz )
+ {
+ cz = new MEDPARTITIONER::ConnectZone();
+ cz->setName( "Nodal Connect Zone defined by MEDPARTITIONER" );
+ cz->setLocalDomainNumber ( idomain );
+ cz->setDistantDomainNumber( idomainNear );
+ czVec.push_back(cz);
+ }
+
+ cz->setNodeCorresp( &corresp[0], ToIdType( corresp.size()/2 ));
+ }
+ }
+
+ // ==========================================================
+ // 2) split cell correspondences by cell geometry types
+ // ==========================================================
+
+ for ( size_t i = 0; i < czVec.size(); ++i )
+ {
+ MEDPARTITIONER::ConnectZone* cz = czVec[i];
+ if ( !cz ||
+ cz->getEntityCorrespNumber( 0,0 ) == 0 ||
+ cz->getLocalDomainNumber () > (int)_mesh.size() ||
+ cz->getDistantDomainNumber() > (int)_mesh.size() )
+ continue;
+ MEDCoupling::MEDCouplingUMesh* mesh1 = _mesh[ cz->getLocalDomainNumber () ];
+ MEDCoupling::MEDCouplingUMesh* mesh2 = _mesh[ cz->getDistantDomainNumber() ];
+
+ // separate ids of two domains
+ const MEDCoupling::MEDCouplingSkyLineArray *corrArray = cz->getEntityCorresp( 0, 0 );
+ const DataArrayIdType* ids12 = corrArray->getValuesArray();
+ MCAuto<DataArrayIdType> ids1, ids2, ids12Sorted;
+ ids1 = ids12->selectByTupleIdSafeSlice( 0, corrArray->getLength(), 2 );
+ ids2 = ids12->selectByTupleIdSafeSlice( 1, corrArray->getLength(), 2 );
+
+ // renumber cells according to mesh->sortCellsInMEDFileFrmt()
+ renumber( ids1, o2nRenumber[ cz->getLocalDomainNumber() ]);
+ renumber( ids2, o2nRenumber[ cz->getDistantDomainNumber() ]);
+
+ // check nb cell types
+ std::set<INTERP_KERNEL::NormalizedCellType> types1, types2;
+ types1 = mesh1->getTypesOfPart( ids1->begin(), ids1->end() );
+ types2 = mesh2->getTypesOfPart( ids2->begin(), ids2->end() );
+ if ( types1.size() < 1 || types2.size() < 1 )
+ continue; // parallel mode?
+
+ MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
+ for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
+ if ( czVec[j] &&
+ czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
+ czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
+ cz21 = czVec[j];
+
+ if ( types1.size() == 1 && types2.size() == 1 ) // split not needed, only sort
+ {
+ ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
+ cz->setEntityCorresp( *types1.begin(), *types2.begin(),
+ ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
+
+ if ( cz21 )// set 2->1 correspondence
+ {
+ ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
+ cz21->setEntityCorresp( *types2.begin(), *types1.begin(),
+ ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
+ }
+ }
+ else // split and sort
+ {
+ typedef std::pair< std::vector< mcIdType >, std::vector< mcIdType > > T2Vecs;
+ T2Vecs idsByType[ INTERP_KERNEL::NORM_MAXTYPE ][ INTERP_KERNEL::NORM_MAXTYPE ];
+ int t1, t2;
+
+ const mcIdType nbIds = ids1->getNbOfElems();
+ const mcIdType * p1 = ids1->begin(), * p2 = ids2->begin();
+ for ( mcIdType i = 0; i < nbIds; ++i )
+ {
+ t1 = mesh1->getTypeOfCell( p1[ i ]);
+ t2 = mesh2->getTypeOfCell( p2[ i ]);
+ T2Vecs & ids = idsByType[ t1 ][ t2 ];
+ ids.first .push_back( p1[ i ]);
+ ids.second.push_back( p1[ i ]);
+ }
+
+ const int maxType = int( INTERP_KERNEL::NORM_MAXTYPE );
+ for ( t1 = 0; t1 < maxType; ++t1 )
+ for ( t2 = 0; t2 < maxType; ++t2 )
+ {
+ T2Vecs & ids = idsByType[ t1 ][ t2 ];
+ if ( ids.first.empty() ) continue;
+ p1 = & ids.first[0];
+ p2 = & ids.second[0];
+ ids1->desallocate();
+ ids1->pushBackValsSilent( p1, p1+ids.first.size() );
+ ids2->desallocate();
+ ids2->pushBackValsSilent( p2, p2+ids.first.size() );
+ ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
+
+ cz->setEntityCorresp( t1, t2,
+ ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
+
+ if ( cz21 )// set 2->1 correspondence
+ {
+ ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
+ cz21->setEntityCorresp( t2, t1,
+ ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
+ break;
+ }
+ }
+ }// split and sort
+
+ cz->setEntityCorresp( 0, 0, 0, 0 ); // erase ids computed by _topology
+ if ( cz21 )
+ cz21->setEntityCorresp( 0, 0, 0, 0 );
+
+ } // loop on czVec
+
+
+ // ==========================================
+ // 3) sort node ids to be in ascending order
+ // ==========================================
+
+ const bool removeEqual = ( nbInitialDomains > 1 );
+
+ for ( size_t i = 0; i < czVec.size(); ++i )
+ {
+ MEDPARTITIONER::ConnectZone* cz = czVec[i];
+ if ( !cz || cz->getNodeNumber() < 1 )
+ continue;
+ if ( cz->getDistantDomainNumber() < cz->getLocalDomainNumber() )
+ continue; // treat a pair of domains once
+
+ MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
+ for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
+ if ( czVec[j] &&
+ czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
+ czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
+ cz21 = czVec[j];
+
+ // separate ids of two domains
+ const MEDCoupling::MEDCouplingSkyLineArray *corrArray = cz->getNodeCorresp();
+ const DataArrayIdType *ids12 = corrArray->getValuesArray();
+ MCAuto<DataArrayIdType> ids1, ids2, ids12Sorted;
+ ids1 = ids12->selectByTupleIdSafeSlice( 0, corrArray->getLength(), 2 );
+ ids2 = ids12->selectByTupleIdSafeSlice( 1, corrArray->getLength(), 2 );
+
+ ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/0, removeEqual );
+ cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
+
+ if ( cz21 )// set 2->1 correspondence
+ {
+ ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0, false );
+ cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
}
}
}
*/
//================================================================================
-void MEDPARTITIONER::MeshCollection::buildConnectZones()
+void MEDPARTITIONER::MeshCollection::buildBoundaryFaces()
{
+ if (_topology->nbDomain() < 2 || !_subdomain_boundary_creates )
+ return;
+
if ( getMeshDimension() < 2 )
return;
- using ParaMEDMEM::MEDCouplingUMesh;
- using ParaMEDMEM::DataArrayDouble;
- using ParaMEDMEM::DataArrayInt;
+ using MEDCoupling::MEDCouplingUMesh;
+ using MEDCoupling::DataArrayDouble;
+ using MEDCoupling::DataArrayInt;
std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
- int nbMeshes = faceMeshes.size();
+ std::size_t nbMeshes = faceMeshes.size();
//preparing bounding box trees for accelerating search of coincident faces
std::vector<BBTreeOfDim* > bbTrees(nbMeshes);
std::vector<DataArrayDouble*>bbox (nbMeshes);
- for (int inew = 0; inew < nbMeshes-1; inew++)
+ for (unsigned int inew = 0; inew < nbMeshes-1; inew++)
if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
{
- DataArrayDouble* bcCoords = faceMeshes[inew]->getBarycenterAndOwner();
+ DataArrayDouble* bcCoords = faceMeshes[inew]->computeCellCenterOfMass();
bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
bbox[inew]->getConstPointer(),0,0,
}
// loop on domains to find joint faces between them
- for (int inew1 = 0; inew1 < nbMeshes; inew1++ )
+ for (unsigned int inew1 = 0; inew1 < nbMeshes; inew1++ )
{
- for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
+ for (unsigned int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
{
MEDCouplingUMesh* mesh1 = 0;
MEDCouplingUMesh* mesh2 = 0;
bool mesh1Here = true, mesh2Here = true;
if (isParallelMode())
{
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
mesh1Here = _domain_selector->isMyDomain(inew1);
mesh2Here = _domain_selector->isMyDomain(inew2);
if ( !mesh1Here && mesh2Here )
if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
// find coincident faces
- std::vector< int > faces1, faces2;
+ std::vector< mcIdType > faces1, faces2;
if ( mesh1 && mesh2 )
{
- const DataArrayDouble* coords2 = mesh2->getBarycenterAndOwner();
- const double* c2 = coords2->getConstPointer();
- const int dim = coords2->getNumberOfComponents();
- const int nbFaces2 = mesh2->getNumberOfCells();
- const int nbFaces1 = mesh1->getNumberOfCells();
+ const DataArrayDouble* coords2 = mesh2->computeCellCenterOfMass();
+ const double* c2 = coords2->getConstPointer();
+ const std::size_t dim = coords2->getNumberOfComponents();
+ const mcIdType nbFaces2 = mesh2->getNumberOfCells();
+ const mcIdType nbFaces1 = mesh1->getNumberOfCells();
- for (int i2 = 0; i2 < nbFaces2; i2++)
+ for (mcIdType i2 = 0; i2 < nbFaces2; i2++)
{
- std::vector<int> coincFaces;
+ std::vector<mcIdType> coincFaces;
bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
if (coincFaces.size()!=0)
{
- int i1 = coincFaces[0];
+ mcIdType i1 = coincFaces[0];
// if ( coincFaces.size() > 1 )
// {
// i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
if ( isParallelMode())
{
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if ( mesh1Here && !mesh2Here )
{
//send faces2 to domain of recvMesh
// delete bounding box trees
- for (int inew = 0; inew < nbMeshes-1; inew++)
+ for (unsigned int inew = 0; inew < nbMeshes-1; inew++)
if (isParallelMode() && _domain_selector->isMyDomain(inew))
{
bbox[inew]->decrRef();
*/
//================================================================================
-void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces,
+void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< mcIdType >& faces,
const int inew1,
const int inew2,
const bool is2nd )
_group_info.erase( groupName );
// get family IDs array
- int* famIDs = 0;
+ mcIdType* famIDs = 0;
int inew = (is2nd ? inew2 : inew1 );
- int totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
+ mcIdType totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
std::string cle = Cle1ToStr( "faceFamily_toArray", inew );
if ( !_map_dataarray_int.count(cle) )
{
if ( totalNbFaces > 0 )
{
- ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
+ MEDCoupling::DataArrayIdType* p=MEDCoupling::DataArrayIdType::New();
p->alloc( totalNbFaces, 1 );
p->fillWithZero();
famIDs = p->getPointer();
famIDs = _map_dataarray_int.find(cle)->second->getPointer();
}
// find a family ID of an existing JOINT group
- int familyID = 0;
- std::map<std::string, int>::iterator name2id = _family_info.find( groupName );
+ mcIdType familyID = 0;
+ std::map<std::string, mcIdType>::iterator name2id = _family_info.find( groupName );
if ( name2id != _family_info.end() )
familyID = name2id->second;
// remove faces from the familyID-the family
if ( familyID != 0 && famIDs )
- for ( size_t i = 0; i < totalNbFaces; ++i )
+ for ( mcIdType i = 0; i < totalNbFaces; ++i )
if ( famIDs[i] == familyID )
famIDs[i] = 0;
if ( familyID == 0 ) // generate a family ID for JOINT group
{
- std::set< int > familyIDs;
+ std::set< mcIdType > familyIDs;
for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
familyIDs.insert( name2id->second );
// find the next free family ID
_domain_selector( 0 ),
_i_non_empty_mesh(-1),
_driver_type(MEDPARTITIONER::Undefined),
- _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
+ _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
_family_splitting(false),
_create_empty_groups(false),
_joint_finder(0)
throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
}
}
- for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
+ for ( unsigned int idomain = 0; idomain < _mesh.size(); ++idomain )
if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
_i_non_empty_mesh = idomain;
}
_domain_selector( &domainSelector ),
_i_non_empty_mesh(-1),
_driver_type(MEDPARTITIONER::Undefined),
- _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
+ _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
_family_splitting(false),
_create_empty_groups(false),
_joint_finder(0)
</mesh>\n \
</mapping>\n \
</root>\n";
- std::vector<std::string> meshNames=MEDLoader::GetMeshNames(myfile);
+ std::vector<std::string> meshNames=GetMeshNames(myfile);
xml.replace(xml.find("$fileName"),9,myfile);
xml.replace(xml.find("$meshName"),9,meshNames[0]);
xml.replace(xml.find("$meshName"),9,meshNames[0]);
f<<xml;
f.close();
}
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if (MyGlobals::_World_Size>1)
MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
#endif
}
}
// find non-empty domain mesh
- for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
+ for ( unsigned int idomain = 0; idomain < _mesh.size(); ++idomain )
if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
_i_non_empty_mesh = idomain;
try
{
//check for all proc/file compatibility of _field_descriptions
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
_field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
#else
_field_descriptions=MyGlobals::_Field_Descriptions;
std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
}
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
try
{
//check for all proc/file compatibility of _family_info
_i_non_empty_mesh(-1),
_name(meshname),
_driver_type(MEDPARTITIONER::MedXml),
- _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
+ _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
_family_splitting(false),
_create_empty_groups(false),
_joint_finder(0)
MEDPARTITIONER::MeshCollection::~MeshCollection()
{
- for (int i=0; i<(int)_mesh.size();i++)
+ for (std::size_t i=0; i<_mesh.size();i++)
if (_mesh[i]!=0) _mesh[i]->decrRef();
- for (int i=0; i<(int)_cell_family_ids.size();i++)
+ for (std::size_t i=0; i<_cell_family_ids.size();i++)
if (_cell_family_ids[i]!=0)
_cell_family_ids[i]->decrRef();
- for (int i=0; i<(int)_face_mesh.size();i++)
+ for (std::size_t i=0; i<_face_mesh.size();i++)
if (_face_mesh[i]!=0)
_face_mesh[i]->decrRef();
- for (int i=0; i<(int)_face_family_ids.size();i++)
+ for (std::size_t i=0; i<_face_family_ids.size();i++)
if (_face_family_ids[i]!=0)
_face_family_ids[i]->decrRef();
- for (std::map<std::string, ParaMEDMEM::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
+ for (std::map<std::string, MEDCoupling::DataArrayIdType*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
if ((*it).second!=0)
(*it).second->decrRef();
- for (std::map<std::string, ParaMEDMEM::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
+ for (std::map<std::string, MEDCoupling::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
if ((*it).second!=0)
(*it).second->decrRef();
delete _driver;
if (_topology!=0 && _owns_topology)
delete _topology;
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
delete _joint_finder;
#endif
}
*
* The method creates as many MED-files as there are domains in the
* collection. It also creates a master file that lists all the MED files.
- * The MED files created in ths manner contain joints that describe the
+ * The MED files created in this manner contain joints that describe the
* connectivity between subdomains.
*
* \param filename name of the master file that will contain the list of the MED files
*/
void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
{
- //building the connect zones necessary for writing joints
- if (_topology->nbDomain()>1 && _subdomain_boundary_creates )
- buildConnectZones();
//suppresses link with driver so that it can be changed for writing
delete _driver;
_driver=0;
int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
{
int nb=0;
- for (int i=0; i<_mesh.size(); i++)
+ for (size_t i=0; i<_mesh.size(); i++)
{
if (_mesh[i]) nb++;
}
return nb;
}
-int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
+mcIdType MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
{
- int nb=0;
- for (int i=0; i<_mesh.size(); i++)
+ mcIdType nb=0;
+ for (size_t i=0; i<_mesh.size(); i++)
{
if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
}
return nb;
}
-int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
+mcIdType MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
{
- int nb=0;
- for (int i=0; i<_face_mesh.size(); i++)
+ mcIdType nb=0;
+ for (size_t i=0; i<_face_mesh.size(); i++)
{
if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
}
return nb;
}
-std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
+std::vector<MEDCoupling::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
{
return _mesh;
}
-std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
+std::vector<MEDCoupling::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
{
return _face_mesh;
}
-ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
+MEDCoupling::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
{
return _mesh[idomain];
}
-ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
+MEDCoupling::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
{
return _face_mesh[idomain];
}
std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
{
- return _connect_zones;
+ if ( _topology )
+ return _topology->getCZ();
+
+ static std::vector<MEDPARTITIONER::ConnectZone*> noCZ;
+ return noCZ;
}
MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
return _topology;
}
-void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo)
+void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo, bool takeOwneship)
{
if (_topology!=0)
{
throw INTERP_KERNEL::Exception("topology is already set");
}
else
- _topology = topo;
+ {
+ _topology = topo;
+ _owns_topology = takeOwneship;
+ }
}
-/*! Method creating the cell graph in serial mode
- *
- * \param array returns the pointer to the structure that contains the graph
+/*! Method creating the cell graph in serial mode
+ *
+ * \param array returns the pointer to the structure that contains the graph
* \param edgeweight returns the pointer to the table that contains the edgeweights
* (only used if indivisible regions are required)
*/
-void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
+void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
{
using std::map;
using std::pair;
if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
- const ParaMEDMEM::MEDCouplingUMesh* mesh=_mesh[0];
+ const MEDCoupling::MEDCouplingUMesh* mesh=_mesh[0];
if (MyGlobals::_Verbose>50)
std::cout<<"getting nodal connectivity"<<std::endl;
//looking for reverse nodal connectivity i global numbering
if (isParallelMode() && !_domain_selector->isMyDomain(0))
{
- vector<int> value;
- vector<int> index(1,0);
+ vector<mcIdType> value;
+ vector<mcIdType> index(1,0);
- array=new MEDPARTITIONER::SkyLineArray(index,value);
+ array = MEDCoupling::MEDCouplingSkyLineArray::New(index,value);
return;
}
-
- int meshDim = mesh->getMeshDimension();
-
- ParaMEDMEM::DataArrayInt* indexr=ParaMEDMEM::DataArrayInt::New();
- ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
- int nbNodes=mesh->getNumberOfNodes();
- mesh->getReverseNodalConnectivity(revConn,indexr);
- //problem saturation over 1 000 000 nodes for 1 proc
- if (MyGlobals::_Verbose>100)
- std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
- const int* indexr_ptr=indexr->getConstPointer();
- const int* revConn_ptr=revConn->getConstPointer();
-
- const ParaMEDMEM::DataArrayInt* index;
- const ParaMEDMEM::DataArrayInt* conn;
- conn=mesh->getNodalConnectivity();
- index=mesh->getNodalConnectivityIndex();
- int nbCells=mesh->getNumberOfCells();
-
- if (MyGlobals::_Verbose>100)
- std::cout << "proc " << MyGlobals::_Rank << " : getNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
- const int* index_ptr=index->getConstPointer();
- const int* conn_ptr=conn->getConstPointer();
-
- //creating graph arcs (cell to cell relations)
- //arcs are stored in terms of (index,value) notation
- // 0 3 5 6 6
- // 1 2 3 2 3 3
- // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
- // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
-
- //warning here one node have less than or equal effective number of cell with it
- //but cell could have more than effective nodes
- //because other equals nodes in other domain (with other global inode)
- if (MyGlobals::_Verbose>50)
- std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
-
- vector <int> cell2cell_index(nbCells+1,0);
- vector <int> cell2cell;
- cell2cell.reserve(3*nbCells);
-
- for (int icell=0; icell<nbCells;icell++)
- {
- map<int,int > counter;
- for (int iconn=index_ptr[icell]; iconn<index_ptr[icell+1];iconn++)
- {
- int inode=conn_ptr[iconn];
- for (int iconnr=indexr_ptr[inode]; iconnr<indexr_ptr[inode+1];iconnr++)
- {
- int icell2=revConn_ptr[iconnr];
- map<int,int>::iterator iter=counter.find(icell2);
- if (iter!=counter.end()) (iter->second)++;
- else counter.insert(make_pair(icell2,1));
- }
- }
- for (map<int,int>::const_iterator iter=counter.begin();
- iter!=counter.end();
- iter++)
- if (iter->second >= meshDim)
- {
- cell2cell_index[icell+1]++;
- cell2cell.push_back(iter->first);
- }
-
-
- }
- indexr->decrRef();
- revConn->decrRef();
-
- cell2cell_index[0]=0;
- for (int icell=0; icell<nbCells;icell++)
- cell2cell_index[icell+1]=cell2cell_index[icell]+cell2cell_index[icell+1];
-
-
- if (MyGlobals::_Verbose>50)
- std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
-
- //filling up index and value to create skylinearray structure
- array=new MEDPARTITIONER::SkyLineArray(cell2cell_index,cell2cell);
-
- if (MyGlobals::_Verbose>100)
- {
- std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
- cell2cell_index.size()-1 << " " << cell2cell.size() << std::endl;
- int max=cell2cell_index.size()>15?15:cell2cell_index.size();
- if (cell2cell_index.size()>1)
- {
- for (int i=0; i<max; ++i)
- std::cout<<cell2cell_index[i]<<" ";
- std::cout << "... " << cell2cell_index[cell2cell_index.size()-1] << std::endl;
- for (int i=0; i<max; ++i)
- std::cout<< cell2cell[i] << " ";
- int ll=cell2cell_index[cell2cell_index.size()-1]-1;
- std::cout << "... (" << ll << ") " << cell2cell[ll-1] << " " << cell2cell[ll] << std::endl;
- }
- }
-
+ array=mesh->generateGraph();
}
/*! Method creating the cell graph in multidomain mode
*
* \param edgeweight returns the pointer to the table that contains the edgeweights
* (only used if indivisible regions are required)
*/
-void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
+void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
{
using std::multimap;
using std::vector;
using std::make_pair;
using std::pair;
- std::multimap< int, int > node2cell;
- std::map< pair<int,int>, int > cell2cellcounter;
- std::multimap<int,int> cell2cell;
+ std::multimap< mcIdType, mcIdType > node2cell;
+ std::map< pair<mcIdType,mcIdType>, mcIdType > cell2cellcounter;
+ std::multimap<mcIdType,mcIdType> cell2cell;
- std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
+ std::vector<std::vector<std::multimap<mcIdType,mcIdType> > > commonDistantNodes;
int nbdomain=_topology->nbDomain();
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if (isParallelMode())
{
_joint_finder=new JointFinder(*this);
continue;
meshDim = _mesh[idomain]->getMeshDimension();
- ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
- ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
- int nbNodes=_mesh[idomain]->getNumberOfNodes();
+ MEDCoupling::DataArrayIdType* index=MEDCoupling::DataArrayIdType::New();
+ MEDCoupling::DataArrayIdType* revConn=MEDCoupling::DataArrayIdType::New();
+ mcIdType nbNodes=_mesh[idomain]->getNumberOfNodes();
_mesh[idomain]->getReverseNodalConnectivity(revConn,index);
//problem saturation over 1 000 000 nodes for 1 proc
if (MyGlobals::_Verbose>100)
std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
- int* index_ptr=index->getPointer();
- int* revConnPtr=revConn->getPointer();
- for (int i=0; i<nbNodes; i++)
+ mcIdType* index_ptr=index->getPointer();
+ mcIdType* revConnPtr=revConn->getPointer();
+ for (mcIdType i=0; i<nbNodes; i++)
{
- for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
+ for (mcIdType icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
{
- int globalNode=_topology->convertNodeToGlobal(idomain,i);
- int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
+ mcIdType globalNode=_topology->convertNodeToGlobal(idomain,i);
+ mcIdType globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
node2cell.insert(make_pair(globalNode, globalCell));
}
}
revConn->decrRef();
index->decrRef();
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
for (int iother=0; iother<nbdomain; iother++)
{
- std::multimap<int,int>::iterator it;
+ std::multimap<mcIdType,mcIdType>::iterator it;
int isource=idomain;
int itarget=iother;
for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
{
- int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
- int globalCell=(*it).second;
+ mcIdType globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
+ mcIdType globalCell=(*it).second;
node2cell.insert(make_pair(globalNode, globalCell));
}
}
if (MyGlobals::_Verbose>50)
std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
- for (int inode=0;inode<_topology->nbNodes();inode++)
+ for (mcIdType inode=0;inode<_topology->nbNodes();inode++)
{
- typedef multimap<int,int>::const_iterator MI;
+ typedef multimap<mcIdType,mcIdType>::const_iterator MI;
std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
{
- int icell1=cell1->second;
- int icell2=cell2->second;
- if (icell1>icell2) {int tmp=icell1; icell1=icell2; icell2=tmp;}
- std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
+ mcIdType icell1=cell1->second;
+ mcIdType icell2=cell2->second;
+ if (icell1>icell2) std::swap(icell1,icell2);
+ std::map<pair<mcIdType,mcIdType>,mcIdType>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
else (it->second)++;
}
}
- // for (int icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
+ // for (mcIdType icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
// {
- // typedef multimap<int,int>::const_iterator MI;
+ // typedef multimap<int,mcIdType>::const_iterator MI;
// std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
// for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++) //on nodes with icell
// {
// std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
// for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++) //on one of these cell
// {
- // int icell2=cell2->second;
- // std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
+ // mcIdType icell2=cell2->second;
+ // std::map<pair<int,mcIdType>,mcIdType>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
// if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
// else (it->second)++;
// }
//converting the counter to a multimap structure
- for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
+ for (std::map<pair<mcIdType,mcIdType>,mcIdType>::const_iterator it=cell2cellcounter.begin();
it!=cell2cellcounter.end();
it++)
if (it->second>=meshDim)
if (MyGlobals::_Verbose>50)
std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
//filling up index and value to create skylinearray structure
- std::vector <int> index,value;
+ std::vector <mcIdType> index,value;
index.push_back(0);
- int idep=0;
+ mcIdType idep=0;
for (int idomain=0; idomain<nbdomain; idomain++)
{
if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
- int nbCells=_mesh[idomain]->getNumberOfCells();
- for (int icell=0; icell<nbCells; icell++)
+ mcIdType nbCells=_mesh[idomain]->getNumberOfCells();
+ for (mcIdType icell=0; icell<nbCells; icell++)
{
- int size=0;
- int globalCell=_topology->convertCellToGlobal(idomain,icell);
- multimap<int,int>::iterator it;
- pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
+ mcIdType size=0;
+ mcIdType globalCell=_topology->convertCellToGlobal(idomain,icell);
+ multimap<mcIdType,mcIdType>::iterator it;
+ pair<multimap<mcIdType,mcIdType>::iterator,multimap<mcIdType,mcIdType>::iterator> ret;
ret=cell2cell.equal_range(globalCell);
for (it=ret.first; it!=ret.second; ++it)
{
- int ival=(*it).second; //no adding one existing yet
- for (int i=idep ; i<idep+size ; i++)
+ mcIdType ival=(*it).second; //no adding one existing yet
+ for (mcIdType i=idep ; i<idep+size ; i++)
{
if (value[i]==ival)
{
}
}
- array=new MEDPARTITIONER::SkyLineArray(index,value);
+ array=MEDCoupling::MEDCouplingSkyLineArray::New(index,value);
if (MyGlobals::_Verbose>100)
{
std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
index.size()-1 << " " << value.size() << std::endl;
- int max=index.size()>15?15:index.size();
+ std::size_t max=index.size()>15?15:index.size();
if (index.size()>1)
{
- for (int i=0; i<max; ++i)
+ for (std::size_t i=0; i<max; ++i)
std::cout<<index[i]<<" ";
std::cout << "... " << index[index.size()-1] << std::endl;
- for (int i=0; i<max; ++i)
+ for (std::size_t i=0; i<max; ++i)
std::cout<< value[i] << " ";
- int ll=index[index.size()-1]-1;
+ mcIdType ll=index[index.size()-1]-1;
std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
}
}
* returns a topology based on the new graph
*/
MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
- Graph::splitter_type split,
+ Graph::splitter_type split,
const std::string& options_string,
int *user_edge_weights,
int *user_vertices_weights)
{
if (MyGlobals::_Verbose>10)
std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
-
+
if (nbdomain <1)
throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
- MEDPARTITIONER::SkyLineArray* array=0;
+ MEDCoupling::MEDCouplingSkyLineArray* array=0;
int* edgeweights=0;
if (_topology->nbDomain()>1 || isParallelMode())
buildParallelCellGraph(array,edgeweights);
else
buildCellGraph(array,edgeweights);
-
+
Graph* cellGraph = 0;
switch (split)
{
case Graph::METIS:
if ( isParallelMode() && MyGlobals::_World_Size > 1 )
- {
+ {
#ifdef MED_ENABLE_PARMETIS
- if (MyGlobals::_Verbose>10)
- std::cout << "ParMETISGraph" << std::endl;
- cellGraph=new ParMETISGraph(array,edgeweights);
+ if (MyGlobals::_Verbose>10)
+ std::cout << "ParMETISGraph" << std::endl;
+ cellGraph=new ParMETISGraph(array,edgeweights);
#endif
- }
+ }
if ( !cellGraph )
- {
+ {
#ifdef MED_ENABLE_METIS
- if (MyGlobals::_Verbose>10)
- std::cout << "METISGraph" << std::endl;
- cellGraph=new METISGraph(array,edgeweights);
+ if (MyGlobals::_Verbose>10)
+ std::cout << "METISGraph" << std::endl;
+ cellGraph=new METISGraph(array,edgeweights);
#endif
- }
+ }
if ( !cellGraph )
- throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
+ throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available (or you're not running in //). Check your products, please.");
break;
case Graph::SCOTCH:
*/
MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
{
- MEDPARTITIONER::SkyLineArray* array=0;
+ MEDCoupling::MEDCouplingSkyLineArray* array=0;
int* edgeweights=0;
if ( _topology->nbDomain()>1)
Graph* cellGraph;
std::set<int> domains;
- for (int i=0; i<_topology->nbCells(); i++)
+ for (mcIdType i=0; i<_topology->nbCells(); i++)
{
domains.insert(partition[i]);
}
//cellGraph is a shared pointer
Topology *topology=0;
- int nbdomain=domains.size();
+ int nbdomain=(int)domains.size();
topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
// if (array!=0) delete array;
delete cellGraph;
}
}
-ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
+MEDCoupling::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
//getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
-//something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1);
+//something like MCAuto<MEDCouplingFieldDouble> f2=ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1);
{
int rank=MyGlobals::_Rank;
std::string tag="ioldFieldDouble="+IntToStr(iold);
{
if (MyGlobals::_Verbose>300)
std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
- ParaMEDMEM::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
+ MEDCoupling::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
return res;
}
if (MyGlobals::_Verbose>200)
std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
meshName=MyGlobals::_Mesh_Names[iold];
+
+ MCAuto<MEDCoupling::MEDCouplingField> f2Tmp(ReadField((MEDCoupling::TypeOfField) typeField, fileName, meshName, 0, fieldName, DT, IT));
+ MCAuto<MEDCoupling::MEDCouplingFieldDouble> f2(MEDCoupling::DynamicCast<MEDCouplingField,MEDCouplingFieldDouble>(f2Tmp));
- ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
- fileName, meshName, 0, fieldName, DT, IT);
-
- ParaMEDMEM::DataArrayDouble* res=f2->getArray();
+ MEDCoupling::DataArrayDouble* res=f2->getArray();
//to know names of components
std::vector<std::string> browse=BrowseFieldDouble(f2);
std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
MyGlobals::_General_Informations.push_back(localFieldInformation);
res->incrRef();
- f2->decrRef();
_map_dataarray_double[descriptionIold]=res;
return res;
}
//to have unique valid fields names/pointers/descriptions for partitionning
//filter _field_descriptions to be in all procs compliant and equal
{
- int nbfiles=MyGlobals::_File_Names.size(); //nb domains
+ std::size_t nbfiles=MyGlobals::_File_Names.size(); //nb domains
+ if (nbfiles==0)
+ {
+ nbfiles=_topology->nbDomain();
+ }
std::vector<std::string> r2;
//from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
- for (int i=0; i<(int)_field_descriptions.size(); i++)
+ for (std::size_t i=0; i<_field_descriptions.size(); i++)
{
std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
- for (int ii=0; ii<(int)r1.size(); ii++)
+ for (std::size_t ii=0; ii<r1.size(); ii++)
r2.push_back(r1[ii]);
}
//here vector(procs*fields) of serialised vector(description) data
_field_descriptions=r2;
- int nbfields=_field_descriptions.size(); //on all domains
+ std::size_t nbfields=_field_descriptions.size(); //on all domains
if ((nbfields%nbfiles)!=0)
{
if (MyGlobals::_Rank==0)
throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
}
_field_descriptions.resize(nbfields/nbfiles);
- for (int i=0; i<(int)_field_descriptions.size(); i++)
+ for (std::size_t i=0; i<_field_descriptions.size(); i++)
{
std::string str=_field_descriptions[i];
str=EraseTagSerialized(str,"idomain=");
}
//returns true if inodes of a face are in inodes of a cell
-bool isFaceOncell(std::vector< int >& inodesFace, std::vector< int >& inodesCell)
+bool isFaceOncell(std::vector< mcIdType >& inodesFace, std::vector< mcIdType >& inodesCell)
{
- int ires=0;
- int nbok=inodesFace.size();
- for (int i=0; i<nbok; i++)
+ std::size_t ires=0;
+ std::size_t nbok=inodesFace.size();
+ for (std::size_t i=0; i<nbok; i++)
{
- int ii=inodesFace[i];
+ mcIdType ii=inodesFace[i];
if (ii<0)
std::cout << "isFaceOncell problem inodeface<0" << std::endl;
- for (int j=0; j<(int)inodesCell.size(); j++)
+ for (std::size_t j=0; j<inodesCell.size(); j++)
{
if (ii==inodesCell[j])
{
{
if (MyGlobals::_Verbose>200)
std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
- ParaMEDMEM::MEDCouplingUMesh* mcel=_mesh[inew];
- ParaMEDMEM::MEDCouplingUMesh* mfac=_face_mesh[inew];
+ MEDCoupling::MEDCouplingUMesh* mcel=_mesh[inew];
+ MEDCoupling::MEDCouplingUMesh* mfac=_face_mesh[inew];
//to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
- std::vector<int> nodeIds;
+ std::vector<mcIdType> nodeIds;
getNodeIds(*mcel, *mfac, nodeIds);
if (nodeIds.size()==0)
continue; //one empty mesh nothing to do
- ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New();
- ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::DataArrayInt::New();
+ MEDCoupling::DataArrayIdType *revNodalCel=MEDCoupling::DataArrayIdType::New();
+ MEDCoupling::DataArrayIdType *revNodalIndxCel=MEDCoupling::DataArrayIdType::New();
mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
- int *revC=revNodalCel->getPointer();
- int *revIndxC=revNodalIndxCel->getPointer();
+ mcIdType *revC=revNodalCel->getPointer();
+ mcIdType *revIndxC=revNodalIndxCel->getPointer();
- std::vector< int > faceOnCell;
- std::vector< int > faceNotOnCell;
- int nbface=mfac->getNumberOfCells();
- for (int iface=0; iface<nbface; iface++)
+ std::vector< mcIdType > faceOnCell;
+ std::vector< mcIdType > faceNotOnCell;
+ mcIdType nbface=mfac->getNumberOfCells();
+ for (mcIdType iface=0; iface<nbface; iface++)
{
bool ok;
- std::vector< int > inodesFace;
+ std::vector< mcIdType > inodesFace;
mfac->getNodeIdsOfCell(iface, inodesFace);
- int nbnodFace=inodesFace.size();
+ mcIdType nbnodFace=ToIdType(inodesFace.size());
if ( nbnodFace != mfac->getNumberOfNodesInCell( iface ))
continue; // invalid node ids
//set inodesFace in mcel
nbok += (( inodesFace[i]=nodeIds[inodesFace[i]] ) >= 0 );
if ( nbok != nbnodFace )
continue;
- int inod=inodesFace[0];
+ mcIdType inod=inodesFace[0];
if (inod<0)
{
std::cout << "filterFaceOnCell problem 1" << std::endl;
continue;
}
- int nbcell=revIndxC[inod+1]-revIndxC[inod];
- for (int j=0; j<nbcell; j++) //look for each cell with inod
+ mcIdType nbcell=revIndxC[inod+1]-revIndxC[inod];
+ for (mcIdType j=0; j<nbcell; j++) //look for each cell with inod
{
- int icel=revC[revIndxC[inod]+j];
- std::vector< int > inodesCell;
+ mcIdType icel=revC[revIndxC[inod]+j];
+ std::vector< mcIdType > inodesCell;
mcel->getNodeIdsOfCell(icel, inodesCell);
ok=isFaceOncell(inodesFace, inodesCell);
if (ok) break;
if ( faceOnCell.empty() )
_face_mesh[inew] = CreateEmptyMEDCouplingUMesh();
else
- _face_mesh[inew] = (ParaMEDMEM::MEDCouplingUMesh *)
+ _face_mesh[inew] = (MEDCoupling::MEDCouplingUMesh *)
mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true);
mfac->decrRef();
std::string key = Cle1ToStr("faceFamily_toArray",inew);
if ( getMapDataArrayInt().count( key ))
{
- ParaMEDMEM::DataArrayInt * & fam = getMapDataArrayInt()[ key ];
- ParaMEDMEM::DataArrayInt * famFilter = ParaMEDMEM::DataArrayInt::New();
+ MEDCoupling::DataArrayIdType * & fam = getMapDataArrayInt()[ key ];
+ MEDCoupling::DataArrayIdType * famFilter = MEDCoupling::DataArrayIdType::New();
famFilter->alloc(faceOnCell.size(),1);
- int* pfamFilter = famFilter->getPointer();
- int* pfam = fam->getPointer();
+ mcIdType* pfamFilter = famFilter->getPointer();
+ mcIdType* pfam = fam->getPointer();
for ( size_t i=0; i<faceOnCell.size(); i++ )
pfamFilter[i]=pfam[faceOnCell[i]];
fam->decrRef();
