From c18a317702569e09f78108e124b8e73623647b73 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 13 Sep 2022 15:12:42 +0200 Subject: [PATCH] register 1,2,3d elements to numpy to be used in arrays --- libsrc/interface/nginterface_v2.cpp | 4 +- libsrc/meshing/curvedelems.cpp | 14 ++--- libsrc/meshing/hprefinement.cpp | 8 +-- libsrc/meshing/improve2.hpp | 2 +- libsrc/meshing/improve3.cpp | 96 ++++++++++++++--------------- libsrc/meshing/meshclass.cpp | 22 +++---- libsrc/meshing/meshclass.hpp | 2 +- libsrc/meshing/meshtype.cpp | 56 ----------------- libsrc/meshing/meshtype.hpp | 86 ++++++++++++++------------ libsrc/meshing/python_mesh.cpp | 55 +++++++++++++++++ libsrc/meshing/refine.cpp | 14 ++--- libsrc/meshing/secondorder.cpp | 12 ++-- libsrc/meshing/topology.cpp | 2 +- libsrc/visualization/vsmesh.cpp | 6 +- tests/pytest/test_array.py | 16 +++++ 15 files changed, 210 insertions(+), 185 deletions(-) diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index 3f7b001e..e6f3bce2 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -701,9 +701,9 @@ namespace netgen { int hpelnr = -1; if (mesh->GetDimension() == 2) - hpelnr = mesh->SurfaceElement(ei).hp_elnr; + hpelnr = mesh->SurfaceElement(ei).GetHpElnr(); else - hpelnr = mesh->VolumeElement(ei).hp_elnr; + hpelnr = mesh->VolumeElement(ei).GetHpElnr(); if (hpelnr < 0) throw NgException("Ngx_Mesh::GetHPElementLevel: Wrong hp-element number!"); diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index e32a9eb6..50d3fa35 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -1632,7 +1632,7 @@ namespace netgen if (mesh.coarsemesh) { const HPRefElement & hpref_el = - (*mesh.hpelements) [mesh[elnr].hp_elnr]; + (*mesh.hpelements) [mesh[elnr].GetHpElnr()]; return mesh.coarsemesh->GetCurvedElements().IsSurfaceElementCurved (hpref_el.coarse_elnr); } @@ -1679,7 +1679,7 @@ namespace netgen if (mesh.coarsemesh) { const HPRefElement & hpref_el = - (*mesh.hpelements) [mesh[elnr].hp_elnr]; + (*mesh.hpelements) [mesh[elnr].GetHpElnr()]; // xi umrechnen double lami[4]; @@ -2428,7 +2428,7 @@ namespace netgen if (mesh.coarsemesh) { const HPRefElement & hpref_el = - (*mesh.hpelements) [mesh[elnr].hp_elnr]; + (*mesh.hpelements) [mesh[elnr].GetHpElnr()]; return mesh.coarsemesh->GetCurvedElements().IsElementCurved (hpref_el.coarse_elnr); } @@ -2480,7 +2480,7 @@ namespace netgen if (mesh.coarsemesh) { const HPRefElement & hpref_el = - (*mesh.hpelements) [mesh[elnr].hp_elnr]; + (*mesh.hpelements) [mesh[elnr].GetHpElnr()]; return mesh.coarsemesh->GetCurvedElements().IsElementHighOrder (hpref_el.coarse_elnr); } @@ -2519,7 +2519,7 @@ namespace netgen if (mesh.coarsemesh) { const HPRefElement & hpref_el = - (*mesh.hpelements) [mesh[elnr].hp_elnr]; + (*mesh.hpelements) [mesh[elnr].GetHpElnr()]; // xi umrechnen double lami[8]; @@ -4065,7 +4065,7 @@ namespace netgen if (mesh.coarsemesh) { const HPRefElement & hpref_el = - (*mesh.hpelements) [mesh[elnr].hp_elnr]; + (*mesh.hpelements) [mesh[elnr].GetHpElnr()]; // xi umrechnen T lami[4]; @@ -4529,7 +4529,7 @@ namespace netgen if (mesh.coarsemesh) { const HPRefElement & hpref_el = - (*mesh.hpelements) [mesh[elnr].hp_elnr]; + (*mesh.hpelements) [mesh[elnr].GetHpElnr()]; // xi umrechnen T lami[8]; diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index a188ca79..d6819c44 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -1403,7 +1403,7 @@ namespace netgen Element2d el(hpel.np); for(int j=0;j SplitElement (Element old, PointIndex pi0, PointInde ArrayMem new_elements; // split element by cutting edge pi0,pi1 at pinew auto np = old.GetNP(); - old.flags.illegal_valid = 0; + old.Flags().illegal_valid = 0; if(np == 4) { // Split tet into two tets @@ -232,7 +232,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh, break; } - elem.flags.illegal_valid = 0; + elem.Flags().illegal_valid = 0; if (!mesh.LegalTet(elem)) badness_new += 1e4; } @@ -255,7 +255,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh, if (elem[l] == pi1) elem[l] = pi0; - elem.flags.illegal_valid = 0; + elem.Flags().illegal_valid = 0; if (!mesh.LegalTet (elem)) (*testout) << "illegal tet " << ei << endl; } @@ -265,7 +265,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh, for (auto ei : has_both_points) { - mesh[ei].flags.illegal_valid = 0; + mesh[ei].Flags().illegal_valid = 0; mesh[ei].Delete(); } } @@ -505,9 +505,9 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table Element newel1 = oldel; Element newel2 = oldel; - oldel.flags.illegal_valid = 0; - newel1.flags.illegal_valid = 0; - newel2.flags.illegal_valid = 0; + oldel.Flags().illegal_valid = 0; + newel1.Flags().illegal_valid = 0; + newel2.Flags().illegal_valid = 0; for (int l = 0; l < 4; l++) { @@ -536,11 +536,11 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table Element newel1 = oldel; Element newel2 = oldel; - oldel.flags.illegal_valid = 0; + oldel.Flags().illegal_valid = 0; oldel.Delete(); - newel1.flags.illegal_valid = 0; - newel2.flags.illegal_valid = 0; + newel1.Flags().illegal_valid = 0; + newel2.Flags().illegal_valid = 0; 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(), el33, 0); - el31.flags.illegal_valid = 0; - el32.flags.illegal_valid = 0; - el33.flags.illegal_valid = 0; + el31.Flags().illegal_valid = 0; + el32.Flags().illegal_valid = 0; + el33.Flags().illegal_valid = 0; if (!mesh.LegalTet(el31) || !mesh.LegalTet(el32) || @@ -823,8 +823,8 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, bad2 = CalcBad (mesh.Points(), el21, 0) + CalcBad (mesh.Points(), el22, 0); - el21.flags.illegal_valid = 0; - el22.flags.illegal_valid = 0; + el21.Flags().illegal_valid = 0; + el22.Flags().illegal_valid = 0; if (!mesh.LegalTet(el21) || !mesh.LegalTet(el22)) @@ -866,8 +866,8 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, mesh[hasbothpoints[1]].Delete(); mesh[hasbothpoints[2]].Delete(); - el21.flags.illegal_valid = 0; - el22.flags.illegal_valid = 0; + el21.Flags().illegal_valid = 0; + el22.Flags().illegal_valid = 0; mesh.AddVolumeElement(el21); mesh.AddVolumeElement(el22); } @@ -942,10 +942,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, CalcBad (mesh.Points(), el4, 0); - el1.flags.illegal_valid = 0; - el2.flags.illegal_valid = 0; - el3.flags.illegal_valid = 0; - el4.flags.illegal_valid = 0; + el1.Flags().illegal_valid = 0; + el2.Flags().illegal_valid = 0; + el3.Flags().illegal_valid = 0; + el4.Flags().illegal_valid = 0; if (goal != OPT_CONFORM) @@ -978,10 +978,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, CalcBad (mesh.Points(), el3, 0) + CalcBad (mesh.Points(), el4, 0); - el1.flags.illegal_valid = 0; - el2.flags.illegal_valid = 0; - el3.flags.illegal_valid = 0; - el4.flags.illegal_valid = 0; + el1.Flags().illegal_valid = 0; + el2.Flags().illegal_valid = 0; + el3.Flags().illegal_valid = 0; + el4.Flags().illegal_valid = 0; if (goal != OPT_CONFORM) { @@ -1014,10 +1014,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, CalcBad (mesh.Points(), el3b, 0) + CalcBad (mesh.Points(), el4b, 0); - el1b.flags.illegal_valid = 0; - el2b.flags.illegal_valid = 0; - el3b.flags.illegal_valid = 0; - el4b.flags.illegal_valid = 0; + el1b.Flags().illegal_valid = 0; + el2b.Flags().illegal_valid = 0; + el3b.Flags().illegal_valid = 0; + el4b.Flags().illegal_valid = 0; if (goal != OPT_CONFORM) { @@ -1054,10 +1054,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, for (auto i : IntRange(4)) mesh[hasbothpoints[i]].Delete(); - el1.flags.illegal_valid = 0; - el2.flags.illegal_valid = 0; - el3.flags.illegal_valid = 0; - el4.flags.illegal_valid = 0; + el1.Flags().illegal_valid = 0; + el2.Flags().illegal_valid = 0; + el3.Flags().illegal_valid = 0; + el4.Flags().illegal_valid = 0; mesh.AddVolumeElement (el1); mesh.AddVolumeElement (el2); mesh.AddVolumeElement (el3); @@ -1068,10 +1068,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, for (auto i : IntRange(4)) mesh[hasbothpoints[i]].Delete(); - el1b.flags.illegal_valid = 0; - el2b.flags.illegal_valid = 0; - el3b.flags.illegal_valid = 0; - el4b.flags.illegal_valid = 0; + el1b.Flags().illegal_valid = 0; + el2b.Flags().illegal_valid = 0; + el3b.Flags().illegal_valid = 0; + el4b.Flags().illegal_valid = 0; mesh.AddVolumeElement (el1b); mesh.AddVolumeElement (el2b); mesh.AddVolumeElement (el3b); @@ -1174,7 +1174,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, hel[3] = pi2; bad2 += CalcBad (mesh.Points(), hel, 0); - hel.flags.illegal_valid = 0; + hel.Flags().illegal_valid = 0; if (!mesh.LegalTet(hel)) bad2 += 1e4; hel[2] = suroundpts[k % nsuround]; @@ -1183,7 +1183,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, bad2 += CalcBad (mesh.Points(), hel, 0); - hel.flags.illegal_valid = 0; + hel.Flags().illegal_valid = 0; if (!mesh.LegalTet(hel)) bad2 += 1e4; } // (*testout) << "bad2," << l << " = " << bad2 << endl; @@ -1253,7 +1253,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, hel[1] = suroundpts[k % nsuround]; hel[2] = suroundpts[(k+1) % nsuround]; hel[3] = pi2; - hel.flags.illegal_valid = 0; + hel.Flags().illegal_valid = 0; /* (*testout) << nsuround << "-swap, new el,top = " @@ -1309,7 +1309,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, for (ElementIndex eli : myrange) { const auto & el = mesh[eli]; - if(el.flags.fixed) + if(el.Flags().fixed) continue; for (auto pi : el.PNums()) @@ -2412,9 +2412,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI CalcBad (mesh.Points(), el33, 0); - el31.flags.illegal_valid = 0; - el32.flags.illegal_valid = 0; - el33.flags.illegal_valid = 0; + el31.Flags().illegal_valid = 0; + el32.Flags().illegal_valid = 0; + el33.Flags().illegal_valid = 0; if (!mesh.LegalTet(el31) || !mesh.LegalTet(el32) || @@ -2433,9 +2433,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI if (d_badness<0.0) { - el31.flags.illegal_valid = 0; - el32.flags.illegal_valid = 0; - el33.flags.illegal_valid = 0; + el31.Flags().illegal_valid = 0; + el32.Flags().illegal_valid = 0; + el33.Flags().illegal_valid = 0; mesh[eli1].Delete(); mesh[eli2].Delete(); @@ -2655,7 +2655,7 @@ double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh, if(badness_after eltyps.Size()) // eltyps.Append (FREEELEMENT); @@ -622,9 +622,9 @@ namespace netgen */ volelements[ei] = el; - volelements[ei].flags.illegal_valid = 0; - volelements[ei].flags.fixed = 0; - volelements[ei].flags.deleted = 0; + volelements[ei].Flags().illegal_valid = 0; + volelements[ei].Flags().fixed = 0; + volelements[ei].Flags().deleted = 0; } @@ -3241,7 +3241,7 @@ namespace netgen if (dist[el[j]] < elmin) elmin = dist[el[j]]; - el.flags.fixed = elmin > layers; + el.Flags().fixed = elmin > layers; // eltyps.Elem(i) = (elmin <= layers) ? // FREEELEMENT : FIXEDELEMENT; if (elmin <= layers) @@ -4384,7 +4384,7 @@ namespace netgen for (i = 1; i <= ne; i++) { Element & el = (Element&) VolumeElement(i); - el.flags.badel = 0; + el.Flags().badel = 0; int nip = el.GetNIP(); for (j = 1; j <= nip; j++) { @@ -4393,7 +4393,7 @@ namespace netgen if (det > 0) { PrintError ("Element ", i , " has wrong orientation"); - el.flags.badel = 1; + el.Flags().badel = 1; } } } @@ -6472,7 +6472,7 @@ namespace netgen if (el.GetType() != TET) { - VolumeElement(i).flags.badel = 0; + VolumeElement(i).Flags().badel = 0; continue; } @@ -6554,7 +6554,7 @@ namespace netgen } - VolumeElement(i).flags.badel = badel; + VolumeElement(i).Flags().badel = badel; if (badel) badtets++; } diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index c1ea8e50..ff0e6e26 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -365,7 +365,7 @@ namespace netgen Element & operator[] (ElementIndex ei) { return volelements[ei]; } 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; } auto & VolumeElements() { return volelements; } diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 1ea9d220..d49d99cd 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -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) { string * bcname_dummy = nullptr; diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 4ca6c596..6722bde0 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -128,15 +128,7 @@ namespace netgen EdgePointGeomInfo () : edgenr(-1), body(0), dist(0.0), u(0.0), v(0.0) { ; } - - EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) - { - edgenr = gi2.edgenr; - body = gi2.body; - dist = gi2.dist; - u = gi2.u; v = gi2.v; - return *this; - } + EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) = default; }; inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi) @@ -430,8 +422,19 @@ namespace netgen /// a linked list for all segments in the same face SurfaceElementIndex next; + /// + int hp_elnr; public: + static auto GetDataLayout() + { + return std::map({ + { "pnum", offsetof(Element2d, pnum)}, + { "index", offsetof(Element2d, index) }, + { "np", offsetof(Element2d, np) } + }); + } + /// DLL_HEADER Element2d (); 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) { 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; @@ -679,10 +684,6 @@ namespace netgen bool operator==(const Element2d & el2) const; int HasFace(const Element2d& el) const; - /// - int meshdocval; - /// - int hp_elnr; }; ostream & operator<<(ostream & s, const Element2d & el); @@ -719,20 +720,6 @@ namespace netgen ELEMENT_TYPE typ; /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) 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 int index; @@ -747,8 +734,32 @@ namespace netgen float badness; 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; + int hp_elnr; + public: + + static auto GetDataLayout() + { + return std::map({ + { "pnum", offsetof(Element, pnum)}, + { "index", offsetof(Element, index) }, + { "np", offsetof(Element, np) } + }); + } /// DLL_HEADER Element () = default; @@ -763,6 +774,9 @@ namespace netgen DLL_HEADER Element (ELEMENT_TYPE type); /// // Element & operator= (const Element & el2); + + const flagstruct& Flags() const { return flags; } + flagstruct& Flags() { return flags; } /// DLL_HEADER void SetNP (int anp); @@ -909,6 +923,9 @@ namespace netgen /// DLL_HEADER void Invert (); + int GetHpElnr() const { return hp_elnr; } + void SetHpElnr(int _hp_elnr) { hp_elnr = _hp_elnr; } + /// split into 4 node tets void GetTets (NgArray & locels) const; /// split into 4 node tets, local point nrs @@ -993,7 +1010,6 @@ namespace netgen bool IsCurved () const { return is_curved; } void SetCurved (bool acurved) { is_curved = acurved; } - int hp_elnr; }; ostream & operator<<(ostream & s, const Element & el); @@ -1011,10 +1027,7 @@ namespace netgen public: /// DLL_HEADER Segment(); - DLL_HEADER Segment (const Segment& other); - - ~Segment() - { ; } + Segment (const Segment& other) = default; // friend ostream & operator<<(ostream & s, const Segment & seg); @@ -1050,10 +1063,8 @@ namespace netgen /// int meshdocval; - private: bool is_curved; - - public: + int hp_elnr; /* PointIndex operator[] (int i) const { return (i == 0) ? p1 : p2; } @@ -1062,10 +1073,9 @@ namespace netgen { return (i == 0) ? p1 : p2; } */ - Segment& operator=(const Segment & other); + Segment& operator=(const Segment & other) = default; - int hp_elnr; int GetNP() const { diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index ebcbf578..c79013a3 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -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::register_dtype({ + py::detail::field_descriptor { + "nodes", data_layout["pnum"], + ELEMENT_MAXPOINTS * sizeof(PointIndex), + py::format_descriptor::format(), + py::detail::npy_format_descriptor::dtype() }, + py::detail::field_descriptor { + "index", data_layout["index"], sizeof(int), + py::format_descriptor::format(), + py::detail::npy_format_descriptor::dtype() }, + py::detail::field_descriptor { + "np", data_layout["np"], sizeof(int8_t), + py::format_descriptor::format(), + pybind11::dtype("int8") } + }); + } + py::class_(m, "Element2D") .def(py::init ([](int index, std::vector 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::register_dtype({ + py::detail::field_descriptor { + "nodes", data_layout["pnum"], + ELEMENT2D_MAXPOINTS * sizeof(PointIndex), + py::format_descriptor::format(), + py::detail::npy_format_descriptor::dtype() }, + py::detail::field_descriptor { + "index", data_layout["index"], sizeof(int), + py::format_descriptor::format(), + py::detail::npy_format_descriptor::dtype() }, + py::detail::field_descriptor { + "np", data_layout["np"], sizeof(int8_t), + py::format_descriptor::format(), + pybind11::dtype("int8") } + }); + } + py::class_(m, "Element1D") .def(py::init([](py::list vertices, py::list surfaces, int index, int edgenr, py::list trignums) @@ -557,6 +598,20 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) })) ; + if(ngcore_have_numpy) + { + py::detail::npy_format_descriptor::register_dtype({ + py::detail::field_descriptor { + "nodes", offsetof(Segment, pnums), + 3 * sizeof(PointIndex), + py::format_descriptor::format(), + py::detail::npy_format_descriptor::dtype() }, + py::detail::field_descriptor { + "index", offsetof(Segment, edgenr), sizeof(int), + py::format_descriptor::format(), + py::detail::npy_format_descriptor::dtype() }, + }); + } py::class_(m, "Element0D") .def(py::init([](PointIndex vertex, int index) diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp index b745bf43..7f9ed6ad 100644 --- a/libsrc/meshing/refine.cpp +++ b/libsrc/meshing/refine.cpp @@ -414,7 +414,7 @@ namespace netgen { 2, 4, 9 }, { 3, 4, 10 } }; - int elrev = el.flags.reverse; + int elrev = el.Flags().reverse; for (int j = 1; j <= 4; j++) pnums.Elem(j) = el.PNum(j); @@ -480,10 +480,10 @@ namespace netgen for (int k = 1; k <= 4; k++) nel.PNum(k) = pnums.Get(reftab[j][k-1]); nel.SetIndex(ind); - nel.flags.reverse = reverse[j]; + nel.Flags().reverse = reverse[j]; if (elrev) { - nel.flags.reverse = !nel.flags.reverse; + nel.Flags().reverse = !nel.Flags().reverse; swap (nel.PNum(3), nel.PNum(4)); } @@ -794,10 +794,10 @@ namespace netgen if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0) { wrongels++; - mesh.VolumeElement(i).flags.badel = 1; + mesh.VolumeElement(i).Flags().badel = 1; } else - mesh.VolumeElement(i).flags.badel = 0; + mesh.VolumeElement(i).Flags().badel = 0; if (wrongels) { @@ -900,14 +900,14 @@ namespace netgen if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0) { wrongels++; - mesh.VolumeElement(i).flags.badel = 1; + mesh.VolumeElement(i).Flags().badel = 1; (*testout) << "wrong el: "; for (int j = 1; j <= 4; j++) (*testout) << mesh.VolumeElement(i).PNum(j) << " "; (*testout) << endl; } else - mesh.VolumeElement(i).flags.badel = 0; + mesh.VolumeElement(i).Flags().badel = 0; } cout << "wrongels = " << wrongels << endl; } diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 9e03d5fc..b69b1474 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -473,10 +473,10 @@ namespace netgen if (mesh.VolumeElement(i).CalcJacobianBadness (mesh.Points()) > 1e10) { wrongels++; - mesh.VolumeElement(i).flags.badel = 1; + mesh.VolumeElement(i).Flags().badel = 1; } else - mesh.VolumeElement(i).flags.badel = 0; + mesh.VolumeElement(i).Flags().badel = 0; double facok = 0; double factry; @@ -559,7 +559,7 @@ namespace netgen { wrongels++; Element & el = mesh.VolumeElement(i); - el.flags.badel = 1; + el.Flags().badel = 1; if (lam < 1e-4) @@ -574,7 +574,7 @@ namespace netgen */ } else - mesh.VolumeElement(i).flags.badel = 0; + mesh.VolumeElement(i).Flags().badel = 0; } cout << "wrongels = " << wrongels << endl; } @@ -597,10 +597,10 @@ namespace netgen if (illegalels.Test(i)) { cout << "illegal element: " << i << endl; - mesh.VolumeElement(i).flags.badel = 1; + mesh.VolumeElement(i).Flags().badel = 1; } else - mesh.VolumeElement(i).flags.badel = 0; + mesh.VolumeElement(i).Flags().badel = 0; } /* diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 4fa3ac8a..00108838 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -1320,7 +1320,7 @@ namespace netgen if (mesh->coarsemesh && mesh->hpelements->Size() == mesh->GetNE() ) { 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; } diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index 5b6f5e4f..a4fe826d 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -597,8 +597,8 @@ namespace netgen for (int i = 1; i <= mesh->GetNE(); i++) { - if (mesh->VolumeElement(i).flags.badel || - mesh->VolumeElement(i).flags.illegal || + if (mesh->VolumeElement(i).Flags().badel || + mesh->VolumeElement(i).Flags().illegal || (i == vispar.drawelement)) { // copy to be thread-safe @@ -636,7 +636,7 @@ namespace netgen for (ElementIndex ei : mesh->VolumeElements().Range()) { - if ((*mesh)[ei].flags.badel) + if ((*mesh)[ei].Flags().badel) { // copy to be thread-safe Element el = (*mesh)[ei]; diff --git a/tests/pytest/test_array.py b/tests/pytest/test_array.py index db73452e..fc3efd09 100644 --- a/tests/pytest/test_array.py +++ b/tests/pytest/test_array.py @@ -13,5 +13,21 @@ def test_array_numpy(): a.NumPy().sort() 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__": test_array_numpy() + test_mesh_elements_numpy_array_access()