ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-09-23 02:02:49 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master'

Browse Source
This commit is contained in:
Julius Zimmermann 2019-07-31 23:00:20 +02:00
commit 39b34a92cb
41 changed files with 1915 additions and 1663 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

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

16
libsrc/core/utils.hpp
View File

@ -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

34
libsrc/csg/python_csg.cpp
View File

@ -3,6 +3,7 @@
#include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include <csg.hpp>
#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<py::gil_scoped_release>())
;
m.def("GenerateMesh", FunctionPointer
([](shared_ptr<CSGeometry> geo, MeshingParameters & param)
.def("GenerateMesh", [](shared_ptr<CSGeometry> geo,
MeshingParameters* pars, py::kwargs kwargs)
{
auto dummy = make_shared<Mesh>();
SetGlobalMesh (dummy);
dummy->SetGeometry(geo);
MeshingParameters mp;
if(pars) mp = *pars;
{
py::gil_scoped_acquire aq;
CreateMPfromKwargs(mp, kwargs);
}
auto mesh = make_shared<Mesh>();
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<py::gil_scoped_release>())
geo->GenerateMesh (mesh, mp);
return mesh;
}, py::arg("mp") = nullptr,
meshingparameter_description.c_str(),
py::call_guard<py::gil_scoped_release>())
;
m.def("Save", FunctionPointer

5
libsrc/general/ngarray.hpp
View File

@ -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)
{

20
libsrc/geom2d/python_geom2d.cpp
View File

@ -2,6 +2,7 @@
#include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include "../meshing/python_mesh.hpp"
#include <meshing.hpp>
#include <geometry2d.hpp>
@ -362,16 +363,25 @@ 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<SplineGeometry2d> self, MeshingParameters & mparam)
// If we change to c++17 this can become optional<MeshingParameters>
.def("GenerateMesh", [](shared_ptr<SplineGeometry2d> self,
MeshingParameters* pars, py::kwargs kwargs)
{
shared_ptr<Mesh> mesh = make_shared<Mesh> ();
MeshingParameters mp;
if(pars) mp = *pars;
{
py::gil_scoped_acquire aq;
CreateMPfromKwargs(mp, kwargs);
}
auto mesh = make_shared<Mesh>();
mesh->SetGeometry(self);
SetGlobalMesh (mesh);
ng_geometry = self;
self->GenerateMesh(mesh, mparam);
self->GenerateMesh(mesh, mp);
return mesh;
},py::call_guard<py::gil_scoped_release>())
}, py::arg("mp") = nullptr,
py::call_guard<py::gil_scoped_release>(),
meshingparameter_description.c_str())
;
}

1
libsrc/include/nginterface_v2.hpp
View File

@ -284,6 +284,7 @@ namespace netgen
int GetDimension() const;
int GetNLevels() const;
size_t GetNVLevel (int level) const;
int GetNElements (int dim) const;
int GetNNodes (int nt) const;

4
libsrc/interface/nginterface.cpp
View File

@ -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;
}

14
libsrc/interface/nginterface_v2.cpp
View File

@ -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
@ -1156,9 +1164,7 @@ namespace netgen
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);

24
libsrc/meshing/bisect.cpp
View File

@ -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++)
@ -3986,7 +3991,8 @@ 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");

2
libsrc/meshing/curvedelems.hpp
View File

@ -49,7 +49,7 @@ public:
int GetOrder () { return order; }
virtual void DoArchive(Archive& ar)
void DoArchive(Archive& ar)
{
if(ar.Input())
buildJacPols();

52
libsrc/meshing/delaunay2d.cpp
View File

@ -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]]);

2
libsrc/meshing/meshclass.cpp
View File

@ -20,7 +20,7 @@ namespace netgen
segmentht = NULL;
lochfunc = NULL;
mglevels = 1;
// mglevels = 1;
elementsearchtree = NULL;
elementsearchtreets = NextTimeStamp();
majortimestamp = timestamp = NextTimeStamp();

7
libsrc/meshing/meshclass.hpp
View File

