diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index 03df4df3..5e891639 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -413,7 +413,7 @@ namespace netgen geom.GetSurface((mesh.GetFaceDescriptor(k).SurfNr())); - Meshing2Surfaces meshing(*surf, geom.BoundingBox()); + Meshing2Surfaces meshing(*surf, mparam, geom.BoundingBox()); meshing.SetStartTime (starttime); double eps = 1e-8 * geom.MaxSize(); @@ -482,7 +482,7 @@ namespace netgen mparam.checkoverlap = 0; MESHING2_RESULT res = - meshing.GenerateMesh (mesh, maxh, k); + meshing.GenerateMesh (mesh, mparam, maxh, k); if (res != MESHING2_OK) { @@ -531,7 +531,7 @@ namespace netgen meshopt.SetMetricWeight (mparam.elsizeweight); meshopt.SetWriteStatus (0); - meshopt.ImproveMesh (mesh); + meshopt.ImproveMesh (mesh, mparam); } { @@ -553,7 +553,7 @@ namespace netgen meshopt.SetMetricWeight (mparam.elsizeweight); meshopt.SetWriteStatus (0); - meshopt.ImproveMesh (mesh); + meshopt.ImproveMesh (mesh, mparam); } } } @@ -723,7 +723,7 @@ namespace netgen if (mparam.uselocalh) { - mesh->CalcLocalH(); + mesh->CalcLocalH(mparam.grading); mesh->DeleteMesh(); FindPoints (geom, *mesh); @@ -757,7 +757,7 @@ namespace netgen if (mparam.uselocalh && 0) { - mesh->CalcLocalH(); + mesh->CalcLocalH(mparam.grading); mesh->DeleteMesh(); FindPoints (geom, *mesh); diff --git a/libsrc/csg/meshsurf.cpp b/libsrc/csg/meshsurf.cpp index 0921ea4e..19e12eb2 100644 --- a/libsrc/csg/meshsurf.cpp +++ b/libsrc/csg/meshsurf.cpp @@ -15,8 +15,9 @@ Meshing2Surfaces :: Meshing2Surfaces (const Surface & asurface) } */ Meshing2Surfaces :: Meshing2Surfaces (const Surface & asurf, + const MeshingParameters & mp, const Box<3> & abb) - : Meshing2(abb), surface(asurf) + : Meshing2(mp, abb), surface(asurf) { ; } diff --git a/libsrc/csg/meshsurf.hpp b/libsrc/csg/meshsurf.hpp index 72a40eaf..e6c69225 100644 --- a/libsrc/csg/meshsurf.hpp +++ b/libsrc/csg/meshsurf.hpp @@ -14,7 +14,8 @@ namespace netgen /// // Meshing2Surfaces (const Surface & asurf); /// - Meshing2Surfaces (const Surface & asurf, const Box<3> & aboundingbox); + Meshing2Surfaces (const Surface & asurf, const MeshingParameters & mp, + const Box<3> & aboundingbox); protected: /// diff --git a/libsrc/general/flags.cpp b/libsrc/general/flags.cpp index ba6b8356..e5916f87 100644 --- a/libsrc/general/flags.cpp +++ b/libsrc/general/flags.cpp @@ -125,8 +125,8 @@ namespace netgen return *strlistflags.Get(name); else { - static Array hstra(0); - return hstra; + static Array dummy_array(0); + return dummy_array; } } @@ -137,8 +137,8 @@ namespace netgen return *numlistflags.Get(name); else { - static Array hnuma(0); - return hnuma; + static Array dummy_array(0); + return dummy_array; } } diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 024ad588..5b6d43e0 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -350,7 +350,7 @@ namespace netgen mesh->SetLocalH (pmin, pmax, mparam.grading); mesh->SetGlobalH (h); - mesh->CalcLocalH(); + mesh->CalcLocalH(mparam.grading); int bnp = mesh->GetNP(); // boundary points @@ -452,7 +452,7 @@ namespace netgen mparam.quad = hquad || geometry.GetDomainQuadMeshing (domnr); - Meshing2 meshing (Box<3> (pmin, pmax)); + Meshing2 meshing (mparam, Box<3> (pmin, pmax)); Array compress(bnp); compress = -1; @@ -491,7 +491,7 @@ namespace netgen meshing.Delaunay (*mesh, domnr, mparam); else */ - meshing.GenerateMesh (*mesh, h, domnr); + meshing.GenerateMesh (*mesh, mparam, h, domnr); for (SurfaceElementIndex sei = oldnf; sei < mesh->GetNSE(); sei++) (*mesh)[sei].SetIndex (domnr); diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index ac05082f..8c937760 100644 --- a/libsrc/interface/nginterface.cpp +++ b/libsrc/interface/nginterface.cpp @@ -2152,7 +2152,7 @@ int Ng_Bisect_WithInfo ( const char * refinementfile, double ** qualityloss, int } if(!mesh->LocalHFunctionGenerated()) - mesh->CalcLocalH(); + mesh->CalcLocalH(mparam.grading); mesh->LocalHFunction().SetGrading (mparam.grading); diff --git a/libsrc/linalg/densemat.cpp b/libsrc/linalg/densemat.cpp index 11596654..a0066e8f 100644 --- a/libsrc/linalg/densemat.cpp +++ b/libsrc/linalg/densemat.cpp @@ -160,39 +160,12 @@ namespace netgen return *this; } - - - - /* - double & DenseMatrix :: operator() (int i, int j) - { - if (i >= 1 && j >= 1 && i <= height && j <= width) - return Elem(i,j); - else (*myerr) << "DenseMatrix: index (" << i << "," << j << ") out of range (1.." - << height << ",1.." << width << ")\n"; - static double dummy = 0; - return dummy; - } - - double DenseMatrix :: operator() (int i, int j) const - { - if (i >= 1 && j >= 1 && i <= height && j <= width) - return Get(i,j); - else (*myerr) << "DenseMatrix: index (" << i << "," << j << ") out of range (1.." - << height << ",1.." << width << ")\n"; - - static double dummy = 0; - return dummy; - } - */ - DenseMatrix & DenseMatrix :: operator= (double v) { - int i; double * p = data; if (data) - for (i = width*height; i > 0; i--, p++) + for (int i = width*height; i > 0; i--, p++) *p = v; return *this; @@ -202,11 +175,10 @@ namespace netgen DenseMatrix & DenseMatrix :: operator*= (double v) { - int i; double * p = data; if (data) - for (i = width*height; i > 0; i--, p++) + for (int i = width*height; i > 0; i--, p++) *p *= v; return *this; diff --git a/libsrc/linalg/linsearch.cpp b/libsrc/linalg/linsearch.cpp index dcdb6644..a2dd38aa 100644 --- a/libsrc/linalg/linsearch.cpp +++ b/libsrc/linalg/linsearch.cpp @@ -47,10 +47,8 @@ double MinFunction :: FuncGrad (const Vector & x, Vector & g) const /* int n = x.Size(); - static Vector xr; - static Vector xl; - xr.SetSize(n); - xl.SetSize(n); + Vector xr(n); + Vector xl(n); double eps = 1e-6; double fl, fr; diff --git a/libsrc/meshing/improve2.hpp b/libsrc/meshing/improve2.hpp index aa049c5f..4a6e49ea 100644 --- a/libsrc/meshing/improve2.hpp +++ b/libsrc/meshing/improve2.hpp @@ -15,8 +15,8 @@ public: /// MeshOptimize2d (); /// - void ImproveMesh (Mesh & mesh2d); - void ImproveMeshJacobian (Mesh & mesh2d); + void ImproveMesh (Mesh & mesh2d, const MeshingParameters & mp); + void ImproveMeshJacobian (Mesh & mesh2d, const MeshingParameters & mp); void ImproveVolumeMesh (Mesh & mesh); void ProjectBoundaryPoints(Array & surfaceindex, const Array* > & from, Array* > & dest); diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index d7cf86d7..8fe8f916 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -161,7 +161,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh, { const Element & elem = mesh[hasonepi[k]]; double err = CalcTetBadness (mesh[elem[0]], mesh[elem[1]], - mesh[elem[2]], mesh[elem[3]], 0); + mesh[elem[2]], mesh[elem[3]], 0, mparam); bad2 += err; oneperr[k] = err; } diff --git a/libsrc/meshing/improve3.hpp b/libsrc/meshing/improve3.hpp index 100ef54a..6dd52881 100644 --- a/libsrc/meshing/improve3.hpp +++ b/libsrc/meshing/improve3.hpp @@ -22,11 +22,11 @@ public: inline double CalcBad (const Mesh::T_POINTS & points, const Element & elem, - double h) + double h, const MeshingParameters & mp = mparam) { if (elem.GetType() == TET) return CalcTetBadness (points[elem[0]], points[elem[1]], - points[elem[2]], points[elem[3]], h); + points[elem[2]], points[elem[3]], h, mp); return 0; } @@ -34,7 +34,8 @@ CalcBad (const Mesh::T_POINTS & points, const Element & elem, extern double CalcTotalBad (const Mesh::T_POINTS & points, - const Mesh::T_VOLELEMENTS & elements); + const Mesh::T_VOLELEMENTS & elements, + const MeshingParameters & mp = mparam); extern int WrongOrientation (const Mesh::T_POINTS & points, const Element & el); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 0e3bbbc6..28117b2e 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -1731,12 +1731,6 @@ namespace netgen void Mesh :: FindOpenElements (int dom) { static int timer = NgProfiler::CreateTimer ("Mesh::FindOpenElements"); - static int timera = NgProfiler::CreateTimer ("Mesh::FindOpenElements A"); - static int timerb = NgProfiler::CreateTimer ("Mesh::FindOpenElements B"); - static int timerc = NgProfiler::CreateTimer ("Mesh::FindOpenElements C"); - static int timerd = NgProfiler::CreateTimer ("Mesh::FindOpenElements D"); - static int timere = NgProfiler::CreateTimer ("Mesh::FindOpenElements E"); - NgProfiler::RegionTimer reg (timer); int np = GetNP(); @@ -1747,7 +1741,6 @@ namespace netgen numonpoint = 0; - NgProfiler::StartTimer (timera); for (ElementIndex ei = 0; ei < ne; ei++) { const Element & el = (*this)[ei]; @@ -1784,13 +1777,6 @@ namespace netgen elsonpoint.Add (el[j], ei); } } - NgProfiler::StopTimer (timera); - - - - - NgProfiler::StartTimer (timerb); - Array hasface(GetNFD()); @@ -1865,15 +1851,11 @@ namespace netgen } - NgProfiler::StopTimer (timerb); - int ii; PointIndex pi; SurfaceElementIndex sei; Element2d hel; - NgProfiler::RegionTimer regc (timerc); - INDEX_3_CLOSED_HASHTABLE faceht(100); openelements.SetSize(0); @@ -2016,8 +1998,6 @@ namespace netgen BuildBoundaryEdges(); - NgProfiler::RegionTimer regd (timerd); - for (int i = 1; i <= openelements.Size(); i++) { const Element2d & sel = openelements.Get(i); @@ -2041,7 +2021,6 @@ namespace netgen } - NgProfiler::RegionTimer rege (timere); /* for (i = 1; i <= GetNSeg(); i++) @@ -2542,13 +2521,14 @@ namespace netgen - void Mesh :: CalcLocalH () + void Mesh :: CalcLocalH (double grading) { if (!lochfunc) { Point3d pmin, pmax; GetBox (pmin, pmax); - SetLocalH (pmin, pmax, mparam.grading); + // SetLocalH (pmin, pmax, mparam.grading); + SetLocalH (pmin, pmax, grading); } PrintMessage (3, @@ -2675,7 +2655,7 @@ namespace netgen } - void Mesh :: CalcLocalHFromPointDistances(void) + void Mesh :: CalcLocalHFromPointDistances(double grading) { PrintMessage (3, "Calculating local h from point distances"); @@ -2684,7 +2664,8 @@ namespace netgen Point3d pmin, pmax; GetBox (pmin, pmax); - SetLocalH (pmin, pmax, mparam.grading); + // SetLocalH (pmin, pmax, mparam.grading); + SetLocalH (pmin, pmax, grading); } PointIndex i,j; @@ -2709,7 +2690,7 @@ namespace netgen } - void Mesh :: CalcLocalHFromSurfaceCurvature (double elperr) + void Mesh :: CalcLocalHFromSurfaceCurvature (double grading, double elperr) { PrintMessage (3, "Calculating local h from surface curvature"); @@ -2718,7 +2699,8 @@ namespace netgen Point3d pmin, pmax; GetBox (pmin, pmax); - SetLocalH (pmin, pmax, mparam.grading); + // SetLocalH (pmin, pmax, mparam.grading); + SetLocalH (pmin, pmax, grading); } @@ -3047,11 +3029,11 @@ namespace netgen - double Mesh :: ElementError (int eli) const + double Mesh :: ElementError (int eli, const MeshingParameters & mp) const { const Element & el = volelements.Get(eli); return CalcTetBadness (points.Get(el[0]), points.Get(el[1]), - points.Get(el[2]), points.Get(el[3]), -1); + points.Get(el[2]), points.Get(el[3]), -1, mp); } void Mesh :: AddLockedPoint (PointIndex pi) @@ -4145,11 +4127,11 @@ namespace netgen const int element, bool consider3D) const { - static Vec3d col1, col2, col3; - static Vec3d rhs, sol; + Vec3d col1, col2, col3; + Vec3d rhs, sol; const double eps = 1e-6; - static Array loctrigs; + Array loctrigs; //SZ @@ -4425,12 +4407,11 @@ namespace netgen double lami[3], const int element) const { - - static Vec3d col1, col2, col3; - static Vec3d rhs, sol; + Vec3d col1, col2, col3; + Vec3d rhs, sol; const double eps = 1.e-4; - static Array loctets; + Array loctets; VolumeElement(element).GetTets (loctets); @@ -5003,8 +4984,8 @@ namespace netgen void Mesh :: GetSurfaceElementsOfFace (int facenr, Array & sei) const { - static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace"); - NgProfiler::RegionTimer reg (timer); + static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace"); + NgProfiler::RegionTimer reg (timer); /* diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 39f52788..29bb5b32 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -311,7 +311,7 @@ namespace netgen /// - DLL_HEADER double ElementError (int eli) const; + DLL_HEADER double ElementError (int eli, const MeshingParameters & mp) const; /// DLL_HEADER void AddLockedPoint (PointIndex pi); @@ -382,7 +382,7 @@ namespace netgen */ DLL_HEADER double AverageH (int surfnr = 0) const; /// Calculates localh - DLL_HEADER void CalcLocalH (); + DLL_HEADER void CalcLocalH (double grading); /// DLL_HEADER void SetLocalH (const Point3d & pmin, const Point3d & pmax, double grading); /// @@ -391,9 +391,9 @@ namespace netgen DLL_HEADER void RestrictLocalHLine (const Point3d & p1, const Point3d & p2, double hloc); /// number of elements per radius - DLL_HEADER void CalcLocalHFromSurfaceCurvature(double elperr); + DLL_HEADER void CalcLocalHFromSurfaceCurvature(double grading, double elperr); /// - DLL_HEADER void CalcLocalHFromPointDistances(void); + DLL_HEADER void CalcLocalHFromPointDistances(double grading); /// DLL_HEADER void RestrictLocalH (resthtype rht, int nr, double loch); /// diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 6d2fe6df..c19fc795 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -22,7 +22,7 @@ namespace netgen Array connectednodes; - if (&mesh3d.LocalHFunction() == NULL) mesh3d.CalcLocalH(); + if (&mesh3d.LocalHFunction() == NULL) mesh3d.CalcLocalH(mp.grading); mesh3d.Compress(); diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp index 4d01ceef..495f8777 100644 --- a/libsrc/meshing/meshfunc2d.cpp +++ b/libsrc/meshing/meshfunc2d.cpp @@ -41,7 +41,7 @@ namespace netgen { MeshOptimize2d meshopt; meshopt.SetMetricWeight (1); - meshopt.ImproveMesh(mesh); + meshopt.ImproveMesh(mesh, mp); break; } diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp index 0aa6dae1..22ffc31f 100644 --- a/libsrc/meshing/meshing2.cpp +++ b/libsrc/meshing/meshing2.cpp @@ -7,24 +7,23 @@ namespace netgen // global variable for visualization +// static Array locpoints; +// static Array legalpoints; +// static Array plainpoints; +// static Array plainzones; +// static Array loclines; +// // static int geomtrig; +// //static const char * rname; +// static int cntelem, trials, nfaces; +// static int oldnl; +// static int qualclass; - static Array locpoints; - static Array legalpoints; - static Array plainpoints; - static Array plainzones; - static Array loclines; - // static int geomtrig; - //static const char * rname; - static int cntelem, trials, nfaces; - static int oldnl; - static int qualclass; - - Meshing2 :: Meshing2 (const Box<3> & aboundingbox) + Meshing2 :: Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox) { boundingbox = aboundingbox; - - LoadRules (NULL); + + LoadRules (NULL, mp.quad); // LoadRules ("rules/quad.rls"); // LoadRules ("rules/triangle.rls"); @@ -73,8 +72,8 @@ namespace netgen canuse = 0; ruleused = 0; - cntelem = 0; - trials = 0; + // cntelem = 0; + // trials = 0; } void Meshing2 :: EndMesh () @@ -194,7 +193,7 @@ namespace netgen - MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, double gh, int facenr) + MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr) { static int timer = NgProfiler::CreateTimer ("surface meshing"); @@ -227,6 +226,17 @@ namespace netgen double h, his, hshould; + Array locpoints; + Array legalpoints; + Array plainpoints; + Array plainzones; + Array loclines; + int cntelem = 0, trials = 0, nfaces = 0; + int oldnl = 0; + int qualclass; + + + // test for 3d overlaps Box3dTree surfeltree (boundingbox.PMin(), boundingbox.PMax()); @@ -419,7 +429,7 @@ namespace netgen //(*testout) << "locpoints " << locpoints << endl; - if (qualclass > mparam.giveuptol2d) + if (qualclass > mp.giveuptol2d) { PrintMessage (3, "give up with qualclass ", qualclass); PrintMessage (3, "number of frontlines = ", adfront->GetNFL()); @@ -704,7 +714,7 @@ namespace netgen - if (mparam.checkchartboundary) + if (mp.checkchartboundary) { for (int i = 1; i <= chartboundpoints.Size(); i++) { @@ -961,7 +971,7 @@ namespace netgen } - if (found && mparam.checkoverlap) + if (found && mp.checkoverlap) { // cout << "checkoverlap" << endl; // test for overlaps diff --git a/libsrc/meshing/meshing2.hpp b/libsrc/meshing/meshing2.hpp index f859395e..1171f97f 100644 --- a/libsrc/meshing/meshing2.hpp +++ b/libsrc/meshing/meshing2.hpp @@ -43,16 +43,16 @@ class Meshing2 public: /// - DLL_HEADER Meshing2 (const Box<3> & aboundingbox); + DLL_HEADER Meshing2 (const MeshingParameters & mp, const Box<3> & aboundingbox); /// DLL_HEADER virtual ~Meshing2 (); /// Load rules, either from file, or compiled rules - void LoadRules (const char * filename); + void LoadRules (const char * filename, bool quad); /// - DLL_HEADER MESHING2_RESULT GenerateMesh (Mesh & mesh, double gh, int facenr); + DLL_HEADER MESHING2_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr); DLL_HEADER void Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp); DLL_HEADER void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp); diff --git a/libsrc/meshing/meshtool.cpp b/libsrc/meshing/meshtool.cpp index 73d2ad8f..4a349ac9 100644 --- a/libsrc/meshing/meshtool.cpp +++ b/libsrc/meshing/meshtool.cpp @@ -192,7 +192,8 @@ namespace netgen // static double teterrpow = 2; double CalcTetBadness (const Point3d & p1, const Point3d & p2, - const Point3d & p3, const Point3d & p4, double h) + const Point3d & p3, const Point3d & p4, double h, + const MeshingParameters & mp) { double vol, l, ll, lll, ll1, ll2, ll3, ll4, ll5, ll6; double err; @@ -224,7 +225,7 @@ namespace netgen h * h * ( 1 / ll1 + 1 / ll2 + 1 / ll3 + 1 / ll4 + 1 / ll5 + 1 / ll6 ) - 12; - double teterrpow = mparam.opterrpow; + double teterrpow = mp.opterrpow; if(teterrpow < 1) teterrpow = 1; if (teterrpow == 1) return err; @@ -235,7 +236,8 @@ namespace netgen double CalcTetBadnessGrad (const Point3d & p1, const Point3d & p2, const Point3d & p3, const Point3d & p4, double h, - int pi, Vec<3> & grad) + int pi, Vec<3> & grad, + const MeshingParameters & mp) { double vol, l, ll, lll; double err; @@ -350,7 +352,7 @@ namespace netgen double errpow; - double teterrpow = mparam.opterrpow; + double teterrpow = mp.opterrpow; if(teterrpow < 1) teterrpow = 1; if (teterrpow == 1) diff --git a/libsrc/meshing/meshtool.hpp b/libsrc/meshing/meshtool.hpp index 6b7b7526..1b027222 100644 --- a/libsrc/meshing/meshtool.hpp +++ b/libsrc/meshing/meshtool.hpp @@ -46,18 +46,16 @@ extern int CheckCode (); /// -extern double CalcTetBadness (const Point3d & p1, - const Point3d & p2, - const Point3d & p3, - const Point3d & p4, - double h); +extern double CalcTetBadness (const Point3d & p1, const Point3d & p2, + const Point3d & p3, const Point3d & p4, + double h, + const MeshingParameters & mp = mparam); /// -extern double CalcTetBadnessGrad (const Point3d & p1, - const Point3d & p2, - const Point3d & p3, - const Point3d & p4, - double h, int pi, - Vec<3> & grad); +extern double CalcTetBadnessGrad (const Point3d & p1, const Point3d & p2, + const Point3d & p3, const Point3d & p4, + double h, int pi, + Vec<3> & grad, + const MeshingParameters & mp = mparam); /** Calculates volume of an element. diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index a8d40fc6..d279968c 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -452,7 +452,7 @@ namespace netgen DenseMatrix & trans) const { int np = GetNP(); - static DenseMatrix pmat(2, np), dshape(2, np); + DenseMatrix pmat(2, np), dshape(2, np); pmat.SetSize (2, np); dshape.SetSize (2, np); @@ -665,8 +665,8 @@ namespace netgen { int i, j; int nip = GetNIP(); - static DenseMatrix trans(2,2); - static DenseMatrix pmat; + DenseMatrix trans(2,2); + DenseMatrix pmat; pmat.SetSize (2, GetNP()); GetPointMatrix (points, pmat); @@ -783,8 +783,8 @@ namespace netgen } int nip = GetNIP(); - static DenseMatrix trans(2,2), dtrans(2,2); - static DenseMatrix pmat, vmat; + DenseMatrix trans(2,2), dtrans(2,2); + DenseMatrix pmat, vmat; pmat.SetSize (2, GetNP()); vmat.SetSize (2, GetNP()); @@ -846,8 +846,8 @@ namespace netgen { int i, j; int nip = GetNIP(); - static DenseMatrix trans(2,2); - static DenseMatrix pmat; + DenseMatrix trans(2,2); + DenseMatrix pmat; pmat.SetSize (2, GetNP()); @@ -1727,7 +1727,7 @@ namespace netgen DenseMatrix & trans) const { int np = GetNP(); - static DenseMatrix pmat(3, np), dshape(3, np); + DenseMatrix pmat(3, np), dshape(3, np); pmat.SetSize (3, np); dshape.SetSize (3, np); @@ -2037,8 +2037,8 @@ namespace netgen double Element :: CalcJacobianBadness (const T_POINTS & points) const { int nip = GetNIP(); - static DenseMatrix trans(3,3); - static DenseMatrix pmat; + DenseMatrix trans(3,3); + DenseMatrix pmat; pmat.SetSize (3, GetNP()); GetPointMatrix (points, pmat); @@ -2073,8 +2073,8 @@ namespace netgen { int i, j, k; int nip = GetNIP(); - static DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); - static DenseMatrix pmat, vmat; + DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); + DenseMatrix pmat, vmat; pmat.SetSize (3, GetNP()); vmat.SetSize (3, GetNP()); @@ -2148,8 +2148,8 @@ namespace netgen int pi, Vec<3> & grad) const { int nip = GetNIP(); - static DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); - static DenseMatrix pmat, vmat; + DenseMatrix trans(3,3), dtrans(3,3), hmat(3,3); + DenseMatrix pmat, vmat; pmat.SetSize (3, GetNP()); vmat.SetSize (3, GetNP()); diff --git a/libsrc/meshing/netrule3.cpp b/libsrc/meshing/netrule3.cpp index 03771409..de0e35e4 100644 --- a/libsrc/meshing/netrule3.cpp +++ b/libsrc/meshing/netrule3.cpp @@ -1,8 +1,6 @@ #include #include "meshing.hpp" -// #define MARK -// #include namespace netgen @@ -203,7 +201,7 @@ int vnetrule :: IsTriangleInFreeZone (const Point3d & p1, int infreeset, cannot = 0; - static Array pfi(3), pfi2(3); + ArrayMem pfi(3), pfi2(3); // convert from local index to freeset index int i, j; @@ -253,7 +251,7 @@ int vnetrule :: IsTriangleInFreeSet (const Point3d & p1, const Point3d & p2, double hpx, hpy, hpz, v1x, v1y, v1z, v2x, v2y, v2z; int act1, act2, act3, it; int cntout; - static Array activefaces; + Array activefaces; int isin; @@ -870,7 +868,7 @@ int vnetrule :: IsQuadInFreeZone (const Point3d & p1, int infreeset, cannot = 0; - static Array pfi(4), pfi2(4); + ArrayMem pfi(4), pfi2(4); // convert from local index to freeset index int i, j; @@ -933,7 +931,7 @@ int vnetrule :: IsQuadInFreeSet (const Point3d & p1, const Point3d & p2, return 1; } - static Array pi3(3); + ArrayMem pi3(3); int res; pi3.Elem(1) = pi.Get(1); diff --git a/libsrc/meshing/parser2.cpp b/libsrc/meshing/parser2.cpp index 92fd828d..55bd1a5d 100644 --- a/libsrc/meshing/parser2.cpp +++ b/libsrc/meshing/parser2.cpp @@ -493,7 +493,7 @@ void netrule :: LoadRule (istream & ist) extern const char * triarules[]; extern const char * quadrules[]; -void Meshing2 :: LoadRules (const char * filename) +void Meshing2 :: LoadRules (const char * filename, bool quad) { char buf[256]; istream * ist; @@ -520,7 +520,8 @@ void Meshing2 :: LoadRules (const char * filename) /* connect tetrules to one string */ const char ** hcp; - if (!mparam.quad) + // if (!mparam.quad) + if (!quad) { hcp = triarules; PrintMessage (3, "load internal triangle rules"); @@ -544,7 +545,8 @@ void Meshing2 :: LoadRules (const char * filename) tr1.reserve(len+1); - if (!mparam.quad) + // if (!mparam.quad) + if (!quad) hcp = triarules; else hcp = quadrules; diff --git a/libsrc/meshing/ruler3.cpp b/libsrc/meshing/ruler3.cpp index 69a63914..3b6f7a91 100644 --- a/libsrc/meshing/ruler3.cpp +++ b/libsrc/meshing/ruler3.cpp @@ -77,26 +77,25 @@ int Meshing3 :: ApplyRules int loktestmode; - static Array pused; // point is already mapped - static Array fused; // face is already mapped - static Array pmap; // map of reference point to local point - static Array pfixed; // point mapped by face-map - static Array fmapi; // face in reference is mapped to face nr ... - static Array fmapr; // face in reference is rotated to map - static Array transfreezone; // transformed free-zone - static int cnt = 0; - static INDEX_2_CLOSED_HASHTABLE ledges(100); // edges in local environment + Array pused; // point is already mapped + Array fused; // face is already mapped + Array pmap; // map of reference point to local point + Array pfixed; // point mapped by face-map + Array fmapi; // face in reference is mapped to face nr ... + Array fmapr; // face in reference is rotated to map + Array transfreezone; // transformed free-zone + INDEX_2_CLOSED_HASHTABLE ledges(100); // edges in local environment - static Array tempnewpoints; - static Array tempnewfaces; - static Array tempdelfaces; - static Array tempelements; - static Array triboxes; // bounding boxes of local faces + Array tempnewpoints; + Array tempnewfaces; + Array tempdelfaces; + Array tempelements; + Array triboxes; // bounding boxes of local faces + Array pnearness; + Array fnearness; - static Array pnearness; - static Array fnearness; - + static int cnt = 0; cnt++; delfaces.SetSize (0); diff --git a/libsrc/meshing/smoothing2.cpp b/libsrc/meshing/smoothing2.cpp index 0d25efbd..7bcbaffd 100644 --- a/libsrc/meshing/smoothing2.cpp +++ b/libsrc/meshing/smoothing2.cpp @@ -8,84 +8,6 @@ namespace netgen static const MeshOptimize2d * meshthis; - -#ifdef OLD - - void CalcTriangleBadness (double x2, double x3, double y3, double metricweight, - double h, double & badness, double & g1x, double & g1y) - { - // badness = sqrt(3.0) /36 * circumference^2 / area - 1 - // p1 = (0, 0), p2 = (x2, 0), p3 = (x3, y3); - - Vec2d v23; - double l12, l13, l23, cir, area; - static const double c = sqrt(3.0) / 36; - double c1, c2, c3, c4; - - v23.X() = x3 - x2; - v23.Y() = y3; - - l12 = x2; - l13 = sqrt (x3*x3 + y3*y3); - l23 = v23.Length(); - - cir = l12 + l13 + l23; - area = 0.5 * x2 * y3; - - if (area <= 1e-24 * cir * cir) - { - g1x = 0; - g1y = 0; - badness = 1e10; - return; - } - - badness = c * cir * cir / area - 1; - - c1 = 2 * c * cir / area; - c2 = 0.5 * c * cir * cir / (area * area); - - g1x = c1 * ( - 1 - x3 / l13) - c2 * (-v23.Y()); - g1y = c1 * ( - y3 / l13) - c2 * ( v23.X()); - - // metricweight = 0.1; - if (metricweight > 0) - { - // area = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1); - // add: metricweight * (area / h^2 + h^2 / area - 2) - - const double area = x2 * y3; - const double dareax1 = -y3; - const double dareay1 = x3 - x2; - - const double areahh = area / (h * h); - const double fac = metricweight * (areahh - 1 / areahh) / area; - - badness += metricweight * (areahh + 1 / areahh - 2); - g1x += fac * dareax1; - g1y += fac * dareay1; - - /* - // add: metricweight * (l1^2/h^2 + l2^2/h^2 + l3^2/h2 + h^2/l1^2 + h^2/l2^2 + h^2/l3^2 - 6) - double h2 = h*h; - double l1 = x2*x2; - double l2 = x3*x3+y3*y3; - double l3 = (x2-x3)*(x2-x3)+y3*y3; - double dl1dx = 2*(-x2); - double dl1dy = 0; - double dl2dx = -2*x3; - double dl2dy = -2*y3; - - badness += (l1/h2 + l2/h2 + l3/h2 +h2/l1 + h2/l2 + h2/l3-6) * metricweight; - - g1x += metricweight * (dl1dx/h2-h2/(l1*l1)*dl1dx + dl2dx/h2-h2/(l2*l2)*dl2dx); - g1y += metricweight * (dl1dy/h2-h2/(l1*l1)*dl1dy + dl2dy/h2-h2/(l2*l2)*dl2dy); - */ - } - } - -#endif - static const double c_trig = 0.14433756; // sqrt(3.0) / 12 static const double c_trig4 = 0.57735026; // sqrt(3.0) / 3 @@ -162,38 +84,6 @@ namespace netgen - - - - -#ifdef OLD - double CalcTriangleBadness (const Point3d & p1, - const Point3d & p2, - const Point3d & p3, - double metricweight, - double h) - { - double badness; - double g1x, g1y; - - Vec3d e1 (p1, p2); - Vec3d e2 (p1, p3); - - double e1l = e1.Length() + 1e-24; - e1 /= e1l; - double e1e2 = (e1 * e2); - e2.Add (-e1e2, e1); - double e2l = e2.Length(); - - CalcTriangleBadness ( e1l, e1e2, e2l, - metricweight, h, badness, g1x, g1y); - return badness; - } -#endif - - - - double CalcTriangleBadness (const Point3d & p1, const Point3d & p2, const Point3d & p3, @@ -274,6 +164,7 @@ namespace netgen class Opti2SurfaceMinFunction : public MinFunction { const Mesh & mesh; + public: Opti2SurfaceMinFunction (const Mesh & amesh) : mesh(amesh) @@ -416,6 +307,7 @@ namespace netgen class Opti2EdgeMinFunction : public MinFunction { const Mesh & mesh; + public: Opti2EdgeMinFunction (const Mesh & amesh) : mesh(amesh) { } ; @@ -650,7 +542,7 @@ namespace netgen ; } - void MeshOptimize2d :: ImproveMesh (Mesh & mesh) + void MeshOptimize2d :: ImproveMesh (Mesh & mesh, const MeshingParameters & mp) { if (!faceindex) { @@ -658,7 +550,7 @@ namespace netgen for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++) { - ImproveMesh (mesh); + ImproveMesh (mesh, mp); if (multithread.terminate) throw NgException ("Meshing stopped"); } @@ -702,7 +594,7 @@ namespace netgen Array savepoints(mesh.GetNP()); - uselocalh = mparam.uselocalh; + uselocalh = mp.uselocalh; Array nelementsonpoint(mesh.GetNP()); @@ -724,7 +616,7 @@ namespace netgen } - loch = mparam.maxh; + loch = mp.maxh; locmetricweight = metricweight; meshthis = this; diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 94df4d03..47da61cc 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -23,7 +23,6 @@ namespace netgen void MinFunctionSum :: Grad (const Vector & x, Vector & g) const { g = 0.; - // static Vector gi(3); VectorMem<3> gi; for(int i=0; i gi; for(int i=0; iFuncGrad(x,gi); @@ -121,10 +120,6 @@ namespace netgen double PointFunction1 :: FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const { - /* - static Vector hx(3); - static double eps = 1e-6; - */ VectorMem<3> hx; const double eps = 1e-6; @@ -150,8 +145,8 @@ namespace netgen double PointFunction1 :: FuncGrad (const Vector & x, Vector & g) const { - static Vector hx(3); - static double eps = 1e-6; + VectorMem<3> hx; + double eps = 1e-6; hx = x; for (int i = 0; i < 3; i++) @@ -247,9 +242,8 @@ namespace netgen int i; double badness = 0; - static Vector hv(4); - static Vector res; - res.SetSize (m.Height()); + VectorMem<4> hv; + Vector res(m.Height()); for (i = 0;i < 3; i++) hv(i) = vp(i); @@ -270,8 +264,8 @@ namespace netgen double CheapPointFunction1 :: FuncGrad (const Vector & x, Vector & g) const { - static Vector hx(3); - static double eps = 1e-6; + VectorMem<3> hx; + double eps = 1e-6; hx = x; for (int i = 0; i < 3; i++) @@ -546,8 +540,8 @@ namespace netgen double CheapPointFunction :: PointFunctionValue (const Point<3> & pp) const { - static Vector p4(4); - static Vector di; + VectorMem<4> p4; + Vector di; int n = m.Height(); p4(0) = pp(0); @@ -574,8 +568,8 @@ namespace netgen double CheapPointFunction :: PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const { - static Vector p4(4); - static Vector di; + VectorMem<4> p4; + Vector di; int n = m.Height(); @@ -686,7 +680,7 @@ namespace netgen { int n = x.Size(); - static Vector hx; + Vector hx; hx.SetSize(n); double eps = 1e-8; @@ -800,7 +794,7 @@ namespace netgen { Vec3d n, vgrad; Point3d pp1; - static Vector freegrad(3); + VectorMem<3> freegrad; CalcNewPoint (x, pp1); @@ -893,7 +887,7 @@ double Opti3EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const Vec3d n1, n2, v1, vgrad; Point3d pp1; double badness; - static Vector freegrad(3); + VectorMem<3> freegrad; CalcNewPoint (x, pp1); @@ -923,7 +917,8 @@ double Opti3EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const double CalcTotalBad (const Mesh::T_POINTS & points, - const Mesh::T_VOLELEMENTS & elements) + const Mesh::T_VOLELEMENTS & elements, + const MeshingParameters & mp) { double sum = 0; double elbad; @@ -931,7 +926,7 @@ double CalcTotalBad (const Mesh::T_POINTS & points, tets_in_qualclass.SetSize(20); tets_in_qualclass = 0; - double teterrpow = mparam.opterrpow; + double teterrpow = mp.opterrpow; for (int i = 1; i <= elements.Size(); i++) { diff --git a/libsrc/meshing/specials.cpp b/libsrc/meshing/specials.cpp index 977b08f9..a227e268 100644 --- a/libsrc/meshing/specials.cpp +++ b/libsrc/meshing/specials.cpp @@ -156,7 +156,7 @@ void CutOffAndCombine (Mesh & mesh, const Mesh & othermesh) mesh.AddLockedPoint (pmat.Elem(i)); mesh.CalcSurfacesOfNode(); - mesh.CalcLocalH(); + mesh.CalcLocalH(mparam.grading); } diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index a73def71..2a42352b 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -112,8 +112,8 @@ namespace netgen gp_Pnt2d parmid; - parmid.SetX(0.3*(par0.X()+par1.X()+par2.X())); - parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y())); + parmid.SetX( (par0.X()+par1.X()+par2.X()) / 3 ); + parmid.SetY( (par0.Y()+par1.Y()+par2.Y()) / 3 ); if (depth%3 == 0) { diff --git a/libsrc/stlgeom/meshstlsurface.cpp b/libsrc/stlgeom/meshstlsurface.cpp index 137b225e..fcbb1e02 100644 --- a/libsrc/stlgeom/meshstlsurface.cpp +++ b/libsrc/stlgeom/meshstlsurface.cpp @@ -320,7 +320,7 @@ int STLSurfaceMeshing (STLGeometry & geom, mesh.CalcSurfacesOfNode(); optmesh.EdgeSwapping (mesh, 0); mesh.CalcSurfacesOfNode(); - optmesh.ImproveMesh (mesh); + optmesh.ImproveMesh (mesh, mparam); } mesh.Compress(); @@ -606,7 +606,7 @@ void STLSurfaceMeshing1 (STLGeometry & geom, return; PrintMessage(5,"Meshing surface ", fnr, "/", mesh.GetNFD()); - MeshingSTLSurface meshing (geom); + MeshingSTLSurface meshing (geom, mparam); meshing.SetStartTime (starttime); @@ -666,7 +666,7 @@ void STLSurfaceMeshing1 (STLGeometry & geom, /* (*testout) << "start meshing with h = " << h << endl; */ - meshing.GenerateMesh (mesh, h, fnr); // face index + meshing.GenerateMesh (mesh, mparam, h, fnr); // face index extern void Render(); Render(); @@ -727,7 +727,7 @@ void STLSurfaceOptimization (STLGeometry & geom, } case 'm': { - optmesh.ImproveMesh(mesh); + optmesh.ImproveMesh(mesh, mparam); break; } case 'c': @@ -749,8 +749,9 @@ void STLSurfaceOptimization (STLGeometry & geom, -MeshingSTLSurface :: MeshingSTLSurface (STLGeometry & ageom) - : Meshing2(ageom.GetBoundingBox()), geom(ageom) +MeshingSTLSurface :: MeshingSTLSurface (STLGeometry & ageom, + const MeshingParameters & mp) + : Meshing2(mp, ageom.GetBoundingBox()), geom(ageom) { ; } diff --git a/libsrc/stlgeom/meshstlsurface.hpp b/libsrc/stlgeom/meshstlsurface.hpp index be15dbb3..4e678439 100644 --- a/libsrc/stlgeom/meshstlsurface.hpp +++ b/libsrc/stlgeom/meshstlsurface.hpp @@ -23,7 +23,7 @@ class MeshingSTLSurface : public Meshing2 int transformationtrig; public: /// - MeshingSTLSurface (STLGeometry & ageom); + MeshingSTLSurface (STLGeometry & ageom, const MeshingParameters & mp); protected: /// diff --git a/libsrc/stlgeom/stlgeommesh.cpp b/libsrc/stlgeom/stlgeommesh.cpp index ad430305..c20c6acf 100644 --- a/libsrc/stlgeom/stlgeommesh.cpp +++ b/libsrc/stlgeom/stlgeommesh.cpp @@ -1451,7 +1451,8 @@ int STLMeshingDummy (STLGeometry* stlgeometry, Mesh*& mesh, MeshingParameters & stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10), mparam.grading); mesh -> LoadLocalMeshSize (mparam.meshsizefilename); - mesh -> CalcLocalHFromSurfaceCurvature (stlparam.resthsurfmeshcurvfac); + mesh -> CalcLocalHFromSurfaceCurvature (mparam.grading, + stlparam.resthsurfmeshcurvfac); mparam.optimize2d = "cmsmSm"; STLSurfaceOptimization (*stlgeometry, *mesh, mparam); #ifdef STAT_STREAM @@ -1502,7 +1503,7 @@ int STLMeshingDummy (STLGeometry* stlgeometry, Mesh*& mesh, MeshingParameters & stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10), mparam.grading); mesh -> LoadLocalMeshSize (mparam.meshsizefilename); - mesh -> CalcLocalH (); + mesh -> CalcLocalH (mparam.grading); } diff --git a/libsrc/stlgeom/stlpkg.cpp b/libsrc/stlgeom/stlpkg.cpp index 79c416dc..4282633a 100644 --- a/libsrc/stlgeom/stlpkg.cpp +++ b/libsrc/stlgeom/stlpkg.cpp @@ -571,7 +571,8 @@ namespace netgen stlgeometry -> RestrictLocalH(*mesh, mparam.maxh); if (stlparam.resthsurfmeshcurvenable) - mesh -> CalcLocalHFromSurfaceCurvature (stlparam.resthsurfmeshcurvfac); + mesh -> CalcLocalHFromSurfaceCurvature (mparam.grading, + stlparam.resthsurfmeshcurvfac); } return TCL_OK; diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index f15ae281..b989cd03 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -4454,7 +4454,7 @@ namespace netgen int argc, tcl_const char *argv[]) { int i; - static char buf[1000]; + char buf[1000]; buf[0] = 0; if (argc >= 2) diff --git a/ng/ngpkg.cpp b/ng/ngpkg.cpp index 7a5d3b3f..26a8de6f 100644 --- a/ng/ngpkg.cpp +++ b/ng/ngpkg.cpp @@ -379,7 +379,7 @@ namespace netgen mesh->GetNE(), " Elements."); mesh->SetGlobalH (mparam.maxh); - mesh->CalcLocalH(); + mesh->CalcLocalH(mparam.grading); return TCL_OK; } @@ -1024,7 +1024,7 @@ namespace netgen } mesh->SetGlobalH (mparam.maxh); - mesh->CalcLocalH(); + mesh->CalcLocalH(mparam.grading); return TCL_OK; } @@ -1083,7 +1083,7 @@ namespace netgen Mesh othermesh; othermesh.Load (argv[1]); othermesh.SetGlobalH (mparam.maxh); - othermesh.CalcLocalH(); + othermesh.CalcLocalH(mparam.grading); CutOffAndCombine (*mesh, othermesh); return TCL_OK; @@ -1644,7 +1644,7 @@ namespace netgen */ if(!mesh->LocalHFunctionGenerated()) - mesh->CalcLocalH(); + mesh->CalcLocalH(mparam.grading); mesh->LocalHFunction().SetGrading (mparam.grading); ref . Bisect (*mesh, biopt); diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index ecde9c91..de46b2ab 100644 --- a/nglib/nglib.cpp +++ b/nglib/nglib.cpp @@ -366,7 +366,7 @@ namespace nglib //MeshingParameters mparam; mp->Transfer_Parameters(); - m->CalcLocalH(); + m->CalcLocalH(mparam.grading); MeshVolume (mparam, *m); RemoveIllegalElements (*m);