+namespace
+{
+ using namespace ParaMEDMEM;
+ //================================================================================
+ /*!
+ * \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<ParaMEDMEM::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;
+ ParaMEDMEM::MEDCouplingUMesh* mesh1 = _mesh[ cz->getLocalDomainNumber () ];
+ ParaMEDMEM::MEDCouplingUMesh* mesh2 = _mesh[ cz->getDistantDomainNumber() ];
+
+ // separate ids of two domains
+ const ParaMEDMEM::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 ParaMEDMEM::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 ParaMEDMEM::MEDCouplingUMesh;
+ using ParaMEDMEM::DataArrayDouble;
+ using ParaMEDMEM::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 )
+ {
+ ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::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 );
+}
+