@ -175,7 +175,9 @@ namespace netgen
/// number of refinement levels
int mglevels;
// int mglevels;
// number of vertices on each refinement level:
NgArray<size_t> level_nv;
/// refinement hierarchy
NgArray<PointIndices<2>,PointIndex::BASE> mlbetweennodes;
/// parent element of volume element
@ -796,7 +798,8 @@ namespace netgen
shared_ptr<NetgenGeometry> GetGeometry() const
{
return geometry;
static auto global_geometry = make_shared<NetgenGeometry>();
return geometry ? geometry : global_geometry;
}
void SetGeometry (shared_ptr<NetgenGeometry> geom)
{

2
libsrc/meshing/meshfunc.cpp
View File

@ -690,7 +690,7 @@ namespace netgen
case 'j': mesh3d.ImproveMeshJacobian(mp); break;
}
}
mesh3d.mglevels = 1;
// mesh3d.mglevels = 1;
MeshQuality3d (mesh3d);
}

3
libsrc/meshing/meshfunc2d.cpp
View File

@ -63,10 +63,7 @@ namespace netgen
}
if (secondorder)
{
if (mesh.GetGeometry())
mesh.GetGeometry()->GetRefinement().MakeSecondOrder(mesh);
else
Refinement().MakeSecondOrder(mesh);
}
}

168
libsrc/meshing/meshing2.cpp
View File

