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

register 1,2,3d elements to numpy to be used in arrays

Browse Source
This commit is contained in:
Christopher Lackner 2022-09-13 15:12:42 +02:00
parent 79b266fac3
commit c18a317702
15 changed files with 210 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
libsrc/interface/nginterface_v2.cpp
View File

@ -701,9 +701,9 @@ namespace netgen
{ {
int hpelnr = -1; int hpelnr = -1;
if (mesh->GetDimension() == 2) if (mesh->GetDimension() == 2)
hpelnr = mesh->SurfaceElement(ei).hp_elnr; hpelnr = mesh->SurfaceElement(ei).GetHpElnr();
else else
hpelnr = mesh->VolumeElement(ei).hp_elnr; hpelnr = mesh->VolumeElement(ei).GetHpElnr();
if (hpelnr < 0) if (hpelnr < 0)
throw NgException("Ngx_Mesh::GetHPElementLevel: Wrong hp-element number!"); throw NgException("Ngx_Mesh::GetHPElementLevel: Wrong hp-element number!");

14
libsrc/meshing/curvedelems.cpp
View File

@ -1632,7 +1632,7 @@ namespace netgen
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
return mesh.coarsemesh->GetCurvedElements().IsSurfaceElementCurved (hpref_el.coarse_elnr); return mesh.coarsemesh->GetCurvedElements().IsSurfaceElementCurved (hpref_el.coarse_elnr);
} }
@ -1679,7 +1679,7 @@ namespace netgen
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
// xi umrechnen // xi umrechnen
double lami[4]; double lami[4];
@ -2428,7 +2428,7 @@ namespace netgen
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
return mesh.coarsemesh->GetCurvedElements().IsElementCurved (hpref_el.coarse_elnr); return mesh.coarsemesh->GetCurvedElements().IsElementCurved (hpref_el.coarse_elnr);
} }
@ -2480,7 +2480,7 @@ namespace netgen
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
return mesh.coarsemesh->GetCurvedElements().IsElementHighOrder (hpref_el.coarse_elnr); return mesh.coarsemesh->GetCurvedElements().IsElementHighOrder (hpref_el.coarse_elnr);
} }
@ -2519,7 +2519,7 @@ namespace netgen
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
// xi umrechnen // xi umrechnen
double lami[8]; double lami[8];
@ -4065,7 +4065,7 @@ namespace netgen
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
// xi umrechnen // xi umrechnen
T lami[4]; T lami[4];
@ -4529,7 +4529,7 @@ namespace netgen
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
// xi umrechnen // xi umrechnen
T lami[8]; T lami[8];

8
libsrc/meshing/hprefinement.cpp
View File

