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

further work on boundarylayers

Browse Source 
better calculation of growthvector, fix bug with addsegment
This commit is contained in:
Christopher Lackner 2020-07-03 19:51:06 +02:00
parent 7da5cfd3de
commit fdd718739f
1 changed files with 17 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
libsrc/meshing/boundarylayer.cpp
View File

@ -247,14 +247,13 @@ namespace netgen
if(growthvectors[pi].Length2() == 0.) if(growthvectors[pi].Length2() == 0.)
continue; continue;
auto& g = growthvectors[pi]; auto& g = growthvectors[pi];
auto gn = g * n; auto ng = n * g;
auto gg = g * g; auto gg = g * g;
auto nn = n * n; auto nn = n * n;
auto l2 = -2*gn/(gn*gn/gg + nn); if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
auto l1 = l2 * gn/gg; auto a = -ng*ng/(ng*ng-nn * gg);
auto new_g = g + 0.5 * (l1 * g + l2 * n); auto b = ng*gg/(ng*ng-nn*gg);
if(new_g * g > 0) g += a*g + b*n;
g = new_g;
} }
} }
@ -296,11 +295,11 @@ namespace netgen
// Set them all to "1" to initially activate all segments // Set them all to "1" to initially activate all segments
segsel.Set(); segsel.Set();
Array<Array<SegmentIndex>> segmap(nseg);
// remove double segments (if multiple surfaces come together // remove double segments (if more than 2 surfaces come together
// in one edge. If one of them is mapped, keep that one and // in one edge. If one of them is mapped, keep that one and
// map the others to it. // map the others to it.
Array<Array<SegmentIndex>> segmap(nseg);
for(SegmentIndex sei = 0; sei < nseg; sei++) for(SegmentIndex sei = 0; sei < nseg; sei++)
{ {
if(!segsel.Test(sei)) continue; if(!segsel.Test(sei)) continue;
@ -320,6 +319,7 @@ namespace netgen
segmap[main].Append(s); segmap[main].Append(s);
segmap[other].SetSize(0); segmap[other].SetSize(0);
segmap[main].Append(other); segmap[main].Append(other);
if(other == sei) sej = nseg;
} }
} }
} }
@ -333,7 +333,9 @@ namespace netgen
// surface index is part of the "hit-list" // surface index is part of the "hit-list"
if(segsel.Test(sei)) if(segsel.Test(sei))
{ {
auto& segi = mesh[sei]; // copy here since we will add segments and this would
// invalidate a reference!
auto segi = mesh[sei];
if(blp.surfid.Contains(segi.si)) if(blp.surfid.Contains(segi.si))
{ {
// clear the bit to indicate that this segment has been processed // clear the bit to indicate that this segment has been processed
@ -342,7 +344,9 @@ namespace netgen
// Find matching segment pair on other surface // Find matching segment pair on other surface
for(SegmentIndex sej = 0; sej < nseg; sej++) for(SegmentIndex sej = 0; sej < nseg; sej++)
{ {
auto& segj = mesh[sej]; // copy here since we will add segments and this would
// invalidate a reference!
auto segj = mesh[sej];
// Find the segment pair on the neighbouring surface element // Find the segment pair on the neighbouring surface element
// Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0] // Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0]
if(segsel.Test(sej) && ((segi[0] == segj[1]) && (segi[1] == segj[0]))) if(segsel.Test(sej) && ((segi[0] == segj[1]) && (segi[1] == segj[0])))
@ -456,7 +460,9 @@ namespace netgen
mesh.AddSegment(s1); mesh.AddSegment(s1);
mesh.AddSegment(s2); mesh.AddSegment(s2);
} }
// segi[0] = mapto[segi[0]] not working somehow? // do not use segi (not even with reference, since
// mesh.AddSegment will resize segment array and
// invalidate reference), this is why we copy it!!!
mesh[sei][0] = mapto[segi[0]]; mesh[sei][0] = mapto[segi[0]];
mesh[sei][1] = mapto[segi[1]]; mesh[sei][1] = mapto[segi[1]];
mesh[sej][0] = mapto[segj[0]]; mesh[sej][0] = mapto[segj[0]];