Merge branch 'fix_point_type_of_vertices' into 'master'

fix point type of geo vertices (FIXEDPOINT) -> locked points

See merge request jschoeberl/netgen!486
This commit is contained in:
Joachim Schöberl 2022-03-02 10:47:37 +00:00
commit 3918990b0e
7 changed files with 174 additions and 94 deletions

View File

@ -529,14 +529,12 @@ namespace netgen
tree.Insert(mesh[pi], pi); tree.Insert(mesh[pi], pi);
vert2meshpt[vert->GetHash()] = pi; vert2meshpt[vert->GetHash()] = pi;
mesh[pi].Singularity(vert->properties.hpref); mesh[pi].Singularity(vert->properties.hpref);
mesh[pi].SetType(FIXEDPOINT);
if(vert->properties.name) Element0d el(pi, pi);
{ el.name = vert->properties.GetName();
Element0d el(pi, pi); mesh.SetCD3Name(pi, el.name);
el.name = vert->properties.GetName(); mesh.pointelements.Append (el);
mesh.SetCD3Name(pi, el.name);
mesh.pointelements.Append (el);
}
} }
for(auto & vert : vertices) for(auto & vert : vertices)

View File

@ -180,12 +180,12 @@ namespace netgen
} }
} }
void InterpolateSurfaceGrowthVectors(const Mesh & mesh, const BoundaryLayerParameters& blp, int fd_old, FlatArray<Vec<3>, PointIndex> growthvectors) void InterpolateSurfaceGrowthVectors(const Mesh & mesh, const BoundaryLayerParameters& blp, int fd_old, FlatArray<Vec<3>, PointIndex> growthvectors, const Table<SurfaceElementIndex, PointIndex> & p2sel)
{ {
auto np = mesh.GetNP();
// interpolate growth vectors at inner surface points from surrounding edge points // interpolate growth vectors at inner surface points from surrounding edge points
Array<Point<2>, PointIndex> delaunay_points(mesh.GetNP()); Array<Point<2>, PointIndex> delaunay_points(np);
Array<int, PointIndex> p2face(mesh.GetNP()); Array<PointIndex, PointIndex> pmap(np); // maps duplicated points
p2face = 0;
Array<SurfaceElementIndex> surface_els; Array<SurfaceElementIndex> surface_els;
Array<PointIndex> edge_points; Array<PointIndex> edge_points;
@ -195,33 +195,75 @@ namespace netgen
if(!blp.surfid.Contains(facei)) if(!blp.surfid.Contains(facei))
continue; continue;
p2face = 0;
edge_points.SetSize(0); edge_points.SetSize(0);
surface_points.SetSize(0); surface_points.SetSize(0);
surface_els.SetSize(0); surface_els.SetSize(0);
delaunay_points.SetSize(np);
pmap.SetSize(np);
mesh.GetSurfaceElementsOfFace (facei, surface_els); mesh.GetSurfaceElementsOfFace (facei, surface_els);
Box<2> bbox ( Box<2>::EMPTY_BOX ); Box<2> bbox ( Box<2>::EMPTY_BOX );
for(auto sei : surface_els) for(auto sei : surface_els)
{ {
const auto & sel = mesh[sei]; auto sel = mesh[sei];
for (auto i : Range(sel.GetNP())) for (auto i : Range(sel.GetNP()))
{ {
auto pi = sel[i]; auto & gi = sel.GeomInfo()[i];
if(p2face[pi] != 0) Point<2> p = {gi.u, gi.v};
bbox.Add(p);
}
}
BoxTree<2> tree(bbox);
for(auto pi : mesh.Points().Range())
{
auto n_surf_els = p2sel[pi].Size();
bool has_relevant_sel = false;
for(auto sei : p2sel[pi])
if(mesh[sei].GetIndex()==facei)
{
has_relevant_sel = true;
break;
}
if(!has_relevant_sel)
continue;
if(mesh[pi].Type() <= EDGEPOINT)
edge_points.Append(pi);
else
surface_points.Append(pi);
// the same point might have different uv coordinates (closed edges for instance)
// duplicate these points for the delaunay tree
bool inserted = false;
for(auto sei : p2sel[pi])
{
auto sel = mesh[sei];
if(sel.GetIndex()!=facei)
continue; continue;
p2face[pi] = facei; PointGeomInfo gi = sel.GeomInfo()[sel.PNums().Pos(pi)];
Point<2> p = {gi.u, gi.v};
bool found = false;
tree.GetFirstIntersecting( p, p, [&] (const auto pi_found) { return found = true; });
if(found)
continue;
if(mesh[pi].Type() <= EDGEPOINT) auto pi_new = pi;
edge_points.Append(pi); if(inserted)
{
pi_new = delaunay_points.Append(p);
pmap.Append(pi);
edge_points.Append(pi_new);
}
else else
surface_points.Append(pi); {
delaunay_points[pi] = p;
auto & gi = sel.GeomInfo()[i]; pmap[pi] = pi;
// TODO: project to plane if u,v not available? }
delaunay_points[pi] = {gi.u, gi.v}; tree.Insert(p, pi_new);
bbox.Add(delaunay_points[pi]); inserted = true;
} }
} }
@ -231,23 +273,24 @@ namespace netgen
DelaunayMesh dmesh( delaunay_points, bbox ); DelaunayMesh dmesh( delaunay_points, bbox );
for(auto pi : edge_points) for(auto pi : edge_points)
{
p2face[pi] = 0;
dmesh.AddPoint(pi); dmesh.AddPoint(pi);
}
std::map<PointIndex, double> weights; std::map<PointIndex, double> weights;
for(auto pi : surface_points) for(auto pi : surface_points)
{ {
dmesh.AddPoint(pi, &weights); dmesh.CalcIntersecting(pi);
dmesh.CalcWeights(pi, weights);
auto & v = growthvectors[pi]; auto & v = growthvectors[pi];
auto n = 1./v.Length() * v;
for(auto & [pi_other, weight] : weights) for(auto & [pi_other, weight] : weights)
v += weight * growthvectors[pi_other]; {
// interpolate only the tangential part of the growth vector
auto t = weight * growthvectors[pmap[pi_other]];
t -= (t * n) * n;
v += t;
}
} }
} }
} }
void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp) void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp)
@ -271,6 +314,8 @@ namespace netgen
if(!blp.outside) if(!blp.outside)
domains.Invert(); domains.Invert();
auto& meshtopo = mesh.GetTopology();
meshtopo.SetBuildVertex2Element(true);
mesh.UpdateTopology(); mesh.UpdateTopology();
bool have_single_segments = HaveSingleSegments(mesh); bool have_single_segments = HaveSingleSegments(mesh);
@ -280,8 +325,6 @@ namespace netgen
else else
segments = mesh.LineSegments(); segments = mesh.LineSegments();
auto& meshtopo = mesh.GetTopology();
int np = mesh.GetNP(); int np = mesh.GetNP();
int ne = mesh.GetNE(); int ne = mesh.GetNE();
int nse = mesh.GetNSE(); int nse = mesh.GetNSE();
@ -329,26 +372,50 @@ namespace netgen
} }
} }
// mark points for remapping const auto & p2sel = mesh.CreatePoint2SurfaceElementTable();
for(const auto& sel : mesh.SurfaceElements())
{ for(auto pi : mesh.Points().Range())
auto n = surfacefacs[sel.GetIndex()] * getSurfaceNormal(sel); {
if(n.Length2() != 0.) const auto & p = mesh[pi];
{ if(p.Type() == INNERPOINT)
for(auto pi : sel.PNums()) continue;
{
if(interpolate_growth_vectors && mesh[pi].Type() >= SURFACEPOINT) std::map<int, Vec<3>> normals;
continue;
auto & np = growthvectors[pi]; // calculate one normal vector per face (average with angles as weights for multiple surface elements within a face)
if(np.Length() == 0) { np = n; continue; } for(auto sei : p2sel[pi])
auto npn = np * n; {
auto npnp = np * np; const auto & sel = mesh[sei];
auto nn = n * n; auto facei = sel.GetIndex();
if(nn-npn*npn/npnp == 0) { np = n; continue; } if(!blp.surfid.Contains(facei))
np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np); continue;
}
} auto n = surfacefacs[sel.GetIndex()] * getSurfaceNormal(sel);
}
int itrig = sel.PNums().Pos(pi);
itrig += sel.GetNP();
auto v0 = (mesh[sel.PNumMod(itrig+1)] - mesh[pi]).Normalize();
auto v1 = (mesh[sel.PNumMod(itrig-1)] - mesh[pi]).Normalize();
if(normals.count(facei)==0)
normals[facei] = {0.,0.,0.};
normals[facei] += acos(v0*v1)*n;
}
for(auto & [facei, n] : normals)
n *= 1.0/n.Length();
// combine normal vectors for each face to keep uniform distances
auto & np = growthvectors[pi];
for(auto & [facei, n] : normals)
{
if(np.Length() == 0) { np = n; continue; }
auto npn = np * n;
auto npnp = np * np;
auto nn = n * n;
if(nn-npn*npn/npnp == 0) { np = n; continue; }
np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np);
}
}
// Bit array to keep track of segments already processed // Bit array to keep track of segments already processed
BitArray segs_done(nseg+1); BitArray segs_done(nseg+1);
@ -432,8 +499,6 @@ namespace netgen
for(auto i : Range(sel.PNums())) for(auto i : Range(sel.PNums()))
{ {
auto pi = sel.PNums()[i]; auto pi = sel.PNums()[i];
if(interpolate_growth_vectors && mesh[pi].Type() >= SURFACEPOINT)
continue;
if(growthvectors[pi].Length2() == 0.) if(growthvectors[pi].Length2() == 0.)
continue; continue;
auto next = sel.PNums()[(i+1)%sel.GetNV()]; auto next = sel.PNums()[(i+1)%sel.GetNV()];
@ -492,25 +557,29 @@ namespace netgen
if(seg.edgenr-1 == edgenr && mesh[seg[0]].Type() == FIXEDPOINT) if(seg.edgenr-1 == edgenr && mesh[seg[0]].Type() == FIXEDPOINT)
{ {
points.Append(seg[0]); points.Append(seg[0]);
points.Append(seg[1]);
break; break;
} }
} }
while(true) while(true)
{ {
bool point_found = false;
for(auto si : meshtopo.GetVertexSegments(points.Last())) for(auto si : meshtopo.GetVertexSegments(points.Last()))
{ {
const auto& seg = mesh[si]; const auto& seg = mesh[si];
if(seg.edgenr-1 != edgenr) if(seg.edgenr-1 != edgenr)
continue; continue;
if(seg[0] == points.Last() && if(seg[0] == points.Last() && points[points.Size()-2] !=seg[1])
(points.Size() < 2 || points[points.Size()-2] !=seg[1]))
{ {
points.Append(seg[1]); points.Append(seg[1]);
point_found = true;
break; break;
} }
} }
if(mesh[points.Last()].Type() == FIXEDPOINT) if(mesh[points.Last()].Type() == FIXEDPOINT)
break; break;
if(!point_found)
throw Exception(string("Could not find connected list of line segements for edge ") + edgenr);
} }
// tangential part of growth vectors // tangential part of growth vectors
@ -541,7 +610,7 @@ namespace netgen
} }
if(interpolate_growth_vectors) if(interpolate_growth_vectors)
InterpolateSurfaceGrowthVectors(mesh, blp, fd_old, growthvectors); InterpolateSurfaceGrowthVectors(mesh, blp, fd_old, growthvectors, p2sel);
// insert new points // insert new points
for (PointIndex pi = 1; pi <= np; pi++) for (PointIndex pi = 1; pi <= np; pi++)
@ -879,6 +948,7 @@ namespace netgen
} }
mesh.GetTopology().ClearEdges(); mesh.GetTopology().ClearEdges();
mesh.UpdateTopology(); mesh.UpdateTopology();
mesh.SetGeometry(nullptr);
} }
void AddDirection( Vec<3> & a, Vec<3> b ) void AddDirection( Vec<3> & a, Vec<3> b )

