netgen/libsrc/meshing/python_mesh.cpp

571 lines
21 KiB
C++
Raw Normal View History

2014-08-31 19:05:24 +06:00
#ifdef NG_PYTHON
2014-08-30 06:15:59 +06:00
#include <boost/python.hpp>
2014-08-31 15:14:18 +06:00
#include <boost/python/slice.hpp>
2014-10-08 21:48:48 +06:00
#include <../general/ngpython.hpp>
2014-08-30 06:15:59 +06:00
#include <mystdlib.h>
#include "meshing.hpp"
using namespace netgen;
namespace bp = boost::python;
2014-12-19 19:03:36 +05:00
namespace netgen
{
extern shared_ptr<NetgenGeometry> ng_geometry;
}
2014-08-30 06:15:59 +06:00
2014-08-31 15:14:18 +06:00
template <typename T, int BASE = 0, typename TIND = int>
2014-08-30 06:15:59 +06:00
void ExportArray ()
{
string name = string("Array_") + typeid(T).name();
2014-08-31 15:14:18 +06:00
bp::class_<Array<T,BASE,TIND>,boost::noncopyable>(name.c_str())
.def ("__len__", &Array<T,BASE,TIND>::Size)
2014-08-30 06:15:59 +06:00
.def ("__getitem__",
2014-08-31 15:14:18 +06:00
FunctionPointer ([](Array<T,BASE,TIND> & self, TIND i) -> T&
2014-08-30 06:15:59 +06:00
{
2014-08-31 15:14:18 +06:00
if (i < BASE || i >= BASE+self.Size())
2014-08-30 06:15:59 +06:00
bp::exec("raise IndexError()\n");
return self[i];
}),
bp::return_value_policy<bp::reference_existing_object>())
2014-08-31 15:14:18 +06:00
.def ("__iter__",
bp::range (FunctionPointer([](Array<T,BASE,TIND> & self) { return &self[BASE]; }),
FunctionPointer([](Array<T,BASE,TIND> & self) { return &self[BASE+self.Size()]; })))
2014-08-30 06:15:59 +06:00
;
}
2015-10-22 20:26:43 +05:00
void TranslateException (const NgException & ex)
{
string err = string("Netgen exception: ")+ex.What();
PyErr_SetString(PyExc_RuntimeError, err.c_str());
}
2014-08-30 06:15:59 +06:00
2015-01-16 15:29:25 +05:00
DLL_HEADER void ExportNetgenMeshing()
2014-08-30 06:15:59 +06:00
{
2014-10-08 21:48:48 +06:00
ModuleScope module("meshing");
2014-08-30 06:15:59 +06:00
2015-10-22 20:26:43 +05:00
bp::register_exception_translator<NgException>(&TranslateException);
2014-08-31 15:14:18 +06:00
bp::class_<PointIndex>("PointId", bp::init<int>())
2014-08-30 06:15:59 +06:00
.def("__repr__", &ToString<PointIndex>)
.def("__str__", &ToString<PointIndex>)
.add_property("nr", &PointIndex::operator int)
2015-05-20 13:39:55 +05:00
.def("__eq__" , FunctionPointer( [](PointIndex &self, PointIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](PointIndex &self ) { return static_cast<int>(self); }) )
2014-08-30 06:15:59 +06:00
;
2014-09-26 02:23:31 +06:00
2015-05-18 19:19:38 +05:00
bp::class_<ElementIndex>("ElementId3D", bp::init<int>())
.def("__repr__", &ToString<ElementIndex>)
.def("__str__", &ToString<ElementIndex>)
.add_property("nr", &ElementIndex::operator int)
2015-05-20 13:39:55 +05:00
.def("__eq__" , FunctionPointer( [](ElementIndex &self, ElementIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](ElementIndex &self ) { return static_cast<int>(self); }) )
2015-05-18 19:19:38 +05:00
;
bp::class_<SurfaceElementIndex>("ElementId2D", bp::init<int>())
.def("__repr__", &ToString<SurfaceElementIndex>)
.def("__str__", &ToString<SurfaceElementIndex>)
.add_property("nr", &SurfaceElementIndex::operator int)
2015-05-20 13:39:55 +05:00
.def("__eq__" , FunctionPointer( [](SurfaceElementIndex &self, SurfaceElementIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](SurfaceElementIndex &self ) { return static_cast<int>(self); }) )
2015-05-18 19:19:38 +05:00
;
bp::class_<SegmentIndex>("ElementId1D", bp::init<int>())
.def("__repr__", &ToString<SegmentIndex>)
.def("__str__", &ToString<SegmentIndex>)
.add_property("nr", &SegmentIndex::operator int)
2015-05-20 13:39:55 +05:00
.def("__eq__" , FunctionPointer( [](SegmentIndex &self, SegmentIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](SegmentIndex &self ) { return static_cast<int>(self); }) )
2015-05-18 19:19:38 +05:00
;
2014-09-26 02:23:31 +06:00
/*
2014-08-31 15:14:18 +06:00
bp::class_<Point<3>> ("Point")
.def(bp::init<double,double,double>())
;
2014-09-26 02:23:31 +06:00
*/
bp::class_<MeshPoint /* ,bp::bases<Point<3>> */ >("MeshPoint")
2015-05-18 19:19:38 +05:00
.def(bp::init<Point<3>>())
.def("__str__", &ToString<MeshPoint>)
.def("__repr__", &ToString<MeshPoint>)
2014-08-31 15:14:18 +06:00
.add_property("p", FunctionPointer([](const MeshPoint & self)
{
bp::list l;
2014-11-24 21:23:34 +05:00
l.append ( self[0] );
l.append ( self[1] );
l.append ( self[2] );
return bp::tuple(l);
2014-08-31 15:14:18 +06:00
}))
;
2015-01-20 22:41:16 +05:00
2014-08-30 06:15:59 +06:00
bp::class_<Element>("Element3D")
2015-05-18 19:19:38 +05:00
.def("__init__", bp::make_constructor
(FunctionPointer ([](int index, bp::list vertices)
{
2015-09-24 03:35:18 +05:00
if (bp::len(vertices) == 4)
{
Element * tmp = new Element(TET);
for (int i = 0; i < 4; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
tmp->SetIndex(index);
return tmp;
}
if (bp::len(vertices) == 8)
{
Element * tmp = new Element(HEX);
for (int i = 0; i < 8; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
tmp->SetIndex(index);
return tmp;
}
2015-05-18 19:19:38 +05:00
}),
bp::default_call_policies(), // need it to use arguments
(bp::arg("index")=1,bp::arg("vertices"))),
"create volume element"
)
.def("__repr__", &ToString<Element>)
2014-08-31 18:12:31 +06:00
.add_property("index", &Element::GetIndex, &Element::SetIndex)
2014-08-30 06:15:59 +06:00
.add_property("vertices",
FunctionPointer ([](const Element & self) -> bp::list
{
bp::list li;
for (int i = 0; i < self.GetNV(); i++)
li.append (self[i]);
return li;
}))
;
2014-10-01 19:16:34 +06:00
2015-01-20 22:41:16 +05:00
bp::class_<Element2d>("Element2D")
2015-05-18 19:19:38 +05:00
.def("__init__", bp::make_constructor
(FunctionPointer ([](int index, bp::list vertices)
{
2015-09-24 03:35:18 +05:00
if (bp::len(vertices) == 3)
{
Element2d * tmp = new Element2d(TRIG);
for (int i = 0; i < 3; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
tmp->SetIndex(index);
return tmp;
}
else
{
Element2d * tmp = new Element2d(QUAD);
for (int i = 0; i < 4; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
tmp->SetIndex(index);
return tmp;
}
2015-05-18 19:19:38 +05:00
}),
bp::default_call_policies(), // need it to use arguments
(bp::arg("index")=1,bp::arg("vertices"))),
"create surface element"
)
2015-01-20 22:41:16 +05:00
.add_property("index", &Element2d::GetIndex, &Element2d::SetIndex)
.add_property("vertices",
FunctionPointer([](const Element2d & self) -> bp::list
{
bp::list li;
for (int i = 0; i < self.GetNV(); i++)
li.append(self[i]);
return li;
}))
;
2015-05-18 19:19:38 +05:00
bp::class_<Segment>("Element1D")
.def("__init__", bp::make_constructor
2015-08-31 20:41:26 +05:00
(FunctionPointer ([](bp::list vertices, bp::list surfaces, int index)
2015-05-18 19:19:38 +05:00
{
Segment * tmp = new Segment;
for (int i = 0; i < 2; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
2015-08-31 20:41:26 +05:00
tmp -> si = index;
if (len(surfaces))
{
tmp->surfnr1 = bp::extract<int>(surfaces[0]);
tmp->surfnr2 = bp::extract<int>(surfaces[1]);
}
2015-05-18 19:19:38 +05:00
return tmp;
}),
bp::default_call_policies(),
2015-08-31 20:41:26 +05:00
(bp::arg("vertices"),
bp::arg("surfaces")=bp::list(),
bp::arg("index")=1
)),
2015-05-18 19:19:38 +05:00
"create segment element"
)
2015-11-11 22:46:18 +05:00
.def("__repr__", &ToString<Segment>)
2015-05-18 19:19:38 +05:00
.add_property("vertices",
FunctionPointer ([](const Segment & self) -> bp::list
{
bp::list li;
for (int i = 0; i < 2; i++)
li.append (self[i]);
return li;
}))
.add_property("surfaces",
FunctionPointer ([](const Segment & self) -> bp::list
{
bp::list li;
li.append (self.surfnr1);
li.append (self.surfnr2);
return li;
}))
;
2015-11-11 22:46:18 +05:00
bp::class_<Element0d>("Element0D")
.def("__init__", bp::make_constructor
(FunctionPointer ([](PointIndex vertex, int index)
{
Element0d * tmp = new Element0d;
tmp->pnum = vertex;
tmp->index = index;
return tmp;
}),
bp::default_call_policies(),
(bp::arg("vertex"),
bp::arg("index")=1
)),
"create point element"
)
.def("__repr__", &ToString<Element0d>)
.add_property("vertices",
FunctionPointer ([](const Element0d & self) -> bp::list
{
bp::list li;
li.append (self.pnum);
return li;
}))
;
2015-05-18 19:19:38 +05:00
bp::class_<FaceDescriptor>("FaceDescriptor")
.def(bp::init<const FaceDescriptor&>())
2015-08-29 16:54:00 +05:00
.def("__init__", bp::make_constructor
2015-09-01 13:50:15 +05:00
(FunctionPointer ([](int surfnr, int domin, int domout, int bc)
2015-08-29 16:54:00 +05:00
{
auto fd = new FaceDescriptor();
fd->SetSurfNr(surfnr);
fd->SetDomainIn(domin);
fd->SetDomainOut(domout);
2015-09-01 13:50:15 +05:00
fd->SetBCProperty(bc);
2015-08-29 16:54:00 +05:00
return fd;
}),
bp::default_call_policies(), // need it to use arguments
(bp::arg("surfnr")=1,
bp::arg("domin")=1,
2015-09-01 13:50:15 +05:00
bp::arg("domout")=0,
bp::arg("bc")=0
)),
2015-08-29 16:54:00 +05:00
"create facedescriptor")
2015-05-18 19:19:38 +05:00
.def("__str__", &ToString<FaceDescriptor>)
.def("__repr__", &ToString<FaceDescriptor>)
.add_property("surfnr", &FaceDescriptor::SurfNr, &FaceDescriptor::SetSurfNr)
2015-08-29 16:54:00 +05:00
.add_property("domin", &FaceDescriptor::DomainIn, &FaceDescriptor::SetDomainIn)
.add_property("domout", &FaceDescriptor::DomainOut, &FaceDescriptor::SetDomainOut)
2015-09-01 13:50:15 +05:00
.add_property("bc", &FaceDescriptor::BCProperty, &FaceDescriptor::SetBCProperty)
.add_property("bcname", FunctionPointer ([](FaceDescriptor & self) -> string { return self.GetBCName(); }))
2015-05-18 19:19:38 +05:00
;
2014-08-30 06:15:59 +06:00
ExportArray<Element>();
2014-08-31 15:14:18 +06:00
ExportArray<Element2d>();
2015-05-18 19:19:38 +05:00
ExportArray<Segment>();
2015-11-11 22:46:18 +05:00
ExportArray<Element0d>();
2014-08-31 15:14:18 +06:00
ExportArray<MeshPoint,PointIndex::BASE,PointIndex>();
2015-05-18 19:19:38 +05:00
ExportArray<FaceDescriptor>();
2014-08-30 06:15:59 +06:00
;
2016-01-29 19:29:03 +05:00
#if BOOST_VERSION >= 106000
2016-01-29 19:04:57 +05:00
bp::register_ptr_to_python<shared_ptr<Mesh>>();
#endif
2014-12-11 21:06:48 +05:00
bp::class_<Mesh,shared_ptr<Mesh>,boost::noncopyable>("Mesh", bp::no_init)
2015-10-22 19:32:58 +05:00
// .def(bp::init<>("create empty mesh"))
.def("__init__", bp::make_constructor
(FunctionPointer ([](int dim)
{
auto mesh = make_shared<Mesh>();
mesh->SetDimension(dim);
return mesh;
}),
bp::default_call_policies(), // need it to use named arguments
(
bp::arg("dim")=3
)
))
2014-08-30 06:15:59 +06:00
.def("__str__", &ToString<Mesh>)
2014-12-19 19:03:36 +05:00
.def("Load", FunctionPointer
([](Mesh & self, const string & filename)
{
2015-08-31 20:41:26 +05:00
istream * infile;
if (filename.find(".vol.gz") != string::npos)
infile = new igzstream (filename.c_str());
else
infile = new ifstream (filename.c_str());
// ifstream input(filename);
2016-02-27 00:35:09 +05:00
#ifdef PARALLEL
// int id;
MPI_Comm_rank(MPI_COMM_WORLD, &id);
MPI_Comm_size(MPI_COMM_WORLD, &ntasks);
if (id == 0)
{
self.Load(*infile);
self.Distribute();
}
else
{
self.SendRecvMesh();
}
#else
2015-08-31 20:41:26 +05:00
self.Load(*infile);
2016-02-27 00:35:09 +05:00
#endif
2014-12-19 19:03:36 +05:00
for (int i = 0; i < geometryregister.Size(); i++)
{
2015-08-31 20:41:26 +05:00
NetgenGeometry * hgeom = geometryregister[i]->LoadFromMeshFile (*infile);
2014-12-19 19:03:36 +05:00
if (hgeom)
{
ng_geometry.reset (hgeom);
break;
}
}
}))
// static_cast<void(Mesh::*)(const string & name)>(&Mesh::Load))
2014-08-30 06:15:59 +06:00
.def("Save", static_cast<void(Mesh::*)(const string & name)const>(&Mesh::Save))
2015-08-31 20:41:26 +05:00
.add_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
2014-08-30 06:15:59 +06:00
.def("Elements3D",
2014-09-03 15:07:10 +06:00
static_cast<Array<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
2014-08-30 06:15:59 +06:00
bp::return_value_policy<bp::reference_existing_object>())
2014-08-31 15:14:18 +06:00
.def("Elements2D",
2014-09-03 15:07:10 +06:00
static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements),
2014-08-31 15:14:18 +06:00
bp::return_value_policy<bp::reference_existing_object>())
2015-05-18 19:19:38 +05:00
.def("Elements1D",
static_cast<Array<Segment>&(Mesh::*)()> (&Mesh::LineSegments),
bp::return_value_policy<bp::reference_existing_object>())
2015-11-11 22:46:18 +05:00
.def("Elements0D", FunctionPointer([] (Mesh & self) -> Array<Element0d>&
{
return self.pointelements;
} ),
bp::return_value_policy<bp::reference_existing_object>())
2015-05-18 19:19:38 +05:00
2014-08-31 15:14:18 +06:00
.def("Points",
2014-09-03 15:07:10 +06:00
static_cast<Mesh::T_POINTS&(Mesh::*)()> (&Mesh::Points),
2014-08-30 06:15:59 +06:00
bp::return_value_policy<bp::reference_existing_object>())
2015-05-18 19:19:38 +05:00
.def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor),
bp::return_value_policy<bp::reference_existing_object>())
2015-12-11 18:31:28 +05:00
.def("GetNFaceDescriptors", &Mesh::GetNFD)
2014-08-30 06:15:59 +06:00
2014-08-31 15:14:18 +06:00
.def("__getitem__", FunctionPointer ([](const Mesh & self, PointIndex pi)
{
return self[pi];
}))
.def ("Add", FunctionPointer ([](Mesh & self, MeshPoint p)
{
return self.AddPoint (Point3d(p));
}))
2015-05-18 19:19:38 +05:00
.def ("Add", FunctionPointer ([](Mesh & self, const Element & el)
{
return self.AddVolumeElement (el);
}))
.def ("Add", FunctionPointer ([](Mesh & self, const Element2d & el)
{
return self.AddSurfaceElement (el);
}))
.def ("Add", FunctionPointer ([](Mesh & self, const Segment & el)
{
return self.AddSegment (el);
}))
2015-11-11 22:46:18 +05:00
.def ("Add", FunctionPointer ([](Mesh & self, const Element0d & el)
{
return self.pointelements.Append (el);
}))
2015-05-18 19:19:38 +05:00
.def ("Add", FunctionPointer ([](Mesh & self, const FaceDescriptor & fd)
{
return self.AddFaceDescriptor (fd);
}))
2015-09-01 13:50:15 +05:00
.def ("SetBCName", &Mesh::SetBCName)
2015-09-01 22:21:52 +05:00
.def ("GetBCName", FunctionPointer([](Mesh & self, int bc)->string
{ return self.GetBCName(bc); }))
.def ("SetMaterial", &Mesh::SetMaterial)
.def ("GetMaterial", FunctionPointer([](Mesh & self, int domnr)
{ return string(self.GetMaterial(domnr)); }))
2015-08-29 16:54:00 +05:00
.def ("GenerateVolumeMesh", FunctionPointer
([](Mesh & self)
{
cout << "generate vol mesh" << endl;
MeshingParameters mp;
2015-12-11 18:31:28 +05:00
mp.optsteps3d = 5;
2015-08-29 16:54:00 +05:00
MeshVolume (mp, self);
2015-09-02 22:01:49 +05:00
OptimizeVolume (mp, self);
2015-08-29 16:54:00 +05:00
}))
2015-12-11 18:31:28 +05:00
.def ("OptimizeVolumeMesh", FunctionPointer
([](Mesh & self)
{
MeshingParameters mp;
mp.optsteps3d = 5;
OptimizeVolume (mp, self);
}))
2015-08-29 16:54:00 +05:00
.def ("Refine", FunctionPointer
([](Mesh & self)
{
2016-02-26 16:19:10 +05:00
if (self.GetGeometry())
self.GetGeometry()->GetRefinement().Refine(self);
else
Refinement().Refine(self);
2015-08-29 16:54:00 +05:00
}))
2016-05-06 00:26:33 +05:00
.def ("BuildSearchTree", &Mesh::BuildElementSearchTree)
2014-12-18 19:00:58 +05:00
.def ("BoundaryLayer", FunctionPointer
([](Mesh & self, int bc, bp::list thicknesses, int volnr, bp::list materials)
{
int n = bp::len(thicknesses);
BoundaryLayerParameters blp;
for (int i = 1; i <= self.GetNFD(); i++)
if (self.GetFaceDescriptor(i).BCProperty() == bc)
blp.surfid.Append (i);
cout << "add layer at surfaces: " << blp.surfid << endl;
blp.prismlayers = n;
blp.growthfactor = 1.0;
// find max domain nr
int maxind = 0;
for (ElementIndex ei = 0; ei < self.GetNE(); ei++)
maxind = max (maxind, self[ei].GetIndex());
cout << "maxind = " << maxind << endl;
for ( int i=0; i<n; i++ )
{
blp.heights.Append( bp::extract<double>(thicknesses[i])()) ;
blp.new_matnrs.Append( maxind+1+i );
self.SetMaterial (maxind+1+i, bp::extract<string>(materials[i])().c_str());
}
blp.bulk_matnr = volnr;
GenerateBoundaryLayer (self, blp);
}
))
.def ("BoundaryLayer", FunctionPointer
2015-01-09 02:18:33 +05:00
([](Mesh & self, int bc, double thickness, int volnr, string material)
2014-12-18 19:00:58 +05:00
{
BoundaryLayerParameters blp;
2014-12-18 21:46:54 +05:00
for (int i = 1; i <= self.GetNFD(); i++)
if (self.GetFaceDescriptor(i).BCProperty() == bc)
blp.surfid.Append (i);
cout << "add layer at surfaces: " << blp.surfid << endl;
2014-12-18 19:00:58 +05:00
blp.prismlayers = 1;
blp.hfirst = thickness;
blp.growthfactor = 1.0;
2014-12-18 21:46:54 +05:00
// find max domain nr
int maxind = 0;
for (ElementIndex ei = 0; ei < self.GetNE(); ei++)
maxind = max (maxind, self[ei].GetIndex());
cout << "maxind = " << maxind << endl;
self.SetMaterial (maxind+1, material.c_str());
2015-01-09 02:18:33 +05:00
blp.new_matnr = maxind+1;
blp.bulk_matnr = volnr;
2014-12-18 19:00:58 +05:00
GenerateBoundaryLayer (self, blp);
}
))
2014-08-31 15:14:18 +06:00
;
2014-08-30 06:15:59 +06:00
2014-08-31 15:14:18 +06:00
typedef MeshingParameters MP;
2014-10-05 22:54:16 +06:00
bp::class_<MP> ("MeshingParameters", bp::init<>())
.def("__init__", bp::make_constructor
2015-11-25 15:25:03 +05:00
(FunctionPointer ([](double maxh, bool quad_dominated)
2014-10-05 22:54:16 +06:00
{
auto tmp = new MeshingParameters;
tmp->maxh = maxh;
2015-11-25 15:25:03 +05:00
tmp->quad = int(quad_dominated);
2014-10-05 22:54:16 +06:00
return tmp;
}),
bp::default_call_policies(), // need it to use arguments
2015-11-25 15:25:03 +05:00
(
bp::arg("maxh")=1000,
bp::arg("quad_dominated")=false
)),
2014-10-05 22:54:16 +06:00
"create meshing parameters"
2015-11-25 15:25:03 +05:00
)
2014-08-31 15:14:18 +06:00
.def("__str__", &ToString<MP>)
.add_property("maxh",
FunctionPointer ([](const MP & mp ) { return mp.maxh; }),
FunctionPointer ([](MP & mp, double maxh) { return mp.maxh = maxh; }))
2014-08-30 06:15:59 +06:00
;
2014-12-05 21:00:01 +05:00
bp::def("SetTestoutFile", FunctionPointer ([] (const string & filename)
{
delete testout;
testout = new ofstream (filename);
}));
2014-12-09 19:58:26 +05:00
bp::def("SetMessageImportance", FunctionPointer ([] (int importance)
{
int old = printmessage_importance;
printmessage_importance = importance;
return old;
}));
2014-08-30 06:15:59 +06:00
}
BOOST_PYTHON_MODULE(libmesh) {
ExportNetgenMeshing();
}
2014-08-31 19:05:24 +06:00
#endif
2014-08-30 06:15:59 +06:00