-// Copyright (C) 2007-2012 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2015 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.
+// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
//
#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 "MEDCouplingAutoRefCountObjectPtr.hxx"
+#include "MEDCouplingFieldDouble.hxx"
#include "MEDCouplingMemArray.hxx"
-#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.hxx"
-#include "MEDCouplingFieldDouble.hxx"
-#include "PointLocator3DIntersectorP0P0.hxx"
-
-#include "MEDCouplingAutoRefCountObjectPtr.hxx"
-#include "BBTree.txx"
+#include "MEDCouplingSkyLineArray.hxx"
+#include "MEDCouplingUMesh.hxx"
+#include "MEDLoader.hxx"
+#include "MEDLoaderBase.hxx"
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
#include <mpi.h>
#endif
-#if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS)
+#ifdef MED_ENABLE_PARMETIS
+#include "MEDPARTITIONER_ParMetisGraph.hxx"
+#endif
+#ifdef MED_ENABLE_METIS
#include "MEDPARTITIONER_MetisGraph.hxx"
#endif
-
#ifdef MED_ENABLE_SCOTCH
#include "MEDPARTITIONER_ScotchGraph.hxx"
#endif
_domain_selector( 0 ),
_i_non_empty_mesh(-1),
_driver_type(MEDPARTITIONER::MedXml),
- _subdomain_boundary_creates(false),
+ _subdomain_boundary_creates( MyGlobals::_Creates_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(false),
+ _subdomain_boundary_creates(MyGlobals::_Creates_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<MEDCoupling::DataArrayInt*> 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
////////////////////
//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
std::cout<<"treating cell and face families"<<std::endl;
castIntField(initialCollection.getMesh(),
- this->getMesh(),
- initialCollection.getCellFamilyIds(),
- "cellFamily");
+ this->getMesh(),
+ initialCollection.getCellFamilyIds(),
+ "cellFamily");
castIntField(initialCollection.getFaceMesh(),
- this->getFaceMesh(),
- initialCollection.getFaceFamilyIds(),
- "faceFamily");
+ 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 (size_t i=0; i<_mesh.size(); 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<int> > >& new2oldIds,
+ std::vector<MEDCoupling::DataArrayInt*> & 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++)
}
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 (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;
if (meshes.size()>1)
{
- ParaMEDMEM::DataArrayInt* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
+ MEDCoupling::DataArrayInt* 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++)
if (splitMeshes[inew][i]!=0)
splitMeshes[inew][i]->decrRef();
- }
+ }
if (MyGlobals::_Verbose>300)
std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl;
}
/*!
- \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)
{
double* bbox;
- BBTree<3>* tree;
+ BBTreeOfDim* tree = 0;
+ int dim = 3;
if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
{
// std::map<pair<double,pair<double, double> >, int > nodeClassifier;
- int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
- bbox=new double[nvertices*2*3];
- 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();
+ bbox=new double[nvertices*2*dim];
- for (int i=0; i<nvertices*3;i++)
+ for (int i=0; i<nvertices*dim;i++)
{
bbox[i*2]=coordsPtr[i]-1e-8;
bbox[i*2+1]=coordsPtr[i]+1e-8;
}
- tree=new BBTree<3>(bbox,0,0,nvertices,1e-9);
+ tree=new BBTreeOfDim( dim, bbox,0,0,nvertices,1e-9);
}
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*3;
+ double* coordsPtr=coords->getPointer()+inode*dim;
vector<int> elems;
tree->getElementsAroundPoint(coordsPtr,elems);
if (elems.size()==0) continue;
}
else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
#else
- if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
+ if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
#endif
- {
- ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords();
- for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
- {
- double* coordsPtr=coords->getPointer()+inode*3;
- vector<int> elems;
- tree->getElementsAroundPoint(coordsPtr,elems);
- if (elems.size()==0) continue;
- nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
- }
- }
+ {
+ MEDCoupling::DataArrayDouble* coords = getMesh(inew)->getCoords();
+ for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
+ {
+ double* coordsPtr=coords->getPointer()+inode*dim;
+ vector<int> 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(iold)))
{
delete tree;
delete[] bbox;
}
- }
+ }
}
-void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
+void getNodeIds(MEDCoupling::MEDCouplingUMesh& meshOne, MEDCoupling::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
{
using std::vector;
+ using MEDPARTITIONER::BBTreeOfDim;
if (!&meshOne || !&meshTwo) return; //empty or not existing
double* bbox;
- BBTree<3>* tree;
+ BBTreeOfDim* tree = 0;
int nv1=meshOne.getNumberOfNodes();
- bbox=new double[nv1*6];
- ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords();
+ MEDCoupling::DataArrayDouble* coords=meshOne.getCoords();
+ int dim = coords->getNumberOfComponents();
+
+ bbox=new double[nv1*2*dim];
double* coordsPtr=coords->getPointer();
- for (int i=0; i<nv1*3; i++)
+ for (int i=0; i<nv1*dim; i++)
{
bbox[i*2]=coordsPtr[i]-1e-8;
bbox[i*2+1]=coordsPtr[i]+1e-8;
}
- tree=new BBTree<3>(bbox,0,0,nv1,1e-9);
-
+ tree=new BBTreeOfDim( dim, bbox,0,0,nv1,1e-9);
+
int nv2=meshTwo.getNumberOfNodes();
nodeIds.resize(nv2,-1);
coords=meshTwo.getCoords();
for (int inode=0; inode<nv2; inode++)
{
- double* coordsPtr2=coords->getPointer()+inode*3;
+ double* coordsPtr2=coords->getPointer()+inode*dim;
vector<int> elems;
tree->getElementsAroundPoint(coordsPtr2,elems);
if (elems.size()==0) continue;
const std::multimap<std::pair<int,int>, std::pair<int,int> >& nodeMapping,
std::vector<std::vector<std::vector<int> > >& new2oldIds)
{
-
//splitMeshes structure will contain the partition of
//the old faces on the new ones
//splitMeshes[4][2] contains the faces from old domain 2
//that have to be added to domain 4
-
+
using std::vector;
using std::map;
using std::multimap;
using std::pair;
using std::make_pair;
-
+
if (MyGlobals::_Verbose>10)
std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl;
if (_topology==0)
throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes");
-
+
int nbNewDomain=_topology->nbDomain();
int nbOldDomain=initialCollection.getTopology()->nbDomain();
-
- vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
- vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
-
- vector< vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
-
+
+ vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
+ vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
+
+ 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);
-
+
//init null pointer for empty meshes
for (int inew=0; inew<nbNewDomain; inew++)
{
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)
std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes empty mesh from iodDomain "<<iold<<std::endl;
}
}
-
-#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++)
{
}
if (!initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
_domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold));
- int nb=0;
- if (splitMeshes[inew][iold])
- nb=splitMeshes[inew][iold]->getNumberOfCells();
+ //int nb=0;
+ //if (splitMeshes[inew][iold])
+ // nb=splitMeshes[inew][iold]->getNumberOfCells();
//std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<<inew<<" "<<iold<<" "<<nb<<std::endl;//" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
}
empty->decrRef();
meshesCastTo.resize(nbNewDomain);
for (int inew=0; inew<nbNewDomain; inew++)
{
- vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
+ vector<const 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);
}
-
+
+ MEDCoupling::MEDCouplingUMesh *bndMesh = 0;
+ if ( _subdomain_boundary_creates &&
+ _mesh[inew] &&
+ _mesh[inew]->getNumberOfCells()>0 )
+ {
+ bndMesh =
+ ((MEDCoupling::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
+ if (bndMesh->getNumberOfCells()>0)
+ myMeshes.push_back( bndMesh );
+ }
+
if (myMeshes.size()>0)
{
- meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
+ meshesCastTo[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshes(myMeshes);
+ meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef();
}
else
{
- ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
+ MEDCoupling::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
meshesCastTo[inew]=empty;
}
for (int iold=0; iold<nbOldDomain; iold++)
if (splitMeshes[inew][iold]!=0)
splitMeshes[inew][iold]->decrRef();
+ if ( bndMesh )
+ bndMesh->decrRef();
}
if (MyGlobals::_Verbose>300)
std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
-void MEDPARTITIONER::MeshCollection::castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
- std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
- std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
- std::string nameArrayTo)
+void MEDPARTITIONER::MeshCollection::castIntField(std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastFrom,
+ std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastTo,
+ std::vector<MEDCoupling::DataArrayInt*>& arrayFrom,
+ std::string nameArrayTo)
{
using std::vector;
//preparing bounding box trees for accelerating source-target node identifications
if (MyGlobals::_Verbose>99)
std::cout<<"making accelerating structures"<<std::endl;
- std::vector<BBTree<3,int>* > acceleratingStructures(ioldMax);
- std::vector<ParaMEDMEM::DataArrayDouble*>bbox(ioldMax);
+ std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
+ std::vector<MEDCoupling::DataArrayDouble*>bbox(ioldMax);
for (int iold =0; iold< ioldMax; iold++)
if (isParallelMode() && _domain_selector->isMyDomain(iold))
{
- ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
+ MEDCoupling::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
- acceleratingStructures[iold]=new BBTree<3,int> (bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
+ acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
+ sourceCoords->decrRef();
}
-
// send-recv operations
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
for (int inew=0; inew<inewMax; inew++)
{
for (int iold=0; iold<ioldMax; iold++)
{
//receive mesh
vector<int> recvIds;
- ParaMEDMEM::MEDCouplingUMesh* recvMesh;
+ 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;
}
void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
- const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
- const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
+ const MEDCoupling::MEDCouplingUMesh& sourceMesh,
+ const MEDCoupling::MEDCouplingUMesh& targetMesh,
const int* fromArray,
std::string nameArrayTo,
- const BBTree<3,int>* myTree)
+ const BBTreeOfDim* myTree)
{
if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
- ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
+ MEDCoupling::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
const double* tc=targetCoords->getConstPointer();
int targetSize=targetMesh.getNumberOfCells();
int sourceSize=sourceMesh.getNumberOfCells();
-if (MyGlobals::_Verbose>200)
-std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
+ if (MyGlobals::_Verbose>200)
+ std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
std::vector<int> ccI;
std::string str,cle;
str=nameArrayTo+"_toArray";
cle=Cle1ToStr(str,inew);
int* toArray;
- const BBTree<3>* tree;
+ const BBTreeOfDim* tree;
bool cleantree=false;
- ParaMEDMEM::DataArrayDouble* sourceBBox=0;
+ MEDCoupling::DataArrayDouble* sourceBBox=0;
+ int dim = targetCoords->getNumberOfComponents();
if (myTree==0)
{
sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8);
- tree=new BBTree<3> (sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
+ tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
cleantree=true;
}
else tree=myTree;
{
if (MyGlobals::_Is0verbose>100)
std::cout << "create " << cle << " size " << targetSize << std::endl;
- ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
+ MEDCoupling::DataArrayInt* p=MEDCoupling::DataArrayInt::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;
+
for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
{
std::vector<int> intersectingElems;
- tree->getElementsAroundPoint(tc+itargetnode*3,intersectingElems); // to be changed in 2D
+ tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
if (intersectingElems.size()!=0)
{
int isourcenode=intersectingElems[0];
- toArray[itargetnode]=fromArray[isourcenode];
- ccI.push_back(itargetnode);
- ccI.push_back(isourcenode);
+ if ( intersectingElems.size() > 1 )
+ {
+ i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
+ isourcenode = intersectingElems[ i2nb->second++ ];
+ }
+ if ( isourcenode < sourceSize ) // protection from invalid elements
+ {
+ toArray[itargetnode]=fromArray[isourcenode];
+ ccI.push_back(itargetnode);
+ ccI.push_back(isourcenode);
+ }
}
}
if (MyGlobals::_Verbose>200)
std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
_map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
-
+
targetCoords->decrRef();
if (cleantree) delete tree;
if (sourceBBox !=0) sourceBBox->decrRef();
-
- }
+}
void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
{
std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
if (descriptionField.find(nameTo)==std::string::npos)
continue; //only nameTo accepted in Fields name description
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
for (int inew=0; inew<inewMax; inew++)
{
for (int iold=0; iold<ioldMax; iold++)
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);
}
}
for (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::DataArrayInt*>::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::DataArrayInt *ccI=it1->second;
if (MyGlobals::_Verbose>300)
std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
" 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();
}
}
+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
+ */
+ //================================================================================
+
+ DataArrayInt* sortCorrespondences( DataArrayInt* ids1,
+ DataArrayInt* ids2,
+ int delta,
+ bool removeEqual = false)
+ {
+ // sort
+ MEDCouplingAutoRefCountObjectPtr< DataArrayInt > renumN2O = ids1->buildPermArrPerLevel();
+ ids1->renumberInPlaceR( renumN2O->begin() );
+ ids2->renumberInPlaceR( renumN2O->begin() );
+
+ if ( removeEqual )
+ {
+ ids1 = ids1->buildUnique();
+ ids2 = ids2->buildUnique();
+ }
+ if ( delta != 0 )
+ {
+ int * id = ids1->getPointer();
+ for ( ; id < ids1->end(); ++id )
+ ++(*id);
+ id = ids2->getPointer();
+ for ( ; id < ids2->end(); ++id )
+ ++(*id);
+ }
+
+ // join
+ DataArrayInt* ids12 = DataArrayInt::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( DataArrayInt* ids, const DataArrayInt* o2nRenumber )
+ {
+ if ( !ids || !o2nRenumber )
+ return;
+ int * id = ids->getPointer();
+ const int * 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::DataArrayInt*> & 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< int > > > 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;
+ int node1 = nmIt1->second.second;
+ for ( ; nmIt2 != nmEnd; ++nmIt2 )
+ {
+ int dom2 = nmIt2->second.first;
+ int 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< int > & 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], 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 DataArrayInt* ids12 = corrArray->getValueArray();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids1, ids2, ids12Sorted;
+ ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 );
+ ids2 = ids12->selectByTupleId2( 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< int >, std::vector< int > > T2Vecs;
+ T2Vecs idsByType[ INTERP_KERNEL::NORM_MAXTYPE ][ INTERP_KERNEL::NORM_MAXTYPE ];
+ int t1, t2;
+
+ const int nbIds = ids1->getNbOfElems();
+ const int * p1 = ids1->begin(), * p2 = ids2->begin();
+ for ( int 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 DataArrayInt *ids12 = corrArray->getValueArray();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids1, ids2, ids12Sorted;
+ ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 );
+ ids2 = ids12->selectByTupleId2( 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 );
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces"
+ * group (where "n" and "p" are domain IDs)
+ */
+//================================================================================
+
+void MEDPARTITIONER::MeshCollection::buildBoundaryFaces()
+{
+ if (_topology->nbDomain() < 2 || !_subdomain_boundary_creates )
+ return;
+
+ if ( getMeshDimension() < 2 )
+ return;
+
+ using MEDCoupling::MEDCouplingUMesh;
+ using MEDCoupling::DataArrayDouble;
+ using MEDCoupling::DataArrayInt;
+
+ std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
+ int 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++)
+ if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
+ {
+ DataArrayDouble* bcCoords = faceMeshes[inew]->getBarycenterAndOwner();
+ bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
+ bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
+ bbox[inew]->getConstPointer(),0,0,
+ bbox[inew]->getNumberOfTuples());
+ bcCoords->decrRef();
+ }
+
+ // loop on domains to find joint faces between them
+ for (int inew1 = 0; inew1 < nbMeshes; inew1++ )
+ {
+ for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
+ {
+ MEDCouplingUMesh* mesh1 = 0;
+ MEDCouplingUMesh* mesh2 = 0;
+ //MEDCouplingUMesh* recvMesh = 0;
+ bool mesh1Here = true, mesh2Here = true;
+ if (isParallelMode())
+ {
+#ifdef HAVE_MPI
+ mesh1Here = _domain_selector->isMyDomain(inew1);
+ mesh2Here = _domain_selector->isMyDomain(inew2);
+ if ( !mesh1Here && mesh2Here )
+ {
+ //send mesh2 to domain of mesh1
+ _domain_selector->sendMesh(*faceMeshes[inew2],
+ _domain_selector->getProcessorID(inew1));
+ }
+ else if ( mesh1Here && !mesh2Here )
+ {
+ //receiving mesh2 from a distant domain
+ _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2));
+ if ( faceMeshes[ inew2 ] )
+ faceMeshes[ inew2 ]->decrRef();
+ faceMeshes[ inew2 ] = mesh2;
+ }
+#endif
+ }
+ if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ];
+ if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
+
+ // find coincident faces
+ std::vector< int > 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();
+
+ for (int i2 = 0; i2 < nbFaces2; i2++)
+ {
+ std::vector<int> coincFaces;
+ bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
+ if (coincFaces.size()!=0)
+ {
+ int i1 = coincFaces[0];
+ // if ( coincFaces.size() > 1 )
+ // {
+ // i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
+ // i1 = coincFaces[ i2nb->second++ ];
+ // }
+ if ( i1 < nbFaces1 ) // protection from invalid elements
+ {
+ faces1.push_back( i1 );
+ faces2.push_back( i2 );
+ }
+ }
+ }
+ coords2->decrRef();
+ }
+
+ if ( isParallelMode())
+ {
+#ifdef HAVE_MPI
+ if ( mesh1Here && !mesh2Here )
+ {
+ //send faces2 to domain of recvMesh
+ SendIntVec(faces2, _domain_selector->getProcessorID(inew2));
+ }
+ else if ( !mesh1Here && mesh2Here )
+ {
+ //receiving ids of faces from a domain of mesh1
+ RecvIntVec(faces2, _domain_selector->getProcessorID(inew1));
+ }
+#endif
+ }
+ // if ( recvMesh )
+ // recvMesh->decrRef();
+
+ // Create group "JOINT_inew1_inew2_Faces" and corresponding families
+ for ( int is2nd = 0; is2nd < 2; ++is2nd )
+ {
+ createJointGroup( is2nd ? faces2 : faces1,
+ inew1 , inew2, is2nd );
+ }
+
+ } // loop on the 2nd domains (inew2)
+ } // loop on the 1st domains (inew1)
+
+
+ // delete bounding box trees
+ for (int inew = 0; inew < nbMeshes-1; inew++)
+ if (isParallelMode() && _domain_selector->isMyDomain(inew))
+ {
+ bbox[inew]->decrRef();
+ delete bbTrees[inew];
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families
+ * \param faces - face ids to include into the group
+ * \param inew1 - index of the 1st domain
+ * \param inew2 - index of the 2nd domain
+ * \param is2nd - in which (1st or 2nd) domain to create the group
+ */
+//================================================================================
+
+void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces,
+ const int inew1,
+ const int inew2,
+ const bool is2nd )
+{
+ // get the name of JOINT group
+ std::string groupName;
+ {
+ std::ostringstream oss;
+ oss << "JOINT_"
+ << (is2nd ? inew2 : inew1 ) << "_"
+ << (is2nd ? inew1 : inew2 ) << "_"
+ << ( getMeshDimension()==2 ? "Edge" : "Face" );
+ groupName = oss.str();
+ }
+
+ // remove existing "JOINT_*" group
+ _group_info.erase( groupName );
+
+ // get family IDs array
+ int* famIDs = 0;
+ int inew = (is2nd ? inew2 : inew1 );
+ int 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 )
+ {
+ MEDCoupling::DataArrayInt* p=MEDCoupling::DataArrayInt::New();
+ p->alloc( totalNbFaces, 1 );
+ p->fillWithZero();
+ famIDs = p->getPointer();
+ _map_dataarray_int[cle]=p;
+ }
+ }
+ else
+ {
+ 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 );
+ if ( name2id != _family_info.end() )
+ familyID = name2id->second;
+
+ // remove faces from the familyID-the family
+ if ( familyID != 0 && famIDs )
+ for ( int i = 0; i < totalNbFaces; ++i )
+ if ( famIDs[i] == familyID )
+ famIDs[i] = 0;
+
+ if ( faces.empty() )
+ return;
+
+ if ( familyID == 0 ) // generate a family ID for JOINT group
+ {
+ std::set< int > familyIDs;
+ for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
+ familyIDs.insert( name2id->second );
+ // find the next free family ID
+ int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd;
+ do
+ {
+ if ( !familyIDs.count( ++familyID ))
+ --freeIdCount;
+ }
+ while ( freeIdCount > 0 );
+ }
+
+ // push faces to familyID-th group
+ if ( faces.back() >= totalNbFaces )
+ throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID");
+ for ( size_t i = 0; i < faces.size(); ++i )
+ famIDs[ faces[i] ] = familyID;
+
+ // register JOINT group and family
+ _family_info[ groupName ] = familyID; // name of the group and family is same
+ _group_info [ groupName ].push_back( groupName );
+}
+
/*! constructing the MESH collection from a distributed file
*
* \param filename name of the master file containing the list of all the MED files
_domain_selector( 0 ),
_i_non_empty_mesh(-1),
_driver_type(MEDPARTITIONER::Undefined),
- _subdomain_boundary_creates(false),
+ _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
_family_splitting(false),
_create_empty_groups(false),
_joint_finder(0)
_domain_selector( &domainSelector ),
_i_non_empty_mesh(-1),
_driver_type(MEDPARTITIONER::Undefined),
- _subdomain_boundary_creates(false),
+ _subdomain_boundary_creates(MyGlobals::_Creates_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.c_str());
+ 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
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(false),
+ _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
_family_splitting(false),
_create_empty_groups(false),
_joint_finder(0)
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::DataArrayInt*>::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
}
*/
void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
{
- //building the connect zones necessary for writing joints
- // if (_topology->nbDomain()>1)
- // 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++;
}
int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() 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=nb+_mesh[i]->getNumberOfCells();
}
int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
{
int nb=0;
- for (int i=0; i<_face_mesh.size(); i++)
+ 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
+/*! 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(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
+{
+
+ using std::map;
+ using std::vector;
+ using std::make_pair;
+ using std::pair;
+
+ if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
+ 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);
+
+ array=new MEDCoupling::MEDCouplingSkyLineArray(index,value);
+ return;
+ }
+ array=mesh->generateGraph();
+}
+/*! Method creating the cell graph in multidomain 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::buildParallelCellGraph(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
{
using std::multimap;
using std::vector;
std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
int nbdomain=_topology->nbDomain();
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
if (isParallelMode())
{
_joint_finder=new JointFinder(*this);
if (MyGlobals::_Verbose>50)
std::cout<<"getting nodal connectivity"<<std::endl;
//looking for reverse nodal connectivity i global numbering
+ int meshDim = 3;
for (int idomain=0; idomain<nbdomain; idomain++)
{
if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
continue;
+ meshDim = _mesh[idomain]->getMeshDimension();
- ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
- ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
+ MEDCoupling::DataArrayInt* index=MEDCoupling::DataArrayInt::New();
+ MEDCoupling::DataArrayInt* revConn=MEDCoupling::DataArrayInt::New();
int nbNodes=_mesh[idomain]->getNumberOfNodes();
_mesh[idomain]->getReverseNodalConnectivity(revConn,index);
//problem saturation over 1 000 000 nodes for 1 proc
}
revConn->decrRef();
index->decrRef();
-#ifdef HAVE_MPI2
+#ifdef HAVE_MPI
for (int iother=0; iother<nbdomain; iother++)
{
std::multimap<int,int>::iterator it;
for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
it!=cell2cellcounter.end();
it++)
- if (it->second>=3)
+ if (it->second>=meshDim)
{
- cell2cell.insert(std::make_pair(it->first.first,it->first.second)); //should be adapted for 2D!
+ cell2cell.insert(std::make_pair(it->first.first,it->first.second));
cell2cell.insert(std::make_pair(it->first.second, it->first.first));
}
index.push_back(idep);
}
}
-
- array=new MEDPARTITIONER::SkyLineArray(index,value);
+
+ array=new MEDCoupling::MEDCouplingSkyLineArray(index,value);
if (MyGlobals::_Verbose>100)
{
* 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;
- buildCellGraph(array,edgeweights);
-
- Graph* cellGraph;
+
+ if (_topology->nbDomain()>1 || isParallelMode())
+ buildParallelCellGraph(array,edgeweights);
+ else
+ buildCellGraph(array,edgeweights);
+
+ Graph* cellGraph = 0;
switch (split)
{
case Graph::METIS:
-#if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS)
- if (MyGlobals::_Verbose>10)
- std::cout << "METISGraph" << std::endl;
- cellGraph=new METISGraph(array,edgeweights);
-#else
- throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
+ if ( isParallelMode() && MyGlobals::_World_Size > 1 )
+ {
+#ifdef MED_ENABLE_PARMETIS
+ 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);
#endif
+ }
+ if ( !cellGraph )
+ throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
break;
+
case Graph::SCOTCH:
#ifdef MED_ENABLE_SCOTCH
if (MyGlobals::_Verbose>10)
*/
MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
{
- MEDPARTITIONER::SkyLineArray* array=0;
+ MEDCoupling::MEDCouplingSkyLineArray* array=0;
int* edgeweights=0;
- buildCellGraph(array,edgeweights);
+ if ( _topology->nbDomain()>1)
+ buildParallelCellGraph(array,edgeweights);
+ else
+ buildCellGraph(array,edgeweights);
+
Graph* cellGraph;
std::set<int> domains;
for (int i=0; i<_topology->nbCells(); i++)
std::ostringstream oss;
oss<<name<<"_"<<i;
if (!isParallelMode() || _domain_selector->isMyDomain(i))
- _mesh[i]->setName(oss.str().c_str());
+ _mesh[i]->setName(oss.str());
}
}
-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.c_str(),f1->getMesh()->getName(),0,f1->getName(),0,1);
+//something like 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)
FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
meshName=MyGlobals::_Mesh_Names[iold];
- ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
- fileName.c_str(), meshName.c_str(), 0, fieldName.c_str(), DT, IT);
+ MEDCoupling::MEDCouplingFieldDouble* f2=ReadField((MEDCoupling::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);
//filter _field_descriptions to be in all procs compliant and equal
{
int 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 (int inew=0; inew<_topology->nbDomain(); inew++)
{
- if (isParallelMode() && _domain_selector->isMyDomain(inew))
+ if (!isParallelMode() || _domain_selector->isMyDomain(inew))
{
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;
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::DataArrayInt *revNodalCel=MEDCoupling::DataArrayInt::New();
+ MEDCoupling::DataArrayInt *revNodalIndxCel=MEDCoupling::DataArrayInt::New();
mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
int *revC=revNodalCel->getPointer();
int *revIndxC=revNodalIndxCel->getPointer();
-
+
std::vector< int > faceOnCell;
std::vector< int > faceNotOnCell;
int nbface=mfac->getNumberOfCells();
std::vector< int > inodesFace;
mfac->getNodeIdsOfCell(iface, inodesFace);
int nbnodFace=inodesFace.size();
+ if ( nbnodFace != mfac->getNumberOfNodesInCell( iface ))
+ continue; // invalid node ids
//set inodesFace in mcel
- for (int i=0; i<nbnodFace; i++) inodesFace[i]=nodeIds[inodesFace[i]];
+ int nbok = 0;
+ for (int i=0; i<nbnodFace; i++)
+ nbok += (( inodesFace[i]=nodeIds[inodesFace[i]] ) >= 0 );
+ if ( nbok != nbnodFace )
+ continue;
int inod=inodesFace[0];
if (inod<0)
- std::cout << "filterFaceOnCell problem 1" << std::endl;
+ {
+ 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
{
std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
}
}
-
+
revNodalCel->decrRef();
revNodalIndxCel->decrRef();
-
- std::string keyy;
- keyy=Cle1ToStr("filterFaceOnCell",inew);
- _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
- keyy=Cle1ToStr("filterNotFaceOnCell",inew);
- _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
+
+ // std::string keyy;
+ // keyy=Cle1ToStr("filterFaceOnCell",inew);
+ // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
+ // keyy=Cle1ToStr("filterNotFaceOnCell",inew);
+ // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
+
+ // filter the face mesh
+ if ( faceOnCell.empty() )
+ _face_mesh[inew] = CreateEmptyMEDCouplingUMesh();
+ else
+ _face_mesh[inew] = (MEDCoupling::MEDCouplingUMesh *)
+ mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true);
+ mfac->decrRef();
+
+ // filter the face families
+ std::string key = Cle1ToStr("faceFamily_toArray",inew);
+ if ( getMapDataArrayInt().count( key ))
+ {
+ MEDCoupling::DataArrayInt * & fam = getMapDataArrayInt()[ key ];
+ MEDCoupling::DataArrayInt * famFilter = MEDCoupling::DataArrayInt::New();
+ famFilter->alloc(faceOnCell.size(),1);
+ int* pfamFilter = famFilter->getPointer();
+ int* pfam = fam->getPointer();
+ for ( size_t i=0; i<faceOnCell.size(); i++ )
+ pfamFilter[i]=pfam[faceOnCell[i]];
+ fam->decrRef();
+ fam = famFilter;
+ }
}
}
}
