ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-15 10:28:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/master' into boundarylayer_fixes

Browse Source
This commit is contained in:
Matthias Hochsteger 2024-01-08 10:38:12 +01:00
commit d09c00796b
16 changed files with 264 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
libsrc/core/exception.cpp
View File

@ -80,7 +80,7 @@ namespace ngcore
// 1 libngcore.dylib 0x000000010ddb298c _ZL21ngcore_signal_handleri + 316
constexpr char reset_shell[] = "\033[0m";
constexpr char green[] = "\033[32m";
constexpr char yellow[] = "\033[33m";
[[maybe_unused]] constexpr char yellow[] = "\033[33m";
std::istringstream in(s);

1
libsrc/core/exception.hpp
View File

@ -10,7 +10,6 @@
namespace ngcore
{
NGCORE_API std::string GetBackTrace();
// Exception for code that shouldn't be executed

24
libsrc/interface/rw_medit.cpp
View File

@ -50,7 +50,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int nvert;
fin >> nvert;
Point<3> p{0.,0.,0.};
for(auto k : Range(nvert)) {
for([[maybe_unused]] auto k : Range(nvert)) {
for(auto i : Range(dim))
fin >> p[i];
fin >> index;
@ -61,7 +61,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int nedge;
fin >> nedge;
Segment seg;
for(auto k : Range(nedge)) {
for([[maybe_unused]] auto k : Range(nedge)) {
for(auto i : Range(2))
fin >> seg[i];
fin >> seg.edgenr;
@ -74,7 +74,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int ntrig, index;
fin >> ntrig;
Element2d sel;
for(auto k : Range(ntrig)) {
for([[maybe_unused]] auto k : Range(ntrig)) {
for(auto i : Range(3))
fin >> sel[i];
fin >> index;
@ -86,7 +86,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int ntet;
fin >> ntet;
Element el(4);
for(auto k : Range(ntet)) {
for([[maybe_unused]] auto k : Range(ntet)) {
for(auto i : Range(4))
fin >> el[i];
fin >> index;
@ -99,7 +99,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int ncorners;
fin >> ncorners;
Element0d el;
for(auto k : Range(ncorners)) {
for([[maybe_unused]] auto k : Range(ncorners)) {
fin >> el.pnum;
}
}
@ -107,7 +107,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int nverts;
fin >> nverts;
int vert;
for(auto k : Range(nverts)) {
for([[maybe_unused]] auto k : Range(nverts)) {
fin >> vert;
}
}
@ -115,7 +115,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int nnormals;
fin >> nnormals;
Vec<3> normal;
for(auto k : Range(nnormals)) {
for([[maybe_unused]] auto k : Range(nnormals)) {
fin >> normal[0];
fin >> normal[1];
fin >> normal[2];
@ -126,7 +126,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
fin >> nnormals;
int vert;
int normal;
for(auto k : Range(nnormals)) {
for([[maybe_unused]] auto k : Range(nnormals)) {
fin >> normal;
fin >> vert;
}
@ -135,7 +135,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int ntangents;
fin >> ntangents;
Vec<3> tangent;
for(auto k : Range(ntangents)) {
for([[maybe_unused]] auto k : Range(ntangents)) {
fin >> tangent[0];
fin >> tangent[1];
fin >> tangent[2];
@ -146,7 +146,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
fin >> ntangents;
int vert;
int tangent;
for(auto k : Range(ntangents)) {
for([[maybe_unused]] auto k : Range(ntangents)) {
fin >> tangent;
fin >> vert;
}
@ -155,7 +155,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int nridges;
fin >> nridges;
int ridge;
for(auto k : Range(nridges)) {
for([[maybe_unused]] auto k : Range(nridges)) {
fin >> ridge;
}
}
@ -164,7 +164,7 @@ void ReadMeditFormat (Mesh & mesh, const filesystem::path & filename, map<tuple<
int nitems;
fin >> nitems;
string s;
for(auto i : Range(nitems))
for([[maybe_unused]] auto i : Range(nitems))
fin >> s; // read one line
}
}

2
libsrc/meshing/basegeom.cpp
View File

@ -829,7 +829,7 @@ namespace netgen
if(face.primary == &face)
{
// check if this face connects two identified closesurfaces
auto & idents = mesh.GetIdentifications();
// auto & idents = mesh.GetIdentifications();
std::set<int> relevant_edges;
auto segments = face.GetBoundary(mesh);
for(const auto &s : segments)

201
libsrc/meshing/boundarylayer.cpp
View File

@ -572,6 +572,7 @@ struct GrowthVectorLimiter {
limiter.LimitGrowthVector(pi_to, sei, trig_shift, seg_shift);
});
<<<<<<< HEAD
// for (auto i : Range(limits))
// if(limits[i] < 1.0)
// cout << i << ": " << limits[i] << endl;
@ -580,6 +581,200 @@ struct GrowthVectorLimiter {
auto pi_from = limiter.map_from[pi_to];
if(pi_from.IsValid())
limits[pi_from] = min(limits[pi_from], limits[pi_to]);
=======
// Calculate parallel projections
Vec<3> ab_base_norm = (b_base - a_base).Normalize();
double a_vec_x = Dot(a_vec, ab_base_norm);
double b_vec_x = Dot(b_vec, -ab_base_norm);
// double ratio_parallel = (a_vec_x + b_vec_x) / ab_base;
// Calculate surface normal at point si
Vec<3> surface_normal = getNormal(mesh[si]);
double a_vec_y = abs(Dot(a_vec, surface_normal));
double b_vec_y = abs(Dot(b_vec, surface_normal));
double diff_perpendicular = abs(a_vec_y - b_vec_y);
double tan_alpha = diff_perpendicular / (ab_base - a_vec_x - b_vec_x);
double TAN_ALPHA_LIMIT = 0.36397; // Approximately 20 degrees in radians
if (tan_alpha > TAN_ALPHA_LIMIT) {
if (a_vec_y > b_vec_y) {
double correction = (TAN_ALPHA_LIMIT / tan_alpha * diff_perpendicular + b_vec_y) / a_vec_y;
limits[pi1] *= correction;
}
else {
double correction = (TAN_ALPHA_LIMIT / tan_alpha * diff_perpendicular + a_vec_y) / b_vec_y;
limits[pi2] *= correction;
}
}
};
auto neighbour_limiter = [&](PointIndex pi1, PointIndex pi2, SurfaceElementIndex si) {
parallel_limiter(pi1, pi2, si);
perpendicular_limiter(pi1, pi2, si);
};
auto modifiedsmooth = [&](size_t nsteps) {
for ([[maybe_unused]] auto i : Range(nsteps))
for (SurfaceElementIndex sei : mesh.SurfaceElements().Range())
{
// assuming triangle
neighbour_limiter(mesh[sei].PNum(1), mesh[sei].PNum(2), sei);
neighbour_limiter(mesh[sei].PNum(2), mesh[sei].PNum(3), sei);
neighbour_limiter(mesh[sei].PNum(3), mesh[sei].PNum(1), sei);
}
};
/*
auto smooth = [&] (size_t nsteps) {
for([[maybe_unused]] auto i : Range(nsteps))
for(const auto & sel : mesh.SurfaceElements())
{
double min_limit = 999;
for(auto pi : sel.PNums())
min_limit = min(min_limit, limits[pi]);
for(auto pi : sel.PNums())
limits[pi] = min(limits[pi], 1.4*min_limit);
}
};
*/
// check for self-intersection within new elements (prisms/hexes)
auto self_intersection = [&] () {
for(SurfaceElementIndex sei : mesh.SurfaceElements().Range())
{
auto facei = mesh[sei].GetIndex();
if(facei < nfd_old && !params.surfid.Contains(facei))
continue;
auto sel = mesh[sei];
auto np = sel.GetNP();
// check if a new edge intesects the plane of any opposing face
double lam;
for(auto i : Range(np))
for(auto fi : Range(np-2))
if(isIntersectingPlane(GetMappedSeg(sel[i]), GetMappedFace(sei, i+fi+1), lam))
if(lam < 1.0)
limits[sel[i]] *= lam;
}
};
// first step: intersect with other surface elements that are boundary of domain the layer is grown into
// second (and subsequent) steps: intersect with other boundary layers, allow restriction by 20% in each step
auto changed_domains = domains;
if(!params.outside)
changed_domains.Invert();
bool limit_reached = true;
double lam_lower_limit = 1.0;
int step = 0;
while(limit_reached || step<3)
{
Array<double, PointIndex> new_limits;
new_limits.SetSize(np);
new_limits = 1.0;
if(step>1)
lam_lower_limit *= 0.8;
limit_reached = false;
// build search tree with all surface elements (bounding box of a surface element also covers the generated boundary layer)
Box<3> bbox(Box<3>::EMPTY_BOX);
for(auto pi : mesh.Points().Range())
{
bbox.Add(mesh[pi]);
bbox.Add(mesh[pi]+limits[pi]*height*growthvectors[pi]);
}
BoxTree<3> tree(bbox);
for(auto sei : mesh.SurfaceElements().Range())
{
const auto & sel = mesh[sei];
Box<3> box(Box<3>::EMPTY_BOX);
const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
if(!changed_domains.Test(fd.DomainIn()) &&
!changed_domains.Test(fd.DomainOut()))
continue;
for(auto pi : sel.PNums())
box.Add(mesh[pi]);
// also add moved points to bounding box
if(params.surfid.Contains(sel.GetIndex()))
for(auto pi : sel.PNums())
box.Add(mesh[pi]+limits[pi]*height*growthvectors[pi]);
tree.Insert(box, sei);
}
for(auto pi : mesh.Points().Range())
{
if(mesh[pi].Type() == INNERPOINT)
continue;
if(growthvectors[pi].Length2() == 0.0)
continue;
Box<3> box(Box<3>::EMPTY_BOX);
auto seg = GetMappedSeg(pi);
box.Add(seg[0]);
box.Add(seg[1]);
double lam = 1.0;
tree.GetFirstIntersecting(box.PMin(), box.PMax(), [&](SurfaceElementIndex sei)
{
const auto & sel = mesh[sei];
if(sel.PNums().Contains(pi))
return false;
auto face = GetMappedFace(sei, -2);
double lam_ = 999;
bool is_bl_sel = params.surfid.Contains(sel.GetIndex());
if (step == 0)
{
face = GetMappedFace(sei, -1);
if (isIntersectingFace(seg, face, lam_))
{
if (is_bl_sel)
lam_ *= params.limit_safety;
lam = min(lam, lam_);
}
}
if(step==1)
{
if(isIntersectingFace(seg, face, lam_))
{
if(is_bl_sel) // allow only half the distance if the opposing surface element has a boundary layer too
lam_ *= params.limit_safety;
lam = min(lam, lam_);
}
}
// if the opposing surface element has a boundary layer, we need to additionally intersect with the new faces
if(step>1 && is_bl_sel)
{
for(auto facei : Range(-1, sel.GetNP()))
{
auto face = GetMappedFace(sei, facei);
if(isIntersectingFace(seg, face, lam_)) // && lam_ > other_limit)
{
lam = min(lam, lam_);
}
}
}
return false;
});
if(lam<1)
{
if(lam<lam_lower_limit && step>1)
{
limit_reached = true;
lam = lam_lower_limit;
}
}
new_limits[pi] = min(limits[pi], lam* limits[pi]);
}
step++;
limits = new_limits;
if (step > 0)
modifiedsmooth(1);
>>>>>>> origin/master
}
for(auto i : Range(growthvectors))
@ -1012,8 +1207,12 @@ struct GrowthVectorLimiter {
for(auto pi : points)
{
auto sels = p2sel[pi];
<<<<<<< HEAD
Vec<3> new_gw = getGW(pi);
new_gw = 0.;
=======
Vec<3> new_gw = growthvectors[pi];
>>>>>>> origin/master
// int cnt = 1;
std::set<PointIndex> suround;
suround.insert(pi);
@ -2113,7 +2312,7 @@ struct GrowthVectorLimiter {
points.Append(pi);
auto p2el = mesh.CreatePoint2ElementTable(is_inner_point);
// smooth growth vectors to shift additional element layers to the inside and fix flipped tets
for([[maybe_unused]] auto step : Range(10))
{

18
libsrc/meshing/boundarylayer2d.cpp
View File

@ -165,11 +165,12 @@ namespace netgen
}
auto oldnf = mesh.GetNSE();
auto res = meshing.GenerateMesh (mesh, mp, mp.maxh, domain);
// auto res =
meshing.GenerateMesh (mesh, mp, mp.maxh, domain);
for (SurfaceElementIndex sei : Range(oldnf, mesh.GetNSE()))
mesh[sei].SetIndex (domain);
int hsteps = mp.optsteps2d;
// int hsteps = mp.optsteps2d;
const char * optstr = mp.optimize2d.c_str();
MeshOptimize2d meshopt(mesh);
@ -219,7 +220,7 @@ namespace netgen
int np = mesh.GetNP();
int nseg = line_segments.Size();
int ne = mesh.GetNSE();
// int ne = mesh.GetNSE();
mesh.UpdateTopology();
double total_thickness = 0.0;
@ -246,7 +247,7 @@ namespace netgen
Array<int> si_map(mesh.GetNFD()+2);
si_map = -1;
int fd_old = mesh.GetNFD();
// int fd_old = mesh.GetNFD();
int max_edge_nr = -1;
int max_domain = -1;
@ -284,7 +285,8 @@ namespace netgen
{
FaceDescriptor new_fd(0, 0, 0, -1);
new_fd.SetBCProperty(new_domain);
int new_fd_index = mesh.AddFaceDescriptor(new_fd);
// int new_fd_index =
mesh.AddFaceDescriptor(new_fd);
if(should_make_new_domain)
mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(domain-1));
}
@ -625,7 +627,7 @@ namespace netgen
auto & pm0 = mapto[seg[0]];
auto & pm1 = mapto[seg[1]];
auto newindex = si_map[domain];
// auto newindex = si_map[domain];
Segment s = seg;
s.geominfo[0] = {};
@ -663,7 +665,7 @@ namespace netgen
auto p0 = mesh[pi0];
auto p1 = mesh[pi1];
auto q0 = mesh[pi2];
auto q1 = mesh[pi3];
// auto q1 = mesh[pi3];
Vec<2> n = {-p1[1]+p0[1], p1[0]-p0[0]};
Vec<2> v = { q0[0]-p0[0], q0[1]-p0[1]};

2
libsrc/meshing/meshclass.cpp
View File

@ -7416,6 +7416,7 @@ namespace netgen
auto v = pmax-pmin;
double eps = v.Length()*1e-8;
/*
auto onPlane = [&] (const MeshPoint & p) -> bool
{
auto v = p_plane-p;
@ -7425,6 +7426,7 @@ namespace netgen
// auto ret = fabs(v*n_plane)/l;
return fabs(v*n_plane) < eps;
};
*/
/*
auto mirror = [&] (PointIndex pi) -> PointIndex

13
libsrc/meshing/parallelmesh.cpp
View File

@ -1856,13 +1856,15 @@ namespace netgen
void Mesh :: PartHybridMesh ()
{
#ifdef METIS
throw Exception("PartHybridMesh not supported");
#ifdef METISxxx
int ne = GetNE();
int nn = GetNP();
int nedges = topology.GetNEdges();
idxtype *xadj, * adjacency, *v_weights = NULL, *e_weights = NULL;
idxtype *xadj, * adjacency;
// idxtype *v_weights = NULL, *e_weights = NULL;
int weightflag = 0;
int numflag = 0;
@ -1953,6 +1955,8 @@ namespace netgen
void Mesh :: PartDualHybridMesh ( ) // NgArray<int> & neloc )
{
throw Exception("PartDualHybridMesh not supported");
#ifdef OLD
#ifdef METIS
int ne = GetNE();
@ -2065,7 +2069,8 @@ namespace netgen
#else
cout << "partdualmesh not available" << endl;
#endif
#endif
}
@ -2131,7 +2136,6 @@ namespace netgen
idxtype *v_weights = NULL, *e_weights = NULL;
idxtype weightflag = 0;
// int numflag = 0;
idxtype nparts = ntasks - 1;
@ -2142,6 +2146,7 @@ namespace netgen
BubbleSort (adjacency.Range (xadj[el], xadj[el+1]));
#ifdef METIS4
idxtype weightflag = 0;
int options[5];
options[0] = 0;
METIS_PartGraphKway ( &ne, &xadj[0], &adjacency[0], v_weights, e_weights, &weightflag,

28
libsrc/meshing/topology.cpp
View File

@ -201,24 +201,24 @@ namespace netgen
INDEX_4 face4(el[elfaces[j][0]], el[elfaces[j][1]],
el[elfaces[j][2]], el[elfaces[j][3]]);
int facedir = 0;
// int facedir = 0;
if (min2 (face4.I1(), face4.I2()) >
min2 (face4.I4(), face4.I3()))
{ // z - flip
facedir += 1;
// facedir += 1;
swap (face4.I1(), face4.I4());
swap (face4.I2(), face4.I3());
}
if (min2 (face4.I1(), face4.I4()) >
min2 (face4.I2(), face4.I3()))
{ // x - flip
facedir += 2;
// facedir += 2;
swap (face4.I1(), face4.I2());
swap (face4.I3(), face4.I4());
}
if (face4.I2() > face4.I4())
{ // diagonal flip
facedir += 4;
// facedir += 4;
swap (face4.I2(), face4.I4());
}
@ -263,27 +263,27 @@ namespace netgen
{ // triangle
// int facenum;
int facedir;
// int facedir;
INDEX_4 face(el.PNum(elfaces[0][0]),
el.PNum(elfaces[0][1]),
el.PNum(elfaces[0][2]),0);
facedir = 0;
// facedir = 0;
if (face.I1() > face.I2())
{
swap (face.I1(), face.I2());
facedir += 1;
// facedir += 1;
}
if (face.I2() > face.I3())
{
swap (face.I2(), face.I3());
facedir += 2;
// facedir += 2;
}
if (face.I1() > face.I2())
{
swap (face.I1(), face.I2());
facedir += 4;
// facedir += 4;
}
if (face.I1() != v) continue;
@ -312,31 +312,31 @@ namespace netgen
{
// quad
// int facenum;
int facedir;
// int facedir;
INDEX_4 face4(el.PNum(elfaces[0][0]),
el.PNum(elfaces[0][1]),
el.PNum(elfaces[0][2]),
el.PNum(elfaces[0][3]));
facedir = 0;
// facedir = 0;
if (min2 (face4.I1(), face4.I2()) >
min2 (face4.I4(), face4.I3()))
{ // z - orientation
facedir += 1;
// facedir += 1;
swap (face4.I1(), face4.I4());
swap (face4.I2(), face4.I3());
}
if (min2 (face4.I1(), face4.I4()) >
min2 (face4.I2(), face4.I3()))
{ // x - orientation
facedir += 2;
// facedir += 2;
swap (face4.I1(), face4.I2());
swap (face4.I3(), face4.I4());
}
if (face4.I2() > face4.I4())
{
facedir += 4;
// facedir += 4;
swap (face4.I2(), face4.I4());
}

2
libsrc/occ/Partition_Inter2d.cxx
View File

@ -516,7 +516,7 @@ static void EdgesPartition(const TopoDS_Face& F,
}
}
Standard_Boolean AffichPurge = Standard_False;
// Standard_Boolean AffichPurge = Standard_False;
if ( LV1.IsEmpty()) return;

6
libsrc/occ/Partition_Loop.cxx
View File

@ -67,7 +67,7 @@
#include <TopExp_Explorer.hxx>
static char* name = new char[100];
static int nbe = 0;
// static int nbe = 0;
#ifdef WIN32
#define M_PI 3.14159265358979323846
@ -198,7 +198,7 @@ static Standard_Boolean SelectEdge(const TopoDS_Face& F,
Cc->D2(uc, PC, CTg1, CTg2);
C->D2(u, P, Tg1, Tg2);
Standard_Real angle;
Standard_Real angle = 0.0;
if (CE.Orientation () == TopAbs_REVERSED && E.Orientation () == TopAbs_FORWARD) {
angle = CTg1.Angle(Tg1.Reversed());
@ -363,7 +363,7 @@ void Partition_Loop::Perform()
}
}
int i = 0;
// int i = 0;
while (!End) {
//-------------------------------
// Construction of a wire.

4
libsrc/occ/Partition_Spliter.cxx
View File

@ -1459,7 +1459,7 @@ void Partition_Spliter::MakeEdges (const TopoDS_Edge& E,
if (VOnE.Extent() < 3) { // do not rebuild not cut edge
if (( VF.IsSame( VOnE.First() ) && VL.IsSame( VOnE.Last() )) ||
VL.IsSame( VOnE.First() ) && VF.IsSame( VOnE.Last() ) ) {
(VL.IsSame( VOnE.First() ) && VF.IsSame( VOnE.Last() )) ) {
NE.Append( E );
return;
}
@ -1476,7 +1476,7 @@ void Partition_Spliter::MakeEdges (const TopoDS_Edge& E,
if (SV.Length() < 3) { // do not rebuild not cut edge
if (( VF.IsSame( SV.First() ) && VL.IsSame( SV.Last() )) ||
VL.IsSame( SV.First() ) && VF.IsSame( SV.Last() ) ) {
( VL.IsSame( SV.First() ) && VF.IsSame( SV.Last() )) ) {
NE.Append( E );
return;
}

2
libsrc/occ/occgeom.cpp
View File

@ -1714,11 +1714,13 @@ namespace netgen
BRepTools::Read(shape, ss, builder);
}
/*
// enumerate shapes and archive only integers
auto my_hash = [](const TopoDS_Shape & key) {
auto occ_hash = key.HashCode(1<<31UL);
return std::hash<decltype(occ_hash)>()(occ_hash);
};
*/
TopTools_IndexedMapOfShape shape_map;
Array<TopoDS_Shape> shape_list;

4
libsrc/occ/python_occ_shapes.cpp
View File

@ -1935,7 +1935,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
optional<TopoDS_Wire> auxspine) {
if (twist)
{
auto [pnt, angle] = *twist;
// auto [pnt, angle] = *twist;
/*
cyl = Cylinder((0,0,0), Z, r=1, h=1).faces[0]
@ -2230,7 +2230,7 @@ tangents : Dict[int, gp_Vec2d]
TColgp_Array1OfPnt poles(0, vpoles.size()-1);
TColStd_Array1OfReal knots(0, vpoles.size()+degree);
TColStd_Array1OfInteger mult(0, vpoles.size()+degree);
int cnt = 0;
// int cnt = 0;
for (int i = 0; i < vpoles.size(); i++)
{

2
libsrc/stlgeom/python_stl.cpp
View File

@ -278,7 +278,7 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
throw Exception("Point not found in triangle");
self->SetSelectTrig(trignr);
self->SetNodeOfSelTrig(point_in_trig);
for(auto i : Range(count))
for([[maybe_unused]] auto i : Range(count))
self->MoveSelectedPointToMiddle();
})
.def("Save", &STLGeometry::Save)

6
nglib/nglib_occ.cpp
View File

@ -148,13 +148,15 @@ namespace nglib
numpoints = me->GetNP();
// Initially set up only for surface meshing without any optimisation
int perfstepsend = MESHCONST_MESHSURFACE;
// int perfstepsend = MESHCONST_MESHSURFACE;
// Check and if required, enable surface mesh optimisation step
/*
if(mp->optsurfmeshenable)
{
perfstepsend = MESHCONST_OPTSURFACE;
perfstepsend = MESHCONST_OPTSURFACE;
}
*/
occgeom->MeshSurface(*me, mparam);
occgeom->OptimizeSurface(*me, mparam);