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