diff --git a/CMakeLists.txt b/CMakeLists.txt index d573504d..6e599d3e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,6 +34,8 @@ option( USE_SUPERBUILD "use ccache" ON) set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_modules") +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + if(APPLE) set(INSTALL_DIR_DEFAULT /Applications/Netgen.app) else(APPLE) diff --git a/libsrc/core/utils.hpp b/libsrc/core/utils.hpp index 35877d42..3645ea18 100644 --- a/libsrc/core/utils.hpp +++ b/libsrc/core/utils.hpp @@ -103,6 +103,22 @@ namespace ngcore void operator delete[] (void * p) { aligned_free(p); } }; + // checks if string starts with sequence + inline bool StartsWith(const std::string& str, const std::string& start) + { + if(start.size() > str.size()) + return false; + return std::equal(start.begin(), start.end(), str.begin()); + } + + // checks if string ends with sequence + inline bool EndsWith(const std::string& str, const std::string& end) + { + if(end.size() > str.size()) + return false; + return std::equal(end.rbegin(), end.rend(), str.rbegin()); + } + } // namespace ngcore #endif // NETGEN_CORE_UTILS_HPP diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 603ff98c..7319d8e7 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -3,6 +3,7 @@ #include <../general/ngpython.hpp> #include #include +#include "../meshing/python_mesh.hpp" using namespace netgen; @@ -693,26 +694,25 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails! res["max"] = MoveToNumpy(max); return res; }, py::call_guard()) - ; - - m.def("GenerateMesh", FunctionPointer - ([](shared_ptr geo, MeshingParameters & param) + .def("GenerateMesh", [](shared_ptr geo, + MeshingParameters* pars, py::kwargs kwargs) { - auto dummy = make_shared(); - SetGlobalMesh (dummy); - dummy->SetGeometry(geo); + MeshingParameters mp; + if(pars) mp = *pars; + { + py::gil_scoped_acquire aq; + CreateMPfromKwargs(mp, kwargs); + } + auto mesh = make_shared(); + SetGlobalMesh (mesh); + mesh->SetGeometry(geo); ng_geometry = geo; geo->FindIdenticSurfaces(1e-8 * geo->MaxSize()); - try - { - geo->GenerateMesh (dummy, param); - } - catch (NgException ex) - { - cout << "Caught NgException: " << ex.What() << endl; - } - return dummy; - }),py::call_guard()) + geo->GenerateMesh (mesh, mp); + return mesh; + }, py::arg("mp") = nullptr, + meshingparameter_description.c_str(), + py::call_guard()) ; m.def("Save", FunctionPointer diff --git a/libsrc/general/ngarray.hpp b/libsrc/general/ngarray.hpp index 4b66d65d..5312edf2 100644 --- a/libsrc/general/ngarray.hpp +++ b/libsrc/general/ngarray.hpp @@ -293,6 +293,11 @@ namespace netgen size = nsize; } + void SetSize0() + { + size = 0; + } + /// Change physical size. Keeps logical size. Keeps contents. void SetAllocSize (size_t nallocsize) { diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index d04b8cd7..f0e56598 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -2,6 +2,7 @@ #include <../general/ngpython.hpp> #include +#include "../meshing/python_mesh.hpp" #include #include @@ -362,17 +363,26 @@ DLL_HEADER void ExportGeom2d(py::module &m) //cout << i << " : " << self.splines[i]->GetPoint(0.1) << " , " << self.splines[i]->GetPoint(0.5) << endl; } })) - .def("GenerateMesh", [](shared_ptr self, MeshingParameters & mparam) + // If we change to c++17 this can become optional + .def("GenerateMesh", [](shared_ptr self, + MeshingParameters* pars, py::kwargs kwargs) { - shared_ptr mesh = make_shared (); + MeshingParameters mp; + if(pars) mp = *pars; + { + py::gil_scoped_acquire aq; + CreateMPfromKwargs(mp, kwargs); + } + auto mesh = make_shared(); mesh->SetGeometry(self); SetGlobalMesh (mesh); ng_geometry = self; - self->GenerateMesh(mesh, mparam); + self->GenerateMesh(mesh, mp); return mesh; - },py::call_guard()) - - ; + }, py::arg("mp") = nullptr, + py::call_guard(), + meshingparameter_description.c_str()) + ; } diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp index 6e06fc82..081422d3 100644 --- a/libsrc/include/nginterface_v2.hpp +++ b/libsrc/include/nginterface_v2.hpp @@ -284,7 +284,8 @@ namespace netgen int GetDimension() const; int GetNLevels() const; - + size_t GetNVLevel (int level) const; + int GetNElements (int dim) const; int GetNNodes (int nt) const; diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index 1896583a..7c647119 100644 --- a/libsrc/interface/nginterface.cpp +++ b/libsrc/interface/nginterface.cpp @@ -1724,7 +1724,9 @@ void Ng_SetSurfaceElementOrders (int enr, int ox, int oy) int Ng_GetNLevels () { - return (mesh) ? mesh->mglevels : 0; + if (!mesh) return 0; + return max(size_t(1), mesh -> level_nv.Size()); + // return (mesh) ? mesh->mglevels : 0; } diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index 1d8f9366..a961e999 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -130,7 +130,15 @@ namespace netgen int Ngx_Mesh :: GetNLevels() const { - return mesh -> mglevels; + return max(size_t(1), mesh -> level_nv.Size()); + } + + size_t Ngx_Mesh :: GetNVLevel(int level) const + { + if (level >= mesh->level_nv.Size()) + return mesh->GetNV(); + else + return mesh->level_nv[level]; } int Ngx_Mesh :: GetNElements (int dim) const @@ -1155,10 +1163,8 @@ namespace netgen biopt.refine_hp = 1; biopt.task_manager = task_manager; biopt.tracer = tracer; - - const Refinement & ref = mesh->GetGeometry()->GetRefinement(); - ref.Bisect (*mesh, biopt); + mesh->GetGeometry()->GetRefinement().Bisect (*mesh, biopt); (*tracer)("call updatetop", false); mesh -> UpdateTopology(task_manager, tracer); (*tracer)("call updatetop", true); diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 8cb57bcd..6390fff0 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -2755,14 +2755,15 @@ namespace netgen inf.close(); } - - - if (mesh.mglevels == 1 || idmaps.Size() > 0) - BisectTetsCopyMesh(mesh, NULL, opt, idmaps, refelementinfofileread); - - mesh.ComputeNVertices(); - + + // if (mesh.mglevels == 1 || idmaps.Size() > 0) + if (mesh.level_nv.Size() == 0) // || idmaps.Size() ???? + { + BisectTetsCopyMesh(mesh, NULL, opt, idmaps, refelementinfofileread); + mesh.level_nv.Append (mesh.GetNV()); + } + int np = mesh.GetNV(); mesh.SetNP(np); @@ -2773,7 +2774,7 @@ namespace netgen // int initnp = np; // int maxsteps = 3; - mesh.mglevels++; + // mesh.mglevels++; /* if (opt.refinementfilename || opt.usemarkedelements) @@ -3807,7 +3808,8 @@ namespace netgen // write multilevel hierarchy to mesh: np = mesh.GetNP(); mesh.mlbetweennodes.SetSize(np); - if (mesh.mglevels <= 2) + // if (mesh.mglevels <= 2) + if (mesh.level_nv.Size() <= 1) { PrintMessage(4,"RESETTING mlbetweennodes"); for (int i = 1; i <= np; i++) @@ -3817,6 +3819,9 @@ namespace netgen } } + mesh.level_nv.Append (np); + + /* for (i = 1; i <= cutedges.GetNBags(); i++) for (j = 1; j <= cutedges.GetBagSize(i); j++) @@ -3982,11 +3987,12 @@ namespace netgen } } - + // Check/Repair static bool repaired_once; - if(mesh.mglevels == 1) + // if(mesh.mglevels == 1) + if(mesh.level_nv.Size() == 1) repaired_once = false; //mesh.Save("before.vol"); diff --git a/libsrc/meshing/curvedelems.hpp b/libsrc/meshing/curvedelems.hpp index f0d5d21c..f1a732a0 100644 --- a/libsrc/meshing/curvedelems.hpp +++ b/libsrc/meshing/curvedelems.hpp @@ -49,7 +49,7 @@ public: int GetOrder () { return order; } - virtual void DoArchive(Archive& ar) + void DoArchive(Archive& ar) { if(ar.Input()) buildJacPols(); diff --git a/libsrc/meshing/delaunay2d.cpp b/libsrc/meshing/delaunay2d.cpp index 0fc535dc..a61016e2 100644 --- a/libsrc/meshing/delaunay2d.cpp +++ b/libsrc/meshing/delaunay2d.cpp @@ -81,12 +81,12 @@ namespace netgen Box<3> bbox ( Box<3>::EMPTY_BOX ); double maxh = 0; - for (int i = 0; i < adfront->GetNFL(); i++) + for (int i = 0; i < adfront.GetNFL(); i++) { - const FrontLine & line = adfront->GetLine (i); + const FrontLine & line = adfront.GetLine (i); - const Point<3> & p1 = adfront->GetPoint(line.L().I1()); - const Point<3> & p2 = adfront->GetPoint(line.L().I2()); + const Point<3> & p1 = adfront.GetPoint(line.L().I1()); + const Point<3> & p2 = adfront.GetPoint(line.L().I2()); maxh = max (maxh, Dist (p1, p2)); @@ -115,12 +115,12 @@ namespace netgen { mesh.LocalHFunction().ClearFlags(); - for (int i = 0; i < adfront->GetNFL(); i++) + for (int i = 0; i < adfront.GetNFL(); i++) { - const FrontLine & line = adfront->GetLine(i); + const FrontLine & line = adfront.GetLine(i); - Box<3> bbox (adfront->GetPoint (line.L().I1())); - bbox.Add (adfront->GetPoint (line.L().I2())); + Box<3> bbox (adfront.GetPoint (line.L().I1())); + bbox.Add (adfront.GetPoint (line.L().I2())); double filld = filldist * bbox.Diam(); @@ -130,7 +130,7 @@ namespace netgen } - mesh.LocalHFunction().FindInnerBoxes (adfront, NULL); + mesh.LocalHFunction().FindInnerBoxes (&adfront, NULL); npoints.SetSize(0); mesh.LocalHFunction().GetInnerPoints (npoints); @@ -162,14 +162,14 @@ namespace netgen if (meshbox.IsIn (npoints.Get(i))) { PointIndex gpnum = mesh.AddPoint (npoints.Get(i)); - adfront->AddPoint (npoints.Get(i), gpnum); + adfront.AddPoint (npoints.Get(i), gpnum); if (debugparam.slowchecks) { (*testout) << npoints.Get(i) << endl; Point<2> p2d (npoints.Get(i)(0), npoints.Get(i)(1)); - if (!adfront->Inside(p2d)) + if (!adfront.Inside(p2d)) { cout << "add outside point" << endl; (*testout) << "outside" << endl; @@ -187,29 +187,29 @@ namespace netgen loch2.ClearFlags(); - for (int i = 0; i < adfront->GetNFL(); i++) + for (int i = 0; i < adfront.GetNFL(); i++) { - const FrontLine & line = adfront->GetLine(i); + const FrontLine & line = adfront.GetLine(i); - Box<3> bbox (adfront->GetPoint (line.L().I1())); - bbox.Add (adfront->GetPoint (line.L().I2())); + Box<3> bbox (adfront.GetPoint (line.L().I1())); + bbox.Add (adfront.GetPoint (line.L().I2())); loch2.SetH (bbox.Center(), bbox.Diam()); } - for (int i = 0; i < adfront->GetNFL(); i++) + for (int i = 0; i < adfront.GetNFL(); i++) { - const FrontLine & line = adfront->GetLine(i); + const FrontLine & line = adfront.GetLine(i); - Box<3> bbox (adfront->GetPoint (line.L().I1())); - bbox.Add (adfront->GetPoint (line.L().I2())); + Box<3> bbox (adfront.GetPoint (line.L().I1())); + bbox.Add (adfront.GetPoint (line.L().I2())); bbox.Increase (filldist * bbox.Diam()); loch2.CutBoundary (bbox); } - loch2.FindInnerBoxes (adfront, NULL); + loch2.FindInnerBoxes (&adfront, NULL); // outer points : smooth mesh-grading npoints.SetSize(0); @@ -220,7 +220,7 @@ namespace netgen if (meshbox.IsIn (npoints.Get(i))) { PointIndex gpnum = mesh.AddPoint (npoints.Get(i)); - adfront->AddPoint (npoints.Get(i), gpnum); + adfront.AddPoint (npoints.Get(i), gpnum); } } @@ -257,11 +257,11 @@ namespace netgen // face bounding box: Box<3> bbox (Box<3>::EMPTY_BOX); - for (int i = 0; i < adfront->GetNFL(); i++) + for (int i = 0; i < adfront.GetNFL(); i++) { - const FrontLine & line = adfront->GetLine(i); - bbox.Add (Point<3> (adfront->GetPoint (line.L()[0]))); - bbox.Add (Point<3> (adfront->GetPoint (line.L()[1]))); + const FrontLine & line = adfront.GetLine(i); + bbox.Add (Point<3> (adfront.GetPoint (line.L()[0]))); + bbox.Add (Point<3> (adfront.GetPoint (line.L()[1]))); } for (int i = 0; i < mesh.LockedPoints().Size(); i++) @@ -402,7 +402,7 @@ namespace netgen if (trig[0] < 0) continue; Point<3> c = Center (mesh[trig[0]], mesh[trig[1]], mesh[trig[2]]); - if (!adfront->Inside (Point<2> (c(0),c(1)))) continue; + if (!adfront.Inside (Point<2> (c(0),c(1)))) continue; Vec<3> n = Cross (mesh[trig[1]]-mesh[trig[0]], mesh[trig[2]]-mesh[trig[0]]); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 695ea807..8d95f222 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -20,7 +20,7 @@ namespace netgen segmentht = NULL; lochfunc = NULL; - mglevels = 1; + // mglevels = 1; elementsearchtree = NULL; elementsearchtreets = NextTimeStamp(); majortimestamp = timestamp = NextTimeStamp(); diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 78ea8035..09d17753 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -175,7 +175,9 @@ namespace netgen /// number of refinement levels - int mglevels; + // int mglevels; + // number of vertices on each refinement level: + NgArray level_nv; /// refinement hierarchy NgArray,PointIndex::BASE> mlbetweennodes; /// parent element of volume element @@ -795,8 +797,9 @@ namespace netgen shared_ptr GetGeometry() const - { - return geometry; + { + static auto global_geometry = make_shared(); + return geometry ? geometry : global_geometry; } void SetGeometry (shared_ptr geom) { diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index b3630df6..b7d91a14 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -690,7 +690,7 @@ namespace netgen case 'j': mesh3d.ImproveMeshJacobian(mp); break; } } - mesh3d.mglevels = 1; + // mesh3d.mglevels = 1; MeshQuality3d (mesh3d); } diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp index a5b0afef..9c794da2 100644 --- a/libsrc/meshing/meshfunc2d.cpp +++ b/libsrc/meshing/meshfunc2d.cpp @@ -63,10 +63,7 @@ namespace netgen } if (secondorder) { - if (mesh.GetGeometry()) - mesh.GetGeometry()->GetRefinement().MakeSecondOrder(mesh); - else - Refinement().MakeSecondOrder(mesh); + mesh.GetGeometry()->GetRefinement().MakeSecondOrder(mesh); } } diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp index 91b008e3..685d4865 100644 --- a/libsrc/meshing/meshing2.cpp +++ b/libsrc/meshing/meshing2.cpp @@ -19,15 +19,33 @@ namespace netgen // static int qualclass; + static Array> global_trig_rules; + static Array> global_quad_rules; + + Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox) + : adfront(aboundingbox), boundingbox(aboundingbox) { - boundingbox = aboundingbox; - - LoadRules (NULL, mp.quad); + static Timer t("Mesing2::Meshing2"); RegionTimer r(t); + + auto & globalrules = mp.quad ? global_quad_rules : global_trig_rules; + if (!globalrules.Size()) + { + LoadRules (NULL, mp.quad); + for (auto * rule : rules) + globalrules.Append (unique_ptr(rule)); + } + else + { + for (auto i : globalrules.Range()) + rules.Append (globalrules[i].get()); + } // LoadRules ("rules/quad.rls"); // LoadRules ("rules/triangle.rls"); - adfront = new AdFront2(boundingbox); + + + // adfront = new AdFront2(boundingbox); starttime = GetTime(); maxarea = -1; @@ -35,18 +53,14 @@ namespace netgen Meshing2 :: ~Meshing2 () - { - delete adfront; - for (int i = 0; i < rules.Size(); i++) - delete rules[i]; - } + { ; } void Meshing2 :: AddPoint (const Point3d & p, PointIndex globind, MultiPointGeomInfo * mgi, bool pointonsurface) { //(*testout) << "add point " << globind << endl; - adfront ->AddPoint (p, globind, mgi, pointonsurface); + adfront.AddPoint (p, globind, mgi, pointonsurface); } void Meshing2 :: AddBoundaryElement (int i1, int i2, @@ -57,7 +71,7 @@ namespace netgen { PrintSysError ("addboundaryelement: illegal geominfo"); } - adfront -> AddLine (i1-1, i2-1, gi1, gi2); + adfront.AddLine (i1-1, i2-1, gi1, gi2); } @@ -219,7 +233,6 @@ namespace netgen int z1, z2, oldnp(-1); bool found; int rulenr(-1); - Point<3> p1, p2; const PointGeomInfo * blgeominfo1; const PointGeomInfo * blgeominfo2; @@ -227,9 +240,8 @@ namespace netgen bool morerisc; bool debugflag; - double h, his, hshould; - - + // double h; + NgArray locpoints; NgArray legalpoints; NgArray plainpoints; @@ -340,7 +352,7 @@ namespace netgen const char * savetask = multithread.task; multithread.task = "Surface meshing"; - adfront ->SetStartFront (); + adfront.SetStartFront (); int plotnexttrial = 999; @@ -349,7 +361,11 @@ namespace netgen NgProfiler::StopTimer (ts3); - while (!adfront ->Empty() && !multithread.terminate) + static Timer tloop("surfacemeshing mainloop"); + static Timer tgetlocals("surfacemeshing getlocals"); + { + RegionTimer rloop(tloop); + while (!adfront.Empty() && !multithread.terminate) { NgProfiler::RegionTimer reg1 (timer1); @@ -364,13 +380,13 @@ namespace netgen multithread.percent = 0; */ - locpoints.SetSize(0); - loclines.SetSize(0); - pindex.SetSize(0); - lindex.SetSize(0); - delpoints.SetSize(0); - dellines.SetSize(0); - locelements.SetSize(0); + locpoints.SetSize0(); + loclines.SetSize0(); + pindex.SetSize0(); + lindex.SetSize0(); + delpoints.SetSize0(); + dellines.SetSize0(); + locelements.SetSize0(); @@ -388,11 +404,11 @@ namespace netgen // unique-pgi, multi-pgi - upgeominfo.SetSize(0); - mpgeominfo.SetSize(0); + upgeominfo.SetSize0(); + mpgeominfo.SetSize0(); - nfaces = adfront->GetNFL(); + nfaces = adfront.GetNFL(); trials ++; @@ -408,27 +424,27 @@ namespace netgen (*testout) << "\n"; } - - int baselineindex = adfront -> SelectBaseLine (p1, p2, blgeominfo1, blgeominfo2, qualclass); - + Point<3> p1, p2; + int baselineindex = adfront.SelectBaseLine (p1, p2, blgeominfo1, blgeominfo2, qualclass); found = 1; - his = Dist (p1, p2); + double his = Dist (p1, p2); - Point3d pmid = Center (p1, p2); - hshould = CalcLocalH (pmid, mesh.GetH (pmid)); + Point<3> pmid = Center (p1, p2); + double hshould = CalcLocalH (pmid, mesh.GetH (pmid)); if (gh < hshould) hshould = gh; mesh.RestrictLocalH (pmid, hshould); - h = hshould; + double h = hshould; double hinner = (3 + qualclass) * max2 (his, hshould); - adfront ->GetLocals (baselineindex, locpoints, mpgeominfo, loclines, + tgetlocals.Start(); + adfront.GetLocals (baselineindex, locpoints, mpgeominfo, loclines, pindex, lindex, 2*hinner); - + tgetlocals.Stop(); NgProfiler::RegionTimer reg2 (timer2); @@ -440,7 +456,7 @@ namespace netgen if (qualclass > mp.giveuptol2d) { PrintMessage (3, "give up with qualclass ", qualclass); - PrintMessage (3, "number of frontlines = ", adfront->GetNFL()); + PrintMessage (3, "number of frontlines = ", adfront.GetNFL()); // throw NgException ("Give up 2d meshing"); break; } @@ -456,8 +472,8 @@ namespace netgen morerisc = 0; - PointIndex gpi1 = adfront -> GetGlobalIndex (pindex.Get(loclines[0].I1())); - PointIndex gpi2 = adfront -> GetGlobalIndex (pindex.Get(loclines[0].I2())); + PointIndex gpi1 = adfront.GetGlobalIndex (pindex.Get(loclines[0].I1())); + PointIndex gpi2 = adfront.GetGlobalIndex (pindex.Get(loclines[0].I2())); debugflag = @@ -515,6 +531,12 @@ namespace netgen *testout << "3d points: " << endl << locpoints << endl; } + + for (size_t i = 0; i < locpoints.Size(); i++) + TransformToPlain (locpoints[i], mpgeominfo[i], + plainpoints[i], h, plainzones[i]); + + /* for (int i = 1; i <= locpoints.Size(); i++) { // (*testout) << "pindex(i) = " << pindex[i-1] << endl; @@ -525,6 +547,7 @@ namespace netgen // (*testout) << plainpoints.Get(i).X() << " " << plainpoints.Get(i).Y() << endl; //(*testout) << "transform " << locpoints.Get(i) << " to " << plainpoints.Get(i).X() << " " << plainpoints.Get(i).Y() << endl; } + */ // (*testout) << endl << endl << endl; @@ -579,7 +602,7 @@ namespace netgen if (IsLineVertexOnChart (locpoints.Get(loclines.Get(i).I1()), locpoints.Get(loclines.Get(i).I2()), innerp, - adfront->GetLineGeomInfo (lindex.Get(i), innerp))) + adfront.GetLineGeomInfo (lindex.Get(i), innerp))) // pgeominfo.Get(loclines.Get(i).I(innerp)))) { @@ -651,31 +674,34 @@ namespace netgen legalpoints.SetSize(plainpoints.Size()); + legalpoints = 1; + /* for (int i = 1; i <= legalpoints.Size(); i++) legalpoints.Elem(i) = 1; - + */ + double avy = 0; - for (int i = 1; i <= plainpoints.Size(); i++) - avy += plainpoints.Elem(i).Y(); + for (size_t i = 0; i < plainpoints.Size(); i++) + avy += plainpoints[i].Y(); avy *= 1./plainpoints.Size(); - for (int i = 1; i <= plainpoints.Size(); i++) + for (auto i : Range(plainpoints)) { - if (plainzones.Elem(i) < 0) + if (plainzones[i] < 0) { - plainpoints.Elem(i) = Point2d (1e4, 1e4); - legalpoints.Elem(i) = 0; + plainpoints[i] = Point2d (1e4, 1e4); + legalpoints[i] = 0; } - if (pindex.Elem(i) == -1) + if (pindex[i] == -1) { - legalpoints.Elem(i) = 0; + legalpoints[i] = 0; } - if (plainpoints.Elem(i).Y() < -1e-10*avy) // changed + if (plainpoints[i].Y() < -1e-10*avy) // changed { - legalpoints.Elem(i) = 0; + legalpoints[i] = 0; } } /* @@ -758,12 +784,14 @@ namespace netgen { multithread.drawing = 1; glrender(1); - cout << "qualclass 100, nfl = " << adfront->GetNFL() << endl; + cout << "qualclass 100, nfl = " << adfront.GetNFL() << endl; } */ if (found) { + static Timer t("ApplyRules"); + RegionTimer r(t); rulenr = ApplyRules (plainpoints, legalpoints, maxlegalpoint, loclines, maxlegalline, locelements, dellines, qualclass, mp); @@ -818,7 +846,7 @@ namespace netgen // for (i = 1; i <= oldnl; i++) - // adfront -> ResetClass (lindex[i]); + // adfront.ResetClass (lindex[i]); /* @@ -947,7 +975,7 @@ namespace netgen for (j = 1; j <= 2; j++) { upgeominfo.Elem(loclines.Get(dellines.Get(i)).I(j)) = - adfront -> GetLineGeomInfo (lindex.Get(dellines.Get(i)), j); + adfront.GetLineGeomInfo (lindex.Get(dellines.Get(i)), j); } */ @@ -1145,7 +1173,7 @@ namespace netgen // cout << "overlap !!!" << endl; #endif for (int k = 1; k <= 5; k++) - adfront -> IncrementClass (lindex.Get(1)); + adfront.IncrementClass (lindex.Get(1)); found = 0; @@ -1179,10 +1207,10 @@ namespace netgen int nlgpi2 = loclines.Get(i).I2(); if (nlgpi1 <= pindex.Size() && nlgpi2 <= pindex.Size()) { - nlgpi1 = adfront->GetGlobalIndex (pindex.Get(nlgpi1)); - nlgpi2 = adfront->GetGlobalIndex (pindex.Get(nlgpi2)); + nlgpi1 = adfront.GetGlobalIndex (pindex.Get(nlgpi1)); + nlgpi2 = adfront.GetGlobalIndex (pindex.Get(nlgpi2)); - int exval = adfront->ExistsLine (nlgpi1, nlgpi2); + int exval = adfront.ExistsLine (nlgpi1, nlgpi2); if (exval) { cout << "ERROR: new line exits, val = " << exval << endl; @@ -1211,8 +1239,8 @@ namespace netgen int tpi2 = locelements.Get(i).PNumMod (j+1); if (tpi1 <= pindex.Size() && tpi2 <= pindex.Size()) { - tpi1 = adfront->GetGlobalIndex (pindex.Get(tpi1)); - tpi2 = adfront->GetGlobalIndex (pindex.Get(tpi2)); + tpi1 = adfront.GetGlobalIndex (pindex.Get(tpi1)); + tpi2 = adfront.GetGlobalIndex (pindex.Get(tpi2)); if (doubleedge.Used (INDEX_2(tpi1, tpi2))) { @@ -1241,7 +1269,7 @@ namespace netgen for (int i = oldnp+1; i <= locpoints.Size(); i++) { PointIndex globind = mesh.AddPoint (locpoints.Get(i)); - pindex.Elem(i) = adfront -> AddPoint (locpoints.Get(i), globind); + pindex.Elem(i) = adfront.AddPoint (locpoints.Get(i), globind); } for (int i = oldnl+1; i <= loclines.Size(); i++) @@ -1271,7 +1299,7 @@ namespace netgen cout << "new el: illegal geominfo" << endl; } - adfront -> AddLine (pindex.Get(loclines.Get(i).I1()), + adfront.AddLine (pindex.Get(loclines.Get(i).I1()), pindex.Get(loclines.Get(i).I2()), upgeominfo.Get(loclines.Get(i).I1()), upgeominfo.Get(loclines.Get(i).I2())); @@ -1296,7 +1324,7 @@ namespace netgen { mtri.PNum(j) = locelements.Elem(i).PNum(j) = - adfront -> GetGlobalIndex (pindex.Get(locelements.Get(i).PNum(j))); + adfront.GetGlobalIndex (pindex.Get(locelements.Get(i).PNum(j))); } @@ -1375,7 +1403,7 @@ namespace netgen } for (int i = 1; i <= dellines.Size(); i++) - adfront -> DeleteLine (lindex.Get(dellines.Get(i))); + adfront.DeleteLine (lindex.Get(dellines.Get(i))); // rname = rules.Get(rulenr)->Name(); #ifdef MYGRAPH @@ -1398,7 +1426,7 @@ namespace netgen if ( debugparam.haltsuccess || debugflag ) { - // adfront -> PrintOpenSegments (*testout); + // adfront.PrintOpenSegments (*testout); cout << "success of rule" << rules.Get(rulenr)->Name() << endl; multithread.drawing = 1; multithread.testmode = 1; @@ -1420,7 +1448,7 @@ namespace netgen (*testout) << "locpoints " << endl; for (int i = 1; i <= pindex.Size(); i++) - (*testout) << adfront->GetGlobalIndex (pindex.Get(i)) << endl; + (*testout) << adfront.GetGlobalIndex (pindex.Get(i)) << endl; (*testout) << "old number of lines = " << oldnl << endl; for (int i = 1; i <= loclines.Size(); i++) @@ -1431,7 +1459,7 @@ namespace netgen int hi = 0; if (loclines.Get(i).I(j) >= 1 && loclines.Get(i).I(j) <= pindex.Size()) - hi = adfront->GetGlobalIndex (pindex.Get(loclines.Get(i).I(j))); + hi = adfront.GetGlobalIndex (pindex.Get(loclines.Get(i).I(j))); (*testout) << hi << " "; } @@ -1450,7 +1478,7 @@ namespace netgen } else { - adfront -> IncrementClass (lindex.Get(1)); + adfront.IncrementClass (lindex.Get(1)); if ( debugparam.haltnosuccess || debugflag ) { @@ -1483,7 +1511,7 @@ namespace netgen int hi = 0; if (loclines.Get(i).I(j) >= 1 && loclines.Get(i).I(j) <= pindex.Size()) - hi = adfront->GetGlobalIndex (pindex.Get(loclines.Get(i).I(j))); + hi = adfront.GetGlobalIndex (pindex.Get(loclines.Get(i).I(j))); (*testout) << hi << " "; } @@ -1518,11 +1546,11 @@ namespace netgen } } - + } PrintMessage (3, "Surface meshing done"); - adfront->PrintOpenSegments (*testout); + adfront.PrintOpenSegments (*testout); multithread.task = savetask; @@ -1530,7 +1558,7 @@ namespace netgen EndMesh (); - if (!adfront->Empty()) + if (!adfront.Empty()) return MESHING2_GIVEUP; return MESHING2_OK; diff --git a/libsrc/meshing/meshing2.hpp b/libsrc/meshing/meshing2.hpp index a23e06d2..47fd5b01 100644 --- a/libsrc/meshing/meshing2.hpp +++ b/libsrc/meshing/meshing2.hpp @@ -29,7 +29,7 @@ derive from Meshing2, and replace transformation. class Meshing2 { /// the current advancing front - AdFront2 * adfront; + AdFront2 adfront; /// rules for mesh generation NgArray rules; /// statistics diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 68fee518..bbc5b527 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1191,7 +1191,7 @@ namespace netgen /// power of error (to approximate max err optimization) double opterrpow = 2; /// do block filling ? - int blockfill = 1; + bool blockfill = true; /// block filling up to distance double filldist = 0.1; /// radius of local environment (times h) @@ -1199,11 +1199,11 @@ namespace netgen /// radius of active environment (times h) double relinnersafety = 3; /// use local h ? - int uselocalh = 1; + bool uselocalh = true; /// grading for local h double grading = 0.3; /// use delaunay meshing - int delaunay = 1; + bool delaunay = true; /// maximal mesh size double maxh = 1e10; /// minimal mesh size @@ -1211,19 +1211,19 @@ namespace netgen /// file for meshsize string meshsizefilename = ""; /// start surfacemeshing from everywhere in surface - int startinsurface = 0; + bool startinsurface = false; /// check overlapping surfaces (debug) - int checkoverlap = 1; + bool checkoverlap = true; /// check overlapping surface mesh before volume meshing - int checkoverlappingboundary = 1; + bool checkoverlappingboundary = true; /// check chart boundary (sometimes too restrictive) - int checkchartboundary = 1; + bool checkchartboundary = true; /// safety factor for curvatures (elements per radius) double curvaturesafety = 2; /// minimal number of segments per edge double segmentsperedge = 1; /// use parallel threads - int parthread = 0; + bool parthread = 0; /// weight of element size w.r.t element shape double elsizeweight = 0.2; /// init with default values @@ -1246,29 +1246,29 @@ namespace netgen /// if non-zero, baseelement must have baseelnp points int baseelnp = 0; /// quality tolerances are handled less careful - int sloppy = 1; + bool sloppy = true; /// limit for max element angle (150-180) double badellimit = 175; - bool check_impossible = 0; + bool check_impossible = false; int only3D_domain_nr = 0; /// - int secondorder = 0; + bool secondorder = false; /// high order element curvature int elementorder = 1; /// quad-dominated surface meshing - int quad = 0; + bool quad = false; /// bool try_hexes = false; /// - int inverttets = 0; + bool inverttets = false; /// - int inverttrigs = 0; + bool inverttrigs = false; /// - int autozrefine = 0; + bool autozrefine = false; /// MeshingParameters (); /// diff --git a/libsrc/meshing/parser2.cpp b/libsrc/meshing/parser2.cpp index 3f4892e6..b0fb1bc7 100644 --- a/libsrc/meshing/parser2.cpp +++ b/libsrc/meshing/parser2.cpp @@ -578,7 +578,9 @@ void Meshing2 :: LoadRules (const char * filename, bool quad) delete ist; exit (1); } - + + Timer t("Parsing rules"); + t.Start(); while (!ist->eof()) { buf[0] = 0; @@ -597,7 +599,8 @@ void Meshing2 :: LoadRules (const char * filename, bool quad) //(*testout) << "loop" << endl; } //(*testout) << "POS3" << endl; - + t.Stop(); + delete ist; //delete [] tr1; } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 448775d8..a58faa70 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -2,6 +2,7 @@ #include <../general/ngpython.hpp> #include +#include "python_mesh.hpp" #include #include "meshing.hpp" @@ -856,24 +857,20 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) self.SetMaxHDomain(maxh); }) .def ("GenerateVolumeMesh", - [](Mesh & self, py::object pymp) + [](Mesh & self, MeshingParameters* pars, + py::kwargs kwargs) { - cout << "generate vol mesh" << endl; - MeshingParameters mp; + if(pars) mp = *pars; { py::gil_scoped_acquire acquire; - if (py::extract(pymp).check()) - mp = py::extract(pymp)(); - else - { - mp.optsteps3d = 5; - } + CreateMPfromKwargs(mp, kwargs); } MeshVolume (mp, self); OptimizeVolume (mp, self); - }, - py::arg("mp")=NGDummyArgument(),py::call_guard()) + }, py::arg("mp")=nullptr, + meshingparameter_description.c_str(), + py::call_guard()) .def ("OptimizeVolumeMesh", [](Mesh & self) { @@ -893,20 +890,14 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def ("Refine", FunctionPointer ([](Mesh & self) { - if (self.GetGeometry()) - self.GetGeometry()->GetRefinement().Refine(self); - else - Refinement().Refine(self); + self.GetGeometry()->GetRefinement().Refine(self); self.UpdateTopology(); }),py::call_guard()) .def ("SecondOrder", FunctionPointer ([](Mesh & self) { - if (self.GetGeometry()) - self.GetGeometry()->GetRefinement().MakeSecondOrder(self); - else - Refinement().MakeSecondOrder(self); + self.GetGeometry()->GetRefinement().MakeSecondOrder(self); })) .def ("GetGeometry", [] (Mesh& self) { return self.GetGeometry(); }) @@ -1026,42 +1017,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) ; typedef MeshingParameters MP; - py::class_ (m, "MeshingParameters") + auto mp = py::class_ (m, "MeshingParameters") .def(py::init<>()) - .def(py::init([](double maxh, bool quad_dominated, int optsteps2d, int optsteps3d, - MESHING_STEP perfstepsend, int only3D_domain, const string & meshsizefilename, - double grading, double curvaturesafety, double segmentsperedge) + .def(py::init([](py::kwargs kwargs) { - MP * instance = new MeshingParameters; - instance->maxh = maxh; - instance->quad = int(quad_dominated); - instance->optsteps2d = optsteps2d; - instance->optsteps3d = optsteps3d; - instance->only3D_domain_nr = only3D_domain; - instance->perfstepsend = perfstepsend; - instance->meshsizefilename = meshsizefilename; - - instance->grading = grading; - instance->curvaturesafety = curvaturesafety; - instance->segmentsperedge = segmentsperedge; - return instance; - }), - py::arg("maxh")=1000, - py::arg("quad_dominated")=false, - py::arg("optsteps2d") = 3, - py::arg("optsteps3d") = 3, - py::arg("perfstepsend") = MESHCONST_OPTVOLUME, - py::arg("only3D_domain") = 0, - py::arg("meshsizefilename") = "", - py::arg("grading")=0.3, - py::arg("curvaturesafety")=2, - py::arg("segmentsperedge")=1, - "create meshing parameters" - ) + MeshingParameters mp; + CreateMPfromKwargs(mp, kwargs); + return mp; + }), meshingparameter_description.c_str()) .def("__str__", &ToString) - .def_property("maxh", - FunctionPointer ([](const MP & mp ) { return mp.maxh; }), - FunctionPointer ([](MP & mp, double maxh) { return mp.maxh = maxh; })) + .def_property("maxh", [](const MP & mp ) { return mp.maxh; }, + [](MP & mp, double maxh) { return mp.maxh = maxh; }) + .def_property("grading", [](const MP & mp ) { return mp.grading; }, + [](MP & mp, double grading) { return mp.grading = grading; }) .def("RestrictH", FunctionPointer ([](MP & mp, double x, double y, double z, double h) { diff --git a/libsrc/meshing/python_mesh.hpp b/libsrc/meshing/python_mesh.hpp new file mode 100644 index 00000000..bbe75a97 --- /dev/null +++ b/libsrc/meshing/python_mesh.hpp @@ -0,0 +1,168 @@ + +#include +#include "meshing.hpp" + +namespace netgen +{ + // TODO: Clarify a lot of these parameters + static string meshingparameter_description = R"delimiter( +Meshing Parameters +------------------- + +maxh: float = 1e10 + Global upper bound for mesh size. + +grading: float = 0.3 + Mesh grading how fast the local mesh size can change. + +meshsizefilename: str = None + Load meshsize from file. Can set local mesh size for points + and along edges. File must have the format: + + nr_points + x1, y1, z1, meshsize + x2, y2, z2, meshsize + ... + xn, yn, zn, meshsize + + nr_edges + x11, y11, z11, x12, y12, z12, meshsize + ... + xn1, yn1, zn1, xn2, yn2, zn2, meshsize + +segmentsperedge: float = 1. + Minimal number of segments per edge. + +quad: bool = False + Quad-dominated surface meshing. + +blockfill: bool = True + Do fast blockfilling. + +filldist: float = 0.1 + Block fill up to distance + +delaunay: bool = True + Use delaunay meshing. + +Optimization Parameters +----------------------- + +optimize3d: str = "cmdmustm" + 3d optimization strategy: + m .. move nodes + M .. move nodes, cheap functional + s .. swap faces + c .. combine elements + d .. divide elements + p .. plot, no pause + P .. plot, Pause + h .. Histogramm, no pause + H .. Histogramm, pause + +optsteps3d: int = 3 + Number of 3d optimization steps. + +optimize2d: str = "smsmsmSmSmSm" + 2d optimization strategy: + s .. swap, opt 6 lines/node + S .. swap, optimal elements + m .. move nodes + p .. plot, no pause + P .. plot, pause + c .. combine + +optsteps2d: int = 3 + Number of 2d optimization steps. + +elsizeweight: float = 0.2 + Weight of element size w.r.t. element shape in optimization. + +)delimiter"; + + inline void CreateMPfromKwargs(MeshingParameters& mp, py::kwargs kwargs) + { + if(kwargs.contains("optimize3d")) + mp.optimize3d = py::cast(kwargs["optimize3d"]); + if(kwargs.contains("optsteps3d")) + mp.optsteps3d = py::cast(kwargs["optsteps3d"]); + if(kwargs.contains("optimize2d")) + mp.optimize2d = py::cast(kwargs["optimize2d"]); + if(kwargs.contains("optsteps2d")) + mp.optsteps2d = py::cast(kwargs["optsteps2d"]); + if(kwargs.contains("opterrpow")) + mp.opterrpow = py::cast(kwargs["opterrpow"]); + if(kwargs.contains("blockfill")) + mp.blockfill = py::cast(kwargs["blockfill"]); + if(kwargs.contains("filldist")) + mp.filldist = py::cast(kwargs["filldist"]); + if(kwargs.contains("safety")) + mp.safety = py::cast(kwargs["safety"]); + if(kwargs.contains("relinnersafety")) + mp.relinnersafety = py::cast(kwargs["relinnersafety"]); + if(kwargs.contains("uselocalh")) + mp.uselocalh = py::cast(kwargs["uselocalh"]); + if(kwargs.contains("grading")) + mp.grading = py::cast(kwargs["grading"]); + if(kwargs.contains("delaunay")) + mp.delaunay = py::cast(kwargs["delaunay"]); + if(kwargs.contains("maxh")) + mp.maxh = py::cast(kwargs["maxh"]); + if(kwargs.contains("minh")) + mp.minh = py::cast(kwargs["minh"]); + if(kwargs.contains("meshsizefilename")) + mp.meshsizefilename = py::cast(kwargs["meshsizefilename"]); + if(kwargs.contains("startinsurface")) + mp.startinsurface = py::cast(kwargs["startinsurface"]); + if(kwargs.contains("checkoverlap")) + mp.checkoverlap = py::cast(kwargs["checkoverlap"]); + if(kwargs.contains("checkoverlappingboundary")) + mp.checkoverlappingboundary = py::cast(kwargs["checkoverlappingboundary"]); + if(kwargs.contains("checkchartboundary")) + mp.checkchartboundary = py::cast(kwargs["checkchartboundary"]); + if(kwargs.contains("curvaturesafety")) + mp.curvaturesafety = py::cast(kwargs["curvaturesafety"]); + if(kwargs.contains("segmentsperedge")) + mp.segmentsperedge = py::cast(kwargs["segmentsperedge"]); + if(kwargs.contains("parthread")) + mp.parthread = py::cast(kwargs["parthread"]); + if(kwargs.contains("elsizeweight")) + mp.elsizeweight = py::cast(kwargs["elsizeweight"]); + if(kwargs.contains("perfstepsstart")) + mp.perfstepsstart = py::cast(kwargs["perfstepsstart"]); + if(kwargs.contains("perfstepsend")) + mp.perfstepsend = py::cast(kwargs["perfstepsend"]); + if(kwargs.contains("giveuptol2d")) + mp.giveuptol2d = py::cast(kwargs["giveuptol2d"]); + if(kwargs.contains("giveuptol")) + mp.giveuptol = py::cast(kwargs["giveuptol"]); + if(kwargs.contains("maxoutersteps")) + mp.maxoutersteps = py::cast(kwargs["maxoutersteps"]); + if(kwargs.contains("starshapeclass")) + mp.starshapeclass = py::cast(kwargs["starshapeclass"]); + if(kwargs.contains("baseelnp")) + mp.baseelnp = py::cast(kwargs["baseelnp"]); + if(kwargs.contains("sloppy")) + mp.sloppy = py::cast(kwargs["sloppy"]); + if(kwargs.contains("badellimit")) + mp.badellimit = py::cast(kwargs["badellimit"]); + if(kwargs.contains("check_impossible")) + mp.check_impossible = py::cast(kwargs["check_impossible"]); + if(kwargs.contains("only3D_domain_nr")) + mp.only3D_domain_nr = py::cast(kwargs["only3D_domain_nr"]); + if(kwargs.contains("secondorder")) + mp.secondorder = py::cast(kwargs["secondorder"]); + if(kwargs.contains("elementorder")) + mp.elementorder = py::cast(kwargs["elementorder"]); + if(kwargs.contains("quad")) + mp.quad = py::cast(kwargs["quad"]); + if(kwargs.contains("try_hexes")) + mp.try_hexes = py::cast(kwargs["try_hexes"]); + if(kwargs.contains("inverttets")) + mp.inverttets = py::cast(kwargs["inverttets"]); + if(kwargs.contains("inverttrigs")) + mp.inverttrigs = py::cast(kwargs["inverttrigs"]); + if(kwargs.contains("autozrefine")) + mp.autozrefine = py::cast(kwargs["autozrefine"]); + } +} // namespace netgen diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index fd814fec..18e7077f 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -1357,8 +1357,6 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) { static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t); - int typ = 1; - (*testout) << "Improve Mesh" << "\n"; PrintMessage (3, "ImproveMesh"); @@ -1366,36 +1364,9 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) int ne = GetNE(); - NgArray perrs(np); - perrs = 1.0; - - double bad1 = 0; - double badmax = 0; - if (goal == OPT_QUALITY) { - for (int i = 1; i <= ne; i++) - { - const Element & el = VolumeElement(i); - if (el.GetType() != TET) - continue; - - double hbad = CalcBad (points, el, 0, mp); - for (int j = 0; j < 4; j++) - perrs[el[j]] += hbad; - - bad1 += hbad; - } - - for (int i = perrs.Begin(); i < perrs.End(); i++) - if (perrs[i] > badmax) - badmax = perrs[i]; - badmax = 0; - } - - if (goal == OPT_QUALITY) - { - bad1 = CalcTotalBad (points, volelements, mp); + double bad1 = CalcTotalBad (points, volelements, mp); (*testout) << "Total badness = " << bad1 << endl; PrintMessage (5, "Total badness = ", bad1); } @@ -1407,16 +1378,9 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) //int uselocalh = mparam.uselocalh; - PointFunction * pf; - - if (typ == 1) - pf = new PointFunction(points, volelements, mp); - else - pf = new CheapPointFunction(points, volelements, mp); - - // pf->SetLocalH (h); + PointFunction pf(points, volelements, mp); - Opti3FreeMinFunction freeminf(*pf); + Opti3FreeMinFunction freeminf(pf); OptiParameters par; par.maxit_linsearch = 20; @@ -1460,7 +1424,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) multithread.task = "Smooth Mesh"; for (PointIndex pi : points.Range()) - if ( (*this)[pi].Type() == INNERPOINT && perrs[pi] > 0.01 * badmax) + if ( (*this)[pi].Type() == INNERPOINT ) { if (multithread.terminate) throw NgException ("Meshing stopped"); @@ -1470,11 +1434,11 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) if ( (pi+1-PointIndex::BASE) % printmod == 0) PrintDot (printdot); double lh = pointh[pi]; - pf->SetLocalH (lh); + pf.SetLocalH (lh); par.typx = lh; freeminf.SetPoint (points[pi]); - pf->SetPointIndex (pi); + pf.SetPointIndex (pi); x = 0; int pok; @@ -1482,10 +1446,10 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) if (!pok) { - pok = pf->MovePointToInner (); + pok = pf.MovePointToInner (); freeminf.SetPoint (points[pi]); - pf->SetPointIndex (pi); + pf.SetPointIndex (pi); } if (pok) @@ -1500,13 +1464,11 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) } PrintDot ('\n'); - delete pf; - multithread.task = savetask; if (goal == OPT_QUALITY) { - bad1 = CalcTotalBad (points, volelements, mp); + double bad1 = CalcTotalBad (points, volelements, mp); (*testout) << "Total badness = " << bad1 << endl; PrintMessage (5, "Total badness = ", bad1); } diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index 74d0328c..d80f35e0 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -18,712 +18,722 @@ namespace netgen #define VSMALL 1e-10 - DLL_HEADER bool merge_solids = 1; + DLL_HEADER bool merge_solids = 1; // can you please explain what you intend to compute here (JS) !!! - double Line :: Dist (Line l) - { - Vec<3> n = p1-p0; - Vec<3> q = l.p1-l.p0; - double nq = n*q; + double Line :: Dist (Line l) + { + Vec<3> n = p1-p0; + Vec<3> q = l.p1-l.p0; + double nq = n*q; - Point<3> p = p0 + 0.5*n; - double lambda = (p-l.p0)*n / (nq + VSMALL); + Point<3> p = p0 + 0.5*n; + double lambda = (p-l.p0)*n / (nq + VSMALL); - if (lambda >= 0 && lambda <= 1) + if (lambda >= 0 && lambda <= 1) { - double d = (p-l.p0-lambda*q).Length(); - // if (d < 1e-3) d = 1e99; - return d; + double d = (p-l.p0-lambda*q).Length(); + // if (d < 1e-3) d = 1e99; + return d; } - else - return 1e99; - } + else + return 1e99; + } + + inline Point<3> occ2ng (const gp_Pnt & p) + { + return Point<3> (p.X(), p.Y(), p.Z()); + } + + double ComputeH (double kappa, const MeshingParameters & mparam) + { + /* + double hret; + kappa *= mparam.curvaturesafety; + + if (mparam.maxh * kappa < 1) + hret = mparam.maxh; + else + hret = 1 / (kappa + VSMALL); + + if (mparam.maxh < hret) + hret = mparam.maxh; + + return hret; + */ + // return min(mparam.maxh, 1/kappa); + return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa; + } - double Line :: Length () - { - return (p1-p0).Length(); - } + void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, + BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h, + const MeshingParameters & mparam) + { + int ls = -1; + gp_Pnt pnt0,pnt1,pnt2; + prop->SetParameters (par0.X(), par0.Y()); + pnt0 = prop->Value(); - inline Point<3> occ2ng (const gp_Pnt & p) - { - return Point<3> (p.X(), p.Y(), p.Z()); - } + prop->SetParameters (par1.X(), par1.Y()); + pnt1 = prop->Value(); + prop->SetParameters (par2.X(), par2.Y()); + pnt2 = prop->Value(); - - double ComputeH (double kappa) - { - double hret; - kappa *= mparam.curvaturesafety; - - if (mparam.maxh * kappa < 1) - hret = mparam.maxh; - else - hret = 1 / (kappa + VSMALL); - - if (mparam.maxh < hret) - hret = mparam.maxh; - - return (hret); - } - - - - - void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, - BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0) - { - int ls = -1; - - gp_Pnt pnt0,pnt1,pnt2; - - prop->SetParameters (par0.X(), par0.Y()); - pnt0 = prop->Value(); - - prop->SetParameters (par1.X(), par1.Y()); - pnt1 = prop->Value(); - - prop->SetParameters (par2.X(), par2.Y()); - pnt2 = prop->Value(); - - double aux; - double maxside = pnt0.Distance(pnt1); - ls = 2; - aux = pnt1.Distance(pnt2); - if(aux > maxside) + double aux; + double maxside = pnt0.Distance(pnt1); + ls = 2; + aux = pnt1.Distance(pnt2); + if(aux > maxside) { - maxside = aux; - ls = 0; + maxside = aux; + ls = 0; } - aux = pnt2.Distance(pnt0); - if(aux > maxside) + aux = pnt2.Distance(pnt0); + if(aux > maxside) { - maxside = aux; - ls = 1; + maxside = aux; + ls = 1; } - gp_Pnt2d parmid; + gp_Pnt2d parmid; - parmid.SetX( (par0.X()+par1.X()+par2.X()) / 3 ); - parmid.SetY( (par0.Y()+par1.Y()+par2.Y()) / 3 ); + parmid.SetX( (par0.X()+par1.X()+par2.X()) / 3 ); + parmid.SetY( (par0.Y()+par1.Y()+par2.Y()) / 3 ); - if (depth%3 == 0) + if (depth%3 == 0) { - double curvature = 0; + double curvature = 0; - prop->SetParameters (parmid.X(), parmid.Y()); - if (!prop->IsCurvatureDefined()) - { + prop->SetParameters (parmid.X(), parmid.Y()); + if (!prop->IsCurvatureDefined()) + { (*testout) << "curvature not defined!" << endl; return; - } - curvature = max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature())); + } + curvature = max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature())); - prop->SetParameters (par0.X(), par0.Y()); - if (!prop->IsCurvatureDefined()) - { + prop->SetParameters (par0.X(), par0.Y()); + if (!prop->IsCurvatureDefined()) + { (*testout) << "curvature not defined!" << endl; return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); + } + curvature = max(curvature,max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature()))); - prop->SetParameters (par1.X(), par1.Y()); - if (!prop->IsCurvatureDefined()) - { + prop->SetParameters (par1.X(), par1.Y()); + if (!prop->IsCurvatureDefined()) + { (*testout) << "curvature not defined!" << endl; return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); + } + curvature = max(curvature,max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature()))); - prop->SetParameters (par2.X(), par2.Y()); - if (!prop->IsCurvatureDefined()) - { + prop->SetParameters (par2.X(), par2.Y()); + if (!prop->IsCurvatureDefined()) + { (*testout) << "curvature not defined!" << endl; return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); + } + curvature = max(curvature,max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature()))); - //(*testout) << "curvature " << curvature << endl; + //(*testout) << "curvature " << curvature << endl; - if (curvature < 1e-3) - { + if (curvature < 1e-3) + { //(*testout) << "curvature too small (" << curvature << ")!" << endl; return; // return war bis 10.2.05 auskommentiert - } + } - h = ComputeH (curvature+1e-10); + h = ComputeH (curvature+1e-10, mparam); - if(h < 1e-4*maxside) - return; + if(h < 1e-4*maxside) + return; - if (h > 30) return; + if (h > 30) return; } - if (h < maxside && depth < 10) + if (h < maxside && depth < 10) { - //cout << "\r h " << h << flush; - gp_Pnt2d pm; + //cout << "\r h " << h << flush; + gp_Pnt2d pm; - //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl; - //cout << "par0 " << par0.X() << " " << par0.Y() - //<< " par1 " << par1.X() << " " << par1.Y() - // << " par2 " << par2.X() << " " << par2.Y()<< endl; + //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl; + //cout << "par0 " << par0.X() << " " << par0.Y() + //<< " par1 " << par1.X() << " " << par1.Y() + // << " par2 " << par2.X() << " " << par2.Y()<< endl; - if(ls == 0) - { + if(ls == 0) + { pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y())); - RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); - RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); - } - else if(ls == 1) - { + RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h, mparam); + RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h, mparam); + } + else if(ls == 1) + { pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y())); - RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); - RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); - } - else if(ls == 2) - { + RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h, mparam); + RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h, mparam); + } + else if(ls == 2) + { pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y())); - RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); - RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); - } + RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h, mparam); + RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h, mparam); + } } - else + else { - gp_Pnt pnt; - Point3d p3d; + gp_Pnt pnt; + Point3d p3d; - prop->SetParameters (parmid.X(), parmid.Y()); - pnt = prop->Value(); - p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); - mesh.RestrictLocalH (p3d, h); + prop->SetParameters (parmid.X(), parmid.Y()); + pnt = prop->Value(); + p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); + mesh.RestrictLocalH (p3d, h); - p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z()); - mesh.RestrictLocalH (p3d, h); + p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z()); + mesh.RestrictLocalH (p3d, h); - p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z()); - mesh.RestrictLocalH (p3d, h); + p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z()); + mesh.RestrictLocalH (p3d, h); - p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z()); - mesh.RestrictLocalH (p3d, h); + p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z()); + mesh.RestrictLocalH (p3d, h); - //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl; + //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl; } - } + } - void DivideEdge (TopoDS_Edge & edge, NgArray & ps, - NgArray & params, Mesh & mesh) - { - double s0, s1; - double maxh = mparam.maxh; - int nsubedges = 1; - gp_Pnt pnt, oldpnt; - double svalue[DIVIDEEDGESECTIONS]; + void DivideEdge (TopoDS_Edge & edge, NgArray & ps, + NgArray & params, Mesh & mesh, + const MeshingParameters & mparam) + { + double s0, s1; + double maxh = mparam.maxh; + int nsubedges = 1; + gp_Pnt pnt, oldpnt; + double svalue[DIVIDEEDGESECTIONS]; - GProp_GProps system; - BRepGProp::LinearProperties(edge, system); - double L = system.Mass(); + GProp_GProps system; + BRepGProp::LinearProperties(edge, system); + double L = system.Mass(); - Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); + Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); - double hvalue[DIVIDEEDGESECTIONS+1]; - hvalue[0] = 0; - pnt = c->Value(s0); + double hvalue[DIVIDEEDGESECTIONS+1]; + hvalue[0] = 0; + pnt = c->Value(s0); - double olddist = 0; - double dist = 0; + double olddist = 0; + double dist = 0; - int tmpVal = (int)(DIVIDEEDGESECTIONS); + int tmpVal = (int)(DIVIDEEDGESECTIONS); - for (int i = 1; i <= tmpVal; i++) + for (int i = 1; i <= tmpVal; i++) { - oldpnt = pnt; - pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0)); - hvalue[i] = hvalue[i-1] + - 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* - pnt.Distance(oldpnt); + oldpnt = pnt; + pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0)); + hvalue[i] = hvalue[i-1] + + 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* + pnt.Distance(oldpnt); - //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) - // << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl; + //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) + // << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl; - olddist = dist; - dist = pnt.Distance(oldpnt); + olddist = dist; + dist = pnt.Distance(oldpnt); } - // nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS])); - nsubedges = max (1, int(floor(hvalue[DIVIDEEDGESECTIONS]+0.5))); + // nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS])); + nsubedges = max (1, int(floor(hvalue[DIVIDEEDGESECTIONS]+0.5))); - ps.SetSize(nsubedges-1); - params.SetSize(nsubedges+1); + ps.SetSize(nsubedges-1); + params.SetSize(nsubedges+1); - int i = 1; - int i1 = 0; - do + int i = 1; + int i1 = 0; + do { - if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i) - { + if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i) + { params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0); pnt = c->Value(params[i]); ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z())); i++; - } - i1++; - if (i1 > DIVIDEEDGESECTIONS) - { + } + i1++; + if (i1 > DIVIDEEDGESECTIONS) + { nsubedges = i; ps.SetSize(nsubedges-1); params.SetSize(nsubedges+1); cout << "divide edge: local h too small" << endl; - } + } } while (i < nsubedges); - params[0] = s0; - params[nsubedges] = s1; + params[0] = s0; + params[nsubedges] = s1; - if (params[nsubedges] <= params[nsubedges-1]) + if (params[nsubedges] <= params[nsubedges-1]) { - cout << "CORRECTED" << endl; - ps.SetSize (nsubedges-2); - params.SetSize (nsubedges); - params[nsubedges] = s1; + cout << "CORRECTED" << endl; + ps.SetSize (nsubedges-2); + params.SetSize (nsubedges); + params[nsubedges] = s1; } - } + } - void OCCFindEdges (OCCGeometry & geom, Mesh & mesh) - { - const char * savetask = multithread.task; - multithread.task = "Edge meshing"; + void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam) + { + static Timer t("OCCFindEdges"); RegionTimer r(t); + static Timer tsearch("OCCFindEdges - search point"); + const char * savetask = multithread.task; + multithread.task = "Edge meshing"; - (*testout) << "edge meshing" << endl; + (*testout) << "edge meshing" << endl; - int nvertices = geom.vmap.Extent(); - int nedges = geom.emap.Extent(); + int nvertices = geom.vmap.Extent(); + int nedges = geom.emap.Extent(); - (*testout) << "nvertices = " << nvertices << endl; - (*testout) << "nedges = " << nedges << endl; + (*testout) << "nvertices = " << nvertices << endl; + (*testout) << "nedges = " << nedges << endl; - double eps = 1e-6 * geom.GetBoundingBox().Diam(); + double eps = 1e-6 * geom.GetBoundingBox().Diam(); - for (int i = 1; i <= nvertices; i++) + tsearch.Start(); + for (int i = 1; i <= nvertices; i++) { - gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i))); - MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) ); + gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i))); + MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) ); - bool exists = 0; - if (merge_solids) - for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++) - if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps) - { - exists = 1; - break; - } + bool exists = 0; + if (merge_solids) + for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++) + if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps) + { + exists = 1; + break; + } - if (!exists) - mesh.AddPoint (mp); + if (!exists) + mesh.AddPoint (mp); + } + tsearch.Stop(); + + (*testout) << "different vertices = " << mesh.GetNP() << endl; + + + int first_ep = mesh.GetNP()+1; + + NgArray face2solid[2]; + for (int i = 0; i<2; i++) + { + face2solid[i].SetSize (geom.fmap.Extent()); + face2solid[i] = 0; } - (*testout) << "different vertices = " << mesh.GetNP() << endl; - - - int first_ep = mesh.GetNP()+1; - - NgArray face2solid[2]; - for (int i = 0; i<2; i++) + int solidnr = 0; + for (TopExp_Explorer exp0(geom.shape, TopAbs_SOLID); exp0.More(); exp0.Next()) { - face2solid[i].SetSize (geom.fmap.Extent()); - face2solid[i] = 0; - } - - int solidnr = 0; - for (TopExp_Explorer exp0(geom.shape, TopAbs_SOLID); exp0.More(); exp0.Next()) - { - solidnr++; - for (TopExp_Explorer exp1(exp0.Current(), TopAbs_FACE); exp1.More(); exp1.Next()) - { + solidnr++; + for (TopExp_Explorer exp1(exp0.Current(), TopAbs_FACE); exp1.More(); exp1.Next()) + { TopoDS_Face face = TopoDS::Face(exp1.Current()); int facenr = geom.fmap.FindIndex(face); if (face2solid[0][facenr-1] == 0) - face2solid[0][facenr-1] = solidnr; + face2solid[0][facenr-1] = solidnr; else - face2solid[1][facenr-1] = solidnr; - } + face2solid[1][facenr-1] = solidnr; + } } - int total = 0; - for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) - for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next()) - for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next()) - total++; + int total = 0; + for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) + for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next()) + for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next()) + total++; - int facenr = 0; - int edgenr = 0; + int facenr = 0; + int edgenr = 0; - (*testout) << "faces = " << geom.fmap.Extent() << endl; - int curr = 0; + (*testout) << "faces = " << geom.fmap.Extent() << endl; + int curr = 0; - for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) + for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) { - TopoDS_Face face = TopoDS::Face(geom.fmap(i3)); - facenr = geom.fmap.FindIndex (face); // sollte doch immer == i3 sein ??? JS + TopoDS_Face face = TopoDS::Face(geom.fmap(i3)); + facenr = geom.fmap.FindIndex (face); // sollte doch immer == i3 sein ??? JS - int solidnr0 = face2solid[0][i3-1]; - int solidnr1 = face2solid[1][i3-1]; + int solidnr0 = face2solid[0][i3-1]; + int solidnr1 = face2solid[1][i3-1]; - /* auskommentiert am 3.3.05 von robert - for (exp2.Init (geom.somap(solidnr0), TopAbs_FACE); exp2.More(); exp2.Next()) - { - TopoDS_Face face2 = TopoDS::Face(exp2.Current()); - if (geom.fmap.FindIndex(face2) == facenr) - { - // if (face.Orientation() != face2.Orientation()) swap (solidnr0, solidnr1); - } - } - */ + /* auskommentiert am 3.3.05 von robert + for (exp2.Init (geom.somap(solidnr0), TopAbs_FACE); exp2.More(); exp2.Next()) + { + TopoDS_Face face2 = TopoDS::Face(exp2.Current()); + if (geom.fmap.FindIndex(face2) == facenr) + { + // if (face.Orientation() != face2.Orientation()) swap (solidnr0, solidnr1); + } + } + */ - mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0)); + mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0)); - // Philippose - 06/07/2009 - // Add the face colour to the mesh data - Quantity_Color face_colour; + // Philippose - 06/07/2009 + // Add the face colour to the mesh data + Quantity_Color face_colour; - if(!(geom.face_colours.IsNull()) - && (geom.face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour))) - { + if(!(geom.face_colours.IsNull()) + && (geom.face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour))) + { mesh.GetFaceDescriptor(facenr).SetSurfColour(Vec3d(face_colour.Red(),face_colour.Green(),face_colour.Blue())); - } - else - { + } + else + { mesh.GetFaceDescriptor(facenr).SetSurfColour(Vec3d(0.0,1.0,0.0)); - } + } - if(geom.fnames.Size()>=facenr) - mesh.GetFaceDescriptor(facenr).SetBCName(&geom.fnames[facenr-1]); - mesh.GetFaceDescriptor(facenr).SetBCProperty(facenr); - // ACHTUNG! STIMMT NICHT ALLGEMEIN (RG) + if(geom.fnames.Size()>=facenr) + mesh.GetFaceDescriptor(facenr).SetBCName(&geom.fnames[facenr-1]); + mesh.GetFaceDescriptor(facenr).SetBCProperty(facenr); + // ACHTUNG! STIMMT NICHT ALLGEMEIN (RG) - Handle(Geom_Surface) occface = BRep_Tool::Surface(face); + Handle(Geom_Surface) occface = BRep_Tool::Surface(face); - for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next()) - { + for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next()) + { TopoDS_Shape wire = exp2.Current(); for (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next()) - { - curr++; - (*testout) << "edge nr " << curr << endl; + { + curr++; + (*testout) << "edge nr " << curr << endl; - multithread.percent = 100 * curr / double (total); - if (multithread.terminate) return; + multithread.percent = 100 * curr / double (total); + if (multithread.terminate) return; - TopoDS_Edge edge = TopoDS::Edge (exp3.Current()); - if (BRep_Tool::Degenerated(edge)) - { - //(*testout) << "ignoring degenerated edge" << endl; - continue; - } - - if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) == - geom.vmap.FindIndex(TopExp::LastVertex (edge))) - { - GProp_GProps system; - BRepGProp::LinearProperties(edge, system); - - if (system.Mass() < eps) + TopoDS_Edge edge = TopoDS::Edge (exp3.Current()); + if (BRep_Tool::Degenerated(edge)) { - cout << "ignoring edge " << geom.emap.FindIndex (edge) - << ". closed edge with length < " << eps << endl; - continue; + //(*testout) << "ignoring degenerated edge" << endl; + continue; + } + + if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) == + geom.vmap.FindIndex(TopExp::LastVertex (edge))) + { + GProp_GProps system; + BRepGProp::LinearProperties(edge, system); + + if (system.Mass() < eps) + { + cout << "ignoring edge " << geom.emap.FindIndex (edge) + << ". closed edge with length < " << eps << endl; + continue; + } } - } - Handle(Geom2d_Curve) cof; - double s0, s1; - cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); + Handle(Geom2d_Curve) cof; + double s0, s1; + cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); - int geomedgenr = geom.emap.FindIndex(edge); + int geomedgenr = geom.emap.FindIndex(edge); - NgArray mp; - NgArray params; + NgArray mp; + NgArray params; - DivideEdge (edge, mp, params, mesh); + DivideEdge (edge, mp, params, mesh, mparam); - NgArray pnums; - pnums.SetSize (mp.Size()+2); + NgArray pnums; + pnums.SetSize (mp.Size()+2); - if (!merge_solids) - { - pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)); - pnums[pnums.Size()-1] = geom.vmap.FindIndex (TopExp::LastVertex (edge)); - } - else - { - Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge))); - Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge))); - - pnums[0] = -1; - pnums.Last() = -1; - for (PointIndex pi = 1; pi < first_ep; pi++) + if (!merge_solids) { - if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi; - if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi; + pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)); + pnums[pnums.Size()-1] = geom.vmap.FindIndex (TopExp::LastVertex (edge)); + } + else + { + Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge))); + Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge))); + + pnums[0] = -1; + pnums.Last() = -1; + for (PointIndex pi = 1; pi < first_ep; pi++) + { + if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi; + if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi; + } } - } - for (int i = 1; i <= mp.Size(); i++) - { - bool exists = 0; - int j; - for (j = first_ep; j <= mesh.GetNP(); j++) - if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) - { - exists = 1; - break; - } - - if (exists) - pnums[i] = j; - else - { + for (int i = 1; i <= mp.Size(); i++) + { + bool exists = 0; + tsearch.Start(); + int j; + for (j = first_ep; j <= mesh.GetNP(); j++) + if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) + { + exists = 1; + break; + } + tsearch.Stop(); + + if (exists) + pnums[i] = j; + else + { mesh.AddPoint (mp[i-1]); (*testout) << "add meshpoint " << mp[i-1] << endl; pnums[i] = mesh.GetNP(); - } - } - (*testout) << "NP = " << mesh.GetNP() << endl; - - //(*testout) << pnums[pnums.Size()-1] << endl; - - for (int i = 1; i <= mp.Size()+1; i++) - { - edgenr++; - Segment seg; - - seg[0] = pnums[i-1]; - seg[1] = pnums[i]; - seg.edgenr = edgenr; - seg.si = facenr; - seg.epgeominfo[0].dist = params[i-1]; - seg.epgeominfo[1].dist = params[i]; - seg.epgeominfo[0].edgenr = geomedgenr; - seg.epgeominfo[1].edgenr = geomedgenr; - - gp_Pnt2d p2d; - p2d = cof->Value(params[i-1]); - // if (i == 1) p2d = cof->Value(s0); - seg.epgeominfo[0].u = p2d.X(); - seg.epgeominfo[0].v = p2d.Y(); - p2d = cof->Value(params[i]); - // if (i == mp.Size()+1) p2d = cof -> Value(s1); - seg.epgeominfo[1].u = p2d.X(); - seg.epgeominfo[1].v = p2d.Y(); - - /* - if (occface->IsUPeriodic()) - { - cout << "U Periodic" << endl; - if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > - fabs(seg.epgeominfo[1].u- - (seg.epgeominfo[0].u-occface->UPeriod()))) - seg.epgeominfo[0].u = p2d.X()+occface->UPeriod(); - - if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > - fabs(seg.epgeominfo[1].u- - (seg.epgeominfo[0].u+occface->UPeriod()))) - seg.epgeominfo[0].u = p2d.X()-occface->UPeriod(); + } } + (*testout) << "NP = " << mesh.GetNP() << endl; - if (occface->IsVPeriodic()) + //(*testout) << pnums[pnums.Size()-1] << endl; + + for (int i = 1; i <= mp.Size()+1; i++) { - cout << "V Periodic" << endl; - if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > - fabs(seg.epgeominfo[1].v- - (seg.epgeominfo[0].v-occface->VPeriod()))) - seg.epgeominfo[0].v = p2d.Y()+occface->VPeriod(); + edgenr++; + Segment seg; + + seg[0] = pnums[i-1]; + seg[1] = pnums[i]; + seg.edgenr = edgenr; + seg.si = facenr; + seg.epgeominfo[0].dist = params[i-1]; + seg.epgeominfo[1].dist = params[i]; + seg.epgeominfo[0].edgenr = geomedgenr; + seg.epgeominfo[1].edgenr = geomedgenr; + + gp_Pnt2d p2d; + p2d = cof->Value(params[i-1]); + // if (i == 1) p2d = cof->Value(s0); + seg.epgeominfo[0].u = p2d.X(); + seg.epgeominfo[0].v = p2d.Y(); + p2d = cof->Value(params[i]); + // if (i == mp.Size()+1) p2d = cof -> Value(s1); + seg.epgeominfo[1].u = p2d.X(); + seg.epgeominfo[1].v = p2d.Y(); + + /* + if (occface->IsUPeriodic()) + { + cout << "U Periodic" << endl; + if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > + fabs(seg.epgeominfo[1].u- + (seg.epgeominfo[0].u-occface->UPeriod()))) + seg.epgeominfo[0].u = p2d.X()+occface->UPeriod(); + + if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > + fabs(seg.epgeominfo[1].u- + (seg.epgeominfo[0].u+occface->UPeriod()))) + seg.epgeominfo[0].u = p2d.X()-occface->UPeriod(); + } + + if (occface->IsVPeriodic()) + { + cout << "V Periodic" << endl; + if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > + fabs(seg.epgeominfo[1].v- + (seg.epgeominfo[0].v-occface->VPeriod()))) + seg.epgeominfo[0].v = p2d.Y()+occface->VPeriod(); + + if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > + fabs(seg.epgeominfo[1].v- + (seg.epgeominfo[0].v+occface->VPeriod()))) + seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod(); + } + */ + + if (edge.Orientation() == TopAbs_REVERSED) + { + swap (seg[0], seg[1]); + swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist); + swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); + swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); + } + + mesh.AddSegment (seg); + + //edgesegments[geomedgenr-1]->Append(mesh.GetNSeg()); - if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > - fabs(seg.epgeominfo[1].v- - (seg.epgeominfo[0].v+occface->VPeriod()))) - seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod(); } - */ - - if (edge.Orientation() == TopAbs_REVERSED) - { - swap (seg[0], seg[1]); - swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist); - swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); - swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); - } - - mesh.AddSegment (seg); - - //edgesegments[geomedgenr-1]->Append(mesh.GetNSeg()); - - } - } - } + } + } } - // for(i=1; i<=mesh.GetNSeg(); i++) - // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si - // << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl; - // exit(10); + // for(i=1; i<=mesh.GetNSeg(); i++) + // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si + // << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl; + // exit(10); - mesh.CalcSurfacesOfNode(); - multithread.task = savetask; - } + mesh.CalcSurfacesOfNode(); + multithread.task = savetask; + } - void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend) - { - int i, j, k; - int changed; + void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend, + MeshingParameters & mparam) + { + static Timer t("OCCMeshSurface"); RegionTimer r(t); + + // int i, j, k; + // int changed; - const char * savetask = multithread.task; - multithread.task = "Surface meshing"; + const char * savetask = multithread.task; + multithread.task = "Surface meshing"; - geom.facemeshstatus = 0; + geom.facemeshstatus = 0; - int noldp = mesh.GetNP(); + int noldp = mesh.GetNP(); - double starttime = GetTime(); + double starttime = GetTime(); - NgArray glob2loc(noldp); + NgArray glob2loc(noldp); - //int projecttype = PARAMETERSPACE; + //int projecttype = PARAMETERSPACE; - int projecttype = PARAMETERSPACE; + int projecttype = PARAMETERSPACE; - int notrys = 1; + int notrys = 1; - int surfmesherror = 0; + int surfmesherror = 0; - for (k = 1; k <= mesh.GetNFD(); k++) + for (int k = 1; k <= mesh.GetNFD(); k++) { - if(1==0 && !geom.fvispar[k-1].IsDrawable()) - { + if(1==0 && !geom.fvispar[k-1].IsDrawable()) + { (*testout) << "ignoring face " << k << endl; cout << "ignoring face " << k << endl; continue; - } + } - (*testout) << "mesh face " << k << endl; - multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); - geom.facemeshstatus[k-1] = -1; + (*testout) << "mesh face " << k << endl; + multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); + geom.facemeshstatus[k-1] = -1; - /* - if (k != 42) - { - cout << "skipped" << endl; - continue; - } - */ + /* + if (k != 42) + { + cout << "skipped" << endl; + continue; + } + */ + + FaceDescriptor & fd = mesh.GetFaceDescriptor(k); + + int oldnf = mesh.GetNSE(); + + Box<3> bb = geom.GetBoundingBox(); + + // int projecttype = PLANESPACE; + + static Timer tinit("init"); + tinit.Start(); + Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam); + tinit.Stop(); - FaceDescriptor & fd = mesh.GetFaceDescriptor(k); + static Timer tprint("print"); + tprint.Start(); + if (meshing.GetProjectionType() == PLANESPACE) + PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); + else + PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)"); + tprint.Stop(); + if (surfmesherror) + cout << "Surface meshing error occurred before (in " << surfmesherror << " faces)" << endl; - int oldnf = mesh.GetNSE(); - - Box<3> bb = geom.GetBoundingBox(); - - // int projecttype = PLANESPACE; - - Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype); - - if (meshing.GetProjectionType() == PLANESPACE) - PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); - else - PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)"); - - if (surfmesherror) - cout << "Surface meshing error occurred before (in " << surfmesherror << " faces)" << endl; - - // Meshing2OCCSurfaces meshing(f2, bb); - meshing.SetStartTime (starttime); - - //(*testout) << "Face " << k << endl << endl; + // Meshing2OCCSurfaces meshing(f2, bb); + meshing.SetStartTime (starttime); + //(*testout) << "Face " << k << endl << endl; - if (meshing.GetProjectionType() == PLANESPACE) - { + if (meshing.GetProjectionType() == PLANESPACE) + { + static Timer t("MeshSurface: Find edges and points - Physical"); RegionTimer r(t); int cntp = 0; glob2loc = 0; - for (i = 1; i <= mesh.GetNSeg(); i++) - { - Segment & seg = mesh.LineSegment(i); - if (seg.si == k) - { - for (j = 1; j <= 2; j++) + for (int i = 1; i <= mesh.GetNSeg(); i++) + { + Segment & seg = mesh.LineSegment(i); + if (seg.si == k) { - int pi = (j == 1) ? seg[0] : seg[1]; - if (!glob2loc.Get(pi)) - { - meshing.AddPoint (mesh.Point(pi), pi); - cntp++; - glob2loc.Elem(pi) = cntp; - } + for (int j = 1; j <= 2; j++) + { + int pi = (j == 1) ? seg[0] : seg[1]; + if (!glob2loc.Get(pi)) + { + meshing.AddPoint (mesh.Point(pi), pi); + cntp++; + glob2loc.Elem(pi) = cntp; + } + } } - } - } + } - for (i = 1; i <= mesh.GetNSeg(); i++) - { - Segment & seg = mesh.LineSegment(i); - if (seg.si == k) - { - PointGeomInfo gi0, gi1; - gi0.trignum = gi1.trignum = k; - gi0.u = seg.epgeominfo[0].u; - gi0.v = seg.epgeominfo[0].v; - gi1.u = seg.epgeominfo[1].u; - gi1.v = seg.epgeominfo[1].v; + for (int i = 1; i <= mesh.GetNSeg(); i++) + { + Segment & seg = mesh.LineSegment(i); + if (seg.si == k) + { + PointGeomInfo gi0, gi1; + gi0.trignum = gi1.trignum = k; + gi0.u = seg.epgeominfo[0].u; + gi0.v = seg.epgeominfo[0].v; + gi1.u = seg.epgeominfo[1].u; + gi1.v = seg.epgeominfo[1].v; - meshing.AddBoundaryElement (glob2loc.Get(seg[0]), glob2loc.Get(seg[1]), gi0, gi1); - //(*testout) << gi0.u << " " << gi0.v << endl; - //(*testout) << gi1.u << " " << gi1.v << endl; - } - } - } - else - { + meshing.AddBoundaryElement (glob2loc.Get(seg[0]), glob2loc.Get(seg[1]), gi0, gi1); + //(*testout) << gi0.u << " " << gi0.v << endl; + //(*testout) << gi1.u << " " << gi1.v << endl; + } + } + } + else + { + static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t); + int cntp = 0; - for (i = 1; i <= mesh.GetNSeg(); i++) - if (mesh.LineSegment(i).si == k) - cntp+=2; + for (int i = 1; i <= mesh.GetNSeg(); i++) + if (mesh.LineSegment(i).si == k) + cntp+=2; NgArray< PointGeomInfo > gis; @@ -731,304 +741,310 @@ namespace netgen gis.SetAllocSize (cntp); gis.SetSize (0); - for (i = 1; i <= mesh.GetNSeg(); i++) - { - Segment & seg = mesh.LineSegment(i); - if (seg.si == k) - { - PointGeomInfo gi0, gi1; - gi0.trignum = gi1.trignum = k; - gi0.u = seg.epgeominfo[0].u; - gi0.v = seg.epgeominfo[0].v; - gi1.u = seg.epgeominfo[1].u; - gi1.v = seg.epgeominfo[1].v; - - int locpnum[2] = {0, 0}; - - for (j = 0; j < 2; j++) + for (int i = 1; i <= mesh.GetNSeg(); i++) + { + Segment & seg = mesh.LineSegment(i); + if (seg.si == k) { - PointGeomInfo gi = (j == 0) ? gi0 : gi1; + PointGeomInfo gi0, gi1; + gi0.trignum = gi1.trignum = k; + gi0.u = seg.epgeominfo[0].u; + gi0.v = seg.epgeominfo[0].v; + gi1.u = seg.epgeominfo[1].u; + gi1.v = seg.epgeominfo[1].v; - int l; - for (l = 0; l < gis.Size() && locpnum[j] == 0; l++) - { - double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); + int locpnum[2] = {0, 0}; - if (dist < 1e-10) - locpnum[j] = l+1; - } + for (int j = 0; j < 2; j++) + { + PointGeomInfo gi = (j == 0) ? gi0 : gi1; - if (locpnum[j] == 0) - { - int pi = (j == 0) ? seg[0] : seg[1]; - meshing.AddPoint (mesh.Point(pi), pi); + int l; + for (l = 0; l < gis.Size() && locpnum[j] == 0; l++) + { + double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); + + if (dist < 1e-10) + locpnum[j] = l+1; + } + + if (locpnum[j] == 0) + { + int pi = (j == 0) ? seg[0] : seg[1]; + meshing.AddPoint (mesh.Point(pi), pi); + + gis.SetSize (gis.Size()+1); + gis[l] = gi; + locpnum[j] = l+1; + } + } + + meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1); + //(*testout) << gi0.u << " " << gi0.v << endl; + //(*testout) << gi1.u << " " << gi1.v << endl; - gis.SetSize (gis.Size()+1); - gis[l] = gi; - locpnum[j] = l+1; - } } - - meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1); - //(*testout) << gi0.u << " " << gi0.v << endl; - //(*testout) << gi1.u << " " << gi1.v << endl; - - } - } - } + } + } - // Philippose - 15/01/2009 - double maxh = geom.face_maxh[k-1]; - //double maxh = mparam.maxh; - mparam.checkoverlap = 0; - // int noldpoints = mesh->GetNP(); - int noldsurfel = mesh.GetNSE(); + // Philippose - 15/01/2009 + double maxh = geom.face_maxh[k-1]; + //double maxh = mparam.maxh; + mparam.checkoverlap = 0; + // int noldpoints = mesh->GetNP(); + int noldsurfel = mesh.GetNSE(); - GProp_GProps sprops; - BRepGProp::SurfaceProperties(TopoDS::Face(geom.fmap(k)),sprops); - meshing.SetMaxArea(2.*sprops.Mass()); + static Timer tsurfprop("surfprop"); + tsurfprop.Start(); + GProp_GProps sprops; + BRepGProp::SurfaceProperties(TopoDS::Face(geom.fmap(k)),sprops); + tsurfprop.Stop(); + meshing.SetMaxArea(2.*sprops.Mass()); - MESHING2_RESULT res; + MESHING2_RESULT res; - try { - res = meshing.GenerateMesh (mesh, mparam, maxh, k); - } + try { + static Timer t("GenerateMesh"); RegionTimer reg(t); + res = meshing.GenerateMesh (mesh, mparam, maxh, k); + } - catch (SingularMatrixException) - { + catch (SingularMatrixException) + { (*myerr) << "Singular Matrix" << endl; res = MESHING2_GIVEUP; - } + } - catch (UVBoundsException) - { + catch (UVBoundsException) + { (*myerr) << "UV bounds exceeded" << endl; res = MESHING2_GIVEUP; - } + } - projecttype = PARAMETERSPACE; - - if (res != MESHING2_OK) - { + projecttype = PARAMETERSPACE; + static Timer t1("rest of loop"); RegionTimer reg1(t1); + + if (res != MESHING2_OK) + { if (notrys == 1) - { - for (int i = noldsurfel+1; i <= mesh.GetNSE(); i++) + { + for (int i = noldsurfel+1; i <= mesh.GetNSE(); i++) mesh.DeleteSurfaceElement (i); - mesh.Compress(); + mesh.Compress(); - cout << "retry Surface " << k << endl; + cout << "retry Surface " << k << endl; - k--; - projecttype*=-1; - notrys++; - continue; - } + k--; + projecttype*=-1; + notrys++; + continue; + } else - { - geom.facemeshstatus[k-1] = -1; - PrintError ("Problem in Surface mesh generation"); - surfmesherror++; - // throw NgException ("Problem in Surface mesh generation"); - } - } - else - { + { + geom.facemeshstatus[k-1] = -1; + PrintError ("Problem in Surface mesh generation"); + surfmesherror++; + // throw NgException ("Problem in Surface mesh generation"); + } + } + else + { geom.facemeshstatus[k-1] = 1; - } + } - notrys = 1; + notrys = 1; - for (i = oldnf+1; i <= mesh.GetNSE(); i++) - mesh.SurfaceElement(i).SetIndex (k); + for (int i = oldnf+1; i <= mesh.GetNSE(); i++) + mesh.SurfaceElement(i).SetIndex (k); } -// ofstream problemfile("occmesh.rep"); + // ofstream problemfile("occmesh.rep"); -// problemfile << "SURFACEMESHING" << endl << endl; + // problemfile << "SURFACEMESHING" << endl << endl; - if (surfmesherror) + if (surfmesherror) { - cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl; - cout << "SURFACE MESHING ERROR OCCURRED IN " << surfmesherror << " FACES:" << endl; - for (int i = 1; i <= geom.fmap.Extent(); i++) - if (geom.facemeshstatus[i-1] == -1) + cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl; + cout << "SURFACE MESHING ERROR OCCURRED IN " << surfmesherror << " FACES:" << endl; + for (int i = 1; i <= geom.fmap.Extent(); i++) + if (geom.facemeshstatus[i-1] == -1) { - cout << "Face " << i << endl; -// problemfile << "problem with face " << i << endl; -// problemfile << "vertices: " << endl; - TopExp_Explorer exp0,exp1,exp2; - for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() ) - { + cout << "Face " << i << endl; + // problemfile << "problem with face " << i << endl; + // problemfile << "vertices: " << endl; + TopExp_Explorer exp0,exp1,exp2; + for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() ) + { TopoDS_Wire wire = TopoDS::Wire(exp0.Current()); for ( exp1.Init(wire,TopAbs_EDGE); exp1.More(); exp1.Next() ) - { - TopoDS_Edge edge = TopoDS::Edge(exp1.Current()); - for ( exp2.Init(edge,TopAbs_VERTEX); exp2.More(); exp2.Next() ) - { - TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current()); - gp_Pnt point = BRep_Tool::Pnt(vertex); -// problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl; - } - } - } -// problemfile << endl; + { + TopoDS_Edge edge = TopoDS::Edge(exp1.Current()); + for ( exp2.Init(edge,TopAbs_VERTEX); exp2.More(); exp2.Next() ) + { + TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current()); + gp_Pnt point = BRep_Tool::Pnt(vertex); + // problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl; + } + } + } + // problemfile << endl; } - cout << endl << endl; - cout << "for more information open IGES/STEP Topology Explorer" << endl; -// problemfile.close(); - throw NgException ("Problem in Surface mesh generation"); + cout << endl << endl; + cout << "for more information open IGES/STEP Topology Explorer" << endl; + // problemfile.close(); + throw NgException ("Problem in Surface mesh generation"); } - else + else { -// problemfile << "OK" << endl << endl; -// problemfile.close(); + // problemfile << "OK" << endl << endl; + // problemfile.close(); } - if (multithread.terminate || perfstepsend < MESHCONST_OPTSURFACE) - return; + if (multithread.terminate || perfstepsend < MESHCONST_OPTSURFACE) + return; - multithread.task = "Optimizing surface"; + multithread.task = "Optimizing surface"; - static Timer timer_opt2d("Optimization 2D"); - timer_opt2d.Start(); + static Timer timer_opt2d("Optimization 2D"); + timer_opt2d.Start(); - for (k = 1; k <= mesh.GetNFD(); k++) + for (int k = 1; k <= mesh.GetNFD(); k++) { - // if (k != 42) continue; - // if (k != 36) continue; + // if (k != 42) continue; + // if (k != 36) continue; - // (*testout) << "optimize face " << k << endl; - multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); + // (*testout) << "optimize face " << k << endl; + multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); - FaceDescriptor & fd = mesh.GetFaceDescriptor(k); + FaceDescriptor & fd = mesh.GetFaceDescriptor(k); - PrintMessage (1, "Optimize Surface ", k); - for (i = 1; i <= mparam.optsteps2d; i++) - { + PrintMessage (1, "Optimize Surface ", k); + for (int i = 1; i <= mparam.optsteps2d; i++) + { // (*testout) << "optstep " << i << endl; if (multithread.terminate) return; { - MeshOptimize2dOCCSurfaces meshopt(geom); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - //meshopt.SetMetricWeight (0.2); - meshopt.SetWriteStatus (0); + MeshOptimize2dOCCSurfaces meshopt(geom); + meshopt.SetFaceIndex (k); + meshopt.SetImproveEdges (0); + meshopt.SetMetricWeight (mparam.elsizeweight); + //meshopt.SetMetricWeight (0.2); + meshopt.SetWriteStatus (0); - // (*testout) << "EdgeSwapping (mesh, (i > mparam.optsteps2d/2))" << endl; - meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2)); + // (*testout) << "EdgeSwapping (mesh, (i > mparam.optsteps2d/2))" << endl; + meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2)); } if (multithread.terminate) return; { - MeshOptimize2dOCCSurfaces meshopt(geom); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - //meshopt.SetMetricWeight (0.2); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); + MeshOptimize2dOCCSurfaces meshopt(geom); + meshopt.SetFaceIndex (k); + meshopt.SetImproveEdges (0); + //meshopt.SetMetricWeight (0.2); + meshopt.SetMetricWeight (mparam.elsizeweight); + meshopt.SetWriteStatus (0); - // (*testout) << "ImproveMesh (mesh)" << endl; - meshopt.ImproveMesh (mesh, mparam); + // (*testout) << "ImproveMesh (mesh)" << endl; + meshopt.ImproveMesh (mesh, mparam); } { - MeshOptimize2dOCCSurfaces meshopt(geom); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - //meshopt.SetMetricWeight (0.2); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); + MeshOptimize2dOCCSurfaces meshopt(geom); + meshopt.SetFaceIndex (k); + meshopt.SetImproveEdges (0); + //meshopt.SetMetricWeight (0.2); + meshopt.SetMetricWeight (mparam.elsizeweight); + meshopt.SetWriteStatus (0); - // (*testout) << "CombineImprove (mesh)" << endl; - meshopt.CombineImprove (mesh); + // (*testout) << "CombineImprove (mesh)" << endl; + meshopt.CombineImprove (mesh); } if (multithread.terminate) return; { - MeshOptimize2dOCCSurfaces meshopt(geom); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - //meshopt.SetMetricWeight (0.2); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); + MeshOptimize2dOCCSurfaces meshopt(geom); + meshopt.SetFaceIndex (k); + meshopt.SetImproveEdges (0); + //meshopt.SetMetricWeight (0.2); + meshopt.SetMetricWeight (mparam.elsizeweight); + meshopt.SetWriteStatus (0); - // (*testout) << "ImproveMesh (mesh)" << endl; - meshopt.ImproveMesh (mesh, mparam); + // (*testout) << "ImproveMesh (mesh)" << endl; + meshopt.ImproveMesh (mesh, mparam); } - } + } } - mesh.CalcSurfacesOfNode(); - mesh.Compress(); - timer_opt2d.Stop(); + mesh.CalcSurfacesOfNode(); + mesh.Compress(); + timer_opt2d.Stop(); - multithread.task = savetask; + multithread.task = savetask; - // Gerhard BEGIN - for(int i = 0; i maxhdom; - maxhdom.SetSize (geom.NrSolids()); - maxhdom = mparam.maxh; + NgArray maxhdom; + maxhdom.SetSize (geom.NrSolids()); + maxhdom = mparam.maxh; - mesh.SetMaxHDomain (maxhdom); + mesh.SetMaxHDomain (maxhdom); - Box<3> bb = geom.GetBoundingBox(); - bb.Increase (bb.Diam()/10); + Box<3> bb = geom.GetBoundingBox(); + bb.Increase (bb.Diam()/10); - mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); + mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); - if (mparam.uselocalh) + if (mparam.uselocalh) { - const char * savetask = multithread.task; - multithread.percent = 0; + const char * savetask = multithread.task; + multithread.percent = 0; - mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); + mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); - int nedges = geom.emap.Extent(); + int nedges = geom.emap.Extent(); - double mincurvelength = IGNORECURVELENGTH; - double maxedgelen = 0; - double minedgelen = 1e99; + double mincurvelength = IGNORECURVELENGTH; + double maxedgelen = 0; + double minedgelen = 1e99; - if(occparam.resthminedgelenenable) - { - mincurvelength = occparam.resthminedgelen; - if(mincurvelength < IGNORECURVELENGTH) mincurvelength = IGNORECURVELENGTH; - } + if(occparam.resthminedgelenenable) + { + mincurvelength = occparam.resthminedgelen; + if(mincurvelength < IGNORECURVELENGTH) mincurvelength = IGNORECURVELENGTH; + } - multithread.task = "Setting local mesh size (elements per edge)"; + multithread.task = "Setting local mesh size (elements per edge)"; - // setting elements per edge + // setting elements per edge - for (int i = 1; i <= nedges && !multithread.terminate; i++) - { + for (int i = 1; i <= nedges && !multithread.terminate; i++) + { TopoDS_Edge e = TopoDS::Edge (geom.emap(i)); multithread.percent = 100 * (i-1)/double(nedges); if (BRep_Tool::Degenerated(e)) continue; @@ -1038,10 +1054,10 @@ namespace netgen double len = system.Mass(); if (len < mincurvelength) - { - (*testout) << "ignored" << endl; - continue; - } + { + (*testout) << "ignored" << endl; + continue; + } double localh = len/mparam.segmentsperedge; double s0, s1; @@ -1059,13 +1075,13 @@ namespace netgen TopTools_ListIteratorOfListOfShape parent_face_list; for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next()) - { - TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value()); + { + TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value()); - int face_index = geom.fmap.FindIndex(parent_face); + int face_index = geom.fmap.FindIndex(parent_face); - if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); - } + if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); + } Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); @@ -1082,20 +1098,20 @@ namespace netgen int maxj = max((int) ceil(len/localh), 2); for (int j = 0; j <= maxj; j++) - { - gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); - mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh); - } - } + { + gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); + mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh); + } + } - multithread.task = "Setting local mesh size (edge curvature)"; + multithread.task = "Setting local mesh size (edge curvature)"; - // setting edge curvature + // setting edge curvature - int nsections = 20; + int nsections = 20; - for (int i = 1; i <= nedges && !multithread.terminate; i++) - { + for (int i = 1; i <= nedges && !multithread.terminate; i++) + { double maxcur = 0; multithread.percent = 100 * (i-1)/double(nedges); TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); @@ -1106,30 +1122,30 @@ namespace netgen BRepLProp_CLProps prop(brepc, 2, 1e-5); for (int j = 1; j <= nsections; j++) - { - double s = s0 + j/(double) nsections * (s1-s0); - prop.SetParameter (s); - double curvature = prop.Curvature(); - if(curvature> maxcur) maxcur = curvature; + { + double s = s0 + j/(double) nsections * (s1-s0); + prop.SetParameter (s); + double curvature = prop.Curvature(); + if(curvature> maxcur) maxcur = curvature; - if (curvature >= 1e99) + if (curvature >= 1e99) continue; - gp_Pnt pnt = c->Value (s); + gp_Pnt pnt = c->Value (s); - mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature))); - } + mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature), mparam)); + } // (*testout) << "edge " << i << " max. curvature: " << maxcur << endl; - } + } - multithread.task = "Setting local mesh size (face curvature)"; + multithread.task = "Setting local mesh size (face curvature)"; - // setting face curvature + // setting face curvature - int nfaces = geom.fmap.Extent(); + int nfaces = geom.fmap.Extent(); - for (int i = 1; i <= nfaces && !multithread.terminate; i++) - { + for (int i = 1; i <= nfaces && !multithread.terminate; i++) + { multithread.percent = 100 * (i-1)/double(nfaces); TopoDS_Face face = TopoDS::Face(geom.fmap(i)); TopLoc_Location loc; @@ -1148,32 +1164,32 @@ namespace netgen int ntriangles = triangulation -> NbTriangles(); for (int j = 1; j <= ntriangles; j++) - { - gp_Pnt p[3]; - gp_Pnt2d par[3]; + { + gp_Pnt p[3]; + gp_Pnt2d par[3]; - for (int k = 1; k <=3; k++) - { - int n = triangulation->Triangles()(j)(k); - p[k-1] = triangulation->Nodes()(n).Transformed(loc); - par[k-1] = triangulation->UVNodes()(n); - } + for (int k = 1; k <=3; k++) + { + int n = triangulation->Triangles()(j)(k); + p[k-1] = triangulation->Nodes()(n).Transformed(loc); + par[k-1] = triangulation->UVNodes()(n); + } - //double maxside = 0; - //maxside = max (maxside, p[0].Distance(p[1])); - //maxside = max (maxside, p[0].Distance(p[2])); - //maxside = max (maxside, p[1].Distance(p[2])); - //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; + //double maxside = 0; + //maxside = max (maxside, p[0].Distance(p[1])); + //maxside = max (maxside, p[0].Distance(p[2])); + //maxside = max (maxside, p[1].Distance(p[2])); + //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; - RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0); - //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; - } - } + RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0, 0, mparam); + //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; + } + } - // setting close edges + // setting close edges - if (occparam.resthcloseedgeenable) - { + if (occparam.resthcloseedgeenable) + { multithread.task = "Setting local mesh size (close edges)"; int sections = 100; @@ -1185,300 +1201,300 @@ namespace netgen int nlines = 0; for (int i = 1; i <= nedges && !multithread.terminate; i++) - { - TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); - if (BRep_Tool::Degenerated(edge)) continue; + { + TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); + if (BRep_Tool::Degenerated(edge)) continue; - double s0, s1; - Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); - BRepAdaptor_Curve brepc(edge); - BRepLProp_CLProps prop(brepc, 1, 1e-5); - prop.SetParameter (s0); + double s0, s1; + Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); + BRepAdaptor_Curve brepc(edge); + BRepLProp_CLProps prop(brepc, 1, 1e-5); + prop.SetParameter (s0); - gp_Vec d0 = prop.D1().Normalized(); - double s_start = s0; - int count = 0; - for (int j = 1; j <= sections; j++) - { - double s = s0 + (s1-s0)*(double)j/(double)sections; - prop.SetParameter (s); - gp_Vec d1 = prop.D1().Normalized(); - double cosalpha = fabs(d0*d1); - if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI))) + gp_Vec d0 = prop.D1().Normalized(); + double s_start = s0; + int count = 0; + for (int j = 1; j <= sections; j++) { - count++; - gp_Pnt p0 = c->Value (s_start); - gp_Pnt p1 = c->Value (s); - lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z()); - lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z()); + double s = s0 + (s1-s0)*(double)j/(double)sections; + prop.SetParameter (s); + gp_Vec d1 = prop.D1().Normalized(); + double cosalpha = fabs(d0*d1); + if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI))) + { + count++; + gp_Pnt p0 = c->Value (s_start); + gp_Pnt p1 = c->Value (s); + lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z()); + lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z()); - Box3d box; - box.SetPoint (Point3d(lines[nlines].p0)); - box.AddPoint (Point3d(lines[nlines].p1)); + Box3d box; + box.SetPoint (Point3d(lines[nlines].p0)); + box.AddPoint (Point3d(lines[nlines].p1)); - searchtree->Insert (box.PMin(), box.PMax(), nlines+1); - nlines++; + searchtree->Insert (box.PMin(), box.PMax(), nlines+1); + nlines++; - s_start = s; - d0 = d1; + s_start = s; + d0 = d1; + } } - } - } + } NgArray linenums; for (int i = 0; i < nlines; i++) - { - multithread.percent = (100*i)/double(nlines); - Line & line = lines[i]; + { + multithread.percent = (100*i)/double(nlines); + Line & line = lines[i]; - Box3d box; - box.SetPoint (Point3d(line.p0)); - box.AddPoint (Point3d(line.p1)); - double maxhline = max (mesh.GetH(box.PMin()), - mesh.GetH(box.PMax())); - box.Increase(maxhline); + Box3d box; + box.SetPoint (Point3d(line.p0)); + box.AddPoint (Point3d(line.p1)); + double maxhline = max (mesh.GetH(box.PMin()), + mesh.GetH(box.PMax())); + box.Increase(maxhline); - double mindist = 1e99; - linenums.SetSize(0); - searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums); + double mindist = 1e99; + linenums.SetSize(0); + searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums); - for (int j = 0; j < linenums.Size(); j++) - { - int num = linenums[j]-1; - if (i == num) continue; - if ((line.p0-lines[num].p0).Length2() < 1e-15) continue; - if ((line.p0-lines[num].p1).Length2() < 1e-15) continue; - if ((line.p1-lines[num].p0).Length2() < 1e-15) continue; - if ((line.p1-lines[num].p1).Length2() < 1e-15) continue; - mindist = min (mindist, line.Dist(lines[num])); - } + for (int j = 0; j < linenums.Size(); j++) + { + int num = linenums[j]-1; + if (i == num) continue; + if ((line.p0-lines[num].p0).Length2() < 1e-15) continue; + if ((line.p0-lines[num].p1).Length2() < 1e-15) continue; + if ((line.p1-lines[num].p0).Length2() < 1e-15) continue; + if ((line.p1-lines[num].p1).Length2() < 1e-15) continue; + mindist = min (mindist, line.Dist(lines[num])); + } - mindist /= (occparam.resthcloseedgefac + VSMALL); + mindist /= (occparam.resthcloseedgefac + VSMALL); - if (mindist < 1e-3) - { - (*testout) << "extremely small local h: " << mindist - << " --> setting to 1e-3" << endl; - (*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl; - mindist = 1e-3; - } + if (mindist < 1e-3) + { + (*testout) << "extremely small local h: " << mindist + << " --> setting to 1e-3" << endl; + (*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl; + mindist = 1e-3; + } - mesh.RestrictLocalHLine(line.p0, line.p1, mindist); - } - } + mesh.RestrictLocalHLine(line.p0, line.p1, mindist); + } + } - multithread.task = savetask; + multithread.task = savetask; } - // Philippose - 09/03/2009 - // Added the capability to load the mesh size from a - // file also for OpenCascade Geometry - // Note: - // ** If the "uselocalh" option is ticked in - // the "mesh options...insider" menu, the mesh - // size will be further modified by the topology - // analysis routines. - // ** To use the mesh size file as the sole source - // for defining the mesh size, uncheck the "uselocalh" - // option. - mesh.LoadLocalMeshSize (mparam.meshsizefilename); - } + // Philippose - 09/03/2009 + // Added the capability to load the mesh size from a + // file also for OpenCascade Geometry + // Note: + // ** If the "uselocalh" option is ticked in + // the "mesh options...insider" menu, the mesh + // size will be further modified by the topology + // analysis routines. + // ** To use the mesh size file as the sole source + // for defining the mesh size, uncheck the "uselocalh" + // option. + mesh.LoadLocalMeshSize (mparam.meshsizefilename); + } int OCCGenerateMesh (OCCGeometry & geom, shared_ptr & mesh, MeshingParameters & mparam) - { - multithread.percent = 0; + { + multithread.percent = 0; - if (mparam.perfstepsstart <= MESHCONST_ANALYSE) + if (mparam.perfstepsstart <= MESHCONST_ANALYSE) { - if(mesh.get() == nullptr) - mesh = make_shared(); - mesh->geomtype = Mesh::GEOM_OCC; + if(mesh.get() == nullptr) + mesh = make_shared(); + mesh->geomtype = Mesh::GEOM_OCC; - OCCSetLocalMeshSize(geom,*mesh); + OCCSetLocalMeshSize(geom,*mesh, mparam); } - if (multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE) - return TCL_OK; - - if (mparam.perfstepsstart <= MESHCONST_MESHEDGES) - { - OCCFindEdges (geom, *mesh); - - /* - cout << "Removing redundant points" << endl; - - int i, j; - int np = mesh->GetNP(); - NgArray equalto; - - equalto.SetSize (np); - equalto = 0; - - for (i = 1; i <= np; i++) - { - for (j = i+1; j <= np; j++) - { - if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12)) - equalto[j-1] = i; - } - } - - for (i = 1; i <= np; i++) - if (equalto[i-1]) - { - cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl; - for (j = 1; j <= mesh->GetNSeg(); j++) - { - Segment & seg = mesh->LineSegment(j); - if (seg[0] == i) seg[0] = equalto[i-1]; - if (seg[1] == i) seg[1] = equalto[i-1]; - } - } - - cout << "Removing degenerated segments" << endl; - for (j = 1; j <= mesh->GetNSeg(); j++) - { - Segment & seg = mesh->LineSegment(j); - if (seg[0] == seg[1]) - { - mesh->DeleteSegment(j); - cout << "Deleting Segment " << j << endl; - } - } - - mesh->Compress(); - */ - - /* - for (int i = 1; i <= geom.fmap.Extent(); i++) - { - Handle(Geom_Surface) hf1 = - BRep_Tool::Surface(TopoDS::Face(geom.fmap(i))); - for (int j = i+1; j <= geom.fmap.Extent(); j++) - { - Handle(Geom_Surface) hf2 = - BRep_Tool::Surface(TopoDS::Face(geom.fmap(j))); - if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl; - } - } - */ - -#ifdef LOG_STREAM - (*logout) << "Edges meshed" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; -#endif - } - - if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES) - return TCL_OK; - - if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE) - { - OCCMeshSurface (geom, *mesh, mparam.perfstepsend); - if (multithread.terminate) return TCL_OK; - -#ifdef LOG_STREAM - (*logout) << "Surfaces meshed" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; -#endif - -#ifdef STAT_STREAM - (*statout) << mesh->GetNSeg() << " & " - << mesh->GetNSE() << " & - &" - << GetTime() << " & " << endl; -#endif - - // MeshQuality2d (*mesh); - mesh->CalcSurfacesOfNode(); - } - - if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE) - return TCL_OK; - - if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME) - { - multithread.task = "Volume meshing"; - - MESHING3_RESULT res = MeshVolume (mparam, *mesh); - -/* - ofstream problemfile("occmesh.rep",ios_base::app); - - problemfile << "VOLUMEMESHING" << endl << endl; - if(res != MESHING3_OK) - problemfile << "ERROR" << endl << endl; - else - problemfile << "OK" << endl - << mesh->GetNE() << " elements" << endl << endl; - - problemfile.close(); -*/ - - if (res != MESHING3_OK) return TCL_ERROR; - - if (multithread.terminate) return TCL_OK; - - RemoveIllegalElements (*mesh); - if (multithread.terminate) return TCL_OK; - - MeshQuality3d (*mesh); - -#ifdef STAT_STREAM - (*statout) << GetTime() << " & "; -#endif - -#ifdef LOG_STREAM - (*logout) << "Volume meshed" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; -#endif - } - - if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME) - return TCL_OK; - - if (mparam.perfstepsstart <= MESHCONST_OPTVOLUME) - { - multithread.task = "Volume optimization"; - - OptimizeVolume (mparam, *mesh); - if (multithread.terminate) return TCL_OK; - -#ifdef STAT_STREAM - (*statout) << GetTime() << " & " - << mesh->GetNE() << " & " - << mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl; -#endif - -#ifdef LOG_STREAM - (*logout) << "Volume optimized" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; -#endif - - // cout << "Optimization complete" << endl; - - } - - (*testout) << "NP: " << mesh->GetNP() << endl; - for (int i = 1; i <= mesh->GetNP(); i++) - (*testout) << mesh->Point(i) << endl; - - (*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl; - for (int i = 1; i <= mesh->GetNSeg(); i++) - (*testout) << mesh->LineSegment(i) << endl; - - for (int i = 0; i < mesh->GetNDomains(); i++) - if(geom.snames.Size()) - mesh->SetMaterial( i+1, geom.snames[i] ); + if (multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE) return TCL_OK; - } + + if (mparam.perfstepsstart <= MESHCONST_MESHEDGES) + { + OCCFindEdges (geom, *mesh, mparam); + + /* + cout << "Removing redundant points" << endl; + + int i, j; + int np = mesh->GetNP(); + NgArray equalto; + + equalto.SetSize (np); + equalto = 0; + + for (i = 1; i <= np; i++) + { + for (j = i+1; j <= np; j++) + { + if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12)) + equalto[j-1] = i; + } + } + + for (i = 1; i <= np; i++) + if (equalto[i-1]) + { + cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl; + for (j = 1; j <= mesh->GetNSeg(); j++) + { + Segment & seg = mesh->LineSegment(j); + if (seg[0] == i) seg[0] = equalto[i-1]; + if (seg[1] == i) seg[1] = equalto[i-1]; + } + } + + cout << "Removing degenerated segments" << endl; + for (j = 1; j <= mesh->GetNSeg(); j++) + { + Segment & seg = mesh->LineSegment(j); + if (seg[0] == seg[1]) + { + mesh->DeleteSegment(j); + cout << "Deleting Segment " << j << endl; + } + } + + mesh->Compress(); + */ + + /* + for (int i = 1; i <= geom.fmap.Extent(); i++) + { + Handle(Geom_Surface) hf1 = + BRep_Tool::Surface(TopoDS::Face(geom.fmap(i))); + for (int j = i+1; j <= geom.fmap.Extent(); j++) + { + Handle(Geom_Surface) hf2 = + BRep_Tool::Surface(TopoDS::Face(geom.fmap(j))); + if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl; + } + } + */ + +#ifdef LOG_STREAM + (*logout) << "Edges meshed" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; +#endif + } + + if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES) + return TCL_OK; + + if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE) + { + OCCMeshSurface (geom, *mesh, mparam.perfstepsend, mparam); + if (multithread.terminate) return TCL_OK; + +#ifdef LOG_STREAM + (*logout) << "Surfaces meshed" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; +#endif + +#ifdef STAT_STREAM + (*statout) << mesh->GetNSeg() << " & " + << mesh->GetNSE() << " & - &" + << GetTime() << " & " << endl; +#endif + + // MeshQuality2d (*mesh); + mesh->CalcSurfacesOfNode(); + } + + if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE) + return TCL_OK; + + if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME) + { + multithread.task = "Volume meshing"; + + MESHING3_RESULT res = MeshVolume (mparam, *mesh); + + /* + ofstream problemfile("occmesh.rep",ios_base::app); + + problemfile << "VOLUMEMESHING" << endl << endl; + if(res != MESHING3_OK) + problemfile << "ERROR" << endl << endl; + else + problemfile << "OK" << endl + << mesh->GetNE() << " elements" << endl << endl; + + problemfile.close(); + */ + + if (res != MESHING3_OK) return TCL_ERROR; + + if (multithread.terminate) return TCL_OK; + + RemoveIllegalElements (*mesh); + if (multithread.terminate) return TCL_OK; + + MeshQuality3d (*mesh); + +#ifdef STAT_STREAM + (*statout) << GetTime() << " & "; +#endif + +#ifdef LOG_STREAM + (*logout) << "Volume meshed" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; +#endif + } + + if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHVOLUME) + return TCL_OK; + + if (mparam.perfstepsstart <= MESHCONST_OPTVOLUME) + { + multithread.task = "Volume optimization"; + + OptimizeVolume (mparam, *mesh); + if (multithread.terminate) return TCL_OK; + +#ifdef STAT_STREAM + (*statout) << GetTime() << " & " + << mesh->GetNE() << " & " + << mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl; +#endif + +#ifdef LOG_STREAM + (*logout) << "Volume optimized" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; +#endif + + // cout << "Optimization complete" << endl; + + } + + (*testout) << "NP: " << mesh->GetNP() << endl; + for (int i = 1; i <= mesh->GetNP(); i++) + (*testout) << mesh->Point(i) << endl; + + (*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl; + for (int i = 1; i <= mesh->GetNSeg(); i++) + (*testout) << mesh->LineSegment(i) << endl; + + for (int i = 0; i < mesh->GetNDomains(); i++) + if(geom.snames.Size()) + mesh->SetMaterial( i+1, geom.snames[i] ); + return TCL_OK; + } } #endif diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 2191aaaf..2209a2b5 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -166,7 +166,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a for (exp0.Init(shape, TopAbs_COMPSOLID); exp0.More(); exp0.Next()) nrcs++; double surfacecont = 0; - + { Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; rebuild->Apply(shape); @@ -869,7 +869,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a // Philippose - 15/01/2009 face_maxh.DeleteAll(); face_maxh.SetSize (fmap.Extent()); - face_maxh = mparam.maxh; + face_maxh = 1e99; // mparam.maxh; // Philippose - 15/01/2010 face_maxh_modified.DeleteAll(); diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index f461effd..fc2c68b9 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -114,7 +114,7 @@ namespace netgen { #include "occmeshsurf.hpp" - extern DLL_HEADER MeshingParameters mparam; + // extern DLL_HEADER MeshingParameters mparam; #define PROJECTION_TOLERANCE 1e-10 @@ -128,327 +128,323 @@ namespace netgen - class EntityVisualizationCode - { - int code; + class EntityVisualizationCode + { + int code; - public: + public: - EntityVisualizationCode() - { code = ENTITYISVISIBLE + !ENTITYISHIGHLIGHTED + ENTITYISDRAWABLE;} + EntityVisualizationCode() + { code = ENTITYISVISIBLE + !ENTITYISHIGHLIGHTED + ENTITYISDRAWABLE;} - int IsVisible () - { return code & ENTITYISVISIBLE;} + int IsVisible () + { return code & ENTITYISVISIBLE;} - int IsHighlighted () - { return code & ENTITYISHIGHLIGHTED;} + int IsHighlighted () + { return code & ENTITYISHIGHLIGHTED;} - int IsDrawable () - { return code & ENTITYISDRAWABLE;} + int IsDrawable () + { return code & ENTITYISDRAWABLE;} - void Show () - { code |= ENTITYISVISIBLE;} + void Show () + { code |= ENTITYISVISIBLE;} - void Hide () - { code &= ~ENTITYISVISIBLE;} + void Hide () + { code &= ~ENTITYISVISIBLE;} - void Highlight () - { code |= ENTITYISHIGHLIGHTED;} + void Highlight () + { code |= ENTITYISHIGHLIGHTED;} - void Lowlight () - { code &= ~ENTITYISHIGHLIGHTED;} + void Lowlight () + { code &= ~ENTITYISHIGHLIGHTED;} - void SetDrawable () - { code |= ENTITYISDRAWABLE;} + void SetDrawable () + { code |= ENTITYISDRAWABLE;} - void SetNotDrawable () - { code &= ~ENTITYISDRAWABLE;} - }; + void SetNotDrawable () + { code &= ~ENTITYISDRAWABLE;} + }; - class Line - { - public: - Point<3> p0, p1; + class Line + { + public: + Point<3> p0, p1; + double Dist (Line l); + double Length () { return (p1-p0).Length(); } + }; + - double Dist (Line l); - double Length (); - }; + inline double Det3 (double a00, double a01, double a02, + double a10, double a11, double a12, + double a20, double a21, double a22) + { + return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00; + } + - inline double Det3 (double a00, double a01, double a02, - double a10, double a11, double a12, - double a20, double a21, double a22) - { - return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00; - } + class OCCGeometry : public NetgenGeometry + { + Point<3> center; + public: + TopoDS_Shape shape; + TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; + NgArray fsingular, esingular, vsingular; + Box<3> boundingbox; + NgArray fnames, enames, snames; + // Philippose - 29/01/2009 + // OpenCascade XDE Support + // XCAF Handle to make the face colours available to the rest of + // the system + Handle_XCAFDoc_ColorTool face_colours; + + mutable int changed; + NgArray facemeshstatus; + // Philippose - 15/01/2009 + // Maximum mesh size for a given face + // (Used to explicitly define mesh size limits on individual faces) + NgArray face_maxh; + + // Philippose - 14/01/2010 + // Boolean array to detect whether a face has been explicitly modified + // by the user or not + NgArray face_maxh_modified; + + // Philippose - 15/01/2009 + // Indicates which faces have been selected by the user in geometry mode + // (Currently handles only selection of one face at a time, but an array would + // help to extend this to multiple faces) + NgArray face_sel_status; + + NgArray fvispar, evispar, vvispar; + + double tolerance; + bool fixsmalledges; + bool fixspotstripfaces; + bool sewfaces; + bool makesolids; + bool splitpartitions; + + OCCGeometry() + { + somap.Clear(); + shmap.Clear(); + fmap.Clear(); + wmap.Clear(); + emap.Clear(); + vmap.Clear(); + } + + + DLL_HEADER virtual void Save (string filename) const; + + void DoArchive(Archive& ar); + + DLL_HEADER void BuildFMap(); + + Box<3> GetBoundingBox() const + { return boundingbox; } + int NrSolids() const + { return somap.Extent(); } - class OCCGeometry : public NetgenGeometry - { - Point<3> center; + // Philippose - 17/01/2009 + // Total number of faces in the geometry + int NrFaces() const + { return fmap.Extent(); } - public: - TopoDS_Shape shape; - TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; - NgArray fsingular, esingular, vsingular; - Box<3> boundingbox; - NgArray fnames, enames, snames; - // Philippose - 29/01/2009 - // OpenCascade XDE Support - // XCAF Handle to make the face colours available to the rest of - // the system - Handle_XCAFDoc_ColorTool face_colours; + void SetCenter() + { center = boundingbox.Center(); } - mutable int changed; - NgArray facemeshstatus; + Point<3> Center() const + { return center; } - // Philippose - 15/01/2009 - // Maximum mesh size for a given face - // (Used to explicitly define mesh size limits on individual faces) - NgArray face_maxh; - - // Philippose - 14/01/2010 - // Boolean array to detect whether a face has been explicitly modified - // by the user or not - NgArray face_maxh_modified; + void Project (int surfi, Point<3> & p) const; + bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; - // Philippose - 15/01/2009 - // Indicates which faces have been selected by the user in geometry mode - // (Currently handles only selection of one face at a time, but an array would - // help to extend this to multiple faces) - NgArray face_sel_status; + OCCSurface GetSurface (int surfi) + { + cout << "OCCGeometry::GetSurface using PLANESPACE" << endl; + return OCCSurface (TopoDS::Face(fmap(surfi)), PLANESPACE); + } - NgArray fvispar, evispar, vvispar; + DLL_HEADER void CalcBoundingBox (); + DLL_HEADER void BuildVisualizationMesh (double deflection); + + void RecursiveTopologyTree (const TopoDS_Shape & sh, + stringstream & str, + TopAbs_ShapeEnum l, + bool free, + const char * lname); - double tolerance; - bool fixsmalledges; - bool fixspotstripfaces; - bool sewfaces; - bool makesolids; - bool splitpartitions; + DLL_HEADER void GetTopologyTree (stringstream & str); - OCCGeometry() - { - somap.Clear(); - shmap.Clear(); - fmap.Clear(); - wmap.Clear(); - emap.Clear(); - vmap.Clear(); - } + DLL_HEADER void PrintNrShapes (); + DLL_HEADER void CheckIrregularEntities (stringstream & str); - DLL_HEADER virtual void Save (string filename) const; + DLL_HEADER void SewFaces(); - void DoArchive(Archive& ar); + DLL_HEADER void MakeSolid(); - DLL_HEADER void BuildFMap(); + DLL_HEADER void HealGeometry(); - Box<3> GetBoundingBox() - { return boundingbox;} - - int NrSolids() - { return somap.Extent();} - - // Philippose - 17/01/2009 - // Total number of faces in the geometry - int NrFaces() - { return fmap.Extent();} - - void SetCenter() - { center = boundingbox.Center();} - - Point<3> Center() - { return center;} - - void Project (int surfi, Point<3> & p) const; - bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; - - OCCSurface GetSurface (int surfi) - { - cout << "OCCGeometry::GetSurface using PLANESPACE" << endl; - return OCCSurface (TopoDS::Face(fmap(surfi)), PLANESPACE); - } - - DLL_HEADER void CalcBoundingBox (); - DLL_HEADER void BuildVisualizationMesh (double deflection); - - void RecursiveTopologyTree (const TopoDS_Shape & sh, - stringstream & str, - TopAbs_ShapeEnum l, - bool free, - const char * lname); - - DLL_HEADER void GetTopologyTree (stringstream & str); - - DLL_HEADER void PrintNrShapes (); - - DLL_HEADER void CheckIrregularEntities (stringstream & str); - - DLL_HEADER void SewFaces(); - - DLL_HEADER void MakeSolid(); - - DLL_HEADER void HealGeometry(); - - // Philippose - 15/01/2009 - // Sets the maximum mesh size for a given face - // (Note: Local mesh size limited by the global max mesh size) - void SetFaceMaxH(int facenr, double faceh) - { - if((facenr> 0) && (facenr <= fmap.Extent())) - { - face_maxh[facenr-1] = min(mparam.maxh,faceh); + // Philippose - 15/01/2009 + // Sets the maximum mesh size for a given face + // (Note: Local mesh size limited by the global max mesh size) + void SetFaceMaxH(int facenr, double faceh, const MeshingParameters & mparam) + { + if((facenr> 0) && (facenr <= fmap.Extent())) + { + face_maxh[facenr-1] = min(mparam.maxh,faceh); - // Philippose - 14/01/2010 - // If the face maxh is greater than or equal to the - // current global maximum, then identify the face as - // not explicitly controlled by the user any more - if(faceh >= mparam.maxh) + // Philippose - 14/01/2010 + // If the face maxh is greater than or equal to the + // current global maximum, then identify the face as + // not explicitly controlled by the user any more + if(faceh >= mparam.maxh) { - face_maxh_modified[facenr-1] = 0; + face_maxh_modified[facenr-1] = 0; } - else + else { - face_maxh_modified[facenr-1] = 1; + face_maxh_modified[facenr-1] = 1; } - } - } + } + } - // Philippose - 15/01/2009 - // Returns the local mesh size of a given face - double GetFaceMaxH(int facenr) - { - if((facenr> 0) && (facenr <= fmap.Extent())) - { - return face_maxh[facenr-1]; - } - else - { - return 0.0; - } - } + // Philippose - 15/01/2009 + // Returns the local mesh size of a given face + double GetFaceMaxH(int facenr) + { + if((facenr> 0) && (facenr <= fmap.Extent())) + { + return face_maxh[facenr-1]; + } + else + { + return 0.0; + } + } - // Philippose - 14/01/2010 - // Returns the flag whether the given face - // has a mesh size controlled by the user or not - bool GetFaceMaxhModified(int facenr) - { - return face_maxh_modified[facenr-1]; - } + // Philippose - 14/01/2010 + // Returns the flag whether the given face + // has a mesh size controlled by the user or not + bool GetFaceMaxhModified(int facenr) + { + return face_maxh_modified[facenr-1]; + } - // Philippose - 17/01/2009 - // Returns the index of the currently selected face - int SelectedFace() - { - int i; - - for(i = 1; i <= fmap.Extent(); i++) - { - if(face_sel_status[i-1]) + // Philippose - 17/01/2009 + // Returns the index of the currently selected face + int SelectedFace() + { + for(int i = 1; i <= fmap.Extent(); i++) + { + if(face_sel_status[i-1]) { - return i; + return i; } - } + } - return 0; - } + return 0; + } - // Philippose - 17/01/2009 - // Sets the currently selected face - void SetSelectedFace(int facenr) - { - face_sel_status = 0; + // Philippose - 17/01/2009 + // Sets the currently selected face + void SetSelectedFace(int facenr) + { + face_sel_status = 0; - if((facenr >= 1) && (facenr <= fmap.Extent())) - { - face_sel_status[facenr-1] = 1; - } - } + if((facenr >= 1) && (facenr <= fmap.Extent())) + { + face_sel_status[facenr-1] = 1; + } + } - void LowLightAll() - { - for (int i = 1; i <= fmap.Extent(); i++) - fvispar[i-1].Lowlight(); - for (int i = 1; i <= emap.Extent(); i++) - evispar[i-1].Lowlight(); - for (int i = 1; i <= vmap.Extent(); i++) - vvispar[i-1].Lowlight(); - } + void LowLightAll() + { + for (int i = 1; i <= fmap.Extent(); i++) + fvispar[i-1].Lowlight(); + for (int i = 1; i <= emap.Extent(); i++) + evispar[i-1].Lowlight(); + for (int i = 1; i <= vmap.Extent(); i++) + vvispar[i-1].Lowlight(); + } - DLL_HEADER void GetUnmeshedFaceInfo (stringstream & str); - DLL_HEADER void GetNotDrawableFaces (stringstream & str); - DLL_HEADER bool ErrorInSurfaceMeshing (); + DLL_HEADER void GetUnmeshedFaceInfo (stringstream & str); + DLL_HEADER void GetNotDrawableFaces (stringstream & str); + DLL_HEADER bool ErrorInSurfaceMeshing (); -// void WriteOCC_STL(char * filename); + // void WriteOCC_STL(char * filename); - DLL_HEADER virtual int GenerateMesh (shared_ptr & mesh, MeshingParameters & mparam); + DLL_HEADER virtual int GenerateMesh (shared_ptr & mesh, MeshingParameters & mparam); - DLL_HEADER virtual const Refinement & GetRefinement () const; - }; + DLL_HEADER virtual const Refinement & GetRefinement () const; + }; - class DLL_HEADER OCCParameters - { - public: + class DLL_HEADER OCCParameters + { + public: - /// Factor for meshing close edges - double resthcloseedgefac; + /// Factor for meshing close edges + double resthcloseedgefac; - /// Enable / Disable detection of close edges - int resthcloseedgeenable; + /// Enable / Disable detection of close edges + int resthcloseedgeenable; - /// Minimum edge length to be used for dividing edges to mesh points - double resthminedgelen; + /// Minimum edge length to be used for dividing edges to mesh points + double resthminedgelen; - /// Enable / Disable use of the minimum edge length (by default use 1e-4) - int resthminedgelenenable; + /// Enable / Disable use of the minimum edge length (by default use 1e-4) + int resthminedgelenenable; - /*! - Default Constructor for the OpenCascade - Mesh generation parameter set - */ - OCCParameters(); + /*! + Default Constructor for the OpenCascade + Mesh generation parameter set + */ + OCCParameters(); - /*! - Dump all the OpenCascade specific meshing parameters - to console - */ - void Print (ostream & ost) const; - }; + /*! + Dump all the OpenCascade specific meshing parameters + to console + */ + void Print (ostream & ost) const; + }; - void PrintContents (OCCGeometry * geom); + void PrintContents (OCCGeometry * geom); - DLL_HEADER OCCGeometry * LoadOCC_IGES (const char * filename); - DLL_HEADER OCCGeometry * LoadOCC_STEP (const char * filename); - DLL_HEADER OCCGeometry * LoadOCC_BREP (const char * filename); + DLL_HEADER OCCGeometry * LoadOCC_IGES (const char * filename); + DLL_HEADER OCCGeometry * LoadOCC_STEP (const char * filename); + DLL_HEADER OCCGeometry * LoadOCC_BREP (const char * filename); - DLL_HEADER extern OCCParameters occparam; + DLL_HEADER extern OCCParameters occparam; - // Philippose - 31.09.2009 - // External access to the mesh generation functions within the OCC - // subsystem (Not sure if this is the best way to implement this....!!) - DLL_HEADER extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr & mesh, - MeshingParameters & mparam); + // Philippose - 31.09.2009 + // External access to the mesh generation functions within the OCC + // subsystem (Not sure if this is the best way to implement this....!!) + DLL_HEADER extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr & mesh, + MeshingParameters & mparam); - DLL_HEADER extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); + DLL_HEADER extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam); - DLL_HEADER extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend); + DLL_HEADER extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend, MeshingParameters & mparam); - DLL_HEADER extern void OCCFindEdges (OCCGeometry & geom, Mesh & mesh); + DLL_HEADER extern void OCCFindEdges (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam); } #endif diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp index 56fbf07b..b1d65fdd 100644 --- a/libsrc/occ/occmeshsurf.cpp +++ b/libsrc/occ/occmeshsurf.cpp @@ -173,6 +173,53 @@ namespace netgen gp_Vec du, dv; occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv); + // static Timer t("occ-defintangplane calculations"); + // RegionTimer reg(t); + + Mat<3,2> D1_; + D1_(0,0) = du.X(); D1_(1,0) = du.Y(); D1_(2,0) = du.Z(); + D1_(0,1) = dv.X(); D1_(1,1) = dv.Y(); D1_(2,1) = dv.Z(); + auto D1T_ = Trans(D1_); + auto D1TD1_ = D1T_*D1_; + if (Det (D1TD1_) == 0) throw SingularMatrixException(); + Mat<2,2> DDTinv_; + CalcInverse (D1TD1_, DDTinv_); + + Mat<3,2> Y_; + Vec<3> y1_ = (ap2-ap1).Normalize(); + Vec<3> y2_ = Cross(n, y1_).Normalize(); + for (int i = 0; i < 3; i++) + { + Y_(i,0) = y1_(i); + Y_(i,1) = y2_(i); + } + + auto A_ = DDTinv_ * D1T_ * Y_; + Mat<2,2> Ainv_; + if (Det(A_) == 0) throw SingularMatrixException(); + CalcInverse (A_, Ainv_); + + Vec<2> temp_ = Ainv_ * (psp2-psp1); + double r_ = temp_.Length(); + Mat<2,2> R_; + R_(0,0) = temp_(0)/r_; + R_(1,0) = temp_(1)/r_; + R_(0,1) = -R_(1,0); + R_(1,1) = R_(0,0); + + Ainv_ = Trans(R_) * Ainv_; + + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + { + Amat(i,j) = A_(i,j); + Amatinv(i,j) = Ainv_(i,j); + } + + // temp = Amatinv * (psp2-psp1); + + +#ifdef OLD DenseMatrix D1(3,2), D1T(2,3), DDTinv(2,2); D1(0,0) = du.X(); D1(1,0) = du.Y(); D1(2,0) = du.Z(); D1(0,1) = dv.X(); D1(1,1) = dv.Y(); D1(2,1) = dv.Z(); @@ -190,6 +237,8 @@ namespace netgen if (D1TD1.Det() == 0) throw SingularMatrixException(); CalcInverse (D1TD1, DDTinv); + // cout << " =?= inv = " << DDTinv << endl; + DenseMatrix Y(3,2); Vec<3> y1 = (ap2-ap1).Normalize(); Vec<3> y2 = Cross(n, y1).Normalize(); @@ -226,6 +275,7 @@ namespace netgen R(0,1) = sin (alpha); R(1,1) = cos (alpha); + // cout << "=?= R = " << R << endl; A = A*R; @@ -240,9 +290,10 @@ namespace netgen Amat(i,j) = A(i,j); Amatinv(i,j) = Ainv(i,j); } - + // cout << "=?= Ainv = " << endl << Ainv << endl; temp = Amatinv * (psp2-psp1); - + cout << " =?= Amatinv = " << Amatinv << endl; +#endif }; } @@ -376,7 +427,8 @@ namespace netgen Meshing2OCCSurfaces :: Meshing2OCCSurfaces (const TopoDS_Shape & asurf, - const Box<3> & abb, int aprojecttype) + const Box<3> & abb, int aprojecttype, + const MeshingParameters & mparam) : Meshing2(mparam, Box<3>(abb.PMin(), abb.PMax())), surface(TopoDS::Face(asurf), aprojecttype) { ; diff --git a/libsrc/occ/occmeshsurf.hpp b/libsrc/occ/occmeshsurf.hpp index 37eb230b..d10cc41d 100644 --- a/libsrc/occ/occmeshsurf.hpp +++ b/libsrc/occ/occmeshsurf.hpp @@ -55,6 +55,7 @@ protected: public: OCCSurface (const TopoDS_Face & aface, int aprojecttype) { + static Timer t("occurface ctor"); RegionTimer r(t); topods_face = aface; occface = BRep_Tool::Surface(topods_face); orient = topods_face.Orientation(); @@ -112,7 +113,8 @@ class Meshing2OCCSurfaces : public Meshing2 public: /// - Meshing2OCCSurfaces (const TopoDS_Shape & asurf, const Box<3> & aboundingbox, int aprojecttype); + Meshing2OCCSurfaces (const TopoDS_Shape & asurf, const Box<3> & aboundingbox, + int aprojecttype, const MeshingParameters & mparam); /// int GetProjectionType () diff --git a/libsrc/occ/occpkg.cpp b/libsrc/occ/occpkg.cpp index a5f64dcc..51673451 100644 --- a/libsrc/occ/occpkg.cpp +++ b/libsrc/occ/occpkg.cpp @@ -27,6 +27,7 @@ namespace netgen { extern DLL_HEADER shared_ptr ng_geometry; extern DLL_HEADER shared_ptr mesh; + extern DLL_HEADER MeshingParameters mparam; char * err_needsoccgeometry = (char*) "This operation needs an OCC geometry"; extern char * err_needsmesh; @@ -588,7 +589,7 @@ namespace netgen { if(!occgeometry->GetFaceMaxhModified(i)) { - occgeometry->SetFaceMaxH(i, mparam.maxh); + occgeometry->SetFaceMaxH(i, mparam.maxh, mparam); } } @@ -597,7 +598,7 @@ namespace netgen int facenr = atoi (argv[2]); double surfms = atof (argv[3]); if (occgeometry && facenr >= 1 && facenr <= occgeometry->NrFaces()) - occgeometry->SetFaceMaxH(facenr, surfms); + occgeometry->SetFaceMaxH(facenr, surfms, mparam); } @@ -608,7 +609,7 @@ namespace netgen { int nrFaces = occgeometry->NrFaces(); for (int i = 1; i <= nrFaces; i++) - occgeometry->SetFaceMaxH(i, surfms); + occgeometry->SetFaceMaxH(i, surfms, mparam); } } diff --git a/libsrc/occ/python_occ.cpp b/libsrc/occ/python_occ.cpp index 18514c8e..311fd78c 100644 --- a/libsrc/occ/python_occ.cpp +++ b/libsrc/occ/python_occ.cpp @@ -3,6 +3,7 @@ #include <../general/ngpython.hpp> #include +#include "../meshing/python_mesh.hpp" #include #include @@ -19,6 +20,21 @@ DLL_HEADER void ExportNgOCC(py::module &m) { py::class_, NetgenGeometry> (m, "OCCGeometry", R"raw_string(Use LoadOCCGeometry to load the geometry from a *.step file.)raw_string") .def(py::init<>()) + .def(py::init([] (const string& filename) + { + shared_ptr geo; + if(EndsWith(filename, ".step") || EndsWith(filename, ".stp")) + geo.reset(LoadOCC_STEP(filename.c_str())); + else if(EndsWith(filename, ".brep")) + geo.reset(LoadOCC_BREP(filename.c_str())); + else if(EndsWith(filename, ".iges")) + geo.reset(LoadOCC_IGES(filename.c_str())); + else + throw Exception("Cannot load file " + filename + "\nValid formats are: step, stp, brep, iges"); + ng_geometry = geo; + return geo; + }), py::arg("filename"), + "Load OCC geometry from step, brep or iges file") .def(NGSPickle()) .def("Heal",[](OCCGeometry & self, double tolerance, bool fixsmalledges, bool fixspotstripfaces, bool sewfaces, bool makesolids, bool splitpartitions) { @@ -108,34 +124,35 @@ DLL_HEADER void ExportNgOCC(py::module &m) res["max"] = MoveToNumpy(max); return res; }, py::call_guard()) - ; - m.def("LoadOCCGeometry",FunctionPointer([] (const string & filename) + .def("GenerateMesh", [](shared_ptr geo, + MeshingParameters* pars, py::kwargs kwargs) + { + MeshingParameters mp; + if(pars) mp = *pars; + { + py::gil_scoped_acquire aq; + CreateMPfromKwargs(mp, kwargs); + } + auto mesh = make_shared(); + SetGlobalMesh(mesh); + mesh->SetGeometry(geo); + ng_geometry = geo; + geo->GenerateMesh(mesh,mp); + return mesh; + }, py::arg("mp") = nullptr, + py::call_guard(), + meshingparameter_description.c_str()) + ; + + m.def("LoadOCCGeometry",[] (const string & filename) { - cout << "load OCC geometry"; + cout << "WARNING: LoadOCCGeometry is deprecated! Just use the OCCGeometry(filename) constructor. It is able to read brep and iges files as well!" << endl; ifstream ist(filename); OCCGeometry * instance = new OCCGeometry(); instance = LoadOCC_STEP(filename.c_str()); ng_geometry = shared_ptr(instance, NOOP_Deleter); return ng_geometry; - }),py::call_guard()); - m.def("GenerateMesh", FunctionPointer([] (shared_ptr geo, MeshingParameters ¶m) - { - auto mesh = make_shared(); - SetGlobalMesh(mesh); - mesh->SetGeometry(geo); - ng_geometry = geo; - - try - { - geo->GenerateMesh(mesh,param); - } - catch (NgException ex) - { - cout << "Caught NgException: " << ex.What() << endl; - } - return mesh; - }),py::call_guard()) - ; + },py::call_guard()); } PYBIND11_MODULE(libNgOCC, m) { diff --git a/libsrc/stlgeom/python_stl.cpp b/libsrc/stlgeom/python_stl.cpp index 0d2bb32f..39a52273 100644 --- a/libsrc/stlgeom/python_stl.cpp +++ b/libsrc/stlgeom/python_stl.cpp @@ -4,6 +4,7 @@ #include <../general/ngpython.hpp> #include #include +#include "../meshing/python_mesh.hpp" #ifdef WIN32 #define DLL_HEADER __declspec(dllexport) @@ -21,6 +22,12 @@ DLL_HEADER void ExportSTL(py::module & m) { py::class_, NetgenGeometry> (m,"STLGeometry") .def(py::init<>()) + .def(py::init<>([](const string& filename) + { + ifstream ist(filename); + return shared_ptr(STLGeometry::Load(ist)); + }), py::arg("filename"), + py::call_guard()) .def(NGSPickle()) .def("_visualizationData", [](shared_ptr stl_geo) { @@ -71,29 +78,31 @@ DLL_HEADER void ExportSTL(py::module & m) res["max"] = MoveToNumpy(max); return res; }, py::call_guard()) + .def("GenerateMesh", [] (shared_ptr geo, + MeshingParameters* pars, py::kwargs kwargs) + { + MeshingParameters mp; + if(pars) mp = *pars; + { + py::gil_scoped_acquire aq; + CreateMPfromKwargs(mp, kwargs); + } + auto mesh = make_shared(); + SetGlobalMesh(mesh); + mesh->SetGeometry(geo); + ng_geometry = geo; + geo->GenerateMesh(mesh,mp); + return mesh; + }, py::arg("mp") = nullptr, + py::call_guard(), + meshingparameter_description.c_str()) ; - m.def("LoadSTLGeometry", FunctionPointer([] (const string & filename) - { - ifstream ist(filename); - return shared_ptr(STLGeometry::Load(ist)); - }),py::call_guard()); - m.def("GenerateMesh", FunctionPointer([] (shared_ptr geo, MeshingParameters ¶m) - { - auto mesh = make_shared(); - SetGlobalMesh(mesh); - mesh->SetGeometry(geo); - ng_geometry = geo; - try - { - geo->GenerateMesh(mesh,param); - } - catch (NgException ex) - { - cout << "Caught NgException: " << ex.What() << endl; - } - return mesh; - }),py::call_guard()) - ; + m.def("LoadSTLGeometry", [] (const string & filename) + { + cout << "WARNING: LoadSTLGeometry is deprecated, use the STLGeometry(filename) constructor instead!" << endl; + ifstream ist(filename); + return shared_ptr(STLGeometry::Load(ist)); + },py::call_guard()); } PYBIND11_MODULE(libstl, m) { diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index 7bba75f5..a177b9bd 100644 --- a/nglib/nglib.cpp +++ b/nglib/nglib.cpp @@ -857,7 +857,7 @@ namespace nglib // slate me->DeleteMesh(); - OCCSetLocalMeshSize(*occgeom, *me); + OCCSetLocalMeshSize(*occgeom, *me, mparam); return(NG_OK); } @@ -875,7 +875,7 @@ namespace nglib mp->Transfer_Parameters(); - OCCFindEdges(*occgeom, *me); + OCCFindEdges(*occgeom, *me, mparam); if((me->GetNP()) && (me->GetNFD())) { @@ -920,7 +920,7 @@ namespace nglib perfstepsend = MESHCONST_OPTSURFACE; } - OCCMeshSurface(*occgeom, *me, perfstepsend); + OCCMeshSurface(*occgeom, *me, perfstepsend, mparam); me->CalcSurfacesOfNode(); diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 9f4bc635..477befea 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -7,7 +7,7 @@ configure_file(__init__.py ${CMAKE_CURRENT_BINARY_DIR}/__init__.py @ONLY) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/__init__.py - meshing.py csg.py geom2d.py stl.py gui.py NgOCC.py read_gmsh.py + meshing.py csg.py geom2d.py stl.py gui.py NgOCC.py occ.py read_gmsh.py DESTINATION ${NG_INSTALL_DIR_PYTHON}/${NG_INSTALL_SUFFIX} COMPONENT netgen ) diff --git a/python/NgOCC.py b/python/NgOCC.py index 40f0a51d..ebc7fd51 100644 --- a/python/NgOCC.py +++ b/python/NgOCC.py @@ -1,10 +1,7 @@ -from netgen.libngpy._NgOCC import * -from netgen.libngpy._meshing import MeshingParameters -def NgOCC_meshing_func (geom, **args): - if "mp" in args: - return GenerateMesh (geom, args["mp"]) - else: - return GenerateMesh (geom, MeshingParameters (**args)) +import logging +logger = logging.getLogger(__name__) -OCCGeometry.GenerateMesh = NgOCC_meshing_func +logger.warning("This module is deprecated and just a wrapper for netgen.occ, import netgen.occ instead") + +from .occ import * diff --git a/python/csg.py b/python/csg.py index f9333179..f6b37f1c 100644 --- a/python/csg.py +++ b/python/csg.py @@ -1,34 +1,19 @@ -from netgen.libngpy._csg import * -from netgen.libngpy._meshing import MeshingParameters -from netgen.libngpy._meshing import Pnt -from netgen.libngpy._meshing import Vec -from netgen.libngpy._meshing import Trafo - +from .libngpy._csg import * +from .libngpy._meshing import Pnt, Vec, Trafo +from .meshing import meshsize try: - import libngpy.csgvis as csgvis - from libngpy.csgvis import MouseMove + from . import csgvis + from .csgvis import MouseMove CSGeometry.VS = csgvis.VS SetBackGroundColor = csgvis.SetBackGroundColor del csgvis def VS (obj): return obj.VS() - except: pass - -def csg_meshing_func (geom, **args): - if "mp" in args: - return GenerateMesh (geom, args["mp"]) - else: - return GenerateMesh (geom, MeshingParameters (**args)) -# return GenerateMesh (geom, MeshingParameters (**args)) - -CSGeometry.GenerateMesh = csg_meshing_func - - unit_cube = CSGeometry() p1 = Plane(Pnt(0,0,0),Vec(-1,0,0)).bc("back") p2 = Plane(Pnt(1,1,1),Vec(1,0,0)).bc("front") @@ -37,5 +22,4 @@ p4 = Plane(Pnt(1,1,1),Vec(0,1,0)).bc("right") p5 = Plane(Pnt(0,0,0),Vec(0,0,-1)).bc("bottom") p6 = Plane(Pnt(1,1,1),Vec(0,0,1)).bc("top") unit_cube.Add (p1*p2*p3*p4*p5*p6, col=(0,0,1)) -# unit_cube.Add (OrthoBrick(Pnt(0,0,0), Pnt(1,1,1))) diff --git a/python/geom2d.py b/python/geom2d.py index 34b37870..8af508b1 100644 --- a/python/geom2d.py +++ b/python/geom2d.py @@ -1,24 +1,12 @@ -from netgen.libngpy._geom2d import * -from netgen.libngpy._meshing import * - -tmp_generate_mesh = SplineGeometry.GenerateMesh - -def geom2d_meshing_func (geom, **args): - if "mp" in args: - return tmp_generate_mesh (geom, args["mp"]) - else: - return tmp_generate_mesh (geom, MeshingParameters (**args)) - - -SplineGeometry.GenerateMesh = geom2d_meshing_func - +from .libngpy._geom2d import SplineGeometry +from .meshing import meshsize unit_square = SplineGeometry() -pnts = [ (0,0), (1,0), (1,1), (0,1) ] -lines = [ (0,1,1,"bottom"), (1,2,2,"right"), (2,3,3,"top"), (3,0,4,"left") ] -pnums = [unit_square.AppendPoint(*p) for p in pnts] -for l1,l2,bc,bcname in lines: - unit_square.Append( ["line", pnums[l1], pnums[l2]], bc=bcname) +_pnts = [ (0,0), (1,0), (1,1), (0,1) ] +_lines = [ (0,1,1,"bottom"), (1,2,2,"right"), (2,3,3,"top"), (3,0,4,"left") ] +_pnums = [unit_square.AppendPoint(*p) for p in _pnts] +for l1,l2,bc,bcname in _lines: + unit_square.Append( ["line", _pnums[l1], _pnums[l2]], bc=bcname) def MakeRectangle (geo, p1, p2, bc=None, bcs=None, **args): @@ -141,8 +129,4 @@ SplineGeometry.AddSegment = lambda *args, **kwargs : SplineGeometry.Append(*args SplineGeometry.AddPoint = lambda *args, **kwargs : SplineGeometry.AppendPoint(*args, **kwargs) SplineGeometry.CreatePML = CreatePML -__all__ = ['SplineGeometry', 'unit_square'] - - - diff --git a/python/gui.py b/python/gui.py index 35fa1f34..959b2fd3 100644 --- a/python/gui.py +++ b/python/gui.py @@ -19,4 +19,8 @@ def StartGUI(): pass if not netgen.libngpy._meshing._netgen_executable_started: - StartGUI() + # catch exception for building html docu on server without display + try: + StartGUI() + except: + pass diff --git a/python/meshing.py b/python/meshing.py index f6ad65ce..5fbebdf4 100644 --- a/python/meshing.py +++ b/python/meshing.py @@ -1 +1,22 @@ -from netgen.libngpy._meshing import * +from .libngpy._meshing import * + +class _MeshsizeObject: + pass + +meshsize = _MeshsizeObject() + +meshsize.very_coarse = MeshingParameters(curvaturesafety=1, + segmentsperedge=0.3, + grading=0.7) +meshsize.coarse = MeshingParameters(curvaturesafety=1.5, + segmentsperedge=0.5, + grading=0.5) +meshsize.moderate = MeshingParameters(curvaturesafety=2, + segmentsperedge=1, + grading=0.3) +meshsize.fine = MeshingParameters(curvaturesafety=3, + segmentsperedge=2, + grading=0.2) +meshsize.very_fine = MeshingParameters(curvaturesafety=5, + segmentsperedge=3, + grading=0.1) diff --git a/python/occ.py b/python/occ.py new file mode 100644 index 00000000..f10a8611 --- /dev/null +++ b/python/occ.py @@ -0,0 +1,2 @@ +from .libngpy._NgOCC import * +from .meshing import meshsize diff --git a/python/stl.py b/python/stl.py index 032f2b69..40f4421a 100644 --- a/python/stl.py +++ b/python/stl.py @@ -1,12 +1,2 @@ from netgen.libngpy._stl import * -from netgen.libngpy._meshing import MeshingParameters - - -def stl_meshing_func (geom, **args): - if "mp" in args: - return GenerateMesh (geom, args["mp"]) - else: - return GenerateMesh (geom, MeshingParameters (**args)) -# return GenerateMesh (geom, MeshingParameters (**args)) - -STLGeometry.GenerateMesh = stl_meshing_func +from .meshing import meshsize diff --git a/tests/pytest/test_pickling.py b/tests/pytest/test_pickling.py index 13e20f9e..17f22163 100644 --- a/tests/pytest/test_pickling.py +++ b/tests/pytest/test_pickling.py @@ -48,11 +48,11 @@ def test_pickle_stl(): def test_pickle_occ(): try: - import netgen.NgOCC as occ + import netgen.occ as occ except: import pytest pytest.skip("can't import occ") - geo = occ.LoadOCCGeometry("../../tutorials/frame.step") + geo = occ.OCCGeometry("../../tutorials/frame.step") geo_dump = pickle.dumps(geo) geo2 = pickle.loads(geo_dump) vd1 = geo._visualizationData() diff --git a/tests/pytest/test_savemesh.py b/tests/pytest/test_savemesh.py index 6b7b7e2c..77ccd198 100644 --- a/tests/pytest/test_savemesh.py +++ b/tests/pytest/test_savemesh.py @@ -28,7 +28,7 @@ def CreateGeo(): def test_BBNDsave(): mesh = CreateGeo().GenerateMesh(maxh=0.4,perfstepsend = meshing.MeshingStep.MESHSURFACE) for i in range(2): - mesh.GenerateVolumeMesh(mp = MeshingParameters(only3D_domain=i+1,maxh=0.4)) + mesh.GenerateVolumeMesh(only3D_domain=i+1,maxh=0.4) mesh.SetGeometry(None) mesh.Save("test.vol") mesh2 = meshing.Mesh()