From 4987d12556c44664676d9928e12ddfd2f637c360 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 27 Sep 2019 20:49:12 +0200 Subject: [PATCH 1/2] SplitImprove for triangles --- libsrc/meshing/improve2.cpp | 143 ++++++++++++++++++++++++++++++++++- libsrc/meshing/improve2.hpp | 1 + libsrc/meshing/meshclass.cpp | 24 ++++++ libsrc/meshing/meshclass.hpp | 1 + 4 files changed, 168 insertions(+), 1 deletion(-) diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp index be431f87..fda88898 100644 --- a/libsrc/meshing/improve2.cpp +++ b/libsrc/meshing/improve2.cpp @@ -432,6 +432,7 @@ namespace netgen { if (!faceindex) { + SplitImprove(mesh); PrintMessage (3, "Combine improve"); for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++) @@ -696,7 +697,11 @@ namespace netgen if ( (normals[el[l]] * nv) < 0.5) bad2 += 1e10; - illegal2 += 1-mesh.LegalTrig(el); + Element2d el1 = el; + for (auto i : Range(3)) + if(el1[i]==pi2) + el1[i] = pi1; + illegal2 += 1-mesh.LegalTrig(el1); } bad2 /= hasonepi.Size(); @@ -791,6 +796,142 @@ namespace netgen mesh.SetNextTimeStamp(); } + void MeshOptimize2d :: SplitImprove (Mesh & mesh) + { + if (!faceindex) + { + PrintMessage (3, "Split improve"); + + mesh.CalcSurfacesOfNode(); // TODO: needed? + for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++) + { + SplitImprove (mesh); + + if (multithread.terminate) + throw NgException ("Meshing stopped"); + } + + faceindex = 0; + mesh.Compress(); // TODO: needed? + return; + } + + Array elements; + mesh.GetSurfaceElementsOfFace (faceindex, elements); + + // return if we have quads in this surface + for (auto & ei : elements) + if (mesh[ei].GetNP() != 3) + return; + + // maps from edges to adjacent trigs + INDEX_2_HASHTABLE> els_on_edge(2*elements.Size() + 2); + + // build els_on_edge table + for (SurfaceElementIndex sei : elements) + { + const Element2d & sel = mesh[sei]; + + for (int j = 0; j < 3; j++) + { + PointIndex pi1 = sel.PNumMod(j+2); + PointIndex pi2 = sel.PNumMod(j+3); + + if (mesh.IsSegment (pi1, pi2)) continue; + + INDEX_2 ii2 (pi1, pi2); + ii2.Sort(); + if (els_on_edge.Used (ii2)) + { + auto els = els_on_edge.Get(ii2); + get<1>(els) = sei; + els_on_edge.Set(ii2, els); + } + else + { + els_on_edge.Set (ii2, make_tuple(sei, sei)); + } + } + } + + // split edges of illegal trigs + for (SurfaceElementIndex sei : elements) + { + Element2d & sel = mesh[sei]; + + if (sel.IsDeleted()) continue; + + // TODO: split also bad trigs, nut just illegal ones + if (mesh.LegalTrig(sel)) continue; + + for (int j = 0; j < 3; j++) + { + PointIndex pi1 = sel.PNumMod(j+2); + PointIndex pi2 = sel.PNumMod(j+3); + PointIndex pi3 = sel.PNumMod(j+1); + PointIndex pi4; + PointGeomInfo gi1 = sel.GeomInfoPiMod(j+2); + PointGeomInfo gi2 = sel.GeomInfoPiMod(j+3); + PointGeomInfo gi3 = sel.GeomInfoPiMod(j+1); + PointGeomInfo gi4; + + if (mesh.IsSegment (pi1, pi2)) continue; + + // get neighbor element + INDEX_2 ii2 (pi1, pi2); + ii2.Sort(); + auto els = els_on_edge.Get(ii2); + SurfaceElementIndex other_i = get<0>(els); + if(other_i==sei) other_i = get<1>(els); + auto & other = mesh[other_i]; + + // find opposite point of neighbor element + for (int j = 0; j < 3; j++) + if(other[j]!=pi1 && other[j]!=pi2) + { + pi4 = other[j]; + gi4 = other.GeomInfoPi(j); + break; + } + + // split edge pi1,pi2 + Point<3> p5; + PointIndex pi5; + PointGeomInfo gi5; + + mesh.GetGeometry()->GetRefinement().PointBetween (mesh[pi1], mesh[pi2], 0.5, + faceindex, + gi1, gi2, p5, gi5); + + pi5 = mesh.AddPoint(p5); + + Element2d e1(3); + e1.SetIndex(faceindex); + e1={ {pi1,gi1}, {pi5,gi5}, {pi3,gi3} }; + mesh.AddSurfaceElement( e1 ); + + Element2d e2(3); + e2.SetIndex(faceindex); + e2 ={ {pi5,gi5}, {pi2,gi2}, {pi3,gi3} }; + mesh.AddSurfaceElement( e2 ); + + Element2d e3(3); + e3.SetIndex(faceindex); + e3 ={ {pi1,gi1}, {pi4,gi4}, {pi5,gi5} }; + mesh.AddSurfaceElement( e3 ); + + Element2d e4(3); + e4.SetIndex(faceindex); + e4 ={ {pi4,gi4}, {pi2,gi2}, {pi5,gi5} }; + mesh.AddSurfaceElement( e4 ); + + sel.Delete(); + other.Delete(); + } + } + + mesh.SetNextTimeStamp(); + } void MeshOptimize2d :: CheckMeshApproximation (Mesh & mesh) { diff --git a/libsrc/meshing/improve2.hpp b/libsrc/meshing/improve2.hpp index e0e0c70a..c9fbae4b 100644 --- a/libsrc/meshing/improve2.hpp +++ b/libsrc/meshing/improve2.hpp @@ -24,6 +24,7 @@ public: void EdgeSwapping (Mesh & mesh, int usemetric); void CombineImprove (Mesh & mesh); + void SplitImprove (Mesh & mesh); void GenericImprove (Mesh & mesh); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 64b7f2ba..dd89e962 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -3692,6 +3692,30 @@ namespace netgen bool Mesh :: LegalTrig (const Element2d & el) const { + // Search for surface trigs with same vertices ( may happen for instance with close surfaces in stl geometies ) + if(!illegal_trigs) + { + auto & tab = const_cast(illegal_trigs); + tab = make_unique> (3*GetNSE() + 1); + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + { + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; + + INDEX_3 i3(sel[0], sel[1], sel[2]); + i3.Sort(); + if(tab->Used(i3) && tab->Get(i3)!=sei) + tab -> Set (i3, -1); + else + tab -> Set (i3, sei); + } + + } + INDEX_3 i3 (el[0], el[1], el[2]); + i3.Sort(); + if(illegal_trigs->Used(i3) && illegal_trigs->Get(i3)==-1) + return false; + return 1; if ( /* hp */ 1) // needed for old, simple hp-refinement { diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 4b91fab3..d310df67 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -55,6 +55,7 @@ namespace netgen unique_ptr> segmentht; /// unique_ptr> surfelementht; + unique_ptr> illegal_trigs; /// faces of rest-solid NgArray openelements; From 8defe2f8645406cb0f697ea78894af28a615c02a Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 30 Sep 2019 14:19:12 +0200 Subject: [PATCH 2/2] New function Mesh::FindIllegalTrigs() Find illegal trigs after surface meshing and use this information in surface mesh optimization --- libsrc/csg/genmesh.cpp | 2 ++ libsrc/meshing/meshclass.cpp | 56 +++++++++++++++++++++---------- libsrc/meshing/meshclass.hpp | 4 +++ libsrc/occ/occgenmesh.cpp | 3 ++ libsrc/stlgeom/meshstlsurface.cpp | 3 ++ 5 files changed, 50 insertions(+), 18 deletions(-) diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index b00ecc34..3bc1c440 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -504,6 +504,8 @@ namespace netgen for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++) mesh[sei].SetIndex (k); + auto n_illegal_trigs = mesh.FindIllegalTrigs(); + PrintMessage (3, n_illegal_trigs, " illegal triangles"); // mesh.CalcSurfacesOfNode(); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index dd89e962..23157576 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -3689,31 +3689,51 @@ namespace netgen return 0; } + int Mesh :: FindIllegalTrigs () + { + // Temporary table to store the vertex numbers of all triangles + INDEX_3_CLOSED_HASHTABLE temp_tab(3*GetNSE() + 1); + size_t cnt = 0; + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + { + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; + + INDEX_3 i3(sel[0], sel[1], sel[2]); + i3.Sort(); + if(temp_tab.Used(i3)) + { + temp_tab.Set (i3, -1); + cnt++; + } + else + { + temp_tab.Set (i3, sei); + } + } + + illegal_trigs = make_unique> (2*cnt+1); + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) + { + const Element2d & sel = surfelements[sei]; + if (sel.IsDeleted()) continue; + + INDEX_3 i3(sel[0], sel[1], sel[2]); + i3.Sort(); + if(temp_tab.Get(i3)==-1) + illegal_trigs -> Set (i3, 1); + } + return cnt; + } bool Mesh :: LegalTrig (const Element2d & el) const { // Search for surface trigs with same vertices ( may happen for instance with close surfaces in stl geometies ) if(!illegal_trigs) - { - auto & tab = const_cast(illegal_trigs); - tab = make_unique> (3*GetNSE() + 1); - for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) - { - const Element2d & sel = surfelements[sei]; - if (sel.IsDeleted()) continue; - - INDEX_3 i3(sel[0], sel[1], sel[2]); - i3.Sort(); - if(tab->Used(i3) && tab->Get(i3)!=sei) - tab -> Set (i3, -1); - else - tab -> Set (i3, sei); - } - - } + throw Exception("In Mesh::LegalTrig() - illegal_trigs table not built"); INDEX_3 i3 (el[0], el[1], el[2]); i3.Sort(); - if(illegal_trigs->Used(i3) && illegal_trigs->Get(i3)==-1) + if(illegal_trigs->Used(i3)) return false; return 1; diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index d310df67..2e548f28 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -571,6 +571,10 @@ namespace netgen /// + // Find trigs with same vertices + // return: number of illegal trigs + int FindIllegalTrigs (); + bool LegalTrig (const Element2d & el) const; /** if values non-null, return values in 4-double array: diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index 414b52f6..78c72cb2 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -862,6 +862,9 @@ namespace netgen for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++) mesh[sei].SetIndex (k); + + auto n_illegal_trigs = mesh.FindIllegalTrigs(); + PrintMessage (3, n_illegal_trigs, " illegal triangles"); } // ofstream problemfile("occmesh.rep"); diff --git a/libsrc/stlgeom/meshstlsurface.cpp b/libsrc/stlgeom/meshstlsurface.cpp index 4eadc1e0..87ab7f39 100644 --- a/libsrc/stlgeom/meshstlsurface.cpp +++ b/libsrc/stlgeom/meshstlsurface.cpp @@ -279,6 +279,9 @@ int STLSurfaceMeshing (STLGeometry & geom, class Mesh & mesh, const MeshingParam mesh.FindOpenSegments(); nopen = mesh.GetNOpenSegments(); + auto n_illegal_trigs = mesh.FindIllegalTrigs(); + PrintMessage (3, n_illegal_trigs, " illegal triangles"); + if (nopen) { geom.ClearMarkedSegs();