+ 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 );