From 2ad56cd7efe47b5d54f941abda76dec82abc5c1a Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 20 Jun 2023 12:44:18 +0200 Subject: [PATCH] Add edge/face midpoints to bounding box in element search tree --- libsrc/meshing/meshclass.cpp | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index a5595b23..53f92aa5 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -5221,7 +5221,37 @@ namespace netgen auto & el = volelements[ei]; if(el.IsCurved() && curvedelems->IsElementCurved(ei)) + { + // add edge/face midpoints to box + auto eltype = el.GetType(); + const auto verts = topology.GetVertices(eltype); + + + const auto edges = FlatArray(topology.GetNEdges(eltype), topology.GetEdges0(eltype)); + for (const auto & edge: edges) { + netgen::Point<3> lam = netgen::Point<3>(0.5* (verts[edge[0]] + verts[edge[1]])); + auto p = netgen::Point<3>(0.0); + curvedelems->CalcElementTransformation(lam,ei,p); + box.Add(p); + } + + const auto faces = FlatArray(topology.GetNFaces(eltype), topology.GetFaces0(eltype)); + for (const auto & face: faces) { + netgen::Vec<3> lam = netgen::Vec<3>(verts[face[0]] + verts[face[1] + verts[face[2]]]); + if(face[3] != -1) { + lam += netgen::Vec<3>(verts[face[3]]); + lam *= 0.25; + } + else + lam *= 1.0/3; + auto p = netgen::Point<3>(0.0); + curvedelems->CalcElementTransformation(netgen::Point<3>(lam),ei,p); + box.Add(p); + } + + box.Scale(1.2); + } elementsearchtree -> Insert (box, ei+1);