Merge branch 'fix_extrusion' into 'master'

Fix extrusion

See merge request jschoeberl/netgen!379

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<SplineGeometry<2>> spl)
   {
     splinecurves2d.Set(name,spl);
   }
-  void CSGeometry :: SetSplineCurve (const char * name, SplineGeometry<3> * spl)
+  void CSGeometry :: SetSplineCurve (const char * name, shared_ptr<SplineGeometry<3>> spl)
   {
     splinecurves3d.Set(name,spl);
   }
 
 
-  const SplineGeometry<2> * CSGeometry :: GetSplineCurve2d (const string & name) const
+  shared_ptr<SplineGeometry<2>> 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<SplineGeometry<3>> CSGeometry :: GetSplineCurve3d (const string & name) const
   {
     if (splinecurves3d.Used(name))
       return splinecurves3d[name];

libsrc/csg/csgeom.hpp

@@ -115,9 +115,9 @@ namespace netgen
     SymbolTable<Solid*> solids;
 
     /// all 2d splinecurves
-    SymbolTable< SplineGeometry<2>* > splinecurves2d;
+    SymbolTable<shared_ptr<SplineGeometry<2>>> splinecurves2d;
     /// all 3d splinecurves
-    SymbolTable< SplineGeometry<3>* > splinecurves3d;
+    SymbolTable<shared_ptr<SplineGeometry<3>>> splinecurves3d;
 
     /// all top level objects: solids and surfaces
     NgArray<TopLevelObject*> toplevelobjects;
@@ -232,10 +232,10 @@ namespace netgen
     const SymbolTable<Solid*> & GetSolids () const { return solids; }
 
 
-    void SetSplineCurve (const char * name, SplineGeometry<2> * spl);
+    void SetSplineCurve (const char * name, shared_ptr<SplineGeometry<2>> spl);
-    void SetSplineCurve (const char * name, SplineGeometry<3> * spl);
+    void SetSplineCurve (const char * name, shared_ptr<SplineGeometry<3>> spl);
-    const SplineGeometry<2> * GetSplineCurve2d (const string & name) const;
+    shared_ptr<SplineGeometry<2>> GetSplineCurve2d (const string & name) const;
-    const SplineGeometry<3> * GetSplineCurve3d (const string & name) const;
+    shared_ptr<SplineGeometry<3>> GetSplineCurve3d (const string & name) const;
 
     void DoArchive(Archive& archive) override;
     

libsrc/csg/csgparser.cpp

@@ -511,8 +511,8 @@ namespace netgen
 		  break;
 		}
 	      
-	      Primitive * nprim = new Extrusion(*(geom->GetSplineCurve3d(epath)),
+	      Primitive * nprim = new Extrusion(geom->GetSplineCurve3d(epath),
-						*(geom->GetSplineCurve2d(profile)),
+						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<SplineGeometry<2>>();
 		// 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<SplineGeometry<3>>();
 		// newspline->CSGLoad(scan);
 		LoadSpline (*newspline, scan);
 

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,
+  void ExtrusionFace :: DefineTangentialPlane(const Point<3>& ap1,
-			 const SplineGeometry<2> & profile_in,
+                                              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<SplineGeometry<3>> path_in,
+			 shared_ptr<SplineGeometry<2>> 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; j<profile->GetNSplines(); j++)
       {
-	ExtrusionFace * face = new ExtrusionFace(&((*profile).GetSpline(j)),
+	ExtrusionFace * face = new ExtrusionFace(&(profile->GetSpline(j)),
-						 path,
+						 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);
   }

libsrc/csg/extrusion.hpp

@@ -12,8 +12,10 @@ namespace netgen
     const SplineSeg<2> * profile;
     const SplineGeometry<3> * path;
     Vec<3> glob_z_direction;
+    Array<double> 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;
+    shared_ptr<SplineGeometry<3>> path;
-    const SplineGeometry<2>* profile; // closed, clockwise oriented curve
+    shared_ptr<SplineGeometry<2>> 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,
+    Extrusion(shared_ptr<SplineGeometry<3>> path_in,
-	      const SplineGeometry<2> & profile_in,
+	      shared_ptr<SplineGeometry<2>> profile_in,
 	      const Vec<3> & z_dir);
     // default constructor for archive
     Extrusion() {}

libsrc/csg/python_csg.cpp

@@ -170,7 +170,8 @@ namespace netgen
 
 DLL_HEADER void ExportCSG(py::module &m)
 {
-  py::class_<SplineGeometry<2>> (m, "SplineCurve2d")
+  py::class_<SplineGeometry<2>, shared_ptr<SplineGeometry<2>>>
+    (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<SPSolid> (sol);
                                          }));
-  m.def ("Extrusion", FunctionPointer([](const SplineGeometry<3> & path,
+  m.def ("Extrusion", [](shared_ptr<SplineGeometry<3>> path,
-					 const SplineGeometry<2> & profile,
+                         shared_ptr<SplineGeometry<2>> profile,
-					 Vec<3> n)
+                         Vec<3> d)
-                                         {
+  {
-                                           Extrusion * extr = new Extrusion (path,profile,n);
+    Extrusion * extr = new Extrusion (path,profile,d);
-                                           Solid * sol = new Solid(extr);
+    Solid * sol = new Solid(extr);
-                                           return make_shared<SPSolid> (sol);
+    return make_shared<SPSolid> (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)
         {

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": [

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")