ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-19 01:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

some ting

Browse Source
This commit is contained in:
Matthias Hochsteger 2023-10-23 11:40:37 +02:00
parent 35f084e9aa
commit 57be10cbcf
1 changed files with 35 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
libsrc/meshing/boundarylayer.cpp
View File

@ -14,7 +14,7 @@ namespace netgen
{ {
// checks if a segment is intersecting a plane, spanned by three points, lam will be set s.t. p_intersect = seg[0] + lam * (seg[1]-seg[0]) // checks if a segment is intersecting a plane, spanned by three points, lam will be set s.t. p_intersect = seg[0] + lam * (seg[1]-seg[0])
// bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const array<Point<3>, 3> & trig, const array<double, 3> trig_lams, double & lam, double & lam2) // bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const array<Point<3>, 3> & trig, const array<double, 3> trig_lams, double & lam, double & lam2)
bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const array<Point<3>, 3> & trig, double & lam, array<double, 2> & lam2) bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const array<Point<3>, 3> & trig, array<double, 3> & lam)
{ {
auto t1 = trig[1]-trig[0]; auto t1 = trig[1]-trig[0];
auto t2 = trig[2]-trig[0]; auto t2 = trig[2]-trig[0];
@ -24,24 +24,26 @@ namespace netgen
if(v0n * v1n >= 0) if(v0n * v1n >= 0)
return false; return false;
lam = -v0n/(v1n-v0n); lam[0] = -v0n/(v1n-v0n);
auto p = seg[0] + lam*(seg[1]-seg[0]); auto p = seg[0] + lam[0]*(seg[1]-seg[0]);
lam *= 0.9; lam[0] *= 0.9;
if(lam < -1e-8 || lam>1+1e-8) if(lam[0] < -1e-8 || lam[0]>1+1e-8)
return false; return false;
const auto area = 0.5 * Cross(t1,t2).Length(); const auto area = 0.5 * Cross(t1,t2).Length();
const auto area0 = 0.5 * Cross(t1,{p}).Length(); const auto area0 = 0.5 * Cross(t1,{p}).Length();
const auto area1 = 0.5 * Cross(t2,{p}).Length(); const auto area1 = 0.5 * Cross(t2,{p}).Length();
lam2[0] = area0/area; lam[1] = area0/area;
lam2[1] = area1/area; lam[2] = area1/area;
return true; return true;
} }
bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const array<Point<3>, 3> & trig, double & lam ) bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const array<Point<3>, 3> & trig, double & lam )
{ {
array<double,2> lam2; array<double,3> vlam;
return isIntersectingPlane(seg, trig, lam, lam2); auto result = isIntersectingPlane(seg, trig, vlam);
lam = vlam[0];
return result;
} }
bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const ArrayMem<Point<3>, 4> & face, double & lam) bool isIntersectingPlane ( const array<Point<3>, 2> & seg, const ArrayMem<Point<3>, 4> & face, double & lam)
@ -174,6 +176,30 @@ namespace netgen
void BoundaryLayerTool :: LimitGrowthVector(PointIndex pi, SurfaceElementIndex sei) void BoundaryLayerTool :: LimitGrowthVector(PointIndex pi, SurfaceElementIndex sei)
{ {
auto sel = mesh[sei];
auto seg = GetMappedSeg(pi);
struct Face {
ArrayMem<Point<3>, 4> p;
ArrayMem<double, 4> lam;
};
ArrayMem<Face, 6> faces;
faces.Append({GetMappedFace(sei, -2), {0.0, 0.0, 0.0, 0.0}});
faces.Append({GetMappedFace(sei, -1), {1.0, 1.0, 1.0, 1.0}});
for(auto i : Range(sel.GetNP())) {
Face face;
face.p = GetMappedFace(sei, i);
face.lam = {0.0, 0.0, 1.0, 1.0};
faces.Append(face);
}
ArrayMem<tuple<double,double>, 6> intersections;
for(auto & face : faces) {
array<double, 3> lam;
if(isIntersectingPlane(seg, face.p, face.lam))
intersections.Append({face.lam[0], face.lam[1]});
}
} }