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

Merge remote-tracking branch 'gitlab/master' into external_project_metis

Browse Source
This commit is contained in:
Matthias Hochsteger 2017-03-10 09:46:30 +01:00
commit 2ee076da4d
27 changed files with 424 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
.gitlab-ci.yml
View File

@ -138,7 +138,12 @@ test_win64:
.template_test_linux: &test_linux .template_test_linux: &test_linux
stage: test stage: test
script: script:
- docker run netgen_${CI_BUILD_REF_NAME}_installed:${UBUNTU_VERSION} bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V"' - >-
docker run
-e NETGENDIR=/opt/netgen/bin
-e PYTHONPATH=/opt/netgen/lib/python3/dist-packages
netgen_${CI_BUILD_REF_NAME}_installed:${UBUNTU_VERSION}
bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V"'
test_ubuntu_1510: test_ubuntu_1510:
<<: *ubuntu_1510 <<: *ubuntu_1510

16
CMakeLists.txt
View File

@ -42,10 +42,9 @@ set(CPACK_PACKAGE_VERSION "${PACKAGE_VERSION}")
####################################################################### #######################################################################
if(USE_CCACHE) if(USE_CCACHE)
find_program(CCACHE_FOUND ccache) find_program(CCACHE_FOUND NAMES ccache ccache.bat)
if(CCACHE_FOUND) if(CCACHE_FOUND)
set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE ccache) set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE ${CCACHE_FOUND})
set_property(GLOBAL PROPERTY RULE_LAUNCH_LINK ccache)
endif(CCACHE_FOUND) endif(CCACHE_FOUND)
endif(USE_CCACHE) endif(USE_CCACHE)
@ -259,6 +258,11 @@ execute_process(COMMAND hdiutil create -volname Netgen -srcfolder ${INSTALL_DIR}
") ")
add_custom_target(bundle COMMAND ${CMAKE_COMMAND} "-P" "${CMAKE_CURRENT_BINARY_DIR}/fixup.cmake") add_custom_target(bundle COMMAND ${CMAKE_COMMAND} "-P" "${CMAKE_CURRENT_BINARY_DIR}/fixup.cmake")
####################################################################### #######################################################################
# CTest
enable_testing()
include(CTest)
#######################################################################
add_subdirectory(libsrc) add_subdirectory(libsrc)
add_subdirectory(ng) add_subdirectory(ng)
@ -270,13 +274,9 @@ add_subdirectory(py_tutorials)
add_subdirectory(doc) add_subdirectory(doc)
add_subdirectory(windows) add_subdirectory(windows)
add_subdirectory(nglib) add_subdirectory(nglib)
add_subdirectory(tests)
#######################################################################
# CTest
enable_testing()
include(CTest)
####################################################################### #######################################################################
# Debian packager # Debian packager

3
libsrc/csg/csg.hpp
View File

@ -10,7 +10,6 @@
#include <myadt.hpp> #include <myadt.hpp>
#include <gprim.hpp> #include <gprim.hpp>
#include <meshing.hpp> #include <meshing.hpp>
// #include <geometry2d.hpp> // #include <geometry2d.hpp>
#include "../gprim/spline.hpp" #include "../gprim/spline.hpp"
#include "../gprim/splinegeometry.hpp" #include "../gprim/splinegeometry.hpp"
@ -21,12 +20,12 @@
#include "solid.hpp" #include "solid.hpp"
#include "identify.hpp" #include "identify.hpp"
#include "singularref.hpp" #include "singularref.hpp"
#include "splinesurface.hpp"
#include "csgeom.hpp" #include "csgeom.hpp"
#include "csgparser.hpp" #include "csgparser.hpp"
#include "triapprox.hpp" #include "triapprox.hpp"
#include "algprim.hpp" #include "algprim.hpp"
#include "splinesurface.hpp"
#include "brick.hpp" #include "brick.hpp"
#include "spline3d.hpp" #include "spline3d.hpp"
#include "manifold.hpp" #include "manifold.hpp"

19
libsrc/csg/csgeom.cpp
View File

@ -345,7 +345,25 @@ namespace netgen
const ExtrusionFace * ef = dynamic_cast< const ExtrusionFace * > (GetSurface(i)); const ExtrusionFace * ef = dynamic_cast< const ExtrusionFace * > (GetSurface(i));
const RevolutionFace * rf = dynamic_cast< const RevolutionFace * > (GetSurface(i)); const RevolutionFace * rf = dynamic_cast< const RevolutionFace * > (GetSurface(i));
const DummySurface * dummyf = dynamic_cast< const DummySurface * > (GetSurface(i)); const DummySurface * dummyf = dynamic_cast< const DummySurface * > (GetSurface(i));
const SplineSurface * splines = dynamic_cast<const SplineSurface *> (GetSurface(i));
if (splines)
{
splines->GetBase()->GetPrimitiveData(classname,coeffs);
out << classname << " " << coeffs.Size() << "\n";
for (int j=0; j<coeffs.Size(); j++)
out << coeffs[j] << " ";
out << "\n";
for (auto cut : *(splines->GetCuts()))
{
cut->GetPrimitiveData(classname,coeffs);
out << classname << " " << coeffs.Size() << "\n";
for (int j=0; j<coeffs.Size(); j++)
out << coeffs[j] << " ";
out << "\n";
}
return;
}
if(sp) if(sp)
{ {
@ -819,6 +837,7 @@ namespace netgen
int inv; int inv;
int nsurf = GetNSurf(); int nsurf = GetNSurf();
isidenticto.SetSize(nsurf); isidenticto.SetSize(nsurf);
for (int i = 0; i < nsurf; i++) for (int i = 0; i < nsurf; i++)
isidenticto[i] = i; isidenticto[i] = i;

7
libsrc/csg/csgeom.hpp
View File

@ -143,6 +143,9 @@ namespace netgen
/// filename of inputfile /// filename of inputfile
string filename; string filename;
/// store splinesurfaces, such that added ones do not get deleted before geometry does
Array<shared_ptr<SplineSurface>> spline_surfaces;
public: public:
CSGeometry (); CSGeometry ();
CSGeometry (const string & afilename); CSGeometry (const string & afilename);
@ -315,7 +318,9 @@ namespace netgen
virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam); virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
virtual const Refinement & GetRefinement () const; virtual const Refinement & GetRefinement () const;
void AddSplineSurface (shared_ptr<SplineSurface> ss) { spline_surfaces.Append(ss); }
}; };

4
libsrc/csg/edgeflw.cpp
View File

@ -491,7 +491,7 @@ namespace netgen
if(splinesurface) if(splinesurface)
{ {
auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p);
mesh.SetCD2Name(refedges[i].edgenr,*name); mesh.SetCD2Name(refedges[i].edgenr,name);
} }
else else
{ {
@ -499,7 +499,7 @@ namespace netgen
if(splinesurface2) if(splinesurface2)
{ {
auto name = splinesurface2->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); auto name = splinesurface2->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p);
mesh.SetCD2Name(refedges[i].edgenr,*name); mesh.SetCD2Name(refedges[i].edgenr,name);
} }
} }

20
libsrc/csg/python_csg.cpp
View File

