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

base Index template

Browse Source
This commit is contained in:
Joachim Schoeberl 2024-12-31 12:16:53 +01:00
parent 9efaac072e
commit b7b168e265
15 changed files with 227 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
libsrc/core/array.hpp
View File

@ -216,6 +216,10 @@ namespace ngcore
template <typename T>
constexpr T IndexBASE () { return T(0); }
template <typename T>
constexpr T IndexBASE (T ind) { return IndexBASE<T>(); }
class IndexFromEnd
{

30
libsrc/interface/writegmsh2.cpp
View File

@ -58,7 +58,7 @@ namespace netgen
int np = mesh.GetNP(); /// number of points in mesh
int ne = mesh.GetNE(); /// number of 3D elements in mesh
int nse = mesh.GetNSE(); /// number of surface elements (BC)
int i, j, k, l;
// int i, j, k, l;
/*
@ -66,8 +66,8 @@ namespace netgen
*/
if ((ne > 0)
&& (mesh.VolumeElement(1).GetNP() <= 10)
&& (mesh.SurfaceElement(1).GetNP() <= 6))
&& (mesh.VolumeElements().First().GetNP() <= 10)
&& (mesh.SurfaceElements().First().GetNP() <= 6))
{
cout << "Write GMSH v2.xx Format \n";
cout << "The GMSH v2.xx export is currently available for elements upto 2nd Order\n" << endl;
@ -86,7 +86,7 @@ namespace netgen
outfile << "$Nodes\n";
outfile << np << "\n";
for (i = 1; i <= np; i++)
for (int i = 1; i <= np; i++)
{
const Point3d & p = mesh.Point(i);
outfile << i << " "; /// node number
@ -101,11 +101,11 @@ namespace netgen
outfile << "$Elements\n";
outfile << ne + nse << "\n"; //// number of elements + number of surfaces BC
for (i = 1; i <= nse; i++)
for (auto sei : Range(mesh.SurfaceElements()))
{
int elType = 0;
Element2d el = mesh.SurfaceElement(i);
Element2d el = mesh[sei]; // .SurfaceElement(i);
if(invertsurf) el.Invert();
if(el.GetNP() == 3) elType = GMSH_TRIG; //// GMSH Type for a 3 node triangle
@ -116,7 +116,7 @@ namespace netgen
return;
}
outfile << i;
outfile << sei-IndexBASE(sei)+1;
outfile << " ";
outfile << elType;
outfile << " ";
@ -125,7 +125,7 @@ namespace netgen
outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
/// that means that physical entity = elementary entity (arbitrary approach)
outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
for (j = 1; j <= el.GetNP(); j++)
for (int j = 1; j <= el.GetNP(); j++)
{
outfile << " ";
outfile << el.PNum(triGmsh[j]);
@ -133,12 +133,12 @@ namespace netgen
outfile << "\n";
}
for (i = 1; i <= ne; i++)
for (ElementIndex ei : Range(mesh.VolumeElements()))
{
int i = ei-IndexBASE(ei)+1;
int elType = 0;
Element el = mesh.VolumeElement(i);
Element el = mesh[ei];
if (inverttets) el.Invert();
if(el.GetNP() == 4) elType = GMSH_TET; //// GMSH Element type for 4 node tetrahedron
@ -160,7 +160,7 @@ namespace netgen
outfile << " ";
outfile << 100000 + el.GetIndex(); /// volume number
outfile << " ";
for (j = 1; j <= el.GetNP(); j++)
for (int j = 1; j <= el.GetNP(); j++)
{
outfile << " ";
outfile << el.PNum(tetGmsh[j]);
@ -193,7 +193,7 @@ namespace netgen
outfile << "$Nodes\n";
outfile << np << "\n";
for (i = 1; i <= np; i++)
for (int i = 1; i <= np; i++)
{
const Point3d & p = mesh.Point(i);
outfile << i << " "; /// node number
@ -207,7 +207,7 @@ namespace netgen
outfile << "$Elements\n";
outfile << nse << "\n";
for (k = 1; k <= nse; k++)
for (int k = 1; k <= nse; k++)
{
int elType = 0;
@ -232,7 +232,7 @@ namespace netgen
outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
/// that means that physical entity = elementary entity (arbitrary approach)
outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
for (l = 1; l <= el.GetNP(); l++)
for (int l = 1; l <= el.GetNP(); l++)
{
outfile << " ";
if((elType == GMSH_TRIG) || (elType == GMSH_TRIG6))

2
libsrc/meshing/bisect.cpp
View File

@ -5,10 +5,8 @@
#include "meshing.hpp" // quickfix for parallel
#define noDEBUG
namespace netgen
{
class MarkedTet

8
libsrc/meshing/clusters.cpp
View File

@ -143,19 +143,21 @@ namespace netgen
cluster_reps.Elem(nnums[j]) = nnums[j];
}
*/
ngcore::ParallelForRange
(mesh.SurfaceElements().Range(),
[&] (auto myrange)
{
NgArrayMem<int,9> nnums; // , ednums;
for (int i_ : myrange)
for (SurfaceElementIndex i_ : myrange)
{
int i = i_+1;
const Element2d & el = mesh.SurfaceElement(i);
const Element2d & el = mesh[i_]; // .SurfaceElement(i);
ELEMENT_TYPE typ = el.GetType();
// top.GetSurfaceElementEdges (i, ednums);
auto ednums = top.GetEdges (SurfaceElementIndex(i_));
auto ednums = top.GetEdges (i_);
// cout << "ednums = " << ednums << endl;
int fanum = top.GetSurfaceElementFace (i);

4
libsrc/meshing/improve3.cpp
View File

@ -357,7 +357,7 @@ void MeshOptimize3d :: CombineImprove ()
{
static Timer t("MeshOptimize3d::CombineImprove"); RegionTimer reg(t);
static Timer topt("Optimize");
static Timer tsearch("Search");
static Timer tsearch("Search-combine");
static Timer tbuild_elements_table("Build elements table");
mesh.BuildBoundaryEdges(false);
@ -613,7 +613,7 @@ void MeshOptimize3d :: SplitImprove ()
{
static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);
static Timer topt("Optimize");
static Timer tsearch("Search");
static Timer tsearch("Search-split");
// int np = mesh.GetNP();
int ne = mesh.GetNE();

6
libsrc/meshing/meshclass.cpp
View File

@ -3051,6 +3051,7 @@ namespace netgen
// bool buggy = false;
// ofstream bout("buggy.out");
for (int i = 1; i <= GetNSE(); i++)
{
const Element2d & el = SurfaceElement(i);
@ -3060,10 +3061,10 @@ namespace netgen
{
for (int j = 1; j <= el.GetNP(); j++)
{
INDEX_3 seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex());
PointIndices<3> seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex());
// int data;
if (seg.I1() < PointIndex::BASE || seg.I2() < PointIndex::BASE)
if (!seg.I1().IsValid() || !seg.I2().IsValid())
cerr << "seg = " << seg << endl;
if (faceht.Used(seg))
@ -7299,6 +7300,7 @@ namespace netgen
{
int eli0, eli1;
GetTopology().GetSurface2VolumeElement(sei+1, eli0, eli1);
// auto [ei0,ei1] = GetTopology().GetSurface2VolumeElement(sei); // the way to go
auto & sel = (*this)[sei];
int face = sel.GetIndex();
int domin = VolumeElement(eli0).GetIndex();

2
libsrc/meshing/meshing3.cpp
View File

@ -9,8 +9,6 @@ double minwithoutother;
MeshingStat3d :: MeshingStat3d ()
{
cntsucc = cnttrials = cntelem = qualclass = 0;

133
libsrc/meshing/meshtype.hpp
View File

@ -27,7 +27,6 @@ namespace netgen
*/
enum ELEMENT_TYPE : unsigned char {
SEGMENT = 1, SEGMENT3 = 2,
TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14,
@ -149,6 +148,105 @@ namespace netgen
}
template <typename T, typename TIndex, int Base>
class Index
{
public:
T i;
static constexpr int BASE = Base;
public:
class t_invalid { public: constexpr t_invalid() = default; };
static constexpr t_invalid INVALID{};
constexpr Index () = default;
constexpr Index (const Index&) = default;
constexpr Index (Index &&) = default;
Index & operator= (const Index&) = default;
Index & operator= (Index&&) = default;
// private:
constexpr Index (int ai) : i(ai)
{
#ifdef DEBUG
if (ai < Base)
cout << "illegal Index, use Index::INVALID instead" << endl;
#endif
}
// friend constexpr netgen::TIndex ngcore::IndexBASE<netgen::TIndex> ();
// template <int N> friend class PointIndices;
/*
friend auto operator+ (Index, int) -> TIndex;
friend TIndex operator+ (Index, size_t);
friend TIndex operator+ (int, Index);
friend TIndex operator+ (size_t, Index);
friend constexpr TIndex operator- (Index, int);
friend int operator- (Index, Index);
friend bool operator< (Index a, Index b);
friend bool operator> (Index a, Index b);
friend bool operator>= (Index a, Index b);
friend bool operator<= (Index a, Index b);
friend bool operator== (Index a, Index b);
friend bool operator!= (Index a, Index b);
*/
public:
constexpr Index (t_invalid inv) : i(long(BASE)-1) { ; }
// TIndex & operator= (const TIndex &ai) { i = ai.i; return *this; }
// private:
constexpr operator const int& () const { return i; }
explicit constexpr operator int& () { return i; }
public:
constexpr operator TIndex() const { return TIndex(i); }
constexpr operator TIndex&() { return static_cast<TIndex&>(*this); }
TIndex operator++ (int) { TIndex hi(*this); i++; return hi; }
TIndex operator-- (int) { TIndex hi(*this); i--; return hi; }
TIndex & operator++ () { i++; return *this; }
TIndex operator-- () { i--; return *this; }
TIndex operator+= (int add) { i += add; return *this; }
void Invalidate() { i = long(TIndex::BASE)-1; }
bool IsValid() const { return i+1 != TIndex::BASE; }
// operator bool() const { return IsValid(); }
void DoArchive (Archive & ar) { ar & i; }
};
/*
template <typename T, typename TIndex, int Base>
auto operator+ (Index<T,TIndex,Base> pi, int i) -> TIndex { return TIndex(pi.i+i); }
template <typename T, typename TIndex, int Base>
inline TIndex operator+ (Index<T,TIndex,Base> pi, size_t i) { return TIndex(pi.i+i); }
template <typename T, typename TIndex, int Base>
inline TIndex operator+ (int i, Index<T,TIndex,Base> pi) { return TIndex(pi.i+i); }
template <typename T, typename TIndex, int Base>
inline TIndex operator+ (size_t i, Index<T,TIndex,Base> pi) { return TIndex(pi.i+i); }
template <typename T, typename TIndex, int Base>
constexpr inline auto operator- (Index<T,TIndex,Base> pi, int i) -> TIndex { return TIndex(pi.i-i); }
template <typename T, typename TIndex, int Base>
inline int operator- (Index<T,TIndex,Base> pa, Index<T,TIndex,Base> pb) { return pa.i-pb.i; }
template <typename T, typename TIndex, int Base>
inline bool operator< (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i < 0; }
template <typename T, typename TIndex, int Base>
inline bool operator> (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i > 0; }
template <typename T, typename TIndex, int Base>
inline bool operator>= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i >= 0; }
template <typename T, typename TIndex, int Base>
inline bool operator<= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i-b.i <= 0; }
template <typename T, typename TIndex, int Base>
inline bool operator== (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i == b.i; }
template <typename T, typename TIndex, int Base>
inline bool operator!= (Index<T,TIndex,Base> a, Index<T,TIndex,Base> b) { return a.i != b.i; }
*/
/*
@ -197,6 +295,8 @@ namespace netgen
// #define BASE0
/*
class PointIndex
{
int i;
@ -260,6 +360,14 @@ namespace netgen
void DoArchive (Archive & ar) { ar & i; }
};
*/
class PointIndex : public Index<int,PointIndex,1>
{
public:
using Index::Index;
};
constexpr inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); }
constexpr inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); }
constexpr inline PointIndex operator+ (int i, PointIndex pi) { return PointIndex(pi.i+i); }
@ -501,8 +609,8 @@ namespace netgen
inline bool operator> (ElementIndex ei1, ElementIndex ei2) { return int(ei1) > int(ei2); };
inline bool operator>= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) >= int(ei2); };
inline bool operator<= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) <= int(ei2); };
// these should not be needed:
// these should not be needed:
inline bool operator== (ElementIndex ei1, int ei2) { return int(ei1) == int(ei2); };
inline bool operator< (size_t s, ElementIndex ei2) { return int(s) < int(ei2); };
inline bool operator< (ElementIndex ei1, size_t s) { return int(ei1) < int(s); }; // should not need
@ -531,6 +639,27 @@ namespace netgen
void DoArchive (Archive & ar) { ar & i; }
};
inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, int i) { return SurfaceElementIndex { int(ei) + i }; }
inline SurfaceElementIndex operator+ (size_t s, SurfaceElementIndex ei) { return SurfaceElementIndex(int(ei) + s); }
inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, size_t s) { return SurfaceElementIndex(int(ei) + s); }
inline bool operator== (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) == int(ei2); };
inline bool operator!= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) != int(ei2); };
inline bool operator< (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) < int(ei2); };
inline bool operator> (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) > int(ei2); };
inline bool operator>= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) >= int(ei2); };
inline bool operator<= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) <= int(ei2); };
// these should not be needed:
inline bool operator== (SurfaceElementIndex ei1, int ei2) { return int(ei1) == int(ei2); };
inline bool operator== (int ei2, SurfaceElementIndex ei1) { return int(ei1) == int(ei2); };
inline bool operator!= (SurfaceElementIndex ei1, int ei2) { return int(ei1) != int(ei2); };
inline bool operator< (size_t s, SurfaceElementIndex ei2) { return int(s) < int(ei2); };
inline bool operator< (SurfaceElementIndex ei1, size_t s) { return int(ei1) < int(s); }; // should not need
inline bool operator< (SurfaceElementIndex ei1, int s) { return int(ei1) < int(s); }; // should not need
inline bool operator>= (size_t s, SurfaceElementIndex ei2) { return int(s) >= int(ei2); };
inline bool operator>= (SurfaceElementIndex ei1, int s) { return int(ei1) >= int(s); };
inline void SetInvalid (SurfaceElementIndex & id) { id = -1; }
inline bool IsInvalid (SurfaceElementIndex & id) { return id == -1; }