@ -1403,7 +1403,7 @@ namespace netgen
Element2d el(hpel.np); Element2d el(hpel.np);
for(int j=0;j<hpel.np;j++) for(int j=0;j<hpel.np;j++)
el.PNum(j+1) = hpel.pnums[j]; el.PNum(j+1) = hpel.pnums[j];
el.hp_elnr = i; el.SetHpElnr(i);
el.SetIndex(hpel.index); el.SetIndex(hpel.index);
if(setorders) if(setorders)
el.SetOrder(act_ref+1,act_ref+1,0); el.SetOrder(act_ref+1,act_ref+1,0);
@ -1421,7 +1421,7 @@ namespace netgen
for(int j=0;j<hpel.np;j++) for(int j=0;j<hpel.np;j++)
el.PNum(j+1) = hpel.pnums[j]; el.PNum(j+1) = hpel.pnums[j];
el.SetIndex(hpel.index); el.SetIndex(hpel.index);
el.hp_elnr = i; el.SetHpElnr(i);
if(setorders) if(setorders)
el.SetOrder(act_ref+1,act_ref+1,act_ref+1); el.SetOrder(act_ref+1,act_ref+1,act_ref+1);
if((*mesh.coarsemesh)[ElementIndex{hpel.coarse_elnr}].IsCurved()) if((*mesh.coarsemesh)[ElementIndex{hpel.coarse_elnr}].IsCurved())
@ -1451,7 +1451,7 @@ namespace netgen
for(ElementIndex i=0;i<mesh.GetNE(); i++) for(ElementIndex i=0;i<mesh.GetNE(); i++)
{ {
// Element el = mesh[i] ; // Element el = mesh[i] ;
HPRefElement & hpel = hpelements[mesh[i].hp_elnr]; HPRefElement & hpel = hpelements[mesh[i].GetHpElnr()];
const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType()); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
double dist[3] = {0,0,0}; double dist[3] = {0,0,0};
int ord_dir[3] = {0,0,0}; int ord_dir[3] = {0,0,0};
@ -1523,7 +1523,7 @@ namespace netgen
for(SurfaceElementIndex i=0;i<mesh.GetNSE(); i++) for(SurfaceElementIndex i=0;i<mesh.GetNSE(); i++)
{ {
// Element2d el = mesh[i] ; // Element2d el = mesh[i] ;
HPRefElement & hpel = hpelements[mesh[i].hp_elnr]; HPRefElement & hpel = hpelements[mesh[i].GetHpElnr()];
const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType()); const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
double dist[3] = {0,0,0}; double dist[3] = {0,0,0};
int ord_dir[3] = {0,0,0}; int ord_dir[3] = {0,0,0};

2
libsrc/meshing/improve2.hpp
View File

@ -19,7 +19,7 @@ inline void AppendEdges( const Element & elem, PointIndex pi, Array<std::tuple<P
{ { 0, 1 }, { 0, 2 }, { 0, 3 }, { { 0, 1 }, { 0, 2 }, { 0, 3 },
{ 1, 2 }, { 1, 3 }, { 2, 3 } }; { 1, 2 }, { 1, 3 }, { 2, 3 } };
if(elem.flags.fixed) if(elem.Flags().fixed)
return; return;
for (int j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {

96
libsrc/meshing/improve3.cpp
View File

@ -41,7 +41,7 @@ static ArrayMem<Element, 3> SplitElement (Element old, PointIndex pi0, PointInde
ArrayMem<Element, 3> new_elements; ArrayMem<Element, 3> new_elements;
// split element by cutting edge pi0,pi1 at pinew // split element by cutting edge pi0,pi1 at pinew
auto np = old.GetNP(); auto np = old.GetNP();
old.flags.illegal_valid = 0; old.Flags().illegal_valid = 0;
if(np == 4) if(np == 4)
{ {
// Split tet into two tets // Split tet into two tets
@ -232,7 +232,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
break; break;
} }
elem.flags.illegal_valid = 0; elem.Flags().illegal_valid = 0;
if (!mesh.LegalTet(elem)) if (!mesh.LegalTet(elem))
badness_new += 1e4; badness_new += 1e4;
} }
@ -255,7 +255,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
if (elem[l] == pi1) if (elem[l] == pi1)
elem[l] = pi0; elem[l] = pi0;
elem.flags.illegal_valid = 0; elem.Flags().illegal_valid = 0;
if (!mesh.LegalTet (elem)) if (!mesh.LegalTet (elem))
(*testout) << "illegal tet " << ei << endl; (*testout) << "illegal tet " << ei << endl;
} }
@ -265,7 +265,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
for (auto ei : has_both_points) for (auto ei : has_both_points)
{ {
mesh[ei].flags.illegal_valid = 0; mesh[ei].Flags().illegal_valid = 0;
mesh[ei].Delete(); mesh[ei].Delete();
} }
} }
@ -505,9 +505,9 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
Element newel1 = oldel; Element newel1 = oldel;
Element newel2 = oldel; Element newel2 = oldel;
oldel.flags.illegal_valid = 0; oldel.Flags().illegal_valid = 0;
newel1.flags.illegal_valid = 0; newel1.Flags().illegal_valid = 0;
newel2.flags.illegal_valid = 0; newel2.Flags().illegal_valid = 0;
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
{ {
@ -536,11 +536,11 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
Element newel1 = oldel; Element newel1 = oldel;
Element newel2 = oldel; Element newel2 = oldel;
oldel.flags.illegal_valid = 0; oldel.Flags().illegal_valid = 0;
oldel.Delete(); oldel.Delete();
newel1.flags.illegal_valid = 0; newel1.Flags().illegal_valid = 0;
newel2.flags.illegal_valid = 0; newel2.Flags().illegal_valid = 0;
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
{ {
@ -799,9 +799,9 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el32, 0) + CalcBad (mesh.Points(), el32, 0) +
CalcBad (mesh.Points(), el33, 0); CalcBad (mesh.Points(), el33, 0);
el31.flags.illegal_valid = 0; el31.Flags().illegal_valid = 0;
el32.flags.illegal_valid = 0; el32.Flags().illegal_valid = 0;
el33.flags.illegal_valid = 0; el33.Flags().illegal_valid = 0;
if (!mesh.LegalTet(el31) || if (!mesh.LegalTet(el31) ||
!mesh.LegalTet(el32) || !mesh.LegalTet(el32) ||
@ -823,8 +823,8 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
bad2 = CalcBad (mesh.Points(), el21, 0) + bad2 = CalcBad (mesh.Points(), el21, 0) +
CalcBad (mesh.Points(), el22, 0); CalcBad (mesh.Points(), el22, 0);
el21.flags.illegal_valid = 0; el21.Flags().illegal_valid = 0;
el22.flags.illegal_valid = 0; el22.Flags().illegal_valid = 0;
if (!mesh.LegalTet(el21) || if (!mesh.LegalTet(el21) ||
!mesh.LegalTet(el22)) !mesh.LegalTet(el22))
@ -866,8 +866,8 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
mesh[hasbothpoints[1]].Delete(); mesh[hasbothpoints[1]].Delete();
mesh[hasbothpoints[2]].Delete(); mesh[hasbothpoints[2]].Delete();
el21.flags.illegal_valid = 0; el21.Flags().illegal_valid = 0;
el22.flags.illegal_valid = 0; el22.Flags().illegal_valid = 0;
mesh.AddVolumeElement(el21); mesh.AddVolumeElement(el21);
mesh.AddVolumeElement(el22); mesh.AddVolumeElement(el22);
} }
@ -942,10 +942,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el4, 0); CalcBad (mesh.Points(), el4, 0);
el1.flags.illegal_valid = 0; el1.Flags().illegal_valid = 0;
el2.flags.illegal_valid = 0; el2.Flags().illegal_valid = 0;
el3.flags.illegal_valid = 0; el3.Flags().illegal_valid = 0;
el4.flags.illegal_valid = 0; el4.Flags().illegal_valid = 0;
if (goal != OPT_CONFORM) if (goal != OPT_CONFORM)
@ -978,10 +978,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el3, 0) + CalcBad (mesh.Points(), el3, 0) +
CalcBad (mesh.Points(), el4, 0); CalcBad (mesh.Points(), el4, 0);
el1.flags.illegal_valid = 0; el1.Flags().illegal_valid = 0;
el2.flags.illegal_valid = 0; el2.Flags().illegal_valid = 0;
el3.flags.illegal_valid = 0; el3.Flags().illegal_valid = 0;
el4.flags.illegal_valid = 0; el4.Flags().illegal_valid = 0;
if (goal != OPT_CONFORM) if (goal != OPT_CONFORM)
{ {
@ -1014,10 +1014,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el3b, 0) + CalcBad (mesh.Points(), el3b, 0) +
CalcBad (mesh.Points(), el4b, 0); CalcBad (mesh.Points(), el4b, 0);
el1b.flags.illegal_valid = 0; el1b.Flags().illegal_valid = 0;
el2b.flags.illegal_valid = 0; el2b.Flags().illegal_valid = 0;
el3b.flags.illegal_valid = 0; el3b.Flags().illegal_valid = 0;
el4b.flags.illegal_valid = 0; el4b.Flags().illegal_valid = 0;
if (goal != OPT_CONFORM) if (goal != OPT_CONFORM)
{ {
@ -1054,10 +1054,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
for (auto i : IntRange(4)) for (auto i : IntRange(4))
mesh[hasbothpoints[i]].Delete(); mesh[hasbothpoints[i]].Delete();
el1.flags.illegal_valid = 0; el1.Flags().illegal_valid = 0;
el2.flags.illegal_valid = 0; el2.Flags().illegal_valid = 0;
el3.flags.illegal_valid = 0; el3.Flags().illegal_valid = 0;
el4.flags.illegal_valid = 0; el4.Flags().illegal_valid = 0;
mesh.AddVolumeElement (el1); mesh.AddVolumeElement (el1);
mesh.AddVolumeElement (el2); mesh.AddVolumeElement (el2);
mesh.AddVolumeElement (el3); mesh.AddVolumeElement (el3);
@ -1068,10 +1068,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
for (auto i : IntRange(4)) for (auto i : IntRange(4))
mesh[hasbothpoints[i]].Delete(); mesh[hasbothpoints[i]].Delete();
el1b.flags.illegal_valid = 0; el1b.Flags().illegal_valid = 0;
el2b.flags.illegal_valid = 0; el2b.Flags().illegal_valid = 0;
el3b.flags.illegal_valid = 0; el3b.Flags().illegal_valid = 0;
el4b.flags.illegal_valid = 0; el4b.Flags().illegal_valid = 0;
mesh.AddVolumeElement (el1b); mesh.AddVolumeElement (el1b);
mesh.AddVolumeElement (el2b); mesh.AddVolumeElement (el2b);
mesh.AddVolumeElement (el3b); mesh.AddVolumeElement (el3b);
@ -1174,7 +1174,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
hel[3] = pi2; hel[3] = pi2;
bad2 += CalcBad (mesh.Points(), hel, 0); bad2 += CalcBad (mesh.Points(), hel, 0);
hel.flags.illegal_valid = 0; hel.Flags().illegal_valid = 0;
if (!mesh.LegalTet(hel)) bad2 += 1e4; if (!mesh.LegalTet(hel)) bad2 += 1e4;
hel[2] = suroundpts[k % nsuround]; hel[2] = suroundpts[k % nsuround];
@ -1183,7 +1183,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
bad2 += CalcBad (mesh.Points(), hel, 0); bad2 += CalcBad (mesh.Points(), hel, 0);
hel.flags.illegal_valid = 0; hel.Flags().illegal_valid = 0;
if (!mesh.LegalTet(hel)) bad2 += 1e4; if (!mesh.LegalTet(hel)) bad2 += 1e4;
} }
// (*testout) << "bad2," << l << " = " << bad2 << endl; // (*testout) << "bad2," << l << " = " << bad2 << endl;
@ -1253,7 +1253,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
hel[1] = suroundpts[k % nsuround]; hel[1] = suroundpts[k % nsuround];
hel[2] = suroundpts[(k+1) % nsuround]; hel[2] = suroundpts[(k+1) % nsuround];
hel[3] = pi2; hel[3] = pi2;
hel.flags.illegal_valid = 0; hel.Flags().illegal_valid = 0;
/* /*
(*testout) << nsuround << "-swap, new el,top = " (*testout) << nsuround << "-swap, new el,top = "
@ -1309,7 +1309,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
for (ElementIndex eli : myrange) for (ElementIndex eli : myrange)
{ {
const auto & el = mesh[eli]; const auto & el = mesh[eli];
if(el.flags.fixed) if(el.Flags().fixed)
continue; continue;
for (auto pi : el.PNums()) for (auto pi : el.PNums())
@ -2412,9 +2412,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
CalcBad (mesh.Points(), el33, 0); CalcBad (mesh.Points(), el33, 0);
el31.flags.illegal_valid = 0; el31.Flags().illegal_valid = 0;
el32.flags.illegal_valid = 0; el32.Flags().illegal_valid = 0;
el33.flags.illegal_valid = 0; el33.Flags().illegal_valid = 0;
if (!mesh.LegalTet(el31) || if (!mesh.LegalTet(el31) ||
!mesh.LegalTet(el32) || !mesh.LegalTet(el32) ||
@ -2433,9 +2433,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
if (d_badness<0.0) if (d_badness<0.0)
{ {
el31.flags.illegal_valid = 0; el31.Flags().illegal_valid = 0;
el32.flags.illegal_valid = 0; el32.Flags().illegal_valid = 0;
el33.flags.illegal_valid = 0; el33.Flags().illegal_valid = 0;
mesh[eli1].Delete(); mesh[eli1].Delete();
mesh[eli2].Delete(); mesh[eli2].Delete();
@ -2655,7 +2655,7 @@ double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh,
if(badness_after<badness_before) if(badness_after<badness_before)
{ {
PointIndex pinew = mesh.AddPoint (center); PointIndex pinew = mesh.AddPoint (center);
el.flags.illegal_valid = 0; el.Flags().illegal_valid = 0;
el.Delete(); el.Delete();
for (auto ei1 : has_both_points0) for (auto ei1 : has_both_points0)

22
libsrc/meshing/meshclass.cpp
View File

@ -599,9 +599,9 @@ namespace netgen
{ {
volelements.Append (el); volelements.Append (el);
} }
volelements.Last().flags.illegal_valid = 0; volelements.Last().Flags().illegal_valid = 0;
volelements.Last().flags.fixed = 0; volelements.Last().Flags().fixed = 0;
volelements.Last().flags.deleted = 0; volelements.Last().Flags().deleted = 0;
// while (volelements.Size() > eltyps.Size()) // while (volelements.Size() > eltyps.Size())
// eltyps.Append (FREEELEMENT); // eltyps.Append (FREEELEMENT);
@ -622,9 +622,9 @@ namespace netgen
*/ */
volelements[ei] = el; volelements[ei] = el;
volelements[ei].flags.illegal_valid = 0; volelements[ei].Flags().illegal_valid = 0;
volelements[ei].flags.fixed = 0; volelements[ei].Flags().fixed = 0;
volelements[ei].flags.deleted = 0; volelements[ei].Flags().deleted = 0;
} }
@ -3241,7 +3241,7 @@ namespace netgen
if (dist[el[j]] < elmin) if (dist[el[j]] < elmin)
elmin = dist[el[j]]; elmin = dist[el[j]];
el.flags.fixed = elmin > layers; el.Flags().fixed = elmin > layers;
// eltyps.Elem(i) = (elmin <= layers) ? // eltyps.Elem(i) = (elmin <= layers) ?
// FREEELEMENT : FIXEDELEMENT; // FREEELEMENT : FIXEDELEMENT;
if (elmin <= layers) if (elmin <= layers)
@ -4384,7 +4384,7 @@ namespace netgen
for (i = 1; i <= ne; i++) for (i = 1; i <= ne; i++)
{ {
Element & el = (Element&) VolumeElement(i); Element & el = (Element&) VolumeElement(i);
el.flags.badel = 0; el.Flags().badel = 0;
int nip = el.GetNIP(); int nip = el.GetNIP();
for (j = 1; j <= nip; j++) for (j = 1; j <= nip; j++)
{ {
@ -4393,7 +4393,7 @@ namespace netgen
if (det > 0) if (det > 0)
{ {
PrintError ("Element ", i , " has wrong orientation"); PrintError ("Element ", i , " has wrong orientation");
el.flags.badel = 1; el.Flags().badel = 1;
} }
} }
} }
@ -6472,7 +6472,7 @@ namespace netgen
if (el.GetType() != TET) if (el.GetType() != TET)
{ {
VolumeElement(i).flags.badel = 0; VolumeElement(i).Flags().badel = 0;
continue; continue;
} }
@ -6554,7 +6554,7 @@ namespace netgen
} }
VolumeElement(i).flags.badel = badel; VolumeElement(i).Flags().badel = badel;
if (badel) badtets++; if (badel) badtets++;
} }

