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

use ngscore::Array for surface elements

Browse Source
This commit is contained in:
Joachim Schöberl 2019-08-08 08:44:59 +02:00
parent a95185a714
commit e8960ebae1
18 changed files with 107 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
libsrc/interface/writeuser.cpp
View File

@ -427,7 +427,7 @@ void WriteSTLExtFormat (const Mesh & mesh,
for(int faceNr = 1;faceNr <= faceBCMapping.EntrySize(bcInd); faceNr++)
{
NgArray<SurfaceElementIndex> faceSei;
Array<SurfaceElementIndex> faceSei;
mesh.GetSurfaceElementsOfFace(faceBCMapping.Get(bcInd,faceNr),faceSei);
for (int i = 0; i < faceSei.Size(); i++)

6
libsrc/meshing/bcfunctions.cpp
View File

@ -341,9 +341,9 @@ namespace netgen
// Extract the number of surface elements having a given colour
// And save this number into an array for later sorting
for(int face_index = 1; face_index <= nfd; face_index++)
{
NgArray<SurfaceElementIndex> se_face;
{
Array<SurfaceElementIndex> se_face;
mesh.GetSurfaceElementsOfFace(face_index, se_face);
Vec3d face_colour;

57
libsrc/meshing/improve2.cpp
View File

@ -74,16 +74,24 @@ namespace netgen
NgProfiler::StartTimer (timerstart);
NgArray<SurfaceElementIndex> seia;
Array<SurfaceElementIndex> seia;
mesh.GetSurfaceElementsOfFace (faceindex, seia);
/*
for (int i = 0; i < seia.Size(); i++)
if (mesh[seia[i]].GetNP() != 3)
{
GenericImprove (mesh);
return;
}
*/
for (SurfaceElementIndex sei : seia)
if (mesh[sei].GetNP() != 3)
{
GenericImprove (mesh);
return;
}
int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
NgArray<Neighbour> neighbors(mesh.GetNSE());
@ -356,8 +364,8 @@ namespace netgen
Element2d sw2 (pi3, pi4, pi2);
int legal1 =
mesh.LegalTrig (mesh.SurfaceElement (t1)) +
mesh.LegalTrig (mesh.SurfaceElement (t2));
mesh.LegalTrig (mesh[t1]) +
mesh.LegalTrig (mesh[t2]);
int legal2 =
mesh.LegalTrig (sw1) + mesh.LegalTrig (sw2);
@ -438,35 +446,24 @@ namespace netgen
static int timerstart1 = NgProfiler::CreateTimer ("Combineimprove 2D start1");
NgProfiler::StartTimer (timerstart1);
// int i, j, k, l;
// PointIndex pi;
// SurfaceElementIndex sei;
NgArray<SurfaceElementIndex> seia;
Array<SurfaceElementIndex> seia;
mesh.GetSurfaceElementsOfFace (faceindex, seia);
for (int i = 0; i < seia.Size(); i++)
if (mesh[seia[i]].GetNP() != 3)
for (SurfaceElementIndex sei : seia)
if (mesh[sei].GetNP() != 3)
return;
int surfnr = 0;
if (faceindex)
surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
// PointIndex pi1, pi2;
// MeshPoint p1, p2, pnew;
double bad1, bad2;
Vec<3> nv;
int np = mesh.GetNP();
//int nse = mesh.GetNSE();
TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonnode(np);
NgArray<SurfaceElementIndex> hasonepi, hasbothpi;
@ -639,8 +636,9 @@ namespace netgen
hasonepi.Append (elementsonnode[pi2][k]);
}
bad1 = 0;
double bad1 = 0;
int illegal1 = 0, illegal2 = 0;
for (int k = 0; k < hasonepi.Size(); k++)
{
const Element2d & el = mesh[hasonepi[k]];
@ -665,7 +663,7 @@ namespace netgen
mesh[pi1] = pnew;
mesh[pi2] = pnew;
bad2 = 0;
double bad2 = 0;
for (int k = 0; k < hasonepi.Size(); k++)
{
Element2d & el = mesh[hasonepi[k]];
@ -753,24 +751,17 @@ namespace netgen
for (SurfaceElementIndex sei2 : elementsonnode[pi2])
{
Element2d & el = mesh[sei2];
if(el.IsDeleted()) continue;
elementsonnode.Add (pi1, sei2);
/*
bool haspi1 = 0;
for (int l = 0; l < el.GetNP(); l++)
if (el[l] == pi1)
haspi1 = true;
if (haspi1) continue;
*/
if (el.IsDeleted()) continue;
if (el.PNums().Contains(pi1)) continue;
elementsonnode.Add (pi1, sei2);
for (int l = 0; l < el.GetNP(); l++)
for (auto l : Range(el.GetNP()))
{
if (el[l] == pi2)
{
el[l] = pi1;
el.GeomInfoPi (l+1) = gi;
el.GeomInfo()[l] = gi;
}
fixed[el[l]] = true;

2
libsrc/meshing/improve2gen.cpp
View File

@ -430,7 +430,7 @@ namespace netgen
}
for (int j = 0; j < rule.oldels.Size(); j++)
mesh.DeleteSurfaceElement ( elmap[j] );
mesh.Delete (elmap[j]);
for (int j = 1; j <= pmap.Size(); j++)
nelonnode[pmap.Get(j)] += rule.incelsonnode.Get(j);

8
libsrc/meshing/improve3.hpp
View File

@ -3,7 +3,7 @@
extern double CalcTotalBad (const Mesh::T_POINTS & points,
const Mesh::T_VOLELEMENTS & elements,
const NgArray<Element> & elements,
const MeshingParameters & mp);
@ -33,7 +33,7 @@ public:
double CalcTotalBad (const Mesh::T_POINTS & points,
const Mesh::T_VOLELEMENTS & elements)
const NgArray<Element> & elements)
{
return netgen::CalcTotalBad (points, elements, mp);
}
@ -100,7 +100,7 @@ class JacobianPointFunction : public MinFunction
{
public:
Mesh::T_POINTS & points;
const Mesh::T_VOLELEMENTS & elements;
const NgArray<Element> & elements;
TABLE<INDEX> elementsonpoint;
PointIndex actpind;
@ -109,7 +109,7 @@ public:
public:
JacobianPointFunction (Mesh::T_POINTS & apoints,
const Mesh::T_VOLELEMENTS & aelements);
const NgArray<Element> & aelements);
virtual ~JacobianPointFunction () { ; }
virtual void SetPointIndex (PointIndex aactpind);
virtual double Func (const Vector & x) const;

16
libsrc/meshing/meshclass.cpp
View File

@ -2467,7 +2467,7 @@ namespace netgen
for (int i = 1; i <= GetNSE(); i++)
{
Element2d & sel = surfelements.Elem(i);
Element2d & sel = surfelements[i-1];
bool remove = false;
for (int j = 1; j <= sel.GetNP(); j++)
if (frontpoints.Test(sel.PNum(j)))
@ -2478,9 +2478,9 @@ namespace netgen
for (int i = surfelements.Size(); i >= 1; i--)
{
if (!surfelements.Elem(i).PNum(1).IsValid())
if (!surfelements[i-1].PNum(1).IsValid())
{
surfelements.Elem(i) = surfelements.Last();
surfelements[i-1] = surfelements.Last();
surfelements.DeleteLast();
}
}
@ -3270,8 +3270,12 @@ namespace netgen
for (int i = 0; i < surfelements.Size(); i++)
if (surfelements[i].IsDeleted())
{
surfelements[i] = surfelements.Last();
surfelements.DeleteLast();
/*
surfelements.Delete(i);
i--;
*/
}
for (int i = 0; i < segments.Size(); i++)
@ -3907,7 +3911,7 @@ namespace netgen
for (i = 1; i <= nse; i++)
if (!used.Test(i))
{
Element2d & el = surfelements.Elem(i);
Element2d & el = surfelements[i-1];
int found = 0, foundrev = 0;
for (j = 1; j <= 3; j++)
{
@ -5372,7 +5376,7 @@ namespace netgen
int np = GetNP();
BitArray usedp(np);
NgArray<SurfaceElementIndex> els_of_face;
Array<SurfaceElementIndex> els_of_face;
fdi = 1;
while (fdi <= GetNFD())
@ -5546,7 +5550,7 @@ namespace netgen
}
}
void Mesh :: GetSurfaceElementsOfFace (int facenr, NgArray<SurfaceElementIndex> & sei) const
void Mesh :: GetSurfaceElementsOfFace (int facenr, Array<SurfaceElementIndex> & sei) const
{
static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace");
NgProfiler::RegionTimer reg (timer);

52
libsrc/meshing/meshclass.hpp
View File

@ -24,9 +24,9 @@ namespace netgen
{
public:
typedef ::netgen::T_POINTS T_POINTS;
typedef NgArray<Element, 0, size_t> T_VOLELEMENTS;
// typedef NgArray<Element, 0, size_t> T_VOLELEMENTS;
// typedef NgArray<Element2d, 0, SurfaceElementIndex> T_SURFELEMENTS;
typedef NgArray<Element2d, 0, size_t> T_SURFELEMENTS;
// typedef NgArray<Element2d, 0, size_t> T_SURFELEMENTS;
private:
/// point coordinates
@ -38,9 +38,9 @@ namespace netgen
/// line-segments at edges
NgArray<Segment, 0, size_t> segments;
/// surface elements, 2d-inner elements
T_SURFELEMENTS surfelements;
Array<Element2d> surfelements;
/// volume elements
T_VOLELEMENTS volelements;
NgArray<Element> volelements;
/// points will be fixed forever
NgArray<PointIndex> lockedpoints;
@ -263,10 +263,6 @@ namespace netgen
const Segment & operator[] (SegmentIndex si) const { return segments[si]; }
Segment & operator[] (SegmentIndex si) { return segments[si]; }
/*
const NgArray<Segment> & LineSegments() const { return segments; }
NgArray<Segment> & LineSegments() { return segments; }
*/
const auto & LineSegments() const { return segments; }
auto & LineSegments() { return segments; }
@ -276,30 +272,46 @@ namespace netgen
// write to pre-allocated container, thread-safe
DLL_HEADER void SetSurfaceElement (SurfaceElementIndex sei, const Element2d & el);
// [[deprecated("Use DeleteSurfaceElement(SurfaceElementIndex) instead of int !")]]
[[deprecated("Use Delete(SurfaceElementIndex) instead of int !")]]
void DeleteSurfaceElement (int eli)
{
{
/*
surfelements.Elem(eli).Delete();
surfelements.Elem(eli).PNum(1).Invalidate();
surfelements.Elem(eli).PNum(2).Invalidate();
surfelements.Elem(eli).PNum(3).Invalidate();
*/
surfelements[eli-1].Delete();
surfelements[eli-1].PNum(1).Invalidate();
surfelements[eli-1].PNum(2).Invalidate();
surfelements[eli-1].PNum(3).Invalidate();
timestamp = NextTimeStamp();
}
[[deprecated("Use Delete(SurfaceElementIndex) instead !")]]
void DeleteSurfaceElement (SurfaceElementIndex eli)
{
for (auto & p : surfelements[eli].PNums()) p.Invalidate();
surfelements[eli].Delete();
timestamp = NextTimeStamp();
}
void Delete (SurfaceElementIndex eli)
{
for (auto & p : surfelements[eli].PNums()) p.Invalidate();
surfelements[eli].Delete();
timestamp = NextTimeStamp();
}
auto GetNSE () const { return surfelements.Size(); }
// [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]
Element2d & SurfaceElement(int i) { return surfelements.Elem(i); }
Element2d & SurfaceElement(int i) { return surfelements[i-1]; }
// [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]
const Element2d & SurfaceElement(int i) const { return surfelements.Get(i); }
const Element2d & SurfaceElement(int i) const { return surfelements[i-1]; }
[[deprecated("Use mesh[](SurfaceElementIndex) instead !")]]
Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; }
[[deprecated("Use mesh[](SurfaceElementIndex) instead !")]]
const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; }
const Element2d & operator[] (SurfaceElementIndex ei) const
@ -307,12 +319,12 @@ namespace netgen
Element2d & operator[] (SurfaceElementIndex ei)
{ return surfelements[ei]; }
const T_SURFELEMENTS & SurfaceElements() const { return surfelements; }
T_SURFELEMENTS & SurfaceElements() { return surfelements; }
const auto & SurfaceElements() const { return surfelements; }
auto & SurfaceElements() { return surfelements; }
DLL_HEADER void RebuildSurfaceElementLists ();
DLL_HEADER void GetSurfaceElementsOfFace (int facenr, NgArray<SurfaceElementIndex> & sei) const;
DLL_HEADER void GetSurfaceElementsOfFace (int facenr, Array<SurfaceElementIndex> & sei) const;
DLL_HEADER ElementIndex AddVolumeElement (const Element & el);
// write to pre-allocated container, thread-safe
@ -324,15 +336,13 @@ namespace netgen
Element & VolumeElement(int i) { return volelements.Elem(i); }
// [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]]
const Element & VolumeElement(int i) const { return volelements.Get(i); }
[[deprecated("Use mesh[](VolumeElementIndex) instead !")]]
Element & VolumeElement(ElementIndex i) { return volelements[i]; }
[[deprecated("Use mesh[](VolumeElementIndex) instead !")]]
const Element & VolumeElement(ElementIndex i) const { return volelements[i]; }
const Element & operator[] (ElementIndex ei) const
{ return volelements[ei]; }
Element & operator[] (ElementIndex ei)
{ return volelements[ei]; }
const Element & operator[] (ElementIndex ei) const { return volelements[ei]; }
Element & operator[] (ElementIndex ei) { return volelements[ei]; }
ELEMENTTYPE ElementType (ElementIndex i) const
{ return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; }

2
libsrc/meshing/meshing2.cpp
View File

@ -305,7 +305,7 @@ namespace netgen
}
}
*/
NgArray<SurfaceElementIndex> seia;
Array<SurfaceElementIndex> seia;
mesh.GetSurfaceElementsOfFace (facenr, seia);
for (int i = 0; i < seia.Size(); i++)
{

6
libsrc/meshing/parallelmesh.cpp
View File

@ -749,9 +749,9 @@ namespace netgen
auto & self = const_cast<Mesh&>(*this);
self.points = T_POINTS(0);
self.surfelements = T_SURFELEMENTS(0);
self.volelements = T_VOLELEMENTS(0);
self.segments = NgArray<Segment, 0, size_t>(0);
self.surfelements = Array<Element2d>(0);
self.volelements = NgArray<Element>(0);
self.segments = NgArray<Segment>(0);
self.lockedpoints = NgArray<PointIndex>(0);
auto cleanup_ptr = [](auto & ptr) {
if (ptr != nullptr) {

6
libsrc/meshing/python_mesh.cpp
View File

@ -740,11 +740,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
.def("Elements3D",
static_cast<NgArray<Element,0,size_t>&(Mesh::*)()> (&Mesh::VolumeElements),
static_cast<NgArray<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
py::return_value_policy::reference)
.def("Elements2D",
static_cast<NgArray<Element2d,0,size_t>&(Mesh::*)()> (&Mesh::SurfaceElements),
static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements),
py::return_value_policy::reference)
.def("Elements1D",
@ -813,7 +813,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def ("DeleteSurfaceElement",
FunctionPointer ([](Mesh & self, SurfaceElementIndex i)
{
return self.DeleteSurfaceElement (i);
return self.Delete(i);
}))
.def ("Compress", FunctionPointer ([](Mesh & self)

2
libsrc/meshing/smoothing2.5.cpp
View File

@ -51,7 +51,7 @@ namespace netgen
int i, j, k;
SurfaceElementIndex sei;
NgArray<SurfaceElementIndex> seia;
Array<SurfaceElementIndex> seia;
mesh.GetSurfaceElementsOfFace (faceindex, seia);
/*

8
libsrc/meshing/smoothing2.cpp
View File

@ -741,13 +741,13 @@ namespace netgen
Opti2dLocalData ld;
NgArray<SurfaceElementIndex> seia;
Array<SurfaceElementIndex> seia;
mesh.GetSurfaceElementsOfFace (faceindex, seia);
bool mixed = 0;
for (int i = 0; i < seia.Size(); i++)
if (mesh[seia[i]].GetNP() != 3)
for (auto sei : seia)
if (mesh[sei].GetNP() != 3)
{
mixed = 1;
mixed = true;
break;
}

18
libsrc/meshing/smoothing3.cpp
View File

@ -300,7 +300,7 @@ namespace netgen
{
public:
Mesh::T_POINTS & points;
const Mesh::T_VOLELEMENTS & elements;
const NgArray<Element> & elements;
TABLE<int,PointIndex::BASE> elementsonpoint;
const MeshingParameters & mp;
PointIndex actpind;
@ -308,7 +308,7 @@ namespace netgen
public:
PointFunction (Mesh::T_POINTS & apoints,
const Mesh::T_VOLELEMENTS & aelements,
const NgArray<Element> & aelements,
const MeshingParameters & amp);
virtual ~PointFunction () { ; }
virtual void SetPointIndex (PointIndex aactpind);
@ -323,7 +323,7 @@ namespace netgen
PointFunction :: PointFunction (Mesh::T_POINTS & apoints,
const Mesh::T_VOLELEMENTS & aelements,
const NgArray<Element> & aelements,
const MeshingParameters & amp)
: points(apoints), elements(aelements), elementsonpoint(apoints.Size()), mp(amp)
{
@ -477,7 +477,7 @@ namespace netgen
DenseMatrix m;
public:
CheapPointFunction (Mesh::T_POINTS & apoints,
const Mesh::T_VOLELEMENTS & aelements,
const NgArray<Element> & aelements,
const MeshingParameters & amp);
virtual void SetPointIndex (PointIndex aactpind);
virtual double PointFunctionValue (const Point<3> & pp) const;
@ -486,7 +486,7 @@ namespace netgen
CheapPointFunction :: CheapPointFunction (Mesh::T_POINTS & apoints,
const Mesh::T_VOLELEMENTS & aelements,
const NgArray<Element> & aelements,
const MeshingParameters & amp)
: PointFunction (apoints, aelements, amp)
{
@ -920,7 +920,7 @@ double Opti3EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
double CalcTotalBad (const Mesh::T_POINTS & points,
const Mesh::T_VOLELEMENTS & elements,
const NgArray<Element> & elements,
const MeshingParameters & mp)
{
static Timer t("CalcTotalBad"); RegionTimer reg(t);
@ -985,13 +985,13 @@ int WrongOrientation (const Mesh::T_POINTS & points, const Element & el)
// {
// public:
// Mesh::T_POINTS & points;
// const Mesh::T_VOLELEMENTS & elements;
// const NgArray<Element> & elements;
// TABLE<INDEX> elementsonpoint;
// PointIndex actpind;
// public:
// JacobianPointFunction (Mesh::T_POINTS & apoints,
// const Mesh::T_VOLELEMENTS & aelements);
// const NgArray<Element> & aelements);
// virtual void SetPointIndex (PointIndex aactpind);
// virtual double Func (const Vector & x) const;
@ -1002,7 +1002,7 @@ int WrongOrientation (const Mesh::T_POINTS & points, const Element & el)
JacobianPointFunction ::
JacobianPointFunction (Mesh::T_POINTS & apoints,
const Mesh::T_VOLELEMENTS & aelements)
const NgArray<Element> & aelements)
: points(apoints), elements(aelements), elementsonpoint(apoints.Size())
{
INDEX i;

2
libsrc/meshing/topology.cpp
View File

@ -219,7 +219,7 @@ namespace netgen
for (SurfaceElementIndex elnr : top.GetVertexSurfaceElements(v))
{
const Element2d & el = mesh.SurfaceElement (elnr);
const Element2d & el = mesh[elnr];
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (el.GetType());

4
libsrc/occ/occgenmesh.cpp
View File

@ -831,8 +831,8 @@ namespace netgen
{
if (notrys == 1)
{
for (int i = noldsurfel+1; i <= mesh.GetNSE(); i++)
mesh.DeleteSurfaceElement (i);
for (SurfaceElementIndex sei = noldsurfel; sei < mesh.GetNSE(); sei++)
mesh.Delete(sei);
mesh.Compress();

10
libsrc/occ/occpkg.cpp
View File

@ -711,7 +711,7 @@ namespace netgen
for(int i = 1; i <= mesh->GetNFD(); i++)
{
NgArray<SurfaceElementIndex> surfElems;
Array<SurfaceElementIndex> surfElems;
mesh->GetSurfaceElementsOfFace(i,surfElems);
if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))
@ -735,7 +735,7 @@ namespace netgen
for(int i = 1; i <= mesh->GetNFD(); i++)
{
NgArray<SurfaceElementIndex> surfElems;
Array<SurfaceElementIndex> surfElems;
mesh->GetSurfaceElementsOfFace(i,surfElems);
if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))
@ -759,7 +759,7 @@ namespace netgen
for(int i = 1; i <= mesh->GetNFD(); i++)
{
NgArray<SurfaceElementIndex> surfElems;
Array<SurfaceElementIndex> surfElems;
mesh->GetSurfaceElementsOfFace(i,surfElems);
if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))
@ -783,14 +783,14 @@ namespace netgen
if(strcmp(argv[1], "hideonly") == 0)
{
NgArray<Vec3d> face_colours;
NgArray<Vec3d> face_colours;
GetFaceColours(*mesh,face_colours);
int colourind = atoi (argv[2]);
for(int i = 1; i <= mesh->GetNFD(); i++)
{
NgArray<SurfaceElementIndex> surfElems;
Array<SurfaceElementIndex> surfElems;
mesh->GetSurfaceElementsOfFace(i,surfElems);
if(ColourMatch(face_colours[colourind],mesh->GetFaceDescriptor(i).SurfColour()))

6
libsrc/visualization/vsfieldlines.cpp
View File

@ -576,7 +576,7 @@ namespace netgen
shared_ptr<Mesh> mesh = GetMesh();
if (!mesh) return;
NgArray<SurfaceElementIndex> elements_2d;
Array<SurfaceElementIndex> elements_2d;
//cout << "fieldlines_startface " << fieldlines_startface << endl;
mesh->GetSurfaceElementsOfFace(fieldlines_startface,elements_2d);
@ -592,7 +592,7 @@ namespace netgen
int i;
for(i=0; i<elements_2d.Size(); i++)
{
const Element2d & elem = mesh->SurfaceElement(elements_2d[i]);
const Element2d & elem = (*mesh)[elements_2d[i]];
v1 = mesh->Point(elem[1]) - mesh->Point(elem[0]);
v2 = mesh->Point(elem[2]) - mesh->Point(elem[0]);
@ -613,7 +613,7 @@ namespace netgen
while(startpointsp < startpoints.Size())
{
const Element2d & elem = mesh->SurfaceElement(elements_2d[i]);
const Element2d & elem = (*mesh)[elements_2d[i]];
int numtri = (elem.GetNV() == 3) ? 1 : 2;

2
libsrc/visualization/vsmesh.cpp
View File

@ -1005,7 +1005,7 @@ namespace netgen
int hoplotn = 1 << vispar.subdivisions;
NgArray<SurfaceElementIndex> seia;
Array<SurfaceElementIndex> seia;
for (int faceindex = 1; faceindex <= mesh->GetNFD(); faceindex++)