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

PointIndex in Identifications

Browse Source
This commit is contained in:
Joachim Schoeberl 2024-12-31 21:26:04 +01:00
parent fe21b0bb8b
commit bcbd390f7d
13 changed files with 188 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
libsrc/csg/zrefine.cpp
View File

@ -265,13 +265,20 @@ namespace netgen
auto & identpts = auto & identpts =
mesh.GetIdentifications().GetIdentifiedPoints (); mesh.GetIdentifications().GetIdentifiedPoints ();
/*
for (int i = 1; i <= identpts.GetNBags(); i++) for (int i = 1; i <= identpts.GetNBags(); i++)
for (int j = 1; j <= identpts.GetBagSize(i); j++) for (int j = 1; j <= identpts.GetBagSize(i); j++)
{ {
INDEX_3 pair; INDEX_3 pair;
int dummy; int dummy;
identpts.GetData(i, j, pair, dummy); identpts.GetData(i, j, pair, dummy);
auto idnr = pair[2]; */
for (auto [hash, val] : identpts)\
{
auto [hash_pts, idnr] = hash;
auto [pi1, pi2] = hash_pts;
// auto idnr = pair[2];
const CloseSurfaceIdentification * csid = const CloseSurfaceIdentification * csid =
dynamic_cast<const CloseSurfaceIdentification*> dynamic_cast<const CloseSurfaceIdentification*>
(geom->identifications.Get(idnr)); (geom->identifications.Get(idnr));
@ -282,17 +289,25 @@ namespace netgen
if (first_id.Test (idnr)) if (first_id.Test (idnr))
{ {
first_id.Clear(idnr); first_id.Clear(idnr);
/*
ref_uniform.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels())); ref_uniform.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels()));
ref_singular.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels1())); ref_singular.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels1()));
ref_singular.Append (INDEX_3 (pair.I2(), pair.I1(), csid->RefLevels2())); ref_singular.Append (INDEX_3 (pair.I2(), pair.I1(), csid->RefLevels2()));
*/
ref_uniform.Append (INDEX_3 (pi1, pi2, csid->RefLevels()));
ref_singular.Append (INDEX_3 (pi1, pi2, csid->RefLevels1()));
ref_singular.Append (INDEX_3 (pi2, pi1, csid->RefLevels2()));
} }
} }
else else
{ {
//const NgArray<double> & slices = csid->GetSlices(); //const NgArray<double> & slices = csid->GetSlices();
INDEX_4 i4; INDEX_4 i4;
i4[0] = pair.I1(); // i4[0] = pair.I1();
i4[1] = pair.I2(); // i4[1] = pair.I2();
i4[0] = pi1;
i4[1] = pi2;
i4[2] = idnr; i4[2] = idnr;
i4[3] = csid->GetSlices().Size(); i4[3] = csid->GetSlices().Size();
ref_slices.Append (i4); ref_slices.Append (i4);

4
libsrc/include/nginterface_v2_impl.hpp
View File

@ -339,8 +339,8 @@ NGX_INLINE DLL_HEADER Ng_Buffer<int[2]> Ngx_Mesh :: GetPeriodicVertices(int idnr
mesh->GetIdentifications().GetPairs (idnr+1, apairs); mesh->GetIdentifications().GetPairs (idnr+1, apairs);
for(auto& ind : apairs) for(auto& ind : apairs)
{ {
ind.I1()--; ind.I1() -= IndexBASE<PointIndex>();
ind.I2()--; ind.I2() -= IndexBASE<PointIndex>();
} }
typedef int ti2[2]; typedef int ti2[2];
return { apairs.Size(), (ti2*)(void*)apairs.Release() }; return { apairs.Size(), (ti2*)(void*)apairs.Release() };

10
libsrc/interface/nginterface_v2.cpp
View File

@ -1249,9 +1249,15 @@ int Ngx_Mesh::GetElementOrder (int enr) const
void Ngx_Mesh::GetElementOrders (int enr, int * ox, int * oy, int * oz) const void Ngx_Mesh::GetElementOrders (int enr, int * ox, int * oy, int * oz) const
{ {
if (mesh->GetDimension() == 3) if (mesh->GetDimension() == 3)
mesh->VolumeElement(enr).GetOrder(*ox, *oy, *oz); {
ElementIndex ei = IndexBASE<ElementIndex>() + enr-1;
mesh->VolumeElement(ei).GetOrder(*ox, *oy, *oz);
}
else else
mesh->SurfaceElement(enr).GetOrder(*ox, *oy, *oz); {
SurfaceElementIndex sei = IndexBASE<SurfaceElementIndex>() + enr-1;
mesh->SurfaceElement(sei).GetOrder(*ox, *oy, *oz);
}
} }
void Ngx_Mesh::SetElementOrder (int enr, int order) void Ngx_Mesh::SetElementOrder (int enr, int order)

39
libsrc/interface/writeuser.cpp
View File

@ -94,9 +94,13 @@ void WriteNeutralFormat (const Mesh & mesh,
if (mesh.GetDimension() == 3) if (mesh.GetDimension() == 3)
{ {
outfile << ne << "\n"; outfile << ne << "\n";
/*
for (int i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
{ {
Element el = mesh.VolumeElement(i); Element el = mesh.VolumeElement(i);
*/
for (Element el : mesh.VolumeElements())
{
if (inverttets) if (inverttets)
el.Invert(); el.Invert();
outfile.width(4); outfile.width(4);
@ -112,9 +116,13 @@ void WriteNeutralFormat (const Mesh & mesh,
} }
outfile << nse << "\n"; outfile << nse << "\n";
/*
for (int i = 1; i <= nse; i++) for (int i = 1; i <= nse; i++)
{ {
Element2d el = mesh.SurfaceElement(i); Element2d el = mesh.SurfaceElement(i);
*/
for (Element2d el : mesh.SurfaceElements())
{
if (invertsurf) if (invertsurf)
el.Invert(); el.Invert();
outfile.width(4); outfile.width(4);
@ -583,9 +591,13 @@ void WriteFEPPFormat (const Mesh & mesh,
outfile << ne << "\n"; outfile << ne << "\n";
/*
for (i = 1; i <= ne; i++) for (i = 1; i <= ne; i++)
{ {
const Element & el = mesh.VolumeElement(i); const Element & el = mesh.VolumeElement(i);
*/
for (const Element & el : mesh.VolumeElements())
{
outfile.width(4); outfile.width(4);
outfile << el.GetIndex() << " "; outfile << el.GetIndex() << " ";
outfile.width(4); outfile.width(4);
@ -677,7 +689,6 @@ void WriteEdgeElementFormat (const Mesh & mesh,
int nelements = mesh.GetNE(); int nelements = mesh.GetNE();
int nsurfelem = mesh.GetNSE(); int nsurfelem = mesh.GetNSE();
int nedges = top->GetNEdges(); int nedges = top->GetNEdges();
int i, j;
int inverttets = mparam.inverttets; int inverttets = mparam.inverttets;
int invertsurf = mparam.inverttrigs; int invertsurf = mparam.inverttrigs;
@ -692,7 +703,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
// vertices with coordinates // vertices with coordinates
outfile << npoints << "\n"; outfile << npoints << "\n";
for (i = 1; i <= npoints; i++) for (int i = 1; i <= npoints; i++)
{ {
const Point3d & p = mesh.Point(i); const Point3d & p = mesh.Point(i);
@ -706,16 +717,24 @@ void WriteEdgeElementFormat (const Mesh & mesh,
// element - edge - list // element - edge - list
outfile << nelements << " " << nedges << "\n"; outfile << nelements << " " << nedges << "\n";
/*
for (i = 1; i <= nelements; i++) for (i = 1; i <= nelements; i++)
{ {
Element el = mesh.VolumeElement(i); Element el = mesh.VolumeElement(i);
*/
for (ElementIndex ei : Range(mesh.VolumeElements()))
{
int i = ei-IndexBASE(ei)+1;
Element el = mesh.VolumeElement(ei);
if (inverttets) if (inverttets)
el.Invert(); el.Invert();
outfile.width(4); outfile.width(4);
outfile << el.GetIndex() << " "; outfile << el.GetIndex() << " ";
outfile.width(8); outfile.width(8);
outfile << el.GetNP(); outfile << el.GetNP();
for (j = 1; j <= el.GetNP(); j++) for (int j = 1; j <= el.GetNP(); j++)
{ {
outfile << " "; outfile << " ";
outfile.width(8); outfile.width(8);
@ -723,11 +742,11 @@ void WriteEdgeElementFormat (const Mesh & mesh,
} }
// top->GetElementEdges(i,edges); // top->GetElementEdges(i,edges);
auto eledges = top->GetEdges(ElementIndex(i-1)); auto eledges = top->GetEdges(ei);
outfile << endl << " "; outfile << endl << " ";
outfile.width(8); outfile.width(8);
outfile << eledges.Size(); outfile << eledges.Size();
for (j=1; j <= eledges.Size(); j++) for (int j=1; j <= eledges.Size(); j++)
{ {
outfile << " "; outfile << " ";
outfile.width(8); outfile.width(8);
@ -738,7 +757,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
// orientation: // orientation:
top->GetElementEdgeOrientations(i,edges); top->GetElementEdgeOrientations(i,edges);
outfile << " "; outfile << " ";
for (j=1; j <= edges.Size(); j++) for (int j=1; j <= edges.Size(); j++)
{ {
outfile << " "; outfile << " ";
outfile.width(8); outfile.width(8);
@ -749,7 +768,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
// surface element - edge - list (with boundary conditions) // surface element - edge - list (with boundary conditions)
outfile << nsurfelem << "\n"; outfile << nsurfelem << "\n";
for (i = 1; i <= nsurfelem; i++) for (int i = 1; i <= nsurfelem; i++)
{ {
Element2d el = mesh.SurfaceElement(i); Element2d el = mesh.SurfaceElement(i);
if (invertsurf) if (invertsurf)
@ -758,7 +777,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
outfile.width(8); outfile.width(8);
outfile << el.GetNP(); outfile << el.GetNP();
for (j = 1; j <= el.GetNP(); j++) for (int j = 1; j <= el.GetNP(); j++)
{ {
outfile << " "; outfile << " ";
outfile.width(8); outfile.width(8);
@ -769,7 +788,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
outfile << endl << " "; outfile << endl << " ";
outfile.width(8); outfile.width(8);
outfile << edges.Size(); outfile << edges.Size();
for (j=1; j <= edges.Size(); j++) for (int j=1; j <= edges.Size(); j++)
{ {
outfile << " "; outfile << " ";
outfile.width(8); outfile.width(8);
@ -782,7 +801,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
// int v1, v2; // int v1, v2;
// edge - vertex - list // edge - vertex - list
outfile << nedges << "\n"; outfile << nedges << "\n";
for (i=1; i <= nedges; i++) for (int i=1; i <= nedges; i++)
{ {
// top->GetEdgeVertices(i,v1,v2); // top->GetEdgeVertices(i,v1,v2);
auto [v1,v2] = top->GetEdgeVertices(i-1); auto [v1,v2] = top->GetEdgeVertices(i-1);

34
libsrc/meshing/meshclass.cpp
View File

@ -6438,22 +6438,24 @@ namespace netgen
mapped_points = false; mapped_points = false;
// Add new points // Add new points
for(auto [p1p2, dummy] : identpts) for(auto [hash, dummy] : identpts)
{ {
if(p1p2[2] != nr) auto [hash_pts, hash_nr] = hash;
if(hash_nr != nr)
continue; continue;
auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}]; // auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}];
auto p1 = Point(p1p2.I1()); auto& ipts = inserted_points[ { hash_pts[0], hash_pts[1] }];
auto p2 = Point(p1p2.I2()); auto p1 = Point(hash_pts.I1());
ipts.Append(p1p2.I1()); auto p2 = Point(hash_pts.I2());
mapped_points.SetBit(p1p2.I1()); ipts.Append(hash_pts.I1());
mapped_points.SetBit(hash_pts.I1());
for(auto slice : slices) for(auto slice : slices)
{ {
auto np = p1 + slice * (p2-p1); auto np = p1 + slice * (p2-p1);
auto npi = AddPoint(np); auto npi = AddPoint(np);
ipts.Append(npi); ipts.Append(npi);
} }
ipts.Append(p1p2.I2()); ipts.Append(hash_pts.I2());
} }
// Split segments // Split segments
@ -6644,21 +6646,21 @@ namespace netgen
void Mesh :: CalcMinMaxAngle (double badellimit, double * retvalues) void Mesh :: CalcMinMaxAngle (double badellimit, double * retvalues)
{ {
int i, j;
int lpi1, lpi2, lpi3, lpi4; int lpi1, lpi2, lpi3, lpi4;
double phimax = 0, phimin = 10; double phimax = 0, phimin = 10;
double facephimax = 0, facephimin = 10; double facephimax = 0, facephimin = 10;
int illegaltets = 0, negativetets = 0, badtets = 0; int illegaltets = 0, negativetets = 0, badtets = 0;
for (i = 1; i <= GetNE(); i++) // for (int i = 1; i <= GetNE(); i++)
for (ElementIndex ei : Range(VolumeElements()))
{ {
int badel = 0; int badel = 0;
Element & el = VolumeElement(i); Element & el = VolumeElement(ei);
if (el.GetType() != TET) if (el.GetType() != TET)
{ {
VolumeElement(i).Flags().badel = 0; VolumeElement(ei).Flags().badel = 0;
continue; continue;
} }
@ -6673,8 +6675,8 @@ namespace netgen
{ {
badel = 1; badel = 1;
illegaltets++; illegaltets++;
(*testout) << "illegal tet: " << i << " "; (*testout) << "illegal tet: " << ei << " ";
for (j = 1; j <= el.GetNP(); j++) for (int j = 1; j <= el.GetNP(); j++)
(*testout) << el.PNum(j) << " "; (*testout) << el.PNum(j) << " ";
(*testout) << endl; (*testout) << endl;
} }
@ -6713,7 +6715,7 @@ namespace netgen
// angles in faces // angles in faces
for (j = 1; j <= 4; j++) for (int j = 1; j <= 4; j++)
{ {
Element2d face(TRIG); Element2d face(TRIG);
el.GetFace (j, face); el.GetFace (j, face);
@ -6740,7 +6742,7 @@ namespace netgen
} }
VolumeElement(i).Flags().badel = badel; VolumeElement(ei).Flags().badel = badel;
if (badel) badtets++; if (badel) badtets++;
} }

4
libsrc/meshing/meshfunc.cpp
View File

@ -578,8 +578,10 @@ namespace netgen
el.SetIndex(m_.domain); el.SetIndex(m_.domain);
mesh.AddVolumeElement(el); mesh.AddVolumeElement(el);
} }
// for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints())
// mesh.GetIdentifications().Add(pmap[p1p2[0]], pmap[p1p2[1]], p1p2[2]);
for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints()) for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints())
mesh.GetIdentifications().Add(pmap[p1p2[0]], pmap[p1p2[1]], p1p2[2]); mesh.GetIdentifications().Add( pmap[ get<0>(p1p2)[0] ], pmap[ get<0>(p1p2)[1]] , get<1>(p1p2) );
for(auto i : Range(m.GetIdentifications().GetMaxNr())) for(auto i : Range(m.GetIdentifications().GetMaxNr()))
{ {
mesh.GetIdentifications().SetType(i+1, m.GetIdentifications().GetType(i+1)); mesh.GetIdentifications().SetType(i+1, m.GetIdentifications().GetType(i+1));

67
libsrc/meshing/meshtype.cpp
View File

@ -2708,7 +2708,8 @@ namespace netgen
void Identifications :: DoArchive (Archive & ar) void Identifications :: DoArchive (Archive & ar)
{ {
ar & maxidentnr; ar & maxidentnr;
ar & identifiedpoints & identifiedpoints_nr; // ar & identifiedpoints & identifiedpoints_nr;
#pragma message( "Archive CloseHadhTable missing " __FILE__ )
ar & idpoints_table; ar & idpoints_table;
if (ar.Output()) if (ar.Output())
@ -2735,11 +2736,11 @@ namespace netgen
void Identifications :: Add (PointIndex pi1, PointIndex pi2, int identnr) void Identifications :: Add (PointIndex pi1, PointIndex pi2, int identnr)
{ {
// (*testout) << "Identification::Add, pi1 = " << pi1 << ", pi2 = " << pi2 << ", identnr = " << identnr << endl; // (*testout) << "Identification::Add, pi1 = " << pi1 << ", pi2 = " << pi2 << ", identnr = " << identnr << endl;
INDEX_2 pair (pi1, pi2); PointIndices<2> pair (pi1, pi2);
identifiedpoints.Set (pair, identnr); identifiedpoints.Set (pair, identnr);
INDEX_3 tripl (pi1, pi2, identnr); // INDEX_3 tripl (pi1, pi2, identnr);
identifiedpoints_nr.Set (tripl, 1); identifiedpoints_nr.Set ( { { pi1, pi2 }, identnr }, 1);
if (identnr > maxidentnr) maxidentnr = identnr; if (identnr > maxidentnr) maxidentnr = identnr;
names.SetSize(maxidentnr); names.SetSize(maxidentnr);
@ -2762,8 +2763,9 @@ namespace netgen
bool Identifications :: Get (PointIndex pi1, PointIndex pi2, int nr) const bool Identifications :: Get (PointIndex pi1, PointIndex pi2, int nr) const
{ {
INDEX_3 tripl(pi1, pi2, nr); // INDEX_3 tripl(pi1, pi2, nr);
if (identifiedpoints_nr.Used (tripl)) // if (identifiedpoints_nr.Used (tripl))
if (identifiedpoints_nr.Used ( { { pi1, pi1 }, nr } ) )
return 1; return 1;
else else
return 0; return 0;
@ -2803,23 +2805,30 @@ namespace netgen
{ {
cout << "getmap, identnr = " << identnr << endl; cout << "getmap, identnr = " << identnr << endl;
/*
for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++) for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++)
for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++) for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++)
*/
for (auto [hash, val] : identifiedpoints_nr)
{ {
/*
INDEX_3 i3; INDEX_3 i3;
int dummy; int dummy;
identifiedpoints_nr.GetData (i, j, i3, dummy); identifiedpoints_nr.GetData (i, j, i3, dummy);
*/
if (i3.I3() == identnr || !identnr) auto [hash_pts, hash_nr] = hash;
if (hash_nr == identnr || !identnr)
{ {
/* /*
identmap.Elem(i3.I1()) = i3.I2(); identmap.Elem(i3.I1()) = i3.I2();
if(symmetric) if(symmetric)
identmap.Elem(i3.I2()) = i3.I1(); identmap.Elem(i3.I2()) = i3.I1();
*/ */
identmap[i3.I1()] = i3.I2(); identmap[hash_pts.I1()] = hash_pts.I2();
if(symmetric) if(symmetric)
identmap[i3.I2()] = i3.I1(); identmap[hash_pts.I2()] = hash_pts.I1();
} }
} }
} }
@ -2830,8 +2839,12 @@ namespace netgen
Array<INDEX_3> Identifications :: GetPairs () const Array<INDEX_3> Identifications :: GetPairs () const
{ {
Array<INDEX_3> pairs; Array<INDEX_3> pairs;
for(auto [i3, dummy] : identifiedpoints_nr) for(auto [hash, dummy] : identifiedpoints_nr)
pairs.Append(i3); // pairs.Append(i3);
{
auto [pts,nr] = hash;
pairs.Append ( { pts[0], pts[1], nr } );
}
return pairs; return pairs;
} }
@ -2841,6 +2854,8 @@ namespace netgen
identpairs.SetSize(0); identpairs.SetSize(0);
if (identnr == 0) if (identnr == 0)
{
/*
for (int i = 1; i <= identifiedpoints.GetNBags(); i++) for (int i = 1; i <= identifiedpoints.GetNBags(); i++)
for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++) for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++)
{ {
@ -2849,7 +2864,13 @@ namespace netgen
identifiedpoints.GetData (i, j, i2, nr); identifiedpoints.GetData (i, j, i2, nr);
identpairs.Append (i2); identpairs.Append (i2);
} }
*/
for (auto [hash,val] : identifiedpoints)
identpairs.Append (hash);
}
else else
{
/*
for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++) for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++)
for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++) for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++)
{ {
@ -2860,11 +2881,20 @@ namespace netgen
if (i3.I3() == identnr) if (i3.I3() == identnr)
identpairs.Append (INDEX_2(i3.I1(), i3.I2())); identpairs.Append (INDEX_2(i3.I1(), i3.I2()));
} }
*/
for (auto [hash,val] : identifiedpoints_nr)
{
auto [hash_pts, hash_nr] = hash;
if (hash_nr == identnr)
identpairs.Append (hash_pts);
}
}
} }
void Identifications :: SetMaxPointNr (int maxpnum) void Identifications :: SetMaxPointNr (int maxpnum)
{ {
/*
for (int i = 1; i <= identifiedpoints.GetNBags(); i++) for (int i = 1; i <= identifiedpoints.GetNBags(); i++)
for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++) for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++)
{ {
@ -2878,6 +2908,18 @@ namespace netgen
identifiedpoints.SetData (i, j, i2, -1); identifiedpoints.SetData (i, j, i2, -1);
} }
} }
*/
// can we get data by reference ?
for (auto [hash,data] : identifiedpoints)
{
if (hash.I1() > IndexBASE<PointIndex>()+maxpnum-1 ||
hash.I2() > IndexBASE<PointIndex>()+maxpnum-1)
{
identifiedpoints[hash] = -1;
}
}
} }
// Map points in the identifications to new point numbers // Map points in the identifications to new point numbers
@ -2900,7 +2942,8 @@ namespace netgen
{ {
ost << "Identifications:" << endl; ost << "Identifications:" << endl;
ost << "pairs: " << endl << identifiedpoints << endl; ost << "pairs: " << endl << identifiedpoints << endl;
ost << "pairs and nr: " << endl << identifiedpoints_nr << endl; // ost << "pairs and nr: " << endl << identifiedpoints_nr << endl;
#pragma message( "Can't ostream a tuple " __FILE__ )
ost << "table: " << endl << idpoints_table << endl; ost << "table: " << endl << idpoints_table << endl;
} }

12
libsrc/meshing/meshtype.hpp
View File

@ -244,7 +244,7 @@ namespace netgen
class PointIndex : public Index<int,PointIndex,1> class PointIndex : public Index<int,PointIndex,0>
{ {
public: public:
using Index::Index; using Index::Index;
@ -1763,13 +1763,15 @@ namespace netgen
class Mesh & mesh; class Mesh & mesh;
/// identify points (thin layers, periodic b.c.) /// identify points (thin layers, periodic b.c.)
INDEX_2_HASHTABLE<int> identifiedpoints; // INDEX_2_HASHTABLE<int> identifiedpoints;
ClosedHashTable<PointIndices<2>, int> identifiedpoints;
/// the same, with info about the id-nr /// the same, with info about the id-nr
INDEX_3_HASHTABLE<int> identifiedpoints_nr; // INDEX_3_HASHTABLE<int> identifiedpoints_nr;
ClosedHashTable<std::tuple<PointIndices<2>, int>, int> identifiedpoints_nr;
/// sorted by identification nr /// sorted by identification nr
TABLE<INDEX_2> idpoints_table; TABLE<PointIndices<2>> idpoints_table;
NgArray<ID_TYPE> type; NgArray<ID_TYPE> type;
@ -1805,7 +1807,7 @@ namespace netgen
// bool HasIdentifiedPoints() const { return identifiedpoints != nullptr; } // bool HasIdentifiedPoints() const { return identifiedpoints != nullptr; }
/// ///
INDEX_3_HASHTABLE<int> & GetIdentifiedPoints () auto & GetIdentifiedPoints ()
{ {
return identifiedpoints_nr; return identifiedpoints_nr;
} }

5
libsrc/meshing/python_mesh.cpp
View File

@ -1266,7 +1266,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
py::list points; py::list points;
for(const auto& pair : self.GetIdentifications().GetIdentifiedPoints()) for(const auto& pair : self.GetIdentifications().GetIdentifiedPoints())
{ {
py::tuple pnts = py::make_tuple(pair.first.I1(), pair.first.I2()); // py::tuple pnts = py::make_tuple(pair.first.I1(), pair.first.I2());
auto [pi1, pi2] = get<0> (pair.first);
py::tuple pnts = py::make_tuple(pi1, pi2);
points.append(pnts); points.append(pnts);
} }
return points; return points;

1
libsrc/meshing/secondorder.cpp
View File

@ -558,6 +558,7 @@ namespace netgen
// (*testout) << "bad els: " << endl; // (*testout) << "bad els: " << endl;
wrongels = 0; wrongels = 0;
for (int i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
{ {
if (!illegalels.Test(i) && if (!illegalels.Test(i) &&
mesh.VolumeElement(i). mesh.VolumeElement(i).

15
libsrc/meshing/topology.cpp
View File

@ -1829,8 +1829,8 @@ namespace netgen
void MeshTopology :: GetElementFaces (int elnr, NgArray<int> & elfaces, bool withorientation) const void MeshTopology :: GetElementFaces (int elnr, NgArray<int> & elfaces, bool withorientation) const
{ {
int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType());
ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1); ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
int nfa = GetNFaces (mesh->VolumeElement(ei).GetType());
elfaces.SetSize (nfa); elfaces.SetSize (nfa);
@ -1859,7 +1859,8 @@ namespace netgen
void MeshTopology :: GetElementEdgeOrientations (int elnr, NgArray<int> & eorient) const void MeshTopology :: GetElementEdgeOrientations (int elnr, NgArray<int> & eorient) const
{ {
int ned = GetNEdges (mesh->VolumeElement(elnr).GetType()); ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
int ned = GetNEdges (mesh->VolumeElement(ei).GetType());
eorient.SetSize (ned); eorient.SetSize (ned);
for (int i = 1; i <= ned; i++) for (int i = 1; i <= ned; i++)
// eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; // eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1;
@ -1869,7 +1870,8 @@ namespace netgen
void MeshTopology :: GetElementFaceOrientations (int elnr, NgArray<int> & forient) const void MeshTopology :: GetElementFaceOrientations (int elnr, NgArray<int> & forient) const
{ {
int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType()); ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
int nfa = GetNFaces (mesh->VolumeElement(ei).GetType());
forient.SetSize (nfa); forient.SetSize (nfa);
for (int i = 1; i <= nfa; i++) for (int i = 1; i <= nfa; i++)
// forient.Elem(i) = faces.Get(elnr)[i-1].forient; // forient.Elem(i) = faces.Get(elnr)[i-1].forient;
@ -2111,7 +2113,9 @@ namespace netgen
int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const
{ {
const Element & el = mesh->VolumeElement (elnr); ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
const Element & el = mesh->VolumeElement (ei);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());
int k = locedgenr; int k = locedgenr;
@ -2123,7 +2127,8 @@ namespace netgen
int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const
{ {
const Element & el = mesh->VolumeElement (elnr); ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
const Element & el = mesh->VolumeElement (ei);
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType()); const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType());

43
libsrc/visualization/vsmesh.cpp
View File

@ -515,20 +515,21 @@ namespace netgen
if (vispar.drawelementnumbers) if (vispar.drawelementnumbers)
{ {
NgArray<int> v; NgArray<int> v;
for (int i = 1; i <= mesh->GetNE(); i++) // for (int i = 1; i <= mesh->GetNE(); i++)
for (ElementIndex ei : Range(mesh->VolumeElements()))
{ {
// const ELEMENTTYPE & eltype = mesh->ElementType(i); // const ELEMENTTYPE & eltype = mesh->ElementType(i);
NgArray<int> pnums; NgArray<int> pnums;
Point3d p; Point3d p;
const Element & el = mesh->VolumeElement (i); const Element & el = mesh->VolumeElement (ei);
if ( ! el.PNum(5)) // eltype == TET ) if ( ! el.PNum(5)) // eltype == TET )
{ {
pnums.SetSize(4); pnums.SetSize(4);
for( int j = 0; j < pnums.Size(); j++) for( int j = 0; j < pnums.Size(); j++)
pnums[j] = mesh->VolumeElement(i).PNum(j+1); pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p1 = mesh->Point(pnums[0]);
@ -541,7 +542,7 @@ namespace netgen
{ {
pnums.SetSize(5); pnums.SetSize(5);
for( int j = 0; j < pnums.Size(); j++) for( int j = 0; j < pnums.Size(); j++)
pnums[j] = mesh->VolumeElement(i).PNum(j+1); pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p1 = mesh->Point(pnums[0]);
@ -559,7 +560,7 @@ namespace netgen
{ {
pnums.SetSize(6); pnums.SetSize(6);
for( int j = 0; j < pnums.Size(); j++) for( int j = 0; j < pnums.Size(); j++)
pnums[j] = mesh->VolumeElement(i).PNum(j+1); pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p1 = mesh->Point(pnums[0]);
const Point3d & p2 = mesh->Point(pnums[1]); const Point3d & p2 = mesh->Point(pnums[1]);
@ -574,7 +575,7 @@ namespace netgen
{ {
pnums.SetSize(8); pnums.SetSize(8);
for( int j = 0; j < pnums.Size(); j++) for( int j = 0; j < pnums.Size(); j++)
pnums[j] = mesh->VolumeElement(i).PNum(j+1); pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p1 = mesh->Point(pnums[0]);
const Point3d & p2 = mesh->Point(pnums[1]); const Point3d & p2 = mesh->Point(pnums[1]);
@ -589,7 +590,7 @@ namespace netgen
} }
glRasterPos3d (p.X(), p.Y(), p.Z()); glRasterPos3d (p.X(), p.Y(), p.Z());
snprintf (buf, size(buf), "%d", i); snprintf (buf, size(buf), "%d", ei-IndexBASE(ei));
// glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); // glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
MyOpenGLText (buf); MyOpenGLText (buf);
@ -617,14 +618,15 @@ namespace netgen
static float badelcol[] = { 1.0f, 0.0f, 1.0f, 1.0f }; static float badelcol[] = { 1.0f, 0.0f, 1.0f, 1.0f };
glLineWidth (1.0f); glLineWidth (1.0f);
for (int i = 1; i <= mesh->GetNE(); i++) //for (int i = 1; i <= mesh->GetNE(); i++)
for (ElementIndex ei : Range(mesh->VolumeElements()))
{ {
if (mesh->VolumeElement(i).Flags().badel || if (mesh->VolumeElement(ei).Flags().badel ||
mesh->VolumeElement(i).Flags().illegal || mesh->VolumeElement(ei).Flags().illegal ||
(i == vispar.drawelement)) (ei-IndexBASE(ei) == vispar.drawelement))
{ {
// copy to be thread-safe // copy to be thread-safe
Element el = mesh->VolumeElement (i); Element el = mesh->VolumeElement (ei);
el.GetSurfaceTriangles (faces); el.GetSurfaceTriangles (faces);
glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, badelcol); glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, badelcol);
@ -674,9 +676,9 @@ namespace netgen
} }
for (int i = 1; i <= mesh->GetNE(); i++) for (ElementIndex ei : Range(mesh->VolumeElements()))
{ {
Element el = mesh->VolumeElement (i); Element el = mesh->VolumeElement (ei);
int hascp = 0; int hascp = 0;
for (int j = 1; j <= el.GetNP(); j++) for (int j = 1; j <= el.GetNP(); j++)
if (el.PNum(j) == vispar.centerpoint) if (el.PNum(j) == vispar.centerpoint)
@ -684,7 +686,7 @@ namespace netgen
if (hascp) if (hascp)
{ {
(*testout) << "draw el " << i << " : "; (*testout) << "draw el " << ei << " : ";
for (int j = 1; j <= el.GetNP(); j++) for (int j = 1; j <= el.GetNP(); j++)
(*testout) << el.PNum(j) << " "; (*testout) << el.PNum(j) << " ";
(*testout) << endl; (*testout) << endl;
@ -861,16 +863,21 @@ namespace netgen
auto & idpts = auto & idpts =
mesh->GetIdentifications().GetIdentifiedPoints(); mesh->GetIdentifications().GetIdentifiedPoints();
/*
for (int i = 1; i <= idpts.GetNBags(); i++) for (int i = 1; i <= idpts.GetNBags(); i++)
for (int j = 1; j <= idpts.GetBagSize(i); j++) for (int j = 1; j <= idpts.GetBagSize(i); j++)
{ {
INDEX_3 pts; INDEX_3 pts;
int dummy; // , val; int dummy; // , val;
idpts.GetData (i, j, pts, dummy); idpts.GetData (i, j, pts, dummy);
*/
for (auto [hash, val] : idpts)
{
auto [hash_pts, hash_nr] = hash;
auto [pi1, pi2] = hash_pts;
// val = pts[2]; // val = pts[2];
const Point3d & p1 = mesh->Point(pts.I1()); const Point3d & p1 = mesh->Point(pi1);
const Point3d & p2 = mesh->Point(pts.I2()); const Point3d & p2 = mesh->Point(pi2);
glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
identifiedcol); identifiedcol);