ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-14 18:08:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

BLayer - Option to build only volelements without new surface elements/edges/domains

Browse Source
This commit is contained in:
Matthias Hochsteger 2024-10-03 18:06:34 +02:00
parent 2d653b2672
commit 0daeeb20aa
13 changed files with 235 additions and 388 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

387
libsrc/meshing/boundarylayer.cpp
View File

@ -145,7 +145,7 @@ Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr,
cout << pts << endl;
for (auto *p_seg : segs)
cout << *p_seg << endl;
throw Exception("Something went wrong in getEdgeTangent!");
throw NG_EXCEPTION("Something went wrong in getEdgeTangent!");
}
tangent = mesh[pts[1]] - mesh[pts[0]];
return tangent.Normalize();
@ -253,36 +253,8 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh &mesh_,
static Timer timer("BoundaryLayerTool::ctor");
RegionTimer regt(timer);
ProcessParameters();
// for(auto & seg : mesh.LineSegments())
// seg.edgenr = seg.epgeominfo[1].edgenr;
total_height = 0.0;
for (auto h : par_heights)
total_height += h;
max_edge_nr = -1;
for (const auto &seg : mesh.LineSegments())
if (seg.edgenr > max_edge_nr)
max_edge_nr = seg.edgenr;
int ndom = mesh.GetNDomains();
ndom_old = ndom;
new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1);
new_mat_nrs = -1;
for (auto [bcname, matname] : par_new_mat) {
mesh.SetMaterial(++ndom, matname);
regex pattern(bcname);
for (auto i : Range(1, mesh.GetNFD() + 1)) {
auto &fd = mesh.GetFaceDescriptor(i);
if (regex_match(fd.GetBCName(), pattern))
new_mat_nrs[i] = ndom;
}
}
if (!params.outside)
domains.Invert();
if (domains.NumSet() == 0)
return;
topo.SetBuildVertex2Element(true);
mesh.UpdateTopology();
@ -321,15 +293,17 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() {
isIn != domains.Test(fd.DomainOut())) {
int new_si = mesh.GetNFD() + 1;
surfacefacs[i] = isIn ? 1. : -1.;
// -1 surf nr is so that curving does not do anything
FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(),
isIn ? fd.DomainOut() : new_mat_nrs[i], -1);
new_fd.SetBCProperty(new_si);
new_fd.SetSurfColour(fd.SurfColour());
mesh.AddFaceDescriptor(new_fd);
si_map[i] = new_si;
moved_surfaces.SetBit(i);
mesh.SetBCName(new_si - 1, "mapped_" + name);
if (!insert_only_volume_elements) {
// -1 surf nr is so that curving does not do anything
FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(),
isIn ? fd.DomainOut() : new_mat_nrs[i], -1);
new_fd.SetBCProperty(new_si);
new_fd.SetSurfColour(fd.SurfColour());
mesh.AddFaceDescriptor(new_fd);
si_map[i] = new_si;
mesh.SetBCName(new_si - 1, "mapped_" + name);
}
}
// curving of surfaces with boundary layers will often
// result in pushed through elements, since we do not (yet)
@ -347,6 +321,8 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() {
}
void BoundaryLayerTool ::CreateFaceDescriptorsSides() {
if (insert_only_volume_elements)
return;
BitArray face_done(mesh.GetNFD() + 1);
face_done.Clear();
for (const auto &sel : mesh.SurfaceElements()) {
@ -717,7 +693,8 @@ void BoundaryLayerTool ::InsertNewElements(
sel.GeomInfo()[i].v = 0.0;
}
sel.SetIndex(si_map[segj.si]);
mesh.AddSurfaceElement(sel);
new_sels.Append(sel);
new_sels_on_moved_bnd.Append(sel);
// TODO: Too many, would be enough to only add outermost ones
Segment s1;
@ -771,8 +748,6 @@ void BoundaryLayerTool ::InsertNewElements(
BitArray fixed_points(np + 1);
fixed_points.Clear();
BitArray moveboundarypoint(np + 1);
moveboundarypoint.Clear();
auto p2el = mesh.CreatePoint2ElementTable();
for (SurfaceElementIndex si = 0; si < nse; si++) {
// copy because surfaceels array will be resized!
@ -820,44 +795,7 @@ void BoundaryLayerTool ::InsertNewElements(
for (auto i : Range(points))
newel[i] = newPoint(sel[i], -1, groups[i]);
newel.SetIndex(si_map[sel.GetIndex()]);
mesh.AddSurfaceElement(newel);
// also move volume element adjacent to this surface element accordingly
ElementIndex ei = -1;
// if(groups[0] || groups[1] || groups[2])
// for(auto ei_ : p2el[sel.PNums()[0]])
// {
// const auto & el = mesh[ei_];
// // if(!domains.Test(el.GetIndex())) continue;
// cout << "check " << ei_ << "\t" << el << "\t" << sel << endl;
// auto pnums = el.PNums();
// if(pnums.Contains(sel[1]) && pnums.Contains(sel[2])) {
// ei = ei_;
// break;
// }
// }
if (ei != -1) {
auto &el = mesh[ei];
for (auto i : Range(el.GetNP()))
for (auto j : Range(3)) {
if (groups[j] && el[i] == sel[j]) {
el[i] = newel[j];
break;
}
}
}
} else {
bool has_moved = false;
for (auto p : sel.PNums())
has_moved |= hasMoved(p);
if (has_moved)
for (auto p : sel.PNums()) {
if (hasMoved(p)) {
fixed_points.SetBit(p);
if (is_boundary_moved.Test(sel.GetIndex()))
moveboundarypoint.SetBit(p);
}
}
new_sels.Append(newel);
}
if (is_boundary_moved.Test(sel.GetIndex())) {
for (auto &p : mesh[si].PNums())
@ -866,25 +804,27 @@ void BoundaryLayerTool ::InsertNewElements(
}
}
for (SegmentIndex sei = 0; sei < nseg; sei++) {
auto &seg = segments[sei];
if (is_boundary_moved.Test(seg.si))
for (auto &p : seg.PNums())
if (hasMoved(p))
p = newPoint(p);
// else if(hasMoved(seg[0]) || hasMoved(seg[1]))
// {
// auto tangent = mesh[seg[1]] - mesh[seg[0]];
// if(hasMoved(seg[0]) && growthvectors[seg[0]] * tangent > 0)
// seg[0] = newPoint(seg[0]);
// if(hasMoved(seg[1]) && growthvectors[seg[1]] * tangent < 0)
// seg[1] = newPoint(seg[1]);
// }
}
// for (SegmentIndex sei = 0; sei < nseg; sei++) {
// auto &seg = segments[sei];
// if (is_boundary_moved.Test(seg.si)) {
// if(insert_only_volume_elements) throw
// Exception("insert_only_volume_elements and is_boundary_moved are
// incompatible"); for (auto &p : seg.PNums())
// if (hasMoved(p))
// p = newPoint(p);
// }
// }
// fill holes in surface mesh at special boundary points (i.e. points with >=4
// adjacent boundary faces)
auto p2sel = mesh.CreatePoint2SurfaceElementTable();
auto p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>(
new_sels.Range(),
[&](auto &table, SurfaceElementIndex ei) {
for (PointIndex pi : new_sels[ei].PNums())
table.Add(pi, ei);
},
mesh.GetNP());
for (auto &[special_pi, special_point] : special_boundary_points) {
if (special_point.growth_groups.Size() != 2)
throw Exception("special_point.growth_groups.Size() != 2");
@ -921,7 +861,7 @@ void BoundaryLayerTool ::InsertNewElements(
for (auto &pi : sel.PNums())
if (pi != pi_common && pi != new_special_pi0)
pi = new_special_pi1;
mesh.AddSurfaceElement(sel);
new_sels.Append(sel);
}
}
}
@ -960,155 +900,17 @@ void BoundaryLayerTool ::InsertNewElements(
sel[1] = seg[0];
sel[2] = pi_new_other;
sel.SetIndex(face);
mesh.AddSurfaceElement(sel);
new_sels.Append(sel);
}
}
}
}
}
for (ElementIndex ei = 0; ei < ne; ei++) {
auto el = mesh[ei];
ArrayMem<PointIndex, 4> fixed;
ArrayMem<PointIndex, 4> moved;
bool moved_bnd = false;
for (const auto &p : el.PNums()) {
if (fixed_points.Test(p))
fixed.Append(p);
if (hasMoved(p))
moved.Append(p);
if (moveboundarypoint.Test(p))
moved_bnd = true;
}
bool do_move, do_insert;
if (changed_domains.Test(el.GetIndex())) {
do_move = fixed.Size() && moved_bnd;
do_insert = do_move;
} else {
do_move = !fixed.Size() || moved_bnd;
do_insert = !do_move;
}
// if (do_move) {
// for (auto& p : mesh[ei].PNums())
// if (hasMoved(p)) {
// if (special_boundary_points.count(p)) {
// auto& special_point = special_boundary_points[p];
// auto& group = special_point.growth_groups[0];
// p = group.new_points.Last();
// } else
// p = newPoint(p);
// }
// }
if (do_insert) {
if (el.GetType() == TET) {
if (moved.Size() == 3) // inner corner
{
PointIndex p1 = moved[0];
PointIndex p2 = moved[1];
PointIndex p3 = moved[2];
auto v1 = mesh[p1];
auto n = Cross(mesh[p2] - v1, mesh[p3] - v1);
auto d = mesh[newPoint(p1, 0)] - v1;
if (n * d > 0)
Swap(p2, p3);
PointIndex p4 = p1;
PointIndex p5 = p2;
PointIndex p6 = p3;
for (auto i : Range(par_heights)) {
Element nel(PRISM);
nel[0] = p4;
nel[1] = p5;
nel[2] = p6;
p4 = newPoint(p1, i);
p5 = newPoint(p2, i);
p6 = newPoint(p3, i);
nel[3] = p4;
nel[4] = p5;
nel[5] = p6;
nel.SetIndex(el.GetIndex());
mesh.AddVolumeElement(nel);
}
}
if (moved.Size() == 2) {
if (fixed.Size() == 1) {
PointIndex p1 = moved[0];
PointIndex p2 = moved[1];
for (auto i : Range(par_heights)) {
PointIndex p3 = newPoint(moved[1], i);
PointIndex p4 = newPoint(moved[0], i);
Element nel(PYRAMID);
nel[0] = p1;
nel[1] = p2;
nel[2] = p3;
nel[3] = p4;
nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] -
moved[1];
if (Cross(mesh[p2] - mesh[p1], mesh[p4] - mesh[p1]) *
(mesh[nel[4]] - mesh[nel[1]]) >
0)
Swap(nel[1], nel[3]);
nel.SetIndex(el.GetIndex());
mesh.AddVolumeElement(nel);
p1 = p4;
p2 = p3;
}
}
}
if (moved.Size() == 1 && fixed.Size() == 1) {
PointIndex p1 = moved[0];
for (auto i : Range(par_heights)) {
Element nel = el;
PointIndex p2 = newPoint(moved[0], i);
for (auto &p : nel.PNums()) {
if (p == moved[0])
p = p1;
else if (p == fixed[0])
p = p2;
}
p1 = p2;
mesh.AddVolumeElement(nel);
}
}
} else if (el.GetType() == PYRAMID) {
if (moved.Size() == 2) {
if (fixed.Size() != 2)
throw Exception("This case is not implemented yet! Fixed size = " +
ToString(fixed.Size()));
PointIndex p1 = moved[0];
PointIndex p2 = moved[1];
for (auto i : Range(par_heights)) {
PointIndex p3 = newPoint(moved[1], i);
PointIndex p4 = newPoint(moved[0], i);
Element nel(PYRAMID);
nel[0] = p1;
nel[1] = p2;
nel[2] = p3;
nel[3] = p4;
nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] -
fixed[1] - moved[0] - moved[1];
if (Cross(mesh[p2] - mesh[p1], mesh[p4] - mesh[p1]) *
(mesh[nel[4]] - mesh[nel[1]]) >
0)
Swap(nel[1], nel[3]);
nel.SetIndex(el.GetIndex());
mesh.AddVolumeElement(nel);
p1 = p4;
p2 = p3;
}
} else if (moved.Size() == 1)
throw Exception("This case is not implemented yet!");
} else if (do_move) {
throw Exception(
"Boundarylayer only implemented for tets and pyramids outside "
"yet!");
}
}
}
}
void BoundaryLayerTool ::SetDomInOut() {
if (insert_only_volume_elements)
return;
for (auto i : Range(1, nfd_old + 1))
if (moved_surfaces.Test(i)) {
if (auto dom = mesh.GetFaceDescriptor(si_map[i]).DomainIn();
@ -1121,6 +923,8 @@ void BoundaryLayerTool ::SetDomInOut() {
}
void BoundaryLayerTool ::SetDomInOutSides() {
if (insert_only_volume_elements)
return;
BitArray done(mesh.GetNFD() + 1);
done.Clear();
for (auto sei : Range(mesh.SurfaceElements())) {
@ -1146,6 +950,8 @@ void BoundaryLayerTool ::SetDomInOutSides() {
}
void BoundaryLayerTool ::AddSegments() {
if (insert_only_volume_elements)
return;
if (have_single_segments)
MergeAndAddSegments(mesh, segments, new_segments);
else {
@ -1155,45 +961,11 @@ void BoundaryLayerTool ::AddSegments() {
}
}
void BoundaryLayerTool ::FixVolumeElements() {
static Timer timer("BoundaryLayerTool::FixVolumeElements");
RegionTimer rt(timer);
BitArray is_inner_point(mesh.GetNP() + 1);
is_inner_point.Clear();
auto changed_domains = domains;
if (!params.outside)
changed_domains.Invert();
for (ElementIndex ei : Range(ne))
if (changed_domains.Test(mesh[ei].GetIndex()))
for (auto pi : mesh[ei].PNums())
if (mesh[pi].Type() == INNERPOINT)
is_inner_point.SetBit(pi);
Array<PointIndex> points;
for (auto pi : mesh.Points().Range())
if (is_inner_point.Test(pi))
points.Append(pi);
auto p2el = mesh.CreatePoint2ElementTable(is_inner_point);
// smooth growth vectors to shift additional element layers to the inside and
// fix flipped tets
for ([[maybe_unused]] auto step : Range(0)) {
for (auto pi : points) {
Vec<3> average_gw = 0.0;
auto &els = p2el[pi];
size_t cnt = 0;
for (auto ei : els)
if (ei < ne)
for (auto pi1 : mesh[ei].PNums())
if (pi1 <= np) {
average_gw += growthvectors[pi1];
cnt++;
}
growthvectors[pi] = 1.0 / cnt * average_gw;
}
void BoundaryLayerTool ::AddSurfaceElements() {
for (auto &sel :
insert_only_volume_elements ? new_sels_on_moved_bnd : new_sels) {
cout << "add surface element " << sel << endl;
mesh.AddSurfaceElement(sel);
}
}
@ -1242,10 +1014,14 @@ void BoundaryLayerTool ::ProcessParameters() {
for (auto id : surfids)
par_surfid.Append(id);
}
if (string *mat = get_if<string>(&params.new_material); mat)
par_new_mat = {{".*", *mat}};
else
par_new_mat = *get_if<map<string, string>>(&params.new_material);
insert_only_volume_elements = !params.new_material.has_value();
if (params.new_material) {
if (string *mat = get_if<string>(&*params.new_material); mat)
par_new_mat = {{".*", *mat}};
else
par_new_mat = *get_if<map<string, string>>(&*params.new_material);
}
if (params.project_boundaries.has_value()) {
auto proj_bnd = *params.project_boundaries;
@ -1282,6 +1058,53 @@ void BoundaryLayerTool ::ProcessParameters() {
for (auto i : *get_if<std::vector<int>>(&params.domain))
domains.SetBit(i);
}
if (domains.NumSet() == 0)
return;
total_height = 0.0;
for (auto h : par_heights)
total_height += h;
max_edge_nr = -1;
for (const auto &seg : mesh.LineSegments())
if (seg.edgenr > max_edge_nr)
max_edge_nr = seg.edgenr;
int ndom = mesh.GetNDomains();
ndom_old = ndom;
new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1);
new_mat_nrs = -1;
if (insert_only_volume_elements) {
for (auto i : Range(1, mesh.GetNFD() + 1)) {
auto &fd = mesh.GetFaceDescriptor(i);
auto domin = fd.DomainIn();
auto domout = fd.DomainOut();
for (int dom : {domin, domout})
if (domains.Test(dom)) {
if (params.outside) {
dom = domin + domout - dom;
if (dom == 0)
throw NG_EXCEPTION("No new material specified for boundarylayer "
"on the outside of domain");
}
new_mat_nrs[i] = dom;
}
}
} else {
for (auto [bcname, matname] : par_new_mat) {
mesh.SetMaterial(++ndom, matname);
regex pattern(bcname);
for (auto i : Range(1, mesh.GetNFD() + 1)) {
auto &fd = mesh.GetFaceDescriptor(i);
if (regex_match(fd.GetBCName(), pattern))
new_mat_nrs[i] = ndom;
}
}
}
cout << "new mat numbers " << endl << new_mat_nrs << endl;
if (!params.outside)
domains.Invert();
}
void BoundaryLayerTool ::Perform() {
@ -1296,6 +1119,7 @@ void BoundaryLayerTool ::Perform() {
SetDomInOut();
AddSegments();
AddSurfaceElements();
mesh.CalcSurfacesOfNode();
topo.SetBuildVertex2Element(true);
@ -1311,6 +1135,7 @@ void BoundaryLayerTool ::Perform() {
mesh[pi] += height * (*gw);
}
mesh.CalcSurfacesOfNode();
mesh.GetTopology().ClearEdges();
mesh.SetNextMajorTimeStamp();
mesh.UpdateTopology();

5
libsrc/meshing/boundarylayer.hpp
View File

@ -49,6 +49,7 @@ class BoundaryLayerTool
MeshTopology & topo;
BoundaryLayerParameters params;
Array<Vec<3>, PointIndex> growthvectors;
std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
Table<SurfaceElementIndex, PointIndex> p2sel;
BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved;
@ -62,11 +63,13 @@ class BoundaryLayerTool
// These parameters are derived from given BoundaryLayerParameters and the Mesh
Array<double> par_heights;
Array<int> par_surfid;
bool insert_only_volume_elements;
map<string, string> par_new_mat;
Array<size_t> par_project_boundaries;
bool have_single_segments;
Array<Segment> segments, new_segments;
Array<Element2d, SurfaceElementIndex> new_sels, new_sels_on_moved_bnd;
Array<Array<PointIndex>, PointIndex> mapto;
Array<PointIndex, PointIndex> mapfrom;
@ -91,7 +94,7 @@ class BoundaryLayerTool
void SetDomInOut();
void SetDomInOutSides();
void AddSegments();
void FixVolumeElements();
void AddSurfaceElements();
Vec<3> getNormal(const Element2d & el)
{

2
libsrc/meshing/boundarylayer2d.cpp
View File

@ -1,5 +1,5 @@
#include <mystdlib.h>
#include "meshing.hpp"
#include "boundarylayer.hpp"
#include "meshing2.hpp"
#include "../geom2d/csg2d.hpp"

10
libsrc/meshing/boundarylayer_interpolate.cpp
View File

@ -2,10 +2,6 @@
namespace netgen {
// TODO: Hack, move this to the header or restructure the whole growth_vectors
// storage
static std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
void BoundaryLayerTool ::InterpolateGrowthVectors() {
int new_max_edge_nr = max_edge_nr;
for (const auto &seg : segments)
@ -113,10 +109,10 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
if (plast != seg[0] && plast != seg[1])
continue;
auto pnew = plast == seg[0] ? seg[1] : seg[0];
if(pnew == points[0] && points.Size()>1) {
if (pnew == points[0] && points.Size() > 1) {
}
if (points_set.count(pnew) > 0 && (pnew != points[0] || points.Size()==2))
if (points_set.count(pnew) > 0 &&
(pnew != points[0] || points.Size() == 2))
continue;
edge_len += (mesh[points.Last()] - mesh[pnew]).Length();
points.Append(pnew);

158
libsrc/meshing/boundarylayer_limiter.hpp
View File

@ -25,17 +25,31 @@ struct GrowthVectorLimiter {
GrowthVectorLimiter(BoundaryLayerTool &tool_)
: tool(tool_), params(tool_.params), mesh(tool_.mesh),
height(tool_.total_height), growthvectors(tool_.growthvectors),
map_from(mesh.Points().Size()),
p2sel(mesh.CreatePoint2SurfaceElementTable()) {
map_from(mesh.Points().Size()) {
changed_domains = tool.domains;
if (!params.outside)
changed_domains.Invert();
map_from = tool.mapfrom;
p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>(
tool.new_sels.Range(),
[&](auto &table, SurfaceElementIndex ei) {
for (PointIndex pi : tool.new_sels[ei].PNums())
table.Add(pi, ei);
},
mesh.GetNP());
}
auto SurfaceElementsRange() { return Range(tool.nse + tool.new_sels.Size()); }
const auto &Get(SurfaceElementIndex sei) {
if (sei < tool.nse)
return mesh[sei];
return tool.new_sels[sei - tool.nse];
}
std::pair<double, double> GetMinMaxLimit(SurfaceElementIndex sei) {
const auto &sel = mesh[sei];
const auto &sel = Get(sei);
double min_limit = GetLimit(sel[0]);
double max_limit = min_limit;
for (auto i : IntRange(1, sel.GetNP())) {
@ -75,7 +89,7 @@ struct GrowthVectorLimiter {
Point<3> GetPoint(PointIndex pi_to, double shift = 1.,
bool apply_limit = false) {
if (tool.growth_vector_map.count(pi_to) == 0)
if (pi_to <= tool.np || tool.growth_vector_map.count(pi_to) == 0)
return mesh[pi_to];
return mesh[pi_to] + GetVector(pi_to, shift, apply_limit);
@ -97,7 +111,7 @@ struct GrowthVectorLimiter {
auto GetTrig(SurfaceElementIndex sei, double shift = 0.0,
bool apply_limit = false) {
auto sel = mesh[sei];
auto sel = Get(sei);
std::array<Point<3>, 3> trig;
for (auto i : Range(3))
trig[i] = GetPoint(sel[i], shift, apply_limit);
@ -105,7 +119,7 @@ struct GrowthVectorLimiter {
}
auto GetMappedTrig(SurfaceElementIndex sei, double shift = 0.0) {
auto sel = mesh[sei];
auto sel = Get(sei);
std::array<Point<3>, 3> trig;
for (auto i : Range(3))
trig[i] = GetMappedPoint(sel[i], shift);
@ -115,7 +129,7 @@ struct GrowthVectorLimiter {
auto GetSideTrig(SurfaceElementIndex sei, int index, double shift = 0.0,
bool grow_first_vertex = true) {
auto trig = GetMappedTrig(sei, 0.0);
auto sel = mesh[sei];
auto sel = Get(sei);
auto index1 = (index + 1) % 3;
if (!grow_first_vertex)
index1 = (index + 2) % 3;
@ -133,7 +147,7 @@ struct GrowthVectorLimiter {
auto seg = GetSeg(pi_to, seg_shift, true);
for (auto pi : mesh[sei].PNums()) {
for (auto pi : Get(sei).PNums()) {
if (pi == pi_from)
return false;
if (map_from[pi] == pi_from)
@ -172,7 +186,7 @@ struct GrowthVectorLimiter {
double max_limit = max(GetLimit(pi_to), trig_max_limit);
bool result = false;
result = SetLimit(pi_to, s * max_limit);
for (auto pi : mesh[sei].PNums())
for (auto pi : Get(sei).PNums())
result = SetLimit(pi, s * max_limit);
return result;
} else {
@ -198,8 +212,8 @@ struct GrowthVectorLimiter {
return;
for (PointIndex pi : IntRange(tool.np, mesh.GetNP())) {
std::set<PointIndex> pis;
for (auto sel : p2sel[pi])
for (auto pi_ : mesh[sel].PNums())
for (auto sei : p2sel[pi])
for (auto pi_ : Get(sei).PNums())
pis.insert(pi_);
ArrayMem<double, 20> limits;
for (auto pi : pis) {
@ -230,7 +244,7 @@ struct GrowthVectorLimiter {
// from original surface elements
if (sei >= tool.nse)
return false;
const auto sel = mesh[sei];
const auto sel = Get(sei);
auto np = sel.GetNP();
for (auto i : Range(np)) {
if (sel[i] > tool.np)
@ -254,8 +268,6 @@ struct GrowthVectorLimiter {
for (SurfaceElementIndex sei : mesh.SurfaceElements().Range()) {
auto sel = mesh[sei];
const auto &fd = mesh.GetFaceDescriptor(sel.GetIndex());
if (sei >= tool.nse)
continue;
if (sel.GetNP() == 4)
continue;
@ -338,9 +350,9 @@ struct GrowthVectorLimiter {
tree = make_unique<BoxTree<3>>(bbox);
for (auto sei : mesh.SurfaceElements().Range()) {
const auto &sel = mesh[sei];
auto sel_index = mesh[sei].GetIndex();
for (auto sei : SurfaceElementsRange()) {
const auto &sel = Get(sei);
auto sel_index = sel.GetIndex();
Box<3> box(Box<3>::EMPTY_BOX);
for (auto pi : sel.PNums()) {
@ -369,7 +381,7 @@ struct GrowthVectorLimiter {
box.Add(GetPoint(pi_to, GetLimit(pi_from)));
tree->GetFirstIntersecting(box.PMin(), box.PMax(),
[&](SurfaceElementIndex sei) {
const auto &sel = mesh[sei];
const auto &sel = Get(sei);
if (sel.PNums().Contains(pi_from))
return false;
if (sel.PNums().Contains(pi_to))
@ -392,8 +404,8 @@ struct GrowthVectorLimiter {
mesh.GetBox(pmin, pmax);
BoxTree<3, SurfaceElementIndex> setree(pmin, pmax);
for (auto sei : mesh.SurfaceElements().Range()) {
const Element2d &tri = mesh[sei];
for (auto sei : SurfaceElementsRange()) {
const Element2d &tri = Get(sei);
Box<3> box(Box<3>::EMPTY_BOX);
for (PointIndex pi : tri.PNums())
@ -403,61 +415,71 @@ struct GrowthVectorLimiter {
setree.Insert(box, sei);
}
for (auto sei : mesh.SurfaceElements().Range()) {
const Element2d &tri = mesh[sei];
for (auto sei : SurfaceElementsRange()) {
const Element2d &tri = Get(sei);
Box<3> box(Box<3>::EMPTY_BOX);
for (PointIndex pi : tri.PNums())
box.Add(GetPoint(pi, 1.0, true));
setree.GetFirstIntersecting(
box.PMin(), box.PMax(), [&](SurfaceElementIndex sej) {
const Element2d &tri2 = mesh[sej];
setree.GetFirstIntersecting(box.PMin(), box.PMax(), [&](size_t sej) {
const Element2d &tri2 = Get(sej);
if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer())
return false;
if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer())
return false;
netgen::Point<3> tri1_points[3], tri2_points[3];
const netgen::Point<3> *trip1[3], *trip2[3];
for (int k = 0; k < 3; k++) {
trip1[k] = &tri1_points[k];
trip2[k] = &tri2_points[k];
netgen::Point<3> tri1_points[3], tri2_points[3];
const netgen::Point<3> *trip1[3], *trip2[3];
for (int k = 0; k < 3; k++) {
trip1[k] = &tri1_points[k];
trip2[k] = &tri2_points[k];
}
auto set_points = [&]() {
for (int k = 0; k < 3; k++) {
tri1_points[k] = GetPoint(tri[k], 1.0, true);
tri2_points[k] = GetPoint(tri2[k], 1.0, true);
}
};
set_points();
int counter = 0;
while (IntersectTriangleTriangle(&trip1[0], &trip2[0])) {
changed = true;
PointIndex pi_max_limit = PointIndex::INVALID;
for (PointIndex pi :
{tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]})
if (pi > tool.np && (!pi_max_limit.IsValid() ||
GetLimit(pi) > GetLimit(pi_max_limit)))
pi_max_limit = map_from[pi];
if (!pi_max_limit.IsValid())
break;
ScaleLimit(pi_max_limit, 0.9);
set_points();
counter++;
if (counter > 5) {
break;
cerr << "Limit intersecting surface elements: too many "
"limitation steps, sels: "
<< Get(sei) << '\t' << Get(sej) << endl;
for (auto si : {sei, sej}) {
auto sel = Get(si);
cerr << "Limits: ";
for (auto pi : sel.PNums())
cerr << GetLimit(pi) << ",\t";
cerr << endl;
for (auto pi : sel.PNums())
cerr << GetPoint(pi, 1.0, true) << "\t";
cerr << endl;
}
auto set_points = [&]() {
for (int k = 0; k < 3; k++) {
tri1_points[k] = GetPoint(tri[k], 1.0, true);
tri2_points[k] = GetPoint(tri2[k], 1.0, true);
}
};
set_points();
int counter = 0;
while (IntersectTriangleTriangle(&trip1[0], &trip2[0])) {
changed = true;
PointIndex pi_max_limit = PointIndex::INVALID;
for (PointIndex pi :
{tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]})
if (pi > tool.np &&
(!pi_max_limit.IsValid() ||
limits[tool.mapfrom[pi]] > limits[pi_max_limit]))
pi_max_limit = tool.mapfrom[pi];
if (!pi_max_limit.IsValid())
break;
limits[pi_max_limit] *= 0.9;
set_points();
counter++;
if (counter > 30) {
cerr << "Limit intersecting surface elements: too many "
"limitation steps, sels: "
<< mesh[sei] << '\t' << mesh[sej] << endl;
break;
}
}
return false;
});
cerr << "pi_max_limit " << pi_max_limit << endl;
break;
}
}
return false;
});
}
}
}
@ -491,7 +513,7 @@ struct GrowthVectorLimiter {
[&](PointIndex pi_to, SurfaceElementIndex sei) {
if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift))
limit_counter++;
auto sel = mesh[sei];
auto sel = Get(sei);
bool is_mapped = true;
for (auto pi : sel.PNums()) {
if (pi >= tool.np)

2
libsrc/meshing/debugging.cpp
View File

@ -36,13 +36,13 @@ namespace netgen
mesh->SetMaterial(2, "2_points");
mesh->SetMaterial(3, "3_points");
mesh->SetMaterial(4, "4_points");
mesh->Compress();
mesh->ClearSurfaceElements();
for (auto & el : openelements)
mesh->AddSurfaceElement( el );
mesh->Compress();
return mesh;
}

1
libsrc/meshing/improve3.cpp
View File

@ -2637,6 +2637,7 @@ double MeshOptimize3d :: SplitImprove2Element (
Element & elem = mesh[ei0];
if (elem.IsDeleted()) return false;
if (ei0 == ei) continue;
if (elem.GetType() != TET) return false;
if (elem[0] == pi1 || elem[1] == pi1 || elem[2] == pi1 || elem[3] == pi1 || (elem.GetNP()==5 && elem[4]==pi1) )
if(!has_both_points0.Contains(ei0))

3
libsrc/meshing/meshfunc.cpp
View File

@ -3,6 +3,7 @@
#include <mystdlib.h>
#include "meshing.hpp"
#include "debugging.hpp"
#include "boundarylayer.hpp"
namespace netgen
{
@ -421,7 +422,7 @@ namespace netgen
mesh.FindOpenElements(domain);
PrintMessage (5, mesh.GetNOpenElements(), " open faces");
// GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol");
// GetOpenElements( mesh, domain )->Save("open_"+ToString(domain)+"_"+ToString(cntsteps)+".vol");
cntsteps++;

2
libsrc/meshing/meshing.hpp
View File

@ -53,7 +53,7 @@ namespace netgen
#include "bisect.hpp"
#include "hprefinement.hpp"
#include "boundarylayer.hpp"
// #include "boundarylayer.hpp"
#include "specials.hpp"
#include "validate.hpp"
#include "basegeom.hpp"

22
libsrc/meshing/meshtype.cpp
View File

@ -2928,14 +2928,17 @@ namespace netgen
case 1: print_vec(std::get<1>(mp.thickness)); break;
}
ost <<"\n new_material: ";
switch(mp.new_material.index())
{
case 0: ost << std::get<0>(mp.new_material); break;
case 1:
for (const auto & [key, value] : std::get<1>(mp.new_material))
ost << key << " -> " << value << ", ";
break;
}
if(!mp.new_material) ost << "nullopt";
else {
switch(mp.new_material->index())
{
case 0: ost << std::get<0>(*mp.new_material); break;
case 1:
for (const auto & [key, value] : std::get<1>(*mp.new_material))
ost << key << " -> " << value << ", ";
break;
}
}
ost << "\n domain: ";
switch(mp.domain.index())
{
@ -2957,9 +2960,8 @@ namespace netgen
ost << "nullopt";
ost << "\n grow_edges: " << mp.grow_edges;
ost << "\n limit_growth_vectors: " << mp.limit_growth_vectors;
ost << "\n sides_keep_surfaceindex: " << mp.sides_keep_surfaceindex;
ost << "\n sides_keep_surfaceindex: " << (mp.sides_keep_surfaceindex ? ToString(*mp.sides_keep_surfaceindex) : "nullopt");
ost << "\n keep_surfaceindex: " << mp.keep_surfaceindex;
ost << "\n limit_safety: " << mp.limit_safety;
ost << endl;
return ost;
}

18
libsrc/meshing/meshtype.hpp
View File

@ -1276,18 +1276,16 @@ namespace netgen
struct BoundaryLayerParameters
{
std::variant<string, int, std::vector<int>> boundary;
std::variant<double, std::vector<double>> thickness;
std::variant<string, std::map<string, string>> new_material;
std::variant<string, int, std::vector<int>> domain;
bool outside;
std::optional<std::variant<string, std::vector<int>>> project_boundaries;
bool grow_edges;
bool limit_growth_vectors;
bool sides_keep_surfaceindex;
bool keep_surfaceindex;
double limit_safety = 0.3; // alloow only 30% of the growth vector length
std::variant<string, int, std::vector<int>> boundary = ".*";
std::optional<std::variant<string, std::map<string, string>>> new_material = nullopt;
std::optional<std::variant<string, std::vector<int>>> project_boundaries = nullopt;
bool outside = false;
bool grow_edges = true;
bool limit_growth_vectors = true;
std::optional<bool> sides_keep_surfaceindex = nullopt; // !outside by default
bool keep_surfaceindex = false;
};

12
libsrc/meshing/python_mesh.cpp
View File

@ -8,6 +8,7 @@
#include <mystdlib.h>
#include "meshing.hpp"
#include "boundarylayer.hpp"
// #include <csg.hpp>
// #include <geometry2d.hpp>
#include <../interface/rw_medit.hpp>
@ -1730,9 +1731,8 @@ py::arg("point_tolerance") = -1.)
std::optional<std::variant<string, std::vector<int>>> project_boundaries,
bool grow_edges,
bool limit_growth_vectors,
bool sides_keep_surfaceindex,
bool keep_surfaceindex,
double limit_safety)
std::optional<bool> sides_keep_surfaceindex,
bool keep_surfaceindex)
{
BoundaryLayerParameters blp;
blp.boundary = boundary;
@ -1745,14 +1745,13 @@ py::arg("point_tolerance") = -1.)
blp.limit_growth_vectors = limit_growth_vectors;
blp.sides_keep_surfaceindex = sides_keep_surfaceindex;
blp.keep_surfaceindex = keep_surfaceindex;
blp.limit_safety = limit_safety;
return blp;
}),
py::arg("boundary"), py::arg("thickness"), py::arg("new_material"),
py::arg("domain") = ".*", py::arg("outside") = false,
py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true,
py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false,
py::arg("keep_surfaceindex")=false, py::arg("limit_safety")=0.3,
py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=nullopt,
py::arg("keep_surfaceindex")=false,
R"delimiter(
Add boundary layer to mesh.
@ -1805,7 +1804,6 @@ project_boundaries : Optional[str] = None
.def_readwrite("limit_growth_vectors", &BoundaryLayerParameters::limit_growth_vectors)
.def_readwrite("sides_keep_surfaceindex", &BoundaryLayerParameters::sides_keep_surfaceindex)
.def_readwrite("keep_surfaceindex", &BoundaryLayerParameters::keep_surfaceindex)
.def_readwrite("limit_safety", &BoundaryLayerParameters::limit_safety)
;
py::implicitly_convertible<py::dict, BoundaryLayerParameters>();

1
ng/ngpkg.cpp
View File

@ -9,6 +9,7 @@ The interface between the GUI and the netgen library
#include <linalg.hpp>
#include <meshing.hpp>
#include "../libsrc/meshing/boundarylayer.hpp"
#include <inctcl.hpp>