netgen/libsrc/meshing/paralleltop.cpp
2009-01-25 19:37:14 +00:00

1702 lines
40 KiB
C++

#ifdef PARALLEL
#include <meshing.hpp>
#include "paralleltop.hpp"
namespace netgen
{
void ParallelMeshTopology :: Reset ()
{
*testout << "ParallelMeshTopology::Reset" << endl;
if ( ntasks == 1 ) return;
PrintMessage ( 4, "RESET");
int nvold = nv;
int nedold = ned;
int nfaold = nfa;
ne = mesh.GetNE();
nv = mesh.GetNV();
nseg = mesh.GetNSeg();
nsurfel = mesh.GetNSE();
ned = mesh.GetTopology().GetNEdges();
nfa = mesh.GetTopology().GetNFaces();
loc2distedge.ChangeSize (ned);
for (int i = 0; i < ned; i++)
if (loc2distedge[i].Size() == 0)
loc2distedge.Add (i, -1); // will be the global nr
loc2distface.ChangeSize (nfa);
for (int i = 0; i < nfa; i++)
if (loc2distface[i].Size() == 0)
loc2distface.Add (i, -1); // will be the global nr
if ( !isexchangevert )
{
isexchangevert = new BitArray (nv * ( ntasks+1 ));
isexchangevert->Clear();
}
if ( !isexchangeedge )
{
isexchangeedge = new BitArray (ned*(ntasks+1) );
isexchangeedge->Clear();
}
if ( !isexchangeface )
{
isexchangeface = new BitArray (nfa*(ntasks+1) );
isexchangeface->Clear();
}
if ( !isexchangeel )
{
isexchangeel = new BitArray (ne*(ntasks+1) );
isexchangeel->Clear();
}
// if the number of vertices did not change, return
if ( nvold == nv ) return;
// faces and edges get new numbers -> delete
isexchangeface -> SetSize(nfa*(ntasks+1) );
isexchangeedge -> SetSize(ned*(ntasks+1) );
isexchangeface -> Clear();
isexchangeedge -> Clear();
SetNV(nv);
SetNE(ne);
if ( !isghostedge.Size() )
{
isghostedge.SetSize(ned);
isghostedge.Clear();
}
if ( !isghostface.Size() )
{
isghostface.SetSize(nfa);
isghostface.Clear();
}
}
ParallelMeshTopology :: ~ParallelMeshTopology ()
{
delete isexchangeface;
delete isexchangevert;
delete isexchangeedge;
delete isexchangeel;
}
ParallelMeshTopology :: ParallelMeshTopology ( const netgen::Mesh & amesh )
: mesh(amesh)
{
ned = 0; //mesh.GetTopology().GetNEdges();
nfa = 0; //mesh.GetTopology().GetNFaces();
nv = 0;
ne = 0;
np = 0;
nseg = 0;
nsurfel = 0;
neglob = 0;
nvglob = 0;
nparel = 0;
isexchangeface = 0;
isexchangevert = 0;
isexchangeel = 0;
isexchangeedge = 0;
coarseupdate = 0;
isghostedge.SetSize(0);
isghostface.SetSize(0);
overlap = 0;
}
int ParallelMeshTopology :: Glob2Loc_Vert (int globnum )
{
for (int i = 1; i <= nv; i++)
if ( globnum == loc2distvert[i][0] )
return i;
return -1;
}
int ParallelMeshTopology :: Glob2Loc_VolEl (int globnum )
{
int locnum = -1;
for (int i = 0; i < ne; i++)
{
if ( globnum == loc2distel[i][0] )
{
locnum = i+1;
}
}
return locnum;
}
int ParallelMeshTopology :: Glob2Loc_SurfEl (int globnum )
{
int locnum = -1;
for (int i = 0; i < nsurfel; i++)
{
if ( globnum == loc2distsurfel[i][0] )
{
locnum = i+1;
}
}
return locnum;
}
int ParallelMeshTopology :: Glob2Loc_Segm (int globnum )
{
int locnum = -1;
for (int i = 0; i < nseg; i++)
{
if ( globnum == loc2distsegm[i][0] )
{
locnum = i+1;
}
}
return locnum;
}
void ParallelMeshTopology :: Print() const
{
(*testout) << endl << "TOPOLOGY FOR PARALLEL MESHES" << endl << endl;
for ( int i = 1; i <= nv; i++ )
if ( IsExchangeVert (i) )
{
(*testout) << "exchange point " << i << ": global " << GetLoc2Glob_Vert(i) << endl;
for ( int dest = 0; dest < ntasks; dest ++)
if ( dest != id )
if ( GetDistantPNum( dest, i ) > 0 )
(*testout) << " p" << dest << ": " << GetDistantPNum ( dest, i ) << endl;
}
for ( int i = 1; i <= ned; i++ )
if ( IsExchangeEdge ( i ) )
{
int v1, v2;
mesh . GetTopology().GetEdgeVertices(i, v1, v2);
(*testout) << "exchange edge " << i << ": global vertices " << GetLoc2Glob_Vert(v1) << " "
<< GetLoc2Glob_Vert(v2) << endl;
for ( int dest = 0; dest < ntasks; dest++)
if ( GetDistantEdgeNum ( dest, i ) > 0 )
if ( dest != id )
{
(*testout) << " p" << dest << ": " << GetDistantEdgeNum ( dest, i ) << endl;
}
}
for ( int i = 1; i <= nfa; i++ )
if ( IsExchangeFace(i) )
{
Array<int> facevert;
mesh . GetTopology().GetFaceVertices(i, facevert);
(*testout) << "exchange face " << i << ": global vertices " ;
for ( int fi=0; fi < facevert.Size(); fi++)
(*testout) << GetLoc2Glob_Vert(facevert[fi]) << " ";
(*testout) << endl;
for ( int dest = 0; dest < ntasks; dest++)
if ( dest != id )
{
if ( GetDistantFaceNum ( dest, i ) >= 0 )
(*testout) << " p" << dest << ": " << GetDistantFaceNum ( dest, i ) << endl;
}
}
for ( int i = 1; i < mesh.GetNE(); i++)
{
if ( !IsExchangeElement(i) ) continue;
Array<int> vert;
const Element & el = mesh.VolumeElement(i);
(*testout) << "parallel local element " << i << endl;
(*testout) << "vertices " ;
for ( int j = 0; j < el.GetNV(); j++)
(*testout) << el.PNum(j+1) << " ";
(*testout) << "is ghost " << IsGhostEl(i) << endl;
(*testout) << endl;
}
}
int ParallelMeshTopology :: GetDistantPNum ( int proc, int locpnum ) const
{
if ( proc == 0 )
return loc2distvert[locpnum][0];
for (int i = 1; i < loc2distvert[locpnum].Size(); i += 2)
if ( loc2distvert[locpnum][i] == proc )
return loc2distvert[locpnum][i+1];
return -1;
}
int ParallelMeshTopology :: GetDistantFaceNum ( int proc, int locfacenum ) const
{
if ( proc == 0 )
return loc2distface[locfacenum-1][0];
for ( int i = 1; i < loc2distface[locfacenum-1].Size(); i+=2 )
if ( loc2distface[locfacenum-1][i] == proc )
return loc2distface[locfacenum-1][i+1];
return -1;
}
int ParallelMeshTopology :: GetDistantEdgeNum ( int proc, int locedgenum ) const
{
if ( proc == 0 )
return loc2distedge[locedgenum-1][0];
for ( int i = 1; i < loc2distedge[locedgenum-1].Size(); i+=2 )
if ( loc2distedge[locedgenum-1][i] == proc )
return loc2distedge[locedgenum-1][i+1];
return -1;
}
int ParallelMeshTopology :: GetDistantElNum ( int proc, int locelnum ) const
{
if ( proc == 0 )
return loc2distel[locelnum-1][0];
for ( int i = 1; i < loc2distel[locelnum-1].Size(); i+=2 )
if ( loc2distel[locelnum-1][i] == proc )
return loc2distel[locelnum-1][i+1];
return -1;
}
/*
//
// gibt anzahl an distant pnums zurueck
int ParallelMeshTopology :: GetNDistantPNums ( int locpnum ) const
{
return loc2distvert[locpnum].Size() / 2 + 1;
}
int ParallelMeshTopology :: GetNDistantFaceNums ( int locfacenum ) const
{
int size = loc2distface[locfacenum-1].Size() / 2 + 1;
return size;
}
int ParallelMeshTopology :: GetNDistantEdgeNums ( int locedgenum ) const
{
int size = loc2distedge[locedgenum-1].Size() / 2 + 1;
return size;
}
int ParallelMeshTopology :: GetNDistantElNums ( int locelnum ) const
{
int size = loc2distel[locelnum-1].Size() / 2 + 1;
return size;
}
*/
// gibt anzahl an distant pnums zurueck
// * pnums entspricht Array<int[2] >
int ParallelMeshTopology :: GetDistantPNums ( int locpnum, int * distpnums ) const
{
// distpnums[0] = loc2distvert[locpnum][0];
// for (int i = 1; i < loc2distvert[locpnum].Size(); i += 2)
// distpnums[ loc2distvert[locpnum][i] ] = loc2distvert[locpnum][i+1];
distpnums[0] = 0;
distpnums[1] = loc2distvert[locpnum][0];
for ( int i = 1; i < loc2distvert[locpnum].Size(); i++ )
distpnums[i+1] = loc2distvert[locpnum][i];
int size = loc2distvert[locpnum].Size() / 2 + 1;
return size;
}
int ParallelMeshTopology :: GetDistantFaceNums ( int locfacenum, int * distfacenums ) const
{
// distfacenums[0] = loc2distface[locfacenum-1][0];
// for ( int i = 1; i < loc2distface[locfacenum-1].Size(); i+=2 )
// distfacenums[loc2distface[locfacenum-1][i]] = loc2distface[locfacenum-1][i+1];
distfacenums[0] = 0;
distfacenums[1] = loc2distface[locfacenum-1][0];
for ( int i = 1; i < loc2distface[locfacenum-1].Size(); i++ )
distfacenums[i+1] = loc2distface[locfacenum-1][i];
int size = loc2distface[locfacenum-1].Size() / 2 + 1;
return size;
}
int ParallelMeshTopology :: GetDistantEdgeNums ( int locedgenum, int * distedgenums ) const
{
// distedgenums[0] = loc2distedge[locedgenum-1][0];
// for ( int i = 1; i < loc2distedge[locedgenum-1].Size(); i+=2 )
// distedgenums[loc2distedge[locedgenum-1][i]] = loc2distedge[locedgenum-1][i+1];
distedgenums[0] = 0;
distedgenums[1] = loc2distedge[locedgenum-1][0];
for ( int i = 1; i < loc2distedge[locedgenum-1].Size(); i++ )
distedgenums[i+1] = loc2distedge[locedgenum-1][i];
int size = loc2distedge[locedgenum-1].Size() / 2 + 1;
return size;
}
int ParallelMeshTopology :: GetDistantElNums ( int locelnum, int * distelnums ) const
{
// distelnums[0] = loc2distel[locelnum-1][0];
// for ( int i = 1; i < loc2distel[locelnum-1].Size(); i+=2 )
// distelnums[loc2distel[locelnum-1][i]] = loc2distel[locelnum-1][i+1];
distelnums[0] = 0;
distelnums[1] = loc2distel[locelnum-1][0];
for ( int i = 1; i < loc2distel[locelnum-1].Size(); i++ )
distelnums[i+1] = loc2distel[locelnum-1][i];
int size = loc2distel[locelnum-1].Size() / 2 + 1;
return size;
}
void ParallelMeshTopology :: SetDistantFaceNum ( int dest, int locnum, int distnum )
{
if ( dest == 0 )
{
loc2distface[locnum-1][0] = distnum;
return;
}
for ( int i = 1; i < loc2distface[locnum-1].Size(); i+=2 )
if ( loc2distface[locnum-1][i] == dest )
{
loc2distface[locnum-1][i+1] = distnum;
return;
}
loc2distface.Add(locnum-1, dest);
loc2distface.Add(locnum-1, distnum);
}
void ParallelMeshTopology :: SetDistantPNum ( int dest, int locnum, int distnum )
{
if ( dest == 0 )
{
loc2distvert[locnum][0] = distnum; // HERE
return;
}
for ( int i = 1; i < loc2distvert[locnum].Size(); i+=2 )
if ( loc2distvert[locnum][i] == dest )
{
loc2distvert[locnum][i+1] = distnum;
return;
}
loc2distvert.Add (locnum, dest);
loc2distvert.Add (locnum, distnum);
}
void ParallelMeshTopology :: SetDistantEdgeNum ( int dest, int locnum, int distnum )
{
if ( dest == 0 )
{
loc2distedge[locnum-1][0] = distnum;
return;
}
for ( int i = 1; i < loc2distedge[locnum-1].Size(); i+=2 )
if ( loc2distedge[locnum-1][i] == dest )
{
loc2distedge[locnum-1][i+1] = distnum;
return;
}
loc2distedge.Add (locnum-1, dest);
loc2distedge.Add (locnum-1, distnum);
}
void ParallelMeshTopology :: SetDistantEl ( int dest, int locnum, int distnum )
{
if ( dest == 0 )
{
loc2distel[locnum-1][0] = distnum;
return;
}
for ( int i = 1; i < loc2distel[locnum-1].Size(); i+=2 )
if ( loc2distel[locnum-1][i] == dest )
{
loc2distel[locnum-1][i+1] = distnum;
return;
}
loc2distel.Add (locnum-1, dest);
loc2distel.Add (locnum-1, distnum);
}
void ParallelMeshTopology :: SetDistantSurfEl ( int dest, int locnum, int distnum )
{
if ( dest == 0 )
{
loc2distsurfel[locnum-1][0] = distnum;
return;
}
for ( int i = 1; i < loc2distsurfel[locnum-1].Size(); i+=2 )
if ( loc2distsurfel[locnum-1][i] == dest )
{
loc2distsurfel[locnum-1][i+1] = distnum;
return;
}
loc2distsurfel.Add (locnum-1, dest);
loc2distsurfel.Add (locnum-1, distnum);
}
void ParallelMeshTopology :: SetDistantSegm ( int dest, int locnum, int distnum )
{
if ( dest == 0 )
{
loc2distsegm[locnum-1][0] = distnum;
return;
}
for (int i = 1; i < loc2distsegm[locnum-1].Size(); i+=2 )
if ( loc2distsegm[locnum-1][i] == dest )
{
loc2distsegm[locnum-1][i+1] = distnum;
return;
}
loc2distsegm.Add (locnum-1, dest);
loc2distsegm.Add (locnum-1, distnum);
}
void ParallelMeshTopology :: GetVertNeighbours ( int vnum, Array<int> & dests ) const
{
dests.SetSize(0);
int i = 1;
while ( i < loc2distvert[vnum].Size() )
{
dests.Append ( loc2distvert[vnum][i] );
i+=2;
}
}
void ParallelMeshTopology :: Update ()
{
ne = mesh.GetNE();
nv = mesh.GetNV();
nseg = mesh.GetNSeg();
nsurfel = mesh.GetNSE();
ned = mesh.GetTopology().GetNEdges();
nfa = mesh.GetTopology().GetNFaces();
}
void ParallelMeshTopology :: UpdateRefinement ()
{
;
}
void ParallelMeshTopology :: UpdateCoarseGridGlobal ()
{
PrintMessage ( 1, "UPDATE GLOBAL COARSEGRID STARTS" ); // JS
// MPI_Barrier (MPI_COMM_WORLD);
// PrintMessage ( 1, "all friends are here " ); // JS
// MPI_Barrier (MPI_COMM_WORLD);
int timer = NgProfiler::CreateTimer ("UpdateCoarseGridGlobal");
NgProfiler::RegionTimer reg(timer);
*testout << "ParallelMeshTopology :: UpdateCoarseGridGlobal" << endl;
const MeshTopology & topology = mesh.GetTopology();
Array<int> sendarray, recvarray;
nfa = topology . GetNFaces();
ned = topology . GetNEdges();
np = mesh . GetNP();
nv = mesh . GetNV();
ne = mesh . GetNE();
nseg = mesh.GetNSeg();
nsurfel = mesh.GetNSE();
// low order processor - save mesh partition
if ( id == 0 )
{
if ( !isexchangeel )
{
isexchangeel = new BitArray ( (ntasks+1) * ne );
isexchangeel -> Clear();
}
for ( int eli = 1; eli <= ne; eli++ )
{
loc2distel[eli-1][0] = eli;
SetExchangeElement ( eli );
const Element & el = mesh . VolumeElement ( eli );
int dest = el . GetPartition ( );
SetExchangeElement ( dest, eli );
for ( int i = 0; i < el.GetNP(); i++ )
{
SetExchangeVert ( dest, el.PNum(i+1) );
SetExchangeVert ( el.PNum(i+1) );
}
Array<int> edges;
topology . GetElementEdges ( eli, edges );
for ( int i = 0; i < edges.Size(); i++ )
{
SetExchangeEdge ( dest, edges[i] );
SetExchangeEdge ( edges[i] );
}
topology . GetElementFaces ( eli, edges );
for ( int i = 0; i < edges.Size(); i++ )
{
SetExchangeFace ( dest, edges[i] );
SetExchangeFace ( edges[i] );
}
}
// HERE
for ( int i = 1; i <= mesh .GetNV(); i++)
loc2distvert[i][0] = i;
for ( int i = 0; i < mesh . GetNSeg(); i++)
loc2distsegm[i][0] = i+1;
for ( int i = 0; i < mesh . GetNSE(); i++)
loc2distsurfel[i][0] = i+1;
for ( int i = 0; i < topology .GetNEdges(); i++)
loc2distedge[i][0] = i+1;
for ( int i = 0; i < topology .GetNFaces(); i++)
loc2distface[i][0] = i+1;
}
if ( id == 0 )
sendarray.Append (nfa);
BitArray recvface(nfa);
recvface.Clear();
/*
Array<int> edges, pnums, faces;
for ( int el = 1; el <= ne; el++ )
{
topology.GetElementFaces (el, faces);
int globeli = GetLoc2Glob_VolEl(el);
for ( int fai = 0; fai < faces.Size(); fai++)
{
int fa = faces[fai];
topology.GetElementEdges ( el, edges );
topology.GetFaceVertices ( fa, pnums );
// send :
// localfacenum
// np
// ned
// globalpnums
// localpnums
// localedgenums mit globalv1, globalv2
sendarray. Append ( fa );
sendarray. Append ( globeli );
sendarray. Append ( pnums.Size() );
sendarray. Append ( edges.Size() );
if (id == 0)
for ( int i = 0; i < pnums.Size(); i++ )
sendarray. Append( pnums[i] );
else
for ( int i = 0; i < pnums.Size(); i++ )
sendarray. Append( GetLoc2Glob_Vert(pnums[i]) );
for ( int i = 0; i < pnums.Size(); i++ )
sendarray. Append(pnums[i] );
for ( int i = 0; i < edges.Size(); i++ )
{
sendarray. Append(edges[i] );
int v1, v2;
topology . GetEdgeVertices ( edges[i], v1, v2 );
int dv1 = GetLoc2Glob_Vert ( v1 );
int dv2 = GetLoc2Glob_Vert ( v2 );
if (id > 0) if ( dv1 > dv2 ) swap ( dv1, dv2 );
sendarray . Append ( dv1 );
sendarray . Append ( dv2 );
}
}
}
*/
// new version
Array<int> edges, pnums, faces, elpnums;
sendarray.Append (ne);
for ( int el = 1; el <= ne; el++ )
{
topology.GetElementFaces (el, faces);
topology.GetElementEdges ( el, edges );
const Element & volel = mesh.VolumeElement (el);
int globeli = GetLoc2Glob_VolEl(el);
sendarray. Append ( globeli );
sendarray. Append ( faces.Size() );
sendarray. Append ( edges.Size() );
sendarray. Append ( volel.GetNP() );
for ( int i = 0; i < faces.Size(); i++ )
sendarray. Append(faces[i] );
for ( int i = 0; i < edges.Size(); i++ )
sendarray. Append(edges[i] );
for ( int i = 0; i < volel.GetNP(); i++ )
if (id == 0)
sendarray. Append(volel[i] );
else
sendarray. Append(GetLoc2Glob_Vert (volel[i]));
}
// end new version
BitArray edgeisinit(ned), vertisinit(np);
edgeisinit.Clear();
vertisinit.Clear();
// Array for temporary use, to find local from global element fast
Array<int,1> glob2loc_el;
if ( id != 0 )
{
glob2loc_el.SetSize (neglob);
glob2loc_el = -1;
for ( int locel = 1; locel <= mesh.GetNE(); locel++)
glob2loc_el[GetLoc2Glob_VolEl(locel)] = locel;
}
// MPI_Barrier (MPI_COMM_WORLD);
MPI_Request sendrequest;
if (id == 0)
{
PrintMessage (4, "UpdateCoarseGridGlobal : bcast, size = ", int (sendarray.Size()*sizeof(int)) );
MyMPI_Bcast ( sendarray );
}
else
MyMPI_ISend ( sendarray, 0, sendrequest );
int nloops = (id == 0) ? ntasks-1 : 1;
for (int hi = 0; hi < nloops; hi++)
{
int sender;
if (id == 0)
{
sender = MyMPI_Recv ( recvarray );
PrintMessage (4, "have received from ", sender);
}
else
{
MyMPI_Bcast ( recvarray );
sender = 0;
}
// compare received vertices with own ones
int ii = 0;
int cntel = 0;
int volel = 1;
if ( id != 0 )
nfaglob = recvarray[ii++];
Array<int> faces, edges;
Array<int> pnums, globalpnums;
int recv_ne = recvarray[ii++];
for (int hi = 0; hi < recv_ne; hi++)
{
int globvolel = recvarray[ii++];
int distnfa = recvarray[ii++];
int distned = recvarray[ii++];
int distnp = recvarray[ii++];
if ( id > 0 )
volel = glob2loc_el[globvolel];
else
volel = globvolel;
if (volel != -1)
{
topology.GetElementFaces( volel, faces);
topology.GetElementEdges ( volel, edges);
const Element & volelement = mesh.VolumeElement (volel);
for ( int i = 0; i < faces.Size(); i++)
SetDistantFaceNum ( sender, faces[i], recvarray[ii++]);
for ( int i = 0; i < edges.Size(); i++)
SetDistantEdgeNum ( sender, edges[i], recvarray[ii++]);
for ( int i = 0; i < distnp; i++)
SetDistantPNum ( sender, volelement[i], recvarray[ii++]);
}
else
ii += distnfa + distned + distnp;
}
}
coarseupdate = 1;
if (id != 0)
{
MPI_Status status;
MPI_Wait (&sendrequest, &status);
}
#ifdef SCALASCA
#pragma pomp inst end(updatecoarsegrid)
#endif
}
void ParallelMeshTopology :: UpdateCoarseGrid ()
{
int timer = NgProfiler::CreateTimer ("UpdateCoarseGrid");
NgProfiler::RegionTimer reg(timer);
#ifdef SCALASCA
#pragma pomp inst begin(updatecoarsegrid)
#endif
(*testout) << "UPDATE COARSE GRID PARALLEL TOPOLOGY " << endl;
PrintMessage ( 1, "UPDATE COARSE GRID PARALLEL TOPOLOGY " );
// find exchange edges - first send exchangeedges locnum, v1, v2
// receive distant distnum, v1, v2
// find matching
const MeshTopology & topology = mesh.GetTopology();
UpdateCoarseGridGlobal();
if ( id == 0 ) return;
Array<int> * sendarray, *recvarray;
sendarray = new Array<int> (0);
recvarray = new Array<int>;
nfa = topology . GetNFaces();
ned = topology . GetNEdges();
np = mesh . GetNP();
nv = mesh . GetNV();
ne = mesh . GetNE();
nseg = mesh.GetNSeg();
nsurfel = mesh.GetNSE();
sendarray->SetSize (0);
BitArray recvface(nfa);
recvface.Clear();
for ( int fa = 1; fa <= nfa; fa++ )
{
if ( !IsExchangeFace ( fa ) ) continue;
Array<int> edges, pnums;
int globfa = GetDistantFaceNum ( 0, fa );
topology.GetFaceEdges ( fa, edges );
topology.GetFaceVertices ( fa, pnums );
// send :
// localfacenum
// globalfacenum
// np
// ned
// globalpnums
// localpnums
// localedgenums mit globalv1, globalv2
//
sendarray -> Append ( fa );
sendarray -> Append ( globfa );
sendarray -> Append ( pnums.Size() );
sendarray -> Append ( edges.Size() );
for ( int i = 0; i < pnums.Size(); i++ )
{
sendarray -> Append( GetLoc2Glob_Vert(pnums[i]) );
}
for ( int i = 0; i < pnums.Size(); i++ )
{
sendarray -> Append(pnums[i] );
}
for ( int i = 0; i < edges.Size(); i++ )
{
sendarray -> Append(edges[i] );
int v1, v2;
topology . GetEdgeVertices ( edges[i], v1, v2 );
int dv1 = GetLoc2Glob_Vert ( v1 );
int dv2 = GetLoc2Glob_Vert ( v2 );
sendarray -> Append ( dv1 );
sendarray -> Append ( dv2 );
}
}
BitArray edgeisinit(ned), vertisinit(np);
edgeisinit.Clear();
vertisinit.Clear();
// Array for temporary use, to find local from global element fast
// only for not too big meshes
// seems ok, as low-order space is treated on one proc
Array<int,1> * glob2locfa;
glob2locfa = new Array<int,1> ( nfaglob );
(*glob2locfa) = -1;
for ( int locfa = 1; locfa <= nfa; locfa++)
{
if ( !IsExchangeFace ( locfa ) ) continue;
(*glob2locfa)[GetDistantFaceNum(0, locfa) ] = locfa;
}
for ( int sender = 1; sender < ntasks; sender ++ )
{
if ( id == sender )
MyMPI_Bcast ( *sendarray, sender-1, MPI_HIGHORDER_COMM);
if ( id != sender )
{
MyMPI_Bcast ( *recvarray, sender-1, MPI_HIGHORDER_COMM);
// compare received vertices with own ones
int ii = 0;
int cntel = 0;
int locfa = 1;
while ( ii< recvarray -> Size() )
{
// receive list :
// distant facenum
// global facenum
// np
// ned
// globalpnums
// distant pnums
// distant edgenums mit globalv1, globalv2
int distfa = (*recvarray)[ii++];
int globfa = (*recvarray)[ii++];
int distnp = (*recvarray)[ii++];
int distned =(*recvarray)[ii++];
int locfa = (*glob2locfa) [globfa];
if ( locfa == -1 )
{
ii += 2*distnp + 3*distned;
locfa = 1;
continue;
}
Array<int> edges;
int fa = locfa;
Array<int> pnums, globalpnums;
topology.GetFaceEdges ( fa, edges );
topology.GetFaceVertices ( fa, pnums );
globalpnums.SetSize ( distnp );
for ( int i = 0; i < distnp; i++)
globalpnums[i] = GetLoc2Glob_Vert ( pnums[i] );
SetDistantFaceNum ( sender, fa, distfa );
// find exchange points
for ( int i = 0; i < distnp; i++)
{
int distglobalpnum = (*recvarray)[ii+i];
for ( int j = 0; j < distnp; j++ )
if ( globalpnums[j] == distglobalpnum )
{
// set sender -- distpnum ---- locpnum
int distpnum = (*recvarray)[ii + i +distnp];
SetDistantPNum ( sender, pnums[j], distpnum );
}
}
int * distedgenums = new int [distned];
// find exchange edges
for ( int i = 0; i < edges.Size(); i++)
{
int v1, v2;
topology . GetEdgeVertices ( edges[i], v1, v2 );
int dv1 = GetLoc2Glob_Vert ( v1 );
int dv2 = GetLoc2Glob_Vert ( v2 );
if ( dv1 > dv2 ) swap ( dv1, dv2 );
for ( int ed = 0; ed < distned; ed++)
{
distedgenums[ed] = (*recvarray)[ii + 2*distnp + 3*ed];
int ddv1 = (*recvarray)[ii + 2*distnp + 3*ed + 1];
int ddv2 = (*recvarray)[ii + 2*distnp + 3*ed + 2];
if ( ddv1 > ddv2 ) swap ( ddv1, ddv2 );
if ( dv1 == ddv1 && dv2 == ddv2 )
{
// set sender -- distednum -- locednum
SetDistantEdgeNum ( sender, edges[i], distedgenums[ed] );
}
}
}
delete [] distedgenums;
ii += 2*distnp + 3*distned;
}
}
}
// set which elements are where for the master processor
delete sendarray; delete recvarray;
if ( id > 0 )
delete glob2locfa;
coarseupdate = 1;
#ifdef SCALASCA
#pragma pomp inst end(updatecoarsegrid)
#endif
}
void ParallelMeshTopology :: UpdateCoarseGridOverlap ()
{
UpdateCoarseGridGlobal();
#ifdef SCALASCA
#pragma pomp inst begin(updatecoarsegrid)
#endif
(*testout) << "UPDATE COARSE GRID PARALLEL TOPOLOGY, OVERLAP " << endl;
PrintMessage ( 1, "UPDATE COARSE GRID PARALLEL TOPOLOGY, OVERLAP " );
const MeshTopology & topology = mesh.GetTopology();
nfa = topology . GetNFaces();
ned = topology . GetNEdges();
np = mesh . GetNP();
nv = mesh . GetNV();
ne = mesh . GetNE();
nseg = mesh.GetNSeg();
nsurfel = mesh.GetNSE();
if ( id != 0 )
{
// find exchange edges - first send exchangeedges locnum, v1, v2
// receive distant distnum, v1, v2
// find matching
Array<int> * sendarray, *recvarray;
sendarray = new Array<int> (0);
recvarray = new Array<int>;
sendarray->SetSize (0);
BitArray recvface(nfa);
recvface.Clear();
for ( int el = 1; el <= ne; el++ )
{
Array<int> edges, pnums, faces;
topology.GetElementFaces (el, faces);
int globeli = GetLoc2Glob_VolEl(el);
for ( int fai = 0; fai < faces.Size(); fai++)
{
int fa = faces[fai];
topology.GetFaceEdges ( fa, edges );
topology.GetFaceVertices ( fa, pnums );
if ( !IsExchangeElement ( el ) ) continue;
int globfa = GetDistantFaceNum(0, fa) ;
// send :
// localfacenum
// globalfacenum
// globalelnum
// np
// ned
// globalpnums
// localpnums
// localedgenums mit globalelnums mit globalv1, globalv2
//
sendarray -> Append ( fa );
sendarray -> Append ( globfa );
sendarray -> Append ( globeli );
sendarray -> Append ( pnums.Size() );
sendarray -> Append ( edges.Size() );
for ( int i = 0; i < pnums.Size(); i++ )
{
sendarray -> Append( GetLoc2Glob_Vert(pnums[i]) );
}
for ( int i = 0; i < pnums.Size(); i++ )
{
sendarray -> Append(pnums[i] );
}
for ( int i = 0; i < edges.Size(); i++ )
{
int globedge = GetDistantEdgeNum(0, edges[i] );
int v1, v2;
topology . GetEdgeVertices ( edges[i], v1, v2 );
int dv1 = GetLoc2Glob_Vert ( v1 );
int dv2 = GetLoc2Glob_Vert ( v2 );
sendarray -> Append(edges[i] );
sendarray -> Append (globedge);
sendarray -> Append ( dv1 );
sendarray -> Append ( dv2 );
}
}
}
BitArray edgeisinit(ned), vertisinit(np);
edgeisinit.Clear();
vertisinit.Clear();
// Array for temporary use, to find local from global element fast
// only for not too big meshes
// seems ok, as low-order space is treated on one proc
Array<int,1> * glob2loc_el;
glob2loc_el = new Array<int,1> ( neglob );
(*glob2loc_el) = -1;
for ( int locel = 1; locel <= mesh.GetNE(); locel++)
(*glob2loc_el)[GetLoc2Glob_VolEl(locel)] = locel;
for ( int sender = 1; sender < ntasks; sender ++ )
{
if ( id == sender )
MyMPI_Bcast (*sendarray, sender-1, MPI_HIGHORDER_COMM);
// {
// for ( int dest = 1; dest < ntasks; dest ++ )
// if ( dest != id)
// {
// MyMPI_Send (*sendarray, dest);
// }
// }
if ( id != sender )
{
// MyMPI_Recv ( *recvarray, sender);
MyMPI_Bcast (*recvarray, sender-1, MPI_HIGHORDER_COMM);
// compare received vertices with own ones
int ii = 0;
int cntel = 0;
int volel = 1;
while ( ii< recvarray -> Size() )
{
// receive list :
// distant facenum
// np
// ned
// globalpnums
// distant pnums
// distant edgenums mit globalv1, globalv2
int distfa = (*recvarray)[ii++];
int globfa = (*recvarray)[ii++];
int globvolel = (*recvarray)[ii++];
int distnp = (*recvarray)[ii++];
int distned =(*recvarray)[ii++];
if ( id > 0 ) // GetLoc2Glob_VolEl ( volel ) != globvolel )
volel = (*glob2loc_el)[globvolel]; //Glob2Loc_VolEl ( globvolel );
else
volel = globvolel;
if ( volel == -1 )
{
ii += 2*distnp + 4*distned;
volel = 1;
continue;
}
Array<int> faces, edges;
topology.GetElementFaces( volel, faces);
topology.GetElementEdges ( volel, edges);
for ( int fai= 0; fai < faces.Size(); fai++ )
{
int fa = faces[fai];
if ( !IsExchangeFace ( fa ) && sender != 0 ) continue;
// if ( recvface.Test ( fa-1 ) ) continue;
Array<int> pnums, globalpnums;
//topology.GetFaceEdges ( fa, edges );
topology.GetFaceVertices ( fa, pnums );
// find exchange faces ...
// have to be of same type
if ( pnums.Size () != distnp ) continue;
globalpnums.SetSize ( distnp );
for ( int i = 0; i < distnp; i++)
globalpnums[i] = GetLoc2Glob_Vert ( pnums[i] );
// test if 3 vertices match
bool match = 1;
for ( int i = 0; i < distnp; i++)
if ( !globalpnums.Contains ( (*recvarray)[ii+i] ) )
match = 0;
if ( !match ) continue;
// recvface.Set(fa-1);
SetDistantFaceNum ( sender, fa, distfa );
SetDistantFaceNum ( 0, fa, globfa );
// find exchange points
for ( int i = 0; i < distnp; i++)
{
int distglobalpnum = (*recvarray)[ii+i];
for ( int j = 0; j < distnp; j++ )
if ( globalpnums[j] == distglobalpnum )
{
// set sender -- distpnum ---- locpnum
int distpnum = (*recvarray)[ii + i +distnp];
SetDistantPNum ( sender, pnums[j], distpnum );
}
}
int * distedgenums = new int [distned];
// find exchange edges
for ( int i = 0; i < edges.Size(); i++)
{
int v1, v2;
topology . GetEdgeVertices ( edges[i], v1, v2 );
int dv1 = GetLoc2Glob_Vert ( v1 );
int dv2 = GetLoc2Glob_Vert ( v2 );
if ( dv1 > dv2 ) swap ( dv1, dv2 );
for ( int ed = 0; ed < distned; ed++)
{
distedgenums[ed] = (*recvarray)[ii + 2*distnp + 4*ed];
int globedgenum = (*recvarray)[ii + 2*distnp + 4*ed + 1];
int ddv1 = (*recvarray)[ii + 2*distnp + 4*ed + 2];
int ddv2 = (*recvarray)[ii + 2*distnp + 4*ed + 3];
if ( ddv1 > ddv2 ) swap ( ddv1, ddv2 );
if ( dv1 == ddv1 && dv2 == ddv2 )
{
// set sender -- distednum -- locednum
SetDistantEdgeNum ( sender, edges[i], distedgenums[ed] );
SetDistantEdgeNum ( 0, edges[i], globedgenum );
}
}
}
delete [] distedgenums;
}
ii += 2*distnp + 4*distned;
}
}
}
// set which elements are where for the master processor
delete sendarray; delete recvarray;
if ( id > 0 )
delete glob2loc_el;
coarseupdate = 1;
}
// send global-local el/face/edge/vert-info to id 0
// nfa = topology . GetNFaces();
// ned = topology . GetNEdges();
// np = mesh . GetNP();
// nv = mesh . GetNV();
// ne = mesh . GetNE();
// nseg = mesh.GetNSeg();
// nsurfel = mesh.GetNSE();
if ( id != 0 )
{
Array<int> * sendarray;
sendarray = new Array<int> (4);
int sendnfa = 0, sendned = 0;
(*sendarray)[0] = ne;
(*sendarray)[1] = nfa;
(*sendarray)[2] = ned;
(*sendarray)[3] = np;
int ii = 4;
for ( int el = 1; el <= ne; el++ )
(*sendarray).Append ( GetLoc2Glob_VolEl (el ) );
for ( int fa = 1; fa <= nfa; fa++ )
{
if ( !IsExchangeFace (fa) ) continue;
sendnfa++;
(*sendarray).Append ( fa );
(*sendarray).Append ( GetDistantFaceNum (0, fa) );
}
for ( int ed = 1; ed <= ned; ed++ )
{
if ( !IsExchangeEdge (ed) ) continue;
sendned++;
sendarray->Append ( ed );
sendarray->Append ( GetDistantEdgeNum(0, ed) );
}
for ( int vnum = 1; vnum <= np; vnum++ )
sendarray->Append ( GetLoc2Glob_Vert(vnum) );
(*sendarray)[1] = sendnfa;
(*sendarray)[2] = sendned;
MyMPI_Send (*sendarray, 0);
delete sendarray;
}
else
{
Array<int> * recvarray = new Array<int>;
for ( int sender = 1; sender < ntasks; sender++ )
{
MyMPI_Recv ( *recvarray, sender);
int distnel = (*recvarray)[0];
int distnfa = (*recvarray)[1];
int distned = (*recvarray)[2];
int distnp = (*recvarray)[3];
int ii = 4;
for ( int el = 1; el <= distnel; el++ )
SetDistantEl ( sender, (*recvarray)[ii++], el );
for ( int fa = 1; fa <= distnfa; fa++ )
{
int distfa = (*recvarray)[ii++];
SetDistantFaceNum ( sender, (*recvarray)[ii++], distfa );
}
for ( int ed = 1; ed <= distned; ed++ )
{
int disted = (*recvarray)[ii++];
SetDistantEdgeNum ( sender, (*recvarray)[ii++], disted );
}
for ( int vnum = 1; vnum <= distnp; vnum++ )
SetDistantPNum ( sender, (*recvarray)[ii++], vnum );
}
delete recvarray;
}
#ifdef SCALASCA
#pragma pomp inst end(updatecoarsegrid)
#endif
}
void ParallelMeshTopology :: UpdateTopology ()
{
// loop over parallel faces and edges, find new local face/edge number,
const MeshTopology & topology = mesh.GetTopology();
int nfa = topology.GetNFaces();
int ned = topology.GetNEdges();
isghostedge.SetSize(ned);
isghostface.SetSize(nfa);
isghostedge.Clear();
isghostface.Clear();
for ( int ed = 1; ed <= ned; ed++)
{
int v1, v2;
topology.GetEdgeVertices ( ed, v1, v2 );
if ( IsGhostVert(v1) || IsGhostVert(v2) )
SetGhostEdge ( ed );
}
Array<int> pnums;
for ( int fa = 1; fa <= nfa; fa++)
{
topology.GetFaceVertices ( fa, pnums );
for ( int i = 0; i < pnums.Size(); i++)
if ( IsGhostVert( pnums[i] ) )
{
SetGhostFace ( fa );
break;
}
}
}
void ParallelMeshTopology :: UpdateExchangeElements()
{
(*testout) << "UPDATE EXCHANGE ELEMENTS " << endl;
const MeshTopology & topology = mesh.GetTopology();
isexchangeedge->SetSize ( (ntasks+1) * topology.GetNEdges() );
isexchangeface->SetSize ( (ntasks+1) * topology.GetNFaces() );
isexchangeedge->Clear();
isexchangeface->Clear();
for ( int eli = 1; eli <= mesh.GetNE(); eli++)
{
if ( ! IsExchangeElement ( eli ) ) continue;
const Element & el = mesh.VolumeElement(eli);
Array<int> faces, edges;
int np = el.NP();
topology.GetElementEdges ( eli, edges );
topology.GetElementFaces ( eli, faces );
for ( int i = 0; i < edges.Size(); i++)
{
SetExchangeEdge ( edges[i] );
}
for ( int i = 0; i < faces.Size(); i++)
{
SetExchangeFace ( faces[i] );
}
for ( int i = 0; i < np; i++)
{
SetExchangeVert ( el[i] );
}
}
if ( id == 0 ) return;
Array<int> ** elementonproc, ** recvelonproc;
elementonproc = new Array<int>*[ntasks];
recvelonproc = new Array<int>*[ntasks];
for ( int i = 1; i < ntasks; i++ )
{
elementonproc[i] = new Array<int>(0);
recvelonproc[i] = new Array<int>(0);
}
for ( int eli = 1; eli <= mesh.GetNE(); eli++ )
{
if ( !IsExchangeElement(eli) ) continue;
for ( int i = 1; i < ntasks; i++ )
if ( GetDistantElNum(i, eli) != -1 && i != id )
{
elementonproc[i] -> Append(eli);
elementonproc[i] -> Append(GetDistantElNum(i, eli));
}
}
for ( int sender = 1; sender < ntasks; sender ++ )
{
if ( id == sender )
for ( int dest = 1; dest < ntasks; dest ++ )
if ( dest != id)
{
MyMPI_Send ( *(elementonproc[dest]), dest);
elementonproc[dest] -> SetSize(0);
}
if ( id != sender )
{
MyMPI_Recv (*( recvelonproc[sender]), sender);
}
}
int ii = 0;
for ( int sender = 1; sender < ntasks; sender++ )
{
if ( sender == id ) continue;
ii = 0;
while ( recvelonproc[sender]->Size() > ii )
{
int distelnum = (*recvelonproc[sender])[ii++];
int locelnum = (*recvelonproc[sender])[ii++];
SetDistantEl ( sender, locelnum, distelnum);
}
recvelonproc[sender]->SetSize(0);
}
BitArray procs(ntasks);
procs.Clear();
for ( int eli = 1; eli <= mesh.GetNE(); eli++)
{
if ( IsGhostEl(eli) ) continue;
if ( !IsExchangeElement(eli) ) continue;
procs.Clear();
int sumprocs = 0;
for ( int i = 1; i < ntasks; i++ )
if ( GetDistantElNum(i, eli) != -1 && i != id )
{
procs.Set(i);
sumprocs++;
}
for ( int dest = 1; dest < ntasks; dest++)
{
if ( !procs.Test(dest) ) continue;
elementonproc[dest]->Append(GetDistantElNum(dest, eli));
elementonproc[dest]->Append(sumprocs);
for ( int i = 1; i < ntasks; i++ )
if ( procs.Test(i) )
{
elementonproc[dest]->Append(i);
elementonproc[dest]->Append(GetDistantElNum(i, eli));
}
}
}
for ( int sender = 1; sender < ntasks; sender ++ )
{
if ( id == sender )
for ( int dest = 1; dest < ntasks; dest ++ )
if ( dest != id)
{
MyMPI_Send ( *(elementonproc[dest]), dest);
delete elementonproc[dest];
}
if ( id != sender )
{
MyMPI_Recv (*( recvelonproc[sender]), sender);
}
}
for ( int sender = 1; sender < ntasks; sender++ )
{
if ( sender == id ) continue;
ii = 0;
while ( recvelonproc[sender]->Size() > ii )
{
int locelnum = (*recvelonproc[sender])[ii++];
int nprocs = (*recvelonproc[sender])[ii++];
for ( int iproc = 0; iproc < nprocs; iproc++)
{
int proc = (*recvelonproc[sender])[ii++];
int distelnum = (*recvelonproc[sender])[ii++];
if ( id == proc ) continue;
SetExchangeElement (locelnum, proc);
SetDistantEl( proc, locelnum, distelnum );
}
}
delete recvelonproc[sender];
}
delete [] elementonproc;
delete [] recvelonproc;
}
void ParallelMeshTopology :: SetNV ( const int anv )
{
*testout << "called setnv" << endl
<< "old size: " << loc2distvert.Size() << endl
<< "new size: " << anv << endl;
loc2distvert.ChangeSize (anv);
for (int i = 1; i <= anv; i++)
if (loc2distvert.EntrySize(i) == 0)
loc2distvert.Add (i, -1); // will be the global nr
BitArray * isexchangevert2 = new BitArray( (ntasks+1) * anv );
isexchangevert2->Clear();
if ( isexchangevert )
{
for ( int i = 0; i < min2( isexchangevert->Size(), isexchangevert2->Size() ); i++ )
if ( isexchangevert->Test(i) ) isexchangevert2->Set(i);
delete isexchangevert;
}
isexchangevert = isexchangevert2;
nv = anv;
}
void ParallelMeshTopology :: SetNE ( const int ane )
{
loc2distel.ChangeSize (ane);
for (int i = 0; i < ane; i++)
{
if (loc2distel[i].Size() == 0)
loc2distel.Add (i, -1); // will be the global nr
}
BitArray * isexchangeel2 = new BitArray ( (ntasks+1) * ane );
isexchangeel2->Clear();
if ( isexchangeel )
{
for ( int i = 0; i < min2(isexchangeel->Size(), isexchangeel2->Size() ) ; i++ )
if ( isexchangeel->Test(i) ) isexchangeel2->Set(i);
delete isexchangeel;
}
ne = ane;
isexchangeel = isexchangeel2;
}
void ParallelMeshTopology :: SetNSE ( int anse )
{
loc2distsurfel.ChangeSize (anse);
for (int i = 0; i < anse; i++)
if (loc2distsurfel[i].Size() == 0)
loc2distsurfel.Add (i, -1); // will be the global nr
nsurfel = anse;
}
void ParallelMeshTopology :: SetNSegm ( int anseg )
{
loc2distsegm.ChangeSize (anseg);
for (int i = 0; i < anseg; i++)
if (loc2distsegm[i].Size() == 0)
loc2distsegm.Add (i, -1); // will be the global nr
nseg = anseg;
}
}
#endif