parent edges for red refinement (thx Guosheng)

This commit is contained in:
Joachim Schoeberl 2021-03-31 07:50:18 +02:00
parent 44c10f663a
commit 096b419f6e
2 changed files with 259 additions and 183 deletions

View File

@ -34,6 +34,9 @@ namespace netgen
mesh.mlbetweennodes = INDEX_2(PointIndex::BASE-1,PointIndex::BASE-1); mesh.mlbetweennodes = INDEX_2(PointIndex::BASE-1,PointIndex::BASE-1);
} }
if (mesh.level_nv.Size() == 0)
mesh.level_nv.Append (mesh.GetNV());
INDEX_2_HASHTABLE<PointIndex> between(mesh.GetNP() + 5); INDEX_2_HASHTABLE<PointIndex> between(mesh.GetNP() + 5);
@ -739,6 +742,7 @@ namespace netgen
mesh.ComputeNVertices(); mesh.ComputeNVertices();
mesh.RebuildSurfaceElementLists(); mesh.RebuildSurfaceElementLists();
mesh.level_nv.Append (mesh.GetNV());
#ifdef PARALLEL #ifdef PARALLEL
if (mesh.GetCommunicator().Size() > 1) if (mesh.GetCommunicator().Size() > 1)

View File

@ -699,196 +699,268 @@ namespace netgen
} ); } );
if (build_parent_edges) if (build_parent_edges)
{ {
static Timer t("build_hierarchy"); RegionTimer reg(t); static Timer t("build_hierarchy"); RegionTimer reg(t);
cnt = 0; cnt = 0;
for (auto verts : edge2vert) cnt[verts[0]]++; for (auto verts : edge2vert) cnt[verts[0]]++;
TABLE<int,PointIndex::BASE> vert2edge (cnt); TABLE<int,PointIndex::BASE> vert2edge (cnt);
for (auto i : edge2vert.Range()) for (auto i : edge2vert.Range())
vert2edge.AddSave (edge2vert[i][0], i); vert2edge.AddSave (edge2vert[i][0], i);
// build edge hierarchy: // build edge hierarchy:
parent_edges.SetSize (ned); parent_edges.SetSize (ned);
parent_edges = { -1, { -1, -1, -1 } }; parent_edges = { -1, { -1, -1, -1 } };
for (size_t i = 0; i < ned; i++) for (size_t i = 0; i < ned; i++)
{ {
auto verts = edge2vert[i]; // 2 vertices of edge auto verts = edge2vert[i]; // 2 vertices of edge
if (verts[0] >= mesh->mlbetweennodes.Size()+PointIndex::BASE || if (verts[0] >= mesh->mlbetweennodes.Size()+PointIndex::BASE ||
verts[1] >= mesh->mlbetweennodes.Size()+PointIndex::BASE) verts[1] >= mesh->mlbetweennodes.Size()+PointIndex::BASE)
continue; continue;
auto pa0 = mesh->mlbetweennodes[verts[0]]; // two parent vertices of v0 auto pa0 = mesh->mlbetweennodes[verts[0]]; // two parent vertices of v0
auto pa1 = mesh->mlbetweennodes[verts[1]]; // two parent vertices of v1 auto pa1 = mesh->mlbetweennodes[verts[1]]; // two parent vertices of v1
// both vertices are on coarsest mesh // both vertices are on coarsest mesh
if (!pa0[0].IsValid() && !pa1[0].IsValid()) if (!pa0[0].IsValid() && !pa1[0].IsValid())
continue; continue;
int issplitedge = 0; int issplitedge = 0;
if (pa0[0] == verts[1] || pa0[1] == verts[1]) if (pa0[0] == verts[1] || pa0[1] == verts[1])
issplitedge = 1; issplitedge = 1;
if (pa1[0] == verts[0] || pa1[1] == verts[0]) if (pa1[0] == verts[0] || pa1[1] == verts[0])
issplitedge = 2; issplitedge = 2;
if (issplitedge) if (issplitedge)
{ {
// cout << "split edge " << endl; // cout << "split edge " << endl;
// edge is obtained by splitting one edge into two parts: // edge is obtained by splitting one edge into two parts:
auto paedge = issplitedge == 1 ? pa0 : pa1; auto paedge = issplitedge == 1 ? pa0 : pa1;
if (paedge[0] > paedge[1]) if (paedge[0] > paedge[1])
Swap (paedge[0], paedge[1]); Swap (paedge[0], paedge[1]);
for (int ednr : vert2edge[paedge[0]]) for (int ednr : vert2edge[paedge[0]])
if (auto cverts = edge2vert[ednr]; cverts[1] == paedge[1]) if (auto cverts = edge2vert[ednr]; cverts[1] == paedge[1])
{ {
int orient = (paedge[0] == verts[0] || paedge[1] == verts[1]) ? 1 : 0; int orient = (paedge[0] == verts[0] || paedge[1] == verts[1]) ? 1 : 0;
parent_edges[i] = { orient, { ednr, -1, -1 } }; parent_edges[i] = { orient, { ednr, -1, -1 } };
} }
} }
else else
{ {
// edge is splitting edge in middle of triangle: bool bisect_edge = false;
for (int j = 1; j <= 2; j++) // edge is splitting edge in middle of triangle:
{ for (int j = 1; j <= 2; j++)
INT<2> paedge1, paedge2, paedge3; {
int orient_inner = 0; INT<2> paedge1, paedge2, paedge3;
if (j == 1) int orient_inner = 0;
{ if (j == 1)
paedge1 = INT<2> (pa0[0], verts[1]); {
paedge2 = INT<2> (pa0[1], verts[1]); paedge1 = INT<2> (pa0[0], verts[1]);
paedge3 = INT<2> (pa0[0], pa0[1]); paedge2 = INT<2> (pa0[1], verts[1]);
orient_inner = 0; paedge3 = INT<2> (pa0[0], pa0[1]);
} orient_inner = 0;
else }
{ else
paedge1 = INT<2> (pa1[0], verts[0]); {
paedge2 = INT<2> (pa1[1], verts[0]); paedge1 = INT<2> (pa1[0], verts[0]);
paedge3 = INT<2> (pa1[0], pa1[1]); paedge2 = INT<2> (pa1[1], verts[0]);
orient_inner = 1; paedge3 = INT<2> (pa1[0], pa1[1]);
} orient_inner = 1;
if (paedge1[0] > paedge1[1]) }
Swap (paedge1[0], paedge1[1]); if (paedge1[0] > paedge1[1])
if (paedge2[0] > paedge2[1]) Swap (paedge1[0], paedge1[1]);
Swap (paedge2[0], paedge2[1]); if (paedge2[0] > paedge2[1])
if (paedge3[0] > paedge3[1]) Swap (paedge2[0], paedge2[1]);
Swap (paedge3[0], paedge3[1]); if (paedge3[0] > paedge3[1])
Swap (paedge3[0], paedge3[1]);
// if first vertex number is -1, then don't try to find entry in node2edge hash table // if first vertex number is -1, then don't try to find entry in node2edge hash table
if ( paedge1[0] == PointIndex::BASE-1 || paedge2[0] == PointIndex::BASE-1 ) if ( paedge1[0] == PointIndex::BASE-1 || paedge2[0] == PointIndex::BASE-1 )
continue; continue;
int paedgenr1=-1, paedgenr2=-1, paedgenr3=-1, orient1 = 0, orient2 = 0; int paedgenr1=-1, paedgenr2=-1, paedgenr3=-1, orient1 = 0, orient2 = 0;
for (int ednr : vert2edge[paedge1[0]]) for (int ednr : vert2edge[paedge1[0]])
if (auto cverts = edge2vert[ednr]; cverts[1] == paedge1[1]) if (auto cverts = edge2vert[ednr]; cverts[1] == paedge1[1])
{ {
paedgenr1 = ednr; paedgenr1 = ednr;
orient1 = (paedge1[0] == verts[0] || paedge1[1] == verts[1]) ? 1 : 0; orient1 = (paedge1[0] == verts[0] || paedge1[1] == verts[1]) ? 1 : 0;
} }
for (int ednr : vert2edge[paedge2[0]]) for (int ednr : vert2edge[paedge2[0]])
if (auto cverts = edge2vert[ednr]; cverts[1] == paedge2[1]) if (auto cverts = edge2vert[ednr]; cverts[1] == paedge2[1])
{ {
paedgenr2 = ednr; paedgenr2 = ednr;
orient2 = (paedge2[0] == verts[0] || paedge2[1] == verts[1]) ? 1 : 0; orient2 = (paedge2[0] == verts[0] || paedge2[1] == verts[1]) ? 1 : 0;
} }
for (int ednr : vert2edge[paedge3[0]]) for (int ednr : vert2edge[paedge3[0]])
if (auto cverts = edge2vert[ednr]; cverts[1] == paedge3[1]) if (auto cverts = edge2vert[ednr]; cverts[1] == paedge3[1])
paedgenr3 = ednr; paedgenr3 = ednr;
if (paedgenr1 != -1 && paedgenr2 != -1) if (paedgenr1 != -1 && paedgenr2 != -1){
parent_edges[i] = { orient1+2*orient2+4*orient_inner, { paedgenr1, paedgenr2, paedgenr3 } }; bisect_edge = true;
} parent_edges[i] = { orient1+2*orient2+4*orient_inner, { paedgenr1, paedgenr2, paedgenr3 } };
}
}
if (!bisect_edge) // not a bisect edge (then a red edge)
{
INT<2> paedge1, paedge2, paedge3;
int orient1 = 0, orient2 = 0, orient3=0;
int orient_inner = 0;
paedge1 = INT<2> (pa0[0], pa0[1]);
paedge2 = INT<2> (pa1[0], pa1[1]);
// find common vertex and the thrid pa edge
if (pa0[0]==pa1[0]){// 00
//orient1 = 0;
orient2 = 1;
if (pa0[1]<pa1[1]){
orient3 = 1;
paedge3 = INT<2> (pa0[1], pa1[1]);
}else{
//orient3 = 0;
paedge3 = INT<2> (pa1[1], pa0[1]);
}
}
else if (pa0[0]==pa1[1]){//01
//orient1 = 0;
//orient2 = 0;
if (pa0[1]<pa1[0]){
orient3 = 1;
paedge3 = INT<2> (pa0[1], pa1[0]);
}else{
//orient3 = 0;
paedge3 = INT<2> (pa1[0], pa0[1]);
}
}
else if (pa0[1]==pa1[0]){//10
orient1 = 1;
orient2 = 1;
if (pa0[0]<pa1[1]){
orient3 = 1;
paedge3 = INT<2> (pa0[0], pa1[1]);
}else{
//orient3 = 0;
paedge3 = INT<2> (pa1[1], pa0[0]);
}
}
else if (pa0[1]==pa1[1]){//11
orient1 = 1;
//orient2 = 0;
if (pa0[0]<pa1[0]){
orient3 = 1;
paedge3 = INT<2> (pa0[0], pa1[0]);
}else{
//orient3 = 0;
paedge3 = INT<2> (pa1[0], pa0[0]);
}
}
int paedgenr1=-1, paedgenr2=-1, paedgenr3=-1;
for (int ednr : vert2edge[paedge1[0]])
if (auto cverts = edge2vert[ednr]; cverts[1] == paedge1[1])
paedgenr1 = ednr;
for (int ednr : vert2edge[paedge2[0]])
if (auto cverts = edge2vert[ednr]; cverts[1] == paedge2[1])
paedgenr2 = ednr;
for (int ednr : vert2edge[paedge3[0]])
if (auto cverts = edge2vert[ednr]; cverts[1] == paedge3[1])
paedgenr3 = ednr;
parent_edges[i] = { 8+orient1+2*orient2+4*orient3, { paedgenr1, paedgenr2, paedgenr3 } };
//cout <<8+orient1+2*orient2+4*orient3 <<":"<<paedgenr1 <<", "<< paedgenr2 << ", "<< paedgenr3 << endl;
}
// TODO: quad edges // TODO: quad edges
/* /*
if (parentedges[i][0] == -1) if (parentedges[i][0] == -1)
{ {
// quad split // quad split
if (pa1[0] != pa2[0] && if (pa1[0] != pa2[0] &&
pa1[0] != pa2[1] && pa1[0] != pa2[1] &&
pa1[1] != pa2[0] && pa1[1] != pa2[0] &&
pa1[1] != pa2[1]) pa1[1] != pa2[1])
for (int j = 1; j <= 2; j++) for (int j = 1; j <= 2; j++)
{ {
INT<2> paedge1, paedge2; INT<2> paedge1, paedge2;
if (j == 1) if (j == 1)
{ {
paedge1 = INT<2> (pa1[0], pa2[0]); paedge1 = INT<2> (pa1[0], pa2[0]);
paedge2 = INT<2> (pa1[1], pa2[1]); paedge2 = INT<2> (pa1[1], pa2[1]);
} }
else else
{ {
paedge1 = INT<2> (pa1[0], pa2[1]); paedge1 = INT<2> (pa1[0], pa2[1]);
paedge2 = INT<2> (pa1[1], pa2[0]); paedge2 = INT<2> (pa1[1], pa2[0]);
} }
int paedgenr1 = 0, paedgenr2 = 0; int paedgenr1 = 0, paedgenr2 = 0;
int orient1 = 1, orient2 = 1; int orient1 = 1, orient2 = 1;
if (paedge1[0] > paedge1[1]) if (paedge1[0] > paedge1[1])
{ {
Swap (paedge1[0], paedge1[1]); Swap (paedge1[0], paedge1[1]);
orient1 = 0; orient1 = 0;
} }
if (paedge2[0] > paedge2[1]) if (paedge2[0] > paedge2[1])
{ {
Swap (paedge2[0], paedge2[1]); Swap (paedge2[0], paedge2[1]);
orient2 = 0; orient2 = 0;
} }
if ( paedge1[0] == -1 || paedge2[0] == -1 ) if ( paedge1[0] == -1 || paedge2[0] == -1 )
continue; continue;
if (node2edge.Used (paedge1) && node2edge.Used (paedge2)) if (node2edge.Used (paedge1) && node2edge.Used (paedge2))
{ {
paedgenr1 = node2edge.Get (paedge1); paedgenr1 = node2edge.Get (paedge1);
paedgenr2 = node2edge.Get (paedge2); paedgenr2 = node2edge.Get (paedge2);
parentedges[i][0] = 2 * paedgenr1 + orient1; parentedges[i][0] = 2 * paedgenr1 + orient1;
parentedges[i][1] = 2 * paedgenr2 + orient2; parentedges[i][1] = 2 * paedgenr2 + orient2;
} }
} }
} }
if (parentedges[i][0] == -1) if (parentedges[i][0] == -1)
{ {
// triangle split into quad+trig (from anisotropic pyramids) // triangle split into quad+trig (from anisotropic pyramids)
for (int j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
for (int k = 0; k < 2; k++) for (int k = 0; k < 2; k++)
{ {
INT<2> paedge (pa1[1-j], pa2[1-k]); INT<2> paedge (pa1[1-j], pa2[1-k]);
int orientpa = 1; int orientpa = 1;
if (paedge[0] > paedge[1]) if (paedge[0] > paedge[1])
{ {
Swap (paedge[0], paedge[1]); Swap (paedge[0], paedge[1]);
orientpa = 0; orientpa = 0;
} }
if (pa1[j] == pa2[k] && node2edge.Used(paedge)) if (pa1[j] == pa2[k] && node2edge.Used(paedge))
{ {
int paedgenr = node2edge.Get (paedge); int paedgenr = node2edge.Get (paedge);
parentedges[i][0] = 2 * paedgenr + orientpa; parentedges[i][0] = 2 * paedgenr + orientpa;
} }
} }
} }
*/ */
}
} }
} /*
for (int i : Range(parent_edges))
/* {
for (int i : Range(parent_edges)) auto [info, nrs] = parent_edges[i];
{ cout << "edge " << i << " has " << info << ", nrs = " << nrs[0] << " " << nrs[1] << endl;
auto [info, nrs] = parent_edges[i]; }
cout << "edge " << i << " has " << info << ", nrs = " << nrs[0] << " " << nrs[1] << endl; */
} }
*/
}
} }