@ -207,11 +207,11 @@ DLL_HEADER void ExportCSG(py::module &m)
; ;
py::class_<SplineSurface, shared_ptr<SplineSurface>> (m, "SplineSurface", py::class_<SplineSurface, shared_ptr<SplineSurface>> (m, "SplineSurface",
"A surface for co dim 2 integrals on the splines") "A surface for co dim 2 integrals on the splines")
.def("__init__", FunctionPointer ([](SplineSurface* instance, shared_ptr<SPSolid> base, py::list cuts) .def("__init__", FunctionPointer ([](SplineSurface* instance, shared_ptr<SPSolid> base, py::list cuts)
{ {
auto primitive = dynamic_cast<OneSurfacePrimitive*> (base->GetSolid()->GetPrimitive()); auto primitive = dynamic_cast<OneSurfacePrimitive*> (base->GetSolid()->GetPrimitive());
auto acuts = new Array<OneSurfacePrimitive*>(); auto acuts = make_shared<Array<shared_ptr<OneSurfacePrimitive>>>();
for(int i = 0; i<py::len(cuts);i++) for(int i = 0; i<py::len(cuts);i++)
{ {
py::extract<shared_ptr<SPSolid>> sps(cuts[i]); py::extract<shared_ptr<SPSolid>> sps(cuts[i]);
@ -219,13 +219,14 @@ DLL_HEADER void ExportCSG(py::module &m)
throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!"); throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!");
auto sp = dynamic_cast<OneSurfacePrimitive*>(sps()->GetSolid()->GetPrimitive()); auto sp = dynamic_cast<OneSurfacePrimitive*>(sps()->GetSolid()->GetPrimitive());
if(sp) if(sp)
acuts->Append(sp); acuts->Append(shared_ptr<OneSurfacePrimitive>(sp));
else else
throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!"); throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!");
} }
if(!primitive) if(!primitive)
throw NgException("Base is not a SurfacePrimitive in constructor of SplineSurface!"); throw NgException("Base is not a SurfacePrimitive in constructor of SplineSurface!");
new (instance) SplineSurface(primitive,acuts); new (instance) SplineSurface(shared_ptr<OneSurfacePrimitive>(primitive),acuts);
py::object obj = py::cast(instance);
}),py::arg("base"), py::arg("cuts")=py::list()) }),py::arg("base"), py::arg("cuts")=py::list())
.def("AddPoint", FunctionPointer .def("AddPoint", FunctionPointer
([] (SplineSurface & self, double x, double y, double z, bool hpref) ([] (SplineSurface & self, double x, double y, double z, bool hpref)
@ -237,8 +238,8 @@ DLL_HEADER void ExportCSG(py::module &m)
.def("AddSegment", FunctionPointer .def("AddSegment", FunctionPointer
([] (SplineSurface & self, int i1, int i2, string bcname, double maxh) ([] (SplineSurface & self, int i1, int i2, string bcname, double maxh)
{ {
auto str = new string(bcname); auto seg = make_shared<LineSeg<3>>(self.GetPoint(i1),self.GetPoint(i2));
self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str,maxh); self.AppendSegment(seg,bcname,maxh);
}), }),
py::arg("pnt1"),py::arg("pnt2"),py::arg("bcname")="default", py::arg("maxh")=-1.) py::arg("pnt1"),py::arg("pnt2"),py::arg("bcname")="default", py::arg("maxh")=-1.)
; ;
@ -433,14 +434,15 @@ DLL_HEADER void ExportCSG(py::module &m)
([] (CSGeometry & self, shared_ptr<SplineSurface> surf) ([] (CSGeometry & self, shared_ptr<SplineSurface> surf)
{ {
auto cuttings = surf->CreateCuttingSurfaces(); auto cuttings = surf->CreateCuttingSurfaces();
auto spsol = make_shared<SPSolid>(new Solid(&*surf)); auto spsol = make_shared<SPSolid>(new Solid(surf.get()));
for(auto cut : (*cuttings)){ for(auto cut : (*cuttings)){
spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(cut))); spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(cut.get())));
} }
spsol->AddSurfaces(self); spsol->AddSurfaces(self);
int tlonr = self.SetTopLevelObject(spsol->GetSolid(), &*surf); int tlonr = self.SetTopLevelObject(spsol->GetSolid(), surf.get());
for(auto p : surf->GetPoints()) for(auto p : surf->GetPoints())
self.AddUserPoint(p); self.AddUserPoint(p);
self.AddSplineSurface(surf);
}), }),
py::arg("SplineSurface")) py::arg("SplineSurface"))

46
libsrc/csg/splinesurface.cpp
View File

@ -10,14 +10,14 @@ void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const
geompoints.Last().hpref = hpref; geompoints.Last().hpref = hpref;
} }
void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh) void SplineSurface :: AppendSegment(shared_ptr<SplineSeg<3>> sp, string & bcname, double amaxh)
{ {
splines.Append(spline); splines.Append(sp);
bcnames.Append(bcname); bcnames.Append(bcname);
maxh.Append(amaxh); maxh.Append(amaxh);
} }
string* SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const string SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const
{ {
double eps = 1e-5; double eps = 1e-5;
@ -32,36 +32,44 @@ void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const
return bcnames[i]; return bcnames[i];
} }
} }
return new string("default"); return "default";
} }
Array<OneSurfacePrimitive*>* SplineSurface :: CreateCuttingSurfaces() const const shared_ptr<Array<shared_ptr<OneSurfacePrimitive>>> SplineSurface :: CreateCuttingSurfaces()
{ {
auto cuttings = new Array<OneSurfacePrimitive*>(); if(all_cuts)
return all_cuts;
auto cuttings = make_shared<Array<shared_ptr<OneSurfacePrimitive>>>();
for (auto cut : *cuts) for (auto cut : *cuts)
cuttings->Append(cut); cuttings->Append(cut);
for(int i = 0; i<splines.Size(); i++) for(int i = 0; i<splines.Size(); i++)
{ {
auto spline = splines[i]; auto spline = splines[i];
auto lineseg = dynamic_cast<LineSeg<3>*>(spline); auto lineseg = dynamic_cast<LineSeg<3>*>(spline.get());
auto p1 = Point<3>(spline->GetPoint(0)); if(lineseg)
Project(p1); {
auto p2 = Point<3>(spline->GetPoint(1)); auto p1 = Point<3>(spline->GetPoint(0));
Project(p2); Project(p1);
auto vec = Vec<3>(p2)-Vec<3>(p1); auto p2 = Point<3>(spline->GetPoint(1));
auto plane = new Plane(p1,-Cross(vec,baseprimitive->GetNormalVector(p1))); Project(p2);
if(maxh[i]>0) auto vec = Vec<3>(p2)-Vec<3>(p1);
{ auto plane = make_shared<Plane>(p1,-Cross(vec,baseprimitive->GetNormalVector(p1)));
plane->SetMaxH(maxh[i]); if(maxh[i]>0)
} {
cuttings->Append(plane); plane->SetMaxH(maxh[i]);
}
cuttings->Append(plane);
}
else
throw NgException("Spline type not implemented for SplineSurface!");
} }
all_cuts = cuttings;
return cuttings; return cuttings;
} }
void SplineSurface :: Print(ostream & str) const void SplineSurface :: Print(ostream & str) const
{ {
str << "SplineSurface " << endl; str << "SplineSurface with base " << *baseprimitive << endl;
} }
} }

26
libsrc/csg/splinesurface.hpp
View File

@ -8,19 +8,20 @@ namespace netgen
{ {
protected: protected:
Array<GeomPoint<3>> geompoints; Array<GeomPoint<3>> geompoints;
Array<SplineSeg<3>*> splines; Array<shared_ptr<SplineSeg<3>>> splines;
Array<string*> bcnames; Array<string> bcnames;
Array<double> maxh; Array<double> maxh;
OneSurfacePrimitive* baseprimitive; shared_ptr<OneSurfacePrimitive> baseprimitive;
Array<OneSurfacePrimitive*>* cuts; shared_ptr<Array<shared_ptr<OneSurfacePrimitive>>> cuts;
shared_ptr<Array<shared_ptr<OneSurfacePrimitive>>> all_cuts;
public: public:
SplineSurface(OneSurfacePrimitive* abaseprimitive, Array<OneSurfacePrimitive*>* acuts) : SplineSurface(shared_ptr<OneSurfacePrimitive> abaseprimitive, shared_ptr<Array<shared_ptr<OneSurfacePrimitive>>> acuts) :
OneSurfacePrimitive(), baseprimitive(abaseprimitive), cuts(acuts) OneSurfacePrimitive(), baseprimitive(abaseprimitive), cuts(acuts)
{ ; } { ; }
virtual ~SplineSurface() { ; } virtual ~SplineSurface() { ; }
const Array<SplineSeg<3>*> & GetSplines() const { return splines; } const auto & GetSplines() const { return splines; }
int GetNSplines() const { return splines.Size(); } int GetNSplines() const { return splines.Size(); }
const Array<GeomPoint<3>>& GetPoints() const { return geompoints; } const Array<GeomPoint<3>>& GetPoints() const { return geompoints; }
string GetSplineType(const int i) const { return splines[i]->GetType(); } string GetSplineType(const int i) const { return splines[i]->GetType(); }
@ -28,16 +29,15 @@ namespace netgen
const SplineSeg<3> & GetSpline(const int i) const { return *splines[i]; } const SplineSeg<3> & GetSpline(const int i) const { return *splines[i]; }
int GetNP() const { return geompoints.Size(); } int GetNP() const { return geompoints.Size(); }
const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; } const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; }
string* GetBCName(int i) const { return bcnames[i]; } string GetBCName(int i) const { return bcnames[i]; }
string* GetBCNameOf(Point<3> p1, Point<3> p2) const; string GetBCNameOf(Point<3> p1, Point<3> p2) const;
DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false); DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false);
void AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh = -1); void AppendSegment(shared_ptr<SplineSeg<3>> spline, string & bcname, double amaxh = -1);
OneSurfacePrimitive* GetBase() const { return baseprimitive; }
Array<OneSurfacePrimitive*>* CreateCuttingSurfaces() const;
const shared_ptr<Array<shared_ptr<OneSurfacePrimitive>>> CreateCuttingSurfaces();
const shared_ptr<Array<shared_ptr<OneSurfacePrimitive>>> GetCuts() const { return all_cuts; }
const shared_ptr<OneSurfacePrimitive> GetBase() const { return baseprimitive; }
virtual void Project (Point<3> & p3d) const { baseprimitive->Project(p3d); } virtual void Project (Point<3> & p3d) const { baseprimitive->Project(p3d); }
virtual double CalcFunctionValue (const Point<3> & point) const virtual double CalcFunctionValue (const Point<3> & point) const

30
libsrc/general/array.hpp
View File

