From 4f6f77f2f3e13f25c7f016b5135e95e522529164 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 25 Jan 2009 19:37:14 +0000 Subject: [PATCH] move parallel files --- libsrc/Makefile.am | 3 +- libsrc/Makefile.in | 4 +- libsrc/meshing/Makefile.am | 3 +- libsrc/meshing/Makefile.in | 17 +- libsrc/meshing/meshing.hpp | 2 +- libsrc/meshing/parallelmesh.cpp | 1910 +++++++++++++++++++++++++++++++ libsrc/meshing/paralleltop.cpp | 1701 +++++++++++++++++++++++++++ libsrc/meshing/paralleltop.hpp | 269 +++++ 8 files changed, 3900 insertions(+), 9 deletions(-) create mode 100644 libsrc/meshing/parallelmesh.cpp create mode 100644 libsrc/meshing/paralleltop.cpp create mode 100644 libsrc/meshing/paralleltop.hpp diff --git a/libsrc/Makefile.am b/libsrc/Makefile.am index 8395087d..2b7b9866 100644 --- a/libsrc/Makefile.am +++ b/libsrc/Makefile.am @@ -2,5 +2,4 @@ AM_CPPFLAGS = METASOURCES = AUTO -SUBDIRS = csg general geom2d gprim include interface linalg meshing occ \ - parallel stlgeom visualization +SUBDIRS = csg general geom2d gprim include interface linalg meshing occ stlgeom visualization diff --git a/libsrc/Makefile.in b/libsrc/Makefile.in index acdb70df..67a585b9 100644 --- a/libsrc/Makefile.in +++ b/libsrc/Makefile.in @@ -216,9 +216,7 @@ top_builddir = @top_builddir@ top_srcdir = @top_srcdir@ AM_CPPFLAGS = METASOURCES = AUTO -SUBDIRS = csg general geom2d gprim include interface linalg meshing occ \ - parallel stlgeom visualization - +SUBDIRS = csg general geom2d gprim include interface linalg meshing occ stlgeom visualization all: all-recursive .SUFFIXES: diff --git a/libsrc/meshing/Makefile.am b/libsrc/meshing/Makefile.am index e9f4760b..0dec170e 100644 --- a/libsrc/meshing/Makefile.am +++ b/libsrc/meshing/Makefile.am @@ -23,4 +23,5 @@ libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp \ pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp \ ruler2.cpp ruler3.cpp secondorder.cpp smoothing2.5.cpp \ smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp \ - topology.cpp triarls.cpp validate.cpp zrefine.cpp \ No newline at end of file + topology.cpp triarls.cpp validate.cpp zrefine.cpp \ + parallelmesh.cpp paralleltop.cpp paralleltop.hpp diff --git a/libsrc/meshing/Makefile.in b/libsrc/meshing/Makefile.in index 54de88dc..dd6ce836 100644 --- a/libsrc/meshing/Makefile.in +++ b/libsrc/meshing/Makefile.in @@ -57,7 +57,8 @@ am_libmesh_la_OBJECTS = adfront2.lo adfront3.lo bisect.lo \ prism2rls.lo pyramid2rls.lo pyramidrls.lo quadrls.lo refine.lo \ ruler2.lo ruler3.lo secondorder.lo smoothing2.5.lo \ smoothing2.lo smoothing3.lo specials.lo tetrarls.lo \ - topology.lo triarls.lo validate.lo zrefine.lo + topology.lo triarls.lo validate.lo zrefine.lo parallelmesh.lo \ + paralleltop.lo libmesh_la_OBJECTS = $(am_libmesh_la_OBJECTS) DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir) depcomp = $(SHELL) $(top_srcdir)/depcomp @@ -71,6 +72,15 @@ CXXLD = $(CXX) CXXLINK = $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \ --mode=link $(CXXLD) $(AM_CXXFLAGS) $(CXXFLAGS) $(AM_LDFLAGS) \ $(LDFLAGS) -o $@ +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \ + --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) \ + $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +CCLD = $(CC) +LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \ + --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \ + $(LDFLAGS) -o $@ SOURCES = $(libmesh_la_SOURCES) DIST_SOURCES = $(libmesh_la_SOURCES) HEADERS = $(noinst_HEADERS) @@ -254,7 +264,8 @@ libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp \ pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp \ ruler2.cpp ruler3.cpp secondorder.cpp smoothing2.5.cpp \ smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp \ - topology.cpp triarls.cpp validate.cpp zrefine.cpp + topology.cpp triarls.cpp validate.cpp zrefine.cpp \ + parallelmesh.cpp paralleltop.cpp paralleltop.hpp all: all-am @@ -331,6 +342,8 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/msghandler.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/netrule2.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/netrule3.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/parallelmesh.Plo@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/paralleltop.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/parser2.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/parser3.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/prism2rls.Plo@am__quote@ diff --git a/libsrc/meshing/meshing.hpp b/libsrc/meshing/meshing.hpp index c695c693..23b2c69a 100644 --- a/libsrc/meshing/meshing.hpp +++ b/libsrc/meshing/meshing.hpp @@ -67,7 +67,7 @@ namespace netgen #include "validate.hpp" #ifdef PARALLEL -#include "../parallel/paralleltop.hpp" +#include "paralleltop.hpp" // #include "../parallel/parallelmesh.hpp" #endif diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp new file mode 100644 index 00000000..8591fc3b --- /dev/null +++ b/libsrc/meshing/parallelmesh.cpp @@ -0,0 +1,1910 @@ +#ifdef PARALLEL + +#include +#include "paralleltop.hpp" + + +#ifdef METIS +namespace metis { extern "C" { +#include +} } +#endif + + +using namespace metis; + +namespace netgen +{ + + + // slaves receive the mesh from the master + void Mesh :: ReceiveParallelMesh ( ) + { + int timer = NgProfiler::CreateTimer ("ReceiveParallelMesh"); + NgProfiler::RegionTimer reg(timer); + + int timer_pts = NgProfiler::CreateTimer ("Receive Points"); + int timer_els = NgProfiler::CreateTimer ("Receive Elements"); + int timer_sels = NgProfiler::CreateTimer ("Receive Surface Elements"); + int timer_edges = NgProfiler::CreateTimer ("Receive Edges"); + + +#ifdef SCALASCA +#pragma pomp inst begin(loadmesh) +#endif + + // PrintMessage (1, "LOAD PARALLEL MESH"); + // MPI_Barrier (MPI_COMM_WORLD); + + string st; + double doublebuf[100]; + int i, n; + int tag_dim = 10, tag_token = 100, tag_n=11, tag_pnum=12, tag_point=13; + int tag_index = 101, tag_facedescr = 102; + MPI_Status status; + + bool endmesh = false; + + int dim; + int nelglob, nelloc, nvglob, nedglob, nfaglob; + + // receive global values + MPI_Bcast( &nelglob, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &nvglob, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &nedglob, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &nfaglob, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &dimension, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MyMPI_Recv ( nelloc, 0 ); + + paralleltop -> SetNVGlob ( nvglob ); + paralleltop -> SetNEGlob ( nelglob ); + + INDEX_CLOSED_HASHTABLE glob2loc_vert_ht (1); + + PrintMessage (1, "Receive mesh"); + + int ve = 0; + while (!endmesh) + { + MyMPI_Recv ( st, 0 ); + + // receive vertices + if (st == "vertex") + { + NgProfiler::RegionTimer reg(timer_pts); + + Array pointarray; + MyMPI_Recv ( pointarray, 0 ); + + int numvert = pointarray.Size() / 5; + paralleltop -> SetNV (numvert); + + glob2loc_vert_ht.SetSize (3*numvert+1); + + for ( int vert=0; vertSetLoc2Glob_Vert ( vert+1, globvert ); + glob2loc_vert_ht.Set (globvert, vert+1); + + netgen::Point<3> p; + p(0) = pointarray[vert*5+1]; + p(1) = pointarray[vert*5+2]; + p(2) = pointarray[vert*5+3]; + AddPoint (p); + (*this)[PointIndex(vert+1)] .Singularity ( pointarray[vert*5+4] ); + } + + Array dist_pnums; + MyMPI_Recv ( dist_pnums, 0); + + for (int hi = 0; hi < dist_pnums.Size(); hi += 3) + { + paralleltop -> + SetDistantPNum ( dist_pnums[hi+1], dist_pnums[hi], dist_pnums[hi+2]); + } + } + + if ( strcmp (st.c_str(), "volumeelements" ) == 0 ) + { + NgProfiler::RegionTimer reg(timer_els); + + *testout << "receiving elements" << endl; + + Element el; + + Array elarray; + MyMPI_Recv ( elarray, 0); + + int ind = 0; + int elnum = 1; + int nelloc = elarray[ind++]; + + paralleltop -> SetNE (nelloc); + + while ( ind < elarray.Size() ) + { + paralleltop->SetLoc2Glob_VolEl ( elnum, elarray[ind++]); + + el.SetIndex(elarray[ind++]); + el.SetNP(elarray[ind++]); + + for ( int j = 0; j < el.GetNP(); j++) + el[j] = glob2loc_vert_ht.Get (elarray[ind++]); + + AddVolumeElement (el); + elnum++; + } + } + + if (strcmp (st.c_str(), "facedescriptor") == 0) + { + Array doublebuf; + MyMPI_Recv( doublebuf, 0 ); + int faceind = AddFaceDescriptor (FaceDescriptor(int(doublebuf[0]), int(doublebuf[1]), int(doublebuf[2]), 0)); + GetFaceDescriptor(faceind).SetBCProperty (int(doublebuf[3])); + GetFaceDescriptor(faceind).domin_singular = doublebuf[4]; + GetFaceDescriptor(faceind).domout_singular = doublebuf[5]; + } + + + + if (strcmp (st.c_str(), "surfaceelementsgi") == 0) + { + NgProfiler::RegionTimer reg(timer_sels); + int j; + int surfnr, bcp, domin, domout, nep, faceind = 0; + int globsel; + int * selbuf; + selbuf = 0; + int bufsize; + // receive: + // faceindex + // nep + // tri.pnum + // tri.geominfopi.trignum + int nlocsel; + MyMPI_Recv ( selbuf, bufsize, 0); + int ii= 0; + int sel = 0; + + nlocsel = selbuf[ii++]; + paralleltop -> SetNSE ( nlocsel ); + + while ( ii < bufsize-1 ) + { + globsel = selbuf[ii++]; + faceind = selbuf[ii++]; + bool isghost = selbuf[ii++]; + nep = selbuf[ii++]; + Element2d tri(nep); + tri.SetIndex(faceind); + for( j=1; j<=nep; j++) + { + tri.PNum(j) = glob2loc_vert_ht.Get (selbuf[ii++]); + + // tri.GeomInfoPi(j).trignum = paralleltop->Glob2Loc_SurfEl(selbuf[ii++]); + tri.GeomInfoPi(j).trignum = selbuf[ii++]; + // Frage JS->AS: Brauchst du die trignum irgendwo ? + // sie sonst nur bei STL - Geometrien benötigt + // die Umrechnung war ein bottleneck ! + } + + tri.SetGhost(isghost); + + paralleltop->SetLoc2Glob_SurfEl ( sel+1, globsel ); + AddSurfaceElement (tri); + sel ++; + } + + delete [] selbuf ; + } + + if (strcmp (st.c_str(), "edgesegmentsgi") == 0) + { + NgProfiler::RegionTimer reg(timer_edges); + int endtime, starttime; + starttime = clock(); + + double * segmbuf = 0; + int bufsize; + MyMPI_Recv ( segmbuf, bufsize, 0); + Segment seg; + int hi; + int globsegi; + int ii = 0; + int segi = 1; + int nsegloc = int ( bufsize / 14 ) ; + paralleltop -> SetNSegm ( nsegloc ); + while ( ii < bufsize ) + { + globsegi = int (segmbuf[ii++]); + seg.si = int (segmbuf[ii++]); + + seg.p1 = glob2loc_vert_ht.Get (int(segmbuf[ii++])); + seg.p2 = glob2loc_vert_ht.Get (int(segmbuf[ii++])); + seg.geominfo[0].trignum = int( segmbuf[ii++] ); + seg.geominfo[1].trignum = int ( segmbuf[ii++]); + seg.surfnr1 = int ( segmbuf[ii++]); + seg.surfnr2 = int ( segmbuf[ii++]); + seg.edgenr = int ( segmbuf[ii++]); + seg.epgeominfo[0].dist = segmbuf[ii++]; + seg.epgeominfo[1].edgenr = int (segmbuf[ii++]); + seg.epgeominfo[1].dist = segmbuf[ii++]; + + seg.singedge_left = segmbuf[ii++]; + seg.singedge_right = segmbuf[ii++]; + + seg.epgeominfo[0].edgenr = seg.epgeominfo[1].edgenr; + + seg.domin = seg.surfnr1; + seg.domout = seg.surfnr2; + if ( seg.p1 >0 && seg.p2 > 0 ) + { + paralleltop-> SetLoc2Glob_Segm ( segi, globsegi ); + + AddSegment (seg); + segi++; + } + + } + delete [] segmbuf; + endtime = clock(); + (*testout) << "Receiving Time fde = " << double(endtime - starttime)/CLOCKS_PER_SEC << endl; + } + + /* + for ( int eli = 1; eli < GetNE(); eli++ ) + { + Element & el = VolumeElement(eli); + } + */ + + + if (strcmp (st.c_str(), "endmesh") == 0) + { + endmesh = true; + } + } + + + // ohne diesem Zusammenwarten gibts Abstürze, und ich weiß nicht warum ??? + // MPI_Barrier (MPI_COMM_WORLD); + // PrintMessage (1, "Have recevied the mesh"); + // MPI_Barrier (MPI_COMM_WORLD); + + + + // paralleltop -> SetNV ( this -> GetNV() ); +// for ( int i = 0; i < GetNV(); i++ ) +// paralleltop -> SetLoc2Glob_Vert ( i+1, (*loc2globvert)[i] ); + + + int timerloc = NgProfiler::CreateTimer ("Update local mesh"); + int timerloc2 = NgProfiler::CreateTimer ("CalcSurfacesOfNode"); + + NgProfiler::RegionTimer regloc(timerloc); + PrintMessage (2, "Got ", GetNE(), " elements"); + + NgProfiler::StartTimer (timerloc2); + + CalcSurfacesOfNode (); + + NgProfiler::StopTimer (timerloc2); + + // BuildConnectedNodes (); + + topology -> Update(); + +// UpdateOverlap(); + clusters -> Update(); + + SetNextMajorTimeStamp(); + + // paralleltop->Print(); + +#ifdef SCALASCA +#pragma pomp inst end(loadmesh) +#endif + } + + + + + + + + // distribute the mesh to the slave processors + // call it only for the master ! + void Mesh :: Distribute () + { + if ( id != 0 || ntasks == 1 ) return; + // metis partition of mesh, only if more than one proc + +#ifdef SCALASCA +#pragma pomp inst begin(metis) +#endif + + + // partition mesh + ParallelMetis (); + +#ifdef SCALASCA +#pragma pomp inst end(metis) +#endif + + // send partition + SendMesh (); + + + paralleltop -> UpdateCoarseGrid(); + +#ifdef SCALASCA +#pragma pomp inst end(loadmesh_seq) +#endif + + // paralleltop -> Print(); + } + + + + + + + void Mesh :: ParallelMetis ( ) + { + int timer = NgProfiler::CreateTimer ("Mesh::Partition"); + NgProfiler::RegionTimer reg(timer); + + PrintMessage (3, "Metis called"); + + if (GetDimension() == 2) + { + PartDualHybridMesh2D ( ); // neloc ); + return; + } + + + int ne = GetNE(); + int nn = GetNP(); + + if (ntasks <= 2) + { + if (ntasks == 1) return; + + for (int i=1; i<=ne; i++) + VolumeElement(i).SetPartition(1); + + return; + } + + + bool uniform_els = true; + ELEMENT_TYPE elementtype = TET; // VolumeElement(1).GetType(); + // metis works only for uniform tet/hex meshes + for ( int el = 2; el <= GetNE(); el++ ) + 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; + break; + } + + // uniform_els = false; + if (!uniform_els) + { + PartHybridMesh ( ); // neloc ); + return; + } + + + + // uniform (TET) mesh, JS + int npe = VolumeElement(1).GetNP(); + + idxtype *elmnts; + elmnts = new idxtype[ne*npe]; + + int etype; + if (elementtype == TET) + etype = 2; + else if (elementtype == HEX) + etype = 3; + + + for (int i=1; i<=ne; i++) + for (int j=1; j<=npe; j++) + elmnts[(i-1)*npe+(j-1)] = VolumeElement(i).PNum(j)-1; + + int numflag = 0; + int nparts = ntasks-1; + + int edgecut; + idxtype *epart, *npart; + epart = new idxtype[ne]; + npart = new idxtype[nn]; + + +// if ( ntasks == 1 ) +// { +// (*this) = *mastermesh; +// nparts = 4; +// metis :: METIS_PartMeshDual (&ne, &nn, elmnts, &etype, &numflag, &nparts, +// &edgecut, epart, npart); +// cout << "done" << endl; + +// cout << "edge-cut: " << edgecut << ", balance: " << metis :: ComputeElementBalance(ne, nparts, epart) << endl; + +// for (int i=1; i<=ne; i++) +// { +// mastermesh->VolumeElement(i).SetPartition(epart[i-1]); +// } + +// return; +// } + + + cout << "call metis ... " << flush; + + int timermetis = NgProfiler::CreateTimer ("Metis itself"); + NgProfiler::StartTimer (timermetis); + + metis :: METIS_PartMeshDual (&ne, &nn, elmnts, &etype, &numflag, &nparts, + &edgecut, epart, npart); + + NgProfiler::StopTimer (timermetis); + + cout << "complete" << endl; + cout << "edge-cut: " << edgecut << ", balance: " << metis :: ComputeElementBalance(ne, nparts, epart) << endl; + + + // partition numbering by metis : 0 ... ntasks - 1 + // we want: 1 ... ntasks + + for (int i=1; i<=ne; i++) + VolumeElement(i).SetPartition(epart[i-1] + 1); + + + delete [] elmnts; + delete [] epart; + delete [] npart; + } + + + + void Mesh :: PartHybridMesh () // Array & neloc ) + { + + int ne = GetNE(); + + int nn = GetNP(); + int nedges = topology->GetNEdges(); + + idxtype *xadj, * adjacency, *v_weights = NULL, *e_weights = NULL; + + int weightflag = 0; + int numflag = 0; + int nparts = ntasks - 1; + + int options[5]; + options[0] = 0; + int edgecut; + idxtype * part; + + xadj = new idxtype[nn+1]; + part = new idxtype[nn]; + + Array cnt(nn+1); + cnt = 0; + + for ( int edge = 1; edge <= nedges; edge++ ) + { + int v1, v2; + topology->GetEdgeVertices ( edge, v1, v2); + cnt[v1-1] ++; + cnt[v2-1] ++; + } + + xadj[0] = 0; + for ( int n = 1; n <= nn; n++ ) + { + xadj[n] = idxtype(xadj[n-1] + cnt[n-1]); + } + + adjacency = new idxtype[xadj[nn]]; + cnt = 0; + + for ( int edge = 1; edge <= nedges; edge++ ) + { + int v1, v2; + topology->GetEdgeVertices ( edge, v1, v2); + adjacency[ xadj[v1-1] + cnt[v1-1] ] = v2-1; + adjacency[ xadj[v2-1] + cnt[v2-1] ] = v1-1; + cnt[v1-1]++; + cnt[v2-1]++; + } + + for ( int vert = 0; vert < nn; vert++ ) + { + FlatArray array ( cnt[vert], &adjacency[ xadj[vert] ] ); + BubbleSort(array); + } + + metis :: METIS_PartGraphKway ( &nn, xadj, adjacency, v_weights, e_weights, &weightflag, + &numflag, &nparts, options, &edgecut, part ); + + Array nodesinpart(ntasks); + + for ( int el = 1; el <= ne; el++ ) + { + Element & volel = VolumeElement(el); + nodesinpart = 0; + + // VolumeElement(el).SetPartition(part[ volel[1] ] + 1); + + int el_np = volel.GetNP(); + int partition = 0; + for ( int i = 0; i < el_np; i++ ) + nodesinpart[ part[volel[i]-1]+1 ] ++; + + for ( int i = 1; i < ntasks; i++ ) + if ( nodesinpart[i] > nodesinpart[partition] ) + partition = i; + + volel.SetPartition(partition); + + } + + /* + for ( int i=1; i<=ne; i++) + { + neloc[ VolumeElement(i).GetPartition() ] ++; + } + */ + + delete [] xadj; + delete [] part; + delete [] adjacency; + } + + + void Mesh :: PartDualHybridMesh ( ) // Array & neloc ) + { + int ne = GetNE(); + + int nn = GetNP(); + int nedges = topology->GetNEdges(); + int nfaces = topology->GetNFaces(); + + idxtype *xadj, * adjacency, *v_weights = NULL, *e_weights = NULL; + + int weightflag = 0; + int numflag = 0; + int nparts = ntasks - 1; + + int options[5]; + options[0] = 0; + int edgecut; + idxtype * part; + + Array facevolels1(nfaces), facevolels2(nfaces); + facevolels1 = -1; + facevolels2 = -1; + + Array elfaces; + xadj = new idxtype[ne+1]; + part = new idxtype[ne]; + + Array cnt(ne+1); + cnt = 0; + + for ( int el=1; el <= ne; el++ ) + { + Element volel = VolumeElement(el); + topology->GetElementFaces(el, elfaces); + for ( int i = 0; i < elfaces.Size(); i++ ) + { + if ( facevolels1[elfaces[i]-1] == -1 ) + facevolels1[elfaces[i]-1] = el; + else + { + facevolels2[elfaces[i]-1] = el; + cnt[facevolels1[elfaces[i]-1]-1]++; + cnt[facevolels2[elfaces[i]-1]-1]++; + } + } + } + + xadj[0] = 0; + for ( int n = 1; n <= ne; n++ ) + { + xadj[n] = idxtype(xadj[n-1] + cnt[n-1]); + } + + adjacency = new idxtype[xadj[ne]]; + cnt = 0; + + for ( int face = 1; face <= nfaces; face++ ) + { + int e1, e2; + e1 = facevolels1[face-1]; + e2 = facevolels2[face-1]; + if ( e2 == -1 ) continue; + adjacency[ xadj[e1-1] + cnt[e1-1] ] = e2-1; + adjacency[ xadj[e2-1] + cnt[e2-1] ] = e1-1; + cnt[e1-1]++; + cnt[e2-1]++; + } + + for ( int el = 0; el < ne; el++ ) + { + FlatArray array ( cnt[el], &adjacency[ xadj[el] ] ); + BubbleSort(array); + } + + int timermetis = NgProfiler::CreateTimer ("Metis itself"); + NgProfiler::StartTimer (timermetis); + + metis :: METIS_PartGraphKway ( &ne, xadj, adjacency, v_weights, e_weights, &weightflag, + &numflag, &nparts, options, &edgecut, part ); + + NgProfiler::StopTimer (timermetis); + + Array nodesinpart(ntasks); + + for ( int el = 1; el <= ne; el++ ) + { + Element & volel = VolumeElement(el); + nodesinpart = 0; + + VolumeElement(el).SetPartition(part[el-1 ] + 1); + + } + + /* + for ( int i=1; i<=ne; i++) + { + neloc[ VolumeElement(i).GetPartition() ] ++; + } + */ + + delete [] xadj; + delete [] part; + delete [] adjacency; + } + + + + + + void Mesh :: PartDualHybridMesh2D ( ) + { + int ne = GetNSE(); + int nv = GetNV(); + + Array xadj(ne+1); + Array adjacency(ne*4); + + // first, build the vertex 2 element table: + Array cnt(nv); + cnt = 0; + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + for (int j = 0; j < (*this)[sei].GetNP(); j++) + cnt[ (*this)[sei][j] ] ++; + + TABLE vert2els(cnt); + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + for (int j = 0; j < (*this)[sei].GetNP(); j++) + vert2els.Add ((*this)[sei][j], sei); + + + // find all neighbour elements + int cntnb = 0; + Array marks(ne); // to visit each neighbour just once + marks = -1; + for (SurfaceElementIndex sei = 0; sei < ne; sei++) + { + xadj[sei] = cntnb; + for (int j = 0; j < (*this)[sei].GetNP(); j++) + { + PointIndex vnr = (*this)[sei][j]; + + // all elements with at least one common vertex + for (int k = 0; k < vert2els[vnr].Size(); k++) + { + SurfaceElementIndex sei2 = vert2els[vnr][k]; + if (sei == sei2) continue; + if (marks[sei2] == sei) continue; + + // neighbour, if two common vertices + int common = 0; + for (int m1 = 0; m1 < (*this)[sei].GetNP(); m1++) + for (int m2 = 0; m2 < (*this)[sei2].GetNP(); m2++) + if ( (*this)[sei][m1] == (*this)[sei2][m2]) + common++; + + if (common >= 2) + { + marks[sei2] = sei; // mark as visited + adjacency[cntnb++] = sei2; + } + } + } + } + xadj[ne] = cntnb; + + idxtype *v_weights = NULL, *e_weights = NULL; + + int weightflag = 0; + int numflag = 0; + int nparts = ntasks - 1; + + int options[5]; + options[0] = 0; + int edgecut; + Array part(ne); + + for ( int el = 0; el < ne; el++ ) + BubbleSort (adjacency.Range (xadj[el], xadj[el+1])); + + metis :: METIS_PartGraphKway ( &ne, &xadj[0], &adjacency[0], v_weights, e_weights, &weightflag, + &numflag, &nparts, options, &edgecut, &part[0] ); + + for (SurfaceElementIndex sei = 0; sei < ne; sei++) + (*this) [sei].SetPartition (part[sei]+1); + } + + + + + + + void Mesh :: SendMesh () const // Mesh * mastermesh, Array & neloc ) const + { + const Mesh * mastermesh = this; // the original plan was different + + PrintMessage ( 1, "Sending Mesh to local processors" ); + + int timer = NgProfiler::CreateTimer ("SendMesh"); + int timer2 = NgProfiler::CreateTimer ("SM::Prepare Points"); + int timer3 = NgProfiler::CreateTimer ("SM::Send Points"); + int timer4 = NgProfiler::CreateTimer ("SM::Send Elements"); + // int timer5 = NgProfiler::CreateTimer ("SM::Prepare Poins"); + + NgProfiler::RegionTimer reg(timer); + +#ifdef SCALASCA +#pragma pomp inst begin(sendmesh) +#endif + + NgProfiler::StartTimer (timer2); + + clock_t starttime, endtime, soltime; + starttime = clock(); + + + paralleltop -> SetNV ( GetNV() ); + paralleltop -> SetNE ( GetNE() ); + paralleltop -> SetNSegm ( GetNSeg() ); + paralleltop -> SetNSE ( GetNSE() ); + + + MPI_Request sendrequest[ntasks]; + + for ( int dest = 1; dest < ntasks; dest++) + MyMPI_Send ("mesh", dest); + + // MPI_Barrier (MPI_COMM_WORLD); + + + Array num_els_on_proc(ntasks); + num_els_on_proc = 0; + for (ElementIndex ei = 0; ei < mastermesh->GetNE(); ei++) + num_els_on_proc[(*this)[ei].GetPartition()]++; + + int nel = GetNE(); + int nv = GetNV(); + int nedges = (GetTopology().GetNEdges()); + int nfaces = GetTopology().GetNFaces(); + int dim = dimension; + MPI_Bcast( &nel, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &nv, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &nedges, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &nfaces, 1, MPI_INT, 0, MPI_COMM_WORLD ); + MPI_Bcast( &dim, 1, MPI_INT, 0, MPI_COMM_WORLD ); + for ( int dest = 1; dest < ntasks; dest++ ) + MyMPI_Send (num_els_on_proc[dest], dest); + + + // get number of vertices for each processor + Array elarraysize(ntasks); + Array nelloc ( ntasks ); + + nelloc = 0; + elarraysize = 1; + + PrintMessage ( 3, "Sending vertices"); + + TABLE els_of_proc (num_els_on_proc); + for (ElementIndex ei = 0; ei < mastermesh->GetNE(); ei++) + els_of_proc.Add ( (*this)[ei].GetPartition(), ei); + + + Array vert_flag ( GetNV() ); + + Array num_verts_on_proc (ntasks); + Array num_procs_on_vert ( GetNV() ); + + num_verts_on_proc = 0; + num_procs_on_vert = 0; + + vert_flag = -1; + for (int dest = 1; dest < ntasks; dest++) + { + FlatArray els = els_of_proc[dest]; + + for (int hi = 0; hi < els.Size(); hi++) + { + const Element & el = (*this) [ els[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]); + } + } + + elarraysize[dest] += 3 + el.GetNP(); + nelloc[dest] ++; + paralleltop -> SetDistantEl ( dest, els[hi]+1, nelloc[dest] ); + } + } + + TABLE verts_of_proc (num_verts_on_proc); + TABLE procs_of_vert (num_procs_on_vert); + TABLE loc_num_of_vert (num_procs_on_vert); + + vert_flag = -1; + for (int dest = 1; dest < ntasks; dest++) + { + FlatArray els = els_of_proc[dest]; + + for (int hi = 0; hi < els.Size(); hi++) + { + const Element & el = (*mastermesh) [ els[hi] ]; + + for (int i = 0; i < el.GetNP(); i++) + { + PointIndex epi = el.PNum(i+1); + if (vert_flag[epi] < dest) + { + vert_flag[epi] = dest; + procs_of_vert.Add (epi, dest); + } + } + } + } + + for (int vert = 1; vert <= mastermesh->GetNP(); vert++ ) + { + FlatArray procs = procs_of_vert[vert]; + for (int j = 0; j < procs.Size(); j++) + { + int dest = procs[j]; + verts_of_proc.Add (dest, vert); + loc_num_of_vert.Add (vert, verts_of_proc[dest].Size()); + } + } + + + Array nvi5(ntasks); + for (int i = 0; i < ntasks; i++) + nvi5[i] = 5 * num_verts_on_proc[i]; + + TABLE pointarrays(nvi5); + + NgProfiler::StopTimer (timer2); + NgProfiler::StartTimer (timer3); + + for (int dest = 1; dest < ntasks; dest++) + { + FlatArray verts = verts_of_proc[dest]; + + for ( int j = 0, ii = 0; j < verts.Size(); j++) + { + const MeshPoint & hp = mastermesh -> Point (verts[j]); + pointarrays.Add (dest, double(verts[j])); + pointarrays.Add (dest, hp(0)); + pointarrays.Add (dest, hp(1)); + pointarrays.Add (dest, hp(2)); + pointarrays.Add (dest, hp.Singularity()); + } + + MyMPI_Send ( "vertex", dest ); + MyMPI_ISend ( pointarrays[dest], dest, sendrequest[dest] ); + MPI_Request_free (&sendrequest[dest]); + } + + NgProfiler::StopTimer (timer3); + NgProfiler::StartTimer (timer4); + + Array num_distpnums(ntasks); + num_distpnums = 0; + + for (int vert = 1; vert <= mastermesh -> GetNP(); vert++) + { + FlatArray procs = procs_of_vert[vert]; + for (int j = 0; j < procs.Size(); j++) + num_distpnums[procs[j]] += 3 * (procs.Size()-1); + } + + TABLE distpnums (num_distpnums); + + for (int vert = 1; vert <= mastermesh -> GetNP(); vert++) + { + FlatArray procs = procs_of_vert[vert]; + for (int j = 0; j < procs.Size(); j++) + for (int k = 0; k < procs.Size(); k++) + if (j != k) + { + distpnums.Add (procs[j], loc_num_of_vert[vert][j]); + distpnums.Add (procs[j], procs_of_vert[vert][k]); + distpnums.Add (procs[j], loc_num_of_vert[vert][k]); + } + } + + for ( int dest = 1; dest < ntasks; dest ++ ) + { + MyMPI_ISend ( distpnums[dest], dest, sendrequest[dest] ); + MPI_Request_free (&sendrequest[dest]); + } + + + endtime = clock(); + (*testout) << "Sending Time verts = " << double(endtime - starttime)/CLOCKS_PER_SEC << endl; + + + + PrintMessage ( 3, "Sending elements" ); + TABLE elementarrays(elarraysize); + + starttime = clock(); + + // for ( int dest = 1; dest < ntasks; dest ++ ) + // MyMPI_Send ( "volumeelements", dest); + + for ( int dest = 1; dest < ntasks; dest++ ) + elementarrays.Add (dest, nelloc[dest]); + + for ( int ei = 1; ei <= mastermesh->GetNE(); ei++) + { + const Element & el = mastermesh -> VolumeElement (ei); + int dest = el.GetPartition(); + + if ( dest > 0 ) + { + // send volume element + + elementarrays.Add (dest, ei); // + elementarrays.Add (dest, el.GetIndex()); + elementarrays.Add (dest, el.GetNP()); + for ( int ii=0; ii double6(6); + for ( int dest = 1; dest < ntasks; dest++) + for ( int fdi = 1; fdi <= mastermesh->GetNFD(); fdi++) + { + MyMPI_Send("facedescriptor", dest); + + double6[0] = GetFaceDescriptor(fdi).SurfNr(); + double6[1] = GetFaceDescriptor(fdi).DomainIn(); + double6[2] = GetFaceDescriptor(fdi).DomainOut(); + double6[3] = GetFaceDescriptor(fdi).BCProperty(); + double6[4] = GetFaceDescriptor(fdi).domin_singular; + double6[5] = GetFaceDescriptor(fdi).domout_singular; + + MyMPI_Send ( double6, dest); + } + + endtime = clock(); + (*testout) << "Sending Time fdi = " << double(endtime - starttime)/CLOCKS_PER_SEC << endl; + starttime = clock(); + + // hasglobalsurf( edgenr * ntasks + dest ) .... global edge edgenr is in mesh at dest +// BitArray hasglobaledge; + +// hasglobaledge.SetSize(ntasks * nedglob); +// hasglobaledge.Clear(); +// cout << "mmf " << mastermesh->GetTopology().GetNFaces() << endl; + // determine sizes of local surface element arrays + + PrintMessage ( 3, "Sending Surface elements" ); + + Array nlocsel(ntasks), bufsize ( ntasks), seli(ntasks); + for ( int i = 0; i < ntasks; i++) + { + nlocsel[i] = 0; + bufsize[i] = 1; + seli[i] = 1; + } + + for ( int sei = 1; sei <= mastermesh -> GetNSE(); sei ++ ) + { + int ei1, ei2; + mastermesh -> GetTopology().GetSurface2VolumeElement (sei, ei1, ei2); + const Element2d & sel = mastermesh -> SurfaceElement (sei); + int dest; + // first element + + for (int j = 0; j < 2; j++) + { + int ei = (j == 0) ? ei1 : ei2; + + if ( ei > 0 && ei <= mastermesh->GetNE() ) + { + const Element & el = mastermesh -> VolumeElement (ei); + dest = el.GetPartition(); + nlocsel[dest] ++; + bufsize[dest] += 4 + 2*sel.GetNP(); + } + } + } + + int ** selbuf = 0; + selbuf = new int*[ntasks]; + for ( int i = 0; i < ntasks; i++) + if ( bufsize[i] > 0 ) + {*(selbuf+i) = new int[bufsize[i]];} + else + selbuf[i] = 0; + + + + Array nselloc (ntasks); + nselloc = 0; + + for ( int dest = 1; dest < ntasks; dest++ ) + { + MyMPI_Send ( "surfaceelementsgi", dest); + selbuf[dest][0] = nlocsel[dest]; + } + + for ( int sei = 1; sei <= mastermesh -> GetNSE(); sei ++ ) + { + int ei1, ei2; + mastermesh -> GetTopology().GetSurface2VolumeElement (sei, ei1, ei2); + const Element2d & sel = mastermesh -> SurfaceElement (sei); + int dest; + + + + int isghost = 0; + + for (int j = 0; j < 2; j++) + { + int ei = (j == 0) ? ei1 : ei2; + + if ( ei > 0 && ei <= mastermesh->GetNE() ) + { + const Element & el = mastermesh -> VolumeElement (ei); + dest = el.GetPartition(); + if (dest > 0) + { + // send: + // sei + // faceind + // nep + // tri.pnums, tri.geominfopi.trignums + + selbuf[dest][seli[dest]++] = sei; + selbuf[dest][seli[dest]++] = sel.GetIndex(); + selbuf[dest][seli[dest]++] = isghost; + selbuf[dest][seli[dest]++] = sel.GetNP(); + + for ( int ii = 1; ii <= sel.GetNP(); ii++) + { + selbuf[dest][seli[dest]++] = sel.PNum(ii); + selbuf[dest][seli[dest]++] = sel.GeomInfoPi(ii).trignum; + } + nselloc[dest] ++; + paralleltop -> SetDistantSurfEl ( dest, sei, nselloc[dest] ); + isghost = 1; + } + + } + } + } + + + for ( int dest = 1; dest < ntasks; dest++) + MyMPI_Send( selbuf[dest], bufsize[dest], dest); + + for ( int dest = 0; dest < ntasks; dest++ ) + { + if (selbuf[dest]) + delete [] *(selbuf+dest); + } + delete [] selbuf; + + + endtime = clock(); + (*testout) << "Sending Time surfels = " << double(endtime - starttime)/CLOCKS_PER_SEC << endl; + starttime = clock(); + + + PrintMessage ( 3, "Sending Edge Segments"); + for ( int dest = 1; dest < ntasks; dest++ ) + MyMPI_Send ( "edgesegmentsgi", dest); + + + Array nlocseg(ntasks), segi(ntasks); + for ( int i = 0; i < ntasks; i++) + { + nlocseg[i] = 0; + bufsize[i] = 0; + segi[i] = 0; + } + + for ( int segi = 1; segi <= mastermesh -> GetNSeg(); segi ++ ) + { + Array volels; + const MeshTopology & topol = mastermesh -> GetTopology(); + topol . GetSegmentVolumeElements ( segi, volels ); + const Segment & segm = mastermesh -> LineSegment (segi); + int dest; + for (int j = 0; j < volels.Size(); j++) + { + int ei = volels[j]; + int dest; + if ( ei > 0 && ei <= mastermesh->GetNE() ) + { + const Element & el = mastermesh -> VolumeElement (ei); + dest = el.GetPartition(); + nlocseg[dest] ++; + bufsize[dest] += 14; + } + } + + } + + + double ** segmbuf; + segmbuf = new double*[ntasks]; + for ( int i = 0; i < ntasks; i++) + if ( bufsize[i] > 0 ) + segmbuf[i] = new double[bufsize[i]]; + else + segmbuf[i] = 0; + + // cout << "mastermesh " << mastermesh -> GetNSeg() << " lineseg " << mastermesh -> LineSegment (1) << endl; + for ( int ls=1; ls <= mastermesh -> GetNSeg(); ls++) + { + Array volels; + mastermesh -> GetTopology().GetSegmentVolumeElements ( ls, volels ); + const Segment & seg = mastermesh -> LineSegment (ls); + int dest; + + for (int j = 0; j < volels.Size(); j++) + { + int ei = volels[j]; + int dest; + if ( ei > 0 && ei <= mastermesh->GetNE() ) + { + const Element & el = mastermesh -> VolumeElement (ei); + dest = el.GetPartition(); + + + if ( dest > 0 ) + { + segmbuf[dest][segi[dest]++] = ls; + segmbuf[dest][segi[dest]++] = seg.si; + segmbuf[dest][segi[dest]++] = seg.p1; + segmbuf[dest][segi[dest]++] = seg.p2; + segmbuf[dest][segi[dest]++] = seg.geominfo[0].trignum; + segmbuf[dest][segi[dest]++] = seg.geominfo[1].trignum; + segmbuf[dest][segi[dest]++] = seg.surfnr1; + segmbuf[dest][segi[dest]++] = seg.surfnr2; + segmbuf[dest][segi[dest]++] = seg.edgenr; + segmbuf[dest][segi[dest]++] = seg.epgeominfo[0].dist; + segmbuf[dest][segi[dest]++] = seg.epgeominfo[1].edgenr; + segmbuf[dest][segi[dest]++] = seg.epgeominfo[1].dist; + segmbuf[dest][segi[dest]++] = seg.singedge_right; + segmbuf[dest][segi[dest]++] = seg.singedge_left; + } + paralleltop -> SetDistantSegm ( dest, ls, int ( segi[dest] / 14 ) ); + } + } + } + + PrintMessage ( 3, "sending segments" ); + + for ( int dest = 1; dest < ntasks; dest++) + { + MyMPI_Send( segmbuf[dest], bufsize[dest], dest); + } + + for ( int dest = 0; dest < ntasks; dest++ ) + { + if ( segmbuf[dest] ) + delete [] segmbuf[dest]; + } + delete [] segmbuf; + + PrintMessage ( 3, "segments sent"); + + endtime = clock(); + (*testout) << "Sending Time segments = " << double(endtime - starttime)/CLOCKS_PER_SEC << endl; + starttime = clock(); + + + for ( int dest = 1; dest < ntasks; dest++ ) + MyMPI_Send("endmesh", dest); + + + for ( int dest = 1; dest < ntasks; dest ++ ) + { + MPI_Status status; + MPI_Wait (&sendrequest[dest], &status); + } + +#ifdef SCALASCA +#pragma pomp inst end(sendmesh) +#endif + } + + + + + + + + + void Mesh :: UpdateOverlap() + { + (*testout) << "UPDATE OVERLAP" << endl; + Array * globelnums; + +#ifdef SCALASCA +#pragma pomp inst begin(updateoverlap) +#endif + + paralleltop->IncreaseOverlap(); + + if ( id > 0 ) + { + int nvglob = paralleltop->GetNVGlob (); + int nelglob = paralleltop->GetNEGlob(); + + int nv = GetNV(); + int ned = topology -> GetNEdges(); + int nfa = topology -> GetNFaces(); + int nel = GetNE(); + + Array glob2loc_vert(nvglob); + glob2loc_vert = -1; + + for ( int locv = 1; locv <= GetNV(); locv++) + { + int globv = paralleltop->GetLoc2Glob_Vert(locv); + glob2loc_vert[globv] = locv; + } + + BitArray addedpoint ( paralleltop -> GetNVGlob () ); + BitArray addedel ( paralleltop -> GetNEGlob () ); + addedpoint.Clear(); + addedel.Clear(); + + Array distvert(ntasks), distel(ntasks), nsenddistel(ntasks); + + nsenddistel = 0; + + MyMPI_Allgather (GetNV(), distvert.Range(1, ntasks), MPI_HIGHORDER_COMM); + MyMPI_Allgather (GetNE(), distel.Range(1, ntasks), MPI_HIGHORDER_COMM); + + BitArray appendedpoint ( GetNP() * ntasks ); + appendedpoint.Clear(); + + TABLE sendpoints(ntasks); + TABLE sendelements(ntasks); + TABLE sendsel(ntasks); + + /* + TABLE cluster_neighbors(nv+ned+nfa+nel); + + for ( int node = 1; node <= nv+ned+nfa+nel; node++ ) + { + int cluster_rep; + cluster_rep = clusters->GetVertexRepresentant(node); + + if ( node == cluster_rep ) continue; + + Array dests; + int nneigh = 0; + if ( node - GetNV() <= 0 ) // cluster representant is vertex + { + int vert = node; + if ( paralleltop -> IsExchangeVert( vert ) ) + { + nneigh = paralleltop -> GetNDistantPNums(vert); + dests.SetSize(2*nneigh); + paralleltop -> GetDistantPNums ( vert, &dests[0] ); + } + } + else if ( node - GetNV() - ned <= 0 ) // cluster representant is edge + { + int edge = node - GetNV(); + if ( paralleltop -> IsExchangeEdge( edge ) ) + { + nneigh = paralleltop -> GetNDistantEdgeNums(edge); + dests.SetSize(2*nneigh); + paralleltop -> GetDistantEdgeNums ( edge, &dests[0] ); + } + } + else if ( node - GetNV() - ned - nfa <= 0 ) // cluster representant is face + { + int face = node - GetNV() - ned; + if ( paralleltop -> IsExchangeFace( face ) ) + { + nneigh = paralleltop -> GetNDistantFaceNums(face); + dests.SetSize(2*nneigh); + paralleltop -> GetDistantFaceNums ( face, &dests[0] ); + } + } + else // cluster representant is element + { + int el = node - GetNV() - ned - nfa; + if ( paralleltop -> IsExchangeElement( el ) ) + { + nneigh = paralleltop -> GetNDistantElNums(el); + dests.SetSize(2*nneigh); + paralleltop -> GetDistantElNums ( el, &dests[0] ); + } + } + + for ( int j = 1; j < nneigh; j++ ) + if ( !cluster_neighbors[cluster_rep].Contains ( dests[2*j] ) ) + { + cluster_neighbors.Add( cluster_rep-1, dests[2*j] ); + } + } + */ + + for ( int i = 0; i < ntasks; i++ ) + { + sendelements.Add (i, 0); + sendsel.Add(i, 0); + } + + Array nsentsel (ntasks); + nsentsel = 0; + + for ( int seli = 1; seli <= GetNSE(); seli++ ) + { + const Element2d & sel = SurfaceElement(seli); + int selnp = sel.GetNP(); + Array vert (selnp); + + Array alldests (0), dests; + + bool isparsel = false; + for ( int i = 0; i < selnp; i++ ) + { + vert[i] = sel.PNum(i+1); + if ( paralleltop -> IsExchangeVert ( vert[i] ) ) + { + isparsel = true; + paralleltop -> GetVertNeighbours ( vert[i], dests ); + for ( int j = 0; j < dests.Size(); j++ ) + if ( !alldests.Contains ( dests[j] ) ) + alldests.Append( dests[j] ); + } + } + + /* + int face = topology->GetSurfaceElementFace(seli); + int cluster_rep = clusters->GetFaceRepresentant(face); + if ( cluster_neighbors[cluster_rep-1].Size() > 0 ) + { + isparsel = true; + for ( int j = 0; j < cluster_neighbors[cluster_rep-1].Size(); j++ ) + if ( !alldests.Contains ( cluster_neighbors[cluster_rep-1][j] ) ) + alldests.Append( cluster_neighbors[cluster_rep-1][j] ); + + } + */ + + if ( !isparsel ) continue; + + for ( int i = 0; i < alldests.Size(); i ++ ) + { + // send the surface element to all distant procs: + + // loc number, + // number of points + // global vert numbers + // surface_element_index + + int dest = alldests[i]; + + // ***************** MISSING id = 0 + if ( dest == 0 ) continue; + + sendsel.Add (dest, seli); + sendsel.Add (dest, selnp); + + for ( int ii=0; ii GetLoc2Glob_Vert (vert[ii]) ); + } + sendsel.Add (dest, sel.GetIndex() ); + nsentsel[dest] ++; + } + } + for ( int dest = 1; dest < ntasks; dest++ ) + sendsel[dest][0] = nsentsel[dest]; + + for ( int eli = 1; eli <= GetNE(); eli++ ) + { + const Element & el = VolumeElement(eli); + int elnp = el.GetNP(); + Array vert (elnp); + + Array alldests (0), dests; + + for ( int i = 0; i < elnp; i++ ) + { + vert[i] = el.PNum(i+1); + if ( paralleltop -> IsExchangeVert ( vert[i] ) ) + { + paralleltop -> GetVertNeighbours ( vert[i], dests ); + for ( int j = 0; j < dests.Size(); j++ ) + if ( !alldests.Contains ( dests[j] ) ) + { + alldests.Append( dests[j] ); + paralleltop->SetExchangeElement ( dests[j], eli ); + } + paralleltop->SetExchangeElement ( eli ); + } + } + + /* + int cluster_rep = clusters->GetElementRepresentant(eli); + if ( cluster_neighbors[cluster_rep-1].Size() > 0 ) + { + for ( int j = 0; j < cluster_neighbors[cluster_rep-1].Size(); j++ ) + if ( !alldests.Contains ( cluster_neighbors[cluster_rep-1][j] ) ) + { + alldests.Append( cluster_neighbors[cluster_rep-1][j] ); + paralleltop->SetExchangeElement ( cluster_neighbors[cluster_rep-1][j], eli ); + } + paralleltop->SetExchangeElement ( eli ); + + } + */ + } + + for ( int eli = 1; eli <= GetNE(); eli++ ) + { + const Element & el = VolumeElement(eli); + int elnp = el.GetNP(); + Array vert (elnp); + + // append to point list: + // local pnum + // global pnum + // point coordinates + + Array points(elnp); + for ( int i = 0; i < elnp; i++ ) + { + vert[i] = el.PNum(i+1); + points[i] = Point(vert[i]); + Array knowndests; + // send point to all dests which get the volume element + for ( int dest = 0; dest < ntasks; dest ++ ) + { + // nur die neuen verts + if ( !paralleltop -> IsExchangeElement ( dest, eli ) ) continue; + + // jeder vertex nur ein mal + if ( appendedpoint.Test( (vert[i]-1) * ntasks + dest ) ) continue; + + appendedpoint.Set( (vert[i]-1) * ntasks + dest ); + paralleltop -> SetExchangeVert (dest, vert[i]); + paralleltop -> SetExchangeVert ( vert[i] ); + + + // append vertex to be sent + // loc pnum + // glob pnum + // coords + + // local pnum + sendpoints.Add (dest, vert[i]); + + // global pnum + sendpoints.Add (dest, paralleltop -> GetLoc2Glob_Vert ( vert[i] ) ); + // coordinates + sendpoints.Add (dest, points[i].X() ); + sendpoints.Add (dest, points[i].Y() ); + sendpoints.Add (dest, points[i].Z() ); + } + } + + + for ( int dest = 1; dest < ntasks; dest ++ ) + { + // send the volume element to all distant procs: + + // loc number, + // glob number + // number of points + // glob vertices + // element_index + if ( !paralleltop -> IsExchangeElement ( dest, eli ) ) continue; + + // loc number + sendelements.Add (dest, eli); + // glob number + sendelements.Add (dest, paralleltop -> GetLoc2Glob_VolEl(eli)); + + sendelements.Add (dest, elnp); + + for ( int j = 0; j < elnp; j++ ) + sendelements.Add (dest, paralleltop -> GetLoc2Glob_Vert(vert[j]) ); + sendelements.Add (dest, el.GetIndex() ); + + distel[dest]++; + nsenddistel[dest] ++; + paralleltop -> SetDistantEl ( dest, eli, -1);// distel[dest] ); + } + + + } + for ( int dest = 1; dest < ntasks; dest++ ) + if ( dest != id ) + sendelements[dest][0] = nsenddistel[dest]; + // find parallel surface elements, if there, append to sendsel - list + + + distel = 0; + + // sizes of sendpoints, sendelements, sendsels + Array sendsize_pts(ntasks), recvsize_pts(ntasks); + Array sendsize_els(ntasks), recvsize_els(ntasks); + Array sendsize_sels(ntasks), recvsize_sels(ntasks); + + for (int i = 0; i < ntasks; i++) + { + sendsize_pts[i] = sendpoints[i].Size(); + sendsize_els[i] = sendelements[i].Size(); + sendsize_sels[i] = sendsel[i].Size(); + } + + MyMPI_Alltoall (sendsize_pts.Range(1, ntasks), recvsize_pts.Range(1, ntasks), MPI_HIGHORDER_COMM); + MyMPI_Alltoall (sendsize_els.Range(1, ntasks), recvsize_els.Range(1, ntasks), MPI_HIGHORDER_COMM); + MyMPI_Alltoall (sendsize_sels.Range(1, ntasks), recvsize_sels.Range(1, ntasks), MPI_HIGHORDER_COMM); + recvsize_pts[0] = 0; + recvsize_els[0] = 0; + recvsize_sels[0] = 0; + + TABLE recvpoints(recvsize_pts); + TABLE recvelements(recvsize_els); + TABLE recvsel(recvsize_sels); + + recvpoints.SetElementSizesToMaxSizes (); + recvelements.SetElementSizesToMaxSizes (); + recvsel.SetElementSizesToMaxSizes (); + + /* + Array sendrequest(3*ntasks), recvrequest(3*ntasks); + Array status(3*ntasks); + + for ( int proc = 1; proc < ntasks; proc++) + { + MyMPI_ISend ( sendpoints[proc], proc, sendrequest[proc] ); + MyMPI_IRecv ( recvpoints[proc], proc, recvrequest[proc] ); + } + + for ( int proc = 1; proc < ntasks; proc++) + { + MPI_Wait(&sendrequest[proc], &status[proc]); + MPI_Wait(&recvrequest[proc], &status[proc]); + MyMPI_ISend ( sendelements[proc], proc, sendrequest[proc+1]); + MyMPI_IRecv ( recvelements[proc], proc, recvrequest[proc+1]); + } + for ( int proc = 1; proc < ntasks; proc++) + { + MPI_Wait(&sendrequest[proc+1], &status[proc+1]); + MPI_Wait(&recvrequest[proc+1], &status[proc+1]); + MyMPI_ISend ( sendsel[proc], proc, sendrequest[proc+2] ); + MyMPI_IRecv ( recvsel[proc], proc, recvrequest[proc+2] ); + } + for ( int proc = 1; proc < ntasks; proc++) + { + MPI_Wait(&sendrequest[proc+2], &status[proc+2]); + MPI_Wait(&recvrequest[proc+2], &status[proc+2]); + } + */ + + Array requests; + + for ( int proc = 1; proc < ntasks; proc++) + { + requests.Append (MyMPI_ISend ( sendpoints[proc], proc )); + requests.Append (MyMPI_IRecv ( recvpoints[proc], proc )); + } + + for ( int proc = 1; proc < ntasks; proc++) + { + requests.Append (MyMPI_ISend ( sendelements[proc], proc )); + requests.Append (MyMPI_IRecv ( recvelements[proc], proc )); + } + for ( int proc = 1; proc < ntasks; proc++) + { + requests.Append (MyMPI_ISend ( sendsel[proc], proc )); + requests.Append (MyMPI_IRecv ( recvsel[proc], proc )); + } + + MPI_Status stat; + for (int i = 0; i < requests.Size(); i++) + MPI_Wait (&requests[i], &stat); + + + + + + Array * distpnum2parpnum; + distpnum2parpnum = new Array [2]; + distpnum2parpnum[0].SetSize(0); + distpnum2parpnum[1].SetSize(0); + + Array firstdistpnum (ntasks); + + + for ( int sender = 1; sender < ntasks; sender++) + { + firstdistpnum[sender] = distpnum2parpnum[0].Size(); + + if ( sender == id ) continue; + + int ii = 0; + // receiving points + // dist pnum + // glob pnum + // coords + int numrecvpts = int ( recvpoints[sender].Size() / 5 ); + + paralleltop -> SetNV ( GetNV() + numrecvpts ); + int expectnp = GetNV () + numrecvpts; + + // received points + while ( ii < recvpoints[sender].Size() ) + { + int distpnum = int ( (recvpoints[sender])[ii++] ); + int globpnum = int ( (recvpoints[sender])[ii++] ); + + Point3d point; + point.X() = (recvpoints[sender])[ii++]; + point.Y() = (recvpoints[sender])[ii++]; + point.Z() = (recvpoints[sender])[ii++]; + + // append point as ghost + // if not already there + int pnum= glob2loc_vert[globpnum];//paralleltop -> Glob2Loc_Vert ( globpnum ); + if ( pnum <= 0 ) + { + pnum = AddPoint ( point, true ); + } + paralleltop -> SetDistantPNum ( 0, pnum, globpnum ); + glob2loc_vert[globpnum] = pnum; + paralleltop -> SetDistantPNum ( sender, pnum, distpnum ); + paralleltop -> SetExchangeVert ( pnum ); + } + + ii = 0; + + int recvnel = (recvelements[sender])[ii++]; + + paralleltop -> SetNE ( recvnel + GetNE() ); + + while ( ii < recvelements[sender].Size() ) + { + // receive list: + // distant number, + // glob number + // number of points + // glob vertices + // element_index + + int distelnum = (recvelements[sender])[ii++]; + int globelnum = (recvelements[sender])[ii++] ; + int elnp = (recvelements[sender])[ii++] ; + Array pnums(elnp), globpnums(elnp); + + // append volel + ELEMENT_TYPE eltype; + switch ( elnp ) + { + case 4: eltype = TET; break; + case 5: eltype = PYRAMID; break; + case 6: eltype = PRISM; break; + case 8: eltype = HEX; break; + } + + Element el ( eltype ) ; + + for ( int i = 0; i < elnp; i++ ) + { + globpnums[i] = int ( (recvelements[sender])[ii++] ); + pnums[i] = glob2loc_vert[globpnums[i]]; //paralleltop -> Glob2Loc_Vert(globpnums[i]); + } + + el.SetIndex ( (recvelements[sender])[ii++] ); + el.SetGhost ( 1 ); + + + for ( int i = 0; i < elnp; i++) + { + (int&) el[i] = pnums[i]; + } + + int eli = AddVolumeElement (el) + 1; + + paralleltop -> SetDistantEl ( sender, eli, distelnum); + paralleltop -> SetDistantEl ( 0, eli, globelnum ); + paralleltop -> SetExchangeElement ( eli ); + } + + ii = 0; + int nrecvsendsel = 0; + if ( recvsel[sender].Size() > 0 ) + nrecvsendsel = (recvsel[sender])[ii++]; + + paralleltop -> SetNSE ( nrecvsendsel + GetNSE() ); + + while ( ii < recvsel[sender] . Size() ) + { + // receive list: + // distant number, + // number of points + // global vert numbers + // surface_element_index + + int distselnum = (recvsel[sender])[ii++]; + int selnp = (recvsel[sender])[ii++] ; + + Array globpnums(selnp); + Array pnums(selnp); + + // append volel + ELEMENT_TYPE eltype; + switch ( selnp ) + { + case 4: eltype = QUAD; break; + case 3: eltype = TRIG; break; + } + + Element2d sel ( eltype ) ; + for ( int i = 0; i < selnp; i++ ) + { + globpnums[i] = int ( (recvsel[sender])[ii++] ); + pnums[i] = glob2loc_vert[globpnums[i]]; + } + + sel.SetIndex ( (recvsel[sender])[ii++] ); + sel.SetGhost ( 1 ); + + + for ( int i = 0; i < selnp; i++) + { + (int&) sel[i] = pnums[i]; + } + + int seli = AddSurfaceElement (sel); + + } + + + } + + delete [] distpnum2parpnum; + + } + + globelnums = new Array; + if ( id == 0 ) + { + for ( int dest = 1; dest < ntasks; dest++) + { + MyMPI_Recv ( *globelnums, dest ); + for ( int i = 0; i < globelnums->Size(); i++ ) + { + paralleltop -> SetDistantEl ( dest, (*globelnums)[i],i+1 ); + paralleltop -> SetExchangeElement ( dest, (*globelnums)[i] ); + } + } + } + else + { + globelnums -> SetSize(GetNE()); + for ( int i = 0; i < GetNE(); i++ ) + { + (*globelnums)[i] = paralleltop -> GetLoc2Glob_VolEl ( i+1 ); + } + MyMPI_Send ( *globelnums, 0 ); + } + + delete globelnums; + // send which elements are where + + topology -> Update(); + + // edge, facenums have probably changed as elements were added + // paralleltop has to be updated + + + paralleltop -> UpdateExchangeElements(); + + + paralleltop -> UpdateCoarseGridOverlap(); + //paralleltop -> UpdateTopology(); + +// *testout << "############################################" << endl << endl; +// paralleltop -> Print(); + + clusters -> Update(); + ; +#ifdef SCALASCA +#pragma pomp inst end(updateoverlap) +#endif + + } + + + + + + + + +} + + + +#endif diff --git a/libsrc/meshing/paralleltop.cpp b/libsrc/meshing/paralleltop.cpp new file mode 100644 index 00000000..4cd36b70 --- /dev/null +++ b/libsrc/meshing/paralleltop.cpp @@ -0,0 +1,1701 @@ +#ifdef PARALLEL + + +#include +#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 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 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 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 & 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 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 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 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 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 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 faces, edges; + Array 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 * sendarray, *recvarray; + sendarray = new Array (0); + recvarray = new Array; + + 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 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 * glob2locfa; + glob2locfa = new Array ( 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 edges; + int fa = locfa; + + Array 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 * sendarray, *recvarray; + sendarray = new Array (0); + recvarray = new Array; + + sendarray->SetSize (0); + + BitArray recvface(nfa); + recvface.Clear(); + + for ( int el = 1; el <= ne; el++ ) + { + Array 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 * glob2loc_el; + + glob2loc_el = new Array ( 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 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 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 * sendarray; + sendarray = new Array (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 * recvarray = new Array; + + 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 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 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 ** elementonproc, ** recvelonproc; + elementonproc = new Array*[ntasks]; + recvelonproc = new Array*[ntasks]; + + for ( int i = 1; i < ntasks; i++ ) + { + elementonproc[i] = new Array(0); + recvelonproc[i] = new Array(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 diff --git a/libsrc/meshing/paralleltop.hpp b/libsrc/meshing/paralleltop.hpp new file mode 100644 index 00000000..971394ca --- /dev/null +++ b/libsrc/meshing/paralleltop.hpp @@ -0,0 +1,269 @@ +#ifndef FILE_PARALLELTOP +#define FILE_PARALLELTOP + +#include + +extern int ntasks; + +class ParallelMeshTopology +{ + const Mesh & mesh; + + // number of local elements, vertices, points (?), edges, faces + int ne, nv, np, ned, nfa; + + // number of local segments and surface elements + int nseg, nsurfel; + + // number of global elements, vertices, ???, faces + int neglob, nvglob; + int nparel; + int nfaglob; + + /** + mapping from local to distant vertex number + each row of the table corresponds to one vertex + each row contains a list of pairs (procnr, dist_vnum) + */ + TABLE loc2distvert; + TABLE loc2distedge, loc2distface, loc2distel; + TABLE loc2distsegm, loc2distsurfel; + + BitArray * isexchangeface, * isexchangevert, * isexchangeedge, * isexchangeel; + + BitArray isghostedge, isghostface; + + bool coarseupdate; + int overlap; + +public: + + ParallelMeshTopology (const Mesh & amesh); + ~ParallelMeshTopology (); + + /// set number of local vertices, reset sizes of loc2dist_vert, isexchangevert... + void SetNV ( int anv ); + void SetNE ( int ane ); + void SetNSE ( int anse ); + void SetNSegm ( int anseg ); + + void Reset (); + + void SetLoc2Glob_Vert ( int locnum, int globnum ) { loc2distvert[locnum][0] = globnum; } + void SetLoc2Glob_VolEl ( int locnum, int globnum ) { loc2distel[locnum-1][0] = globnum; } + void SetLoc2Glob_SurfEl ( int locnum, int globnum ) { loc2distsurfel[locnum-1][0] = globnum; } + void SetLoc2Glob_Segm ( int locnum, int globnum ) { loc2distsegm[locnum-1][0] = globnum; } + + int GetLoc2Glob_Vert ( int locnum ) const { return loc2distvert[locnum][0]; } + int GetLoc2Glob_VolEl ( int locnum ) const { return loc2distel[locnum-1][0]; } + + void GetVertNeighbours ( int vnum, Array & dests ) const; + + int Glob2Loc_SurfEl ( int globnum ); + int Glob2Loc_VolEl ( int globnum ); + int Glob2Loc_Segm ( int globnum ); + int Glob2Loc_Vert ( int globnum ); + + int GetNDistantPNums ( int locpnum ) const { return loc2distvert[locpnum].Size() / 2 + 1; } + int GetNDistantFaceNums ( int locfacenum ) const { return loc2distface[locfacenum-1].Size() / 2 + 1; } + int GetNDistantEdgeNums ( int locedgenum ) const { return loc2distedge[locedgenum-1].Size() / 2 + 1; } + int GetNDistantElNums ( int locelnum ) const { return loc2distel[locelnum-1].Size() / 2 + 1; } + + int GetDistantPNum ( int proc, int locpnum ) const; + int GetDistantEdgeNum ( int proc, int locedgenum ) const; + int GetDistantFaceNum ( int proc, int locedgenum ) const; + int GetDistantElNum ( int proc, int locelnum ) const; + + int GetDistantPNums ( int locpnum, int * distpnums ) const; + int GetDistantEdgeNums ( int locedgenum, int * distedgenums ) const; + int GetDistantFaceNums ( int locedgenum, int * distfacenums ) const; + int GetDistantElNums ( int locelnum, int * distfacenums ) const; + + void Print() const; + + void SetExchangeFace ( int fnr ) { isexchangeface->Set( (fnr-1)*(ntasks+1) ); } + void SetExchangeVert ( int vnum ) { isexchangevert->Set( (vnum-1)*(ntasks+1) ); } + void SetExchangeElement ( int elnum ) { isexchangeel->Set((elnum-1)*(ntasks+1) ); } + void SetExchangeEdge ( int ednum ) { isexchangeedge->Set ((ednum-1)*(ntasks+1)); } + + + // fuer diese Fkts kann man ja O(N) - bitarrays lassen + bool IsExchangeVert ( PointIndex vnum ) const + { + return loc2distvert[vnum].Size() > 1; + // return isexchangevert->Test((vnum-1)*(ntasks+1)); + } + + bool IsExchangeEdge ( int ednum ) const + { + /* + if ( isexchangeedge->Test( (ednum-1)*(ntasks+1) ) != + ( loc2distedge[ednum-1].Size() > 1) ) + { + cerr << "isexedge differs, id = " << id << ", ednum = " << ednum << endl; + } + */ + // return loc2distedge[ednum-1].Size() > 1; + int i = (ednum-1)*(ntasks+1); + return isexchangeedge->Test(i); + } + + bool IsExchangeFace ( int fnum ) const + { + /* + if ( isexchangeface->Test( (fnum-1)*(ntasks+1) ) != + ( loc2distface[fnum-1].Size() > 1) ) + { + cerr << "it differs, id = " << id << ", fnum = " << fnum << endl; + } + */ + // nur hier funktioniert's so nicht ???? (JS) + // return loc2distface[fnum-1].Size() > 1; + return isexchangeface->Test( (fnum-1)*(ntasks+1) ); + } + + bool IsExchangeElement ( int elnum ) const + { + // return loc2distel[elnum-1].Size() > 1; + return isexchangeel->Test((elnum-1)*(ntasks+1)); + } + + bool IsExchangeSEl ( int selnum ) const + { + return loc2distsurfel[selnum-1].Size() > 1; + } + + + void SetExchangeFace (int dest, int fnr ) + { + // isexchangeface->Set((fnr-1)*(ntasks+1) + (dest+1)); + } + + void SetExchangeVert (int dest, int vnum ) + { + // isexchangevert->Set((vnum-1)*(ntasks+1) + (dest+1) ); + } + + void SetExchangeElement (int dest, int vnum ) + { + isexchangeel->Set( (vnum-1)*(ntasks+1) + (dest+1) ); + } + + void SetExchangeEdge (int dest, int ednum ) + { + // isexchangeedge->Set ( (ednum-1)*(ntasks+1) + (dest+1) ); + } + + + bool IsExchangeVert (int dest, int vnum ) const + { + FlatArray exchange = loc2distvert[vnum]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; + // return isexchangevert->Test((vnum-1)*(ntasks+1) + (dest+1) ); + } + + bool IsExchangeEdge (int dest, int ednum ) const + { + FlatArray exchange = loc2distedge[ednum-1]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; + + // return isexchangeedge->Test((ednum-1)*(ntasks+1) + (dest+1)); + } + + bool IsExchangeFace (int dest, int fnum ) const + { + FlatArray exchange = loc2distface[fnum-1]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; + + // return isexchangeface->Test((fnum-1)*(ntasks+1) + (dest+1) ); + } + + bool IsExchangeElement (int dest, int elnum ) const + { + // das geht auch nicht + /* + FlatArray exchange = loc2distel[elnum-1]; + for (int i = 1; i < exchange.Size(); i += 2) + if (exchange[i] == dest) return true; + return false; + */ + return isexchangeel->Test((elnum-1)*(ntasks+1) + (dest+1)); + } + + int Overlap() const + { return overlap; } + + int IncreaseOverlap () + { overlap++; return overlap; } + + void Update(); + + void UpdateCoarseGrid(); + void UpdateCoarseGridOverlap(); + void UpdateRefinement (); + void UpdateTopology (); + void UpdateExchangeElements(); + + void UpdateCoarseGridGlobal(); + + bool DoCoarseUpdate() const + { return !coarseupdate; } + + + void SetDistantFaceNum ( int dest, int locnum, int distnum ); + void SetDistantPNum ( int dest, int locnum, int distnum ); + void SetDistantEdgeNum ( int dest, int locnum, int distnum ); + void SetDistantEl ( int dest, int locnum, int distnum ); + void SetDistantSurfEl ( int dest, int locnum, int distnum ); + void SetDistantSegm ( int dest, int locnum, int distnum ); + + bool IsGhostEl ( int elnum ) const { return mesh.VolumeElement(elnum).IsGhost(); } + bool IsGhostFace ( int fanum ) const { return isghostface.Test(fanum-1); } + bool IsGhostEdge ( int ednum ) const { return isghostedge.Test(ednum-1); } + bool IsGhostVert ( int pnum ) const { return mesh.Point(pnum).IsGhost(); } + +// inline void SetGhostVert ( const int pnum ) +// { isghostvert.Set(pnum-1); } + + void SetGhostEdge ( int ednum ) + { isghostedge.Set(ednum-1); } + + void ClearGhostEdge ( int ednum ) + { isghostedge.Clear(ednum-1); } + + void SetGhostFace ( int fanum ) + { isghostface.Set(fanum-1); } + + void ClearGhostFace ( int fanum ) + { isghostface.Clear(fanum-1); } + + bool IsElementInPartition ( int elnum, int dest ) const + { + return IsExchangeElement ( dest, elnum); + } + + void SetElementInPartition ( int elnum, int dest ) + { + SetExchangeElement ( dest, elnum ); + } + + void SetNVGlob ( int anvglob ) { nvglob = anvglob; } + void SetNEGlob ( int aneglob ) { neglob = aneglob; } + + int GetNVGlob () { return nvglob; } + int GetNEGlob () { return neglob; } +}; + + + + + + + +#endif