diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 4eac1c45..745adaef 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -62,7 +62,63 @@ namespace netgen delete vert2pointelement; } + template + void LoopOverEdges (const Mesh & mesh, MeshTopology & top, PointIndex v, + FUNC func) + { + for (ElementIndex elnr : top.GetVertexElements(v)) + { + const Element & el = mesh[elnr]; + int neledges = MeshTopology::GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + // edge.Sort(); + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + if (edge.I1() != v) continue; + + func (edge, elnr, k, 3, edgedir); + } + } + + for (SurfaceElementIndex elnr : top.GetVertexSurfaceElements(v)) + { + const Element2d & el = mesh[elnr]; + + int neledges = MeshTopology::GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + // edge.Sort(); + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + if (edge.I1() != v) continue; + + func (edge, elnr, k, 2, edgedir); + } + } + + for (SegmentIndex elnr : top.GetVertexSegments(v)) + { + const Segment & el = mesh[elnr]; + INDEX_2 edge(el[0], el[1]); + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + edge.Sort(); + if (edge.I1() != v) continue; + + func (edge, elnr, 0, 1, edgedir); + } + } + template void LoopOverFaces (const Mesh & mesh, MeshTopology & top, PointIndex v, FUNC func) @@ -435,163 +491,83 @@ namespace netgen } Array edgenr(nv); - Array edgeflag(nv); - Array vertex2; - - edgeflag = PointIndex::BASE-1; ned = edge2vert.Size(); + int max_edge_on_vertex = 0; for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++) { + int onv = vert2edge[i].Size() + vert2vertcoarse[i].Size() + + 4*(*vert2element)[i].Size() + 2*(*vert2surfelement)[i].Size() + (*vert2segment)[i].Size(); + max_edge_on_vertex = max (onv, max_edge_on_vertex); + } + + + INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); + Array vertex2; + + for (PointIndex v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) + { + v2eht.DeleteData(); vertex2.SetSize (0); - for (int j = 0; j < vert2edge[i].Size(); j++) + for (int j = 0; j < vert2edge[v].Size(); j++) { - int ednr = vert2edge[i][j]; + int ednr = vert2edge[v][j]; int i2 = edge2vert.Get(ednr)[1]; - edgeflag[i2] = i; + v2eht.Set (i2, v); edgenr[i2] = ednr; } - for (int j = 0; j < vert2vertcoarse[i].Size(); j++) + for (int j = 0; j < vert2vertcoarse[v].Size(); j++) { - int v2 = vert2vertcoarse[i][j]; - if (edgeflag[v2] < i) + int v2 = vert2vertcoarse[v][j]; + if (!v2eht.Used(v2)) { - edgeflag[v2] = i; + v2eht.Set (v2, v); vertex2.Append (v2); } } - FlatArray v2els = (*vert2element)[i]; - for (int j = 0; j < v2els.Size(); j++) - { - const Element & el = mesh[v2els[j]]; - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - edge.Sort(); - if (edge.I1() != i) continue; - - if (edgeflag[edge.I2()] < i) - { - vertex2.Append (edge.I2()); - edgeflag[edge.I2()] = i; - } - } - } - - for (int j = 0; j < (*vert2surfelement)[i].Size(); j++) - { - SurfaceElementIndex elnr = (*vert2surfelement)[i][j]; - const Element2d & el = mesh.SurfaceElement (elnr); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - edge.Sort(); - if (edge.I1() != i) continue; - - if (edgeflag[edge.I2()] < i) - { - vertex2.Append (edge.I2()); - edgeflag[edge.I2()] = i; - } - } - } - - for (int j = 0; j < (*vert2segment)[i].Size(); j++) - { - SegmentIndex elnr = (*vert2segment)[i][j]; - const Segment & el = mesh.LineSegment (elnr); - - INDEX_2 edge(el[0], el[1]); - edge.Sort(); - if (edge.I1() != i) continue; - - if (edgeflag[edge.I2()] < i) - { - vertex2.Append (edge.I2()); - edgeflag[edge.I2()] = i; - } - } + LoopOverEdges (mesh, *this, v, + [&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) + { + if (!v2eht.Used(edge.I2())) + { + vertex2.Append (edge.I2()); + v2eht.Set (edge.I2(), v); + } + }); QuickSort (vertex2); + for (int j = 0; j < vertex2.Size(); j++) { edgenr[vertex2[j]] = ++ned; - edge2vert.Append (INDEX_2 (i, vertex2[j])); + edge2vert.Append (INDEX_2 (v, vertex2[j])); } + LoopOverEdges (mesh, *this, v, + [&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) + { + int edgenum = edgenr[edge.I2()]; + switch (element_dim) + { + case 3: + edges[elnr][loc_edge].nr = edgenum-1; + edges[elnr][loc_edge].orient = edgedir; + break; + case 2: + surfedges[elnr][loc_edge].nr = edgenum-1; + surfedges[elnr][loc_edge].orient = edgedir; + break; + case 1: + segedges[elnr].nr = edgenum-1; + segedges[elnr].orient = edgedir; + break; + } + }); - for (int j = 0; j < (*vert2element)[i].Size(); j++) - { - ElementIndex elnr = (*vert2element)[i][j]; - const Element & el = mesh[elnr]; - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) continue; - - int edgenum = edgenr[edge.I2()]; - edges[elnr][k].nr = edgenum-1; - edges[elnr][k].orient = edgedir; - } - } - - for (int j = 0; j < (*vert2surfelement)[i].Size(); j++) - { - SurfaceElementIndex elnr = (*vert2surfelement)[i][j]; - const Element2d & el = mesh.SurfaceElement (elnr); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) continue; - - int edgenum = edgenr[edge.I2()]; - surfedges[elnr][k].nr = edgenum-1; - surfedges[elnr][k].orient = edgedir; - } - } - - for (int j = 0; j < (*vert2segment)[i].Size(); j++) - { - SegmentIndex elnr = (*vert2segment)[i][j]; - const Segment & el = mesh.LineSegment (elnr); - - INDEX_2 edge(el[0], el[1]); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) continue; - - int edgenum = edgenr[edge.I2()]; - segedges[elnr].nr = edgenum-1; - segedges[elnr].orient = edgedir; - } } }