fix periodic identifications for meshes where edges touch

This commit is contained in:
Christopher Lackner 2024-08-28 11:07:11 +02:00
parent 69f2ea5635
commit 1497bf36cc

View File

@ -8,25 +8,26 @@ namespace netgen
{ {
struct PointTree struct PointTree
{ {
BoxTree<3> tree; std::map<int, BoxTree<3>> tree;
Box<3> bounding_box;
PointTree( Box<3> bb ) : tree(bb) {} PointTree( Box<3> bb ) : bounding_box(bb) {}
void Insert(Point<3> p, PointIndex n) void Insert(Point<3> p, PointIndex n, int index)
{ {
tree.Insert(p, p, n); if(tree.count(index) == 0)
tree.emplace(index, bounding_box);
tree.at(index).Insert(p, p, n);
} }
PointIndex Find(Point<3> p) const PointIndex Find(Point<3> p, int index) const
{ {
ArrayMem<int, 1> points; ArrayMem<int, 1> points;
tree.GetIntersecting(p, p, points); tree.at(index).GetIntersecting(p, p, points);
if(points.Size()==0) if(points.Size()==0)
throw Exception("cannot find mapped point " + ToString(p)); throw Exception("cannot find mapped point " + ToString(p));
return points[0]; return points[0];
} }
double GetTolerance() { return tree.GetTolerance(); }
}; };
DLL_HEADER GeometryRegisterArray geometryregister; DLL_HEADER GeometryRegisterArray geometryregister;
@ -584,7 +585,6 @@ namespace netgen
for(auto & vert : vertices) for(auto & vert : vertices)
{ {
auto pi = mesh.AddPoint(vert->GetPoint(), vert->properties.layer); auto pi = mesh.AddPoint(vert->GetPoint(), vert->properties.layer);
tree.Insert(mesh[pi], pi);
vert2meshpt[vert->nr] = pi; vert2meshpt[vert->nr] = pi;
mesh[pi].Singularity(vert->properties.hpref); mesh[pi].Singularity(vert->properties.hpref);
mesh[pi].SetType(FIXEDPOINT); mesh[pi].SetType(FIXEDPOINT);
@ -716,7 +716,8 @@ namespace netgen
for(auto i : Range(edge_points)) for(auto i : Range(edge_points))
{ {
auto pi = mesh.AddPoint(edge_points[i], edge->properties.layer); auto pi = mesh.AddPoint(edge_points[i], edge->properties.layer);
tree.Insert(mesh[pi], pi); if(edge->identifications.Size())
tree.Insert(mesh[pi], pi, edge->nr);
pnums[i+1] = pi; pnums[i+1] = pi;
} }
@ -757,7 +758,7 @@ namespace netgen
p_other = (*ident.trafo)(mesh[pi]); p_other = (*ident.trafo)(mesh[pi]);
else else
static_cast<GeometryEdge*>(ident.to)->ProjectPoint(p_other, nullptr); static_cast<GeometryEdge*>(ident.to)->ProjectPoint(p_other, nullptr);
auto pi_other = tree.Find(p_other); auto pi_other = tree.Find(p_other, ident.to->nr);
identifications.Add(pi, pi_other, ident.name, ident.type); identifications.Add(pi, pi_other, ident.name, ident.type);
} }
} }
@ -866,7 +867,7 @@ namespace netgen
for(auto pi : s.PNums()) for(auto pi : s.PNums())
if(!is_point_in_tree[pi]) if(!is_point_in_tree[pi])
{ {
tree.Insert(mesh[pi], pi); tree.Insert(mesh[pi], pi, -1);
is_point_in_tree[pi] = true; is_point_in_tree[pi] = true;
} }
@ -936,8 +937,8 @@ namespace netgen
edges[mapped_edges[edgenr]]->ProjectPoint(p2, nullptr); edges[mapped_edges[edgenr]]->ProjectPoint(p2, nullptr);
edges[mapped_edges[edgenr]]->ProjectPoint(p3, nullptr); edges[mapped_edges[edgenr]]->ProjectPoint(p3, nullptr);
} }
sel[2] = tree.Find(p2); sel[2] = tree.Find(p2, -1);
sel[3] = tree.Find(p3); sel[3] = tree.Find(p3, -1);
// find mapped segment to set PointGeomInfo correctly // find mapped segment to set PointGeomInfo correctly
Segment s_other; Segment s_other;
@ -1017,7 +1018,7 @@ namespace netgen
if(mesh[pi].Type() == SURFACEPOINT && pi_to_face[pi]==-1) if(mesh[pi].Type() == SURFACEPOINT && pi_to_face[pi]==-1)
{ {
pi_to_face[pi] = face->nr; pi_to_face[pi] = face->nr;
tree.Insert(mesh[pi], pi); tree.Insert(mesh[pi], pi, -1);
pi_of_face[face->nr].Append(pi); pi_of_face[face->nr].Append(pi);
} }
} }
@ -1085,7 +1086,7 @@ namespace netgen
else else
edge->ProjectPoint(p, nullptr); edge->ProjectPoint(p, nullptr);
} }
tree.Insert(p, pi); tree.Insert(p, pi, -1);
is_point_in_tree[pi] = true; is_point_in_tree[pi] = true;
} }
} }
@ -1097,7 +1098,7 @@ namespace netgen
{ {
auto pi = seg[i]; auto pi = seg[i];
if(!pmap[pi].IsValid()) if(!pmap[pi].IsValid())
pmap[tree.Find(mesh[pi])] = pi; pmap[tree.Find(mesh[pi], -1)] = pi;
// store uv values (might be different values for same point in case of internal edges) // store uv values (might be different values for same point in case of internal edges)
double u = seg.epgeominfo[i].u; double u = seg.epgeominfo[i].u;