@ -267,6 +267,15 @@ namespace netgen
(*this)[i] = a2[i]; (*this)[i] = a2[i];
} }
/// array move
Array (Array && a2)
: FlatArray<T,BASE,TIND> (a2.size, a2.data), allocsize(a2.allocsize), ownmem(a2.ownmem)
{
a2.size = 0;
a2.data = nullptr;
a2.allocsize = 0;
a2.ownmem = false;
}
/// if responsible, deletes memory /// if responsible, deletes memory
@ -377,7 +386,15 @@ namespace netgen
return *this; return *this;
} }
Array & operator= (Array && a2)
{
Swap (data, a2.data);
Swap (size, a2.size);
Swap (allocsize, a2.allocsize);
Swap (ownmem, a2.ownmem);
return *this;
}
private: private:
/// resize array, at least to size minsize. copy contents /// resize array, at least to size minsize. copy contents
@ -391,7 +408,16 @@ namespace netgen
T * p = new T[nsize]; T * p = new T[nsize];
int mins = (nsize < size) ? nsize : size; int mins = (nsize < size) ? nsize : size;
memcpy (p, data, mins * sizeof(T)); // memcpy (p, data, mins * sizeof(T));
#if defined(__GNUG__) && __GNUC__ < 5 && !defined(__clang__)
for (size_t i = 0; i < mins; i++) p[i] = move(data[i]);
#else
if (std::is_trivially_copyable<T>::value)
memcpy (p, data, sizeof(T)*mins);
else
for (size_t i = 0; i < mins; i++) p[i] = move(data[i]);
#endif
if (ownmem) if (ownmem)
delete [] data; delete [] data;

16
libsrc/general/table.hpp
View File

@ -34,11 +34,25 @@ protected:
public: public:
/// ///
BASE_TABLE (BASE_TABLE && table2)
: data(move(table2.data)), oneblock(table2.oneblock)
{
table2.oneblock = nullptr;
}
BASE_TABLE (int size); BASE_TABLE (int size);
/// ///
BASE_TABLE (const FlatArray<int> & entrysizes, int elemsize); BASE_TABLE (const FlatArray<int> & entrysizes, int elemsize);
/// ///
~BASE_TABLE (); ~BASE_TABLE ();
BASE_TABLE & operator= (BASE_TABLE && table2)
{
data = move(table2.data);
Swap (oneblock, table2.oneblock);
return *this;
}
/// ///
void SetSize (int size); void SetSize (int size);
/// ///
@ -94,7 +108,7 @@ class TABLE : public BASE_TABLE
public: public:
/// Creates table. /// Creates table.
inline TABLE () : BASE_TABLE(0) { ; } inline TABLE () : BASE_TABLE(0) { ; }
/// Creates table of size size /// Creates table of size size
inline TABLE (int size) : BASE_TABLE (size) { ; } inline TABLE (int size) : BASE_TABLE (size) { ; }

23
libsrc/geom2d/genmesh2d.cpp
View File

