ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-27 21:30:35 +05:00
Code Issues Packages Projects Releases Wiki Activity

load 2d meshes in parallel

Browse Source
This commit is contained in:
Joachim Schoeberl 2011-07-16 23:13:26 +00:00
parent 496d33ff17
commit 40eae2fbae
6 changed files with 326 additions and 367 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

243
libsrc/meshing/parallelmesh.cpp
View File

@ -31,10 +31,9 @@ namespace netgen
if (id == 0) if (id == 0)
PrintMessage (1, "Send/Receive mesh"); PrintMessage (1, "Send/Receive mesh");
{ {
// distribute global information // distribute global information
int nelglob, nvglob; int nelglob, nseglob, nvglob;
if (id == 0) if (id == 0)
{ {
paralleltop -> SetNV (GetNV()); paralleltop -> SetNV (GetNV());
@ -43,15 +42,19 @@ namespace netgen
paralleltop -> SetNSE (GetNSE()); paralleltop -> SetNSE (GetNSE());
nelglob = GetNE(); nelglob = GetNE();
nseglob = GetNSE();
nvglob = GetNV(); nvglob = GetNV();
} }
MyMPI_Bcast (dimension);
MyMPI_Bcast (nelglob); MyMPI_Bcast (nelglob);
MyMPI_Bcast (nseglob);
MyMPI_Bcast (nvglob); MyMPI_Bcast (nvglob);
if (id > 0) if (id > 0)
{ {
paralleltop -> SetNEGlob (nelglob); paralleltop -> SetNEGlob (nelglob);
paralleltop -> SetNSEGlob (nseglob);
paralleltop -> SetNVGlob (nvglob); paralleltop -> SetNVGlob (nvglob);
} }
} }
@ -98,25 +101,31 @@ namespace netgen
{ {
Array<MPI_Request> sendrequests; Array<MPI_Request> sendrequests;
PrintMessage ( 3, "Sending vertices");
Array<int> num_els_on_proc(ntasks); Array<int> num_els_on_proc(ntasks);
num_els_on_proc = 0; num_els_on_proc = 0;
for (ElementIndex ei = 0; ei < GetNE(); ei++) for (ElementIndex ei = 0; ei < GetNE(); ei++)
num_els_on_proc[(*this)[ei].GetPartition()]++; num_els_on_proc[(*this)[ei].GetPartition()]++;
// get number of vertices for each processor
Array<int> nelloc (ntasks);
nelloc = 0;
// elarraysize = 1;
PrintMessage ( 3, "Sending vertices");
TABLE<ElementIndex> els_of_proc (num_els_on_proc); TABLE<ElementIndex> els_of_proc (num_els_on_proc);
for (ElementIndex ei = 0; ei < GetNE(); ei++) for (ElementIndex ei = 0; ei < GetNE(); ei++)
els_of_proc.Add ( (*this)[ei].GetPartition(), ei); els_of_proc.Add ( (*this)[ei].GetPartition(), ei);
Array<int> num_sels_on_proc(ntasks);
num_sels_on_proc = 0;
for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
num_sels_on_proc[(*this)[ei].GetPartition()]++;
TABLE<SurfaceElementIndex> sels_of_proc (num_sels_on_proc);
for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
sels_of_proc.Add ( (*this)[ei].GetPartition(), ei);
Array<int, PointIndex::BASE> vert_flag (GetNV()); Array<int, PointIndex::BASE> vert_flag (GetNV());
Array<int, PointIndex::BASE> num_procs_on_vert (GetNV()); Array<int, PointIndex::BASE> num_procs_on_vert (GetNV());
Array<int> num_verts_on_proc (ntasks); Array<int> num_verts_on_proc (ntasks);
@ -125,6 +134,12 @@ namespace netgen
num_procs_on_vert = 0; num_procs_on_vert = 0;
vert_flag = -1; vert_flag = -1;
Array<int> nelloc (ntasks);
nelloc = 0;
Array<int> nselloc (ntasks);
nselloc = 0;
for (int dest = 1; dest < ntasks; dest++) for (int dest = 1; dest < ntasks; dest++)
{ {
FlatArray<ElementIndex> els = els_of_proc[dest]; FlatArray<ElementIndex> els = els_of_proc[dest];
@ -145,11 +160,32 @@ namespace netgen
paralleltop -> SetDistantPNum ( dest, epi, num_verts_on_proc[dest]); paralleltop -> SetDistantPNum ( dest, epi, num_verts_on_proc[dest]);
} }
} }
// elarraysize[dest] += 3 + el.GetNP();
nelloc[dest] ++; nelloc[dest] ++;
paralleltop -> SetDistantEl ( dest, els[hi]+1, nelloc[dest] ); paralleltop -> SetDistantEl ( dest, els[hi]+1, nelloc[dest] );
} }
FlatArray<SurfaceElementIndex> sels = sels_of_proc[dest];
for (int hi = 0; hi < sels.Size(); hi++)
{
const Element2d & el = (*this) [ sels[hi] ];
for (int i = 0; i < el.GetNP(); i++)
{
PointIndex epi = el[i];
if (vert_flag[epi] < dest)
{
vert_flag[epi] = dest;
num_verts_on_proc[dest]++;
num_procs_on_vert[epi]++;
paralleltop -> SetDistantPNum ( dest, epi, num_verts_on_proc[dest]);
}
}
nselloc[dest] ++;
paralleltop -> SetDistantSurfEl ( dest, sels[hi]+1, nselloc[dest] );
}
} }
TABLE<PointIndex> verts_of_proc (num_verts_on_proc); TABLE<PointIndex> verts_of_proc (num_verts_on_proc);
@ -163,7 +199,6 @@ namespace netgen
for (int hi = 0; hi < els.Size(); hi++) for (int hi = 0; hi < els.Size(); hi++)
{ {
const Element & el = (*this) [ els[hi] ]; const Element & el = (*this) [ els[hi] ];
for (int i = 0; i < el.GetNP(); i++) for (int i = 0; i < el.GetNP(); i++)
{ {
PointIndex epi = el[i]; PointIndex epi = el[i];
@ -174,6 +209,22 @@ namespace netgen
} }
} }
} }
FlatArray<SurfaceElementIndex> sels = sels_of_proc[dest];
for (int hi = 0; hi < sels.Size(); hi++)
{
const Element2d & el = (*this) [ sels[hi] ];
for (int i = 0; i < el.GetNP(); i++)
{
PointIndex epi = el[i];
if (vert_flag[epi] < dest)
{
vert_flag[epi] = dest;
procs_of_vert.Add (epi, dest);
}
}
}
} }
for (int vert = 1; vert <= GetNP(); vert++ ) for (int vert = 1; vert <= GetNP(); vert++ )
@ -241,6 +292,8 @@ namespace netgen
for ( int dest = 1; dest < ntasks; dest ++ ) for ( int dest = 1; dest < ntasks; dest ++ )
sendrequests.Append (MyMPI_ISend (distpnums[dest], dest, MPI_TAG_MESH+1)); sendrequests.Append (MyMPI_ISend (distpnums[dest], dest, MPI_TAG_MESH+1));
PrintMessage ( 3, "Sending elements" ); PrintMessage ( 3, "Sending elements" );
Array<int> elarraysize (ntasks); Array<int> elarraysize (ntasks);
@ -271,7 +324,6 @@ namespace netgen
PrintMessage ( 3, "Sending Face Descriptors" ); PrintMessage ( 3, "Sending Face Descriptors" );
Array<double> fddata (6 * GetNFD()); Array<double> fddata (6 * GetNFD());
for (int fdi = 1; fdi <= GetNFD(); fdi++) for (int fdi = 1; fdi <= GetNFD(); fdi++)
{ {
@ -288,59 +340,31 @@ namespace netgen
PrintMessage ( 3, "Sending Surface elements" ); PrintMessage ( 3, "Sending Surface elements" );
Array <int> nlocsel(ntasks), bufsize(ntasks); // , seli(ntasks); Array <int> nlocsel(ntasks), bufsize(ntasks);
for ( int i = 0; i < ntasks; i++) nlocsel = 0;
{ bufsize = 1;
nlocsel[i] = 0;
bufsize[i] = 1;
}
for (int sei = 1; sei <= GetNSE(); sei++ ) for (int sei = 1; sei <= GetNSE(); sei++ )
{ {
int ei1, ei2;
GetTopology().GetSurface2VolumeElement (sei, ei1, ei2);
const Element2d & sel = SurfaceElement (sei); const Element2d & sel = SurfaceElement (sei);
int dest = sel.GetPartition();
for (int j = 0; j < 2; j++)
{
int ei = (j == 0) ? ei1 : ei2;
if ( ei > 0 && ei <= GetNE() )
{
const Element & el = VolumeElement (ei);
int dest = el.GetPartition();
nlocsel[dest] ++; nlocsel[dest] ++;
bufsize[dest] += 4 + 2*sel.GetNP(); bufsize[dest] += 4 + 2*sel.GetNP();
} }
}
}
TABLE<int> selbuf(bufsize); TABLE<int> selbuf(bufsize);
Array<int> nselloc (ntasks);
nselloc = 0;
for (int dest = 1; dest < ntasks; dest++ ) for (int dest = 1; dest < ntasks; dest++ )
selbuf.Add (dest, nlocsel[dest]); selbuf.Add (dest, nlocsel[dest]);
for (int sei = 1; sei <= GetNSE(); sei ++ ) for (int sei = 1; sei <= GetNSE(); sei ++ )
{ {
int ei1, ei2;
GetTopology().GetSurface2VolumeElement (sei, ei1, ei2);
const Element2d & sel = SurfaceElement (sei); const Element2d & sel = SurfaceElement (sei);
int dest = sel.GetPartition();
int isghost = 0;
for (int j = 0; j < 2; j++)
{
int ei = (j == 0) ? ei1 : ei2;
if ( ei > 0 && ei <= GetNE() )
{
const Element & el = VolumeElement (ei);
int dest = el.GetPartition();
if (dest > 0)
{
selbuf.Add (dest, sei); selbuf.Add (dest, sei);
selbuf.Add (dest, sel.GetIndex()); selbuf.Add (dest, sel.GetIndex());
selbuf.Add (dest, isghost); selbuf.Add (dest, 0);
selbuf.Add (dest, sel.GetNP()); selbuf.Add (dest, sel.GetNP());
for ( int ii = 1; ii <= sel.GetNP(); ii++) for ( int ii = 1; ii <= sel.GetNP(); ii++)
@ -348,12 +372,6 @@ namespace netgen
selbuf.Add (dest, sel.PNum(ii)); selbuf.Add (dest, sel.PNum(ii));
selbuf.Add (dest, sel.GeomInfoPi(ii).trignum); selbuf.Add (dest, sel.GeomInfoPi(ii).trignum);
} }
nselloc[dest] ++;
paralleltop -> SetDistantSurfEl ( dest, sei, nselloc[dest] );
isghost = 1;
}
}
}
} }
for (int dest = 1; dest < ntasks; dest++) for (int dest = 1; dest < ntasks; dest++)
@ -373,13 +391,13 @@ namespace netgen
{ {
Array<int> volels; Array<int> volels;
const MeshTopology & topol = GetTopology(); const MeshTopology & topol = GetTopology();
topol . GetSegmentVolumeElements ( segi, volels ); topol . GetSegmentSurfaceElements ( segi, volels );
for (int j = 0; j < volels.Size(); j++) for (int j = 0; j < volels.Size(); j++)
{ {
int ei = volels[j]; int ei = volels[j];
if ( ei > 0 && ei <= GetNE() ) if ( ei > 0 && ei <= GetNSE() )
{ {
const Element & el = VolumeElement (ei); const Element2d & el = SurfaceElement (ei);
int dest = el.GetPartition(); int dest = el.GetPartition();
nlocseg[dest] ++; nlocseg[dest] ++;
bufsize[dest] += 14; bufsize[dest] += 14;
@ -392,15 +410,15 @@ namespace netgen
for ( int ls=1; ls <= GetNSeg(); ls++) for ( int ls=1; ls <= GetNSeg(); ls++)
{ {
Array<int> volels; Array<int> volels;
GetTopology().GetSegmentVolumeElements ( ls, volels ); GetTopology().GetSegmentSurfaceElements ( ls, volels );
const Segment & seg = LineSegment (ls); const Segment & seg = LineSegment (ls);
for (int j = 0; j < volels.Size(); j++) for (int j = 0; j < volels.Size(); j++)
{ {
int ei = volels[j]; int ei = volels[j];
if ( ei > 0 && ei <= GetNE() ) if ( ei > 0 && ei <= GetNSE() )
{ {
const Element & el = VolumeElement (ei); const Element2d & el = SurfaceElement (ei);
int dest = el.GetPartition(); int dest = el.GetPartition();
if ( dest > 0 ) if ( dest > 0 )
@ -430,6 +448,8 @@ namespace netgen
sendrequests.Append (MyMPI_ISend(segmbuf[dest], dest, MPI_TAG_MESH+5)); sendrequests.Append (MyMPI_ISend(segmbuf[dest], dest, MPI_TAG_MESH+5));
MPI_Waitall (sendrequests.Size(), &sendrequests[0], MPI_STATUS_IGNORE); MPI_Waitall (sendrequests.Size(), &sendrequests[0], MPI_STATUS_IGNORE);
MPI_Barrier(MPI_COMM_WORLD);
} }
@ -451,8 +471,6 @@ namespace netgen
// string st; // string st;
// receive vertices // receive vertices
NgProfiler::StartTimer (timer_pts); NgProfiler::StartTimer (timer_pts);
@ -463,7 +481,7 @@ namespace netgen
paralleltop -> SetNV (numvert); paralleltop -> SetNV (numvert);
// INDEX_CLOSED_HASHTABLE<int> glob2loc_vert_ht (3*numvert+1); // INDEX_CLOSED_HASHTABLE<int> glob2loc_vert_ht (3*numvert+1);
INDEX_HASHTABLE<int> glob2loc_vert_ht (100); INDEX_HASHTABLE<int> glob2loc_vert_ht (3*numvert+1);
for (int vert = 0; vert < numvert; vert++) for (int vert = 0; vert < numvert; vert++)
{ {
@ -489,7 +507,6 @@ namespace netgen
NgProfiler::StopTimer (timer_pts); NgProfiler::StopTimer (timer_pts);
*testout << "got " << numvert << " vertices" << endl; *testout << "got " << numvert << " vertices" << endl;
{ {
Element el; Element el;
@ -498,7 +515,6 @@ namespace netgen
NgProfiler::RegionTimer reg(timer_els); NgProfiler::RegionTimer reg(timer_els);
// double sumcosts = 0, sum = 0;
for (int ind = 0, elnum = 1; ind < elarray.Size(); elnum++) for (int ind = 0, elnum = 1; ind < elarray.Size(); elnum++)
{ {
paralleltop->SetLoc2Glob_VolEl ( elnum, elarray[ind++]); paralleltop->SetLoc2Glob_VolEl ( elnum, elarray[ind++]);
@ -506,27 +522,13 @@ namespace netgen
el.SetIndex(elarray[ind++]); el.SetIndex(elarray[ind++]);
el.SetNP(elarray[ind++]); el.SetNP(elarray[ind++]);
/*
for ( int j = 0; j < el.GetNP(); j++)
{
int costs = 1; // glob2loc_vert_ht.CalcPositionCosts (elarray[ind+j]);
// *testout << "costs to hash " << elarray[ind+j] << ", hashvalue = "
// << glob2loc_vert_ht.HashValue(elarray[ind+j]) << ": " << costs << endl;
sumcosts += costs;
sum++;
}
*/
for ( int j = 0; j < el.GetNP(); j++) for ( int j = 0; j < el.GetNP(); j++)
el[j] = glob2loc_vert_ht.Get (elarray[ind++]); el[j] = glob2loc_vert_ht.Get (elarray[ind++]);
AddVolumeElement (el); AddVolumeElement (el);
} }
// *testout << "average hashing costs: " << sumcosts / sum << endl;
// cout << "average hashing costs: " << sumcosts / sum << endl;
} }
{ {
Array<double> fddata; Array<double> fddata;
MyMPI_Recv (fddata, 0, MPI_TAG_MESH+3); MyMPI_Recv (fddata, 0, MPI_TAG_MESH+3);
@ -541,7 +543,6 @@ namespace netgen
{ {
NgProfiler::RegionTimer reg(timer_sels); NgProfiler::RegionTimer reg(timer_sels);
// surface elements
Array<int> selbuf; Array<int> selbuf;
MyMPI_Recv ( selbuf, 0, MPI_TAG_MESH+4); MyMPI_Recv ( selbuf, 0, MPI_TAG_MESH+4);
@ -620,12 +621,14 @@ namespace netgen
} }
} }
MPI_Barrier(MPI_COMM_WORLD);
int timerloc = NgProfiler::CreateTimer ("Update local mesh"); int timerloc = NgProfiler::CreateTimer ("Update local mesh");
int timerloc2 = NgProfiler::CreateTimer ("CalcSurfacesOfNode"); int timerloc2 = NgProfiler::CreateTimer ("CalcSurfacesOfNode");
NgProfiler::RegionTimer regloc(timerloc); NgProfiler::RegionTimer regloc(timerloc);
PrintMessage (2, "Got ", GetNE(), " elements"); PrintMessage (2, "Got ", GetNE(), " elements");
PrintMessage (2, "Got ", GetNSE(), " surface elements");
NgProfiler::StartTimer (timerloc2); NgProfiler::StartTimer (timerloc2);
@ -638,6 +641,7 @@ namespace netgen
SetNextMajorTimeStamp(); SetNextMajorTimeStamp();
// paralleltop->Print(); // paralleltop->Print();
cout << "receive mesh complete, id = " << id << endl;
} }
@ -648,14 +652,10 @@ namespace netgen
// call it only for the master ! // call it only for the master !
void Mesh :: Distribute () void Mesh :: Distribute ()
{ {
if ( id != 0 || ntasks == 1 ) return; if (id != 0 || ntasks == 1 ) return;
// metis partition of mesh, only if more than one proc
// partition mesh
#ifdef METIS #ifdef METIS
ParallelMetis (); ParallelMetis ();
#else #else
for (ElementIndex ei = 0; ei < GetNE(); ei++) for (ElementIndex ei = 0; ei < GetNE(); ei++)
@ -668,17 +668,11 @@ namespace netgen
for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++) for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
*testout << "sel(" << int(ei) << ") is in part " << (*this)[ei].GetPartition() << endl; *testout << "sel(" << int(ei) << ") is in part " << (*this)[ei].GetPartition() << endl;
// send partition // send partition
// for (int dest = 1; dest < ntasks; dest++)
// MyMPI_Send ("mesh", dest, MPI_TAG_CMD);
MyMPI_SendCmd ("mesh"); MyMPI_SendCmd ("mesh");
SendRecvMesh (); SendRecvMesh ();
paralleltop -> UpdateCoarseGrid(); paralleltop -> UpdateCoarseGrid();
// paralleltop -> Print(); // paralleltop -> Print();
} }
@ -704,49 +698,41 @@ namespace netgen
int ne = GetNE(); int ne = GetNE();
int nn = GetNP(); int nn = GetNP();
if (ntasks <= 2) if (ntasks <= 2 || ne <= 1)
{ {
if (ntasks == 1) return; if (ntasks == 1) return;
for (int i=1; i<=ne; i++) for (int i=1; i<=ne; i++)
VolumeElement(i).SetPartition(1); VolumeElement(i).SetPartition(1);
for (int i=1; i<=GetNSE(); i++)
SurfaceElement(i).SetPartition(1);
return; return;
} }
bool uniform_els = true; bool uniform_els = true;
ELEMENT_TYPE elementtype = TET; // VolumeElement(1).GetType();
// metis works only for uniform tet/hex meshes ELEMENT_TYPE elementtype = TET;
for ( int el = 2; el <= GetNE(); el++ ) for ( int el = 1; el <= GetNE(); el++ )
if ( VolumeElement(el).GetType() != elementtype ) if ( VolumeElement(el).GetType() != elementtype )
{ {
// int nelperproc = ne / (ntasks-1);
// for (int i=1; i<=ne; i++)
// {
// int partition = i / nelperproc + 1;
// if ( partition >= ntasks ) partition = ntasks-1;
// VolumeElement(i).SetPartition(partition);
// }
uniform_els = false; uniform_els = false;
break; break;
} }
// uniform_els = false;
if (!uniform_els) if (!uniform_els)
{ {
PartHybridMesh ( ); // neloc ); PartHybridMesh ( );
return; return;
} }
// uniform (TET) mesh, JS // uniform (TET) mesh, JS
int npe = VolumeElement(1).GetNP(); int npe = VolumeElement(1).GetNP();
Array<idxtype> elmnts(ne*npe);
idxtype *elmnts;
elmnts = new idxtype[ne*npe];
int etype; int etype;
if (elementtype == TET) if (elementtype == TET)
@ -763,10 +749,7 @@ namespace netgen
int nparts = ntasks-1; int nparts = ntasks-1;
int edgecut; int edgecut;
idxtype *epart, *npart; Array<idxtype> epart(ne), npart(nn);
epart = new idxtype[ne];
npart = new idxtype[nn];
// if ( ntasks == 1 ) // if ( ntasks == 1 )
// { // {
@ -787,35 +770,49 @@ namespace netgen
// } // }
cout << "call metis ... " << flush; cout << "call metis ... " << flush;
int timermetis = NgProfiler::CreateTimer ("Metis itself"); int timermetis = NgProfiler::CreateTimer ("Metis itself");
NgProfiler::StartTimer (timermetis); NgProfiler::StartTimer (timermetis);
METIS_PartMeshDual (&ne, &nn, elmnts, &etype, &numflag, &nparts, METIS_PartMeshDual (&ne, &nn, &elmnts[0], &etype, &numflag, &nparts,
&edgecut, epart, npart); &edgecut, &epart[0], &npart[0]);
NgProfiler::StopTimer (timermetis); NgProfiler::StopTimer (timermetis);
cout << "complete" << endl; cout << "complete" << endl;
cout << "edge-cut: " << edgecut << ", balance: " << ComputeElementBalance(ne, nparts, epart) << endl; cout << "edge-cut: " << edgecut << ", balance: "
<< ComputeElementBalance(ne, nparts, &epart[0]) << endl;
// partition numbering by metis : 0 ... ntasks - 1 // partition numbering by metis : 0 ... ntasks - 1
// we want: 1 ... ntasks // we want: 1 ... ntasks
for (int i=1; i<=ne; i++) for (int i=1; i<=ne; i++)
VolumeElement(i).SetPartition(epart[i-1] + 1); VolumeElement(i).SetPartition(epart[i-1] + 1);
delete [] elmnts; for (int sei = 1; sei <= GetNSE(); sei++ )
delete [] epart; {
delete [] npart; int ei1, ei2;
GetTopology().GetSurface2VolumeElement (sei, ei1, ei2);
Element2d & sel = SurfaceElement (sei);
for (int j = 0; j < 2; j++)
{
int ei = (j == 0) ? ei1 : ei2;
if ( ei > 0 && ei <= GetNE() )
{
sel.SetPartition (VolumeElement(ei).GetPartition());
break;
}
}
}
} }
#endif #endif
void Mesh :: PartHybridMesh () // Array<int> & neloc ) void Mesh :: PartHybridMesh ()
{ {
#ifdef METIS #ifdef METIS
int ne = GetNE(); int ne = GetNE();

239
libsrc/meshing/paralleltop.cpp
View File

@ -429,16 +429,12 @@ namespace netgen
if (id == 0) if (id == 0)
PrintMessage ( 3, "UPDATE GLOBAL COARSEGRID STARTS" ); // JS PrintMessage ( 3, "UPDATE GLOBAL COARSEGRID STARTS" ); // JS
MPI_Group MPI_GROUP_WORLD; MPI_Group MPI_GROUP_WORLD;
int process_ranks[] = { 0 };
int n_ho = netgen::ntasks - 1; MPI_Comm_group (MPI_COMM_WORLD, &MPI_GROUP_WORLD);
int * process_ranks = new int[netgen::ntasks-1]; MPI_Group_excl (MPI_GROUP_WORLD, 1, process_ranks, &MPI_HIGHORDER_WORLD);
for ( int i = 0; i < netgen::ntasks-1; i++ ) MPI_Comm_create (MPI_COMM_WORLD, MPI_HIGHORDER_WORLD, &MPI_HIGHORDER_COMM);
process_ranks[i] = i+1;
MPI_Comm_group ( MPI_COMM_WORLD, &MPI_GROUP_WORLD);
MPI_Group_incl ( MPI_GROUP_WORLD, n_ho, process_ranks, & MPI_HIGHORDER_WORLD);
MPI_Comm_create ( MPI_COMM_WORLD, MPI_HIGHORDER_WORLD, & MPI_HIGHORDER_COMM);
int timer = NgProfiler::CreateTimer ("UpdateCoarseGridGlobal"); int timer = NgProfiler::CreateTimer ("UpdateCoarseGridGlobal");
@ -446,6 +442,8 @@ namespace netgen
*testout << "ParallelMeshTopology :: UpdateCoarseGridGlobal" << endl; *testout << "ParallelMeshTopology :: UpdateCoarseGridGlobal" << endl;
cout << "update global, id = " << id << endl;
const MeshTopology & topology = mesh.GetTopology(); const MeshTopology & topology = mesh.GetTopology();
Array<int> sendarray, recvarray; Array<int> sendarray, recvarray;
@ -481,67 +479,12 @@ namespace netgen
} }
if ( id == 0 ) if ( id == 0 )
{
sendarray.Append (nfa); sendarray.Append (nfa);
sendarray.Append (ned);
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; Array<int> edges, pnums, faces, elpnums;
sendarray.Append (ne); sendarray.Append (ne);
for ( int el = 1; el <= ne; el++ ) for ( int el = 1; el <= ne; el++ )
@ -551,6 +494,7 @@ namespace netgen
const Element & volel = mesh.VolumeElement (el); const Element & volel = mesh.VolumeElement (el);
int globeli = GetLoc2Glob_VolEl(el); int globeli = GetLoc2Glob_VolEl(el);
// cout << "el = " << el << ", globeli = " << globeli << endl;
sendarray. Append ( globeli ); sendarray. Append ( globeli );
sendarray. Append ( faces.Size() ); sendarray. Append ( faces.Size() );
@ -568,63 +512,47 @@ namespace netgen
else else
sendarray. Append(GetLoc2Glob_Vert (volel[i])); sendarray. Append(GetLoc2Glob_Vert (volel[i]));
} }
// end new version
sendarray.Append (nsurfel);
BitArray edgeisinit(ned), vertisinit(np); for (int el = 1; el <= nsurfel; el++)
edgeisinit.Clear();
vertisinit.Clear();
// Array for temporary use, to find local from global element fast
Array<int,1> glob2loc_el;
if ( id != 0 )
{ {
topology.GetSurfaceElementEdges ( el, edges );
const Element2d & surfel = mesh.SurfaceElement (el);
sendarray. Append ( el );
sendarray. Append ( edges.Size() );
sendarray. Append ( surfel.GetNP() );
for ( int i = 0; i < edges.Size(); i++ )
sendarray. Append(edges[i] );
for ( int i = 0; i < surfel.GetNP(); i++ )
sendarray. Append(surfel[i] );
}
}
if (id == 0)
MyMPI_Bcast ( sendarray );
else
MyMPI_Bcast ( recvarray );
cout << "have done bcast, id = " << id << endl;
if (id != 0)
{
Array<int,1> glob2loc_el;
glob2loc_el.SetSize (neglob); glob2loc_el.SetSize (neglob);
glob2loc_el = -1; glob2loc_el = -1;
for ( int locel = 1; locel <= mesh.GetNE(); locel++) for ( int locel = 1; locel <= mesh.GetNE(); locel++)
glob2loc_el[GetLoc2Glob_VolEl(locel)] = 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
sendrequest = MyMPI_ISend ( sendarray, 0, MPI_TAG_MESH );
int nloops = (id == 0) ? ntasks-1 : 1;
for (int hi = 0; hi < nloops; hi++)
{
int sender;
if (id == 0)
{
sender = MyMPI_Recv ( recvarray, MPI_TAG_MESH );
// PrintMessage (4, "have received from ", sender);
}
else
{
MyMPI_Bcast ( recvarray );
sender = 0;
}
// compare received vertices with own ones
int ii = 0; int ii = 0;
// int cntel = 0;
int volel = 1;
if ( id != 0 )
nfaglob = recvarray[ii++]; nfaglob = recvarray[ii++];
nedglob = recvarray[ii++];
if (id == 0) continue;
Array<int> faces, edges; Array<int> faces, edges;
Array<int> pnums, globalpnums; Array<int> pnums, globalpnums;
@ -637,11 +565,7 @@ namespace netgen
int distned = recvarray[ii++]; int distned = recvarray[ii++];
int distnp = recvarray[ii++]; int distnp = recvarray[ii++];
if ( id > 0 ) int volel = glob2loc_el[globvolel];
volel = glob2loc_el[globvolel];
else
volel = globvolel;
if (volel != -1) if (volel != -1)
{ {
topology.GetElementFaces( volel, faces); topology.GetElementFaces( volel, faces);
@ -649,24 +573,67 @@ namespace netgen
const Element & volelement = mesh.VolumeElement (volel); const Element & volelement = mesh.VolumeElement (volel);
for ( int i = 0; i < faces.Size(); i++) for ( int i = 0; i < faces.Size(); i++)
SetDistantFaceNum ( sender, faces[i], recvarray[ii++]); SetDistantFaceNum ( 0, faces[i], recvarray[ii++]);
for ( int i = 0; i < edges.Size(); i++) for ( int i = 0; i < edges.Size(); i++)
SetDistantEdgeNum ( sender, edges[i], recvarray[ii++]); SetDistantEdgeNum ( 0, edges[i], recvarray[ii++]);
for ( int i = 0; i < distnp; i++) for ( int i = 0; i < volelement.GetNP(); i++)
SetDistantPNum ( sender, volelement[i], recvarray[ii++]); SetDistantPNum ( 0, volelement[i], recvarray[ii++]);
} }
else else
ii += distnfa + distned + distnp; ii += distnfa + distned + distnp;
} }
Array<int,1> glob2loc_sel;
int recv_nse = recvarray[ii++];
nseglob = recv_nse;
glob2loc_sel.SetSize (nseglob);
glob2loc_sel = -1;
for ( int locel = 1; locel <= mesh.GetNSE(); locel++)
glob2loc_sel[GetLoc2Glob_SurfEl(locel)] = locel;
for (int hi = 0; hi < recv_nse; hi++)
{
int globvolel = recvarray[ii++];
int distned = recvarray[ii++];
int distnp = recvarray[ii++];
int surfel = glob2loc_sel[globvolel];
if (surfel != -1)
{
topology.GetSurfaceElementEdges ( surfel, edges);
const Element2d & element = mesh.SurfaceElement (surfel);
for ( int i = 0; i < edges.Size(); i++)
SetDistantEdgeNum ( 0, edges[i], recvarray[ii++]);
for ( int i = 0; i < element.GetNP(); i++)
SetDistantPNum ( 0, element[i], recvarray[ii++]);
}
else
ii += distned + distnp;
} }
coarseupdate = 1; }
if (id != 0) if (id != 0)
MPI_Wait (&sendrequest, MPI_STATUS_IGNORE); {
*testout << "l2d - vert = " << loc2distvert << endl;
*testout << "l2d - edge = " << loc2distedge << endl;
*testout << "l2d - el = " << loc2distel << endl;
*testout << "l2d - sel = " << loc2distsurfel << endl;
}
MPI_Barrier (MPI_COMM_WORLD);
cout << "update global complete, id = " << id << endl;
MPI_Barrier (MPI_COMM_WORLD);
coarseupdate = 1;
} }
@ -676,6 +643,7 @@ namespace netgen
void ParallelMeshTopology :: UpdateCoarseGrid () void ParallelMeshTopology :: UpdateCoarseGrid ()
{ {
cout << "updatecoarsegrid, id = " << id << endl;
static int timer = NgProfiler::CreateTimer ("UpdateCoarseGrid"); static int timer = NgProfiler::CreateTimer ("UpdateCoarseGrid");
NgProfiler::RegionTimer reg(timer); NgProfiler::RegionTimer reg(timer);
@ -684,6 +652,7 @@ namespace netgen
if (id == 0) if (id == 0)
PrintMessage (1, "UPDATE COARSE GRID PARALLEL TOPOLOGY "); PrintMessage (1, "UPDATE COARSE GRID PARALLEL TOPOLOGY ");
// find exchange edges - first send exchangeedges locnum, v1, v2 // find exchange edges - first send exchangeedges locnum, v1, v2
// receive distant distnum, v1, v2 // receive distant distnum, v1, v2
// find matching // find matching
@ -719,8 +688,10 @@ namespace netgen
Array<int> cnt_send(ntasks-1); Array<int> cnt_send(ntasks-1);
MPI_Barrier (MPI_HIGHORDER_COMM);
cout << "before ex edges" << endl;
*testout << "before ex edges" << endl;
MPI_Barrier (MPI_HIGHORDER_COMM);
NgProfiler::StartTimer (timere); NgProfiler::StartTimer (timere);
@ -793,6 +764,15 @@ namespace netgen
NgProfiler::StopTimer (timere); NgProfiler::StopTimer (timere);
MPI_Barrier (MPI_HIGHORDER_COMM);
cout << "before ex faces" << endl;
*testout << "before ex faces" << endl;
MPI_Barrier (MPI_HIGHORDER_COMM);
if (mesh.GetDimension() == 3)
{
NgProfiler::StartTimer (timerf); NgProfiler::StartTimer (timerf);
glob2loc.SetSize (nfaglob); glob2loc.SetSize (nfaglob);
@ -855,7 +835,12 @@ namespace netgen
NgProfiler::StopTimer (timerf); NgProfiler::StopTimer (timerf);
}
MPI_Barrier (MPI_HIGHORDER_COMM);
cout << "after ex faces" << endl;
*testout << "after ex faces" << endl;
MPI_Barrier (MPI_HIGHORDER_COMM);
// set which elements are where for the master processor // set which elements are where for the master processor

5
libsrc/meshing/paralleltop.hpp
View File

@ -17,9 +17,10 @@ namespace netgen
int nseg, nsurfel; int nseg, nsurfel;
// number of global elements, vertices, ???, faces // number of global elements, vertices, ???, faces
int neglob, nvglob; int neglob, nseglob, nvglob;
int nparel; int nparel;
int nfaglob; int nfaglob;
int nedglob;
/** /**
mapping from local to distant vertex number mapping from local to distant vertex number
@ -45,6 +46,7 @@ namespace netgen
void SetNVGlob ( int anvglob ) { nvglob = anvglob; } void SetNVGlob ( int anvglob ) { nvglob = anvglob; }
void SetNEGlob ( int aneglob ) { neglob = aneglob; } void SetNEGlob ( int aneglob ) { neglob = aneglob; }
void SetNSEGlob ( int anseglob ) { nseglob = anseglob; }
int GetNVGlob () { return nvglob; } int GetNVGlob () { return nvglob; }
int GetNEGlob () { return neglob; } int GetNEGlob () { return neglob; }
@ -59,6 +61,7 @@ namespace netgen
int GetLoc2Glob_Vert ( int locnum ) const { return loc2distvert[locnum][0]; } int GetLoc2Glob_Vert ( int locnum ) const { return loc2distvert[locnum][0]; }
int GetLoc2Glob_VolEl ( int locnum ) const { return loc2distel[locnum-1][0]; } int GetLoc2Glob_VolEl ( int locnum ) const { return loc2distel[locnum-1][0]; }
int GetLoc2Glob_SurfEl ( int locnum ) const { return loc2distsurfel[locnum-1][0]; }
void GetVertNeighbours ( int vnum, Array<int> & dests ) const; void GetVertNeighbours ( int vnum, Array<int> & dests ) const;

97
libsrc/meshing/topology.cpp
View File

@ -1114,50 +1114,10 @@ namespace netgen
if ( ntasks > 1 ) if ( ntasks > 1 )
{ {
if ( !paralleltop.DoCoarseUpdate() ) continue; if ( !paralleltop.DoCoarseUpdate() ) continue;
// "illegal" faces are exchange faces
/*
(*testout) << "exchange face : " << i << endl;
(*testout) << "points = " << face2vert[i] << endl;
*/
// (*testout) << "global points = ";
// for ( int j = 0; j < 3; j++ )
// // (*testout) << face2vert[i].I(j+1) << " -> "
// // << paralleltop.GetLoc2Glob_Vert( face2vert[i].I(j+1) ) << ", ";
// (*testout) << endl;
// if ( !paralleltop.IsExchangeFace (i+1) )
// paralleltop.SetRefinementFace (i+1);
/*
paralleltop.SetExchangeFace (i+1);
for (int j = 0; j < 4; j++)
{
if ( face2vert[i].I(j+1) > 0 )
paralleltop.SetExchangeVert(face2vert[i].I(j+1));
}
Array<int> faceedges;
GetFaceEdges (i+1, faceedges);
for ( int j = 0; j < faceedges.Size(); j++)
{
paralleltop.SetExchangeEdge ( faceedges[j] );
int v1, v2;
GetEdgeVertices(faceedges[j], v1, v2 );
}
*/
/*
(*testout) << "pos = ";
for (int j = 0; j < 4; j++)
if (face2vert[i].I(j+1) >= 1)
(*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " ";
(*testout) << endl;
*/
} }
else else
{
#endif #endif
{
(*testout) << "illegal face : " << i << endl; (*testout) << "illegal face : " << i << endl;
(*testout) << "points = " << face2vert[i] << endl; (*testout) << "points = " << face2vert[i] << endl;
(*testout) << "pos = "; (*testout) << "pos = ";
@ -1185,31 +1145,30 @@ namespace netgen
} }
} }
#ifdef PARALLEL
} }
#endif
} }
} }
#ifndef PARALLEL
if (cnt_err)
cout << cnt_err << " elements are not matching !!!" << endl;
#else
if (cnt_err && ntasks == 1) if (cnt_err && ntasks == 1)
cout << cnt_err << " elements are not matching !!!" << endl; cout << cnt_err << " elements are not matching !!!" << endl;
else if (cnt_err && ntasks > 1)
{ if (cnt_err && ntasks > 1)
// cout << "p" << id << ": " << cnt_err << " elements are not local" << endl;
isparallel = 1; isparallel = 1;
} }
else if ( ntasks > 1 )
;
// cout << "p" << id << ": " << "Partition " << id << " is totally local" << endl;
#endif
} }
#ifdef PARALLEL #ifdef PARALLEL
if ( isparallel ) cout << "is par = " << int(isparallel) << ", id = " << id << endl;
if (mesh.GetDimension() == 3)
if (isparallel != (id != 0))
{
cerr << " ****************************** " << endl;
cerr << " **** wrong isparallel ****** " << endl;
cerr << " ****************************** " << endl;
}
if (id != 0) // if (isparallel)
{ {
paralleltop.Update(); paralleltop.Update();
if ( paralleltop.DoCoarseUpdate() ) if ( paralleltop.DoCoarseUpdate() )
@ -1219,11 +1178,6 @@ namespace netgen
}
/* /*
for (i = 1; i <= ne; i++) for (i = 1; i <= ne; i++)
{ {
@ -1737,7 +1691,8 @@ namespace netgen
void MeshTopology :: GetSegmentVolumeElements ( int segnr, Array<int> & volels ) const void MeshTopology ::
GetSegmentVolumeElements ( int segnr, Array<int> & volels ) const
{ {
int v1, v2; int v1, v2;
GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 ); GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 );
@ -1749,6 +1704,24 @@ namespace netgen
for ( int eli1=1; eli1 <= volels1.Size(); eli1++) for ( int eli1=1; eli1 <= volels1.Size(); eli1++)
if ( volels2.Contains( volels1.Elem(eli1) ) ) if ( volels2.Contains( volels1.Elem(eli1) ) )
volels.Append ( volels1.Elem(eli1) ); volels.Append ( volels1.Elem(eli1) );
} }
void MeshTopology ::
GetSegmentSurfaceElements (int segnr, Array<int> & els) const
{
int v1, v2;
GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 );
Array<int> els1, els2;
GetVertexSurfaceElements ( v1, els1 );
GetVertexSurfaceElements ( v2, els2 );
els.SetSize(0);
for ( int eli1=1; eli1 <= els1.Size(); eli1++)
if ( els2.Contains( els1.Elem(eli1) ) )
els.Append ( els1.Elem(eli1) );
}
} }

3
libsrc/meshing/topology.hpp
View File

@ -124,7 +124,8 @@ public:
int GetVerticesEdge ( int v1, int v2) const; int GetVerticesEdge ( int v1, int v2) const;
void GetSegmentVolumeElements ( int segnr, Array<int> & surfels ) const; void GetSegmentVolumeElements ( int segnr, Array<int> & els ) const;
void GetSegmentSurfaceElements ( int segnr, Array<int> & els ) const;
}; };

2
libsrc/visualization/vssolution.cpp
View File

@ -2192,7 +2192,7 @@ namespace netgen
if (minv == maxv) maxv = minv+1e-6; if (minv == maxv) maxv = minv+1e-6;
#ifdef PARALLEL #ifdef PARALLEL
if (id == 0) if ((ntasks > 1) && (id == 0))
{ {
minv = 1e99; minv = 1e99;
maxv = -1e99; maxv = -1e99;