1 #ifdef NG_PYTHON
2 
3 #include <../general/ngpython.hpp>
4 #include <core/python_ngcore.hpp>
5 #include <csg.hpp>
6 #include "../meshing/python_mesh.hpp"
7 
8 
9 using namespace netgen;
10 using namespace pybind11::literals;
11 
12 namespace netgen
13 {
14   extern shared_ptr<NetgenGeometry> ng_geometry;
15 }
16 
17 
18 
19 // a shadow solid tree using shared pointers.
20 
21 class SPSolid
22 {
23   shared_ptr<SPSolid> s1, s2;
24   Solid * solid;
25   int bc = -1;
26   string bcname = "";
27   double maxh = -1;
28   string material;
29   bool owner;
30   double red = 0, green = 0, blue = 1;
31   bool transp = false;
32 public:
33   enum optyp { TERM, SECTION, UNION, SUB, EXISTING };
34 
SPSolid(Solid * as)35   SPSolid (Solid * as) : solid(as), owner(true), op(TERM) { ; }
SPSolid(Solid * as,int)36   SPSolid (Solid * as, int /*dummy*/)
37     : solid(as), owner(false), op(EXISTING) { ; }
~SPSolid()38   ~SPSolid ()
39   {
40     ; // if (owner) delete solid;
41   }
SPSolid(optyp aop,shared_ptr<SPSolid> as1,shared_ptr<SPSolid> as2)42   SPSolid (optyp aop, shared_ptr<SPSolid> as1, shared_ptr<SPSolid> as2)
43     : s1(as1), s2(as2), owner(true), op(aop)
44   {
45     if (aop == UNION)
46       solid = new Solid (Solid::UNION, s1->GetSolid(), s2->GetSolid());
47     else if (aop == SECTION)
48       solid = new Solid (Solid::SECTION, s1->GetSolid(), s2->GetSolid());
49     else if (aop == SUB)
50       solid = new Solid (Solid::SUB, s1->GetSolid()); // , s2->GetSolid());
51   }
52 
GetSolid()53   Solid * GetSolid() { return solid; }
GetSolid() const54   const Solid * GetSolid() const { return solid; }
55 
GiveUpOwner()56   void GiveUpOwner()
57   {
58     owner = false;
59     if (s1) s1 -> GiveUpOwner();
60     if (s2) s2 -> GiveUpOwner();
61   }
62 
AddSurfaces(CSGeometry & geom)63   void AddSurfaces(CSGeometry & geom)
64   {
65     if (op == TERM)
66       geom.AddSurfaces (solid->GetPrimitive());
67     if (s1) s1 -> AddSurfaces (geom);
68     if (s2) s2 -> AddSurfaces (geom);
69   }
70 
SetMaterial(string mat)71   void SetMaterial (string mat)  { material = mat; }
72 
GetMaterial()73   string GetMaterial ()
74   {
75     if (!material.empty()) return material;
76     if (s1)
77       {
78         string s1mat = s1->GetMaterial();
79         if (!s1mat.empty()) return s1mat;
80       }
81     if (s2)
82       {
83         string s2mat = s2->GetMaterial();
84         if (!s2mat.empty()) return s2mat;
85       }
86     return material;
87   }
88 
SetBC(int abc)89   void SetBC(int abc)
90   {
91     if (bc == -1)
92       {
93         bc = abc;
94         if (s1) s1 -> SetBC(bc);
95         if (s2) s2 -> SetBC(bc);
96         if (op == TERM)
97           {
98             Primitive * prim = solid -> GetPrimitive();
99             for (int i = 0; i < prim->GetNSurfaces(); i++)
100               prim->GetSurface(i).SetBCProperty (abc);
101             // cout << "set " << prim->GetNSurfaces() << " surfaces to bc " << bc << endl;
102           }
103       }
104   }
105 
SetBCName(string name)106   void SetBCName(string name)
107   {
108     if (bcname == "")
109       {
110         bcname = name;
111         if (s1) s1 -> SetBCName(name);
112         if (s2) s2 -> SetBCName(name);
113         if (op == TERM)
114           {
115             Primitive * prim = solid -> GetPrimitive();
116             for (int i = 0; i < prim->GetNSurfaces(); i++)
117               prim->GetSurface(i).SetBCName (name);
118             // cout << "set " << prim->GetNSurfaces() << " surfaces to bc " << bc << endl;
119           }
120       }
121   }
122 
123 
124 
SetMaxH(double amaxh)125   void SetMaxH(double amaxh)
126   {
127     if (maxh == -1)
128       {
129         maxh = amaxh;
130         if (s1) s1 -> SetMaxH(maxh);
131         if (s2) s2 -> SetMaxH(maxh);
132         if (op == TERM)
133           {
134             Primitive * prim = solid -> GetPrimitive();
135             for (int i = 0; i < prim->GetNSurfaces(); i++)
136               prim->GetSurface(i).SetMaxH (maxh);
137           }
138       }
139   }
140 
SetColor(double ared,double agreen,double ablue)141   void SetColor(double ared, double agreen, double ablue)
142   {
143     red = ared;
144     green = agreen;
145     blue = ablue;
146   }
147 
GetRed() const148   double GetRed() const { return red; }
GetGreen() const149   double GetGreen() const { return green; }
GetBlue() const150   double GetBlue() const { return blue; }
151 
SetTransparent()152   void SetTransparent() { transp = true; }
IsTransparent()153   bool IsTransparent() { return transp; }
154 
155 private:
156   optyp op;
157 };
158 
operator <<(ostream & ost,const SPSolid & sol)159 inline ostream & operator<< (ostream & ost, const SPSolid & sol)
160 {
161   ost << *sol.GetSolid();
162   return ost;
163 }
164 
165 namespace netgen
166 {
167   extern CSGeometry * ParseCSG (istream & istr, CSGeometry *instance=nullptr);
168 }
169 
170 
171 
ExportCSG(py::module & m)172 DLL_HEADER void ExportCSG(py::module &m)
173 {
174   py::class_<SplineGeometry<2>, shared_ptr<SplineGeometry<2>>>
175     (m, "SplineCurve2d")
176     .def(py::init<>())
177     .def ("AddPoint", FunctionPointer
178           ([] (SplineGeometry<2> & self, double x, double y)
179            {
180              self.geompoints.Append (GeomPoint<2> (Point<2> (x,y)));
181              return self.geompoints.Size()-1;
182            }))
183     .def ("AddSegment", [] (SplineGeometry<2> & self, int i1, int i2,
184                             string bcname, double maxh)
185     {
186       self.splines.Append (new LineSeg<2> (self.geompoints[i1], self.geompoints[i2], maxh, bcname));
187     }, "p1"_a, "p2"_a, "bcname"_a="default", "maxh"_a=1e99)
188     .def ("AddSegment", [] (SplineGeometry<2> & self, int i1, int i2,
189                             int i3, string bcname, double maxh)
190     {
191       self.splines.Append (new SplineSeg3<2> (self.geompoints[i1], self.geompoints[i2], self.geompoints[i3], bcname, maxh));
192     }, "p1"_a, "p2"_a, "p3"_a, "bcname"_a="default", "maxh"_a=1e99)
193     ;
194 
195   py::class_<SplineGeometry<3>,shared_ptr<SplineGeometry<3>>> (m,"SplineCurve3d")
196     .def(py::init<>())
197     .def ("AddPoint", FunctionPointer
198           ([] (SplineGeometry<3> & self, double x, double y, double z)
199            {
200              self.geompoints.Append (GeomPoint<3> (Point<3> (x,y,z)));
201              return self.geompoints.Size()-1;
202            }))
203     .def ("AddSegment", FunctionPointer
204           ([] (SplineGeometry<3> & self, int i1, int i2)
205            {
206              self.splines.Append (new LineSeg<3> (self.geompoints[i1], self.geompoints[i2]));
207            }))
208     .def ("AddSegment", FunctionPointer
209           ([] (SplineGeometry<3> & self, int i1, int i2, int i3)
210            {
211              self.splines.Append (new SplineSeg3<3> (self.geompoints[i1], self.geompoints[i2], self.geompoints[i3]));
212            }))
213     ;
214 
215   py::class_<SplineSurface, shared_ptr<SplineSurface>> (m, "SplineSurface",
216                         "A surface for co dim 2 integrals on the splines")
217     .def(py::init([](shared_ptr<SPSolid> base, py::list cuts)
218 	     {
219 	       auto primitive = dynamic_cast<OneSurfacePrimitive*> (base->GetSolid()->GetPrimitive());
220 	       auto acuts = make_shared<NgArray<shared_ptr<OneSurfacePrimitive>>>();
221 	       for(int i = 0; i<py::len(cuts);i++)
222 		 {
223 		   py::extract<shared_ptr<SPSolid>> sps(cuts[i]);
224 		   if(!sps.check())
225 		     throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!");
226 		   auto sp = dynamic_cast<OneSurfacePrimitive*>(sps()->GetSolid()->GetPrimitive());
227 		   if(sp)
228 		     acuts->Append(shared_ptr<OneSurfacePrimitive>(sp));
229 		   else
230 		     throw Exception("Cut must be SurfacePrimitive in constructor of SplineSurface!");
231 		 }
232 	       if(!primitive)
233 		 throw Exception("Base is not a SurfacePrimitive in constructor of SplineSurface!");
234 	       return make_shared<SplineSurface>(shared_ptr<OneSurfacePrimitive>(primitive),acuts);
235 	     }),py::arg("base"), py::arg("cuts")=py::list())
236     .def("AddPoint", FunctionPointer
237 	 ([] (SplineSurface & self, double x, double y, double z, bool hpref)
238 	  {
239 	    self.AppendPoint(Point<3>(x,y,z),hpref);
240 	    return self.GetNP()-1;
241 	  }),
242 	 py::arg("x"),py::arg("y"),py::arg("z"),py::arg("hpref")=false)
243     .def("AddSegment", [] (SplineSurface & self, int i1, int i2, string bcname, double maxh)
244 	  {
245             auto seg = make_shared<LineSeg<3>>(self.GetPoint(i1),self.GetPoint(i2));
246 	    self.AppendSegment(seg,bcname,maxh);
247 	  },
248 	 py::arg("pnt1"),py::arg("pnt2"),py::arg("bcname")="default", py::arg("maxh")=-1.)
249     .def("AddSegment", [] (SplineSurface& self, int i1, int i2, int i3, string bcname, double maxh)
250          {
251            auto seg = make_shared<SplineSeg3<3>>(self.GetPoint(i1), self.GetPoint(i2), self.GetPoint(i3));
252            self.AppendSegment(seg, bcname, maxh);
253          }, py::arg("pnt1"),py::arg("pnt2"), py::arg("pnt3"),py::arg("bcname")="default", py::arg("maxh")=-1.)
254     ;
255 
256   py::class_<SPSolid, shared_ptr<SPSolid>> (m, "Solid")
257     .def ("__str__", &ToString<SPSolid>)
258     .def ("__add__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::UNION, self, other); }))
259     .def ("__mul__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::SECTION, self, other); }))
260     .def ("__sub__", FunctionPointer
261           ([] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other)
262            { return make_shared<SPSolid> (SPSolid::SECTION, self,
263                                           make_shared<SPSolid> (SPSolid::SUB, other, nullptr)); }))
264 
265     .def ("bc", FunctionPointer([](shared_ptr<SPSolid> & self, int nr) -> shared_ptr<SPSolid>
266                                 { self->SetBC(nr); return self; }))
267     .def ("bc", FunctionPointer([](shared_ptr<SPSolid> & self, string name) -> shared_ptr<SPSolid>
268                                 { self->SetBCName(name); return self; }))
269     .def ("maxh", FunctionPointer([](shared_ptr<SPSolid> & self, double maxh) -> shared_ptr<SPSolid>
270                                 { self->SetMaxH(maxh); return self; }))
271     .def ("mat", FunctionPointer([](shared_ptr<SPSolid> & self, string mat) -> shared_ptr<SPSolid>
272                                  { self->SetMaterial(mat); return self; }))
273     .def ("mat", &SPSolid::GetMaterial)
274     .def("col", FunctionPointer([](shared_ptr<SPSolid> & self, py::list rgb) -> shared_ptr<SPSolid>
275                                 {
276                                   py::extract<double> red(rgb[0]);
277                                   py::extract<double> green(rgb[1]);
278                                   py::extract<double> blue(rgb[2]);
279                                   self->SetColor(red(),green(),blue());
280                                   return self;
281                                 }))
282     .def("transp", FunctionPointer([](shared_ptr<SPSolid> & self)->shared_ptr < SPSolid > { self->SetTransparent(); return self; }))
283     ;
284 
285   m.def ("Sphere", FunctionPointer([](Point<3> c, double r)
286                                      {
287                                        Sphere * sp = new Sphere (c, r);
288                                        Solid * sol = new Solid (sp);
289                                        return make_shared<SPSolid> (sol);
290                                      }));
291   m.def ("Ellipsoid", FunctionPointer([](Point<3> m, Vec<3> a, Vec<3> b, Vec<3> c)
292                                      {
293                                        Ellipsoid * ell = new Ellipsoid (m, a, b, c);
294                                        Solid * sol = new Solid (ell);
295                                        return make_shared<SPSolid> (sol);
296                                      }));
297   m.def ("Plane", FunctionPointer([](Point<3> p, Vec<3> n)
298                                     {
299                                       Plane * sp = new Plane (p,n);
300                                       Solid * sol = new Solid (sp);
301                                       return make_shared<SPSolid> (sol);
302                                     }));
303   m.def ("Cone", FunctionPointer([](Point<3> a, Point<3> b, double ra, double rb)
304                                        {
305                                          Cone * cyl = new Cone (a, b, ra, rb);
306                                          Solid * sol = new Solid (cyl);
307                                          return make_shared<SPSolid> (sol);
308                                        }));
309   m.def ("Cylinder", FunctionPointer([](Point<3> a, Point<3> b, double r)
310                                        {
311                                          Cylinder * cyl = new Cylinder (a, b, r);
312                                          Solid * sol = new Solid (cyl);
313                                          return make_shared<SPSolid> (sol);
314                                        }));
315   m.def ("OrthoBrick", FunctionPointer([](Point<3> p1, Point<3> p2)
316                                          {
317                                            OrthoBrick * brick = new OrthoBrick (p1,p2);
318                                            Solid * sol = new Solid (brick);
319                                            return make_shared<SPSolid> (sol);
320                                          }));
321   m.def ("Torus", FunctionPointer([](Point<3> c, Vec<3> n, double R, double r)
322                                          {
323                                            Torus * torus = new Torus (c,n,R,r);
324                                            Solid * sol = new Solid (torus);
325                                            return make_shared<SPSolid> (sol);
326                                          }));
327   m.def ("Revolution", [](Point<3> p1, Point<3> p2,
328                           shared_ptr<SplineGeometry<2>> spline)
329   {
330     Revolution * rev = new Revolution (p1, p2, spline);
331     Solid * sol = new Solid(rev);
332     return make_shared<SPSolid> (sol);
333   });
334   m.def ("Extrusion", [](shared_ptr<SplineGeometry<3>> path,
335                          shared_ptr<SplineGeometry<2>> profile,
336                          Vec<3> d)
337   {
338     Extrusion * extr = new Extrusion (path,profile,d);
339     Solid * sol = new Solid(extr);
340     return make_shared<SPSolid> (sol);
341   }, py::arg("path"), py::arg("profile"), py::arg("d"),
342      R"delimiter(A body of extrusion is defined by its profile
343 (which has to be a closed, clockwiseoriented 2D curve),
344  by a path (a 3D curve) and a vector d. It is constructed
345  as follows: Take a point p on the path and denote the
346  (unit-)tangent of the path in this point by t. If we cut
347  the body by the plane given by p and t as normal vector,
348  the cut is the profile. The profile is oriented by the
349  (local) y-direction `y:=d−(d·t)t` and the (local) x-direction
350  `x:=t \times y`.
351 The following points have to be noticed:
352  * If the path is not closed, then also the body is NOT closed.
353    In this case e.g. planes or orthobricks have to be used to
354    construct a closed body.
355  * The path has to be smooth, i.e. the tangents at the end- resp.
356    start-point of two consecutive spline or line patches have to
357    have the same directions.
358 )delimiter");
359   m.def("EllipticCone", [](const Point<3>& a, const Vec<3>& v, const Vec<3>& w,
360                             double h, double r)
361         {
362           auto ellcone = new EllipticCone(a,v,w,h,r);
363           auto sol = new Solid(ellcone);
364           return make_shared<SPSolid>(sol);
365         }, py::arg("a"), py::arg("vl"), py::arg("vs"), py::arg("h"), py::arg("r"),
366         R"raw_string(
367 An elliptic cone, given by the point 'a' at the base of the cone along the main axis,
368 the vectors v and w of the long and short axis of the ellipse, respectively,
369 the height of the cone, h, and ratio of base long axis length to top long axis length, r
370 
371 Note: The elliptic cone has to be truncated by planes similar to a cone or an elliptic cylinder.
372 When r =1, the truncated elliptic cone becomes an elliptic cylinder.
373 When r tends to zero, the truncated elliptic cone tends to a full elliptic cone.
374 However, when r = 0, the top part becomes a point(tip) and meshing fails!
375 )raw_string");
376 
377   m.def("Polyhedron", [](py::list points, py::list faces)
378   {
379     auto poly = new Polyhedra();
380     for(auto p : points)
381       poly->AddPoint(py::cast<Point<3>>(p));
382     int fnr = 0;
383     for(auto face : faces)
384       {
385         auto lface = py::cast<py::list>(face);
386         if(py::len(lface) == 3)
387           poly->AddFace(py::cast<int>(lface[0]),
388                         py::cast<int>(lface[1]),
389                         py::cast<int>(lface[2]),
390                         fnr++);
391         else if(py::len(lface) == 4)
392           {
393             poly->AddFace(py::cast<int>(lface[0]),
394                           py::cast<int>(lface[1]),
395                           py::cast<int>(lface[2]),
396                           fnr);
397             poly->AddFace(py::cast<int>(lface[0]),
398                           py::cast<int>(lface[2]),
399                           py::cast<int>(lface[3]),
400                           fnr++);
401           }
402       }
403     return make_shared<SPSolid>(new Solid(poly));
404   });
405 
406   m.def ("Or", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
407                                  {
408                                    return make_shared<SPSolid> (SPSolid::UNION, s1, s2);
409                                  }));
410   m.def ("And", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
411                                   {
412                                     return make_shared<SPSolid> (SPSolid::SECTION, s1, s2);
413                                   }));
414 
415 
416   py::class_<CSGeometry, NetgenGeometry, shared_ptr<CSGeometry>> (m, "CSGeometry")
417     .def(py::init<>())
418     .def(py::init([](const string& filename)
419                   {
420                     ifstream ist (filename);
421                     auto geo = make_shared<CSGeometry>();
422                     ParseCSG(ist, geo.get());
423                     geo->FindIdenticSurfaces(1e-8 * geo->MaxSize());
424                     return geo;
425                   }), py::arg("filename"))
426     .def(NGSPickle<CSGeometry>())
427     .def("Save", FunctionPointer([] (CSGeometry & self, string filename)
428                                  {
429                                    cout << "save geometry to file " << filename << endl;
430                                    self.Save (filename);
431                                  }))
432     .def("Add",
433          [] (CSGeometry & self, shared_ptr<SPSolid> solid, py::list bcmod, double maxh,
434              py::tuple col, bool transparent, int layer)
435           {
436             solid->AddSurfaces (self);
437             solid->GiveUpOwner();
438             int tlonr = self.SetTopLevelObject (solid->GetSolid());
439             self.GetTopLevelObject(tlonr) -> SetMaterial(solid->GetMaterial());
440             self.GetTopLevelObject(tlonr) -> SetRGB(solid->GetRed(),solid->GetGreen(),solid->GetBlue());
441             // self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent());
442             self.GetTopLevelObject(tlonr)->SetTransparent(transparent);
443             self.GetTopLevelObject(tlonr)->SetMaxH(maxh);
444             self.GetTopLevelObject(tlonr)->SetLayer(layer);
445 
446             // cout << "rgb = " << py::len(rgb) << endl;
447             if (py::len(col)==3)
448               self.GetTopLevelObject(tlonr) -> SetRGB(py::cast<double>(col[0]),
449                                                       py::cast<double>(col[1]),
450                                                       py::cast<double>(col[2]));
451 
452             // bcmod is list of tuples ( solid, bcnr )
453             for (int i = 0; i < py::len(bcmod); i++)
454               {
455                 py::tuple tup = py::extract<py::tuple> (bcmod[i]) ();
456                 auto mod_solid = py::extract<shared_ptr<SPSolid>> (tup[0]) ();
457                 int mod_nr = -1;
458                 string * bcname = nullptr;
459                 py::object val = tup[1];
460                 if (py::extract<int>(val).check()) mod_nr = py::extract<int> (val)();
461                 if (py::extract<string>(val).check()) bcname = new string ( py::extract<string> (val)());
462 
463                 NgArray<int> si;
464                 mod_solid -> GetSolid() -> GetSurfaceIndices (si);
465                 // cout << "change bc on surfaces: " << si << " to " << mod_nr << endl;
466 
467                 for (int j = 0; j < si.Size(); j++)
468                   {
469                     CSGeometry::BCModification bcm;
470                     bcm.bcname = bcname ? new string (*bcname) : nullptr;
471                     bcm.tlonr = tlonr;
472                     bcm.si = si[j];
473 		    bcm.bcnr = mod_nr;
474 		    self.bcmodifications.Append (bcm);
475                   }
476                 delete bcname;
477               }
478             return tlonr;
479           },
480          py::arg("solid"), py::arg("bcmod")=py::list(), py::arg("maxh")=1e99,
481          py::arg("col")=py::tuple(), py::arg("transparent")=false, py::arg("layer")=1
482          )
483 
484     .def("AddSurface", FunctionPointer
485          ([] (CSGeometry & self, shared_ptr<SPSolid> surface, shared_ptr<SPSolid> solid)
486           {
487             solid->AddSurfaces (self);
488             solid->GiveUpOwner();
489             Surface & surf = surface->GetSolid()->GetPrimitive()->GetSurface();
490             int tlonr = self.SetTopLevelObject (solid->GetSolid(), &surf);
491             // self.GetTopLevelObject(tlonr) -> SetMaterial(solid->GetMaterial());
492             self.GetTopLevelObject(tlonr) -> SetBCProp(surf.GetBCProperty());
493             self.GetTopLevelObject(tlonr) -> SetBCName(surf.GetBCName());
494 
495             self.GetTopLevelObject(tlonr) -> SetRGB(solid->GetRed(),solid->GetGreen(),solid->GetBlue());
496             self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent());
497           }),
498          py::arg("surface"), py::arg("solid")
499          )
500     .def("AddSplineSurface", FunctionPointer
501 	 ([] (CSGeometry & self, shared_ptr<SplineSurface> surf)
502 	  {
503 	    auto cuttings = surf->CreateCuttingSurfaces();
504 	    auto spsol = make_shared<SPSolid>(new Solid(surf.get()));
505 	    for(auto cut : (*cuttings)){
506 	      spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(cut.get())));
507 	    }
508 	    spsol->AddSurfaces(self);
509 	    int tlonr = self.SetTopLevelObject(spsol->GetSolid(), surf.get());
510 	    self.GetTopLevelObject(tlonr) -> SetBCProp(surf->GetBase()->GetBCProperty());
511 	    self.GetTopLevelObject(tlonr) -> SetBCName(surf->GetBase()->GetBCName());
512 	    self.GetTopLevelObject(tlonr) -> SetMaxH(surf->GetBase()->GetMaxH());
513             NgArray<Point<3>> non_midpoints;
514             for(auto spline : surf->GetSplines())
515               {
516                 non_midpoints.Append(spline->GetPoint(0));
517               }
518 	    for(auto p : non_midpoints)
519 		self.AddUserPoint(p);
520             self.AddSplineSurface(surf);
521 	  }),
522 	  py::arg("SplineSurface"))
523     .def("SingularFace", [] (CSGeometry & self, shared_ptr<SPSolid> sol, shared_ptr<SPSolid> surfaces, double factor)
524          {
525            int tlonum = -1;
526            for (int i = 0; i < self.GetNTopLevelObjects(); i++)
527              if (self.GetTopLevelObject(i)->GetSolid() == sol->GetSolid())
528                tlonum = i;
529            if (tlonum == -1) throw NgException("not a top-level-object");
530            if (!surfaces) surfaces = sol;
531            auto singface = new SingularFace(tlonum+1, surfaces->GetSolid(), factor);
532            self.singfaces.Append(singface);
533          }, py::arg("solid"), py::arg("surfaces")=nullptr, py::arg("factor")=0.25)
534     .def("SingularEdge", [] (CSGeometry & self, shared_ptr<SPSolid> s1,shared_ptr<SPSolid> s2, double factor)
535          {
536            auto singedge = new SingularEdge(1, -1, self, s1->GetSolid(), s2->GetSolid(), factor);
537            self.singedges.Append (singedge);
538          })
539     .def("SingularPoint", [] (CSGeometry & self, shared_ptr<SPSolid> s1,shared_ptr<SPSolid> s2,
540                              shared_ptr<SPSolid> s3, double factor)
541          {
542            auto singpoint = new SingularPoint(1, s1->GetSolid(), s2->GetSolid(), s3->GetSolid(), factor);
543            self.singpoints.Append (singpoint);
544          })
545     .def("CloseSurfaces", FunctionPointer
546          ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, py::list aslices )
547           {
548             NgArray<int> si1, si2;
549             s1->GetSolid()->GetSurfaceIndices (si1);
550             s2->GetSolid()->GetSurfaceIndices (si2);
551             Flags flags;
552 
553             try
554             {
555                 int n = py::len(aslices);
556                 Array<double> slices(n);
557                 for(int i=0; i<n; i++)
558                 {
559                     slices[i]= py::extract<double>(aslices[i])();
560                 }
561                 flags.SetFlag("slices", slices);
562             }
563             catch( py::error_already_set const & ) {
564                 cout << "caught python error:" << endl;
565                 PyErr_Print();
566             }
567 
568             const TopLevelObject * domain = nullptr;
569             self.AddIdentification
570               (new CloseSurfaceIdentification
571                (self.GetNIdentifications()+1, self,
572                 self.GetSurface (si1[0]), self.GetSurface (si2[0]),
573                 domain,
574                 flags));
575           }),
576          py::arg("solid1"), py::arg("solid2"), py::arg("slices")
577          )
578     .def("CloseSurfaces", FunctionPointer
579          ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2,
580               int reflevels, shared_ptr<SPSolid> domain_solid)
581           {
582             NgArray<int> si1, si2;
583             s1->GetSolid()->GetSurfaceIndices (si1);
584             s2->GetSolid()->GetSurfaceIndices (si2);
585             cout << "surface ids1 = " << si1 << endl;
586             cout << "surface ids2 = " << si2 << endl;
587 
588             Flags flags;
589             const TopLevelObject * domain = nullptr;
590             if (domain_solid)
591               domain = self.GetTopLevelObject(domain_solid->GetSolid());
592 
593             self.AddIdentification
594               (new CloseSurfaceIdentification
595                (self.GetNIdentifications()+1, self,
596                 self.GetSurface (si1[0]), self.GetSurface (si2[0]),
597                 domain,
598                 flags));
599           }),
600          py::arg("solid1"), py::arg("solid2"), py::arg("reflevels")=2, py::arg("domain")=nullptr
601          )
602 
603     .def("PeriodicSurfaces", FunctionPointer
604          ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2,
605               Transformation<3> trafo)
606           {
607             NgArray<int> si1, si2;
608             s1->GetSolid()->GetSurfaceIndices (si1);
609             s2->GetSolid()->GetSurfaceIndices (si2);
610             cout << "identify surfaces " << si1[0] << " and " << si2[0] << endl;
611             self.AddIdentification
612               (new PeriodicIdentification
613                (self.GetNIdentifications()+1, self,
614                 self.GetSurface (si1[0]), self.GetSurface (si2[0]),
615                 trafo));
616           }),
617          py::arg("solid1"), py::arg("solid2"),
618          py::arg("trafo")=Transformation<3>(Vec<3>(0,0,0))
619          )
620     .def("NameEdge", [] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, string name)
621          {
622            Array<Surface*> surfs1, surfs2;
623            s1->GetSolid()->ForEachSurface( [&surfs1] (Surface * s, bool inv) { surfs1.Append(s); });
624            s2->GetSolid()->ForEachSurface( [&surfs2] (Surface * s, bool inv) { surfs2.Append(s); });
625            for (auto s1 : surfs1)
626              for (auto s2 : surfs2)
627                self.named_edges[tuple(s1,s2)] = name;
628          })
629 
630     .def("AddPoint", [] (CSGeometry & self, Point<3> p, variant<int,string> index) -> CSGeometry&
631          {
632            if (auto pint = std::get_if<int> (&index))
633              self.AddUserPoint(CSGeometry::UserPoint(p, *pint));
634            if (auto pstr = std::get_if<string> (&index))
635              self.AddUserPoint(CSGeometry::UserPoint(p, *pstr));
636            return self;
637          })
638 
639     .def("GetTransparent", FunctionPointer
640          ([] (CSGeometry & self, int tlonr)
641           {
642             return self.GetTopLevelObject(tlonr)->GetTransparent();
643           }),
644          py::arg("tlonr")
645          )
646     .def("SetTransparent", FunctionPointer
647          ([] (CSGeometry & self, int tlonr, bool transparent)
648           {
649             self.GetTopLevelObject(tlonr)->SetTransparent(transparent);
650           }),
651          py::arg("tlonr"), py::arg("transparent")
652          )
653 
654     .def("GetVisible", FunctionPointer
655          ([] (CSGeometry & self, int tlonr)
656           {
657             return self.GetTopLevelObject(tlonr)->GetVisible();
658           }),
659          py::arg("tlonr")
660          )
661     .def("SetVisible", FunctionPointer
662          ([] (CSGeometry & self, int tlonr, bool visible)
663           {
664             self.GetTopLevelObject(tlonr)->SetVisible(visible);
665           }),
666          py::arg("tlonr"), py::arg("visible")
667          )
668     .def("SetBoundingBox", FunctionPointer
669          ([] (CSGeometry & self, Point<3> pmin, Point<3> pmax)
670           {
671             self.SetBoundingBox(Box<3> (pmin, pmax));
672           }),
673          py::arg("pmin"), py::arg("pmax")
674          )
675     .def("Draw", FunctionPointer
676          ([] (shared_ptr<CSGeometry> self)
677           {
678              self->FindIdenticSurfaces(1e-6);
679              self->CalcTriangleApproximation(0.01, 20);
680              ng_geometry = self;
681           })
682          )
683     .def("GetSolids", [](CSGeometry& self)
684                       {
685                         py::list lst;
686                         for(auto i : Range(self.GetSolids().Size()))
687                           lst.append(make_shared<SPSolid>(self.GetSolids()[i], 1234));
688                         return lst;
689                       })
690     .def_property_readonly ("ntlo", &CSGeometry::GetNTopLevelObjects)
691     .def("_visualizationData", [](shared_ptr<CSGeometry> csg_geo)
692          {
693            std::vector<float> vertices;
694            std::vector<int> trigs;
695            std::vector<float> normals;
696            std::vector<float> min = {std::numeric_limits<float>::max(),
697                                      std::numeric_limits<float>::max(),
698                                      std::numeric_limits<float>::max()};
699            std::vector<float> max = {std::numeric_limits<float>::lowest(),
700                                      std::numeric_limits<float>::lowest(),
701                                      std::numeric_limits<float>::lowest()};
702            std::vector<string> surfnames;
703            for (int i = 0; i < csg_geo->GetNSurf(); i++)
704              {
705                auto surf = csg_geo->GetSurface(i);
706                surfnames.push_back(surf->GetBCName());
707              }
708            csg_geo->FindIdenticSurfaces(1e-6);
709            csg_geo->CalcTriangleApproximation(0.01,100);
710            auto nto = csg_geo->GetNTopLevelObjects();
711            size_t np = 0;
712            size_t ntrig = 0;
713            for (int i = 0; i < nto; i++){
714              np += csg_geo->GetTriApprox(i)->GetNP();
715              ntrig += csg_geo->GetTriApprox(i)->GetNT();
716            }
717            vertices.reserve(np*3);
718            trigs.reserve(ntrig*4);
719            normals.reserve(np*3);
720            int offset_points = 0;
721            for (int i = 0; i < nto; i++)
722              {
723                auto triapprox = csg_geo->GetTriApprox(i);
724                for (int j = 0; j < triapprox->GetNP(); j++)
725                  for(int k = 0; k < 3; k++) {
726                    float val = triapprox->GetPoint(j)[k];
727                    vertices.push_back(val);
728                    min[k] = min2(min[k], val);
729                    max[k] = max2(max[k],val);
730                    normals.push_back(triapprox->GetNormal(j)[k]);
731                  }
732                for (int j = 0; j < triapprox->GetNT(); j++)
733                  {
734                    for(int k = 0; k < 3; k++)
735                      trigs.push_back(triapprox->GetTriangle(j)[k]+offset_points);
736                    trigs.push_back(triapprox->GetTriangle(j).SurfaceIndex());
737                  }
738                offset_points += triapprox->GetNP();
739              }
740            py::gil_scoped_acquire ac;
741            py::dict res;
742            py::list snames;
743            for(auto name : surfnames)
744              snames.append(py::cast(name));
745            res["vertices"] = MoveToNumpy(vertices);
746            res["triangles"] = MoveToNumpy(trigs);
747            res["normals"] = MoveToNumpy(normals);
748            res["surfnames"] = snames;
749            res["min"] = MoveToNumpy(min);
750            res["max"] = MoveToNumpy(max);
751            return res;
752          }, py::call_guard<py::gil_scoped_release>())
753   .def("GenerateMesh", [](shared_ptr<CSGeometry> geo,
754                           MeshingParameters* pars, py::kwargs kwargs)
755            {
756              MeshingParameters mp;
757              if(pars) mp = *pars;
758              {
759                py::gil_scoped_acquire aq;
760                CreateMPfromKwargs(mp, kwargs);
761              }
762              auto mesh = make_shared<Mesh>();
763              SetGlobalMesh (mesh);
764              mesh->SetGeometry(geo);
765 	     ng_geometry = geo;
766              geo->FindIdenticSurfaces(1e-8 * geo->MaxSize());
767              auto result = geo->GenerateMesh (mesh, mp);
768              if(result != 0)
769                throw Exception("Meshing failed!");
770              return mesh;
771            }, py::arg("mp") = nullptr,
772        meshingparameter_description.c_str(),
773     py::call_guard<py::gil_scoped_release>())
774     ;
775 
776   m.def("Save", FunctionPointer
777           ([](const Mesh & self, const string & filename, const CSGeometry & geom)
778            {
779              ostream * outfile;
780              if (filename.substr (filename.length()-3, 3) == ".gz")
781                outfile = new ogzstream (filename.c_str());
782              else
783                outfile = new ofstream (filename.c_str());
784 
785              self.Save (*outfile);
786              *outfile << endl << endl << "endmesh" << endl << endl;
787              geom.SaveToMeshFile (*outfile);
788              delete outfile;
789            }),py::call_guard<py::gil_scoped_release>())
790     ;
791 
792 
793 
794   m.def("ZRefinement", FunctionPointer
795           ([](Mesh & mesh, CSGeometry & geom)
796           {
797             ZRefinementOptions opt;
798             opt.minref = 5;
799             ZRefinement (mesh, &geom, opt);
800           }),py::call_guard<py::gil_scoped_release>())
801     ;
802 }
803 
PYBIND11_MODULE(libcsg,m)804 PYBIND11_MODULE(libcsg, m) {
805   ExportCSG(m);
806 }
807 #endif
808 
809