11
libsrc/meshing/parallelmesh.cpp
View File

@ -593,17 +593,16 @@ namespace netgen
DynamicTable<int> elementarrays(elarraysize);
for (int ei = 1; ei <= GetNE(); ei++)
for (ElementIndex ei : VolumeElements().Range())
{
const Element & el = VolumeElement (ei);
// int dest = el.GetPartition();
int dest = vol_partition[ei-1];
int dest = vol_partition[ei];
elementarrays.Add (dest, ei);
elementarrays.Add (dest, int(ei+1));
elementarrays.Add (dest, el.GetIndex());
elementarrays.Add (dest, el.GetNP());
for (int i = 0; i < el.GetNP(); i++)
elementarrays.Add (dest, el[i]);
for (PointIndex pi : el.PNums())
elementarrays.Add (dest, pi);
}
tbuildelementtable.Stop();

15
libsrc/meshing/refine.cpp
View File

@ -865,9 +865,10 @@ namespace netgen
for (int k = 1; k <= 3; k++)
{
fhelp.Clear();
for (int i = 1; i <= mesh.GetNE(); i++)
// for (int i = 1; i <= mesh.GetNE(); i++)
for (const Element & el : mesh.VolumeElements())
{
const Element & el = mesh.VolumeElement(i);
// const Element & el = mesh.VolumeElement(i);
int freeel = 0;
for (int j = 1; j <= el.GetNP(); j++)
if (free.Test(el.PNum(j)))
@ -897,19 +898,19 @@ namespace netgen
wrongels = 0;
for (int i = 1; i <= mesh.GetNE(); i++)
for (ElementIndex ei : mesh.VolumeElements().Range())
{
if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
if (mesh.VolumeElement(ei).Volume(mesh.Points()) < 0)
{
wrongels++;
mesh.VolumeElement(i).Flags().badel = 1;
mesh.VolumeElement(ei).Flags().badel = 1;
(*testout) << "wrong el: ";
for (int j = 1; j <= 4; j++)
(*testout) << mesh.VolumeElement(i).PNum(j) << " ";
(*testout) << mesh.VolumeElement(ei).PNum(j) << " ";
(*testout) << endl;
}
else
mesh.VolumeElement(i).Flags().badel = 0;
mesh.VolumeElement(ei).Flags().badel = 0;
}
cout << "wrongels = " << wrongels << endl;
}

7
libsrc/meshing/topology.hpp
View File

@ -185,6 +185,13 @@ public:
elnr2 = surf2volelement.Get(selnr)[1];
}
std::array<ElementIndex,2> GetSurface2VolumeElement (SurfaceElementIndex sei)
{
return { ElementIndex( surf2volelement.Get(sei+1)[0] - 1),
ElementIndex( surf2volelement.Get(sei+1)[1] - 1) };
}
int GetFace2SurfaceElement (int fnr) const { return face2surfel[fnr-1]; }
SegmentIndex GetSegmentOfEdge(int edgenr) const { return edge2segment[edgenr-1]; }