2
libsrc/meshing/meshclass.hpp
View File

@ -365,7 +365,7 @@ namespace netgen
Element & operator[] (ElementIndex ei) { return volelements[ei]; } Element & operator[] (ElementIndex ei) { return volelements[ei]; }
ELEMENTTYPE ElementType (ElementIndex i) const ELEMENTTYPE ElementType (ElementIndex i) const
{ return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; } { return (volelements[i].Flags().fixed) ? FIXEDELEMENT : FREEELEMENT; }
const auto & VolumeElements() const { return volelements; } const auto & VolumeElements() const { return volelements; }
auto & VolumeElements() { return volelements; } auto & VolumeElements() { return volelements; }

56
libsrc/meshing/meshtype.cpp
View File

@ -178,62 +178,6 @@ namespace netgen
*/ */
} }
Segment::Segment (const Segment & other)
:
edgenr(other.edgenr),
singedge_left(other.singedge_left),
singedge_right(other.singedge_right),
seginfo(other.seginfo),
si(other.si),
domin(other.domin),
domout(other.domout),
tlosurf(other.tlosurf),
geominfo(),
surfnr1(other.surfnr1),
surfnr2(other.surfnr2),
epgeominfo(),
meshdocval(other.meshdocval),
is_curved(other.is_curved),
hp_elnr(other.hp_elnr)
{
for (int j = 0; j < 3; j++)
pnums[j] = other.pnums[j];
geominfo[0] = other.geominfo[0];
geominfo[1] = other.geominfo[1];
epgeominfo[0] = other.epgeominfo[0];
epgeominfo[1] = other.epgeominfo[1];
}
Segment& Segment::operator=(const Segment & other)
{
if (&other != this)
{
pnums[0] = other[0];
pnums[1] = other[1];
edgenr = other.edgenr;
singedge_left = other.singedge_left;
singedge_right = other.singedge_right;
seginfo = other.seginfo;
si = other.si;
domin = other.domin;
domout = other.domout;
tlosurf = other.tlosurf;
geominfo[0] = other.geominfo[0];
geominfo[1] = other.geominfo[1];
surfnr1 = other.surfnr1;
surfnr2 = other.surfnr2;
epgeominfo[0] = other.epgeominfo[0];
epgeominfo[1] = other.epgeominfo[1];
pnums[2] = other.pnums[2];
meshdocval = other.meshdocval;
hp_elnr = other.hp_elnr;
is_curved = other.is_curved;
}
return *this;
}
void Segment :: DoArchive (Archive & ar) void Segment :: DoArchive (Archive & ar)
{ {
string * bcname_dummy = nullptr; string * bcname_dummy = nullptr;

86
libsrc/meshing/meshtype.hpp
View File

@ -128,15 +128,7 @@ namespace netgen
EdgePointGeomInfo () EdgePointGeomInfo ()
: edgenr(-1), body(0), dist(0.0), u(0.0), v(0.0) { ; } : edgenr(-1), body(0), dist(0.0), u(0.0), v(0.0) { ; }
EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) = default;
EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2)
{
edgenr = gi2.edgenr;
body = gi2.body;
dist = gi2.dist;
u = gi2.u; v = gi2.v;
return *this;
}
}; };
inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi) inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi)
@ -430,8 +422,19 @@ namespace netgen
/// a linked list for all segments in the same face /// a linked list for all segments in the same face
SurfaceElementIndex next; SurfaceElementIndex next;
///
int hp_elnr;
public: public:
static auto GetDataLayout()
{
return std::map<string, int>({
{ "pnum", offsetof(Element2d, pnum)},
{ "index", offsetof(Element2d, index) },
{ "np", offsetof(Element2d, np) }
});
}
/// ///
DLL_HEADER Element2d (); DLL_HEADER Element2d ();
Element2d (const Element2d &) = default; Element2d (const Element2d &) = default;
@ -587,6 +590,8 @@ namespace netgen
void SetOrder (int ox, int oy, int /* oz */) { orderx = ox; ordery = oy;} void SetOrder (int ox, int oy, int /* oz */) { orderx = ox; ordery = oy;}
void SetOrder (int ox, int oy) { orderx = ox; ordery = oy;} void SetOrder (int ox, int oy) { orderx = ox; ordery = oy;}
int GetHpElnr() const { return hp_elnr; }
void SetHpElnr(int _hp_elnr) { hp_elnr = _hp_elnr; }
/// ///
void GetBox (const T_POINTS & points, Box3d & box) const; void GetBox (const T_POINTS & points, Box3d & box) const;
@ -679,10 +684,6 @@ namespace netgen
bool operator==(const Element2d & el2) const; bool operator==(const Element2d & el2) const;
int HasFace(const Element2d& el) const; int HasFace(const Element2d& el) const;
///
int meshdocval;
///
int hp_elnr;
}; };
ostream & operator<<(ostream & s, const Element2d & el); ostream & operator<<(ostream & s, const Element2d & el);
@ -719,20 +720,6 @@ namespace netgen
ELEMENT_TYPE typ; ELEMENT_TYPE typ;
/// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism)
int8_t np; int8_t np;
///
class flagstruct {
public:
bool marked:1; // marked for refinement
bool badel:1; // angles worse then limit
bool reverse:1; // for refinement a la Bey
bool illegal:1; // illegal, will be split or swapped
bool illegal_valid:1; // is illegal-flag valid ?
bool badness_valid:1; // is badness valid ?
bool refflag:1; // mark element for refinement
bool strongrefflag:1;
bool deleted:1; // element is deleted, will be removed from array
bool fixed:1; // don't change element in optimization
};
/// sub-domain index /// sub-domain index
int index; int index;
@ -747,8 +734,32 @@ namespace netgen
float badness; float badness;
bool is_curved:1; // element is (high order) curved bool is_curved:1; // element is (high order) curved
public: class flagstruct {
public:
bool marked:1; // marked for refinement
bool badel:1; // angles worse then limit
bool reverse:1; // for refinement a la Bey
bool illegal:1; // illegal, will be split or swapped
bool illegal_valid:1; // is illegal-flag valid ?
bool badness_valid:1; // is badness valid ?
bool refflag:1; // mark element for refinement
bool strongrefflag:1;
bool deleted:1; // element is deleted, will be removed from array
bool fixed:1; // don't change element in optimization
};
flagstruct flags; flagstruct flags;
int hp_elnr;
public:
static auto GetDataLayout()
{
return std::map<string, int>({
{ "pnum", offsetof(Element, pnum)},
{ "index", offsetof(Element, index) },
{ "np", offsetof(Element, np) }
});
}
/// ///
DLL_HEADER Element () = default; DLL_HEADER Element () = default;
@ -763,6 +774,9 @@ namespace netgen
DLL_HEADER Element (ELEMENT_TYPE type); DLL_HEADER Element (ELEMENT_TYPE type);
/// ///
// Element & operator= (const Element & el2); // Element & operator= (const Element & el2);
const flagstruct& Flags() const { return flags; }
flagstruct& Flags() { return flags; }
/// ///
DLL_HEADER void SetNP (int anp); DLL_HEADER void SetNP (int anp);
@ -909,6 +923,9 @@ namespace netgen
/// ///
DLL_HEADER void Invert (); DLL_HEADER void Invert ();
int GetHpElnr() const { return hp_elnr; }
void SetHpElnr(int _hp_elnr) { hp_elnr = _hp_elnr; }
/// split into 4 node tets /// split into 4 node tets
void GetTets (NgArray<Element> & locels) const; void GetTets (NgArray<Element> & locels) const;
/// split into 4 node tets, local point nrs /// split into 4 node tets, local point nrs
@ -993,7 +1010,6 @@ namespace netgen
bool IsCurved () const { return is_curved; } bool IsCurved () const { return is_curved; }
void SetCurved (bool acurved) { is_curved = acurved; } void SetCurved (bool acurved) { is_curved = acurved; }
int hp_elnr;
}; };
ostream & operator<<(ostream & s, const Element & el); ostream & operator<<(ostream & s, const Element & el);
@ -1011,10 +1027,7 @@ namespace netgen
public: public:
/// ///
DLL_HEADER Segment(); DLL_HEADER Segment();
DLL_HEADER Segment (const Segment& other); Segment (const Segment& other) = default;
~Segment()
{ ; }
// friend ostream & operator<<(ostream & s, const Segment & seg); // friend ostream & operator<<(ostream & s, const Segment & seg);
@ -1050,10 +1063,8 @@ namespace netgen
/// ///
int meshdocval; int meshdocval;
private:
bool is_curved; bool is_curved;
int hp_elnr;
public:
/* /*
PointIndex operator[] (int i) const PointIndex operator[] (int i) const
{ return (i == 0) ? p1 : p2; } { return (i == 0) ? p1 : p2; }
@ -1062,10 +1073,9 @@ namespace netgen
{ return (i == 0) ? p1 : p2; } { return (i == 0) ? p1 : p2; }
*/ */
Segment& operator=(const Segment & other); Segment& operator=(const Segment & other) = default;
int hp_elnr;
int GetNP() const int GetNP() const
{ {

55
libsrc/meshing/python_mesh.cpp
View File

@ -434,6 +434,27 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
})) }))
; ;
if(ngcore_have_numpy)
{
auto data_layout = Element::GetDataLayout();
py::detail::npy_format_descriptor<Element>::register_dtype({
py::detail::field_descriptor {
"nodes", data_layout["pnum"],
ELEMENT_MAXPOINTS * sizeof(PointIndex),
py::format_descriptor<int[ELEMENT_MAXPOINTS]>::format(),
py::detail::npy_format_descriptor<int[ELEMENT_MAXPOINTS]>::dtype() },
py::detail::field_descriptor {
"index", data_layout["index"], sizeof(int),
py::format_descriptor<int>::format(),
py::detail::npy_format_descriptor<int>::dtype() },
py::detail::field_descriptor {
"np", data_layout["np"], sizeof(int8_t),
py::format_descriptor<signed char>::format(),
pybind11::dtype("int8") }
});
}
py::class_<Element2d>(m, "Element2D") py::class_<Element2d>(m, "Element2D")
.def(py::init ([](int index, std::vector<PointIndex> vertices) .def(py::init ([](int index, std::vector<PointIndex> vertices)
{ {
@ -493,6 +514,26 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
})) }))
; ;
if(ngcore_have_numpy)
{
auto data_layout = Element2d::GetDataLayout();
py::detail::npy_format_descriptor<Element2d>::register_dtype({
py::detail::field_descriptor {
"nodes", data_layout["pnum"],
ELEMENT2D_MAXPOINTS * sizeof(PointIndex),
py::format_descriptor<int[ELEMENT2D_MAXPOINTS]>::format(),
py::detail::npy_format_descriptor<int[ELEMENT2D_MAXPOINTS]>::dtype() },
py::detail::field_descriptor {
"index", data_layout["index"], sizeof(int),
py::format_descriptor<int>::format(),
py::detail::npy_format_descriptor<int>::dtype() },
py::detail::field_descriptor {
"np", data_layout["np"], sizeof(int8_t),
py::format_descriptor<signed char>::format(),
pybind11::dtype("int8") }
});
}
py::class_<Segment>(m, "Element1D") py::class_<Segment>(m, "Element1D")
.def(py::init([](py::list vertices, py::list surfaces, int index, int edgenr, .def(py::init([](py::list vertices, py::list surfaces, int index, int edgenr,
py::list trignums) py::list trignums)
@ -557,6 +598,20 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
})) }))
; ;
if(ngcore_have_numpy)
{
py::detail::npy_format_descriptor<Segment>::register_dtype({
py::detail::field_descriptor {
"nodes", offsetof(Segment, pnums),
3 * sizeof(PointIndex),
py::format_descriptor<int[3]>::format(),
py::detail::npy_format_descriptor<int[3]>::dtype() },
py::detail::field_descriptor {
"index", offsetof(Segment, edgenr), sizeof(int),
py::format_descriptor<int>::format(),
py::detail::npy_format_descriptor<int>::dtype() },
});
}
py::class_<Element0d>(m, "Element0D") py::class_<Element0d>(m, "Element0D")
.def(py::init([](PointIndex vertex, int index) .def(py::init([](PointIndex vertex, int index)

14
libsrc/meshing/refine.cpp
View File

@ -414,7 +414,7 @@ namespace netgen
{ 2, 4, 9 }, { 2, 4, 9 },
{ 3, 4, 10 } }; { 3, 4, 10 } };
int elrev = el.flags.reverse; int elrev = el.Flags().reverse;
for (int j = 1; j <= 4; j++) for (int j = 1; j <= 4; j++)
pnums.Elem(j) = el.PNum(j); pnums.Elem(j) = el.PNum(j);
@ -480,10 +480,10 @@ namespace netgen
for (int k = 1; k <= 4; k++) for (int k = 1; k <= 4; k++)
nel.PNum(k) = pnums.Get(reftab[j][k-1]); nel.PNum(k) = pnums.Get(reftab[j][k-1]);
nel.SetIndex(ind); nel.SetIndex(ind);
nel.flags.reverse = reverse[j]; nel.Flags().reverse = reverse[j];
if (elrev) if (elrev)
{ {
nel.flags.reverse = !nel.flags.reverse; nel.Flags().reverse = !nel.Flags().reverse;
swap (nel.PNum(3), nel.PNum(4)); swap (nel.PNum(3), nel.PNum(4));
} }
@ -794,10 +794,10 @@ namespace netgen
if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0) if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
{ {
wrongels++; wrongels++;
mesh.VolumeElement(i).flags.badel = 1; mesh.VolumeElement(i).Flags().badel = 1;
} }
else else
mesh.VolumeElement(i).flags.badel = 0; mesh.VolumeElement(i).Flags().badel = 0;
if (wrongels) if (wrongels)
{ {
@ -900,14 +900,14 @@ namespace netgen
if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0) if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
{ {
wrongels++; wrongels++;
mesh.VolumeElement(i).flags.badel = 1; mesh.VolumeElement(i).Flags().badel = 1;
(*testout) << "wrong el: "; (*testout) << "wrong el: ";
for (int j = 1; j <= 4; j++) for (int j = 1; j <= 4; j++)
(*testout) << mesh.VolumeElement(i).PNum(j) << " "; (*testout) << mesh.VolumeElement(i).PNum(j) << " ";
(*testout) << endl; (*testout) << endl;
} }
else else
mesh.VolumeElement(i).flags.badel = 0; mesh.VolumeElement(i).Flags().badel = 0;
} }
cout << "wrongels = " << wrongels << endl; cout << "wrongels = " << wrongels << endl;
} }

