diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index aa3386a9..f6818eba 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -479,61 +479,82 @@ namespace netgen mesh.LoadLocalMeshSize(mparam.meshsizefilename); } - void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array> & points, Array & params, int layer) + void GeometryEdge :: Divide(const MeshingParameters & mparam, const Mesh & mesh, Array> & points, Array & params) { - static Timer tdivedgesections("Divide edge sections"); - static Timer tdivide("Divide Edges"); - RegionTimer rt(tdivide); - // -------------------- DivideEdge ----------------- - static constexpr size_t divide_edge_sections = 10000; - double hvalue[divide_edge_sections+1]; - hvalue[0] = 0; + static Timer tdivedgesections("Divide edge sections"); + static Timer tdivide("Divide Edges"); + RegionTimer rt(tdivide); + // -------------------- DivideEdge ----------------- + tdivedgesections.Start(); + auto layer = properties.layer; + double safety = 0.5*(1.-mparam.grading); - Point<3> old_pt = edge->GetPoint(0.); - // calc local h for edge - tdivedgesections.Start(); - for(auto i : Range(divide_edge_sections)) + double lam = 0.0; + Point<3> p = GetPoint(0.0); + auto old_p = p; + Array hvalue, fine_params; + hvalue.Append(.0); + + while (lam<1. && hvalue.Size() < 20000) { + fine_params.Append(lam); + auto h = mesh.GetH(old_p, layer); + auto step = safety * h/GetTangent(lam).Length(); + lam += step; + lam = min2(lam, 1.0); + p = GetPoint(lam); + hvalue.Append((hvalue.Size()==0 ? 0.0 : hvalue.Last()) + 1./h * (p-old_p).Length()); + old_p = p; + } + + fine_params.Append(1.0); + + if(hvalue.Size()==20000 && lam<1.0) + cout << "Warning: Could not divide Edge" << endl; + + tdivedgesections.Stop(); + + auto n = hvalue.Size()-1; + int nsubedges = max2(1, int(floor(hvalue.Last()+0.5))); + points.SetSize(nsubedges-1); + params.SetSize(nsubedges+1); + + int i1 = 0; + for(auto i : Range(1,nsubedges)) + { + auto h_target = i*hvalue.Last()/nsubedges; + while(hvalue[i1]GetPoint(double(i+1)/divide_edge_sections); - hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt, layer) * (pt-old_pt).Length(); - old_pt = pt; + points.SetSize(i-1); + params.SetSize(i+1); + cout << "divide edge: local h too small" << endl; + break; } - int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5))); - tdivedgesections.Stop(); - points.SetSize(nsubedges-1); - params.SetSize(nsubedges+1); - int i = 1; - int i1 = 0; - do - { - if (hvalue[i1]/hvalue[divide_edge_sections]*nsubedges >= i) - { - params[i] = (double(i1)/divide_edge_sections); - points[i-1] = MeshPoint(edge->GetPoint(params[i])); - i++; - } - i1++; - if (i1 > divide_edge_sections) - { - nsubedges = i; - points.SetSize(nsubedges-1); - params.SetSize(nsubedges+1); - cout << "divide edge: local h too small" << endl; - } + // interpolate lam between points + auto lam0 = fine_params[i1-1]; + auto lam1 = fine_params[i1]; + auto h0 = hvalue[i1-1]; + auto h1 = hvalue[i1]; - } while(i < nsubedges); + auto fac = (h_target-h0)/(h1-h0); + auto lam = lam0 + fac*(lam1-lam0); + params[i] = lam; + points[i-1] = MeshPoint(GetPoint(params[i])); + } - params[0] = 0.; - params[nsubedges] = 1.; + params[0] = 0.; + params[nsubedges] = 1.; - if(params[nsubedges] <= params[nsubedges-1]) - { - cout << "CORRECTED" << endl; - points.SetSize (nsubedges-2); - params.SetSize (nsubedges); - params[nsubedges-1] = 1.; - } + if(params[nsubedges] <= params[nsubedges-1]) + { + cout << "CORRECTED" << endl; + points.SetSize (nsubedges-2); + params.SetSize (nsubedges); + params[nsubedges-1] = 1.; + } } void NetgenGeometry :: FindEdges(Mesh& mesh, @@ -616,7 +637,7 @@ namespace netgen } else { - DivideEdge(edge, mparam, mesh, edge_points, params, edge->properties.layer); + edge->Divide(mparam, mesh, edge_points, params); } } else diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 6b38f11f..bb9a508c 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -110,6 +110,7 @@ namespace netgen } virtual Vec<3> GetTangent(double t) const = 0; virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override; + virtual void Divide(const MeshingParameters & mparam, const Mesh & mesh, Array> & points, Array & params); }; class DLL_HEADER GeometryFace : public GeometryShape diff --git a/libsrc/occ/occ_edge.cpp b/libsrc/occ/occ_edge.cpp index 1685aa9b..a3a2b06e 100644 --- a/libsrc/occ/occ_edge.cpp +++ b/libsrc/occ/occ_edge.cpp @@ -69,6 +69,7 @@ namespace netgen gp_Pnt p; gp_Vec v; curve->D1(t, p, v); - return occ2ng(v); + return occ2ng(v) * (s1-s0); } + } diff --git a/libsrc/occ/occ_face.cpp b/libsrc/occ/occ_face.cpp index 88fd3fb9..2e60469a 100644 --- a/libsrc/occ/occ_face.cpp +++ b/libsrc/occ/occ_face.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include "occ_edge.hpp" #include "occ_face.hpp" @@ -111,9 +112,16 @@ namespace netgen for(auto i : Range(2)) { + // take uv from CurveOnSurface as start value but project again for better accuracy + // (cof->Value yields wrong values (outside of surface) for complicated faces auto uv = cof->Value(s[i]); - seg.epgeominfo[i].u = uv.X(); - seg.epgeominfo[i].v = uv.Y(); + PointGeomInfo gi; + gi.u = uv.X(); + gi.v = uv.Y(); + Point<3> pproject = mesh[seg[i]]; + ProjectPointGI(pproject, gi); + seg.epgeominfo[i].u = gi.u; + seg.epgeominfo[i].v = gi.v; } bool do_swap = ORIENTATION == REVERSED; @@ -234,7 +242,11 @@ namespace netgen double OCCFace::GetCurvature(const PointGeomInfo& gi) const { - throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__)); + BRepAdaptor_Surface sf(face, Standard_True); + BRepLProp_SLProps prop2(sf, 2, 1e-5); + prop2.SetParameters (gi.u, gi.v); + return max(fabs(prop2.MinCurvature()), + fabs(prop2.MaxCurvature())); } void OCCFace::RestrictH(Mesh& mesh, const MeshingParameters& mparam) const diff --git a/libsrc/occ/occ_utils.hpp b/libsrc/occ/occ_utils.hpp index 91cf78a3..2f55d405 100644 --- a/libsrc/occ/occ_utils.hpp +++ b/libsrc/occ/occ_utils.hpp @@ -276,13 +276,12 @@ namespace netgen } }; - - inline gp_Pnt Center (TopoDS_Shape shape) + inline auto Properties (TopoDS_Shape shape) { GProp_GProps props; double tol; switch (shape.ShapeType()) - { + { case TopAbs_SOLID: case TopAbs_COMPOUND: case TopAbs_COMPSOLID: @@ -298,8 +297,18 @@ namespace netgen BRepGProp::LinearProperties(shape, props, tol); break; default: BRepGProp::LinearProperties(shape, props); - } - return props.CentreOfMass(); + } + return props; + } + + inline gp_Pnt Center (TopoDS_Shape shape) + { + return Properties(shape).CentreOfMass(); + } + + inline double Mass (TopoDS_Shape shape) + { + return Properties(shape).Mass(); } } diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index d2903765..517e1ebb 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -537,9 +537,7 @@ namespace netgen multithread.percent = 100 * (i-1)/double(nedges); if (BRep_Tool::Degenerated(e)) continue; - GProp_GProps system; - BRepGProp::LinearProperties(e, system); - double len = system.Mass(); + double len = Mass(e); if (len < mincurvelength) { @@ -587,7 +585,7 @@ namespace netgen localh = min2(localh, OCCGeometry::GetProperties(e).maxh); maxedgelen = max (maxedgelen, len); minedgelen = min (minedgelen, len); - int maxj = max((int) ceil(len/localh), 2); + int maxj = max((int) ceil(len/localh)*2, 2); for (int j = 0; j <= maxj; j++) { diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index cb17666e..8959a19d 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -715,23 +715,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) .def("Properties", [] (const TopoDS_Shape & shape) { - GProp_GProps props; - switch (shape.ShapeType()) - { - case TopAbs_FACE: - case TopAbs_SHELL: - BRepGProp::SurfaceProperties (shape, props); break; - case TopAbs_SOLID: - case TopAbs_COMPOUND: - case TopAbs_COMPSOLID: - BRepGProp::VolumeProperties (shape, props); break; - default: - BRepGProp::LinearProperties(shape, props); - // throw Exception("Properties implemented only for FACE"); - } - double mass = props.Mass(); - gp_Pnt center = props.CentreOfMass(); - return tuple( py::cast(mass), py::cast(center) ); + auto props = Properties(shape); + return tuple( py::cast(props.Mass()), py::cast(props.CentreOfMass()) ); }, "returns tuple of shape properties, currently ('mass', 'center'") .def_property_readonly("center", [](const TopoDS_Shape & shape) { @@ -739,20 +724,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) }, "returns center of gravity of shape") .def_property_readonly("mass", [](const TopoDS_Shape & shape) { - GProp_GProps props; - switch (shape.ShapeType()) - { - case TopAbs_FACE: - case TopAbs_SHELL: - BRepGProp::SurfaceProperties (shape, props); break; - case TopAbs_SOLID: - case TopAbs_COMPOUND: - case TopAbs_COMPSOLID: - BRepGProp::VolumeProperties (shape, props); break; - default: - BRepGProp::LinearProperties(shape, props); - } - return props.Mass(); + return Mass(shape); }, "returns mass of shape, what is length, face, or volume") .def("Move", [](const TopoDS_Shape & shape, const gp_Vec v)