diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index cc3e4464..e1ada634 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -176,9 +176,8 @@ namespace netgen solids.DeleteAll (); - for (int i = 0; i < splinecurves2d.Size(); i++) - delete splinecurves2d[i]; splinecurves2d.DeleteAll(); + splinecurves3d.DeleteAll(); /* for (int i = 0; i < surfaces.Size(); i++) @@ -712,24 +711,24 @@ namespace netgen - void CSGeometry :: SetSplineCurve (const char * name, SplineGeometry<2> * spl) + void CSGeometry :: SetSplineCurve (const char * name, shared_ptr> spl) { splinecurves2d.Set(name,spl); } - void CSGeometry :: SetSplineCurve (const char * name, SplineGeometry<3> * spl) + void CSGeometry :: SetSplineCurve (const char * name, shared_ptr> spl) { splinecurves3d.Set(name,spl); } - const SplineGeometry<2> * CSGeometry :: GetSplineCurve2d (const string & name) const + shared_ptr> CSGeometry :: GetSplineCurve2d (const string & name) const { if (splinecurves2d.Used(name)) return splinecurves2d[name]; else return NULL; } - const SplineGeometry<3> * CSGeometry :: GetSplineCurve3d (const string & name) const + shared_ptr> CSGeometry :: GetSplineCurve3d (const string & name) const { if (splinecurves3d.Used(name)) return splinecurves3d[name]; diff --git a/libsrc/csg/csgeom.hpp b/libsrc/csg/csgeom.hpp index 35453b30..17d0a3d9 100644 --- a/libsrc/csg/csgeom.hpp +++ b/libsrc/csg/csgeom.hpp @@ -115,9 +115,9 @@ namespace netgen SymbolTable solids; /// all 2d splinecurves - SymbolTable< SplineGeometry<2>* > splinecurves2d; + SymbolTable>> splinecurves2d; /// all 3d splinecurves - SymbolTable< SplineGeometry<3>* > splinecurves3d; + SymbolTable>> splinecurves3d; /// all top level objects: solids and surfaces NgArray toplevelobjects; @@ -232,10 +232,10 @@ namespace netgen const SymbolTable & GetSolids () const { return solids; } - void SetSplineCurve (const char * name, SplineGeometry<2> * spl); - void SetSplineCurve (const char * name, SplineGeometry<3> * spl); - const SplineGeometry<2> * GetSplineCurve2d (const string & name) const; - const SplineGeometry<3> * GetSplineCurve3d (const string & name) const; + void SetSplineCurve (const char * name, shared_ptr> spl); + void SetSplineCurve (const char * name, shared_ptr> spl); + shared_ptr> GetSplineCurve2d (const string & name) const; + shared_ptr> GetSplineCurve3d (const string & name) const; void DoArchive(Archive& archive) override; diff --git a/libsrc/csg/csgparser.cpp b/libsrc/csg/csgparser.cpp index 4e3c1786..05c1bdcd 100644 --- a/libsrc/csg/csgparser.cpp +++ b/libsrc/csg/csgparser.cpp @@ -511,8 +511,8 @@ namespace netgen break; } - Primitive * nprim = new Extrusion(*(geom->GetSplineCurve3d(epath)), - *(geom->GetSplineCurve2d(profile)), + Primitive * nprim = new Extrusion(geom->GetSplineCurve3d(epath), + geom->GetSplineCurve2d(profile), z_dir); geom->AddSurfaces (nprim); return new Solid(nprim); @@ -1186,7 +1186,7 @@ namespace netgen ParseChar (scan, '='); ParseChar (scan, '('); - SplineGeometry<2> * newspline = new SplineGeometry<2>; + auto newspline = make_shared>(); // newspline->CSGLoad(scan); LoadSpline (*newspline, scan); @@ -1212,7 +1212,7 @@ namespace netgen ParseChar (scan, '='); ParseChar (scan, '('); - SplineGeometry<3> * newspline = new SplineGeometry<3>; + auto newspline = make_shared>(); // newspline->CSGLoad(scan); LoadSpline (*newspline, scan); diff --git a/libsrc/csg/extrusion.cpp b/libsrc/csg/extrusion.cpp index c5518ab8..82b2e651 100644 --- a/libsrc/csg/extrusion.cpp +++ b/libsrc/csg/extrusion.cpp @@ -41,6 +41,18 @@ namespace netgen loc_z_dir[i] = glob_z_direction; } } + + double cum_angle = 0.; + for(auto i : Range(path->GetSplines())) + { + const auto& sp = path->GetSpline(i); + auto t1 = sp.GetTangent(0.); + t1.Normalize(); + auto t2 = sp.GetTangent(1.); + t2.Normalize(); + cum_angle += acos(t1 * t2); + angles.Append(cum_angle); + } profile->GetCoeff(profile_spline_coeff); latest_point3d = -1.111e30; @@ -656,20 +668,35 @@ namespace netgen dez /= lenz; dez -= (dez * ez) * ez; } - - Extrusion :: Extrusion(const SplineGeometry<3> & path_in, - const SplineGeometry<2> & profile_in, + void ExtrusionFace :: DefineTangentialPlane(const Point<3>& ap1, + const Point<3>& ap2) + { + Surface::DefineTangentialPlane(ap1, ap2); + tangential_plane_seg = latest_seg; + } + + void ExtrusionFace :: ToPlane(const Point<3>& p3d, Point<2>& p2d, + double h, int& zone) const + { + Surface::ToPlane(p3d, p2d, h, zone); + double angle = angles[tangential_plane_seg] - angles[latest_seg]; + if(fabs(angle) > 3.14/2.) + zone = -1; + } + + Extrusion :: Extrusion(shared_ptr> path_in, + shared_ptr> profile_in, const Vec<3> & z_dir) : - path(&path_in), profile(&profile_in), z_direction(z_dir) + path(path_in), profile(profile_in), z_direction(z_dir) { surfaceactive.SetSize(0); surfaceids.SetSize(0); for(int j=0; jGetNSplines(); j++) { - ExtrusionFace * face = new ExtrusionFace(&((*profile).GetSpline(j)), - path, + ExtrusionFace * face = new ExtrusionFace(&(profile->GetSpline(j)), + path.get(), z_direction); faces.Append(face); surfaceactive.Append(true); @@ -856,7 +883,7 @@ namespace netgen return retval; if(latestfacenum >= 0) - return faces[latestfacenum]->VecInFace(p,v2,0); + return faces[latestfacenum]->VecInFace(p,v2,eps); else return VecInSolid(p,v2,eps); } diff --git a/libsrc/csg/extrusion.hpp b/libsrc/csg/extrusion.hpp index e9680eff..af951c51 100644 --- a/libsrc/csg/extrusion.hpp +++ b/libsrc/csg/extrusion.hpp @@ -12,8 +12,10 @@ namespace netgen const SplineSeg<2> * profile; const SplineGeometry<3> * path; Vec<3> glob_z_direction; + Array angles; bool deletable; + int tangential_plane_seg; NgArray< const SplineSeg3<3> * > spline3_path; NgArray< const LineSeg<3> * > line_path; @@ -114,6 +116,11 @@ namespace netgen Vec<3> & ex, Vec<3> & ey, Vec<3> & ez, Vec<3> & dex, Vec<3> & dey, Vec<3> & dez) const; + void DefineTangentialPlane(const Point<3>& ap1, + const Point<3>& ap2) override; + void ToPlane(const Point<3>& p3d, Point<2>& p2d, + double h, int& zone) const override; + }; @@ -121,8 +128,8 @@ namespace netgen class Extrusion : public Primitive { private: - const SplineGeometry<3>* path; - const SplineGeometry<2>* profile; // closed, clockwise oriented curve + shared_ptr> path; + shared_ptr> profile; // closed, clockwise oriented curve Vec<3> z_direction; @@ -131,8 +138,8 @@ namespace netgen mutable int latestfacenum; public: - Extrusion(const SplineGeometry<3> & path_in, - const SplineGeometry<2> & profile_in, + Extrusion(shared_ptr> path_in, + shared_ptr> profile_in, const Vec<3> & z_dir); // default constructor for archive Extrusion() {} diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 0954b2d6..a650a5ae 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -170,7 +170,8 @@ namespace netgen DLL_HEADER void ExportCSG(py::module &m) { - py::class_> (m, "SplineCurve2d") + py::class_, shared_ptr>> + (m, "SplineCurve2d") .def(py::init<>()) .def ("AddPoint", FunctionPointer ([] (SplineGeometry<2> & self, double x, double y) @@ -329,14 +330,31 @@ DLL_HEADER void ExportCSG(py::module &m) Solid * sol = new Solid(rev); return make_shared (sol); })); - m.def ("Extrusion", FunctionPointer([](const SplineGeometry<3> & path, - const SplineGeometry<2> & profile, - Vec<3> n) - { - Extrusion * extr = new Extrusion (path,profile,n); - Solid * sol = new Solid(extr); - return make_shared (sol); - })); + m.def ("Extrusion", [](shared_ptr> path, + shared_ptr> profile, + Vec<3> d) + { + Extrusion * extr = new Extrusion (path,profile,d); + Solid * sol = new Solid(extr); + return make_shared (sol); + }, py::arg("path"), py::arg("profile"), py::arg("d"), + R"delimiter(A body of extrusion is defined by its profile +(which has to be a closed, clockwiseoriented 2D curve), + by a path (a 3D curve) and a vector d. It is constructed + as follows: Take a point p on the path and denote the + (unit-)tangent of the path in this point by t. If we cut + the body by the plane given by p and t as normal vector, + the cut is the profile. The profile is oriented by the + (local) y-direction `y:=d−(d·t)t` and the (local) x-direction + `x:=t \times y`. +The following points have to be noticed: + * If the path is not closed, then also the body is NOT closed. + In this case e.g. planes or orthobricks have to be used to + construct a closed body. + * The path has to be smooth, i.e. the tangents at the end- resp. + start-point of two consecutive spline or line patches have to + have the same directions. +)delimiter"); m.def("EllipticCone", [](const Point<3>& a, const Vec<3>& v, const Vec<3>& w, double h, double r) { diff --git a/tests/pytest/results.json b/tests/pytest/results.json index 0ece1b09..45b04c68 100644 --- a/tests/pytest/results.json +++ b/tests/pytest/results.json @@ -1272,6 +1272,98 @@ "total_badness": 65897.969985 } ], + "extrusion.geo": [ + { + "angles_tet": [ + 6.6841, + 171.53 + ], + "angles_trig": [ + 11.293, + 152.07 + ], + "ne1d": 172, + "ne2d": 286, + "ne3d": 241, + "quality_histogram": "[0, 0, 7, 56, 39, 22, 0, 0, 0, 0, 2, 0, 5, 18, 18, 29, 21, 10, 9, 5]", + "total_badness": 775.80779693 + }, + { + "angles_tet": [ + 16.097, + 160.17 + ], + "angles_trig": [ + 15.327, + 149.09 + ], + "ne1d": 104, + "ne2d": 152, + "ne3d": 124, + "quality_histogram": "[0, 0, 0, 0, 10, 19, 39, 26, 9, 3, 1, 1, 5, 1, 3, 2, 3, 2, 0, 0]", + "total_badness": 353.53219387 + }, + { + "angles_tet": [ + 11.505, + 165.94 + ], + "angles_trig": [ + 15.054, + 147.98 + ], + "ne1d": 134, + "ne2d": 196, + "ne3d": 167, + "quality_histogram": "[0, 0, 0, 1, 3, 35, 33, 8, 5, 21, 11, 9, 10, 11, 7, 3, 5, 3, 1, 1]", + "total_badness": 417.63980201 + }, + { + "angles_tet": [ + 6.6841, + 171.53 + ], + "angles_trig": [ + 11.293, + 152.07 + ], + "ne1d": 172, + "ne2d": 286, + "ne3d": 241, + "quality_histogram": "[0, 0, 7, 56, 39, 22, 0, 0, 0, 0, 2, 0, 5, 18, 18, 29, 21, 10, 9, 5]", + "total_badness": 775.80779693 + }, + { + "angles_tet": [ + 17.691, + 140.88 + ], + "angles_trig": [ + 18.812, + 116.06 + ], + "ne1d": 276, + "ne2d": 570, + "ne3d": 646, + "quality_histogram": "[0, 0, 0, 0, 0, 0, 4, 17, 30, 51, 51, 55, 81, 69, 74, 73, 72, 52, 12, 5]", + "total_badness": 1020.0235117 + }, + { + "angles_tet": [ + 14.402, + 155.5 + ], + "angles_trig": [ + 24.071, + 119.03 + ], + "ne1d": 442, + "ne2d": 1220, + "ne3d": 2802, + "quality_histogram": "[0, 0, 0, 0, 0, 2, 5, 2, 4, 12, 25, 55, 147, 298, 311, 503, 457, 536, 342, 103]", + "total_badness": 3603.431162 + } + ], "fichera.geo": [ { "angles_tet": [ diff --git a/tests/pytest/test_tutorials.py b/tests/pytest/test_tutorials.py index 90bccc3f..236749da 100644 --- a/tests/pytest/test_tutorials.py +++ b/tests/pytest/test_tutorials.py @@ -63,12 +63,12 @@ def getMeshingparameters(filename): standard = [MeshingParameters()] + [MeshingParameters(ms) for ms in (meshsize.very_coarse, meshsize.coarse, meshsize.moderate, meshsize.fine, meshsize.very_fine)] if filename == "shell.geo": return [] # do not test this example cause it needs so long... - if filename == "extrusion.geo": - return [] # this segfaults right now if filename == "manyholes2.geo": return [standard[1]] # this gets too big for finer meshsizes if filename in ("manyholes.geo", "frame.step"): return standard[:3] # this gets too big for finer meshsizes + if filename == "extrusion.geo": + return standard[:-1] if filename == "screw.step": return standard[3:] # coarser meshes don't work here if filename == "cylsphere.geo": @@ -142,7 +142,11 @@ def generateResultFile(): continue meshdata = [] for mp in mps: - mesh = generateMesh(_file, mp) + try: + mesh = generateMesh(_file, mp) + except Exception as e: + print("Meshingparameters: ", mp) + raise e meshdata.append( getData(mesh, mp) ) data[_file] = meshdata print("needed", time.time() - start, "seconds")