12
libsrc/meshing/secondorder.cpp
View File

@ -473,10 +473,10 @@ namespace netgen
if (mesh.VolumeElement(i).CalcJacobianBadness (mesh.Points()) > 1e10) if (mesh.VolumeElement(i).CalcJacobianBadness (mesh.Points()) > 1e10)
{ {
wrongels++; wrongels++;
mesh.VolumeElement(i).flags.badel = 1; mesh.VolumeElement(i).Flags().badel = 1;
} }
else else
mesh.VolumeElement(i).flags.badel = 0; mesh.VolumeElement(i).Flags().badel = 0;
double facok = 0; double facok = 0;
double factry; double factry;
@ -559,7 +559,7 @@ namespace netgen
{ {
wrongels++; wrongels++;
Element & el = mesh.VolumeElement(i); Element & el = mesh.VolumeElement(i);
el.flags.badel = 1; el.Flags().badel = 1;
if (lam < 1e-4) if (lam < 1e-4)
@ -574,7 +574,7 @@ namespace netgen
*/ */
} }
else else
mesh.VolumeElement(i).flags.badel = 0; mesh.VolumeElement(i).Flags().badel = 0;
} }
cout << "wrongels = " << wrongels << endl; cout << "wrongels = " << wrongels << endl;
} }
@ -597,10 +597,10 @@ namespace netgen
if (illegalels.Test(i)) if (illegalels.Test(i))
{ {
cout << "illegal element: " << i << endl; cout << "illegal element: " << i << endl;
mesh.VolumeElement(i).flags.badel = 1; mesh.VolumeElement(i).Flags().badel = 1;
} }
else else
mesh.VolumeElement(i).flags.badel = 0; mesh.VolumeElement(i).Flags().badel = 0;
} }
/* /*

2
libsrc/meshing/topology.cpp
View File

@ -1320,7 +1320,7 @@ namespace netgen
if (mesh->coarsemesh && mesh->hpelements->Size() == mesh->GetNE() ) if (mesh->coarsemesh && mesh->hpelements->Size() == mesh->GetNE() )
{ {
const HPRefElement & hpref_el = const HPRefElement & hpref_el =
(*mesh->hpelements) [ (*mesh)[vertels[k]].hp_elnr]; (*mesh->hpelements) [ (*mesh)[vertels[k]].GetHpElnr()];
(*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl; (*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl;
} }

6
libsrc/visualization/vsmesh.cpp
View File

@ -597,8 +597,8 @@ namespace netgen
for (int i = 1; i <= mesh->GetNE(); i++) for (int i = 1; i <= mesh->GetNE(); i++)
{ {
if (mesh->VolumeElement(i).flags.badel || if (mesh->VolumeElement(i).Flags().badel ||
mesh->VolumeElement(i).flags.illegal || mesh->VolumeElement(i).Flags().illegal ||
(i == vispar.drawelement)) (i == vispar.drawelement))
{ {
// copy to be thread-safe // copy to be thread-safe
@ -636,7 +636,7 @@ namespace netgen
for (ElementIndex ei : mesh->VolumeElements().Range()) for (ElementIndex ei : mesh->VolumeElements().Range())
{ {
if ((*mesh)[ei].flags.badel) if ((*mesh)[ei].Flags().badel)
{ {
// copy to be thread-safe // copy to be thread-safe
Element el = (*mesh)[ei]; Element el = (*mesh)[ei];

16
tests/pytest/test_array.py
View File

@ -13,5 +13,21 @@ def test_array_numpy():
a.NumPy().sort() a.NumPy().sort()
assert(all(a == array([0,0,2,2,5]))) assert(all(a == array([0,0,2,2,5])))
def test_mesh_elements_numpy_array_access():
from netgen.csg import unit_cube
mesh = unit_cube.GenerateMesh()
np_els = mesh.Elements3D().NumPy()
vol_nodes = np_els["nodes"]
indices = np_els["index"]
nps = np_els["np"]
for nodes, el, index, np in zip(vol_nodes, mesh.Elements3D(), indices, nps):
for n1, n2 in zip(nodes, el.vertices):
assert n1 == n2
for n in nodes[len(el.vertices):]:
assert n == 0
assert el.index == index
assert len(el.vertices) == np
if __name__ == "__main__": if __name__ == "__main__":
test_array_numpy() test_array_numpy()
test_mesh_elements_numpy_array_access()