mirror of
https://github.com/NGSolve/netgen.git
synced 2025-01-27 13:20:34 +05:00
Periodic Mesh somewhat working in 2 and 3 dimensions.
This commit is contained in:
Lukas
parent
6034f6ecc3
commit
29cfd7533c
@ -44,6 +44,7 @@ namespace netgen
|
|||||||
cd2names.SetSize(0);
|
cd2names.SetSize(0);
|
||||||
|
|
||||||
#ifdef PARALLEL
|
#ifdef PARALLEL
|
||||||
|
this->comm = MPI_COMM_WORLD;
|
||||||
paralleltop = new ParallelMeshTopology (*this);
|
paralleltop = new ParallelMeshTopology (*this);
|
||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
|
|||||||
@ -32,6 +32,11 @@ namespace netgen
|
|||||||
/// point coordinates
|
/// point coordinates
|
||||||
T_POINTS points;
|
T_POINTS points;
|
||||||
|
|
||||||
|
#ifdef PARALLEL
|
||||||
|
// The communicator for this mesh. (more or less dummy for now!)
|
||||||
|
MPI_Comm comm;
|
||||||
|
#endif
|
||||||
|
|
||||||
/// line-segments at edges
|
/// line-segments at edges
|
||||||
Array<Segment, 0, size_t> segments;
|
Array<Segment, 0, size_t> segments;
|
||||||
/// surface elements, 2d-inner elements
|
/// surface elements, 2d-inner elements
|
||||||
|
|||||||
@ -35,53 +35,25 @@ namespace netgen
|
|||||||
|
|
||||||
void Mesh :: SendRecvMesh ()
|
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)
|
if (id == 0)
|
||||||
PrintMessage (1, "Send/Receive mesh");
|
PrintMessage (1, "Send/Receive mesh");
|
||||||
|
|
||||||
{
|
// TODO: why is this here??
|
||||||
// distribute global information
|
if (id == 0)
|
||||||
int nelglob, nseglob, nvglob;
|
{
|
||||||
if (id == 0)
|
paralleltop -> SetNV (GetNV());
|
||||||
{
|
paralleltop -> SetNE (GetNE());
|
||||||
paralleltop -> SetNV (GetNV());
|
paralleltop -> SetNSegm (GetNSeg());
|
||||||
paralleltop -> SetNE (GetNE());
|
paralleltop -> SetNSE (GetNSE());
|
||||||
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<int> 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);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if (id == 0)
|
if (id == 0)
|
||||||
SendMesh ();
|
SendMesh ();
|
||||||
@ -101,19 +73,32 @@ namespace netgen
|
|||||||
{
|
{
|
||||||
Array<MPI_Request> sendrequests;
|
Array<MPI_Request> sendrequests;
|
||||||
|
|
||||||
PrintMessage ( 3, "Sending vertices");
|
int dim = GetDimension();
|
||||||
|
MyMPI_Bcast(dim);
|
||||||
|
|
||||||
|
|
||||||
|
const_cast<MeshTopology&>(GetTopology()).Update();
|
||||||
|
|
||||||
|
PrintMessage ( 3, "Sending nr of elements");
|
||||||
|
|
||||||
|
MPI_Comm comm = this->comm;
|
||||||
|
|
||||||
Array<int> num_els_on_proc(ntasks);
|
Array<int> num_els_on_proc(ntasks);
|
||||||
num_els_on_proc = 0;
|
num_els_on_proc = 0;
|
||||||
for (ElementIndex ei = 0; ei < GetNE(); ei++)
|
for (ElementIndex ei = 0; ei < GetNE(); ei++)
|
||||||
num_els_on_proc[(*this)[ei].GetPartition()]++;
|
num_els_on_proc[(*this)[ei].GetPartition()]++;
|
||||||
|
|
||||||
|
MPI_Scatter (&num_els_on_proc[0], 1, MPI_INT,
|
||||||
|
MPI_IN_PLACE, -1, MPI_INT, 0, comm);
|
||||||
|
|
||||||
TABLE<ElementIndex> els_of_proc (num_els_on_proc);
|
TABLE<ElementIndex> els_of_proc (num_els_on_proc);
|
||||||
for (ElementIndex ei = 0; ei < GetNE(); ei++)
|
for (ElementIndex ei = 0; ei < GetNE(); ei++)
|
||||||
els_of_proc.Add ( (*this)[ei].GetPartition(), ei);
|
els_of_proc.Add ( (*this)[ei].GetPartition(), ei);
|
||||||
|
|
||||||
|
PrintMessage ( 3, "Building vertex/proc mapping");
|
||||||
|
|
||||||
|
// TODO: sels/segs also have to consider identifications
|
||||||
|
|
||||||
Array<int> num_sels_on_proc(ntasks);
|
Array<int> num_sels_on_proc(ntasks);
|
||||||
num_sels_on_proc = 0;
|
num_sels_on_proc = 0;
|
||||||
for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
|
for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
|
||||||
@ -133,132 +118,160 @@ namespace netgen
|
|||||||
segs_of_proc.Add ( (*this)[ei].GetPartition(), ei);
|
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<int, PointIndex::BASE> vert_flag (GetNV());
|
Array<int, PointIndex::BASE> vert_flag (GetNV());
|
||||||
Array<int, PointIndex::BASE> num_procs_on_vert (GetNV());
|
Array<int, PointIndex::BASE> num_procs_on_vert (GetNV());
|
||||||
Array<int> num_verts_on_proc (ntasks);
|
Array<int> num_verts_on_proc (ntasks);
|
||||||
|
|
||||||
num_verts_on_proc = 0;
|
num_verts_on_proc = 0;
|
||||||
num_procs_on_vert = 0;
|
num_procs_on_vert = 0;
|
||||||
vert_flag = -1;
|
|
||||||
|
Array<INDEX_2> per_pairs;
|
||||||
|
Array<INDEX_2> pp2;
|
||||||
for (int dest = 1; dest < ntasks; dest++)
|
auto & idents = GetIdentifications();
|
||||||
|
bool has_periodic = false;
|
||||||
|
for (int idnr = 1; idnr < idents.GetMaxNr()+1; idnr++)
|
||||||
{
|
{
|
||||||
FlatArray<ElementIndex> els = els_of_proc[dest];
|
if(idents.GetType(idnr)!=Identifications::PERIODIC) continue;
|
||||||
for (int hi = 0; hi < els.Size(); hi++)
|
has_periodic = true;
|
||||||
{
|
idents.GetPairs(idnr, pp2);
|
||||||
const Element & el = (*this) [ els[hi] ];
|
per_pairs.Append(pp2);
|
||||||
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<SurfaceElementIndex> sels = sels_of_proc[dest];
|
|
||||||
for (int hi = 0; hi < sels.Size(); hi++)
|
|
||||||
{
|
|
||||||
const Element2d & el = (*this) [ sels[hi] ];
|
|
||||||
for (int i = 0; i < el.GetNP(); i++)
|
|
||||||
{
|
|
||||||
PointIndex epi = el[i];
|
|
||||||
if (vert_flag[epi] < dest)
|
|
||||||
{
|
|
||||||
vert_flag[epi] = dest;
|
|
||||||
|
|
||||||
num_verts_on_proc[dest]++;
|
|
||||||
num_procs_on_vert[epi]++;
|
|
||||||
|
|
||||||
paralleltop -> SetDistantPNum (dest, epi);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
FlatArray<SegmentIndex> 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);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
cout << "per_pairs: " << endl << per_pairs << endl;
|
||||||
|
Array<int, PointIndex::BASE> 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<PointIndex, PointIndex::BASE> 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<int> 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<allvs.Size(); j++)
|
||||||
|
{
|
||||||
|
auto pervs2 = per_verts[allvs[j]];
|
||||||
|
for (int l = 0; l < pervs2.Size(); l++)
|
||||||
|
{
|
||||||
|
auto addv = pervs2[l];
|
||||||
|
if (allvs.Contains(addv) || addv==k) continue;
|
||||||
|
changed = true;
|
||||||
|
allvs.Append(addv);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
f(k, allvs);
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
iterate_per_verts_trans([&](auto k, auto & allvs) {
|
||||||
|
npvs[k] = allvs.Size();
|
||||||
|
});
|
||||||
|
/** as per_verts, but TRANSITIVE **/
|
||||||
|
TABLE<PointIndex, PointIndex::BASE> per_verts_trans(npvs);
|
||||||
|
iterate_per_verts_trans([&](auto k, auto & allvs) {
|
||||||
|
for (int j = 0; j<allvs.Size(); j++)
|
||||||
|
per_verts_trans.Add(k, allvs[j]);
|
||||||
|
});
|
||||||
|
cout << "per_verts_trans: " << endl << per_verts_trans << endl;
|
||||||
|
for (int k = PointIndex::BASE; k < GetNV()+PointIndex::BASE; k++) {
|
||||||
|
BubbleSort(per_verts_trans[k]);
|
||||||
|
}
|
||||||
|
cout << "per_verts_trans: " << endl << per_verts_trans << endl;
|
||||||
|
|
||||||
|
|
||||||
|
/** iterates over all vertices, calls lambda function for (vertex, destination) **/
|
||||||
|
auto iterate_vertices = [&](auto f) {
|
||||||
|
vert_flag = -1;
|
||||||
|
for (int dest = 1; dest < ntasks; dest++)
|
||||||
|
{
|
||||||
|
FlatArray<ElementIndex> 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<SurfaceElementIndex> sels = sels_of_proc[dest];
|
||||||
|
for (int hi = 0; hi < sels.Size(); hi++)
|
||||||
|
{
|
||||||
|
const Element2d & el = (*this) [ sels[hi] ];
|
||||||
|
for (int i = 0; i < el.GetNP(); i++)
|
||||||
|
f(el[i], dest);
|
||||||
|
}
|
||||||
|
FlatArray<SegmentIndex> 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<PointIndex> verts_of_proc (num_verts_on_proc);
|
TABLE<PointIndex> verts_of_proc (num_verts_on_proc);
|
||||||
TABLE<int, PointIndex::BASE> procs_of_vert (num_procs_on_vert);
|
TABLE<int, PointIndex::BASE> procs_of_vert (num_procs_on_vert);
|
||||||
TABLE<int, PointIndex::BASE> loc_num_of_vert (num_procs_on_vert);
|
TABLE<int, PointIndex::BASE> loc_num_of_vert (num_procs_on_vert);
|
||||||
|
|
||||||
vert_flag = -1;
|
/** Write vertex/proc mappingfs to tables **/
|
||||||
for (int dest = 1; dest < ntasks; dest++)
|
iterate_vertices([&](auto vertex, auto dest) {
|
||||||
{
|
auto addit = [&] (auto vertex, auto dest) {
|
||||||
FlatArray<ElementIndex> els = els_of_proc[dest];
|
if (vert_flag[vertex] < dest)
|
||||||
for (int hi = 0; hi < els.Size(); hi++)
|
{
|
||||||
{
|
vert_flag[vertex] = dest;
|
||||||
const Element & el = (*this) [ els[hi] ];
|
procs_of_vert.Add (vertex, dest);
|
||||||
for (int i = 0; i < el.GetNP(); i++)
|
}
|
||||||
{
|
};
|
||||||
PointIndex epi = el[i];
|
addit(vertex, dest);
|
||||||
if (vert_flag[epi] < dest)
|
auto pers = per_verts_trans[vertex];
|
||||||
{
|
for(int j = 0; j < pers.Size(); j++)
|
||||||
vert_flag[epi] = dest;
|
addit(pers[j], dest);
|
||||||
procs_of_vert.Add (epi, dest);
|
});
|
||||||
}
|
|
||||||
}
|
/** local vertex numbers on distant procs (I think this was only used for debugging??) **/
|
||||||
}
|
|
||||||
|
|
||||||
FlatArray<SurfaceElementIndex> sels = sels_of_proc[dest];
|
|
||||||
for (int hi = 0; hi < sels.Size(); hi++)
|
|
||||||
{
|
|
||||||
const Element2d & el = (*this) [ sels[hi] ];
|
|
||||||
for (int i = 0; i < el.GetNP(); i++)
|
|
||||||
{
|
|
||||||
PointIndex epi = el[i];
|
|
||||||
if (vert_flag[epi] < dest)
|
|
||||||
{
|
|
||||||
vert_flag[epi] = dest;
|
|
||||||
procs_of_vert.Add (epi, dest);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
FlatArray<SegmentIndex> 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);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int vert = 1; vert <= GetNP(); vert++ )
|
for (int vert = 1; vert <= GetNP(); vert++ )
|
||||||
{
|
{
|
||||||
FlatArray<int> procs = procs_of_vert[vert];
|
FlatArray<int> procs = procs_of_vert[vert];
|
||||||
@ -270,6 +283,8 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
PrintMessage ( 3, "Sending Vertices - vertices");
|
||||||
|
|
||||||
for (int dest = 1; dest < ntasks; dest++)
|
for (int dest = 1; dest < ntasks; dest++)
|
||||||
{
|
{
|
||||||
FlatArray<PointIndex> verts = verts_of_proc[dest];
|
FlatArray<PointIndex> verts = verts_of_proc[dest];
|
||||||
@ -296,6 +311,71 @@ namespace netgen
|
|||||||
Array<int> num_distpnums(ntasks);
|
Array<int> num_distpnums(ntasks);
|
||||||
num_distpnums = 0;
|
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<int> 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<pp2.Size(); j++)
|
||||||
|
{
|
||||||
|
INDEX_2 & pair = pp2[j];
|
||||||
|
// both are on same procs!
|
||||||
|
auto ps = procs_of_vert[pair.I1()];
|
||||||
|
for (int l = 0; l < ps.Size(); l++)
|
||||||
|
{
|
||||||
|
ppd_sizes[ps[l]] += 2;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
TABLE<int> 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<pp2.Size(); j++)
|
||||||
|
{
|
||||||
|
INDEX_2 & pair = pp2[j];
|
||||||
|
auto ps = procs_of_vert[pair.I1()];
|
||||||
|
for (int l = 0; l < ps.Size(); l++)
|
||||||
|
{
|
||||||
|
auto p = ps[l];
|
||||||
|
pp_data[p][maxidentnr + idnr]++;
|
||||||
|
pp_data.Add(p, pair.I1());
|
||||||
|
pp_data.Add(p, pair.I2());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
cout << "periodic data: " << endl << pp_data << endl;
|
||||||
|
|
||||||
|
Array<MPI_Request> 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++)
|
for (int vert = 1; vert <= GetNP(); vert++)
|
||||||
{
|
{
|
||||||
FlatArray<int> procs = procs_of_vert[vert];
|
FlatArray<int> procs = procs_of_vert[vert];
|
||||||
@ -368,90 +448,341 @@ namespace netgen
|
|||||||
for (int dest = 1; dest < ntasks; dest++)
|
for (int dest = 1; dest < ntasks; dest++)
|
||||||
sendrequests.Append (MyMPI_ISend (fddata, dest, MPI_TAG_MESH+3));
|
sendrequests.Append (MyMPI_ISend (fddata, dest, MPI_TAG_MESH+3));
|
||||||
|
|
||||||
|
/** Surface Elements **/
|
||||||
|
|
||||||
PrintMessage ( 3, "Sending Surface elements" );
|
PrintMessage ( 3, "Sending Surface elements" );
|
||||||
|
// build sel-identification
|
||||||
Array <int> nlocsel(ntasks), bufsize(ntasks);
|
size_t nse = GetNSE();
|
||||||
|
Array<SurfaceElementIndex> 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<SurfaceElementIndex, 0> os1, os2;
|
||||||
|
for(SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
|
||||||
|
{
|
||||||
|
if(ided_sel[sei]!=-1) continue;
|
||||||
|
const Element2d & sel = (*this)[sei];
|
||||||
|
FlatArray<const PointIndex> 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; m<os1.Size(); m++) {
|
||||||
|
if(!os2.Contains(os1[m])) {
|
||||||
|
cout << "remove " << os1[m] << endl;
|
||||||
|
os1.Delete(m);
|
||||||
|
m--;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if(!os1.Size()) continue;
|
||||||
|
if(os1.Size()>1) {
|
||||||
|
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<const PointIndex> 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]<sei) continue;
|
||||||
|
cout << 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 <int> nlocsel(ntasks), bufsize(ntasks);
|
||||||
nlocsel = 0;
|
nlocsel = 0;
|
||||||
bufsize = 1;
|
bufsize = 1;
|
||||||
|
iterate_sels([&](SurfaceElementIndex sei, const Element2d & sel, int dest){
|
||||||
for (int sei = 1; sei <= GetNSE(); sei++ )
|
nlocsel[dest]++;
|
||||||
{
|
|
||||||
const Element2d & sel = SurfaceElement (sei);
|
|
||||||
int dest = sel.GetPartition();
|
|
||||||
nlocsel[dest] ++;
|
|
||||||
bufsize[dest] += 4 + 2*sel.GetNP();
|
bufsize[dest] += 4 + 2*sel.GetNP();
|
||||||
}
|
});
|
||||||
|
|
||||||
TABLE<int> selbuf(bufsize);
|
TABLE<int> selbuf(bufsize);
|
||||||
|
|
||||||
for (int dest = 1; dest < ntasks; dest++ )
|
for (int dest = 1; dest < ntasks; dest++ )
|
||||||
selbuf.Add (dest, nlocsel[dest]);
|
selbuf.Add (dest, nlocsel[dest]);
|
||||||
|
iterate_sels([&](SurfaceElementIndex sei, const auto & sel, int dest) {
|
||||||
for (int sei = 1; sei <= GetNSE(); sei ++ )
|
|
||||||
{
|
|
||||||
const Element2d & sel = SurfaceElement (sei);
|
|
||||||
int dest = sel.GetPartition();
|
|
||||||
|
|
||||||
selbuf.Add (dest, sei);
|
selbuf.Add (dest, sei);
|
||||||
selbuf.Add (dest, sel.GetIndex());
|
selbuf.Add (dest, sel.GetIndex());
|
||||||
// selbuf.Add (dest, 0);
|
// selbuf.Add (dest, 0);
|
||||||
selbuf.Add (dest, sel.GetNP());
|
selbuf.Add (dest, sel.GetNP());
|
||||||
|
|
||||||
for ( int ii = 1; ii <= sel.GetNP(); ii++)
|
for ( int ii = 1; ii <= sel.GetNP(); ii++)
|
||||||
{
|
{
|
||||||
selbuf.Add (dest, sel.PNum(ii));
|
selbuf.Add (dest, sel.PNum(ii));
|
||||||
selbuf.Add (dest, sel.GeomInfoPi(ii).trignum);
|
selbuf.Add (dest, sel.GeomInfoPi(ii).trignum);
|
||||||
}
|
}
|
||||||
}
|
});
|
||||||
|
cout << "Surface Element Data: " << endl << selbuf << endl;
|
||||||
|
// distribute sel data
|
||||||
for (int dest = 1; dest < ntasks; dest++)
|
for (int dest = 1; dest < ntasks; dest++)
|
||||||
sendrequests.Append (MyMPI_ISend(selbuf[dest], dest, MPI_TAG_MESH+4));
|
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");
|
PrintMessage ( 3, "Sending Edge Segments");
|
||||||
|
// segment identifications (cannot use transitive
|
||||||
|
// because false identifications)
|
||||||
Array <int> nloc(ntasks);
|
cout << "NSEG TOTAL: " << GetNSeg() << endl;
|
||||||
nloc = 0;
|
for(SegmentIndex segi = 0; segi < GetNSeg(); segi++)
|
||||||
bufsize = 1;
|
|
||||||
|
|
||||||
for (int i = 1; i <= GetNSeg(); i++ )
|
|
||||||
{
|
{
|
||||||
const Segment & seg = LineSegment (i);
|
cout << "segment " << segi << ": " << endl;
|
||||||
int dest = seg.GetPartition();
|
const Segment & seg = (*this)[segi];
|
||||||
nloc[dest] ++;
|
cout << seg[0] << " " << seg[1] << endl;
|
||||||
bufsize[dest] += 14;
|
|
||||||
}
|
}
|
||||||
|
// for(PointIndex pi = PointIndex::BASE; pi<GetNP(); pi++)
|
||||||
|
// {
|
||||||
|
// cout << endl << "---" << endl;
|
||||||
|
// cout << "segments for PI " << pi << ": " << endl;
|
||||||
|
// cout << GetTopology().GetVertexSegments(pi) << endl;
|
||||||
|
// cout << "---" << endl << endl;
|
||||||
|
// }
|
||||||
|
// Array<int, PointIndex::BASE> v2ss(GetNP());
|
||||||
|
// v2ss = 0;
|
||||||
|
// for(SegmentIndex segi = 0; segi < GetNSeg(); segi++)
|
||||||
|
// {
|
||||||
|
// const Segment & seg = (*this)[segi];
|
||||||
|
// v2ss[seg[0]]++;
|
||||||
|
// v2ss[seg[1]]++;
|
||||||
|
// }
|
||||||
|
// TABLE<SegmentIndex,PointIndex::BASE> 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<SegmentIndex> osegs1, osegs2, osegs_both;
|
||||||
|
Array<int> 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<ided1.Size(); l++)
|
||||||
|
{
|
||||||
|
// auto ospart = vert2segment[ided1[l]];
|
||||||
|
auto ospart = GetTopology().GetVertexSegments(ided1[l]);
|
||||||
|
for(int j=0; j<ospart.Size(); j++)
|
||||||
|
{
|
||||||
|
if(osegs1.Contains(ospart[j]))
|
||||||
|
throw NgException("double ospart?? contact Lukas...");
|
||||||
|
osegs1.Append(ospart[j]);
|
||||||
|
type1.Append(idents.GetSymmetric(pi1, ided1[l]));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
osegs2.SetSize(0);
|
||||||
|
type2.SetSize(0);
|
||||||
|
for (int l = 0; l<ided2.Size(); l++)
|
||||||
|
{
|
||||||
|
// auto ospart = vert2segment[ided2[l]];
|
||||||
|
auto ospart = GetTopology().GetVertexSegments(ided2[l]);
|
||||||
|
for(int j=0; j<ospart.Size(); j++)
|
||||||
|
{
|
||||||
|
if(osegs2.Contains(ospart[j]))
|
||||||
|
throw NgException("double ospart?? contact Lukas...");
|
||||||
|
osegs2.Append(ospart[j]);
|
||||||
|
type2.Append(idents.GetSymmetric(pi2, ided2[l]));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
osegs_both.SetSize(0);
|
||||||
|
for (int l = 0; l<osegs1.Size(); l++) {
|
||||||
|
auto pos = osegs2.Pos(osegs1[l]);
|
||||||
|
if (pos == -1) continue;
|
||||||
|
if (type1[l] != type2[pos]) continue;
|
||||||
|
osegs_both.Append(osegs1[l]);
|
||||||
|
}
|
||||||
|
for(int l = 0; l<osegs_both.Size(); l++) {
|
||||||
|
int segnp = (*this)[osegs_both[l]].GetNP();
|
||||||
|
if(segnp!=segnp)
|
||||||
|
throw NgException("Tried to identify non-curved and curved Segment!");
|
||||||
|
}
|
||||||
|
for(int l = 0; l<osegs_both.Size(); l++) {
|
||||||
|
f(segi, osegs_both[l]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
};
|
||||||
|
Array<int> per_seg_size(GetNSeg());
|
||||||
|
per_seg_size = 0;
|
||||||
|
iterate_segs1([&](SegmentIndex segi1, SegmentIndex segi2)
|
||||||
|
{ per_seg_size[segi1]++; });
|
||||||
|
TABLE<SegmentIndex> 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<SegmentIndex> 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<allsegs.Size(); j++)
|
||||||
|
{
|
||||||
|
auto persegs2 = per_seg[allsegs[j]];
|
||||||
|
for (int l = 0; l<persegs2.Size(); l++)
|
||||||
|
{
|
||||||
|
auto addseg = persegs2[l];
|
||||||
|
if (allsegs.Contains(addseg) || addseg==segi) continue;
|
||||||
|
allsegs.Append(addseg);
|
||||||
|
changed = true;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
f(segi, allsegs);
|
||||||
|
}
|
||||||
|
};
|
||||||
|
iterate_per_seg_trans([&](SegmentIndex segi, Array<SegmentIndex> & segs){
|
||||||
|
for (int j = 0; j < segs.Size(); j++)
|
||||||
|
per_seg_size[segi] = segs.Size();
|
||||||
|
});
|
||||||
|
TABLE<SegmentIndex> per_seg_trans(per_seg_size);
|
||||||
|
iterate_per_seg_trans([&](SegmentIndex segi, Array<SegmentIndex> & 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<int> dests;
|
||||||
|
auto iterate_segs2 = [&](auto f)
|
||||||
|
{
|
||||||
|
for (SegmentIndex segi = 0; segi<GetNSeg(); segi++)
|
||||||
|
{
|
||||||
|
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();
|
||||||
|
// cout << "ided with seg " << per_seg_trans[segi][l] << ", partition " << dest2 << endl;
|
||||||
|
if(!dests.Contains(dest2))
|
||||||
|
dests.Append(dest2);
|
||||||
|
}
|
||||||
|
for (int l = 0; l < dests.Size(); l++)
|
||||||
|
f(segi, seg, dests[l]);
|
||||||
|
}
|
||||||
|
};
|
||||||
|
Array<int> 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<double> segm_buf(bufsize);
|
TABLE<double> segm_buf(bufsize);
|
||||||
|
iterate_segs2([&](auto segi, const auto & seg, int dest)
|
||||||
/*
|
{
|
||||||
for (int dest = 1; dest < ntasks; dest++ )
|
segm_buf.Add (dest, segi);
|
||||||
segm_buf.Add (dest, nloc[dest]);
|
segm_buf.Add (dest, seg.si);
|
||||||
*/
|
segm_buf.Add (dest, seg.pnums[0]);
|
||||||
for (int i = 1; i <= GetNSeg(); i++ )
|
segm_buf.Add (dest, seg.pnums[1]);
|
||||||
{
|
segm_buf.Add (dest, seg.geominfo[0].trignum);
|
||||||
const Segment & seg = LineSegment (i);
|
segm_buf.Add (dest, seg.geominfo[1].trignum);
|
||||||
int dest = seg.GetPartition();
|
segm_buf.Add (dest, seg.surfnr1);
|
||||||
|
segm_buf.Add (dest, seg.surfnr2);
|
||||||
segm_buf.Add (dest, i);
|
segm_buf.Add (dest, seg.edgenr);
|
||||||
segm_buf.Add (dest, seg.si);
|
segm_buf.Add (dest, seg.epgeominfo[0].dist);
|
||||||
segm_buf.Add (dest, seg.pnums[0]);
|
segm_buf.Add (dest, seg.epgeominfo[1].edgenr);
|
||||||
segm_buf.Add (dest, seg.pnums[1]);
|
segm_buf.Add (dest, seg.epgeominfo[1].dist);
|
||||||
segm_buf.Add (dest, seg.geominfo[0].trignum);
|
segm_buf.Add (dest, seg.singedge_right);
|
||||||
segm_buf.Add (dest, seg.geominfo[1].trignum);
|
segm_buf.Add (dest, seg.singedge_left);
|
||||||
segm_buf.Add (dest, seg.surfnr1);
|
});
|
||||||
segm_buf.Add (dest, seg.surfnr2);
|
cout << "SEGMENT DATA: " << endl << segm_buf << endl;
|
||||||
segm_buf.Add (dest, seg.edgenr);
|
// distrubute segment data
|
||||||
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++)
|
for (int dest = 1; dest < ntasks; dest++)
|
||||||
sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5));
|
sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5));
|
||||||
|
// Segments done
|
||||||
|
|
||||||
|
|
||||||
|
// Array <int> 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<double> 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 <int> nlocseg(ntasks), segi(ntasks);
|
Array <int> nlocseg(ntasks), segi(ntasks);
|
||||||
@ -595,7 +926,16 @@ namespace netgen
|
|||||||
int timer_sels = NgProfiler::CreateTimer ("Receive surface elements");
|
int timer_sels = NgProfiler::CreateTimer ("Receive surface elements");
|
||||||
NgProfiler::RegionTimer reg(timer);
|
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;
|
// string st;
|
||||||
|
|
||||||
// receive vertices
|
// receive vertices
|
||||||
@ -604,6 +944,9 @@ namespace netgen
|
|||||||
Array<int> verts;
|
Array<int> verts;
|
||||||
MyMPI_Recv (verts, 0, MPI_TAG_MESH+1);
|
MyMPI_Recv (verts, 0, MPI_TAG_MESH+1);
|
||||||
|
|
||||||
|
|
||||||
|
cout << " got vertex data : " << endl << verts << endl;
|
||||||
|
|
||||||
int numvert = verts.Size();
|
int numvert = verts.Size();
|
||||||
paralleltop -> SetNV (numvert);
|
paralleltop -> SetNV (numvert);
|
||||||
|
|
||||||
@ -623,6 +966,34 @@ namespace netgen
|
|||||||
MPI_Datatype mptype = MeshPoint::MyGetMPIType();
|
MPI_Datatype mptype = MeshPoint::MyGetMPIType();
|
||||||
MPI_Status status;
|
MPI_Status status;
|
||||||
MPI_Recv( &points[1], numvert, mptype, 0, MPI_TAG_MESH+1, MPI_COMM_WORLD, &status);
|
MPI_Recv( &points[1], numvert, mptype, 0, MPI_TAG_MESH+1, MPI_COMM_WORLD, &status);
|
||||||
|
|
||||||
|
Array<int> 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<int> pairdata(2*npairs, &pp_data[offset]);
|
||||||
|
offset += 2*npairs;
|
||||||
|
for (int k = 0; k<npairs; k++) {
|
||||||
|
PointIndex loc1 = glob2loc_vert_ht.Get(pairdata[2*k]);
|
||||||
|
PointIndex loc2 = glob2loc_vert_ht.Get(pairdata[2*k+1]);
|
||||||
|
cout << loc1 << " <--> " << loc2 << " || ";
|
||||||
|
cout << pairdata[2*k] << " <--> " << pairdata[2*k+1] << endl;
|
||||||
|
idents.Add(loc1, loc2, idnr);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
Array<int> dist_pnums;
|
Array<int> dist_pnums;
|
||||||
MyMPI_Recv (dist_pnums, 0, MPI_TAG_MESH+1);
|
MyMPI_Recv (dist_pnums, 0, MPI_TAG_MESH+1);
|
||||||
@ -648,9 +1019,20 @@ namespace netgen
|
|||||||
|
|
||||||
el.SetIndex(elarray[ind++]);
|
el.SetIndex(elarray[ind++]);
|
||||||
el.SetNP(elarray[ind++]);
|
el.SetNP(elarray[ind++]);
|
||||||
|
|
||||||
|
int ind0 = ind;
|
||||||
|
int ind1 = ind;
|
||||||
for ( int j = 0; j < el.GetNP(); j++)
|
for ( int j = 0; j < el.GetNP(); j++)
|
||||||
el[j] = glob2loc_vert_ht.Get (elarray[ind++]);
|
el[j] = glob2loc_vert_ht.Get (elarray[ind++]);
|
||||||
|
|
||||||
|
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);
|
AddVolumeElement (el);
|
||||||
}
|
}
|
||||||
@ -689,12 +1071,25 @@ namespace netgen
|
|||||||
int nep = selbuf[ii++];
|
int nep = selbuf[ii++];
|
||||||
Element2d tri(nep);
|
Element2d tri(nep);
|
||||||
tri.SetIndex(faceind);
|
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++)
|
for(int j = 1; j <= nep; j++)
|
||||||
{
|
{
|
||||||
tri.PNum(j) = glob2loc_vert_ht.Get (selbuf[ii++]);
|
tri.PNum(j) = glob2loc_vert_ht.Get (selbuf[ii++]);
|
||||||
tri.GeomInfoPi(j).trignum = 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 );
|
paralleltop->SetLoc2Glob_SurfEl ( sel+1, globsel );
|
||||||
AddSurfaceElement (tri);
|
AddSurfaceElement (tri);
|
||||||
sel ++;
|
sel ++;
|
||||||
@ -718,9 +1113,13 @@ namespace netgen
|
|||||||
{
|
{
|
||||||
globsegi = int (segmbuf[ii++]);
|
globsegi = int (segmbuf[ii++]);
|
||||||
seg.si = 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[0] = glob2loc_vert_ht.Get (int(segmbuf[ii++]));
|
||||||
seg.pnums[1] = 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[0].trignum = int( segmbuf[ii++] );
|
||||||
seg.geominfo[1].trignum = int ( segmbuf[ii++]);
|
seg.geominfo[1].trignum = int ( segmbuf[ii++]);
|
||||||
seg.surfnr1 = int ( segmbuf[ii++]);
|
seg.surfnr1 = int ( segmbuf[ii++]);
|
||||||
@ -811,7 +1210,7 @@ namespace netgen
|
|||||||
SetNextMajorTimeStamp();
|
SetNextMajorTimeStamp();
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@ -241,7 +241,7 @@ namespace netgen
|
|||||||
// update new vertices after mesh-refinement
|
// update new vertices after mesh-refinement
|
||||||
if (mesh.mlbetweennodes.Size() > 0)
|
if (mesh.mlbetweennodes.Size() > 0)
|
||||||
{
|
{
|
||||||
// cout << "UpdateCoarseGrid - vertices" << endl;
|
cout << "UpdateCoarseGrid - vertices" << endl;
|
||||||
int newnv = mesh.mlbetweennodes.Size();
|
int newnv = mesh.mlbetweennodes.Size();
|
||||||
loc2distvert.ChangeSize(mesh.mlbetweennodes.Size());
|
loc2distvert.ChangeSize(mesh.mlbetweennodes.Size());
|
||||||
/*
|
/*
|
||||||
@ -376,7 +376,7 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
|
|
||||||
Array<int> sendarray, recvarray;
|
Array<int> sendarray, recvarray;
|
||||||
// cout << "UpdateCoarseGrid - edges" << endl;
|
cout << "UpdateCoarseGrid - edges" << endl;
|
||||||
|
|
||||||
// static int timerv = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex vertices");
|
// static int timerv = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex vertices");
|
||||||
static int timere = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex edges");
|
static int timere = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex edges");
|
||||||
@ -493,7 +493,7 @@ namespace netgen
|
|||||||
|
|
||||||
// MPI_Barrier (MPI_LocalComm);
|
// MPI_Barrier (MPI_LocalComm);
|
||||||
|
|
||||||
// cout << "UpdateCoarseGrid - faces" << endl;
|
cout << "UpdateCoarseGrid - faces" << endl;
|
||||||
if (mesh.GetDimension() == 3)
|
if (mesh.GetDimension() == 3)
|
||||||
{
|
{
|
||||||
NgProfiler::StartTimer (timerf);
|
NgProfiler::StartTimer (timerf);
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user