diff --git a/libsrc/csg/algprim.cpp b/libsrc/csg/algprim.cpp index 4ab7ee9a..5e2b7fd3 100644 --- a/libsrc/csg/algprim.cpp +++ b/libsrc/csg/algprim.cpp @@ -1493,8 +1493,215 @@ namespace netgen +/// Elliptic Cone +/// Josephat Kalezhi (kalezhi@cbu.ac.zm) +/// February 21st, 2018 +/// + +EllipticCone :: EllipticCone (const Point<3> & aa, const Vec<3> & avl, + const Vec<3> & avs, double ah, double avlr) +{ + a = aa; + h = ah; + vlr = avlr; + + if (avl.Length2() >= avs.Length2()) + { + vl = avl; + vs = avs; + } + else + { + vl = avs; + vs = avl; + } + CalcData(); + // Print (cout); +} + + +Primitive * EllipticCone :: CreateDefault () + { + return new EllipticCone (Point<3> (0,0,0), Vec<3> (1,0,0), Vec<3> (0,1,0), 1, 0.5); + } + + + void EllipticCone :: GetPrimitiveData (const char *& classname, Array & coeffs) const + { + classname = "ellipticcone"; + coeffs.SetSize (15); + coeffs.Elem(1) = a(0); + coeffs.Elem(2) = a(1); + coeffs.Elem(3) = a(2); + coeffs.Elem(4) = vl(0); + coeffs.Elem(5) = vl(1); + coeffs.Elem(6) = vl(2); + coeffs.Elem(7) = vs(0); + coeffs.Elem(8) = vs(1); + coeffs.Elem(9) = vs(2); + coeffs.Elem(10) = h; + coeffs.Elem(11) = vlr; + + } + + + void EllipticCone :: SetPrimitiveData (Array & coeffs) + { + + a(0) = coeffs.Elem(1); + a(1) = coeffs.Elem(2); + a(2) = coeffs.Elem(3); + vl(0) = coeffs.Elem(4); + vl(1) = coeffs.Elem(5); + vl(2) = coeffs.Elem(6); + vs(0) = coeffs.Elem(7); + vs(1) = coeffs.Elem(8); + vs(2) = coeffs.Elem(9); + h = coeffs.Elem(10); + vlr = coeffs.Elem(11); + CalcData(); + } + + + +void EllipticCone :: CalcData () +{ + Vec<3> nh = Cross(vl, vs); + nh.Normalize(); + + double lvl = vl.Length(); + double lvs = vs.Length(); + + Vec<3> t1vec = lvl*(vlr -1)*(1/h)*nh; + + Vec<3> va (a); + double t1 = lvl*(1 - (vlr -1)*(1/h)*(va*nh)); + + Vec<3> nvl = (1.0/lvl)*vl; + Vec<3> nvs = (1.0/lvs)*vs; + double ellipt2 = sqr(lvl/lvs); + + cxx = nvl(0)*nvl(0) + ellipt2*nvs(0)*nvs(0) - t1vec(0)*t1vec(0); + cyy = nvl(1)*nvl(1) + ellipt2*nvs(1)*nvs(1) - t1vec(1)*t1vec(1); + czz = nvl(2)*nvl(2) + ellipt2*nvs(2)*nvs(2) - t1vec(2)*t1vec(2); + + cxy = 2*(nvl(0)*nvl(1) + ellipt2*nvs(0)*nvs(1) - t1vec(0)*t1vec(1)); + cxz = 2*(nvl(0)*nvl(2) + ellipt2*nvs(0)*nvs(2) - t1vec(0)*t1vec(2)); + cyz = 2*(nvl(1)*nvl(2) + ellipt2*nvs(1)*nvs(2) - t1vec(1)*t1vec(2)); + + Vec<3> v = -2*((va*nvl)*nvl + ellipt2*(va*nvs)*nvs + t1*t1vec); + cx = v(0); + cy = v(1); + cz = v(2); + + c1 = pow(va*nvl,2) + ellipt2*pow(va*nvs,2) - t1*t1; + + double lvltop = vlr*lvl; + double minlvl = (lvl < lvltop) ? lvl : lvltop; + double maxlvl = max2( lvl,lvltop); + cxx /= maxlvl; cyy /= maxlvl; czz /= maxlvl; + cxy /= maxlvl; cxz /= maxlvl; cyz /= maxlvl; + cx /= maxlvl; cy /= maxlvl; cz /= maxlvl; + c1 /= maxlvl; +} + + +INSOLID_TYPE EllipticCone :: BoxInSolid (const BoxSphere<3> & box) const +{ + double rp, dist; + + Vec<3> cv( box.Center()); + Vec<3> nh = Cross(vl, vs); + nh.Normalize(); + + double lvl = vl.Length(); + Vec<3> t1vec = lvl*(vlr -1)*(1/h)*nh; + Vec<3> va (a); + double t1 = lvl*(1 - (vlr -1)*(1/h)*(va*nh)); + rp = cv*t1vec + t1; + double lvltop = vlr*lvl; + double maxlvl = max2( lvl,lvltop); + + dist = sqrt( CalcFunctionValue(box.Center())*maxlvl + rp*rp) - rp; + + if (dist - box.Diam() > 0) return IS_OUTSIDE; + if (dist + box.Diam() < 0) return IS_INSIDE; + return DOES_INTERSECT; +} + +double EllipticCone :: HesseNorm () const +{ + return 1.0/min(vs.Length2 (),vl.Length2()); +} + +double EllipticCone :: MaxCurvature () const +{ + double a = vs.Length(); + double b = vl.Length(); + + return max2(b/(a*a),a/(b*b)); +} + +double EllipticCone :: MaxCurvatureLoc (const Point<3> & c, + double rad) const +{ +#ifdef JOACHIMxxx + cout << "max curv local" << endl; + return 0.02; +#endif + + // saubere Loesung wird noch notwendig !!! + double a = vs.Length(); + double b = vl.Length(); + return max2(b/(a*a),a/(b*b)); +} + +Point<3> EllipticCone :: GetSurfacePoint () const +{ + return a + vl; +} + + +void EllipticCone :: GetTriangleApproximation +(TriangleApproximation & tas, + const Box<3> & boundingbox, double facets) const +{ + int i, j; + double lg, bg; + int n = int(facets) + 1; + + Vec<3> nh = Cross(vl, vs); + nh.Normalize(); + Vec<3> vh = h*nh; + + double lvl = vl.Length(); + double lvs = vs.Length(); + Vec<3> nvl = (1.0/lvl)*vl; + Vec<3> nvs = (1.0/lvs)*vs; + + for ( j = 0; j <= n; j++ ) + for (i = 0; i <= n; i++) + { + lg = 2 *M_PI * double (i) /n; + bg = double(j) /n; + + Point<3> p = a + (bg *vh) + + (( lvl*(1 + (vlr -1)*bg) * cos(lg)) * nvl) + + (( lvs*(1 + (vlr -1)*bg)* sin(lg) ) * nvs); + + tas.AddPoint (p); + } + + for ( j = 0; j < n; j++) + for ( i = 0; i < n; i++) + { + int pi = i + (n+1) * j; + tas.AddTriangle (TATriangle (0, pi, pi+1, pi+n+2)); + tas.AddTriangle (TATriangle (0, pi, pi+n+2, pi+n+1)); + } +} diff --git a/libsrc/csg/algprim.hpp b/libsrc/csg/algprim.hpp index 9ca93cdd..98715fd3 100644 --- a/libsrc/csg/algprim.hpp +++ b/libsrc/csg/algprim.hpp @@ -365,10 +365,46 @@ namespace netgen }; +/// +/// Elliptic Cone +/// Josephat Kalezhi (kalezhi@cbu.ac.zm) +/// February 21st, 2018 +/// + /// + class EllipticCone : public QuadraticSurface + { + Point<3> a; + Vec<3> vl, vs; + double h, vlr; + public: + /// + EllipticCone (const Point<3> & aa, const Vec<3> & avl, + const Vec<3> & avs, double ah, double avlr); + static Primitive * CreateDefault (); + virtual void GetPrimitiveData (const char *& classname, Array & coeffs) const; + virtual void SetPrimitiveData (Array & coeffs); + + /// + virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const; + /// + virtual double HesseNorm () const; + virtual double MaxCurvature () const; + virtual double MaxCurvatureLoc (const Point<3> & /* c */ , + double /* rad */) const; + /// + virtual Point<3> GetSurfacePoint () const; + + virtual void GetTriangleApproximation (TriangleApproximation & tas, + const Box<3> & bbox, + double facets) const; + + private: + void CalcData(); + }; /** Torus diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index 892d754f..40dda44e 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -477,6 +477,15 @@ namespace netgen delete_them.Append(cone); } + else if(classname == "ellipticcone") + { + EllipticCone * ellipticcone = new EllipticCone(dummypoint,dummyvec,dummyvec,dummydouble,dummydouble); + + ellipticcone->SetPrimitiveData(coeffs); + + AddSurface(ellipticcone); + delete_them.Append(ellipticcone); + } else if(classname == "extrusionface") { ExtrusionFace * ef = diff --git a/libsrc/csg/csgparser.cpp b/libsrc/csg/csgparser.cpp index 8ac0a6fe..c202c9fd 100644 --- a/libsrc/csg/csgparser.cpp +++ b/libsrc/csg/csgparser.cpp @@ -39,6 +39,7 @@ namespace netgen { TOK_SPHERE, "sphere" }, { TOK_CYLINDER, "cylinder" }, { TOK_CONE, "cone" }, + { TOK_ELLIPTICCONE, "ellipticcone" }, { TOK_ELLIPTICCYLINDER, "ellipticcylinder" }, { TOK_ELLIPSOID, "ellipsoid" }, { TOK_ORTHOBRICK, "orthobrick" }, @@ -327,7 +328,19 @@ namespace netgen return new Solid (surf); } + case TOK_ELLIPTICCONE: + { + Point<3> a; + Vec<3> vl, vs; + double h, vlr; + scan.ReadNext(); + scan >> '(' >> a >> ';' >> vl >> ';' >> vs >> ';' >> h >>';' >> vlr >> ')'; + + OneSurfacePrimitive * surf = new EllipticCone ( a, vl, vs, h, vlr ); + geom->AddSurfaces (surf); + return new Solid (surf); + } case TOK_SPHERE: { diff --git a/libsrc/csg/csgparser.hpp b/libsrc/csg/csgparser.hpp index dfbd24ac..0895190a 100644 --- a/libsrc/csg/csgparser.hpp +++ b/libsrc/csg/csgparser.hpp @@ -27,7 +27,7 @@ namespace netgen enum PRIMITIVE_TYPE { TOK_SPHERE = 1, TOK_CYLINDER, TOK_PLANE, TOK_ELLIPTICCYLINDER, - TOK_ELLIPSOID, TOK_CONE, + TOK_ELLIPSOID, TOK_CONE, TOK_ELLIPTICCONE, TOK_ORTHOBRICK, TOK_POLYHEDRON, TOK_TORUS, TOK_TUBE, TOK_GENCYL, TOK_EXTRUSION, TOK_REVOLUTION, diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index ead2f6b5..d6d60e75 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -330,6 +330,23 @@ DLL_HEADER void ExportCSG(py::module &m) Solid * sol = new Solid(extr); return make_shared (sol); })); + m.def("EllipticCone", [](const Point<3>& a, const Vec<3>& v, const Vec<3>& w, + double h, double r) + { + auto ellcone = new EllipticCone(a,v,w,h,r); + auto sol = new Solid(ellcone); + return make_shared(sol); + }, py::arg("a"), py::arg("vl"), py::arg("vs"), py::arg("h"), py::arg("r"), + R"raw_string( +An elliptic cone, given by the point 'a' at the base of the cone along the main axis, +the vectors v and w of the long and short axis of the ellipse, respectively, +the height of the cone, h, and ratio of base long axis length to top long axis length, r + +Note: The elliptic cone has to be truncated by planes similar to a cone or an elliptic cylinder. +When r =1, the truncated elliptic cone becomes an elliptic cylinder. +When r tends to zero, the truncated elliptic cone tends to a full elliptic cone. +However, when r = 0, the top part becomes a point(tip) and meshing fails! +)raw_string"); m.def ("Or", FunctionPointer([](shared_ptr s1, shared_ptr s2) { diff --git a/tutorials/ellipticcone.geo b/tutorials/ellipticcone.geo new file mode 100644 index 00000000..6ccdd0ac --- /dev/null +++ b/tutorials/ellipticcone.geo @@ -0,0 +1,21 @@ +# +## An elliptic cone +# + +# An elliptic cone, given by the point c = (c x , c y , c z ) at the base of the cone along the main axis, +# the vectors v and w of the long and short axis of the ellipse, respectively, +# the height of the cone, h, and ratio of base long axis length to top long axis length, r: +# ellipticcone (c x , c y , c z ; v x , v y , v z ; w x , w y , w z; h; r) + +# Note: The elliptic cone has to be truncated by planes similar to a cone or an elliptic cylinder. +# When r =1, the truncated elliptic cone becomes an elliptic cylinder. +# When r tends to zero, the truncated elliptic cone tends to a full elliptic cone. +# However, when r = 0, the top part becomes a point(tip) and meshing fails! + +algebraic3d + +solid cutcone = ellipticcone ( 0, 0, 0; 5, 0, 0; 0, 2, 0; 5; 0.5) + and plane (0, 0, 0; 0, 0, -1) + and plane (0, 0, 5; 0, 0, 1); + +tlo cutcone;