17
libsrc/occ/occmeshsurf.cpp
View File

@ -366,8 +366,11 @@ namespace netgen
void OCCSurface :: Project (Point<3> & ap, PointGeomInfo & gi)
{
// static Timer t("OccSurface::Project"); RegionTimer reg(t);
// static Timer t2("OccSurface::Project actual");
static Timer t("OccSurface::Project"); RegionTimer reg(t);
static Timer tanal("OccSurface::Project analysis");
static Timer ttol("OccSurface::Project approximation");
static Timer t2("OccSurface::Project actual");
// try Newton's method ...
@ -473,13 +476,17 @@ namespace netgen
// double u,v;
// JS : shouldn't we move these 2 lines to the constructor ?
// tanal.Start();
Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
// ShapeAnalysis_Surface su( occface );
// tanal.Stop();
ttol.Start();
auto toltool = BRep_Tool::Tolerance( topods_face );
ttol.Stop();
// gp_Pnt2d suval = su->ValueOfUV ( pnt, toltool);
// t2.Start();
t2.Start();
gp_Pnt2d suval = su->NextValueOfUV (gp_Pnt2d(u,v), pnt, toltool);
// t2.Stop();
t2.Stop();
suval.Coord( u, v);
pnt = occface->Value( u, v );

17
libsrc/occ/occmeshsurf.hpp
View File

@ -5,10 +5,13 @@
#include "occgeom.hpp"
#include "mydefs.hpp"
#include <BRep_Tool.hxx>
#include <TopoDS_Face.hxx>
#include <Geom_Surface.hxx>
#include <ShapeAnalysis.hxx>
#include <ShapeAnalysis_Surface.hxx>
#include <GeomLProp_SLProps.hxx>
#define PARAMETERSPACE -1
#define PLANESPACE 1
@ -30,7 +33,8 @@ public:
Handle(Geom_Surface) occface;
TopAbs_Orientation orient;
int projecttype;
ShapeAnalysis_Surface su;
Standard_Real toltool;
protected:
Point<3> p1;
Point<3> p2;
@ -60,10 +64,15 @@ protected:
public:
OCCSurface (const TopoDS_Face & aface, int aprojecttype)
: topods_face(aface),
occface(BRep_Tool::Surface(topods_face)),
su( occface ),
toltool(BRep_Tool::Tolerance(topods_face))
{
static Timer t("occurface ctor"); RegionTimer r(t);
topods_face = aface;
occface = BRep_Tool::Surface(topods_face);
// occface = BRep_Tool::Surface(topods_face);
orient = topods_face.Orientation();
projecttype = aprojecttype;
ShapeAnalysis::GetFaceUVBounds (topods_face, umin, umax, vmin, vmax);
@ -72,6 +81,8 @@ public:
umax += fabs(umax-umin)/100.0;
vmax += fabs(vmax-vmin)/100.0;
// projecttype = PLANESPACE;
// su = ShapeAnalysis_Surface( occface );
/*
TopExp_Explorer exp1;
exp1.Init (topods_face, TopAbs_WIRE);

8
libsrc/occ/occpkg.cpp
View File

@ -817,20 +817,28 @@ namespace netgen
if(strcmp(argv[1], "showall") == 0)
{
/*
for(int i = 1; i <= mesh->GetNSE(); i++)
{
mesh->SurfaceElement(i).Visible(1);
}
*/
for (auto & el : mesh->SurfaceElements())
el.Visible(1);
mesh->SetNextTimeStamp();
}
if(strcmp(argv[1], "hideall") == 0)
{
/*
for(int i = 1; i <= mesh->GetNSE(); i++)
{
mesh->SurfaceElement(i).Visible(0);
}
*/
for (auto & el : mesh->SurfaceElements())
el.Visible(0);
mesh->SetNextTimeStamp();
}

6
nglib/nglib.cpp
View File

@ -270,7 +270,7 @@ namespace nglib
NGLIB_API Ng_Surface_Element_Type
Ng_GetSurfaceElement (Ng_Mesh * mesh, int num, int * pi)
{
const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1));
for (int i = 1; i <= el.GetNP(); i++)
pi[i-1] = el.PNum(i);
Ng_Surface_Element_Type et;
@ -301,7 +301,7 @@ namespace nglib
NGLIB_API Ng_Volume_Element_Type
Ng_GetVolumeElement (Ng_Mesh * mesh, int num, int * pi)
{
const Element & el = ((Mesh*)mesh)->VolumeElement(num);
const Element & el = ((Mesh*)mesh)->VolumeElement(ElementIndex(num-1));
for (int i = 1; i <= el.GetNP(); i++)
pi[i-1] = el.PNum(i);
Ng_Volume_Element_Type et;
@ -439,7 +439,7 @@ namespace nglib
NGLIB_API Ng_Surface_Element_Type
Ng_GetElement_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum)
{
const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1));
for (int i = 1; i <= el.GetNP(); i++)
pi[i-1] = el.PNum(i);