From 29cfd7533ccd6c894404e2019493ae6dce5085e3 Mon Sep 17 00:00:00 2001 From: Lukas Date: Thu, 19 Jul 2018 17:33:26 +0200 Subject: [PATCH 1/6] Periodic Mesh somewhat working in 2 and 3 dimensions. --- libsrc/meshing/meshclass.cpp | 1 + libsrc/meshing/meshclass.hpp | 5 + libsrc/meshing/parallelmesh.cpp | 841 +++++++++++++++++++++++--------- libsrc/meshing/paralleltop.cpp | 6 +- 4 files changed, 629 insertions(+), 224 deletions(-) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index a89a6ff3..d49e37a4 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -44,6 +44,7 @@ namespace netgen cd2names.SetSize(0); #ifdef PARALLEL + this->comm = MPI_COMM_WORLD; paralleltop = new ParallelMeshTopology (*this); #endif } diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 1befedce..7b6ac3cd 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -32,6 +32,11 @@ namespace netgen /// point coordinates T_POINTS points; +#ifdef PARALLEL + // The communicator for this mesh. (more or less dummy for now!) + MPI_Comm comm; +#endif + /// line-segments at edges Array segments; /// surface elements, 2d-inner elements diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 50f64c1c..09dcad7c 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -35,53 +35,25 @@ namespace netgen void Mesh :: SendRecvMesh () { + int id, np; + MPI_Comm_rank(this->comm, &id); + MPI_Comm_size(this->comm, &np); + + if (np == 1) { + throw NgException("SendRecvMesh called, but only one rank in communicator!!"); + } + if (id == 0) PrintMessage (1, "Send/Receive mesh"); - { - // distribute global information - int nelglob, nseglob, nvglob; - if (id == 0) - { - paralleltop -> SetNV (GetNV()); - paralleltop -> SetNE (GetNE()); - paralleltop -> SetNSegm (GetNSeg()); - paralleltop -> SetNSE (GetNSE()); - - nelglob = GetNE(); - nseglob = GetNSE(); - nvglob = GetNV(); - } - - MyMPI_Bcast (dimension); - MyMPI_Bcast (nelglob); - MyMPI_Bcast (nseglob); - MyMPI_Bcast (nvglob); - } - - - - { - // distribute number of local elements - if (id == 0) - { - Array num_els_on_proc(ntasks); - num_els_on_proc = 0; - for (ElementIndex ei = 0; ei < GetNE(); ei++) - num_els_on_proc[(*this)[ei].GetPartition()]++; - - MPI_Scatter (&num_els_on_proc[0], 1, MPI_INT, - MPI_IN_PLACE, -1, MPI_INT, 0, MPI_COMM_WORLD); - } - else - { - int nelloc; - MPI_Scatter (NULL, 0, MPI_INT, - &nelloc, 1, MPI_INT, 0, MPI_COMM_WORLD); - - paralleltop -> SetNE (nelloc); - } - } + // TODO: why is this here?? + if (id == 0) + { + paralleltop -> SetNV (GetNV()); + paralleltop -> SetNE (GetNE()); + paralleltop -> SetNSegm (GetNSeg()); + paralleltop -> SetNSE (GetNSE()); + } if (id == 0) SendMesh (); @@ -101,19 +73,32 @@ namespace netgen { Array sendrequests; - PrintMessage ( 3, "Sending vertices"); + int dim = GetDimension(); + MyMPI_Bcast(dim); + + const_cast(GetTopology()).Update(); + + PrintMessage ( 3, "Sending nr of elements"); + + MPI_Comm comm = this->comm; Array num_els_on_proc(ntasks); num_els_on_proc = 0; for (ElementIndex ei = 0; ei < GetNE(); ei++) num_els_on_proc[(*this)[ei].GetPartition()]++; + MPI_Scatter (&num_els_on_proc[0], 1, MPI_INT, + MPI_IN_PLACE, -1, MPI_INT, 0, comm); + TABLE els_of_proc (num_els_on_proc); for (ElementIndex ei = 0; ei < GetNE(); ei++) els_of_proc.Add ( (*this)[ei].GetPartition(), ei); + PrintMessage ( 3, "Building vertex/proc mapping"); + // TODO: sels/segs also have to consider identifications + Array num_sels_on_proc(ntasks); num_sels_on_proc = 0; for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++) @@ -133,132 +118,160 @@ namespace netgen segs_of_proc.Add ( (*this)[ei].GetPartition(), ei); - + /** + Counts for proc->vertex, vertex->proc mapping. + + Whenever two vertices are identified by periodicity, any proc + that gets one of the vertices actually gets both of them. + This has to be transitive, that is, if + a <-> b and b <-> c, + then any proc that has vertex a also has vertices b and c! + **/ Array vert_flag (GetNV()); Array num_procs_on_vert (GetNV()); Array num_verts_on_proc (ntasks); - num_verts_on_proc = 0; num_procs_on_vert = 0; - vert_flag = -1; - - - for (int dest = 1; dest < ntasks; dest++) + + Array per_pairs; + Array pp2; + auto & idents = GetIdentifications(); + bool has_periodic = false; + for (int idnr = 1; idnr < idents.GetMaxNr()+1; idnr++) { - 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); - } - } - } - - - FlatArray sels = sels_of_proc[dest]; - for (int hi = 0; hi < sels.Size(); hi++) - { - const Element2d & el = (*this) [ sels[hi] ]; - for (int i = 0; i < el.GetNP(); i++) - { - PointIndex epi = el[i]; - if (vert_flag[epi] < dest) - { - vert_flag[epi] = dest; - - num_verts_on_proc[dest]++; - num_procs_on_vert[epi]++; - - paralleltop -> SetDistantPNum (dest, epi); - } - } - } - - FlatArray segs = segs_of_proc[dest]; - for (int hi = 0; hi < segs.Size(); hi++) - { - const Segment & el = (*this) [segs[hi]]; - for (int i = 0; i < 2; 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); - } - } - } + if(idents.GetType(idnr)!=Identifications::PERIODIC) continue; + has_periodic = true; + idents.GetPairs(idnr, pp2); + per_pairs.Append(pp2); } + + cout << "per_pairs: " << endl << per_pairs << endl; + Array npvs(GetNV()); + npvs = 0; + for (int k = 0; k < per_pairs.Size(); k++) { + npvs[per_pairs[k].I1()]++; + npvs[per_pairs[k].I2()]++; + } + /** for each vertex, gives us all identified vertices **/ + TABLE per_verts(npvs); + for (int k = 0; k < per_pairs.Size(); k++) { + per_verts.Add(per_pairs[k].I1(), per_pairs[k].I2()); + per_verts.Add(per_pairs[k].I2(), per_pairs[k].I1()); + } + for (int k = PointIndex::BASE; k < GetNV()+PointIndex::BASE; k++) { + BubbleSort(per_verts[k]); + } + cout << "per_verts: " << endl << per_verts << endl; + + + auto iterate_per_verts_trans = [&](auto f){ + Array allvs; + for (int k = PointIndex::BASE; k < GetNV()+PointIndex::BASE; k++) + { + allvs.SetSize(0); + allvs.Append(per_verts[k]); + bool changed = true; + while(changed) { + changed = false; + for (int j = 0; j per_verts_trans(npvs); + iterate_per_verts_trans([&](auto k, auto & allvs) { + for (int j = 0; j 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++) + f(el[i], dest); + } + FlatArray sels = sels_of_proc[dest]; + for (int hi = 0; hi < sels.Size(); hi++) + { + const Element2d & el = (*this) [ sels[hi] ]; + for (int i = 0; i < el.GetNP(); i++) + f(el[i], dest); + } + FlatArray segs = segs_of_proc[dest]; + for (int hi = 0; hi < segs.Size(); hi++) + { + const Segment & el = (*this) [segs[hi]]; + for (int i = 0; i < 2; i++) + f(el[i], dest); + } + } + }; + + /** count vertices per proc and procs per vertex **/ + iterate_vertices([&](auto vertex, auto dest){ + auto countit = [&] (auto vertex, auto dest) { + if (vert_flag[vertex] < dest) + { + vert_flag[vertex] = dest; + num_verts_on_proc[dest]++; + num_procs_on_vert[vertex]++; + GetParallelTopology().SetDistantPNum (dest, vertex); + } + }; + countit(vertex, dest); + auto pers = per_verts_trans[vertex]; + for(int j = 0; j < pers.Size(); j++) + countit(pers[j], 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 = (*this) [ els[hi] ]; - for (int i = 0; i < el.GetNP(); i++) - { - PointIndex epi = el[i]; - if (vert_flag[epi] < dest) - { - vert_flag[epi] = dest; - procs_of_vert.Add (epi, dest); - } - } - } - - FlatArray sels = sels_of_proc[dest]; - for (int hi = 0; hi < sels.Size(); hi++) - { - const Element2d & el = (*this) [ sels[hi] ]; - for (int i = 0; i < el.GetNP(); i++) - { - PointIndex epi = el[i]; - if (vert_flag[epi] < dest) - { - vert_flag[epi] = dest; - procs_of_vert.Add (epi, dest); - } - } - } - - FlatArray segs = segs_of_proc[dest]; - for (int hi = 0; hi < segs.Size(); hi++) - { - const Segment & el = (*this) [ segs[hi] ]; - for (int i = 0; i < 2; i++) - { - PointIndex epi = el[i]; - if (vert_flag[epi] < dest) - { - vert_flag[epi] = dest; - procs_of_vert.Add (epi, dest); - } - } - } - - } - + /** Write vertex/proc mappingfs to tables **/ + iterate_vertices([&](auto vertex, auto dest) { + auto addit = [&] (auto vertex, auto dest) { + if (vert_flag[vertex] < dest) + { + vert_flag[vertex] = dest; + procs_of_vert.Add (vertex, dest); + } + }; + addit(vertex, dest); + auto pers = per_verts_trans[vertex]; + for(int j = 0; j < pers.Size(); j++) + addit(pers[j], dest); + }); + + /** local vertex numbers on distant procs (I think this was only used for debugging??) **/ for (int vert = 1; vert <= GetNP(); vert++ ) { FlatArray procs = procs_of_vert[vert]; @@ -270,6 +283,8 @@ namespace netgen } } + PrintMessage ( 3, "Sending Vertices - vertices"); + for (int dest = 1; dest < ntasks; dest++) { FlatArray verts = verts_of_proc[dest]; @@ -296,6 +311,71 @@ namespace netgen Array num_distpnums(ntasks); num_distpnums = 0; + PrintMessage ( 3, "Sending Vertices - identifications"); + + /** + Info about periodic identifications sent to each proc is an array of + integers. + - maxidentnr + - type for each identification + - nr of pairs for each identification (each pair is local!) + - pairs for each periodic ident (global numbers) + **/ + int maxidentnr = idents.GetMaxNr(); + Array ppd_sizes(ntasks); + ppd_sizes = 1 + 2*maxidentnr; + for (int idnr = 1; idnr < idents.GetMaxNr()+1; idnr++) + { + if(idents.GetType(idnr)!=Identifications::PERIODIC) continue; + idents.GetPairs(idnr, pp2); + for(int j = 0; j pp_data(ppd_sizes); + for(int dest = 0; dest < ntasks; dest++) + pp_data.Add(dest, maxidentnr); + for (int dest = 0; dest < ntasks; dest++) + { + for (int idnr = 1; idnr < idents.GetMaxNr()+1; idnr++) + pp_data.Add(dest, idents.GetType(idnr)); + for (int idnr = 1; idnr < idents.GetMaxNr()+1; idnr++) + pp_data.Add(dest, 0); + } + for (int idnr = 1; idnr < idents.GetMaxNr()+1; idnr++) + { + if(idents.GetType(idnr)!=Identifications::PERIODIC) continue; + idents.GetPairs(idnr, pp2); + for(int j = 0; j req_per; + for(int dest = 1; dest < ntasks; dest++) + req_per.Append(MyMPI_ISend(pp_data[dest], dest, MPI_TAG_MESH+1)); + MPI_Waitall(req_per.Size(), &req_per[0], MPI_STATUS_IGNORE); + + PrintMessage ( 3, "Sending Vertices - distprocs"); + for (int vert = 1; vert <= GetNP(); vert++) { FlatArray procs = procs_of_vert[vert]; @@ -368,90 +448,341 @@ namespace netgen for (int dest = 1; dest < ntasks; dest++) sendrequests.Append (MyMPI_ISend (fddata, dest, MPI_TAG_MESH+3)); + /** Surface Elements **/ + PrintMessage ( 3, "Sending Surface elements" ); - - Array nlocsel(ntasks), bufsize(ntasks); + // build sel-identification + size_t nse = GetNSE(); + Array ided_sel(nse); + ided_sel = -1; + bool has_ided_sels = false; + if(GetNE() && has_periodic) //we can only have identified surf-els if we have vol-els (right?) + { + Array os1, os2; + for(SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + { + if(ided_sel[sei]!=-1) continue; + const Element2d & sel = (*this)[sei]; + FlatArray points = sel.PNums(); + auto ided1 = per_verts[points[0]]; + os1.SetSize(0); + for (int j = 0; j < ided1.Size(); j++) + os1.Append(GetTopology().GetVertexSurfaceElements(ided1[j])); + for (int j = 1; j < points.Size(); j++) + { + os2.SetSize(0); + auto p2 = points[j]; + auto ided2 = per_verts[p2]; + for (int l = 0; l < ided2.Size(); l++) + os2.Append(GetTopology().GetVertexSurfaceElements(ided2[l])); + for (int m = 0; m1) { + throw NgException("SurfaceElement identified with more than one other??"); + } + cout << "ID it with SEL " << os1[0] << " with points: " << endl; + const Element2d & sel2 = (*this)[sei]; + FlatArray points2 = sel2.PNums(); + for (int l = 0; l < points2.Size(); l++) + cout << points2[l] << " "; + cout << endl << endl; + has_ided_sels = true; + ided_sel[sei] = os1[0]; + ided_sel[os1[0]] = sei; + } + } + cout << endl << "IDED SELS: " << endl; + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++ ) + { + if(ided_sel[sei]==-1 || ided_sel[sei] " << ided_sel[sei] << endl; + } + cout << endl; + // build sel data to send + auto iterate_sels = [&](auto f) { + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++ ) + { + const Element2d & sel = (*this)[sei]; + int dest = (*this)[sei].GetPartition(); + f(sei, sel, dest); + if(ided_sel[sei]!=-1) + { + int dest2 = (*this)[ided_sel[sei]].GetPartition(); + f(sei, sel, dest2); + } + } + }; + Array nlocsel(ntasks), bufsize(ntasks); nlocsel = 0; bufsize = 1; - - for (int sei = 1; sei <= GetNSE(); sei++ ) - { - const Element2d & sel = SurfaceElement (sei); - int dest = sel.GetPartition(); - nlocsel[dest] ++; + iterate_sels([&](SurfaceElementIndex sei, const Element2d & sel, int dest){ + nlocsel[dest]++; bufsize[dest] += 4 + 2*sel.GetNP(); - } - + }); TABLE selbuf(bufsize); - for (int dest = 1; dest < ntasks; dest++ ) selbuf.Add (dest, nlocsel[dest]); - - for (int sei = 1; sei <= GetNSE(); sei ++ ) - { - const Element2d & sel = SurfaceElement (sei); - int dest = sel.GetPartition(); - + iterate_sels([&](SurfaceElementIndex sei, const auto & sel, int dest) { selbuf.Add (dest, sei); selbuf.Add (dest, sel.GetIndex()); // selbuf.Add (dest, 0); selbuf.Add (dest, sel.GetNP()); - for ( int ii = 1; ii <= sel.GetNP(); ii++) { selbuf.Add (dest, sel.PNum(ii)); selbuf.Add (dest, sel.GeomInfoPi(ii).trignum); } - } - + }); + cout << "Surface Element Data: " << endl << selbuf << endl; + // distribute sel data for (int dest = 1; dest < ntasks; dest++) sendrequests.Append (MyMPI_ISend(selbuf[dest], dest, MPI_TAG_MESH+4)); + + /** + ---Segments--- + Segments always have a twin (one for each surface). These will + be counted as identified, but whatever... + **/ PrintMessage ( 3, "Sending Edge Segments"); - - - Array nloc(ntasks); - nloc = 0; - bufsize = 1; - - for (int i = 1; i <= GetNSeg(); i++ ) + // segment identifications (cannot use transitive + // because false identifications) + cout << "NSEG TOTAL: " << GetNSeg() << endl; + for(SegmentIndex segi = 0; segi < GetNSeg(); segi++) { - const Segment & seg = LineSegment (i); - int dest = seg.GetPartition(); - nloc[dest] ++; - bufsize[dest] += 14; + cout << "segment " << segi << ": " << endl; + const Segment & seg = (*this)[segi]; + cout << seg[0] << " " << seg[1] << endl; } - + // for(PointIndex pi = PointIndex::BASE; pi v2ss(GetNP()); + // v2ss = 0; + // for(SegmentIndex segi = 0; segi < GetNSeg(); segi++) + // { + // const Segment & seg = (*this)[segi]; + // v2ss[seg[0]]++; + // v2ss[seg[1]]++; + // } + // TABLE vert2segment(v2ss); + // for(SegmentIndex segi = 0; segi < GetNSeg(); segi++) + // { + // const Segment & seg = (*this)[segi]; + // vert2segment.Add(seg[0], segi); + // vert2segment.Add(seg[1], segi); + // } + auto iterate_segs1 = [&](auto f) { + Array osegs1, osegs2, osegs_both; + Array type1, type2; + for(SegmentIndex segi = 0; segi < GetNSeg(); segi++) + { + const Segment & seg = (*this)[segi]; + int segnp = seg.GetNP(); + PointIndex pi1 = seg[0]; + // auto ided1 = per_verts_trans[pi1]; + auto ided1 = per_verts[pi1]; + PointIndex pi2 = seg[1]; + // auto ided2 = per_verts_trans[pi2]; + auto ided2 = per_verts[pi2]; + if (!(ided1.Size() && ided2.Size())) continue; + osegs1.SetSize(0); + type1.SetSize(0); + for (int l = 0; l per_seg_size(GetNSeg()); + per_seg_size = 0; + iterate_segs1([&](SegmentIndex segi1, SegmentIndex segi2) + { per_seg_size[segi1]++; }); + TABLE per_seg(per_seg_size); + iterate_segs1([&](SegmentIndex segi1, SegmentIndex segi2) + { per_seg.Add(segi1, segi2); }); + cout << endl << "IDED SEGS: " << endl << per_seg << endl; + // make per_seg transitive + auto iterate_per_seg_trans = [&](auto f){ + Array allsegs; + for (SegmentIndex segi = 0; segi < GetNSeg(); segi++) + { + allsegs.SetSize(0); + allsegs.Append(per_seg[segi]); + bool changed = true; + while (changed) + { + changed = false; + for (int j = 0; j & segs){ + for (int j = 0; j < segs.Size(); j++) + per_seg_size[segi] = segs.Size(); + }); + TABLE per_seg_trans(per_seg_size); + iterate_per_seg_trans([&](SegmentIndex segi, Array & segs){ + for (int j = 0; j < segs.Size(); j++) + per_seg_trans.Add(segi, segs[j]); + }); + cout << endl << "TRANS IDED SEGS: " << endl << per_seg_trans << endl; + // build segment data + Array dests; + auto iterate_segs2 = [&](auto f) + { + for (SegmentIndex segi = 0; segi nloc_seg(ntasks); + // bufsize = 1; //was originally this - why?? + bufsize = 0; + nloc_seg = 0; + iterate_segs2([&](auto segi, const auto & seg, int dest) + { + nloc_seg[dest]++; + bufsize[dest] += 14; + }); TABLE segm_buf(bufsize); - - /* - for (int dest = 1; dest < ntasks; dest++ ) - segm_buf.Add (dest, nloc[dest]); - */ - for (int i = 1; i <= GetNSeg(); i++ ) - { - const Segment & seg = LineSegment (i); - int dest = seg.GetPartition(); - - segm_buf.Add (dest, i); - segm_buf.Add (dest, seg.si); - segm_buf.Add (dest, seg.pnums[0]); - segm_buf.Add (dest, seg.pnums[1]); - segm_buf.Add (dest, seg.geominfo[0].trignum); - segm_buf.Add (dest, seg.geominfo[1].trignum); - segm_buf.Add (dest, seg.surfnr1); - segm_buf.Add (dest, seg.surfnr2); - segm_buf.Add (dest, seg.edgenr); - segm_buf.Add (dest, seg.epgeominfo[0].dist); - segm_buf.Add (dest, seg.epgeominfo[1].edgenr); - segm_buf.Add (dest, seg.epgeominfo[1].dist); - segm_buf.Add (dest, seg.singedge_right); - segm_buf.Add (dest, seg.singedge_left); - } - + iterate_segs2([&](auto segi, const auto & seg, int dest) + { + segm_buf.Add (dest, segi); + segm_buf.Add (dest, seg.si); + segm_buf.Add (dest, seg.pnums[0]); + segm_buf.Add (dest, seg.pnums[1]); + segm_buf.Add (dest, seg.geominfo[0].trignum); + segm_buf.Add (dest, seg.geominfo[1].trignum); + segm_buf.Add (dest, seg.surfnr1); + segm_buf.Add (dest, seg.surfnr2); + segm_buf.Add (dest, seg.edgenr); + segm_buf.Add (dest, seg.epgeominfo[0].dist); + segm_buf.Add (dest, seg.epgeominfo[1].edgenr); + segm_buf.Add (dest, seg.epgeominfo[1].dist); + segm_buf.Add (dest, seg.singedge_right); + segm_buf.Add (dest, seg.singedge_left); + }); + cout << "SEGMENT DATA: " << endl << segm_buf << endl; + // distrubute segment data for (int dest = 1; dest < ntasks; dest++) sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5)); + // Segments done + + + // Array nloc(ntasks); + // nloc = 0; + // bufsize = 1; + + // for (int i = 1; i <= GetNSeg(); i++ ) + // { + // const Segment & seg = LineSegment (i); + // int dest = seg.GetPartition(); + // nloc[dest] ++; + // bufsize[dest] += 14; + // } + + // TABLE segm_buf(bufsize); + + // /* + // for (int dest = 1; dest < ntasks; dest++ ) + // segm_buf.Add (dest, nloc[dest]); + // */ + // for (int i = 1; i <= GetNSeg(); i++ ) + // { + // const Segment & seg = LineSegment (i); + // int dest = seg.GetPartition(); + + // segm_buf.Add (dest, i); + // segm_buf.Add (dest, seg.si); + // segm_buf.Add (dest, seg.pnums[0]); + // segm_buf.Add (dest, seg.pnums[1]); + // segm_buf.Add (dest, seg.geominfo[0].trignum); + // segm_buf.Add (dest, seg.geominfo[1].trignum); + // segm_buf.Add (dest, seg.surfnr1); + // segm_buf.Add (dest, seg.surfnr2); + // segm_buf.Add (dest, seg.edgenr); + // segm_buf.Add (dest, seg.epgeominfo[0].dist); + // segm_buf.Add (dest, seg.epgeominfo[1].edgenr); + // segm_buf.Add (dest, seg.epgeominfo[1].dist); + // segm_buf.Add (dest, seg.singedge_right); + // segm_buf.Add (dest, seg.singedge_left); + // } + + // for (int dest = 1; dest < ntasks; dest++) + // sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5)); /* Array nlocseg(ntasks), segi(ntasks); @@ -595,7 +926,16 @@ namespace netgen int timer_sels = NgProfiler::CreateTimer ("Receive surface elements"); NgProfiler::RegionTimer reg(timer); - + int dim; + MyMPI_Bcast(dim); + SetDimension(dim); + + // Receive number of local elements + int nelloc; + MPI_Scatter (NULL, 0, MPI_INT, + &nelloc, 1, MPI_INT, 0, MPI_COMM_WORLD); + paralleltop -> SetNE (nelloc); + // string st; // receive vertices @@ -604,6 +944,9 @@ namespace netgen Array verts; MyMPI_Recv (verts, 0, MPI_TAG_MESH+1); + + cout << " got vertex data : " << endl << verts << endl; + int numvert = verts.Size(); paralleltop -> SetNV (numvert); @@ -623,6 +966,34 @@ namespace netgen MPI_Datatype mptype = MeshPoint::MyGetMPIType(); MPI_Status status; MPI_Recv( &points[1], numvert, mptype, 0, MPI_TAG_MESH+1, MPI_COMM_WORLD, &status); + + Array pp_data; + MyMPI_Recv(pp_data, 0, MPI_TAG_MESH+1); + + + int maxidentnr = pp_data[0]; + auto & idents = GetIdentifications(); + for (int idnr = 1; idnr < maxidentnr+1; idnr++) + { + + idents.SetType(idnr, (Identifications::ID_TYPE)pp_data[idnr]); + } + + cout << "vertex identifications: " << endl; + int offset = 2*maxidentnr+1; + for(int idnr = 1; idnr < maxidentnr+1; idnr++) + { + int npairs = pp_data[maxidentnr+idnr]; + FlatArray pairdata(2*npairs, &pp_data[offset]); + offset += 2*npairs; + for (int k = 0; k " << loc2 << " || "; + cout << pairdata[2*k] << " <--> " << pairdata[2*k+1] << endl; + idents.Add(loc1, loc2, idnr); + } + } Array dist_pnums; MyMPI_Recv (dist_pnums, 0, MPI_TAG_MESH+1); @@ -648,9 +1019,20 @@ namespace netgen el.SetIndex(elarray[ind++]); el.SetNP(elarray[ind++]); - + + int ind0 = ind; + int ind1 = ind; for ( int j = 0; j < el.GetNP(); j++) el[j] = glob2loc_vert_ht.Get (elarray[ind++]); + + cout << "Add element, glob: " << endl; + for ( int j = 0; j < el.GetNP(); j++) + cout << elarray[ind0++] << " "; + cout << endl; + cout << "-> LOC: " << endl; + for ( int j = 0; j < el.GetNP(); j++) + cout << glob2loc_vert_ht.Get (elarray[ind1++]) << " "; + cout << endl; AddVolumeElement (el); } @@ -689,12 +1071,25 @@ namespace netgen int nep = selbuf[ii++]; Element2d tri(nep); tri.SetIndex(faceind); + int ii0 = ii; + int ii1 = ii; + cout << "Add Surfel, glob: " << endl; + for(int j = 1; j <= nep; j++) { + cout << selbuf[ii0++] << " "; + ii0++; + } + cout << endl; for(int j = 1; j <= nep; j++) { tri.PNum(j) = glob2loc_vert_ht.Get (selbuf[ii++]); tri.GeomInfoPi(j).trignum = selbuf[ii++]; } - + cout << "-> LOC: " << endl; + for(int j = 1; j <= nep; j++) { + cout << selbuf[ii1++] << " "; + glob2loc_vert_ht.Get (selbuf[ii1++]); + } + cout << endl; paralleltop->SetLoc2Glob_SurfEl ( sel+1, globsel ); AddSurfaceElement (tri); sel ++; @@ -718,9 +1113,13 @@ namespace netgen { globsegi = int (segmbuf[ii++]); seg.si = int (segmbuf[ii++]); - + + cout << "add segm, glob: " << endl; + cout << segmbuf[ii] << " " << segmbuf[ii+1] << endl; seg.pnums[0] = glob2loc_vert_ht.Get (int(segmbuf[ii++])); seg.pnums[1] = glob2loc_vert_ht.Get (int(segmbuf[ii++])); + cout << "-> LOC: " << endl; + cout << seg.pnums[0] << " " << seg.pnums[1] << endl; seg.geominfo[0].trignum = int( segmbuf[ii++] ); seg.geominfo[1].trignum = int ( segmbuf[ii++]); seg.surfnr1 = int ( segmbuf[ii++]); @@ -811,7 +1210,7 @@ namespace netgen SetNextMajorTimeStamp(); } - + diff --git a/libsrc/meshing/paralleltop.cpp b/libsrc/meshing/paralleltop.cpp index 80c70121..2af5d0fe 100644 --- a/libsrc/meshing/paralleltop.cpp +++ b/libsrc/meshing/paralleltop.cpp @@ -241,7 +241,7 @@ namespace netgen // update new vertices after mesh-refinement if (mesh.mlbetweennodes.Size() > 0) { - // cout << "UpdateCoarseGrid - vertices" << endl; + cout << "UpdateCoarseGrid - vertices" << endl; int newnv = mesh.mlbetweennodes.Size(); loc2distvert.ChangeSize(mesh.mlbetweennodes.Size()); /* @@ -376,7 +376,7 @@ namespace netgen } Array sendarray, recvarray; - // cout << "UpdateCoarseGrid - edges" << endl; + cout << "UpdateCoarseGrid - edges" << endl; // static int timerv = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex vertices"); static int timere = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex edges"); @@ -493,7 +493,7 @@ namespace netgen // MPI_Barrier (MPI_LocalComm); - // cout << "UpdateCoarseGrid - faces" << endl; + cout << "UpdateCoarseGrid - faces" << endl; if (mesh.GetDimension() == 3) { NgProfiler::StartTimer (timerf); From 9fd51602b151a758ca76dc73703e457f8ca7c82e Mon Sep 17 00:00:00 2001 From: Lukas Date: Fri, 20 Jul 2018 13:29:16 +0200 Subject: [PATCH 2/6] less output --- libsrc/meshing/paralleltop.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libsrc/meshing/paralleltop.cpp b/libsrc/meshing/paralleltop.cpp index 2af5d0fe..80c70121 100644 --- a/libsrc/meshing/paralleltop.cpp +++ b/libsrc/meshing/paralleltop.cpp @@ -241,7 +241,7 @@ namespace netgen // update new vertices after mesh-refinement if (mesh.mlbetweennodes.Size() > 0) { - cout << "UpdateCoarseGrid - vertices" << endl; + // cout << "UpdateCoarseGrid - vertices" << endl; int newnv = mesh.mlbetweennodes.Size(); loc2distvert.ChangeSize(mesh.mlbetweennodes.Size()); /* @@ -376,7 +376,7 @@ namespace netgen } Array sendarray, recvarray; - cout << "UpdateCoarseGrid - edges" << endl; + // cout << "UpdateCoarseGrid - edges" << endl; // static int timerv = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex vertices"); static int timere = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex edges"); @@ -493,7 +493,7 @@ namespace netgen // MPI_Barrier (MPI_LocalComm); - cout << "UpdateCoarseGrid - faces" << endl; + // cout << "UpdateCoarseGrid - faces" << endl; if (mesh.GetDimension() == 3) { NgProfiler::StartTimer (timerf); From 7f424a543ac8657c18d45a184e059d6b9ff301d5 Mon Sep 17 00:00:00 2001 From: Lukas Date: Fri, 20 Jul 2018 13:44:46 +0200 Subject: [PATCH 3/6] periodic mpi formatted --- libsrc/meshing/parallelmesh.cpp | 234 +++++--------------------------- 1 file changed, 31 insertions(+), 203 deletions(-) diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 09dcad7c..7a3d1fec 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -97,8 +97,6 @@ namespace netgen PrintMessage ( 3, "Building vertex/proc mapping"); - // TODO: sels/segs also have to consider identifications - Array num_sels_on_proc(ntasks); num_sels_on_proc = 0; for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++) @@ -119,21 +117,24 @@ namespace netgen /** - Counts for proc->vertex, vertex->proc mapping. - + ----- STRATEGY FOR PERIODIC MESHES ----- + Whenever two vertices are identified by periodicity, any proc that gets one of the vertices actually gets both of them. This has to be transitive, that is, if a <-> b and b <-> c, then any proc that has vertex a also has vertices b and c! + + Surfaceelements and Segments that are identified by + periodicity are treated the same way. + + We need to duplicate these so we have containers to + hold the edges/facets. Afaik, a mesh cannot have nodes + that are not part of some sort of element. + **/ - Array vert_flag (GetNV()); - Array num_procs_on_vert (GetNV()); - Array num_verts_on_proc (ntasks); - num_verts_on_proc = 0; - num_procs_on_vert = 0; - + /** First, we build tables for vertex identification. **/ Array per_pairs; Array pp2; auto & idents = GetIdentifications(); @@ -145,14 +146,13 @@ namespace netgen idents.GetPairs(idnr, pp2); per_pairs.Append(pp2); } - - cout << "per_pairs: " << endl << per_pairs << endl; Array npvs(GetNV()); npvs = 0; for (int k = 0; k < per_pairs.Size(); k++) { npvs[per_pairs[k].I1()]++; npvs[per_pairs[k].I2()]++; } + /** for each vertex, gives us all identified vertices **/ TABLE per_verts(npvs); for (int k = 0; k < per_pairs.Size(); k++) { @@ -162,9 +162,8 @@ namespace netgen for (int k = PointIndex::BASE; k < GetNV()+PointIndex::BASE; k++) { BubbleSort(per_verts[k]); } - cout << "per_verts: " << endl << per_verts << endl; - - + + /** The same table as per_verts, but TRANSITIVE!! **/ auto iterate_per_verts_trans = [&](auto f){ Array allvs; for (int k = PointIndex::BASE; k < GetNV()+PointIndex::BASE; k++) @@ -189,24 +188,24 @@ namespace netgen f(k, allvs); } }; - iterate_per_verts_trans([&](auto k, auto & allvs) { npvs[k] = allvs.Size(); }); - /** as per_verts, but TRANSITIVE **/ TABLE per_verts_trans(npvs); iterate_per_verts_trans([&](auto k, auto & allvs) { for (int j = 0; j vert_flag (GetNV()); + Array num_procs_on_vert (GetNV()); + Array num_verts_on_proc (ntasks); + num_verts_on_proc = 0; + num_procs_on_vert = 0; auto iterate_vertices = [&](auto f) { vert_flag = -1; for (int dest = 1; dest < ntasks; dest++) @@ -234,7 +233,6 @@ namespace netgen } } }; - /** count vertices per proc and procs per vertex **/ iterate_vertices([&](auto vertex, auto dest){ auto countit = [&] (auto vertex, auto dest) { @@ -251,11 +249,9 @@ namespace netgen for(int j = 0; j < pers.Size(); j++) countit(pers[j], 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); - /** Write vertex/proc mappingfs to tables **/ iterate_vertices([&](auto vertex, auto dest) { auto addit = [&] (auto vertex, auto dest) { @@ -270,8 +266,10 @@ namespace netgen for(int j = 0; j < pers.Size(); j++) addit(pers[j], dest); }); - - /** local vertex numbers on distant procs (I think this was only used for debugging??) **/ + /** + local vertex numbers on distant procs + (I think this was only used for debugging??) + **/ for (int vert = 1; vert <= GetNP(); vert++ ) { FlatArray procs = procs_of_vert[vert]; @@ -282,7 +280,6 @@ namespace netgen loc_num_of_vert.Add (vert, verts_of_proc[dest].Size()); } } - PrintMessage ( 3, "Sending Vertices - vertices"); for (int dest = 1; dest < ntasks; dest++) @@ -311,9 +308,10 @@ namespace netgen Array num_distpnums(ntasks); num_distpnums = 0; - PrintMessage ( 3, "Sending Vertices - identifications"); /** + Next, we send the identifications themselfs. + Info about periodic identifications sent to each proc is an array of integers. - maxidentnr @@ -321,6 +319,7 @@ namespace netgen - nr of pairs for each identification (each pair is local!) - pairs for each periodic ident (global numbers) **/ + PrintMessage ( 3, "Sending Vertices - identifications"); int maxidentnr = idents.GetMaxNr(); Array ppd_sizes(ntasks); ppd_sizes = 1 + 2*maxidentnr; @@ -366,14 +365,11 @@ namespace netgen } } } - - cout << "periodic data: " << endl << pp_data << endl; - Array req_per; for(int dest = 1; dest < ntasks; dest++) req_per.Append(MyMPI_ISend(pp_data[dest], dest, MPI_TAG_MESH+1)); MPI_Waitall(req_per.Size(), &req_per[0], MPI_STATUS_IGNORE); - + PrintMessage ( 3, "Sending Vertices - distprocs"); for (int vert = 1; vert <= GetNP(); vert++) @@ -477,7 +473,6 @@ namespace netgen os2.Append(GetTopology().GetVertexSurfaceElements(ided2[l])); for (int m = 0; m1) { throw NgException("SurfaceElement identified with more than one other??"); } - cout << "ID it with SEL " << os1[0] << " with points: " << endl; const Element2d & sel2 = (*this)[sei]; FlatArray points2 = sel2.PNums(); for (int l = 0; l < points2.Size(); l++) @@ -498,13 +492,6 @@ namespace netgen ided_sel[os1[0]] = sei; } } - cout << endl << "IDED SELS: " << endl; - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++ ) - { - if(ided_sel[sei]==-1 || ided_sel[sei] " << ided_sel[sei] << endl; - } - cout << endl; // build sel data to send auto iterate_sels = [&](auto f) { for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++ ) @@ -540,50 +527,13 @@ namespace netgen selbuf.Add (dest, sel.GeomInfoPi(ii).trignum); } }); - cout << "Surface Element Data: " << endl << selbuf << endl; // distribute sel data for (int dest = 1; dest < ntasks; dest++) sendrequests.Append (MyMPI_ISend(selbuf[dest], dest, MPI_TAG_MESH+4)); - /** - ---Segments--- - - Segments always have a twin (one for each surface). These will - be counted as identified, but whatever... - **/ + /** Segments **/ PrintMessage ( 3, "Sending Edge Segments"); - // segment identifications (cannot use transitive - // because false identifications) - cout << "NSEG TOTAL: " << GetNSeg() << endl; - for(SegmentIndex segi = 0; segi < GetNSeg(); segi++) - { - cout << "segment " << segi << ": " << endl; - const Segment & seg = (*this)[segi]; - cout << seg[0] << " " << seg[1] << endl; - } - // for(PointIndex pi = PointIndex::BASE; pi v2ss(GetNP()); - // v2ss = 0; - // for(SegmentIndex segi = 0; segi < GetNSeg(); segi++) - // { - // const Segment & seg = (*this)[segi]; - // v2ss[seg[0]]++; - // v2ss[seg[1]]++; - // } - // TABLE vert2segment(v2ss); - // for(SegmentIndex segi = 0; segi < GetNSeg(); segi++) - // { - // const Segment & seg = (*this)[segi]; - // vert2segment.Add(seg[0], segi); - // vert2segment.Add(seg[1], segi); - // } auto iterate_segs1 = [&](auto f) { Array osegs1, osegs2, osegs_both; Array type1, type2; @@ -592,22 +542,19 @@ namespace netgen const Segment & seg = (*this)[segi]; int segnp = seg.GetNP(); PointIndex pi1 = seg[0]; - // auto ided1 = per_verts_trans[pi1]; auto ided1 = per_verts[pi1]; PointIndex pi2 = seg[1]; - // auto ided2 = per_verts_trans[pi2]; auto ided2 = per_verts[pi2]; if (!(ided1.Size() && ided2.Size())) continue; osegs1.SetSize(0); type1.SetSize(0); for (int l = 0; l per_seg(per_seg_size); iterate_segs1([&](SegmentIndex segi1, SegmentIndex segi2) { per_seg.Add(segi1, segi2); }); - cout << endl << "IDED SEGS: " << endl << per_seg << endl; // make per_seg transitive auto iterate_per_seg_trans = [&](auto f){ Array allsegs; @@ -686,7 +631,6 @@ namespace netgen for (int j = 0; j < segs.Size(); j++) per_seg_trans.Add(segi, segs[j]); }); - cout << endl << "TRANS IDED SEGS: " << endl << per_seg_trans << endl; // build segment data Array dests; auto iterate_segs2 = [&](auto f) @@ -700,7 +644,6 @@ namespace netgen for (int l = 0; l < per_seg_trans[segi].Size(); l++) { int dest2 = (*this)[per_seg_trans[segi][l]].GetPartition(); - // cout << "ided with seg " << per_seg_trans[segi][l] << ", partition " << dest2 << endl; if(!dests.Contains(dest2)) dests.Append(dest2); } @@ -735,124 +678,9 @@ namespace netgen segm_buf.Add (dest, seg.singedge_right); segm_buf.Add (dest, seg.singedge_left); }); - cout << "SEGMENT DATA: " << endl << segm_buf << endl; // distrubute segment data for (int dest = 1; dest < ntasks; dest++) sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5)); - // Segments done - - - // Array nloc(ntasks); - // nloc = 0; - // bufsize = 1; - - // for (int i = 1; i <= GetNSeg(); i++ ) - // { - // const Segment & seg = LineSegment (i); - // int dest = seg.GetPartition(); - // nloc[dest] ++; - // bufsize[dest] += 14; - // } - - // TABLE segm_buf(bufsize); - - // /* - // for (int dest = 1; dest < ntasks; dest++ ) - // segm_buf.Add (dest, nloc[dest]); - // */ - // for (int i = 1; i <= GetNSeg(); i++ ) - // { - // const Segment & seg = LineSegment (i); - // int dest = seg.GetPartition(); - - // segm_buf.Add (dest, i); - // segm_buf.Add (dest, seg.si); - // segm_buf.Add (dest, seg.pnums[0]); - // segm_buf.Add (dest, seg.pnums[1]); - // segm_buf.Add (dest, seg.geominfo[0].trignum); - // segm_buf.Add (dest, seg.geominfo[1].trignum); - // segm_buf.Add (dest, seg.surfnr1); - // segm_buf.Add (dest, seg.surfnr2); - // segm_buf.Add (dest, seg.edgenr); - // segm_buf.Add (dest, seg.epgeominfo[0].dist); - // segm_buf.Add (dest, seg.epgeominfo[1].edgenr); - // segm_buf.Add (dest, seg.epgeominfo[1].dist); - // segm_buf.Add (dest, seg.singedge_right); - // segm_buf.Add (dest, seg.singedge_left); - // } - - // for (int dest = 1; dest < ntasks; dest++) - // sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5)); - - /* - 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 <= GetNSeg(); segi ++ ) - { - Array volels; - const MeshTopology & topol = GetTopology(); - topol . GetSegmentSurfaceElements ( segi, volels ); - for (int j = 0; j < volels.Size(); j++) - { - int ei = volels[j]; - if ( ei > 0 && ei <= GetNSE() ) - { - const Element2d & el = SurfaceElement (ei); - int dest = el.GetPartition(); - nlocseg[dest] ++; - bufsize[dest] += 14; - } - } - } - - TABLE segmbuf(bufsize); - - for ( int ls=1; ls <= GetNSeg(); ls++) - { - Array volels; - GetTopology().GetSegmentSurfaceElements ( ls, volels ); - const Segment & seg = LineSegment (ls); - - for (int j = 0; j < volels.Size(); j++) - { - int ei = volels[j]; - if ( ei > 0 && ei <= GetNSE() ) - { - const Element2d & el = SurfaceElement (ei); - int dest = el.GetPartition(); - - if ( dest > 0 ) - { - segmbuf.Add (dest, ls); - segmbuf.Add (dest, seg.si); - segmbuf.Add (dest, seg.pnums[0]); - segmbuf.Add (dest, seg.pnums[1]); - segmbuf.Add (dest, seg.geominfo[0].trignum); - segmbuf.Add (dest, seg.geominfo[1].trignum); - segmbuf.Add (dest, seg.surfnr1); - segmbuf.Add (dest, seg.surfnr2); - segmbuf.Add (dest, seg.edgenr); - segmbuf.Add (dest, seg.epgeominfo[0].dist); - segmbuf.Add (dest, seg.epgeominfo[1].edgenr); - segmbuf.Add (dest, seg.epgeominfo[1].dist); - segmbuf.Add (dest, seg.singedge_right); - segmbuf.Add (dest, seg.singedge_left); - segi[dest] += 14; - } - // paralleltop -> SetDistantSegm ( dest, ls, int ( segi[dest] / 14 ) ); - } - } - } - - for ( int dest = 1; dest < ntasks; dest++) - sendrequests.Append (MyMPI_ISend(segmbuf[dest], dest, MPI_TAG_MESH+5)); - */ PrintMessage ( 3, "now wait ..."); From b2ae6210fc2d3cfcc263b3d07593c4bdf04e8f47 Mon Sep 17 00:00:00 2001 From: Lukas Date: Fri, 20 Jul 2018 13:56:27 +0200 Subject: [PATCH 4/6] removed more output --- libsrc/meshing/parallelmesh.cpp | 38 --------------------------------- 1 file changed, 38 deletions(-) diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 7a3d1fec..c329c3f0 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -484,9 +484,6 @@ namespace netgen } const Element2d & sel2 = (*this)[sei]; FlatArray points2 = sel2.PNums(); - for (int l = 0; l < points2.Size(); l++) - cout << points2[l] << " "; - cout << endl << endl; has_ided_sels = true; ided_sel[sei] = os1[0]; ided_sel[os1[0]] = sei; @@ -640,7 +637,6 @@ namespace netgen const Segment & seg = (*this)[segi]; dests.SetSize(0); dests.Append(seg.GetPartition()); - // cout << "partition of seg " << segi << " is " << dests[0] << endl; for (int l = 0; l < per_seg_trans[segi].Size(); l++) { int dest2 = (*this)[per_seg_trans[segi][l]].GetPartition(); @@ -773,8 +769,6 @@ namespace netgen MyMPI_Recv (verts, 0, MPI_TAG_MESH+1); - cout << " got vertex data : " << endl << verts << endl; - int numvert = verts.Size(); paralleltop -> SetNV (numvert); @@ -807,7 +801,6 @@ namespace netgen idents.SetType(idnr, (Identifications::ID_TYPE)pp_data[idnr]); } - cout << "vertex identifications: " << endl; int offset = 2*maxidentnr+1; for(int idnr = 1; idnr < maxidentnr+1; idnr++) { @@ -817,8 +810,6 @@ namespace netgen for (int k = 0; k " << loc2 << " || "; - cout << pairdata[2*k] << " <--> " << pairdata[2*k+1] << endl; idents.Add(loc1, loc2, idnr); } } @@ -848,19 +839,8 @@ namespace netgen el.SetIndex(elarray[ind++]); el.SetNP(elarray[ind++]); - int ind0 = ind; - int ind1 = ind; for ( int j = 0; j < el.GetNP(); j++) el[j] = glob2loc_vert_ht.Get (elarray[ind++]); - - cout << "Add element, glob: " << endl; - for ( int j = 0; j < el.GetNP(); j++) - cout << elarray[ind0++] << " "; - cout << endl; - cout << "-> LOC: " << endl; - for ( int j = 0; j < el.GetNP(); j++) - cout << glob2loc_vert_ht.Get (elarray[ind1++]) << " "; - cout << endl; AddVolumeElement (el); } @@ -899,25 +879,11 @@ namespace netgen int nep = selbuf[ii++]; Element2d tri(nep); tri.SetIndex(faceind); - int ii0 = ii; - int ii1 = ii; - cout << "Add Surfel, glob: " << endl; - for(int j = 1; j <= nep; j++) { - cout << selbuf[ii0++] << " "; - ii0++; - } - cout << endl; for(int j = 1; j <= nep; j++) { tri.PNum(j) = glob2loc_vert_ht.Get (selbuf[ii++]); tri.GeomInfoPi(j).trignum = selbuf[ii++]; } - cout << "-> LOC: " << endl; - for(int j = 1; j <= nep; j++) { - cout << selbuf[ii1++] << " "; - glob2loc_vert_ht.Get (selbuf[ii1++]); - } - cout << endl; paralleltop->SetLoc2Glob_SurfEl ( sel+1, globsel ); AddSurfaceElement (tri); sel ++; @@ -942,12 +908,8 @@ namespace netgen globsegi = int (segmbuf[ii++]); seg.si = int (segmbuf[ii++]); - cout << "add segm, glob: " << endl; - cout << segmbuf[ii] << " " << segmbuf[ii+1] << endl; seg.pnums[0] = glob2loc_vert_ht.Get (int(segmbuf[ii++])); seg.pnums[1] = glob2loc_vert_ht.Get (int(segmbuf[ii++])); - cout << "-> LOC: " << endl; - cout << seg.pnums[0] << " " << seg.pnums[1] << endl; seg.geominfo[0].trignum = int( segmbuf[ii++] ); seg.geominfo[1].trignum = int ( segmbuf[ii++]); seg.surfnr1 = int ( segmbuf[ii++]); From 53524579b7b4806a328fb4beca998fd19c575653 Mon Sep 17 00:00:00 2001 From: Lukas Date: Fri, 20 Jul 2018 14:00:34 +0200 Subject: [PATCH 5/6] Fixed indention --- libsrc/meshing/parallelmesh.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index c329c3f0..57ad1264 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -46,7 +46,7 @@ namespace netgen if (id == 0) PrintMessage (1, "Send/Receive mesh"); - // TODO: why is this here?? + // Why is this here?? if (id == 0) { paralleltop -> SetNV (GetNV()); @@ -838,7 +838,7 @@ namespace netgen el.SetIndex(elarray[ind++]); el.SetNP(elarray[ind++]); - + for ( int j = 0; j < el.GetNP(); j++) el[j] = glob2loc_vert_ht.Get (elarray[ind++]); @@ -907,7 +907,7 @@ namespace netgen { globsegi = int (segmbuf[ii++]); seg.si = int (segmbuf[ii++]); - + seg.pnums[0] = glob2loc_vert_ht.Get (int(segmbuf[ii++])); seg.pnums[1] = glob2loc_vert_ht.Get (int(segmbuf[ii++])); seg.geominfo[0].trignum = int( segmbuf[ii++] ); @@ -1000,7 +1000,7 @@ namespace netgen SetNextMajorTimeStamp(); } - + From 47e71acf135d0b49565580db1c8e2e0767507eae Mon Sep 17 00:00:00 2001 From: Lukas Date: Wed, 1 Aug 2018 10:35:26 +0200 Subject: [PATCH 6/6] Force segemnts to stick to surface elements in mesh partition. (surf els already stick to cells in 3d) --- libsrc/meshing/parallelmesh.cpp | 137 +++++++++++++++++++++++--------- 1 file changed, 101 insertions(+), 36 deletions(-) diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 57ad1264..f574f1f6 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -1092,51 +1092,87 @@ namespace netgen LineSegment(i+1).SetPartition(epart[i+GetNE()+GetNSE()] + 1); + // surface elements attached to volume elements + Array boundarypoints (GetNP()); + boundarypoints = false; + + if(GetDimension() == 3) + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + { + const Element2d & el = (*this)[sei]; + for (int j = 0; j < el.GetNP(); j++) + boundarypoints[el[j]] = true; + } + else + for (SegmentIndex segi = 0; segi < GetNSeg(); segi++) + { + const Segment & seg = (*this)[segi]; + for (int j = 0; j < 2; j++) + boundarypoints[seg[j]] = true; + } + + + // Build Pnt2Element table, boundary points only + Array cnt(GetNP()); + cnt = 0; + + auto loop_els_2d = [&](auto f) { + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + { + const Element2d & el = (*this)[sei]; + cout << "surf-el " << sei << " verts: " << endl; + for (int j = 0; j < el.GetNP(); j++) { + cout << el[j] << " "; + f(el[j], sei); + } + cout << endl; + } + }; + auto loop_els_3d = [&](auto f) { + for (ElementIndex ei = 0; ei < GetNE(); ei++) + { + const Element & el = (*this)[ei]; + for (int j = 0; j < el.GetNP(); j++) + f(el[j], ei); + } + }; + auto loop_els = [&](auto f) + { + if (GetDimension() == 3 ) + loop_els_3d(f); + else + loop_els_2d(f); + }; + + + loop_els([&](auto vertex, int index) + { + if(boundarypoints[vertex]) + cnt[vertex]++; + }); + cout << "count: " << endl << cnt << endl; + TABLE pnt2el(cnt); + loop_els([&](auto vertex, int index) + { + if(boundarypoints[vertex]) + pnt2el.Add(vertex, index); + }); + if (GetDimension() == 3) { - - // surface elements attached to volume elements - Array boundarypoints (GetNP()); - boundarypoints = false; - - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) - { - const Element2d & el = (*this)[sei]; - for (int j = 0; j < el.GetNP(); j++) - boundarypoints[el[j]] = true; - } - // Build Pnt2Element table, boundary points only - Array cnt(GetNP()); - cnt = 0; - for (ElementIndex ei = 0; ei < GetNE(); ei++) - { - const Element & el = (*this)[ei]; - for (int j = 0; j < el.GetNP(); j++) - if (boundarypoints[el[j]]) - cnt[el[j]]++; - } - TABLE pnt2el(cnt); - cnt = 0; - for (ElementIndex ei = 0; ei < GetNE(); ei++) - { - const Element & el = (*this)[ei]; - for (int j = 0; j < el.GetNP(); j++) - if (boundarypoints[el[j]]) - pnt2el.Add (el[j], ei); - } - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { Element2d & sel = (*this)[sei]; PointIndex pi1 = sel[0]; - FlatArray els = pnt2el[pi1]; + // FlatArray els = pnt2el[pi1]; + FlatArray els = pnt2el[pi1]; sel.SetPartition (-1); for (int j = 0; j < els.Size(); j++) { - const Element & el = (*this)[els[j]]; + const Element & el = (*this)[ElementIndex(els[j])]; bool hasall = true; @@ -1165,13 +1201,13 @@ namespace netgen { Segment & sel = (*this)[si]; PointIndex pi1 = sel[0]; - FlatArray els = pnt2el[pi1]; + FlatArray els = pnt2el[pi1]; sel.SetPartition (-1); for (int j = 0; j < els.Size(); j++) { - const Element & el = (*this)[els[j]]; + const Element & el = (*this)[ElementIndex(els[j])]; bool haspi[9] = { false }; // max surfnp @@ -1193,7 +1229,36 @@ namespace netgen if (sel.GetPartition() == -1) cerr << "no volume element found" << endl; } - + } + else + { + for (SegmentIndex segi = 0; segi < GetNSeg(); segi++) + { + Segment & seg = (*this)[segi]; + seg.SetPartition(-1); + PointIndex pi1 = seg[0]; + + FlatArray sels = pnt2el[pi1]; + for (int j = 0; j < sels.Size(); j++) + { + SurfaceElementIndex sei = sels[j]; + Element2d & se = (*this)[sei]; + bool found = false; + for (int l = 0; l < se.GetNP(); l++ && !found) + found |= (se[l]==seg[1]); + if(found) { + seg.SetPartition(se.GetPartition()); + break; + } + } + + if (seg.GetPartition() == -1) { + cout << endl << "segi: " << segi << endl; + cout << "points: " << seg[0] << " " << seg[1] << endl; + cout << "surfels: " << endl << sels << endl; + throw NgException("no surface element found"); + } + } } }