diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 2eb7a203..10544c39 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -86,6 +86,9 @@ namespace netgen protected: GeometryVertex *start, *end; public: + // Neighboring domains in 2d + // In 3d unused, EXCEPT for free floating edges in a domain, + // then both are pointing to the containing domain int domin=-1, domout=-1; GeometryEdge( GeometryVertex &start_, GeometryVertex &end_ ) @@ -187,7 +190,10 @@ namespace netgen }; class DLL_HEADER GeometrySolid : public GeometryShape - { }; + { + public: + Array free_edges; // edges with no adjacent face + }; class DLL_HEADER NetgenGeometry { diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index d1cb3c0f..a6540cf2 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -564,7 +564,7 @@ namespace netgen - void Delaunay1 (Mesh & mesh, const MeshingParameters & mp, AdFront3 * adfront, + void Delaunay1 (Mesh & mesh, int domainnr, const MeshingParameters & mp, AdFront3 * adfront, NgArray & tempels, int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax) { @@ -623,6 +623,13 @@ namespace netgen for (PointIndex pi : mesh.LockedPoints()) usep[pi] = true; + // mark points of free edge segments (no adjacent face) + for (auto & seg : mesh.LineSegments()) + if(seg.domin == domainnr && seg.domout == domainnr) + { + usep[seg[0]] = true; + usep[seg[1]] = true; + } NgArray freelist; @@ -1554,7 +1561,7 @@ namespace netgen int np = mesh.GetNP(); - Delaunay1 (mesh, mp, adfront, tempels, oldnp, startel, pmin, pmax); + Delaunay1 (mesh, domainnr, mp, adfront, tempels, oldnp, startel, pmin, pmax); { // improve delaunay - mesh by swapping !!!! diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index 573bdeba..7e323dfd 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -14,6 +14,11 @@ namespace netgen { +static constexpr int tetedges[6][2] = + { { 0, 1 }, { 0, 2 }, { 0, 3 }, + { 1, 2 }, { 1, 3 }, { 2, 3 } }; + + static constexpr int IMPROVEMENT_CONFORMING_EDGE = -1e6; static inline bool NotTooBad(double bad1, double bad2) @@ -258,6 +263,26 @@ double MeshOptimize3d :: CombineImproveEdge ( for (auto ei : has_both_points) badness_old += mesh[ei].GetBadness(); + if (goal == OPT_CONFORM && p0.Type() <= EDGEPOINT) { + // check if the optimization improves conformity with free segments + std::set edges_before, edges_after; + + for (auto ei : has_one_point) { + const auto el = mesh[ei]; + for(auto i : Range(6)) { + auto e0 = el[tetedges[i][0]]; + auto e1 = el[tetedges[i][1]]; + if(e0 == pi0 || e1 == pi0) edges_before.insert(e0 == pi0 ? e1 : e0); + if(e0 == pi1 || e1 == pi1) edges_after.insert(e0 == pi1 ? e1 : e0); + } + } + + for(auto new_edge : edges_after) { + if (edges_before.count(new_edge) == 0 && mesh[new_edge].Type() <= EDGEPOINT && mesh.BoundaryEdge (new_edge, pi0)) + badness_old += GetLegalPenalty(); + } + } + MeshPoint pnew = p0; if (p0.Type() == INNERPOINT) pnew = Center (p0, p1); @@ -1548,10 +1573,6 @@ void MeshOptimize3d :: SwapImproveSurface ( // loop over edges - static const int tetedges[6][2] = - { { 0, 1 }, { 0, 2 }, { 0, 3 }, - { 1, 2 }, { 1, 3 }, { 2, 3 } }; - pi1 = elemi[tetedges[j][0]]; pi2 = elemi[tetedges[j][1]]; @@ -2406,6 +2427,8 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face, !mesh.LegalTet(elem2)) bad1 += GetLegalPenalty(); + if(mesh.BoundaryEdge (pi4, pi5)) + bad1 += GetLegalPenalty(); el31.PNum(1) = pi1; el31.PNum(2) = pi2; @@ -2583,10 +2606,6 @@ double MeshOptimize3d :: SplitImprove2Element ( return false; // search for very flat tets, with two disjoint edges nearly crossing, like a rectangle with diagonals - static constexpr int tetedges[6][2] = - { { 0, 1 }, { 0, 2 }, { 0, 3 }, - { 1, 2 }, { 1, 3 }, { 2, 3 } }; - int minedge = -1; double mindist = 1e99; double minlam0, minlam1; diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index b47ce478..050c00e9 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -80,6 +80,16 @@ namespace netgen m.AddFaceDescriptor( mesh.GetFaceDescriptor(i) ); } + // mark interior edge points + for(const auto& seg : mesh.LineSegments()) + { + if(seg.domin > 0 && seg.domin == seg.domout) + { + ipmap[seg.domin-1][seg[0]] = 1; + ipmap[seg.domin-1][seg[1]] = 1; + } + } + // mark used points for each domain, add surface elements (with wrong point numbers) to domain mesh for(const auto & sel : mesh.SurfaceElements()) { @@ -522,6 +532,7 @@ namespace netgen } } RemoveIllegalElements (mesh, domain); + ConformToFreeSegments (mesh, domain); } void MergeMeshes( Mesh & mesh, Array & md ) @@ -615,18 +626,24 @@ namespace netgen { ParallelFor( md.Range(), [&](int i) { - if (mp.checkoverlappingboundary) - if (md[i].mesh->CheckOverlappingBoundary()) - { - if(debugparam.write_mesh_on_error) - md[i].mesh->Save("overlapping_mesh_domain_"+ToString(md[i].domain)+".vol.gz"); - throw NgException ("Stop meshing since boundary mesh is overlapping"); - } + try { + if (mp.checkoverlappingboundary) + if (md[i].mesh->CheckOverlappingBoundary()) + { + if(debugparam.write_mesh_on_error) + md[i].mesh->Save("overlapping_mesh_domain_"+ToString(md[i].domain)+".vol.gz"); + throw NgException ("Stop meshing since boundary mesh is overlapping"); + } - if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC) - FillCloseSurface( md[i] ); - CloseOpenQuads( md[i] ); - MeshDomain(md[i]); + if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC) + FillCloseSurface( md[i] ); + CloseOpenQuads( md[i] ); + MeshDomain(md[i]); + } + catch (const Exception & e) { + cerr << "Meshing of domain " << i+1 << " failed with error: " << e.what() << endl; + throw e; + } }, md.Size()); } catch(...) @@ -734,6 +751,65 @@ namespace netgen } + void ConformToFreeSegments (Mesh & mesh, int domain) + { + auto geo = mesh.GetGeometry(); + if(!geo) return; + auto n_solids = geo->GetNSolids(); + if(!n_solids) return; + if(geo->GetSolid(domain-1).free_edges.Size() == 0) + return; + + Segment bad_seg; + Array free_segs; + for (auto seg : mesh.LineSegments()) + if(seg.domin == domain && seg.domout == domain) + free_segs.Append(seg); + + auto num_nonconforming = [&] () { + size_t count = 0; + + auto p2el = mesh.CreatePoint2ElementTable(); + + for (auto seg : free_segs) { + + auto has_p0 = p2el[seg[0]]; + bool has_both = false; + + for(auto ei : has_p0) { + if(mesh[ei].PNums().Contains(seg[1])) + has_both = true; + } + + if(!has_both) { + bad_seg = seg; + count++; + } + } + return count; + }; + + for (auto i : Range(5)) { + auto num_bad_segs = num_nonconforming(); + PrintMessage(1, "Non-conforming free segments in domain ", domain, ": ", num_bad_segs); + + if(num_bad_segs == 0) + return; + + MeshingParameters dummymp; + MeshOptimize3d optmesh(mesh, dummymp, OPT_CONFORM); + + for (auto i : Range(3)) { + optmesh.SwapImprove2 (); + optmesh.SwapImprove(); + optmesh.CombineImprove(); + } + } + + if(debugparam.write_mesh_on_error) + mesh.Save("free_segment_not_conformed_dom_"+ToString(domain)+"_seg_"+ToString(bad_seg[0])+"_"+ToString(bad_seg[1])+".vol.gz"); + throw Exception("Segment not resolved in volume mesh in domain " + ToString(domain)+ ", seg: " + ToString(bad_seg)); + } void RemoveIllegalElements (Mesh & mesh3d, int domain) diff --git a/libsrc/meshing/meshfunc.hpp b/libsrc/meshing/meshfunc.hpp index 854cdd8f..8d603212 100644 --- a/libsrc/meshing/meshfunc.hpp +++ b/libsrc/meshing/meshfunc.hpp @@ -31,6 +31,7 @@ DLL_HEADER MESHING3_RESULT OptimizeVolume (const MeshingParameters & mp, Mesh& m // const CSGeometry * geometry = NULL); DLL_HEADER void RemoveIllegalElements (Mesh & mesh3d, int domain = 0); +DLL_HEADER void ConformToFreeSegments (Mesh & mesh3d, int domain); enum MESHING_STEP { diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 13156783..033dd645 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1132,6 +1132,21 @@ namespace netgen } } */ + + std::map> free_edges_in_solid; + for(auto i1 : Range(1, somap.Extent()+1)) + { + auto s = somap(i1); + for (auto edge : MyExplorer(s, TopAbs_EDGE, TopAbs_WIRE)) + if (!emap.Contains(edge)) + { + free_edges_in_solid[i1].Append(emap.Add (edge)); + for (auto vertex : MyExplorer(edge, TopAbs_VERTEX)) + if (!vmap.Contains(vertex)) + vmap.Add (vertex); + } + } + for (auto edge : MyExplorer(shape, TopAbs_EDGE, TopAbs_WIRE)) if (!emap.Contains(edge)) { @@ -1256,6 +1271,16 @@ namespace netgen if(face.Shape().Orientation() == TopAbs_INTERNAL) face.domout = k; } + + if(free_edges_in_solid.count(i1)) + for(auto ei : free_edges_in_solid[i1]) + { + auto & edge = GetEdge(emap(ei)); + edge.properties.maxh = min(edge.properties.maxh, occ_solid->properties.maxh); + edge.domin = k; + edge.domout = k; + occ_solid->free_edges.Append(&GetEdge(emap(ei))); + } solids.Append(std::move(occ_solid)); }