@ -238,6 +238,29 @@ namespace netgen
} }
} }
// first add all vertices (for compatible orientation on periodic bnds)
{
double diam2 = Dist2(pmin, pmax);
for (int i = 0; i < splines.Size(); i++)
for (int j : { 1, 2 })
{
Point<2> hnewp = (j == 1) ? splines[i]->StartPI() : splines[i]->EndPI();
Point<3> newp(hnewp(0), hnewp(1), 0);
int layer = GetSpline(i).layer;
int npi = -1;
for (PointIndex pi = PointIndex::BASE;
pi < mesh2d.GetNP()+PointIndex::BASE; pi++)
if (Dist2 (mesh2d.Point(pi), newp) < 1e-12 * diam2 && mesh2d.Point(pi).GetLayer() == layer)
npi = pi;
if (npi == -1)
{
npi = mesh2d.AddPoint (newp, layer);
searchtree.Insert (newp, npi);
}
}
}
for (int i = 0; i < splines.Size(); i++) for (int i = 0; i < splines.Size(); i++)
if (GetSpline(i).copyfrom == -1) if (GetSpline(i).copyfrom == -1)
{ {

39
libsrc/include/nginterface_v2.hpp
View File

@ -36,54 +36,52 @@ namespace netgen
class Ng_Points class Ng_Points
{ {
public: public:
int num; size_t num;
const int * ptr; const int * ptr;
int Size() const { return num; } size_t Size() const { return num; }
int operator[] (int i) const { return ptr[i]-POINTINDEX_BASE; } int operator[] (size_t i) const { return ptr[i]-POINTINDEX_BASE; }
}; };
class Ng_Vertices class Ng_Vertices
{ {
public: public:
int num; size_t num;
const int * ptr; const int * ptr;
int Size() const { return num; } size_t Size() const { return num; }
int operator[] (int i) const { return ptr[i]-POINTINDEX_BASE; } int operator[] (size_t i) const { return ptr[i]-POINTINDEX_BASE; }
}; };
class Ng_Edges class Ng_Edges
{ {
public: public:
int num; size_t num;
const T_EDGE2 * ptr; const T_EDGE2 * ptr;
int Size() const { return num; } size_t Size() const { return num; }
// int operator[] (int i) const { return abs (ptr[i])-1; } int operator[] (size_t i) const { return ptr[i].nr; }
int operator[] (int i) const { return ptr[i].nr; }
}; };
class Ng_Faces class Ng_Faces
{ {
public: public:
int num; size_t num;
const T_FACE2 * ptr; const T_FACE2 * ptr;
int Size() const { return num; } size_t Size() const { return num; }
// int operator[] (int i) const { return (ptr[i]-1) / 8; } int operator[] (size_t i) const { return ptr[i].nr; }
int operator[] (int i) const { return ptr[i].nr; }
}; };
class Ng_Facets class Ng_Facets
{ {
public: public:
int num; size_t num;
const int * ptr; const int * ptr;
int Size() const { return num; } size_t Size() const { return num; }
int operator[] (int i) const { return ptr[i]; } int operator[] (size_t i) const { return ptr[i]; }
}; };
@ -229,11 +227,14 @@ namespace netgen
Ng_Point GetPoint (int nr) const; Ng_Point GetPoint (int nr) const;
template <int DIM> template <int DIM>
Ng_Element GetElement (int nr) const; Ng_Element GetElement (size_t nr) const;
template <int DIM> template <int DIM>
int GetElementIndex (int nr) const; int GetElementIndex (size_t nr) const;
/// material/boundary label of region, template argument is co-dimension
template <int DIM>
const string & GetMaterialCD (int region_nr) const;
/// Curved Elements: /// Curved Elements:
/// elnr .. element nr /// elnr .. element nr

31
libsrc/include/nginterface_v2_impl.hpp
View File

@ -5,13 +5,13 @@ NGX_INLINE DLL_HEADER Ng_Point Ngx_Mesh :: GetPoint (int nr) const
template <> template <>
NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (int nr) const NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (size_t nr) const
{ {
return (*mesh).pointelements[nr].index; return (*mesh).pointelements[nr].index;
} }
template <> template <>
NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (size_t nr) const
{ {
if(mesh->GetDimension()==3) if(mesh->GetDimension()==3)
return (*mesh)[SegmentIndex(nr)].edgenr; return (*mesh)[SegmentIndex(nr)].edgenr;
@ -20,21 +20,21 @@ NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const
} }
template <> template <>
NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<2> (int nr) const NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<2> (size_t nr) const
{ {
int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex(); int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex();
return mesh->GetFaceDescriptor(ind).BCProperty(); return mesh->GetFaceDescriptor(ind).BCProperty();
} }
template <> template <>
NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<3> (int nr) const NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<3> (size_t nr) const
{ {
return (*mesh)[ElementIndex(nr)].GetIndex(); return (*mesh)[ElementIndex(nr)].GetIndex();
} }
template <> template <>
NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (int nr) const NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (size_t nr) const
{ {
const Element0d & el = mesh->pointelements[nr]; const Element0d & el = mesh->pointelements[nr];
@ -60,7 +60,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (int nr) const
template <> template <>
NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const
{ {
const Segment & el = mesh->LineSegment (SegmentIndex(nr)); const Segment & el = mesh->LineSegment (SegmentIndex(nr));
@ -107,7 +107,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
} }
template <> template <>
NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (int nr) const NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (size_t nr) const
{ {
const Element2d & el = mesh->SurfaceElement (SurfaceElementIndex (nr)); const Element2d & el = mesh->SurfaceElement (SurfaceElementIndex (nr));
@ -146,7 +146,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (int nr) const
} }
template <> template <>
NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (int nr) const NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (size_t nr) const
{ {
const Element & el = mesh->VolumeElement (ElementIndex (nr)); const Element & el = mesh->VolumeElement (ElementIndex (nr));
@ -175,8 +175,23 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (int nr) const
template <> NGX_INLINE DLL_HEADER
const string & Ngx_Mesh :: GetMaterialCD<0> (int region_nr) const
{
return mesh->GetMaterial(region_nr+1);
}
template <> NGX_INLINE DLL_HEADER
const string & Ngx_Mesh :: GetMaterialCD<1> (int region_nr) const
{
return mesh->GetBCName(region_nr);
}
template <> NGX_INLINE DLL_HEADER
const string & Ngx_Mesh :: GetMaterialCD<2> (int region_nr) const
{
return mesh->GetCD2Name(region_nr);
}

23
libsrc/meshing/meshclass.cpp
View File

@ -7,7 +7,7 @@ namespace netgen
static mutex buildsearchtree_mutex; static mutex buildsearchtree_mutex;
Mesh :: Mesh () Mesh :: Mesh ()
: surfarea(*this) : surfarea(*this), topology(*this)
{ {
// volelements.SetName ("vol elements"); // volelements.SetName ("vol elements");
// surfelements.SetName ("surf elements"); // surfelements.SetName ("surf elements");
@ -27,7 +27,7 @@ namespace netgen
numvertices = -1; numvertices = -1;
dimension = 3; dimension = 3;
topology = new MeshTopology (*this); // topology = new MeshTopology (*this);
curvedelems = new CurvedElements (*this); curvedelems = new CurvedElements (*this);
clusters = new AnisotropicClusters (*this); clusters = new AnisotropicClusters (*this);
ident = new Identifications (*this); ident = new Identifications (*this);
@ -57,7 +57,7 @@ namespace netgen
delete segmentht; delete segmentht;
delete curvedelems; delete curvedelems;
delete clusters; delete clusters;
delete topology; // delete topology;
delete ident; delete ident;
delete elementsearchtree; delete elementsearchtree;
delete coarsemesh; delete coarsemesh;
@ -126,8 +126,9 @@ namespace netgen
delete ident; delete ident;
ident = new Identifications (*this); ident = new Identifications (*this);
delete topology; // delete topology;
topology = new MeshTopology (*this); // topology = new MeshTopology (*this);
topology = MeshTopology (*this);
delete curvedelems; delete curvedelems;
curvedelems = new CurvedElements (*this); curvedelems = new CurvedElements (*this);
delete clusters; delete clusters;
@ -1241,7 +1242,7 @@ namespace netgen
if (ntasks == 1) // sequential run only if (ntasks == 1) // sequential run only
{ {
topology -> Update(); topology.Update();
clusters -> Update(); clusters -> Update();
} }
@ -1481,7 +1482,7 @@ namespace netgen
CalcSurfacesOfNode (); CalcSurfacesOfNode ();
topology -> Update(); topology.Update();
clusters -> Update(); clusters -> Update();
SetNextMajorTimeStamp(); SetNextMajorTimeStamp();
@ -4878,12 +4879,12 @@ namespace netgen
//(*testout) << "velement " << velement << endl; //(*testout) << "velement " << velement << endl;
Array<int> faces; Array<int> faces;
topology->GetElementFaces(velement,faces); topology.GetElementFaces(velement,faces);
//(*testout) << "faces " << faces << endl; //(*testout) << "faces " << faces << endl;
for(int i=0; i<faces.Size(); i++) for(int i=0; i<faces.Size(); i++)
faces[i] = topology->GetFace2SurfaceElement(faces[i]); faces[i] = topology.GetFace2SurfaceElement(faces[i]);
//(*testout) << "surfel " << faces << endl; //(*testout) << "surfel " << faces << endl;
@ -4910,7 +4911,7 @@ namespace netgen
} }
Array<int> faces2; Array<int> faces2;
topology->GetElementFaces(velement,faces2); topology.GetElementFaces(velement,faces2);
/* /*
cout << "no matching surf element" << endl cout << "no matching surf element" << endl
<< "p = " << p << endl << "p = " << p << endl
@ -5709,7 +5710,7 @@ namespace netgen
void Mesh :: UpdateTopology (TaskManager tm) void Mesh :: UpdateTopology (TaskManager tm)
{ {
topology->Update(tm); topology.Update(tm);
clusters->Update(tm); clusters->Update(tm);
#ifdef PARALLEL #ifdef PARALLEL
if (paralleltop) if (paralleltop)

4
libsrc/meshing/meshclass.hpp
View File

@ -103,7 +103,7 @@ namespace netgen
mutable int elementsearchtreets; mutable int elementsearchtreets;
/// element -> face, element -> edge etc ... /// element -> face, element -> edge etc ...
MeshTopology * topology; MeshTopology topology;
/// methods for high order elements /// methods for high order elements
class CurvedElements * curvedelems; class CurvedElements * curvedelems;
@ -687,7 +687,7 @@ namespace netgen
const MeshTopology & GetTopology () const const MeshTopology & GetTopology () const
{ return *topology; } { return topology; }
DLL_HEADER void UpdateTopology (TaskManager tm = &DummyTaskManager); DLL_HEADER void UpdateTopology (TaskManager tm = &DummyTaskManager);

2
libsrc/meshing/meshfunc.cpp
View File

@ -75,7 +75,7 @@ namespace netgen
(*testout) << "Meshing subdomain " << k << endl; (*testout) << "Meshing subdomain " << k << endl;
mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k)); mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
mesh3d.CalcSurfacesOfNode(); mesh3d.CalcSurfacesOfNode();
mesh3d.FindOpenElements(k); mesh3d.FindOpenElements(k);

16
libsrc/meshing/meshtype.hpp
View File

@ -25,9 +25,25 @@ namespace netgen
HEX = 25 HEX = 25
}; };
/*
typedef int ELEMENT_EDGE[2]; // initial point, end point typedef int ELEMENT_EDGE[2]; // initial point, end point
typedef int ELEMENT_FACE[4]; // points, last one is -1 for trig typedef int ELEMENT_FACE[4]; // points, last one is -1 for trig
*/
struct ELEMENT_EDGE
{
int vals[2];
int & operator[] (size_t i) { return vals[i]; }
int operator[] (size_t i) const { return vals[i]; }
};
struct ELEMENT_FACE
{
int vals[4];
int & operator[] (size_t i) { return vals[i]; }
int operator[] (size_t i) const { return vals[i]; }
};
#define ELEMENT_MAXPOINTS 12 #define ELEMENT_MAXPOINTS 12
#define ELEMENT2D_MAXPOINTS 8 #define ELEMENT2D_MAXPOINTS 8

12
libsrc/meshing/parallelmesh.cpp
View File

@ -718,7 +718,7 @@ namespace netgen
NgProfiler::StopTimer (timerloc2); NgProfiler::StopTimer (timerloc2);
topology -> Update(); topology.Update();
clusters -> Update(); clusters -> Update();
// paralleltop -> UpdateCoarseGrid(); // paralleltop -> UpdateCoarseGrid();
@ -1218,7 +1218,7 @@ namespace netgen
int ne = GetNE(); int ne = GetNE();
int nn = GetNP(); int nn = GetNP();
int nedges = topology->GetNEdges(); int nedges = topology.GetNEdges();
idxtype *xadj, * adjacency, *v_weights = NULL, *e_weights = NULL; idxtype *xadj, * adjacency, *v_weights = NULL, *e_weights = NULL;
@ -1240,7 +1240,7 @@ namespace netgen
for ( int edge = 1; edge <= nedges; edge++ ) for ( int edge = 1; edge <= nedges; edge++ )
{ {
int v1, v2; int v1, v2;
topology->GetEdgeVertices ( edge, v1, v2); topology.GetEdgeVertices ( edge, v1, v2);
cnt[v1-1] ++; cnt[v1-1] ++;
cnt[v2-1] ++; cnt[v2-1] ++;
} }
@ -1257,7 +1257,7 @@ namespace netgen
for ( int edge = 1; edge <= nedges; edge++ ) for ( int edge = 1; edge <= nedges; edge++ )
{ {
int v1, v2; int v1, v2;
topology->GetEdgeVertices ( edge, v1, v2); topology.GetEdgeVertices ( edge, v1, v2);
adjacency[ xadj[v1-1] + cnt[v1-1] ] = v2-1; adjacency[ xadj[v1-1] + cnt[v1-1] ] = v2-1;
adjacency[ xadj[v2-1] + cnt[v2-1] ] = v1-1; adjacency[ xadj[v2-1] + cnt[v2-1] ] = v1-1;
cnt[v1-1]++; cnt[v1-1]++;
@ -1313,7 +1313,7 @@ namespace netgen
// int nn = GetNP(); // int nn = GetNP();
// int nedges = topology->GetNEdges(); // int nedges = topology->GetNEdges();
int nfaces = topology->GetNFaces(); int nfaces = topology.GetNFaces();
idxtype *xadj, * adjacency, *v_weights = NULL, *e_weights = NULL; idxtype *xadj, * adjacency, *v_weights = NULL, *e_weights = NULL;
@ -1340,7 +1340,7 @@ namespace netgen
for ( int el=1; el <= ne; el++ ) for ( int el=1; el <= ne; el++ )
{ {
Element volel = VolumeElement(el); Element volel = VolumeElement(el);
topology->GetElementFaces(el, elfaces); topology.GetElementFaces(el, elfaces);
for ( int i = 0; i < elfaces.Size(); i++ ) for ( int i = 0; i < elfaces.Size(); i++ )
{ {
if ( facevolels1[elfaces[i]-1] == -1 ) if ( facevolels1[elfaces[i]-1] == -1 )

154
libsrc/meshing/topology.cpp
View File

@ -44,23 +44,14 @@ namespace netgen
MeshTopology :: MeshTopology (const Mesh & amesh) MeshTopology :: MeshTopology (const Mesh & amesh)
: mesh(amesh) : mesh(&amesh)
{ {
buildedges = 1; buildedges = 1;
buildfaces = 1; buildfaces = 1;
vert2element = 0;
vert2surfelement = 0;
vert2segment = 0;
timestamp = -1; timestamp = -1;
} }
MeshTopology :: ~MeshTopology () MeshTopology :: ~MeshTopology () { ; }
{
delete vert2element;
delete vert2surfelement;
delete vert2segment;
delete vert2pointelement;
}
template <typename FUNC> template <typename FUNC>
void LoopOverEdges (const Mesh & mesh, MeshTopology & top, PointIndex v, void LoopOverEdges (const Mesh & mesh, MeshTopology & top, PointIndex v,
@ -350,13 +341,13 @@ namespace netgen
#endif #endif
if (timestamp > mesh.GetTimeStamp()) return; if (timestamp > mesh->GetTimeStamp()) return;
int ne = mesh.GetNE(); int ne = mesh->GetNE();
int nse = mesh.GetNSE(); int nse = mesh->GetNSE();
int nseg = mesh.GetNSeg(); int nseg = mesh->GetNSeg();
int np = mesh.GetNP(); int np = mesh->GetNP();
int nv = mesh.GetNV(); int nv = mesh->GetNV();
if (id == 0) if (id == 0)
PrintMessage (3, "Update mesh topology"); PrintMessage (3, "Update mesh topology");
@ -367,12 +358,7 @@ namespace netgen
(*testout) << "nseg = " << nseg << endl; (*testout) << "nseg = " << nseg << endl;
(*testout) << "np = " << np << endl; (*testout) << "np = " << np << endl;
(*testout) << "nv = " << nv << endl; (*testout) << "nv = " << nv << endl;
delete vert2element;
delete vert2surfelement;
delete vert2segment;
delete vert2pointelement;
Array<int,PointIndex::BASE> cnt(nv); Array<int,PointIndex::BASE> cnt(nv);
Array<int> vnums; Array<int> vnums;
@ -382,69 +368,67 @@ namespace netgen
vertex to surface element vertex to surface element
vertex to segment vertex to segment
*/ */
cnt = 0; cnt = 0;
for (ElementIndex ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
{ {
const Element & el = mesh[ei]; const Element & el = (*mesh)[ei];
for (int j = 0; j < el.GetNV(); j++) for (int j = 0; j < el.GetNV(); j++)
cnt[el[j]]++; cnt[el[j]]++;
} }
vert2element = new TABLE<ElementIndex,PointIndex::BASE> (cnt); vert2element = TABLE<ElementIndex,PointIndex::BASE> (cnt);
for (ElementIndex ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
{ {
const Element & el = mesh[ei]; const Element & el = (*mesh)[ei];
for (int j = 0; j < el.GetNV(); j++) for (int j = 0; j < el.GetNV(); j++)
vert2element->AddSave (el[j], ei); vert2element.AddSave (el[j], ei);
} }
cnt = 0; cnt = 0;
for (SurfaceElementIndex sei = 0; sei < nse; sei++) for (SurfaceElementIndex sei = 0; sei < nse; sei++)
{ {
const Element2d & el = mesh[sei]; const Element2d & el = (*mesh)[sei];
for (int j = 0; j < el.GetNV(); j++) for (int j = 0; j < el.GetNV(); j++)
cnt[el[j]]++; cnt[el[j]]++;
} }
vert2surfelement = new TABLE<SurfaceElementIndex,PointIndex::BASE> (cnt); vert2surfelement = TABLE<SurfaceElementIndex,PointIndex::BASE> (cnt);
for (SurfaceElementIndex sei = 0; sei < nse; sei++) for (SurfaceElementIndex sei = 0; sei < nse; sei++)
{ {
const Element2d & el = mesh[sei]; const Element2d & el = (*mesh)[sei];
for (int j = 0; j < el.GetNV(); j++) for (int j = 0; j < el.GetNV(); j++)
vert2surfelement->AddSave (el[j], sei); vert2surfelement.AddSave (el[j], sei);
} }
cnt = 0; cnt = 0;
for (SegmentIndex si = 0; si < nseg; si++) for (SegmentIndex si = 0; si < nseg; si++)
{ {
const Segment & seg = mesh.LineSegment(si); const Segment & seg = mesh->LineSegment(si);
cnt[seg[0]]++; cnt[seg[0]]++;
cnt[seg[1]]++; cnt[seg[1]]++;
} }
vert2segment = new TABLE<SegmentIndex,PointIndex::BASE> (cnt); vert2segment = TABLE<SegmentIndex,PointIndex::BASE> (cnt);
for (SegmentIndex si = 0; si < nseg; si++) for (SegmentIndex si = 0; si < nseg; si++)
{ {
const Segment & seg = mesh.LineSegment(si); const Segment & seg = mesh->LineSegment(si);
vert2segment->AddSave (seg[0], si); vert2segment.AddSave (seg[0], si);
vert2segment->AddSave (seg[1], si); vert2segment.AddSave (seg[1], si);
} }
cnt = 0; cnt = 0;
for (int pei = 0; pei < mesh.pointelements.Size(); pei++) for (int pei = 0; pei < mesh->pointelements.Size(); pei++)
{ {
const Element0d & pointel = mesh.pointelements[pei]; const Element0d & pointel = mesh->pointelements[pei];
cnt[pointel.pnum]++; cnt[pointel.pnum]++;
} }
vert2pointelement = new TABLE<int,PointIndex::BASE> (cnt); vert2pointelement = TABLE<int,PointIndex::BASE> (cnt);
for (int pei = 0; pei < mesh.pointelements.Size(); pei++) for (int pei = 0; pei < mesh->pointelements.Size(); pei++)
{ {
const Element0d & pointel = mesh.pointelements[pei]; const Element0d & pointel = mesh->pointelements[pei];
vert2pointelement->AddSave (pointel.pnum, pei); vert2pointelement.AddSave (pointel.pnum, pei);
} }
@ -477,15 +461,15 @@ namespace netgen
// ensure all coarse grid and intermediate level edges // ensure all coarse grid and intermediate level edges
cnt = 0; cnt = 0;
for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) for (int i = mesh->mlbetweennodes.Begin(); i < mesh->mlbetweennodes.End(); i++)
{ {
INDEX_2 parents = Sort (mesh.mlbetweennodes[i]); INDEX_2 parents = Sort (mesh->mlbetweennodes[i]);
if (parents[0] >= PointIndex::BASE) cnt[parents[0]]++; if (parents[0] >= PointIndex::BASE) cnt[parents[0]]++;
} }
TABLE<int,PointIndex::BASE> vert2vertcoarse (cnt); TABLE<int,PointIndex::BASE> vert2vertcoarse (cnt);
for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) for (int i = mesh->mlbetweennodes.Begin(); i < mesh->mlbetweennodes.End(); i++)
{ {
INDEX_2 parents = Sort (mesh.mlbetweennodes[i]); INDEX_2 parents = Sort (mesh->mlbetweennodes[i]);
if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]); if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]);
} }
@ -495,7 +479,7 @@ namespace netgen
for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++) for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++)
{ {
int onv = vert2edge[i].Size() + vert2vertcoarse[i].Size() + int onv = vert2edge[i].Size() + vert2vertcoarse[i].Size() +
4*(*vert2element)[i].Size() + 2*(*vert2surfelement)[i].Size() + (*vert2segment)[i].Size(); 4*(vert2element)[i].Size() + 2*(vert2surfelement)[i].Size() + (vert2segment)[i].Size();
max_edge_on_vertex = max (onv, max_edge_on_vertex); max_edge_on_vertex = max (onv, max_edge_on_vertex);
} }
@ -504,7 +488,7 @@ namespace netgen
cnt = 0; cnt = 0;
ParallelForRange ParallelForRange
(tm, mesh.Points().Size(), (tm, mesh->Points().Size(),
[&] (size_t begin, size_t end) [&] (size_t begin, size_t end)
{ {
INDEX_CLOSED_HASHTABLE<int> v2eht(2*max_edge_on_vertex+10); INDEX_CLOSED_HASHTABLE<int> v2eht(2*max_edge_on_vertex+10);
@ -527,7 +511,7 @@ namespace netgen
v2eht.Set (v2, 33); // some value v2eht.Set (v2, 33); // some value
} }
LoopOverEdges (mesh, *this, v, LoopOverEdges (*mesh, *this, v,
[&] (INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) [&] (INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir)
{ {
if (!v2eht.Used (edge.I2())) if (!v2eht.Used (edge.I2()))
@ -542,7 +526,7 @@ namespace netgen
// accumulate number of edges // accumulate number of edges
int ned = edge2vert.Size(); int ned = edge2vert.Size();
for (auto v : mesh.Points().Range()) for (auto v : mesh->Points().Range())
{ {
auto hv = cnt[v]; auto hv = cnt[v];
cnt[v] = ned; cnt[v] = ned;
@ -556,7 +540,7 @@ namespace netgen
// for (PointIndex v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) // for (PointIndex v = PointIndex::BASE; v < nv+PointIndex::BASE; v++)
ParallelForRange ParallelForRange
(tm, mesh.Points().Size(), (tm, mesh->Points().Size(),
[&] (size_t begin, size_t end) [&] (size_t begin, size_t end)
{ {
INDEX_CLOSED_HASHTABLE<int> v2eht(2*max_edge_on_vertex+10); INDEX_CLOSED_HASHTABLE<int> v2eht(2*max_edge_on_vertex+10);
@ -581,7 +565,7 @@ namespace netgen
vertex2.Append (v2); vertex2.Append (v2);
} }
LoopOverEdges (mesh, *this, v, LoopOverEdges (*mesh, *this, v,
[&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) [&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir)
{ {
if (!v2eht.Used(edge.I2())) if (!v2eht.Used(edge.I2()))
@ -600,7 +584,7 @@ namespace netgen
ned++; ned++;
} }
LoopOverEdges (mesh, *this, v, LoopOverEdges (*mesh, *this, v,
[&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) [&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir)
{ {
int edgenum = v2eht.Get(edge.I2()); int edgenum = v2eht.Get(edge.I2());
@ -661,7 +645,7 @@ namespace netgen
int max_face_on_vertex = 0; int max_face_on_vertex = 0;
for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++) for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++)
{ {
int onv = vert2oldface[i].Size() + (*vert2element)[i].Size() + (*vert2surfelement)[i].Size(); int onv = vert2oldface[i].Size() + vert2element[i].Size() + vert2surfelement[i].Size();
max_face_on_vertex = max (onv, max_face_on_vertex); max_face_on_vertex = max (onv, max_face_on_vertex);
} }
@ -680,7 +664,7 @@ namespace netgen
// for (auto v : mesh.Points().Range()) // for (auto v : mesh.Points().Range())
NgProfiler::StartTimer (timer2b1); NgProfiler::StartTimer (timer2b1);
ParallelForRange ParallelForRange
(tm, mesh.Points().Size(), (tm, mesh->Points().Size(),
[&] (size_t begin, size_t end) [&] (size_t begin, size_t end)
{ {
INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10); INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10);
@ -698,7 +682,7 @@ namespace netgen
vert2face.Set (face, 33); // something vert2face.Set (face, 33); // something
} }
int cnti = 0; int cnti = 0;
LoopOverFaces (mesh, *this, v, LoopOverFaces (*mesh, *this, v,
[&] (INDEX_4 i4, int elnr, int j, bool volume, int facedir) [&] (INDEX_4 i4, int elnr, int j, bool volume, int facedir)
{ {
INDEX_3 face(i4.I1(), i4.I2(), i4.I3()); INDEX_3 face(i4.I1(), i4.I2(), i4.I3());
@ -715,7 +699,7 @@ namespace netgen
// accumulate number of faces // accumulate number of faces
int nfa = oldnfa; int nfa = oldnfa;
for (auto v : mesh.Points().Range()) for (auto v : mesh->Points().Range())
{ {
auto hv = cnt[v]; auto hv = cnt[v];
cnt[v] = nfa; cnt[v] = nfa;
@ -726,7 +710,7 @@ namespace netgen
// for (auto v : mesh.Points().Range()) // for (auto v : mesh.Points().Range())
ParallelForRange ParallelForRange
(tm, mesh.Points().Size(), (tm, mesh->Points().Size(),
[&] (size_t begin, size_t end) [&] (size_t begin, size_t end)
{ {
INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10); INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10);
@ -746,7 +730,7 @@ namespace netgen
vert2face.Set (face, fnr); vert2face.Set (face, fnr);
} }
LoopOverFaces (mesh, *this, v, LoopOverFaces (*mesh, *this, v,
[&] (INDEX_4 i4, int elnr, int j, bool volume, int facedir) [&] (INDEX_4 i4, int elnr, int j, bool volume, int facedir)
{ {
INDEX_3 face(i4.I1(), i4.I2(), i4.I3()); INDEX_3 face(i4.I1(), i4.I2(), i4.I3());
@ -775,7 +759,7 @@ namespace netgen
} }
LoopOverFaces (mesh, *this, v, LoopOverFaces (*mesh, *this, v,
[&] (INDEX_4 i4, int elnr, int j, bool volume, int facedir) [&] (INDEX_4 i4, int elnr, int j, bool volume, int facedir)
{ {
INDEX_3 face(i4.I1(), i4.I2(), i4.I3()); INDEX_3 face(i4.I1(), i4.I2(), i4.I3());
@ -1149,10 +1133,10 @@ namespace netgen
{ {
// (*testout) << "is face of element " << vertels[k] << endl; // (*testout) << "is face of element " << vertels[k] << endl;
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]].hp_elnr];
(*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl; (*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl;
} }
@ -1306,7 +1290,7 @@ namespace netgen
void MeshTopology :: GetElementEdges (int elnr, Array<int> & eledges) const void MeshTopology :: GetElementEdges (int elnr, Array<int> & eledges) const
{ {
int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); int ned = GetNEdges (mesh->VolumeElement(elnr).GetType());
eledges.SetSize (ned); eledges.SetSize (ned);
for (int i = 0; i < ned; i++) for (int i = 0; i < ned; i++)
eledges[i] = edges.Get(elnr)[i].nr+1; eledges[i] = edges.Get(elnr)[i].nr+1;
@ -1314,7 +1298,7 @@ namespace netgen
} }
void MeshTopology :: GetElementFaces (int elnr, Array<int> & elfaces, bool withorientation) const void MeshTopology :: GetElementFaces (int elnr, Array<int> & elfaces, bool withorientation) const
{ {
int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType());
elfaces.SetSize (nfa); elfaces.SetSize (nfa);
if (!withorientation) if (!withorientation)
@ -1342,7 +1326,7 @@ namespace netgen
void MeshTopology :: GetElementEdgeOrientations (int elnr, Array<int> & eorient) const void MeshTopology :: GetElementEdgeOrientations (int elnr, Array<int> & eorient) const
{ {
int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); int ned = GetNEdges (mesh->VolumeElement(elnr).GetType());
eorient.SetSize (ned); eorient.SetSize (ned);
for (int i = 1; i <= ned; i++) for (int i = 1; i <= ned; i++)
// eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; // eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1;
@ -1352,7 +1336,7 @@ namespace netgen
void MeshTopology :: GetElementFaceOrientations (int elnr, Array<int> & forient) const void MeshTopology :: GetElementFaceOrientations (int elnr, Array<int> & forient) const
{ {
int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType());
forient.SetSize (nfa); forient.SetSize (nfa);
for (int i = 1; i <= nfa; i++) for (int i = 1; i <= nfa; i++)
// forient.Elem(i) = faces.Get(elnr)[i-1].forient; // forient.Elem(i) = faces.Get(elnr)[i-1].forient;
@ -1366,7 +1350,7 @@ namespace netgen
{ {
// int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); // int ned = GetNEdges (mesh.VolumeElement(elnr).GetType());
if (mesh.GetDimension()==3 || 1) if (mesh->GetDimension()==3 || 1)
{ {
if (orient) if (orient)
{ {
@ -1454,7 +1438,7 @@ namespace netgen
void MeshTopology :: GetSurfaceElementEdges (int elnr, Array<int> & eledges) const void MeshTopology :: GetSurfaceElementEdges (int elnr, Array<int> & eledges) const
{ {
int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); int ned = GetNEdges (mesh->SurfaceElement(elnr).GetType());
eledges.SetSize (ned); eledges.SetSize (ned);
for (int i = 0; i < ned; i++) for (int i = 0; i < ned; i++)
// eledges[i] = abs (surfedges.Get(elnr)[i]); // eledges[i] = abs (surfedges.Get(elnr)[i]);
@ -1463,7 +1447,7 @@ namespace netgen
void MeshTopology :: GetEdges (SurfaceElementIndex elnr, Array<int> & eledges) const void MeshTopology :: GetEdges (SurfaceElementIndex elnr, Array<int> & eledges) const
{ {
int ned = GetNEdges (mesh[elnr].GetType()); int ned = GetNEdges ( (*mesh)[elnr].GetType());
eledges.SetSize (ned); eledges.SetSize (ned);
for (int i = 0; i < ned; i++) for (int i = 0; i < ned; i++)
// eledges[i] = abs (surfedges[elnr][i])-1; // eledges[i] = abs (surfedges[elnr][i])-1;
@ -1486,7 +1470,7 @@ namespace netgen
void MeshTopology :: void MeshTopology ::
GetSurfaceElementEdgeOrientations (int elnr, Array<int> & eorient) const GetSurfaceElementEdgeOrientations (int elnr, Array<int> & eorient) const
{ {
int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); int ned = GetNEdges (mesh->SurfaceElement(elnr).GetType());
eorient.SetSize (ned); eorient.SetSize (ned);
for (int i = 0; i < ned; i++) for (int i = 0; i < ned; i++)
// eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1; // eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1;
@ -1504,7 +1488,7 @@ namespace netgen
int MeshTopology :: GetSurfaceElementEdges (int elnr, int * eledges, int * orient) const int MeshTopology :: GetSurfaceElementEdges (int elnr, int * eledges, int * orient) const
{ {
int i; int i;
if (mesh.GetDimension() == 3 || 1) if (mesh->GetDimension() == 3 || 1)
{ {
if (orient) if (orient)
{ {
@ -1554,7 +1538,7 @@ namespace netgen
int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const
{ {
const Element & el = mesh.VolumeElement (elnr); const Element & el = mesh->VolumeElement (elnr);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());
int k = locedgenr; int k = locedgenr;
@ -1566,7 +1550,7 @@ namespace netgen
int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const
{ {
const Element & el = mesh.VolumeElement (elnr); const Element & el = mesh->VolumeElement (elnr);
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType()); const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType());
@ -1622,7 +1606,7 @@ namespace netgen
int MeshTopology :: GetSurfaceElementEdgeOrientation (int elnr, int locedgenr) const int MeshTopology :: GetSurfaceElementEdgeOrientation (int elnr, int locedgenr) const
{ {
const Element2d & el = mesh.SurfaceElement (elnr); const Element2d & el = mesh->SurfaceElement (elnr);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());
int k = locedgenr; int k = locedgenr;
@ -1633,7 +1617,7 @@ namespace netgen
int MeshTopology :: GetSurfaceElementFaceOrientation2 (int elnr) const int MeshTopology :: GetSurfaceElementFaceOrientation2 (int elnr) const
{ {
const Element2d & el = mesh.SurfaceElement (elnr); const Element2d & el = mesh->SurfaceElement (elnr);
const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType()); const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType());
@ -1687,7 +1671,7 @@ namespace netgen
int MeshTopology :: GetSegmentEdgeOrientation (int elnr) const int MeshTopology :: GetSegmentEdgeOrientation (int elnr) const
{ {
const Segment & el = mesh.LineSegment (elnr); const Segment & el = mesh->LineSegment (elnr);
const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());
int k = 0; int k = 0;
@ -1760,7 +1744,7 @@ namespace netgen
// find one element having all vertices of the face // find one element having all vertices of the face
for (int i = 0; i < els.Size(); i++) for (int i = 0; i < els.Size(); i++)
{ {
const Element & el = mesh[els[i]]; const Element & el = (*mesh)[els[i]];
int nref_faces = GetNFaces (el.GetType()); int nref_faces = GetNFaces (el.GetType());
const ELEMENT_FACE * ref_faces = GetFaces1 (el.GetType()); const ELEMENT_FACE * ref_faces = GetFaces1 (el.GetType());
int nfa_ref_edges = GetNEdges (GetFaceType(fnr)); int nfa_ref_edges = GetNEdges (GetFaceType(fnr));
@ -1849,12 +1833,12 @@ namespace netgen
void MeshTopology :: GetVertexElements (int vnr, Array<ElementIndex> & elements) const void MeshTopology :: GetVertexElements (int vnr, Array<ElementIndex> & elements) const
{ {
if (vert2element) if (vert2element.Size())
{ {
int ne = vert2element->EntrySize(vnr); int ne = vert2element.EntrySize(vnr);
elements.SetSize(ne); elements.SetSize(ne);
for (int i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
elements.Elem(i) = vert2element->Get(vnr, i); elements.Elem(i) = vert2element.Get(vnr, i);
} }
} }
@ -1884,13 +1868,13 @@ namespace netgen
void MeshTopology :: GetVertexSurfaceElements( int vnr, void MeshTopology :: GetVertexSurfaceElements( int vnr,
Array<SurfaceElementIndex> & elements ) const Array<SurfaceElementIndex> & elements ) const
{ {
if (vert2surfelement) if (vert2surfelement.Size())
{ {
int i; int i;
int ne = vert2surfelement->EntrySize(vnr); int ne = vert2surfelement.EntrySize(vnr);
elements.SetSize(ne); elements.SetSize(ne);
for (i = 1; i <= ne; i++) for (i = 1; i <= ne; i++)
elements.Elem(i) = vert2surfelement->Get(vnr, i); elements.Elem(i) = vert2surfelement.Get(vnr, i);
} }
} }

91
libsrc/meshing/topology.hpp
View File

@ -27,33 +27,52 @@ struct T_FACE
int fnr; // 0-based int fnr; // 0-based
}; };
template <typename T, int S>
struct FixArray
{
T vals[S];
T & operator[] (size_t i) { return vals[i]; }
T operator[] (size_t i) const { return vals[i]; }
};
class MeshTopology class MeshTopology
{ {
const Mesh & mesh; const Mesh * mesh;
bool buildedges; bool buildedges;
bool buildfaces; bool buildfaces;
Array<INDEX_2> edge2vert; Array<INDEX_2> edge2vert;
Array<INDEX_4> face2vert; Array<INDEX_4> face2vert;
/*
Array<T_EDGE[12]> edges; Array<T_EDGE[12]> edges;
Array<T_FACE[6]> faces; Array<T_FACE[6]> faces;
Array<T_EDGE[4]> surfedges; Array<T_EDGE[4]> surfedges;
*/
Array<FixArray<T_EDGE,12>> edges;
Array<FixArray<T_FACE,6>> faces;
Array<FixArray<T_EDGE,4>> surfedges;
Array<T_EDGE> segedges; Array<T_EDGE> segedges;
Array<T_FACE> surffaces; Array<T_FACE> surffaces;
Array<INDEX_2> surf2volelement; Array<INDEX_2> surf2volelement;
Array<int> face2surfel; Array<int> face2surfel;
TABLE<ElementIndex,PointIndex::BASE> *vert2element; TABLE<ElementIndex,PointIndex::BASE> vert2element;
TABLE<SurfaceElementIndex,PointIndex::BASE> *vert2surfelement; TABLE<SurfaceElementIndex,PointIndex::BASE> vert2surfelement;
TABLE<SegmentIndex,PointIndex::BASE> *vert2segment; TABLE<SegmentIndex,PointIndex::BASE> vert2segment;
TABLE<int,PointIndex::BASE> *vert2pointelement = nullptr; TABLE<int,PointIndex::BASE> vert2pointelement;
int timestamp; int timestamp;
public: public:
int GetNSurfedges() const {return surfedges.Size();} int GetNSurfedges() const {return surfedges.Size();}
MeshTopology () = default;
MeshTopology (const MeshTopology & top) = default;
MeshTopology (MeshTopology && top) = default;
MeshTopology (const Mesh & amesh); MeshTopology (const Mesh & amesh);
~MeshTopology (); ~MeshTopology ();
MeshTopology & operator= (const MeshTopology & top) = default;
MeshTopology & operator= (MeshTopology && top) = default;
void SetBuildEdges (bool be) void SetBuildEdges (bool be)
{ buildedges = be; } { buildedges = be; }
void SetBuildFaces (bool bf) void SetBuildFaces (bool bf)
@ -145,17 +164,17 @@ public:
void GetVertexElements (int vnr, Array<ElementIndex> & elements) const; void GetVertexElements (int vnr, Array<ElementIndex> & elements) const;
FlatArray<ElementIndex> GetVertexElements (int vnr) const FlatArray<ElementIndex> GetVertexElements (int vnr) const
{ return (*vert2element)[vnr]; } { return vert2element[vnr]; }
void GetVertexSurfaceElements( int vnr, Array<SurfaceElementIndex>& elements ) const; void GetVertexSurfaceElements( int vnr, Array<SurfaceElementIndex>& elements ) const;
FlatArray<SurfaceElementIndex> GetVertexSurfaceElements (int vnr) const FlatArray<SurfaceElementIndex> GetVertexSurfaceElements (int vnr) const
{ return (*vert2surfelement)[vnr]; } { return vert2surfelement[vnr]; }
FlatArray<SegmentIndex> GetVertexSegments (int vnr) const FlatArray<SegmentIndex> GetVertexSegments (int vnr) const
{ return (*vert2segment)[vnr]; } { return vert2segment[vnr]; }
FlatArray<int> GetVertexPointElements (int vnr) const FlatArray<int> GetVertexPointElements (int vnr) const
{ return (*vert2pointelement)[vnr]; } { return vert2pointelement[vnr]; }
int GetVerticesEdge ( int v1, int v2) const; int GetVerticesEdge ( int v1, int v2) const;
void GetSegmentVolumeElements ( int segnr, Array<ElementIndex> & els ) const; void GetSegmentVolumeElements ( int segnr, Array<ElementIndex> & els ) const;
@ -337,22 +356,22 @@ inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et) const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
{ {
static int segm_edges[1][2] = static ELEMENT_EDGE segm_edges[1] =
{ { 1, 2 }}; { { 1, 2 }};
static int trig_edges[3][2] = static ELEMENT_EDGE trig_edges[3] =
{ { 3, 1 }, { { 3, 1 },
{ 2, 3 }, { 2, 3 },
{ 1, 2 }}; { 1, 2 }};
static int quad_edges[4][2] = static ELEMENT_EDGE quad_edges[4] =
{ { 1, 2 }, { { 1, 2 },
{ 3, 4 }, { 3, 4 },
{ 4, 1 }, { 4, 1 },
{ 2, 3 }}; { 2, 3 }};
static int tet_edges[6][2] = static ELEMENT_EDGE tet_edges[6] =
{ { 4, 1 }, { { 4, 1 },
{ 4, 2 }, { 4, 2 },
{ 4, 3 }, { 4, 3 },
@ -360,7 +379,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
{ 1, 3 }, { 1, 3 },
{ 2, 3 }}; { 2, 3 }};
static int prism_edges[9][2] = static ELEMENT_EDGE prism_edges[9] =
{ { 3, 1 }, { { 3, 1 },
{ 1, 2 }, { 1, 2 },
{ 3, 2 }, { 3, 2 },
@ -371,7 +390,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
{ 1, 4 }, { 1, 4 },
{ 2, 5 }}; { 2, 5 }};
static int pyramid_edges[8][2] = static ELEMENT_EDGE pyramid_edges[8] =
{ { 1, 2 }, { { 1, 2 },
{ 2, 3 }, { 2, 3 },
{ 1, 4 }, { 1, 4 },
@ -381,7 +400,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
{ 3, 5 }, { 3, 5 },
{ 4, 5 }}; { 4, 5 }};
static int hex_edges[12][2] = static ELEMENT_EDGE hex_edges[12] =
{ {
{ 1, 2 }, { 1, 2 },
{ 3, 4 }, { 3, 4 },
@ -435,22 +454,22 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et) const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
{ {
static int segm_edges[1][2] = static ELEMENT_EDGE segm_edges[1] =
{ { 0, 1 }}; { { 0, 1 }};
static int trig_edges[3][2] = static ELEMENT_EDGE trig_edges[3] =
{ { 2, 0 }, { { 2, 0 },
{ 1, 2 }, { 1, 2 },
{ 0, 1 }}; { 0, 1 }};
static int quad_edges[4][2] = static ELEMENT_EDGE quad_edges[4] =
{ { 0, 1 }, { { 0, 1 },
{ 2, 3 }, { 2, 3 },
{ 3, 0 }, { 3, 0 },
{ 1, 2 }}; { 1, 2 }};
static int tet_edges[6][2] = static ELEMENT_EDGE tet_edges[6] =
{ { 3, 0 }, { { 3, 0 },
{ 3, 1 }, { 3, 1 },
{ 3, 2 }, { 3, 2 },
@ -458,7 +477,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
{ 0, 2 }, { 0, 2 },
{ 1, 2 }}; { 1, 2 }};
static int prism_edges[9][2] = static ELEMENT_EDGE prism_edges[9] =
{ { 2, 0 }, { { 2, 0 },
{ 0, 1 }, { 0, 1 },
{ 2, 1 }, { 2, 1 },
@ -469,7 +488,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
{ 0, 3 }, { 0, 3 },
{ 1, 4 }}; { 1, 4 }};
static int pyramid_edges[8][2] = static ELEMENT_EDGE pyramid_edges[8] =
{ { 0, 1 }, { { 0, 1 },
{ 1, 2 }, { 1, 2 },
{ 0, 3 }, { 0, 3 },
@ -479,7 +498,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
{ 2, 4 }, { 2, 4 },
{ 3, 4 }}; { 3, 4 }};
static int hex_edges[12][2] = static ELEMENT_EDGE hex_edges[12] =
{ {
{ 0, 1 }, { 0, 1 },
{ 2, 3 }, { 2, 3 },
@ -540,18 +559,18 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
{ {
static const int trig_faces[1][4] = static const ELEMENT_FACE trig_faces[1] =
{ { 1, 2, 3, 0 } }; { { 1, 2, 3, 0 } };
static const int quad_faces[1][4] = static const ELEMENT_FACE quad_faces[1] =
{ { 1, 2, 3, 4 } }; { { 1, 2, 3, 4 } };
static const int tet_faces[4][4] = static const ELEMENT_FACE tet_faces[4] =
{ { 4, 2, 3, 0 }, { { 4, 2, 3, 0 },
{ 4, 3, 1, 0 }, { 4, 3, 1, 0 },
{ 4, 1, 2, 0 }, { 4, 1, 2, 0 },
{ 1, 3, 2, 0 } }; { 1, 3, 2, 0 } };
static const int prism_faces[5][4] = static const ELEMENT_FACE prism_faces[5] =
{ {
{ 1, 3, 2, 0 }, { 1, 3, 2, 0 },
{ 4, 5, 6, 0 }, { 4, 5, 6, 0 },
@ -560,7 +579,7 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
{ 2, 3, 6, 5 } { 2, 3, 6, 5 }
}; };
static const int pyramid_faces[5][4] = static const ELEMENT_FACE pyramid_faces[5] =
{ {
{ 1, 2, 5, 0 }, { 1, 2, 5, 0 },
{ 2, 3, 5, 0 }, { 2, 3, 5, 0 },
@ -569,7 +588,7 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
{ 1, 4, 3, 2 } { 1, 4, 3, 2 }
}; };
static const int hex_faces[6][4] = static const ELEMENT_FACE hex_faces[6] =
{ {
{ 1, 4, 3, 2 }, { 1, 4, 3, 2 },
{ 5, 6, 7, 8 }, { 5, 6, 7, 8 },
@ -622,18 +641,18 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et) inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
{ {
static const int trig_faces[1][4] = static const ELEMENT_FACE trig_faces[1] =
{ { 0, 1, 2, -1 } }; { { 0, 1, 2, -1 } };
static const int quad_faces[1][4] = static const ELEMENT_FACE quad_faces[1] =
{ { 0, 1, 2, 3 } }; { { 0, 1, 2, 3 } };
static const int tet_faces[4][4] = static const ELEMENT_FACE tet_faces[4] =
{ { 3, 1, 2, -1 }, { { 3, 1, 2, -1 },
{ 3, 2, 0, -1 }, { 3, 2, 0, -1 },
{ 3, 0, 1, -1 }, { 3, 0, 1, -1 },
{ 0, 2, 1, -1 } }; { 0, 2, 1, -1 } };
static const int prism_faces[5][4] = static const ELEMENT_FACE prism_faces[5] =
{ {
{ 0, 2, 1, -1 }, { 0, 2, 1, -1 },
{ 3, 4, 5, -1 }, { 3, 4, 5, -1 },
@ -642,7 +661,7 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
{ 1, 2, 5, 4 } { 1, 2, 5, 4 }
}; };
static const int pyramid_faces[5][4] = static const ELEMENT_FACE pyramid_faces[5] =
{ {
{ 0, 1, 4, -1 }, { 0, 1, 4, -1 },
{ 1, 2, 4, -1 }, { 1, 2, 4, -1 },
@ -651,7 +670,7 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
{ 0, 3, 2, 1 } { 0, 3, 2, 1 }
}; };
static const int hex_faces[6][4] = static const ELEMENT_FACE hex_faces[6] =
{ {
{ 0, 3, 2, 1 }, { 0, 3, 2, 1 },
{ 4, 5, 6, 7 }, { 4, 5, 6, 7 },

2
ng/CMakeLists.txt
View File

@ -5,6 +5,8 @@ if(USE_GUI)
# add icon to netgen executable # add icon to netgen executable
enable_language(RC) enable_language(RC)
set(netgen_sources ${netgen_sources} ../windows/netgen.rc) set(netgen_sources ${netgen_sources} ../windows/netgen.rc)
# Don't use ccache here due to incompatiblity with the resource compiler
set_directory_properties(PROPERTIES RULE_LAUNCH_COMPILE "")
endif(WIN32) endif(WIN32)
add_executable(netgen ${netgen_sources}) add_executable(netgen ${netgen_sources})

1
tests/CMakeLists.txt Normal file
View File

@ -0,0 +1 @@
add_subdirectory(pytest)

2
tests/docker_15.10
View File

@ -1,4 +1,4 @@
FROM ubuntu:15.10 FROM ubuntu:15.10
MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at> MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev libboost-python-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy
ADD . /root/src/netgen ADD . /root/src/netgen

2
tests/docker_16.04
View File

@ -1,4 +1,4 @@
FROM ubuntu:16.04 FROM ubuntu:16.04
MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at> MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev libboost-python-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy
ADD . /root/src/netgen ADD . /root/src/netgen

4
tests/pytest/CMakeLists.txt Normal file
View File

@ -0,0 +1,4 @@
if(USE_PYTHON)
add_test(NAME pytest COMMAND ${PYTHON_EXECUTABLE} -m pytest WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
add_custom_target(pytest ${PYTHON_EXECUTABLE} -m pytest WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
endif(USE_PYTHON)

52
tests/pytest/test_savemesh.py Normal file
View File

@ -0,0 +1,52 @@
from netgen.csg import *
from netgen import meshing
import filecmp
import difflib
from math import sqrt, cos, sin
def CreateQuad():
base = Plane(Pnt(0,0,0),Vec(0,0,1))
surface = SplineSurface(base)
pts = [(-0.2,-0.2,0),(-0.2,0.2,0),(0.2,0.2,0),(0.2,-0.2,0)]
geopts = [surface.AddPoint(*p) for p in pts]
for p1,p2,bc in [(0,1,"wire"), (1, 2,"contact"),(2,3,"wire"),(3,0,"wire")]:
surface.AddSegment(geopts[p1],geopts[p2],bc)
return surface
Cross = lambda a,b: [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-b[0]*a[1]]
def CreateGeo():
geo = CSGeometry()
air = OrthoBrick(Pnt(-1,-1,-1),Pnt(1,1,1))
geo.Add(air.mat("air"))
surface = CreateQuad()
geo.AddSplineSurface(surface)
return geo
def test_BBNDsave():
mesh = CreateGeo().GenerateMesh(maxh=0.4,perfstepsend = meshing.MeshingStep.MESHSURFACE)
for i in range(2):
mesh.GenerateVolumeMesh(mp = MeshingParameters(only3D_domain=i+1,maxh=0.4))
mesh.SetGeometry(None)
mesh.Save("test.vol")
mesh2 = meshing.Mesh()
mesh2.Load("test.vol")
mesh2.Save("test2.vol")
with open("test.vol","r") as f:
first = f.readlines()
with open("test2.vol","r") as f:
second = f.readlines()
# exclude the face colours section (because they aren't in the same order)
for i,line in enumerate(first):
if line[0:12] == "face_colours":
first = first[0:i]
second = second[0:i]
break
diff = difflib.context_diff(first,second)
print("File diff:")
l = list(diff)
print(*l)
assert len(l)==0