diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index f5f85ea7..fe0c8bd0 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/libsrc/geom2d/csg2d.cpp @@ -1292,20 +1292,38 @@ void CreateResult(Solid2d & sp, Solid2d & sr, bool UNION) } } +// Check if vertex v is not necessary (i.e. is on the line v->prev, v->next and has same info as v->prev and no pinfo +bool canRemoveVertex( Vertex * v ) +{ + return false; + if(v->spline) + return false; + if(v->pinfo.name != POINT_NAME_DEFAULT) + return false; + if(v->pinfo.maxh != MAXH_DEFAULT) + return false; + + if(v->info.bc != v->prev->info.bc || v->info.maxh != v->prev->info.maxh ) + return false; + + if(fabs(Area(*v->prev,*v,*v->next)) >= EPSILON) + return false; + + return true; +} + void CleanUpResult(Solid2d & sr) { auto & RR = sr.polys; for (Loop& R : RR) { - while ( (R.first.get() != NULL) && (fabs(Area(*R.first->prev,*R.first,*R.first->next)) < EPSILON) ) + while ( (R.first.get() != NULL) && canRemoveVertex(R.first.get())) R.Remove(R.first.get()); if (R.first.get() != NULL) for (Vertex* V : R.Vertices(ALL)) - if (!V->spline && !V->prev->spline && fabs(Area(*V->prev,*V,*V->next)) < EPSILON) - { + if (canRemoveVertex(V)) R.Remove(V); - } } for (int i = RR.Size()-1; i>=0; i--) if(RR[i].Size()==0) @@ -1547,6 +1565,39 @@ Solid2d ClipSolids ( Solid2d && s1, Solid2d && s2, char op) return std::move(res); } +Vertex* Loop :: getNonIntersectionVertex() + { + for (Vertex* v : Vertices(ALL)) + if (!v->is_intersection) + return(v); + + // no non-intersection vertex found -> generate and return temporary vertex + for (Vertex* v : Vertices(ALL)) + // make sure that edge from V to V->next is not collinear with other polygon + if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) ) + { + // add edge midpoint as temporary vertex + if(v->spline) + { + auto p = v->spline->GetPoint(0.5); + auto s = *v->spline; + v->spline = Split(s, 0, 0.5); + auto vnew = v->Insert(p); + vnew->info = v->info; + vnew->spline = Split(s, 0.5, 1.0); + return vnew; + } + else + { + auto p = Center(*v, *v->next); + auto vnew = v->Insert(p); + vnew->info = v->info; + return vnew; + } + } + return(NULL); + } + bool Loop :: IsInside( Point<2> r ) const { int w = 0; @@ -1652,11 +1703,14 @@ Solid2d & Solid2d :: Move( Vec<2> v ) return Transform( [v](Point<2> p) -> Point<2> { return p+v; } ); } -Solid2d & Solid2d :: Scale( double sx, double sy ) +Solid2d & Solid2d :: Scale( double s ) { - if(sy==0.0) - sy=sx; - return Transform( [sx,sy](Point<2> p) -> Point<2> { return{p[0]*sx, p[1]*sy}; } ); + return Transform( [s](Point<2> p) -> Point<2> { return{p[0]*s, p[1]*s}; } ); +} + +Solid2d & Solid2d :: Scale( Vec<2> s ) +{ + return Transform( [s](Point<2> p) -> Point<2> { return{p[0]*s[0], p[1]*s[1]}; } ); } Solid2d & Solid2d :: RotateRad( double ang, Point<2> center ) @@ -1669,8 +1723,8 @@ Solid2d & Solid2d :: RotateRad( double ang, Point<2> center ) p -= c; double x = p[0]; double y = p[1]; - p[0] = cosa*x+sina*y; - p[1] = -sina*x+cosa*y; + p[0] = cosa*x-sina*y; + p[1] = sina*x+cosa*y; p += c; return p; } ); diff --git a/libsrc/geom2d/csg2d.hpp b/libsrc/geom2d/csg2d.hpp index 118060f8..afc1e654 100644 --- a/libsrc/geom2d/csg2d.hpp +++ b/libsrc/geom2d/csg2d.hpp @@ -29,10 +29,10 @@ enum IntersectionType T_INTERSECTION_Q, T_INTERSECTION_P, V_INTERSECTION, - X_OVERLAP, - T_OVERLAP_Q, - T_OVERLAP_P, - V_OVERLAP + X_OVERLAP, // Q0 -- P1 -- Q1 -- P0 (different direction) + T_OVERLAP_Q, // same direction or P inside Q + T_OVERLAP_P, // same direction or Q inside P + V_OVERLAP // one common point }; enum IntersectionLabel @@ -588,23 +588,7 @@ struct Loop // // return and insert a non-intersection vertex // - Vertex* getNonIntersectionVertex() - { - for (Vertex* v : Vertices(ALL)) - if (!v->is_intersection) - return(v); - - // no non-intersection vertex found -> generate and return temporary vertex - for (Vertex* v : Vertices(ALL)) - // make sure that edge from V to V->next is not collinear with other polygon - if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) ) - { - // add edge midpoint as temporary vertex - auto p = Center(*v, *v->next); - return v->Insert(p); - } - return(NULL); - } + Vertex* getNonIntersectionVertex(); void SetBC(string bc) { @@ -694,11 +678,12 @@ struct Solid2d } Solid2d & Move( Vec<2> v ); - Solid2d & Scale( double sx, double sy=0.0 ); + Solid2d & Scale( double s ); + Solid2d & Scale( Vec<2> s ); Solid2d & RotateRad( double ang, Point<2> center = {0,0} ); Solid2d & RotateDeg( double ang, Point<2> center = {0,0} ) { - return RotateRad( ang/180.*M_PI ); + return RotateRad( ang/180.*M_PI, center ); } Solid2d & BC(string bc) diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index b8f5a871..03f57398 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -414,16 +414,9 @@ DLL_HEADER void ExportGeom2d(py::module &m) .def("Copy", [](Solid2d & self) -> Solid2d { return self; }) .def("Move", &Solid2d::Move) - .def("Scale", &Solid2d::Scale) - .def("Rotate", [](Solid2d & self, optional rad, optional deg, Point<2> center ) - { - if(rad) - self.RotateRad(*rad, center); - if(deg) - self.RotateDeg(*deg, center); - return self; - }, py::arg("rad")=nullopt, py::arg("deg")=nullopt, py::arg("center")=Point<2>{0,0}) - + .def("Scale", static_cast(&Solid2d::Scale)) + .def("Scale", static_cast)>(&Solid2d::Scale)) + .def("Rotate", &Solid2d::RotateDeg, py::arg("angle"), py::arg("center")=Point<2>{0,0}) ; diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 2dba11ea..9987ed0f 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -211,6 +211,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) py::class_> (m, "Vec2d") .def(py::init()) + .def(py::init( [] (std::pair xy) + { + return Vec<2>{xy.first, xy.second}; + })) .def ("__str__", &ToString>) .def(py::self==py::self) .def(py::self+py::self) @@ -223,6 +227,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def("__len__", [](Vec<2>& /*unused*/) { return 2; }) ; + py::implicitly_convertible>(); + py::class_> (m, "Vec3d") .def(py::init()) .def(py::init([](py::tuple v) @@ -1261,27 +1267,45 @@ project_boundaries : Optional[str] = None }), py::arg("mapping")) .def(NGSPickle()) .def("GenerateMesh", [](shared_ptr geo, - bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames) + bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames, py::list py_hppts, py::list py_hpbnd) { if (py::len(py_bbbpts) != py::len(py_bbbnames)) throw Exception("In SurfaceGeometry::GenerateMesh bbbpts and bbbnames do not have same lengths."); Array> bbbpts(py::len(py_bbbpts)); Array bbbname(py::len(py_bbbpts)); + Array> hppts(py::len(py_hppts)); + Array hpptsfac(py::len(py_hppts)); + Array hpbnd(py::len(py_hpbnd)); + Array hpbndfac(py::len(py_hpbnd)); for(int i = 0; i(py_bbbpts[i])(); bbbpts[i] = Point<3>(py::extract(pnt[0])(),py::extract(pnt[1])(),py::extract(pnt[2])()); bbbname[i] = py::extract(py_bbbnames[i])(); } + for(int i = 0; i(py_hppts[i])(); + hppts[i] = Point<3>(py::extract(pnt[0])(),py::extract(pnt[1])(),py::extract(pnt[2])()); + //hpptsfac[i] = py::len(pnt) > 3 ? py::extract(pnt[3])() : 0.0; + hpptsfac[i] = py::extract(pnt[3])(); + } + + for(int i = 0; i(py_hpbnd[i])(); + hpbnd[i] = py::extract(bnd[0])(); + hpbndfac[i] = py::extract(bnd[1])(); + } auto mesh = make_shared(); SetGlobalMesh (mesh); mesh->SetGeometry(geo); ng_geometry = geo; - auto result = geo->GenerateMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname); + auto result = geo->GenerateMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname, hppts, hpptsfac, hpbnd, hpbndfac); if(result != 0) throw Exception("SurfaceGeometry: Meshing failed!"); return mesh; - }, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list()) + }, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list(), py::arg("hppts")=py::list(), py::arg("hpbnd")=py::list()) ; ; diff --git a/libsrc/meshing/surfacegeom.cpp b/libsrc/meshing/surfacegeom.cpp index 40d272a8..14a265d3 100644 --- a/libsrc/meshing/surfacegeom.cpp +++ b/libsrc/meshing/surfacegeom.cpp @@ -45,10 +45,10 @@ namespace netgen prr = p[1]+2*eps; pll = p[1]-2*eps; - dr = GetTangentVectors(u, pr); - dl = GetTangentVectors(u, pl); - drr = GetTangentVectors(u, prr); - dll = GetTangentVectors(u, pll); + GetTangentVectors(u, pr, dr); + GetTangentVectors(u, pl, dl); + GetTangentVectors(u, prr, drr); + GetTangentVectors(u, pll, dll); f_vv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]); } @@ -74,11 +74,31 @@ namespace netgen return tang; } + void SurfaceGeometry :: GetTangentVectors(double u, double v, Array>& tang) const + { + + Point<2> pru = Point<2>(u+eps,v); + Point<2> plu = Point<2>(u-eps,v); + Point<2> prru = Point<2>(u+2*eps,v); + Point<2> pllu = Point<2>(u-2*eps,v); + + Point<2> prv = Point<2>(u,v+eps); + Point<2> plv = Point<2>(u,v-eps); + Point<2> prrv = Point<2>(u,v+2*eps); + Point<2> pllv = Point<2>(u,v-2*eps); + + + tang[0] = 1/(12.0*eps)*( 8.0*func(pru) - 8.0*func(plu) - func(prru) + func(pllu) ); + tang[1] = 1/(12.0*eps)*( 8.0*func(prv) - 8.0*func(plv) - func(prrv) + func(pllv) ); + } + + Vec<3> SurfaceGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const { Array> tang = GetTangentVectors(gi->u, gi->v); auto normal = Cross(tang[0], tang[1]); - return Cross(tang[0], tang[1]); + normal.Normalize(); + return normal; } @@ -97,14 +117,13 @@ namespace netgen bool SurfaceGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const { - Array> tangs; + Array> tangs(2); Vec<3> diff, f_uu, f_vv, f_uv; Vec<2> r, dx; - double norm_r, det, energy=0.0, new_energy=0.0, alpha=2.0,u=0.0,v=0.0; + double norm_r, det, energy=0.0, new_energy=0.0, alpha=2.0,u=0.0,v=0.0,maxerr=1e-16; Mat<2,2> mat, inv; - int num=0, maxit=20; - double damping=0.2; - + int num=0, maxit=25; + double damping=0.5; //Solve minimization problem // argmin_(u,v) 0.5*\| f(u,v)-p\|^2 @@ -115,7 +134,7 @@ namespace netgen do { num++; - tangs = GetTangentVectors(gi.u, gi.v); + GetTangentVectors(gi.u, gi.v,tangs); diff = func(Point<2>(gi.u, gi.v)) - Vec<3>(p); energy = diff.Length2(); r = Vec<2>( diff*tangs[0], diff*tangs[1] ); @@ -146,12 +165,13 @@ namespace netgen while (alpha > 1e-10 && new_energy > energy+1e-14); if (alpha <= 1e-10) throw Exception("In SurfaceGeometry::ProjectPointGI: Linesearch min alpha reached!"); + gi.u = u; gi.v = v; } - while ( norm_r > 1e-12 && num < maxit); + while ( norm_r > maxerr && num < maxit); //Stay in reference domain [0,1]^2 if (gi.u < 0 || gi.u > 1 || gi.v < 0 || gi.v > 1) @@ -179,19 +199,13 @@ namespace netgen void SurfaceGeometry :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint, int surfi1, int surfi2, const EdgePointGeomInfo & ap1, const EdgePointGeomInfo & ap2, Point<3> & newp, EdgePointGeomInfo & newgi) const { - newp = p1+secpoint*(p2-p1); - - PointGeomInfo pgi; - pgi.u = ap1.u+secpoint*(ap2.u-ap1.u); - pgi.v = ap1.v+secpoint*(ap2.v-ap1.v); - - ProjectPointGI(surfi1, newp, pgi); - - newgi.u = pgi.u; - newgi.v = pgi.v; + newgi.u = ap1.u+secpoint*(ap2.u-ap1.u); + newgi.v = ap1.v+secpoint*(ap2.v-ap1.v); newgi.edgenr = ap1.edgenr; newgi.body = -1; newgi.dist = -1.0; + + newp = Point<3>(func(Point<2>(newgi.u, newgi.v))); } void SurfaceGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, @@ -200,16 +214,16 @@ namespace netgen const PointGeomInfo & gi2, Point<3> & newp, PointGeomInfo & newgi) const { - newp = p1+secpoint*(p2-p1); - newgi.u = gi1.u+secpoint*(gi2.u-gi1.u); newgi.v = gi1.v+secpoint*(gi2.v-gi1.v); newgi.trignum = -1; - ProjectPointGI(surfi, newp, newgi); + newp = Point<3>(func(Point<2>(newgi.u, newgi.v))); + //newp = p1+secpoint*(p2-p1); + //ProjectPointGI(surfi, newp, newgi); } - int SurfaceGeometry :: GenerateMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames) + int SurfaceGeometry :: GenerateMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames, const Array>& hppoints, const Array& hpptsfac, const Array& hpbnd, const Array& hpbndfac) { mesh->SetDimension(3); @@ -241,6 +255,16 @@ namespace netgen indbbbpts[k] = pids[pids.Size()-1]; } } + + for (int k = 0; k < hppoints.Size(); k++) + { + auto diff = pnt - hppoints[k]; + if(diff.Length2() < 1e-14) + { + (*mesh)[pids[pids.Size()-1]].Singularity(hpptsfac[k]); + } + + } } for (bool f : found) @@ -318,8 +342,19 @@ namespace netgen Segment seg; seg.si = 1; seg.edgenr = 1; - seg.epgeominfo[0].edgenr = 1; - seg.epgeominfo[1].edgenr = 1; + seg.epgeominfo[0].edgenr = 0; + seg.epgeominfo[1].edgenr = 0; + //for hp refinement + seg.singedge_left = 0; + seg.singedge_right = 0; + for (size_t i=0; i < hpbnd.Size(); i++) + { + if (hpbnd[i] == "bottom") + { + seg.singedge_left = hpbndfac[i]; + seg.singedge_right = hpbndfac[i]; + } + } // needed for codim2 in 3d seg.edgenr = 1; for(int i=0; i < nx; i++) @@ -341,6 +376,18 @@ namespace netgen seg.si = 2; seg.edgenr = 2; + seg.singedge_left = 0; + seg.singedge_right = 0; + + for (size_t i=0; i < hpbnd.Size(); i++) + { + if (hpbnd[i] == "right") + { + seg.singedge_left = hpbndfac[i]; + seg.singedge_right = hpbndfac[i]; + } + } + for(int i=0; i> GetTangentVectors(double u, double v) const; + void GetTangentVectors(double u, double v, Array>& tang) const; + + virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override; virtual PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override; @@ -59,7 +62,7 @@ namespace netgen const PointGeomInfo & gi2, Point<3> & newp, PointGeomInfo & newgi) const override; - int GenerateMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames); + int GenerateMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames, const Array>& hppoints, const Array& hpptsfac, const Array& hpbnd, const Array& hpbndfac); }; diff --git a/tests/pytest/test_csg2d.py b/tests/pytest/test_csg2d.py index d962085c..06e875a9 100644 --- a/tests/pytest/test_csg2d.py +++ b/tests/pytest/test_csg2d.py @@ -5,7 +5,7 @@ from pytest import approx def test_two_circles(): c1 = Circle(center=(0,0), radius=1) - c2 = c1.Rotate(deg=45) + c2 = c1.Rotate(45) s = c1*c2 geo = CSG2d() geo.Add(s)