@ -19,15 +19,33 @@ namespace netgen
// static int qualclass;
Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox)
{
boundingbox = aboundingbox;
static Array<unique_ptr<netrule>> global_trig_rules;
static Array<unique_ptr<netrule>> global_quad_rules;
Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox)
: adfront(aboundingbox), boundingbox(aboundingbox)
{
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<netrule>(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,8 +240,7 @@ namespace netgen
bool morerisc;
bool debugflag;
double h, his, hshould;
// double h;
NgArray<Point3d> locpoints;
NgArray<int> legalpoints;
@ -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;

2
libsrc/meshing/meshing2.hpp
View File

@ -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<netrule*> rules;
/// statistics

30
libsrc/meshing/meshtype.hpp
View File

@ -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 ();
///

3
libsrc/meshing/parser2.cpp
View File

@ -579,6 +579,8 @@ void Meshing2 :: LoadRules (const char * filename, bool quad)
exit (1);
}
Timer t("Parsing rules");
t.Start();
while (!ist->eof())
{
buf[0] = 0;
@ -597,6 +599,7 @@ void Meshing2 :: LoadRules (const char * filename, bool quad)
//(*testout) << "loop" << endl;
}
//(*testout) << "POS3" << endl;
t.Stop();
delete ist;
//delete [] tr1;

68
libsrc/meshing/python_mesh.cpp
View File

@ -2,6 +2,7 @@
#include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include "python_mesh.hpp"
#include <mystdlib.h>
#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<MeshingParameters>(pymp).check())
mp = py::extract<MeshingParameters>(pymp)();
else
{
mp.optsteps3d = 5;
}
CreateMPfromKwargs(mp, kwargs);
}
MeshVolume (mp, self);
OptimizeVolume (mp, self);
},
py::arg("mp")=NGDummyArgument(),py::call_guard<py::gil_scoped_release>())
}, py::arg("mp")=nullptr,
meshingparameter_description.c_str(),
py::call_guard<py::gil_scoped_release>())
.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.UpdateTopology();
}),py::call_guard<py::gil_scoped_release>())
.def ("SecondOrder", FunctionPointer
([](Mesh & self)
{
if (self.GetGeometry())
self.GetGeometry()->GetRefinement().MakeSecondOrder(self);
else
Refinement().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_<MP> (m, "MeshingParameters")
auto mp = py::class_<MP> (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<MP>)
.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)
{

168
libsrc/meshing/python_mesh.hpp Normal file
View File

@ -0,0 +1,168 @@
#include <pybind11/pybind11.h>
#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<string>(kwargs["optimize3d"]);
if(kwargs.contains("optsteps3d"))
mp.optsteps3d = py::cast<int>(kwargs["optsteps3d"]);
if(kwargs.contains("optimize2d"))
mp.optimize2d = py::cast<string>(kwargs["optimize2d"]);
if(kwargs.contains("optsteps2d"))
mp.optsteps2d = py::cast<int>(kwargs["optsteps2d"]);
if(kwargs.contains("opterrpow"))
mp.opterrpow = py::cast<double>(kwargs["opterrpow"]);
if(kwargs.contains("blockfill"))
mp.blockfill = py::cast<bool>(kwargs["blockfill"]);
if(kwargs.contains("filldist"))
mp.filldist = py::cast<double>(kwargs["filldist"]);
if(kwargs.contains("safety"))
mp.safety = py::cast<double>(kwargs["safety"]);
if(kwargs.contains("relinnersafety"))
mp.relinnersafety = py::cast<double>(kwargs["relinnersafety"]);
if(kwargs.contains("uselocalh"))
mp.uselocalh = py::cast<bool>(kwargs["uselocalh"]);
if(kwargs.contains("grading"))
mp.grading = py::cast<double>(kwargs["grading"]);
if(kwargs.contains("delaunay"))
mp.delaunay = py::cast<bool>(kwargs["delaunay"]);
if(kwargs.contains("maxh"))
mp.maxh = py::cast<double>(kwargs["maxh"]);
if(kwargs.contains("minh"))
mp.minh = py::cast<double>(kwargs["minh"]);
if(kwargs.contains("meshsizefilename"))
mp.meshsizefilename = py::cast<string>(kwargs["meshsizefilename"]);
if(kwargs.contains("startinsurface"))
mp.startinsurface = py::cast<bool>(kwargs["startinsurface"]);
if(kwargs.contains("checkoverlap"))
mp.checkoverlap = py::cast<bool>(kwargs["checkoverlap"]);
if(kwargs.contains("checkoverlappingboundary"))
mp.checkoverlappingboundary = py::cast<bool>(kwargs["checkoverlappingboundary"]);
if(kwargs.contains("checkchartboundary"))
mp.checkchartboundary = py::cast<bool>(kwargs["checkchartboundary"]);
if(kwargs.contains("curvaturesafety"))
mp.curvaturesafety = py::cast<double>(kwargs["curvaturesafety"]);
if(kwargs.contains("segmentsperedge"))
mp.segmentsperedge = py::cast<double>(kwargs["segmentsperedge"]);
if(kwargs.contains("parthread"))
mp.parthread = py::cast<bool>(kwargs["parthread"]);
if(kwargs.contains("elsizeweight"))
mp.elsizeweight = py::cast<double>(kwargs["elsizeweight"]);
if(kwargs.contains("perfstepsstart"))
mp.perfstepsstart = py::cast<int>(kwargs["perfstepsstart"]);
if(kwargs.contains("perfstepsend"))
mp.perfstepsend = py::cast<int>(kwargs["perfstepsend"]);
if(kwargs.contains("giveuptol2d"))
mp.giveuptol2d = py::cast<int>(kwargs["giveuptol2d"]);
if(kwargs.contains("giveuptol"))
mp.giveuptol = py::cast<int>(kwargs["giveuptol"]);
if(kwargs.contains("maxoutersteps"))
mp.maxoutersteps = py::cast<int>(kwargs["maxoutersteps"]);
if(kwargs.contains("starshapeclass"))
mp.starshapeclass = py::cast<int>(kwargs["starshapeclass"]);
if(kwargs.contains("baseelnp"))
mp.baseelnp = py::cast<int>(kwargs["baseelnp"]);
if(kwargs.contains("sloppy"))
mp.sloppy = py::cast<bool>(kwargs["sloppy"]);
if(kwargs.contains("badellimit"))
mp.badellimit = py::cast<double>(kwargs["badellimit"]);
if(kwargs.contains("check_impossible"))
mp.check_impossible = py::cast<bool>(kwargs["check_impossible"]);
if(kwargs.contains("only3D_domain_nr"))
mp.only3D_domain_nr = py::cast<int>(kwargs["only3D_domain_nr"]);
if(kwargs.contains("secondorder"))
mp.secondorder = py::cast<bool>(kwargs["secondorder"]);
if(kwargs.contains("elementorder"))
mp.elementorder = py::cast<int>(kwargs["elementorder"]);
if(kwargs.contains("quad"))
mp.quad = py::cast<bool>(kwargs["quad"]);
if(kwargs.contains("try_hexes"))
mp.try_hexes = py::cast<bool>(kwargs["try_hexes"]);
if(kwargs.contains("inverttets"))
mp.inverttets = py::cast<bool>(kwargs["inverttets"]);
if(kwargs.contains("inverttrigs"))
mp.inverttrigs = py::cast<bool>(kwargs["inverttrigs"]);
if(kwargs.contains("autozrefine"))
mp.autozrefine = py::cast<bool>(kwargs["autozrefine"]);
}
} // namespace netgen

56
libsrc/meshing/smoothing3.cpp
View File

@ -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<double,PointIndex::BASE> 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;
PointFunction pf(points, volelements, mp);
if (typ == 1)
pf = new PointFunction(points, volelements, mp);
else
pf = new CheapPointFunction(points, volelements, mp);
// pf->SetLocalH (h);
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);
}

