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

[occ] ZRefinement

Browse Source
This commit is contained in:
Christopher Lackner 2022-02-23 12:22:19 +01:00
parent 8f8a4a6dc8
commit ea7e980c8d
5 changed files with 201 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
libsrc/general/hashtabl.hpp
View File

@ -150,6 +150,13 @@ public:
return cont.Get (bnr, pos);
}
T & Get (const INDEX_2 & ahash)
{
int bnr = HashValue (ahash);
int pos = Position (bnr, ahash);
return cont.Get (bnr, pos);
}
///
bool Used (const INDEX_2 & ahash) const
{
@ -214,9 +221,14 @@ public:
int BagNr() const { return bagnr; }
int Pos() const { return pos; }
void operator++ (int)
Iterator operator++ (int)
{
Iterator it(ht, bagnr, pos);
++(*this);
return it;
}
Iterator& operator++()
{
// cout << "begin Operator ++: bagnr = " << bagnr << " - pos = " << pos << endl;
pos++;
while (bagnr < ht.GetNBags() &&
pos == ht.GetBagSize(bagnr+1))
@ -224,7 +236,12 @@ public:
pos = 0;
bagnr++;
}
// cout << "end Operator ++: bagnr = " << bagnr << " - pos = " << pos << endl;
return *this;
}
std::pair<INDEX_2, T> operator*()
{
return std::make_pair(ht.hash[bagnr][pos], ht.cont[bagnr][pos]);
}
bool operator != (int i) const
@ -246,6 +263,18 @@ public:
return GetNBags();
}
Iterator begin () const
{
Iterator it(*this, 0, -1);
it++;
return it;
}
int end() const
{
return GetNBags();
}
void GetData (const Iterator & it,
INDEX_2 & ahash, T & acont) const
{

2
libsrc/general/table.hpp
View File

@ -201,6 +201,8 @@ public:
inline const T & Get (int i, int nr) const
{ return ((T*)data.Get(i).col)[nr-1]; }
inline T & Get (int i, int nr)
{ return ((T*)data.Get(i).col)[nr-1]; }
/** Returns pointer to the first element in row i. */
inline const T * GetLine (int i) const

162
libsrc/meshing/meshclass.cpp
View File

@ -6213,7 +6213,169 @@ namespace netgen
*/
}
void Mesh :: ZRefine(const string& name, const Array<double>& slices)
{
auto nr = GetIdentifications().GetNr(name);
auto& identpts = GetIdentifications().GetIdentifiedPoints();
UpdateTopology();
std::map<std::pair<PointIndex, PointIndex>,
Array<PointIndex>> inserted_points;
BitArray mapped_points(GetNV()+1);
mapped_points = false;
// Add new points
for(auto [p1p2, idnr] : identpts)
{
if(idnr != nr)
continue;
auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}];
auto p1 = Point(p1p2.I1());
auto p2 = Point(p1p2.I2());
ipts.Append(p1p2.I1());
mapped_points.SetBit(p1p2.I1());
for(auto slice : slices)
{
auto np = p1 + slice * (p2-p1);
auto npi = AddPoint(np);
ipts.Append(npi);
}
ipts.Append(p1p2.I2());
}
// Split segments
for(auto si : Range(segments))
{
auto& seg = segments[si];
auto p1 = seg[0];
auto p2 = seg[1];
auto c1 = inserted_points.count({p1, p2});
auto c2 = inserted_points.count({p2, p1});
if(c1 == 0 && c2 == 0)
continue;
if(c2)
Swap(p1,p2);
const auto& ipts = inserted_points[{p1,p2}];
if(c2)
seg[1] = ipts[ipts.Size()-2];
else
seg[1] = ipts[1];
for(auto i : Range(size_t(1), ipts.Size()-1))
{
Segment snew = seg;
if(c2)
{
seg[0] = ipts[ipts.Size()-1-i];
seg[1] = ipts[ipts.Size()-2-i];
}
else
{
snew[0] = ipts[i];
snew[1] = ipts[i+1];
}
AddSegment(snew);
}
}
BitArray sel_done(surfelements.Size());
sel_done = false;
// Split surface elements
auto p2sel = CreatePoint2SurfaceElementTable();
for(const auto& [pair, inserted] : inserted_points)
{
for(auto si : p2sel[pair.first])
{
if(sel_done[si])
continue;
sel_done.SetBit(si);
auto sel = surfelements[si];
map<PointIndex, Array<PointIndex>> mapped_points;
int nmapped = 0;
for(auto i : Range(sel.GetNP()))
{
auto p1 = sel[i];
auto p2 = sel[(i+1)%sel.GetNP()];
auto c1 = inserted_points.count({p1, p2});
auto c2 = inserted_points.count({p2, p1});
if(c1 == 0 && c2 == 0)
continue;
if(c2)
Swap(p1, p2);
auto& ipts = inserted_points[{p1, p2}];
auto& a1 = mapped_points[p1];
auto& a2 = mapped_points[p2];
a1 = ipts.Range(0, ipts.Size()-1);
a2 = ipts.Range(1, ipts.Size());
nmapped = ipts.Size()-1;
}
for(auto i : Range(nmapped))
{
Element2d nsel = sel;
for(auto& pi : nsel.PNums())
if(mapped_points.count(pi))
pi = mapped_points[pi][i];
AddSurfaceElement(nsel);
}
if(nmapped)
surfelements[si].Delete();
}
}
// Split volume elements
BitArray vol_done(volelements.Size());
vol_done = false;
auto p2el = CreatePoint2ElementTable(); // mapped_points);
for(const auto& [pair, inserted] : inserted_points)
{
for(auto ei : p2el[pair.first])
{
if(vol_done[ei])
continue;
vol_done.SetBit(ei);
auto el = volelements[ei];
map<PointIndex, Array<PointIndex>> mapped_points;
int nmapped = 0;
NgArray<int> eledges;
topology.GetElementEdges(ei+1, eledges);
for(auto edgei : eledges)
{
int p1, p2;
topology.GetEdgeVertices(edgei, p1, p2);
auto c1 = inserted_points.count({p1, p2});
auto c2 = inserted_points.count({p2, p1});
if(c1 == 0 && c2 == 0)
continue;
if(c2)
Swap(p1, p2);
auto& ipts = inserted_points[{p1, p2}];
auto& a1 = mapped_points[p1];
auto& a2 = mapped_points[p2];
a1 = ipts.Range(0, ipts.Size()-1);
a2 = ipts.Range(1, ipts.Size());
nmapped = ipts.Size()-1;
}
for(auto i : Range(nmapped))
{
Element nel = el;
for(auto& pi : nel.PNums())
if(mapped_points.count(pi))
pi = mapped_points[pi][i];
AddVolumeElement(nel);
}
if(nmapped)
volelements[ei].Delete();
}
}
Compress();
}
void Mesh :: RebuildSurfaceElementLists ()
{

1
libsrc/meshing/meshclass.hpp
View File

@ -499,6 +499,7 @@ namespace netgen
/// Refines mesh and projects points to true surface
// void Refine (int levels, const CSGeometry * geom);
void ZRefine(const string& name, const Array<double>& slices);
bool BoundaryEdge (PointIndex pi1, PointIndex pi2) const
{

2
libsrc/meshing/python_mesh.cpp
View File

@ -1127,6 +1127,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
self.UpdateTopology();
}),py::call_guard<py::gil_scoped_release>())
.def("ZRefine", &Mesh::ZRefine)
.def ("SecondOrder", FunctionPointer
([](Mesh & self)
{