View File

@ -190,26 +190,7 @@ namespace netgen
if(definitive_overlapping_trig==-1) if(definitive_overlapping_trig==-1)
{ {
Mesh m; // GetMesh(pi_new)->Save("error.vol.gz");
m.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
for(auto pi : points.Range())
m.AddPoint(P3(points[pi]));
for (DelaunayTrig & trig : trigs)
{
if (trig[0] < 0) continue;
Vec<3> n = Cross (P3(points[trig[1]])-P3(points[trig[0]]),
P3(points[trig[2]])-P3(points[trig[0]]));
if (n(2) < 0) Swap (trig[1], trig[2]);
Element2d el(trig[0], trig[1], trig[2]);
el.SetIndex (1);
m.AddSurfaceElement (el);
}
m.Compress();
m.AddPoint(P3(points[pi_new]));
m.Save("error.vol.gz");
throw Exception("point not in any circle "+ ToString(pi_new)); throw Exception("point not in any circle "+ ToString(pi_new));
} }
@ -299,15 +280,18 @@ namespace netgen
auto pi_last = *points.Range().end()-3; auto pi_last = *points.Range().end()-3;
for(auto edge : edges) for(auto edge : edges)
{ {
auto v0 = points[edge[0]] - p;
auto v1 = points[edge[1]] - p;
v0.Normalize();
v1.Normalize();
double angle = acos(v0*v1);
for(PointIndex pi : {edge[0], edge[1]}) for(PointIndex pi : {edge[0], edge[1]})
{ {
if(pi>=pi_last) if(pi>=pi_last)
continue; continue;
if(weights.count(pi)) double weight = angle/(eps+Dist(p, points[pi]));
continue; sum += weight;
double weight = 1.0/(eps+Dist(p, points[pi])); weights[pi] += weight;
sum += weight;
weights[pi] = weight;
} }
} }
double isum = 1.0/sum; double isum = 1.0/sum;
@ -315,15 +299,12 @@ namespace netgen
weight *= isum; weight *= isum;
} }
void DelaunayMesh::AddPoint( PointIndex pi_new, std::map<PointIndex, double> * weights ) void DelaunayMesh::AddPoint( PointIndex pi_new)
{ {
static Timer t("AddPoint"); RegionTimer reg(t); static Timer t("AddPoint"); RegionTimer reg(t);
CalcIntersecting(pi_new); CalcIntersecting(pi_new);
if(weights)
CalcWeights(pi_new, *weights);
for (int j : intersecting) for (int j : intersecting)
{ {
UnsetNeighbours(j); UnsetNeighbours(j);
@ -339,6 +320,31 @@ namespace netgen
tree->DeleteElement (j); tree->DeleteElement (j);
} }
unique_ptr<Mesh> DelaunayMesh::GetMesh(PointIndex pi_new)
{
auto mesh = make_unique<Mesh>();
Mesh & m = *mesh;
m.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
for(auto pi : points.Range())
m.AddPoint(P3(points[pi]));
for (DelaunayTrig & trig : trigs)
{
if (trig[0] < 0) continue;
Vec<3> n = Cross (P3(points[trig[1]])-P3(points[trig[0]]),
P3(points[trig[2]])-P3(points[trig[0]]));
if (n(2) < 0) Swap (trig[1], trig[2]);
Element2d el(trig[0], trig[1], trig[2]);
el.SetIndex (1);
m.AddSurfaceElement (el);
}
m.Compress();
m.AddPoint(P3(points[pi_new]));
return mesh;
}
ostream & operator<< (ostream & ost, DelaunayTrig trig) ostream & operator<< (ostream & ost, DelaunayTrig trig)
{ {
ost << trig[0] << "-" << trig[1] << "-" << trig[2] << endl; ost << trig[0] << "-" << trig[1] << "-" << trig[2] << endl;

View File

@ -66,9 +66,9 @@ namespace netgen
void CalcIntersecting( PointIndex pi_new ); void CalcIntersecting( PointIndex pi_new );
void CalcWeights( PointIndex pi_new, std::map<PointIndex, double> & weights ); void CalcWeights( PointIndex pi_new, std::map<PointIndex, double> & weights );
void AddPoint( PointIndex pi_new, std::map<PointIndex, double> * weights = nullptr ); void AddPoint( PointIndex pi_new );
Array<DelaunayTrig> & GetElements() { return trigs; } Array<DelaunayTrig> & GetElements() { return trigs; }
unique_ptr<Mesh> GetMesh(PointIndex pi_new); // for debugging purposes
}; };
} // namespace netgen } // namespace netgen

View File

@ -2352,6 +2352,8 @@ namespace netgen
for (int i = 0; i < lockedpoints.Size(); i++) for (int i = 0; i < lockedpoints.Size(); i++)
points[lockedpoints[i]].SetType(FIXEDPOINT); points[lockedpoints[i]].SetType(FIXEDPOINT);
for(const auto& pointel : pointelements)
points[pointel.pnum].SetType(FIXEDPOINT);
/* /*
for (i = 0; i < openelements.Size(); i++) for (i = 0; i < openelements.Size(); i++)

View File

@ -102,7 +102,7 @@ namespace netgen
// mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points? // mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points?
for(auto pi : mesh.LockedPoints()) for(auto pi : mesh.LockedPoints())
for(auto i : Range(ret)) for(auto i : Range(ret))
ipmap[i][pi] = 1; ipmap[i][pi] = 2;
// add used points to domain mesh, build point mapping // add used points to domain mesh, build point mapping
for(auto i : Range(ret)) for(auto i : Range(ret))
@ -113,7 +113,10 @@ namespace netgen
for(auto pi : Range(ipmap[i])) for(auto pi : Range(ipmap[i]))
if(ipmap[i][pi]) if(ipmap[i][pi])
{ {
auto pi_new = m.AddPoint( mesh[pi] ); const auto& mp = mesh[pi];
auto pi_new = m.AddPoint( mp, mp.GetLayer(), mp.Type() );
if(ipmap[i][pi] == 2)
mesh.AddLockedPoint(pi_new);
ipmap[i][pi] = pi_new; ipmap[i][pi] = pi_new;
pmap.Append( pi ); pmap.Append( pi );
} }

View File

@ -186,6 +186,7 @@ public:
{ return vert2element[vnr]; } { return vert2element[vnr]; }
void GetVertexSurfaceElements( int vnr, Array<SurfaceElementIndex>& elements ) const; void GetVertexSurfaceElements( int vnr, Array<SurfaceElementIndex>& elements ) const;
const auto & GetVertexSurfaceElements( ) const { return vert2surfelement; }
FlatArray<SurfaceElementIndex> GetVertexSurfaceElements(PointIndex vnr) const FlatArray<SurfaceElementIndex> GetVertexSurfaceElements(PointIndex vnr) const
{ return vert2surfelement[vnr]; } { return vert2surfelement[vnr]; }