134
libsrc/occ/occgenmesh.cpp
View File

@ -41,24 +41,14 @@ namespace netgen
return 1e99;
}
double Line :: Length ()
{
return (p1-p0).Length();
}
inline Point<3> occ2ng (const gp_Pnt & p)
{
return Point<3> (p.X(), p.Y(), p.Z());
}
double ComputeH (double kappa)
double ComputeH (double kappa, const MeshingParameters & mparam)
{
/*
double hret;
kappa *= mparam.curvaturesafety;
@ -70,14 +60,17 @@ namespace netgen
if (mparam.maxh < hret)
hret = mparam.maxh;
return (hret);
return hret;
*/
// return min(mparam.maxh, 1/kappa);
return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;
}
void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0)
BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h,
const MeshingParameters & mparam)
{
int ls = -1;
@ -166,7 +159,7 @@ namespace netgen
h = ComputeH (curvature+1e-10);
h = ComputeH (curvature+1e-10, mparam);
if(h < 1e-4*maxside)
return;
@ -188,20 +181,20 @@ namespace netgen
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);
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);
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);
}
}
@ -232,7 +225,8 @@ namespace netgen
void DivideEdge (TopoDS_Edge & edge, NgArray<MeshPoint> & ps,
NgArray<double> & params, Mesh & mesh)
NgArray<double> & params, Mesh & mesh,
const MeshingParameters & mparam)
{
double s0, s1;
double maxh = mparam.maxh;
@ -312,8 +306,10 @@ namespace netgen
void OCCFindEdges (OCCGeometry & geom, Mesh & mesh)
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";
@ -327,6 +323,7 @@ namespace netgen
double eps = 1e-6 * geom.GetBoundingBox().Diam();
tsearch.Start();
for (int i = 1; i <= nvertices; i++)
{
gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i)));
@ -344,6 +341,7 @@ namespace netgen
if (!exists)
mesh.AddPoint (mp);
}
tsearch.Stop();
(*testout) << "different vertices = " << mesh.GetNP() << endl;
@ -474,7 +472,7 @@ namespace netgen
NgArray <MeshPoint> mp;
NgArray <double> params;
DivideEdge (edge, mp, params, mesh);
DivideEdge (edge, mp, params, mesh, mparam);
NgArray <int> pnums;
pnums.SetSize (mp.Size()+2);
@ -502,6 +500,7 @@ namespace netgen
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)
@ -509,6 +508,7 @@ namespace netgen
exists = 1;
break;
}
tsearch.Stop();
if (exists)
pnums[i] = j;
@ -606,10 +606,13 @@ namespace netgen
void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend)
void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend,
MeshingParameters & mparam)
{
int i, j, k;
int changed;
static Timer t("OCCMeshSurface"); RegionTimer r(t);
// int i, j, k;
// int changed;
const char * savetask = multithread.task;
multithread.task = "Surface meshing";
@ -630,7 +633,7 @@ namespace netgen
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())
{
@ -652,7 +655,6 @@ namespace netgen
}
*/
FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
int oldnf = mesh.GetNSE();
@ -661,32 +663,38 @@ namespace netgen
// int projecttype = PLANESPACE;
Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype);
static Timer tinit("init");
tinit.Start();
Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
tinit.Stop();
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;
// Meshing2OCCSurfaces meshing(f2, bb);
meshing.SetStartTime (starttime);
//(*testout) << "Face " << k << endl << endl;
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++)
for (int i = 1; i <= mesh.GetNSeg(); i++)
{
Segment & seg = mesh.LineSegment(i);
if (seg.si == k)
{
for (j = 1; j <= 2; j++)
for (int j = 1; j <= 2; j++)
{
int pi = (j == 1) ? seg[0] : seg[1];
if (!glob2loc.Get(pi))
@ -699,7 +707,7 @@ namespace netgen
}
}
for (i = 1; i <= mesh.GetNSeg(); i++)
for (int i = 1; i <= mesh.GetNSeg(); i++)
{
Segment & seg = mesh.LineSegment(i);
if (seg.si == k)
@ -719,9 +727,11 @@ namespace netgen
}
else
{
static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t);
int cntp = 0;
for (i = 1; i <= mesh.GetNSeg(); i++)
for (int i = 1; i <= mesh.GetNSeg(); i++)
if (mesh.LineSegment(i).si == k)
cntp+=2;
@ -731,7 +741,7 @@ namespace netgen
gis.SetAllocSize (cntp);
gis.SetSize (0);
for (i = 1; i <= mesh.GetNSeg(); i++)
for (int i = 1; i <= mesh.GetNSeg(); i++)
{
Segment & seg = mesh.LineSegment(i);
if (seg.si == k)
@ -745,7 +755,7 @@ namespace netgen
int locpnum[2] = {0, 0};
for (j = 0; j < 2; j++)
for (int j = 0; j < 2; j++)
{
PointGeomInfo gi = (j == 0) ? gi0 : gi1;
@ -788,13 +798,17 @@ namespace netgen
// int noldpoints = mesh->GetNP();
int noldsurfel = mesh.GetNSE();
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;
try {
static Timer t("GenerateMesh"); RegionTimer reg(t);
res = meshing.GenerateMesh (mesh, mparam, maxh, k);
}
@ -811,6 +825,7 @@ namespace netgen
}
projecttype = PARAMETERSPACE;
static Timer t1("rest of loop"); RegionTimer reg1(t1);
if (res != MESHING2_OK)
{
@ -843,14 +858,14 @@ namespace netgen
notrys = 1;
for (i = oldnf+1; i <= mesh.GetNSE(); i++)
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)
{
@ -860,8 +875,8 @@ namespace netgen
if (geom.facemeshstatus[i-1] == -1)
{
cout << "Face " << i << endl;
// problemfile << "problem with face " << i << endl;
// problemfile << "vertices: " << 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() )
{
@ -873,22 +888,22 @@ namespace netgen
{
TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current());
gp_Pnt point = BRep_Tool::Pnt(vertex);
// problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl;
// problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl;
}
}
}
// problemfile << endl;
// problemfile << endl;
}
cout << endl << endl;
cout << "for more information open IGES/STEP Topology Explorer" << endl;
// problemfile.close();
// problemfile.close();
throw NgException ("Problem in Surface mesh generation");
}
else
{
// problemfile << "OK" << endl << endl;
// problemfile.close();
// problemfile << "OK" << endl << endl;
// problemfile.close();
}
@ -902,7 +917,7 @@ namespace netgen
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;
@ -913,7 +928,7 @@ namespace netgen
FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
PrintMessage (1, "Optimize Surface ", k);
for (i = 1; i <= mparam.optsteps2d; i++)
for (int i = 1; i <= mparam.optsteps2d; i++)
{
// (*testout) << "optstep " << i << endl;
if (multithread.terminate) return;
@ -988,7 +1003,8 @@ namespace netgen
void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh)
void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh,
const MeshingParameters & mparam)
{
mesh.SetGlobalH (mparam.maxh);
mesh.SetMinimalH (mparam.minh);
@ -1117,7 +1133,7 @@ namespace netgen
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;
}
@ -1165,7 +1181,7 @@ namespace netgen
//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);
RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0, 0, mparam);
//cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
}
}
@ -1298,7 +1314,7 @@ namespace netgen
mesh = make_shared<Mesh>();
mesh->geomtype = Mesh::GEOM_OCC;
OCCSetLocalMeshSize(geom,*mesh);
OCCSetLocalMeshSize(geom,*mesh, mparam);
}
if (multithread.terminate || mparam.perfstepsend <= MESHCONST_ANALYSE)
@ -1306,7 +1322,7 @@ namespace netgen
if (mparam.perfstepsstart <= MESHCONST_MESHEDGES)
{
OCCFindEdges (geom, *mesh);
OCCFindEdges (geom, *mesh, mparam);
/*
cout << "Removing redundant points" << endl;
@ -1379,7 +1395,7 @@ namespace netgen
if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
{
OCCMeshSurface (geom, *mesh, mparam.perfstepsend);
OCCMeshSurface (geom, *mesh, mparam.perfstepsend, mparam);
if (multithread.terminate) return TCL_OK;
#ifdef LOG_STREAM
@ -1407,7 +1423,7 @@ namespace netgen
MESHING3_RESULT res = MeshVolume (mparam, *mesh);
/*
/*
ofstream problemfile("occmesh.rep",ios_base::app);
problemfile << "VOLUMEMESHING" << endl << endl;
@ -1418,7 +1434,7 @@ namespace netgen
<< mesh->GetNE() << " elements" << endl << endl;
problemfile.close();
*/
*/
if (res != MESHING3_OK) return TCL_ERROR;

2
libsrc/occ/occgeom.cpp
View File

@ -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();

38
libsrc/occ/occgeom.hpp
View File

@ -114,7 +114,7 @@ namespace netgen
{
#include "occmeshsurf.hpp"
extern DLL_HEADER MeshingParameters mparam;
// extern DLL_HEADER MeshingParameters mparam;
#define PROJECTION_TOLERANCE 1e-10
@ -171,10 +171,8 @@ namespace netgen
{
public:
Point<3> p0, p1;
double Dist (Line l);
double Length ();
double Length () { return (p1-p0).Length(); }
};
@ -250,22 +248,22 @@ namespace netgen
DLL_HEADER void BuildFMap();
Box<3> GetBoundingBox()
{ return boundingbox;}
Box<3> GetBoundingBox() const
{ return boundingbox; }
int NrSolids()
{ return somap.Extent();}
int NrSolids() const
{ return somap.Extent(); }
// Philippose - 17/01/2009
// Total number of faces in the geometry
int NrFaces()
{ return fmap.Extent();}
int NrFaces() const
{ return fmap.Extent(); }
void SetCenter()
{ center = boundingbox.Center();}
{ center = boundingbox.Center(); }
Point<3> Center()
{ return center;}
Point<3> Center() const
{ return center; }
void Project (int surfi, Point<3> & p) const;
bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
@ -300,7 +298,7 @@ namespace netgen
// 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)
void SetFaceMaxH(int facenr, double faceh, const MeshingParameters & mparam)
{
if((facenr> 0) && (facenr <= fmap.Extent()))
{
@ -347,9 +345,7 @@ namespace netgen
// Returns the index of the currently selected face
int SelectedFace()
{
int i;
for(i = 1; i <= fmap.Extent(); i++)
for(int i = 1; i <= fmap.Extent(); i++)
{
if(face_sel_status[i-1])
{
@ -386,7 +382,7 @@ namespace netgen
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> & mesh, MeshingParameters & mparam);
@ -444,11 +440,11 @@ namespace netgen
DLL_HEADER extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr<Mesh> & 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

58
libsrc/occ/occmeshsurf.cpp
View File

@ -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)
{
;

4
libsrc/occ/occmeshsurf.hpp
View File

@ -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 ()

7
libsrc/occ/occpkg.cpp
View File

@ -27,6 +27,7 @@ namespace netgen
{
extern DLL_HEADER shared_ptr<NetgenGeometry> ng_geometry;
extern DLL_HEADER shared_ptr<Mesh> 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);
}
}

61
libsrc/occ/python_occ.cpp
View File

@ -3,6 +3,7 @@
#include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include "../meshing/python_mesh.hpp"
#include <meshing.hpp>
#include <occgeom.hpp>
@ -19,6 +20,21 @@ DLL_HEADER void ExportNgOCC(py::module &m)
{
py::class_<OCCGeometry, shared_ptr<OCCGeometry>, 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<OCCGeometry> 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<OCCGeometry>())
.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<py::gil_scoped_release>())
;
m.def("LoadOCCGeometry",FunctionPointer([] (const string & filename)
.def("GenerateMesh", [](shared_ptr<OCCGeometry> geo,
MeshingParameters* pars, py::kwargs kwargs)
{
cout << "load OCC geometry";
MeshingParameters mp;
if(pars) mp = *pars;
{
py::gil_scoped_acquire aq;
CreateMPfromKwargs(mp, kwargs);
}
auto mesh = make_shared<Mesh>();
SetGlobalMesh(mesh);
mesh->SetGeometry(geo);
ng_geometry = geo;
geo->GenerateMesh(mesh,mp);
return mesh;
}, py::arg("mp") = nullptr,
py::call_guard<py::gil_scoped_release>(),
meshingparameter_description.c_str())
;
m.def("LoadOCCGeometry",[] (const string & filename)
{
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<OCCGeometry>(instance, NOOP_Deleter);
return ng_geometry;
}),py::call_guard<py::gil_scoped_release>());
m.def("GenerateMesh", FunctionPointer([] (shared_ptr<OCCGeometry> geo, MeshingParameters &param)
{
auto mesh = make_shared<Mesh>();
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::gil_scoped_release>())
;
},py::call_guard<py::gil_scoped_release>());
}
PYBIND11_MODULE(libNgOCC, m) {

39
libsrc/stlgeom/python_stl.cpp
View File

@ -4,6 +4,7 @@
#include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include <stlgeom.hpp>
#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_<STLGeometry,shared_ptr<STLGeometry>, NetgenGeometry> (m,"STLGeometry")
.def(py::init<>())
.def(py::init<>([](const string& filename)
{
ifstream ist(filename);
return shared_ptr<STLGeometry>(STLGeometry::Load(ist));
}), py::arg("filename"),
py::call_guard<py::gil_scoped_release>())
.def(NGSPickle<STLGeometry>())
.def("_visualizationData", [](shared_ptr<STLGeometry> stl_geo)
{
@ -71,29 +78,31 @@ DLL_HEADER void ExportSTL(py::module & m)
res["max"] = MoveToNumpy(max);
return res;
}, py::call_guard<py::gil_scoped_release>())
;
m.def("LoadSTLGeometry", FunctionPointer([] (const string & filename)
.def("GenerateMesh", [] (shared_ptr<STLGeometry> geo,
MeshingParameters* pars, py::kwargs kwargs)
{
ifstream ist(filename);
return shared_ptr<STLGeometry>(STLGeometry::Load(ist));
}),py::call_guard<py::gil_scoped_release>());
m.def("GenerateMesh", FunctionPointer([] (shared_ptr<STLGeometry> geo, MeshingParameters &param)
MeshingParameters mp;
if(pars) mp = *pars;
{
py::gil_scoped_acquire aq;
CreateMPfromKwargs(mp, kwargs);
}
auto mesh = make_shared<Mesh>();
SetGlobalMesh(mesh);
mesh->SetGeometry(geo);
ng_geometry = geo;
try
{
geo->GenerateMesh(mesh,param);
}
catch (NgException ex)
{
cout << "Caught NgException: " << ex.What() << endl;
}
geo->GenerateMesh(mesh,mp);
return mesh;
}),py::call_guard<py::gil_scoped_release>())
}, py::arg("mp") = nullptr,
py::call_guard<py::gil_scoped_release>(),
meshingparameter_description.c_str())
;
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>(STLGeometry::Load(ist));
},py::call_guard<py::gil_scoped_release>());
}
PYBIND11_MODULE(libstl, m) {

6
nglib/nglib.cpp
View File

@ -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();

2
python/CMakeLists.txt
View File

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

13
python/NgOCC.py
View File

@ -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 *

26
python/csg.py
View File

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

30
python/geom2d.py
View File

@ -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']

4
python/gui.py
View File

@ -19,4 +19,8 @@ def StartGUI():
pass
if not netgen.libngpy._meshing._netgen_executable_started:
# catch exception for building html docu on server without display
try:
StartGUI()
except:
pass

23
python/meshing.py
View File

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

2
python/occ.py Normal file
View File

@ -0,0 +1,2 @@
from .libngpy._NgOCC import *
from .meshing import meshsize

12
python/stl.py
View File

@ -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

4
tests/pytest/test_pickling.py
View File

@ -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()

2
tests/pytest/test_savemesh.py
View File

@ -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()