diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index b792ce39..5897a986 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -227,8 +227,13 @@ namespace netgen mparam.checkoverlap = 0; - - meshing.GenerateMesh (*mesh, h, domnr); + + /* + if (!mparam.quad) + meshing.Delaunay (*mesh, domnr, mparam); + else + */ + meshing.GenerateMesh (*mesh, h, domnr); for (SurfaceElementIndex sei = oldnf; sei < mesh->GetNSE(); sei++) (*mesh)[sei].SetIndex (domnr); diff --git a/libsrc/gprim/geomobjects.hpp b/libsrc/gprim/geomobjects.hpp index 85059c13..2d86f3c3 100644 --- a/libsrc/gprim/geomobjects.hpp +++ b/libsrc/gprim/geomobjects.hpp @@ -197,6 +197,14 @@ protected: Point pmin, pmax; public: Box () { ; } + + Box ( const Point & p1) + { + for (int i = 0; i < D; i++) + pmin(i) = pmax(i) = p1(i); + } + + Box ( const Point & p1, const Point & p2) { for (int i = 0; i < D; i++) diff --git a/libsrc/interface/readuser.cpp b/libsrc/interface/readuser.cpp index 98f325d9..bd711264 100644 --- a/libsrc/interface/readuser.cpp +++ b/libsrc/interface/readuser.cpp @@ -43,6 +43,7 @@ namespace netgen { Point3d p; in >> p.X() >> p.Y() >> p.Z(); + p.Z() *= 10; mesh.AddPoint (p); } diff --git a/libsrc/linalg/densemat.cpp b/libsrc/linalg/densemat.cpp index 768c7d5c..df5ead6a 100644 --- a/libsrc/linalg/densemat.cpp +++ b/libsrc/linalg/densemat.cpp @@ -262,6 +262,7 @@ namespace netgen if (det == 0) { (*myerr) << "CalcInverse: Matrix singular" << endl; + (*testout) << "CalcInverse: Matrix singular" << endl; return; } @@ -488,6 +489,7 @@ namespace netgen if (max < 1e-20) { cerr << "Inverse matrix: matrix singular" << endl; + *testout << "Inverse matrix: matrix singular" << endl; return; } diff --git a/libsrc/meshing/Makefile.am b/libsrc/meshing/Makefile.am index 1d0a9be6..e673802a 100644 --- a/libsrc/meshing/Makefile.am +++ b/libsrc/meshing/Makefile.am @@ -15,14 +15,15 @@ clusters.hpp hprefinement.hpp improve3.hpp meshtype.hpp \ METASOURCES = AUTO noinst_LTLIBRARIES = libmesh.la -libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp \ - clusters.cpp curvedelems.cpp delaunay.cpp geomsearch.cpp global.cpp \ - hprefinement.cpp improve2.cpp improve2gen.cpp improve3.cpp localh.cpp \ - meshclass.cpp meshfunc.cpp meshfunc2d.cpp meshing2.cpp meshing3.cpp \ - meshtool.cpp meshtype.cpp msghandler.cpp netrule2.cpp \ - netrule3.cpp parser2.cpp parser3.cpp prism2rls.cpp \ - pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp \ - ruler2.cpp ruler3.cpp secondorder.cpp smoothing2.5.cpp \ - smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp \ - topology.cpp triarls.cpp validate.cpp zrefine.cpp bcfunctions.cpp \ +libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp \ + clusters.cpp curvedelems.cpp delaunay.cpp delaunay2d.cpp \ + geomsearch.cpp global.cpp hprefinement.cpp improve2.cpp \ + improve2gen.cpp improve3.cpp localh.cpp meshclass.cpp \ + meshfunc.cpp meshfunc2d.cpp meshing2.cpp meshing3.cpp \ + meshtool.cpp meshtype.cpp msghandler.cpp netrule2.cpp \ + netrule3.cpp parser2.cpp parser3.cpp prism2rls.cpp \ + pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp \ + ruler2.cpp ruler3.cpp secondorder.cpp smoothing2.5.cpp \ + smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp \ + topology.cpp triarls.cpp validate.cpp zrefine.cpp bcfunctions.cpp \ parallelmesh.cpp paralleltop.cpp paralleltop.hpp basegeom.cpp diff --git a/libsrc/meshing/adfront2.cpp b/libsrc/meshing/adfront2.cpp index 3b273bc7..4eea0da1 100644 --- a/libsrc/meshing/adfront2.cpp +++ b/libsrc/meshing/adfront2.cpp @@ -8,8 +8,8 @@ namespace netgen { - AdFront2::FrontPoint2 :: FrontPoint2 (const Point<3> & ap, PointIndex agi, - MultiPointGeomInfo * amgi, bool aonsurface) + FrontPoint2 :: FrontPoint2 (const Point<3> & ap, PointIndex agi, + MultiPointGeomInfo * amgi, bool aonsurface) { p = ap; globalindex = agi; @@ -464,4 +464,45 @@ namespace netgen ost << flush; } + + + bool AdFront2 :: Inside (const Point<2> & p) const + { + int cnt; + Vec<2> n; + Vec<3> v1; + DenseMatrix a(2), ainv(2); + Vector b(2), u(2); + + // random numbers: + n(0) = 0.123871; + n(1) = 0.15432; + + cnt = 0; + for (int i = 0; i < lines.Size(); i++) + if (lines[i].Valid()) + { + const Point<3> & p1 = points[lines[i].L().I1()].P(); + const Point<3> & p2 = points[lines[i].L().I2()].P(); + + v1 = p2 - p1; + + a(0, 0) = v1(0); + a(1, 0) = v1(1); + + a(0, 1) = -n(0); + a(1, 1) = -n(1); + + b(0) = p(0) - p1(0); + b(1) = p(1) - p1(1); + + CalcInverse (a, ainv); + ainv.Mult (b, u); + + if (u(0) >= 0 && u(0) <= 1 && u(1) > 0) + cnt++; + } + + return ((cnt % 2) != 0); + } } diff --git a/libsrc/meshing/adfront2.hpp b/libsrc/meshing/adfront2.hpp index 3f71ff7c..fcaffd60 100644 --- a/libsrc/meshing/adfront2.hpp +++ b/libsrc/meshing/adfront2.hpp @@ -13,8 +13,6 @@ Advancing front class for surfaces */ -class AdFront2 -{ /// class FrontPoint2 @@ -163,6 +161,8 @@ class AdFront2 }; +class AdFront2 +{ /// Array points; /// front points @@ -202,6 +202,11 @@ public: } /// int GetNFL () const { return nfl; } + + const FrontLine & GetLine (int nr) { return lines[nr]; } + const FrontPoint2 & GetPoint (int nr) { return points[nr]; } + + /// int SelectBaseLine (Point<3> & p1, Point<3> & p2, const PointGeomInfo *& geominfo1, @@ -250,6 +255,18 @@ public: { return points[pi].GlobalIndex(); } + + + /// is Point p inside Surface (flat geometry only) + bool Inside (const Point<2> & p) const; + + bool SameSide (const Point<2> & lp1, const Point<2> & lp2, + const Array * testfaces = NULL) const + { + return Inside (lp1) == Inside (lp2); + } + + /// void SetStartFront (); /// diff --git a/libsrc/meshing/adfront3.cpp b/libsrc/meshing/adfront3.cpp index c5dbb8d4..99fd6e20 100644 --- a/libsrc/meshing/adfront3.cpp +++ b/libsrc/meshing/adfront3.cpp @@ -772,7 +772,7 @@ void AdFront3 :: SetStartFront (int /* baseelnp */) bool AdFront3 :: Inside (const Point<3> & p) const { - int i, cnt; + int cnt; Vec3d n, v1, v2; DenseMatrix a(3), ainv(3); Vector b(3), u(3); @@ -783,7 +783,7 @@ bool AdFront3 :: Inside (const Point<3> & p) const n.Z() = -0.43989; cnt = 0; - for (i = 1; i <= faces.Size(); i++) + for (int i = 1; i <= faces.Size(); i++) if (faces.Get(i).Valid()) { const Point<3> & p1 = points[faces.Get(i).Face().PNum(1)].P(); @@ -827,36 +827,29 @@ bool AdFront3 :: Inside (const Point<3> & p) const int AdFront3 :: SameSide (const Point<3> & lp1, const Point<3> & lp2, const Array * testfaces) const { - int i, ii, cnt; - const Point<3> *line[2]; line[0] = &lp1; line[1] = &lp2; - cnt = 0; - Point3d pmin(lp1); Point3d pmax(lp1); pmin.SetToMin (lp2); pmax.SetToMax (lp2); - - static Array aprif; + + ArrayMem aprif; aprif.SetSize(0); if (!testfaces) facetree->GetIntersecting (pmin, pmax, aprif); else - { - for (i = 1; i <= testfaces->Size(); i++) - aprif.Append (testfaces->Get(i)); - } + for (int i = 1; i <= testfaces->Size(); i++) + aprif.Append (testfaces->Get(i)); - // (*testout) << "test ss, p1,p2 = " << lp1 << lp2 << ", inters = " << aprif.Size() << endl; - // for (i = 1; i <= faces.Size(); i++) - for (ii = 1; ii <= aprif.Size(); ii++) + int cnt = 0; + for (int ii = 1; ii <= aprif.Size(); ii++) { - i = aprif.Get(ii); + int i = aprif.Get(ii); if (faces.Get(i).Valid()) { @@ -867,7 +860,6 @@ int AdFront3 :: SameSide (const Point<3> & lp1, const Point<3> & lp2, if (IntersectTriangleLine (&tri[0], &line[0])) cnt++; - } } diff --git a/libsrc/meshing/delaunay2d.cpp b/libsrc/meshing/delaunay2d.cpp new file mode 100644 index 00000000..e92ef963 --- /dev/null +++ b/libsrc/meshing/delaunay2d.cpp @@ -0,0 +1,174 @@ +#include +#include "meshing.hpp" + +// not yet working .... + +namespace netgen +{ + + void Meshing2 :: BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp) + { + double filldist = mp.filldist; + + cout << "blockfill local h" << endl; + cout << "rel filldist = " << filldist << endl; + PrintMessage (3, "blockfill local h"); + + Array > npoints; + + // adfront -> CreateTrees(); + + Box<3> bbox ( Box<3>::EMPTY_BOX ); + double maxh = 0; + + for (int i = 0; i < adfront->GetNFL(); 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()); + + double hi = Dist (p1, p2); + if (hi > maxh) maxh = hi; + + bbox.Add (p1); + bbox.Add (p2); + } + + + cout << "bbox = " << bbox << endl; + + + Point<3> mpc = bbox.Center(); + bbox.Increase (bbox.Diam()/2); + Box<3> meshbox = bbox; + + LocalH loch2 (bbox, 1); + + if (mp.maxh < maxh) maxh = mp.maxh; + + bool changed; + do + { + mesh.LocalHFunction().ClearFlags(); + + for (int i = 0; i < adfront->GetNFL(); i++) + { + const FrontLine & line = adfront->GetLine(i); + + Box<3> bbox (adfront->GetPoint (line.L().I1())); + bbox.Add (adfront->GetPoint (line.L().I2())); + + + double filld = filldist * bbox.Diam(); + bbox.Increase (filld); + + mesh.LocalHFunction().CutBoundary (bbox); + } + + + mesh.LocalHFunction().FindInnerBoxes (adfront, NULL); + + npoints.SetSize(0); + mesh.LocalHFunction().GetInnerPoints (npoints); + + changed = false; + for (int i = 0; i < npoints.Size(); i++) + { + if (mesh.LocalHFunction().GetH(npoints[i]) > 1.5 * maxh) + { + mesh.LocalHFunction().SetH (npoints[i], maxh); + changed = true; + } + } + } + while (changed); + + if (debugparam.slowchecks) + (*testout) << "Blockfill with points: " << endl; + *testout << "loch = " << mesh.LocalHFunction() << endl; + + *testout << "npoints = " << endl << npoints << endl; + + for (int i = 1; i <= npoints.Size(); i++) + { + if (meshbox.IsIn (npoints.Get(i))) + { + int gpnum = mesh.AddPoint (npoints.Get(i)); + 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)) + { + cout << "add outside point" << endl; + (*testout) << "outside" << endl; + } + } + + } + } + + + + // find outer points + + loch2.ClearFlags(); + + for (int i = 0; i < adfront->GetNFL(); i++) + { + const FrontLine & line = adfront->GetLine(i); + + 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++) + { + const FrontLine & line = adfront->GetLine(i); + + 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); + + npoints.SetSize(0); + loch2.GetOuterPoints (npoints); + + for (int i = 1; i <= npoints.Size(); i++) + { + if (meshbox.IsIn (npoints.Get(i))) + { + int gpnum = mesh.AddPoint (npoints.Get(i)); + adfront->AddPoint (npoints.Get(i), gpnum); + } + } + + } + + void Meshing2 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp) + { + cout << "2D Delaunay meshing (in progress)" << endl; + + int oldnp = mesh.GetNP(); + + cout << "np, old = " << mesh.GetNP() << endl; + + BlockFillLocalH (mesh, mp); + + + cout << "np, now = " << mesh.GetNP() << endl; + + } + +} diff --git a/libsrc/meshing/localh.cpp b/libsrc/meshing/localh.cpp index ba313f22..4138fbf4 100644 --- a/libsrc/meshing/localh.cpp +++ b/libsrc/meshing/localh.cpp @@ -5,453 +5,446 @@ namespace netgen { -GradingBox :: GradingBox (const double * ax1, const double * ax2) -{ - h2 = 0.5 * (ax2[0] - ax1[0]); - for (int i = 0; i <= 2; i++) - { - /* - x1[i] = ax1[i]; - x2[i] = ax2[i]; - */ + GradingBox :: GradingBox (const double * ax1, const double * ax2) + { + h2 = 0.5 * (ax2[0] - ax1[0]); + for (int i = 0; i < 3; i++) xmid[i] = 0.5 * (ax1[i] + ax2[i]); - } - /* - (*testout) << "new box: " << xmid[0] << "-" << xmid[1] << "-" << xmid[2] - << " h = " << (x2[0] - x1[0]) << endl; - */ + for (int i = 0; i < 8; i++) + childs[i] = NULL; + father = NULL; - for (int i = 0; i < 8; i++) - childs[i] = NULL; - father = NULL; + flags.cutboundary = 0; + flags.isinner = 0; + flags.oldcell = 0; + flags.pinner = 0; - flags.cutboundary = 0; - flags.isinner = 0; - flags.oldcell = 0; - flags.pinner = 0; - - // hopt = x2[0] - x1[0]; - hopt = 2 * h2; -} + hopt = 2 * h2; + } + BlockAllocator GradingBox :: ball(sizeof (GradingBox)); -BlockAllocator GradingBox :: ball(sizeof (GradingBox)); + void * GradingBox :: operator new(size_t) + { + return ball.Alloc(); + } -void * GradingBox :: operator new(size_t) -{ - return ball.Alloc(); -} - -void GradingBox :: operator delete (void * p) -{ - ball.Free (p); -} + void GradingBox :: operator delete (void * p) + { + ball.Free (p); + } + void GradingBox :: DeleteChilds() + { + for (int i = 0; i < 8; i++) + if (childs[i]) + { + childs[i]->DeleteChilds(); + delete childs[i]; + childs[i] = NULL; + } + } + + LocalH :: LocalH (const Point3d & pmin, const Point3d & pmax, double agrading) + { + double x1[3], x2[3]; + double hmax; -void GradingBox :: DeleteChilds() -{ - int i; - for (i = 0; i < 8; i++) - if (childs[i]) + boundingbox = Box3d (pmin, pmax); + grading = agrading; + + // a small enlargement, non-regular points + double val = 0.0879; + for (int i = 1; i <= 3; i++) { - childs[i]->DeleteChilds(); - delete childs[i]; - childs[i] = NULL; + x1[i-1] = (1 + val * i) * pmin.X(i) - val * i * pmax.X(i); + x2[i-1] = 1.1 * pmax.X(i) - 0.1 * pmin.X(i); } -} + + hmax = x2[0] - x1[0]; + for (int i = 1; i <= 2; i++) + if (x2[i] - x1[i] > hmax) + hmax = x2[i] - x1[i]; + + for (int i = 0; i <= 2; i++) + x2[i] = x1[i] + hmax; + + root = new GradingBox (x1, x2); + boxes.Append (root); + } -LocalH :: LocalH (const Point3d & pmin, const Point3d & pmax, double agrading) -{ - double x1[3], x2[3]; - double hmax; - int i; + LocalH :: LocalH (const Box<3> & box, double agrading) + { + Point3d pmin = box.PMin(); + Point3d pmax = box.PMax(); - boundingbox = Box3d (pmin, pmax); - grading = agrading; + double x1[3], x2[3]; + double hmax; - // a small enlargement, non-regular points - double val = 0.0879; - for (i = 1; i <= 3; i++) - { + boundingbox = Box3d (pmin, pmax); + grading = agrading; - x1[i-1] = (1 + val * i) * pmin.X(i) - val * i * pmax.X(i); - x2[i-1] = 1.1 * pmax.X(i) - 0.1 * pmin.X(i); - } + // a small enlargement, non-regular points + double val = 0.0879; + for (int i = 1; i <= 3; i++) + { + x1[i-1] = (1 + val * i) * pmin.X(i) - val * i * pmax.X(i); + x2[i-1] = 1.1 * pmax.X(i) - 0.1 * pmin.X(i); + } - hmax = x2[0] - x1[0]; - for (i = 1; i <= 2; i++) - if (x2[i] - x1[i] > hmax) - hmax = x2[i] - x1[i]; + hmax = x2[0] - x1[0]; + for (int i = 1; i <= 2; i++) + if (x2[i] - x1[i] > hmax) + hmax = x2[i] - x1[i]; - for (i = 0; i <= 2; i++) - x2[i] = x1[i] + hmax; + for (int i = 0; i <= 2; i++) + x2[i] = x1[i] + hmax; - root = new GradingBox (x1, x2); - boxes.Append (root); -} + root = new GradingBox (x1, x2); + boxes.Append (root); + } -LocalH :: ~LocalH () -{ - root->DeleteChilds(); - delete root; -} -void LocalH :: Delete () -{ - root->DeleteChilds(); -} -void LocalH :: SetH (const Point3d & p, double h) -{ - /* - (*testout) << "Set h at " << p << " to " << h << endl; - if (h < 1e-8) - { + + LocalH :: ~LocalH () + { + root->DeleteChilds(); + delete root; + } + + void LocalH :: Delete () + { + root->DeleteChilds(); + } + + void LocalH :: SetH (const Point3d & p, double h) + { + /* + (*testout) << "Set h at " << p << " to " << h << endl; + if (h < 1e-8) + { cout << "do not set h to " << h << endl; return; - } - */ + } + */ - if (fabs (p.X() - root->xmid[0]) > root->h2 || - fabs (p.Y() - root->xmid[1]) > root->h2 || - fabs (p.Z() - root->xmid[2]) > root->h2) - return; + if (fabs (p.X() - root->xmid[0]) > root->h2 || + fabs (p.Y() - root->xmid[1]) > root->h2 || + fabs (p.Z() - root->xmid[2]) > root->h2) + return; - /* - if (p.X() < root->x1[0] || p.X() > root->x2[0] || - p.Y() < root->x1[1] || p.Y() > root->x2[1] || - p.Z() < root->x1[2] || p.Z() > root->x2[2]) - return; - */ + /* + if (p.X() < root->x1[0] || p.X() > root->x2[0] || + p.Y() < root->x1[1] || p.Y() > root->x2[1] || + p.Z() < root->x1[2] || p.Z() > root->x2[2]) + return; + */ - if (GetH(p) <= 1.2 * h) return; + if (GetH(p) <= 1.2 * h) return; - GradingBox * box = root; - GradingBox * nbox = root; - GradingBox * ngb; - int childnr; - double x1[3], x2[3]; + GradingBox * box = root; + GradingBox * nbox = root; + GradingBox * ngb; + int childnr; + double x1[3], x2[3]; - while (nbox) - { - box = nbox; - childnr = 0; - if (p.X() > box->xmid[0]) childnr += 1; - if (p.Y() > box->xmid[1]) childnr += 2; - if (p.Z() > box->xmid[2]) childnr += 4; - nbox = box->childs[childnr]; - }; - - - while (2 * box->h2 > h) - { - childnr = 0; - if (p.X() > box->xmid[0]) childnr += 1; - if (p.Y() > box->xmid[1]) childnr += 2; - if (p.Z() > box->xmid[2]) childnr += 4; - - double h2 = box->h2; - if (childnr & 1) - { - x1[0] = box->xmid[0]; - x2[0] = x1[0]+h2; // box->x2[0]; - } - else - { - x2[0] = box->xmid[0]; - x1[0] = x2[0]-h2; // box->x1[0]; - } - - if (childnr & 2) - { - x1[1] = box->xmid[1]; - x2[1] = x1[1]+h2; // box->x2[1]; - } - else - { - x2[1] = box->xmid[1]; - x1[1] = x2[1]-h2; // box->x1[1]; - } - - if (childnr & 4) - { - x1[2] = box->xmid[2]; - x2[2] = x1[2]+h2; // box->x2[2]; - } - else - { - x2[2] = box->xmid[2]; - x1[2] = x2[2]-h2; // box->x1[2]; - } - - ngb = new GradingBox (x1, x2); - box->childs[childnr] = ngb; - ngb->father = box; - - boxes.Append (ngb); - box = box->childs[childnr]; - } - - box->hopt = h; - - - double hbox = 2 * box->h2; // box->x2[0] - box->x1[0]; - double hnp = h + grading * hbox; - - Point3d np; - int i; - for (i = 1; i <= 3; i++) - { - np = p; - np.X(i) = p.X(i) + hbox; - SetH (np, hnp); - - np.X(i) = p.X(i) - hbox; - SetH (np, hnp); - } - /* - Point3d np; - int i1, i2, i3; - for (i1 = -1; i1 <= 1; i1++) - for (i2 = -1; i2 <= 1; i2++) - for (i3 = -1; i3 <= 1; i3++) - { - np.X() = p.X() + hbox * i1; - np.Y() = p.Y() + hbox * i2; - np.Z() = p.Z() + hbox * i3; - - SetH (np, hnp); - } - */ -} - - - -double LocalH :: GetH (const Point3d & x) const -{ - const GradingBox * box = root; - const GradingBox * nbox; - int childnr; - - while (1) - { - childnr = 0; - if (x.X() > box->xmid[0]) childnr += 1; - if (x.Y() > box->xmid[1]) childnr += 2; - if (x.Z() > box->xmid[2]) childnr += 4; - nbox = box->childs[childnr]; - if (nbox) + while (nbox) + { box = nbox; - else - { - // (*testout) << "diam = " << (box->x2[0] - box->x1[0]) - // << " h = " << box->hopt << endl; + childnr = 0; + if (p.X() > box->xmid[0]) childnr += 1; + if (p.Y() > box->xmid[1]) childnr += 2; + if (p.Z() > box->xmid[2]) childnr += 4; + nbox = box->childs[childnr]; + }; + + + while (2 * box->h2 > h) + { + childnr = 0; + if (p.X() > box->xmid[0]) childnr += 1; + if (p.Y() > box->xmid[1]) childnr += 2; + if (p.Z() > box->xmid[2]) childnr += 4; + + double h2 = box->h2; + if (childnr & 1) + { + x1[0] = box->xmid[0]; + x2[0] = x1[0]+h2; // box->x2[0]; + } + else + { + x2[0] = box->xmid[0]; + x1[0] = x2[0]-h2; // box->x1[0]; + } + + if (childnr & 2) + { + x1[1] = box->xmid[1]; + x2[1] = x1[1]+h2; // box->x2[1]; + } + else + { + x2[1] = box->xmid[1]; + x1[1] = x2[1]-h2; // box->x1[1]; + } + + if (childnr & 4) + { + x1[2] = box->xmid[2]; + x2[2] = x1[2]+h2; // box->x2[2]; + } + else + { + x2[2] = box->xmid[2]; + x1[2] = x2[2]-h2; // box->x1[2]; + } + + ngb = new GradingBox (x1, x2); + box->childs[childnr] = ngb; + ngb->father = box; + + boxes.Append (ngb); + box = box->childs[childnr]; + } + + box->hopt = h; + + + double hbox = 2 * box->h2; // box->x2[0] - box->x1[0]; + double hnp = h + grading * hbox; + + Point3d np; + for (int i = 1; i <= 3; i++) + { + np = p; + np.X(i) = p.X(i) + hbox; + SetH (np, hnp); + + np.X(i) = p.X(i) - hbox; + SetH (np, hnp); + } + } + + + + double LocalH :: GetH (const Point3d & x) const + { + const GradingBox * box = root; + + while (1) + { + int childnr = 0; + if (x.X() > box->xmid[0]) childnr += 1; + if (x.Y() > box->xmid[1]) childnr += 2; + if (x.Z() > box->xmid[2]) childnr += 4; + + if (box->childs[childnr]) + box = box->childs[childnr]; + else return box->hopt; - } - } -} + } + } -/// minimal h in box (pmin, pmax) -double LocalH :: GetMinH (const Point3d & pmin, const Point3d & pmax) const -{ - Point3d pmin2, pmax2; - for (int j = 1; j <= 3; j++) - if (pmin.X(j) < pmax.X(j)) - { pmin2.X(j) = pmin.X(j); pmax2.X(j) = pmax.X(j); } - else - { pmin2.X(j) = pmax.X(j); pmax2.X(j) = pmin.X(j); } + /// minimal h in box (pmin, pmax) + double LocalH :: GetMinH (const Point3d & pmin, const Point3d & pmax) const + { + Point3d pmin2, pmax2; + for (int j = 1; j <= 3; j++) + if (pmin.X(j) < pmax.X(j)) + { pmin2.X(j) = pmin.X(j); pmax2.X(j) = pmax.X(j); } + else + { pmin2.X(j) = pmax.X(j); pmax2.X(j) = pmin.X(j); } - return GetMinHRec (pmin2, pmax2, root); -} + return GetMinHRec (pmin2, pmax2, root); + } -double LocalH :: GetMinHRec (const Point3d & pmin, const Point3d & pmax, - const GradingBox * box) const -{ - double h2 = box->h2; - if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 || - pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 || - pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2) - return 1e8; - /* - if (pmax.X() < box->x1[0] || pmin.X() > box->x2[0] || - pmax.Y() < box->x1[1] || pmin.Y() > box->x2[1] || - pmax.Z() < box->x1[2] || pmin.Z() > box->x2[2]) - return 1e8; - */ - + double LocalH :: GetMinHRec (const Point3d & pmin, const Point3d & pmax, + const GradingBox * box) const + { + double h2 = box->h2; + if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 || + pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 || + pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2) + return 1e8; - double hmin = 2 * box->h2; // box->x2[0] - box->x1[0]; - int i; + double hmin = 2 * box->h2; // box->x2[0] - box->x1[0]; - for (i = 0; i <= 7; i++) - { + for (int i = 0; i < 8; i++) if (box->childs[i]) - { - double hi = GetMinHRec (pmin, pmax, box->childs[i]); - if (hi < hmin) - hmin = hi; - } - } + hmin = min2 (hmin, GetMinHRec (pmin, pmax, box->childs[i])); - return hmin; -} + return hmin; + } -void LocalH :: CutBoundaryRec (const Point3d & pmin, const Point3d & pmax, - GradingBox * box) -{ - double h2 = box->h2; - if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 || - pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 || - pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2) - return; - /* - if (pmax.X() < box->x1[0] || pmin.X() > box->x2[0] || - pmax.Y() < box->x1[1] || pmin.Y() > box->x2[1] || - pmax.Z() < box->x1[2] || pmin.Z() > box->x2[2]) - return; - */ + void LocalH :: CutBoundaryRec (const Point3d & pmin, const Point3d & pmax, + GradingBox * box) + { + double h2 = box->h2; + if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 || + pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 || + pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2) + return; - box->flags.cutboundary = 1; - for (int i = 0; i < 8; i++) - if (box->childs[i]) - CutBoundaryRec (pmin, pmax, box->childs[i]); -} + + box->flags.cutboundary = 1; + for (int i = 0; i < 8; i++) + if (box->childs[i]) + CutBoundaryRec (pmin, pmax, box->childs[i]); + } -void LocalH :: FindInnerBoxes ( // int (*sameside)(const Point3d & p1, const Point3d & p2), - AdFront3 * adfront, - int (*testinner)(const Point3d & p1)) -{ - int i; + void LocalH :: FindInnerBoxes (AdFront3 * adfront, + int (*testinner)(const Point3d & p1)) + { + int nf = adfront->GetNF(); - int nf = adfront->GetNF(); + for (int i = 0; i < boxes.Size(); i++) + boxes[i] -> flags.isinner = 0; - for (i = 0; i < boxes.Size(); i++) - boxes[i] -> flags.isinner = 0; + root->flags.isinner = 0; - root->flags.isinner = 0; - - Point3d rpmid(root->xmid[0], root->xmid[1], root->xmid[2]); - Vec3d rv(root->h2, root->h2, root->h2); - Point3d rx2 = rpmid + rv; - Point3d rx1 = rpmid - rv; + Point3d rpmid(root->xmid[0], root->xmid[1], root->xmid[2]); + Vec3d rv(root->h2, root->h2, root->h2); + Point3d rx2 = rpmid + rv; + Point3d rx1 = rpmid - rv; - root->flags.pinner = !adfront->SameSide (rpmid, rx2); + root->flags.pinner = !adfront->SameSide (rpmid, rx2); - if (testinner) - (*testout) << "inner = " << root->flags.pinner << " =?= " - << testinner(Point3d(root->xmid[0], root->xmid[1], root->xmid[2])) << endl; + if (testinner) + (*testout) << "inner = " << root->flags.pinner << " =?= " + << testinner(Point3d(root->xmid[0], root->xmid[1], root->xmid[2])) << endl; - Array faceinds(nf); - Array faceboxes(nf); + Array faceinds(nf); + Array faceboxes(nf); - for (i = 1; i <= nf; i++) - { - faceinds.Elem(i) = i; - adfront->GetFaceBoundingBox(i, faceboxes.Elem(i)); - } + for (int i = 1; i <= nf; i++) + { + faceinds.Elem(i) = i; + adfront->GetFaceBoundingBox(i, faceboxes.Elem(i)); + } - for (i = 0; i < 8; i++) - FindInnerBoxesRec2 (root->childs[i], adfront, faceboxes, faceinds, nf); -} + for (int i = 0; i < 8; i++) + FindInnerBoxesRec2 (root->childs[i], adfront, faceboxes, faceinds, nf); + } -void LocalH :: -FindInnerBoxesRec2 (GradingBox * box, - class AdFront3 * adfront, - Array & faceboxes, - Array & faceinds, int nfinbox) -{ - if (!box) return; + void LocalH :: + FindInnerBoxesRec2 (GradingBox * box, + class AdFront3 * adfront, + Array & faceboxes, + Array & faceinds, int nfinbox) + { + if (!box) return; - int i, j; + GradingBox * father = box -> father; - GradingBox * father = box -> father; + Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]); + Vec3d v(box->h2, box->h2, box->h2); + Box3d boxc(c-v, c+v); + + Point3d fc(father->xmid[0], father->xmid[1], father->xmid[2]); + Vec3d fv(father->h2, father->h2, father->h2); + Box3d fboxc(fc-fv, fc+fv); + + Box3d boxcfc(c,fc); + + ArrayMem faceused; + ArrayMem faceused2; + ArrayMem facenotused; + + /* + faceused.SetSize(0); + facenotused.SetSize(0); + faceused2.SetSize(0); + */ + + for (int j = 1; j <= nfinbox; j++) + { + // adfront->GetFaceBoundingBox (faceinds.Get(j), facebox); + const Box3d & facebox = faceboxes.Get(faceinds.Get(j)); - Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]); - Vec3d v(box->h2, box->h2, box->h2); - Box3d boxc(c-v, c+v); + if (boxc.Intersect (facebox)) + faceused.Append(faceinds.Get(j)); + else + facenotused.Append(faceinds.Get(j)); - Point3d fc(father->xmid[0], father->xmid[1], father->xmid[2]); - Vec3d fv(father->h2, father->h2, father->h2); - Box3d fboxc(fc-fv, fc+fv); - - Box3d boxcfc(c,fc); - - - static Array faceused; - static Array faceused2; - static Array facenotused; - - faceused.SetSize(0); - facenotused.SetSize(0); - faceused2.SetSize(0); - - for (j = 1; j <= nfinbox; j++) - { - // adfront->GetFaceBoundingBox (faceinds.Get(j), facebox); - const Box3d & facebox = faceboxes.Get(faceinds.Get(j)); + if (boxcfc.Intersect (facebox)) + faceused2.Append (faceinds.Get(j)); + } - if (boxc.Intersect (facebox)) - faceused.Append(faceinds.Get(j)); - else - facenotused.Append(faceinds.Get(j)); - - if (boxcfc.Intersect (facebox)) - faceused2.Append (faceinds.Get(j)); - } - - for (j = 1; j <= faceused.Size(); j++) - faceinds.Elem(j) = faceused.Get(j); - for (j = 1; j <= facenotused.Size(); j++) - faceinds.Elem(j+faceused.Size()) = facenotused.Get(j); + for (int j = 1; j <= faceused.Size(); j++) + faceinds.Elem(j) = faceused.Get(j); + for (int j = 1; j <= facenotused.Size(); j++) + faceinds.Elem(j+faceused.Size()) = facenotused.Get(j); - if (!father->flags.cutboundary) - { - box->flags.isinner = father->flags.isinner; - box->flags.pinner = father->flags.pinner; - } - else - { - Point3d cf(father->xmid[0], father->xmid[1], father->xmid[2]); + if (!father->flags.cutboundary) + { + box->flags.isinner = father->flags.isinner; + box->flags.pinner = father->flags.pinner; + } + else + { + Point3d cf(father->xmid[0], father->xmid[1], father->xmid[2]); - if (father->flags.isinner) - box->flags.pinner = 1; - else - { - if (adfront->SameSide (c, cf, &faceused2)) - box->flags.pinner = father->flags.pinner; - else - box->flags.pinner = 1 - father->flags.pinner; - } + if (father->flags.isinner) + box->flags.pinner = 1; + else + { + if (adfront->SameSide (c, cf, &faceused2)) + box->flags.pinner = father->flags.pinner; + else + box->flags.pinner = 1 - father->flags.pinner; + } - if (box->flags.cutboundary) - box->flags.isinner = 0; - else - box->flags.isinner = box->flags.pinner; - } + if (box->flags.cutboundary) + box->flags.isinner = 0; + else + box->flags.isinner = box->flags.pinner; + } - int nf = faceused.Size(); - for (i = 0; i < 8; i++) - FindInnerBoxesRec2 (box->childs[i], adfront, faceboxes, faceinds, nf); -} + // cout << "faceused: " << faceused.Size() << ", " << faceused2.Size() << ", " << facenotused.Size() << endl; + + int nf = faceused.Size(); + for (int i = 0; i < 8; i++) + FindInnerBoxesRec2 (box->childs[i], adfront, faceboxes, faceinds, nf); + } + + + + void LocalH :: FindInnerBoxesRec ( int (*inner)(const Point3d & p), + GradingBox * box) + { + if (box->flags.cutboundary) + { + for (int i = 0; i < 8; i++) + if (box->childs[i]) + FindInnerBoxesRec (inner, box->childs[i]); + } + else + { + if (inner (box->PMid())) + SetInnerBoxesRec (box); + } + } @@ -464,217 +457,252 @@ FindInnerBoxesRec2 (GradingBox * box, -/* -void LocalH :: FindInnerBoxes ( // int (*sameside)(const Point3d & p1, const Point3d & p2), - AdFront3 * adfront, - int (*testinner)(const Point3d & p1)) -{ - int i; - for (i = 1; i <= boxes.Size(); i++) - boxes.Elem(i)->flags.isinner = 0; - root->flags.isinner = 0; - Point3d rpmid(root->xmid[0], root->xmid[1], root->xmid[2]); - Point3d rx2 = rpmid + Vec3d (root->h2, root->h2, root->h2); - root->flags.pinner = !adfront->SameSide (rpmid, rx2); + + + + + void LocalH :: FindInnerBoxes (AdFront2 * adfront, + int (*testinner)(const Point<2> & p1)) + { + int nf = adfront->GetNFL(); + + for (int i = 0; i < boxes.Size(); i++) + boxes[i] -> flags.isinner = 0; + + root->flags.isinner = 0; + + Point<2> rpmid(root->xmid[0], root->xmid[1], root->xmid[2]); + Vec<2> rv(root->h2, root->h2); + Point<2> rx2 = rpmid + rv; + Point<2> rx1 = rpmid - rv; + + + root->flags.pinner = !adfront->SameSide (rpmid, rx2); - if (testinner) - (*testout) << "inner = " << root->flags.pinner << " =?= " - << testinner(Point3d(root->xmid[0], root->xmid[1], root->xmid[2])) << endl; + if (testinner) + (*testout) << "inner = " << root->flags.pinner << " =?= " + << testinner(rpmid) << endl; + Array faceinds(nf); + Array > faceboxes(nf); - for (i = 2; i <= boxes.Size(); i++) - { - GradingBox * box = boxes.Elem(i); - GradingBox * father = box -> father; + for (int i = 0; i < nf; i++) + { + faceinds[i] = i; + // adfront->GetFaceBoundingBox(i, faceboxes.Elem(i)); - Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]); - Vec3d v(box->h2, box->h2, box->h2); - Point3d x1 = c-v; - Point3d x2 = c+v; + const FrontLine & line = adfront->GetLine(i); + faceboxes[i].Set (adfront->GetPoint (line.L().I1())); + faceboxes[i].Add (adfront->GetPoint (line.L().I2())); + + } + + for (int i = 0; i < 8; i++) + FindInnerBoxesRec2 (root->childs[i], adfront, faceboxes, faceinds, nf); + } - if (!father->flags.cutboundary) - { - box->flags.isinner = father->flags.isinner; - box->flags.pinner = father->flags.pinner; - } - else - { - Point3d cf(father->xmid[0], father->xmid[1], father->xmid[2]); + void LocalH :: + FindInnerBoxesRec2 (GradingBox * box, + class AdFront2 * adfront, + Array > & faceboxes, + Array & faceinds, int nfinbox) + { + if (!box) return; + + GradingBox * father = box -> father; + + Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]); + Vec3d v(box->h2, box->h2, box->h2); + Box3d boxc(c-v, c+v); - if (father->flags.isinner) - box->flags.pinner = 1; - else - { - if (adfront->SameSide (c, cf)) - box->flags.pinner = father->flags.pinner; - else - box->flags.pinner = 1 - father->flags.pinner; - } + Point3d fc(father->xmid[0], father->xmid[1], father->xmid[2]); + Vec3d fv(father->h2, father->h2, father->h2); + Box3d fboxc(fc-fv, fc+fv); - if (box->flags.cutboundary) - box->flags.isinner = 0; - else - box->flags.isinner = box->flags.pinner; - } - } - // FindInnerBoxesRec (inner, root); -} -*/ + Box3d boxcfc(c,fc); + + ArrayMem faceused; + ArrayMem faceused2; + ArrayMem facenotused; -void LocalH :: FindInnerBoxesRec ( int (*inner)(const Point3d & p), - GradingBox * box) -{ - int i; - if (box->flags.cutboundary) - { - for (i = 0; i < 8; i++) - if (box->childs[i]) - FindInnerBoxesRec (inner, box->childs[i]); - } - else - { - if (inner (Point3d (box->xmid[0], box->xmid[1], box->xmid[2]))) - SetInnerBoxesRec (box); - } -} + for (int j = 1; j <= nfinbox; j++) + { + // adfront->GetFaceBoundingBox (faceinds.Get(j), facebox); + const Box3d & facebox = faceboxes.Get(faceinds.Get(j)); + + if (boxc.Intersect (facebox)) + faceused.Append(faceinds.Get(j)); + else + facenotused.Append(faceinds.Get(j)); - -void LocalH :: SetInnerBoxesRec (GradingBox * box) -{ - box->flags.isinner = 1; - for (int i = 0; i < 8; i++) - if (box->childs[i]) - ClearFlagsRec (box->childs[i]); -} - -void LocalH :: ClearFlagsRec (GradingBox * box) -{ - box->flags.cutboundary = 0; - box->flags.isinner = 0; - for (int i = 0; i < 8; i++) - if (box->childs[i]) - ClearFlagsRec (box->childs[i]); -} - - -void LocalH :: WidenRefinement () -{ - int nb = boxes.Size(); - int i; - // (*testout) << "old boxes: " << nb << endl; - for (i = 1; i <= nb; i++) - { - GradingBox * box = boxes.Get(i); - // double h = box->x2[0] - box->x1[0]; - double h = box->hopt; - Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]); - // (*testout) << " i = " << i - // << " c = " << c << " h = " << h << endl; - - for (int i1 = -1; i1 <= 1; i1++) - for (int i2 = -1; i2 <= 1; i2++) - for (int i3 = -1; i3 <= 1; i3++) - SetH (Point3d (c.X() + i1 * h, - c.Y() + i2 * h, - c.Z() + i3 * h), 1.001 * h); - } -} - -void LocalH :: GetInnerPoints (Array & points) -{ - int i, nb = boxes.Size(); - - for (i = 1; i <= nb; i++) - { - GradingBox * box = boxes.Get(i); - /* - if (box->flags.pinner) - points.Append (box->randomip); - */ - // if (box->flags.pinner) - if (box->flags.isinner) - { - Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]); - points.Append (c); - /* - cout << "add point " << c << "; h = " << box->hopt - << "; max-min = " << (box->x2[0]-box->x1[0]) << endl; - */ - } - } -} - - - -void LocalH :: GetOuterPoints (Array & points) -{ - int i, nb = boxes.Size(); - - for (i = 1; i <= nb; i++) - { - GradingBox * box = boxes.Get(i); - if (!box->flags.isinner && - !box->flags.cutboundary) - { - Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]); - points.Append (c); - } - } -} - - - -void LocalH :: Convexify () -{ - ConvexifyRec (root); -} - -void LocalH :: ConvexifyRec (GradingBox * box) -{ - Point3d center(box->xmid[0], box->xmid[1], box->xmid[2]); - Point3d hp; - - double size = 2 * box->h2; // box->x2[0] - box->x1[0]; - double dx = 0.6 * size; - - double maxh = box->hopt; - int i; + if (boxcfc.Intersect (facebox)) + faceused2.Append (faceinds.Get(j)); + } + + for (int j = 1; j <= faceused.Size(); j++) + faceinds.Elem(j) = faceused.Get(j); + for (int j = 1; j <= facenotused.Size(); j++) + faceinds.Elem(j+faceused.Size()) = facenotused.Get(j); - - for (i = 1; i <= 6; i++) - { - hp = center; - switch (i) - { - case 1: hp.X() += dx; break; - case 2: hp.X() -= dx; break; - case 3: hp.Y() += dx; break; - case 4: hp.Y() -= dx; break; - case 5: hp.Z() += dx; break; - case 6: hp.Z() -= dx; break; - } + if (!father->flags.cutboundary) + { + box->flags.isinner = father->flags.isinner; + box->flags.pinner = father->flags.pinner; + } + else + { + Point3d cf(father->xmid[0], father->xmid[1], father->xmid[2]); - double hh = GetH (hp); - if (hh > maxh) maxh = hh; - } + if (father->flags.isinner) + box->flags.pinner = 1; + else + { + Point<2> c2d (c.X(), c.Y()); + Point<2> cf2d (cf.X(), cf.Y()); + if (adfront->SameSide (c2d, cf2d, &faceused2)) + box->flags.pinner = father->flags.pinner; + else + box->flags.pinner = 1 - father->flags.pinner; + } + + if (box->flags.cutboundary) + box->flags.isinner = 0; + else + box->flags.isinner = box->flags.pinner; + } - if (maxh < 0.95 * box->hopt) - SetH (center, maxh); + // cout << "faceused: " << faceused.Size() << ", " << faceused2.Size() << ", " << facenotused.Size() << endl; - for (i = 0; i < 8; i++) - if (box->childs[i]) - ConvexifyRec (box->childs[i]); -} + int nf = faceused.Size(); + for (int i = 0; i < 8; i++) + FindInnerBoxesRec2 (box->childs[i], adfront, faceboxes, faceinds, nf); + } -void LocalH :: PrintMemInfo (ostream & ost) const -{ - ost << "LocalH: " << boxes.Size() << " boxes of " << sizeof(GradingBox) - << " bytes = " << boxes.Size()*sizeof(GradingBox) << " bytes" << endl; -} + + + void LocalH :: FindInnerBoxesRec ( int (*inner)(const Point<2> & p), + GradingBox * box) + { + if (box->flags.cutboundary) + { + for (int i = 0; i < 8; i++) + if (box->childs[i]) + FindInnerBoxesRec (inner, box->childs[i]); + } + else + { + Point<2> p2d(box->PMid()(0), box->PMid()(1)); + if (inner (p2d)) + SetInnerBoxesRec (box); + } + } + + + + + + + + + + + + + + + + + + void LocalH :: SetInnerBoxesRec (GradingBox * box) + { + box->flags.isinner = 1; + for (int i = 0; i < 8; i++) + if (box->childs[i]) + ClearFlagsRec (box->childs[i]); + } + + void LocalH :: ClearFlagsRec (GradingBox * box) + { + box->flags.cutboundary = 0; + box->flags.isinner = 0; + for (int i = 0; i < 8; i++) + if (box->childs[i]) + ClearFlagsRec (box->childs[i]); + } + + + void LocalH :: WidenRefinement () + { + for (int i = 0; i < boxes.Size(); i++) + { + double h = boxes[i]->hopt; + Point3d c = boxes[i]->PMid(); + + for (int i1 = -1; i1 <= 1; i1++) + for (int i2 = -1; i2 <= 1; i2++) + for (int i3 = -1; i3 <= 1; i3++) + SetH (Point3d (c.X() + i1 * h, + c.Y() + i2 * h, + c.Z() + i3 * h), 1.001 * h); + } + } + + void LocalH :: GetInnerPoints (Array > & points) + { + for (int i = 0; i < boxes.Size(); i++) + if (boxes[i] -> flags.isinner) + points.Append ( boxes[i] -> PMid() ); + } + + + void LocalH :: GetOuterPoints (Array > & points) + { + for (int i = 0; i < boxes.Size(); i++) + if (!boxes[i]->flags.isinner && !boxes[i]->flags.cutboundary) + points.Append ( boxes[i] -> PMid()); + } + + + void LocalH :: Convexify () + { + ConvexifyRec (root); + } + + void LocalH :: ConvexifyRec (GradingBox * box) + { + Point<3> center = box -> PMid(); + + double size = 2 * box->h2; // box->x2[0] - box->x1[0]; + double dx = 0.6 * size; + + double maxh = box->hopt; + + for (int i = 0; i < 3; i++) + { + Point<3> hp = center; + hp(i) += dx; + maxh = max2 (maxh, GetH(hp)); + hp(i) = center(i)-dx; + maxh = max2 (maxh, GetH(hp)); + } + + if (maxh < 0.95 * box->hopt) + SetH (center, maxh); + + for (int i = 0; i < 8; i++) + if (box->childs[i]) + ConvexifyRec (box->childs[i]); + } + + void LocalH :: PrintMemInfo (ostream & ost) const + { + ost << "LocalH: " << boxes.Size() << " boxes of " << sizeof(GradingBox) + << " bytes = " << boxes.Size()*sizeof(GradingBox) << " bytes" << endl; + } } diff --git a/libsrc/meshing/localh.hpp b/libsrc/meshing/localh.hpp index 3c1aa0d9..583b6e0a 100644 --- a/libsrc/meshing/localh.hpp +++ b/libsrc/meshing/localh.hpp @@ -13,12 +13,6 @@ /// box for grading class GradingBox { - /* - /// xmin - float x1[3]; - /// xmax - float x2[3]; - */ /// xmid float xmid[3]; /// half edgelength @@ -30,6 +24,8 @@ class GradingBox /// double hopt; /// +public: + struct { unsigned int cutboundary:1; @@ -37,14 +33,17 @@ class GradingBox unsigned int oldcell:1; unsigned int pinner:1; } flags; -public: + /// GradingBox (const double * ax1, const double * ax2); /// void DeleteChilds(); /// - friend class LocalH; + Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); } + double H2() const { return h2; } + + friend class LocalH; static BlockAllocator ball; void * operator new(size_t); @@ -53,6 +52,7 @@ public: + /** Control of 3D mesh grading */ @@ -70,6 +70,8 @@ public: /// LocalH (const Point3d & pmin, const Point3d & pmax, double grading); /// + LocalH (const Box<3> & box, double grading); + /// ~LocalH(); /// void Delete(); @@ -83,14 +85,19 @@ public: double GetMinH (const Point3d & pmin, const Point3d & pmax) const; /// mark boxes intersecting with boundary-box - void CutBoundary (const Point3d & pmin, const Point3d & pmax) - { CutBoundaryRec (pmin, pmax, root); } - + // void CutBoundary (const Point3d & pmin, const Point3d & pmax) + // { CutBoundaryRec (pmin, pmax, root); } + void CutBoundary (const Box<3> & box) + { CutBoundaryRec (box.PMin(), box.PMax(), root); } + /// find inner boxes - void FindInnerBoxes ( // int (*sameside)(const Point3d & p1, const Point3d & p2), - class AdFront3 * adfront, + void FindInnerBoxes (class AdFront3 * adfront, int (*testinner)(const Point3d & p1)); + void FindInnerBoxes (class AdFront2 * adfront, + int (*testinner)(const Point<2> & p1)); + + /// clears all flags void ClearFlags () { ClearFlagsRec(root); } @@ -99,10 +106,10 @@ public: void WidenRefinement (); /// get points in inner elements - void GetInnerPoints (Array & points); + void GetInnerPoints (Array > & points); /// get points in outer closure - void GetOuterPoints (Array & points); + void GetOuterPoints (Array > & points); /// void Convexify (); @@ -131,6 +138,18 @@ private: Array & finds, int nfinbox); + + void FindInnerBoxesRec ( int (*inner)(const Point<2> & p), + GradingBox * box); + + /// + void FindInnerBoxesRec2 (GradingBox * box, + class AdFront2 * adfront, + Array > & faceboxes, + Array & finds, int nfinbox); + + + /// void SetInnerBoxesRec (GradingBox * box); @@ -139,7 +158,26 @@ private: /// void ConvexifyRec (GradingBox * box); + + friend ostream & operator<< (ostream & ost, const LocalH & loch); }; + + +inline ostream & operator<< (ostream & ost, const GradingBox & box) +{ + ost << "gradbox, pmid = " << box.PMid() << ", h2 = " << box.H2() + << " cutbound = " << box.flags.cutboundary << " isinner = " << box.flags.isinner + << endl; + return ost; +} + +inline ostream & operator<< (ostream & ost, const LocalH & loch) +{ + for (int i = 0; i < loch.boxes.Size(); i++) + ost << "box[" << i << "] = " << *(loch.boxes[i]); + return ost; +} + #endif diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index feeca6ae..db36fe1a 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -4365,7 +4365,7 @@ namespace netgen - if(el.GetType() == TET || el.GetType() == TET10) + if (el.GetType() == TET) { retval = (lam(0) > -eps && lam(1) > -eps && diff --git a/libsrc/meshing/meshing2.hpp b/libsrc/meshing/meshing2.hpp index ce80c8a7..278b4e64 100644 --- a/libsrc/meshing/meshing2.hpp +++ b/libsrc/meshing/meshing2.hpp @@ -54,6 +54,10 @@ public: /// MESHING2_RESULT GenerateMesh (Mesh & mesh, double gh, int facenr); + void Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp); + void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp); + + /// void AddPoint (const Point3d & p, PointIndex globind, MultiPointGeomInfo * mgi = NULL, bool pointonsurface = true); diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp index ff683d63..3421c319 100644 --- a/libsrc/meshing/meshing3.cpp +++ b/libsrc/meshing/meshing3.cpp @@ -1082,7 +1082,7 @@ void Meshing3 :: BlockFill (Mesh & mesh, double gh) } - +/* static const AdFront3 * locadfront; static int TestInner (const Point3d & p) { @@ -1092,62 +1092,45 @@ static int TestSameSide (const Point3d & p1, const Point3d & p2) { return locadfront->SameSide (p1, p2); } - +*/ void Meshing3 :: BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp) { - int i, j; - double filldist = mp.filldist; (*testout) << "blockfill local h" << endl; (*testout) << "rel filldist = " << filldist << endl; PrintMessage (3, "blockfill local h"); - /* - (*mycout) << "boxes: " << mesh.LocalHFunction().GetNBoxes() << endl - << "filldist = " << filldist << endl; - */ - Array npoints; + + Array > npoints; adfront -> CreateTrees(); - Point3d mpmin, mpmax; - // mesh.GetBox (mpmin, mpmax); - bool firstp = 1; - + Box<3> bbox ( Box<3>::EMPTY_BOX ); double maxh = 0; - for (i = 1; i <= adfront->GetNF(); i++) + + for (int i = 1; i <= adfront->GetNF(); i++) { const MiniElement2d & el = adfront->GetFace(i); - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { const Point3d & p1 = adfront->GetPoint (el.PNumMod(j)); const Point3d & p2 = adfront->GetPoint (el.PNumMod(j+1)); - double hi = Dist (p1, p2); - if (hi > maxh) - { - maxh = hi; - //(*testout) << "reducing maxh to " << maxh << " because of " << p1 << " and " << p2 << endl; - } - if (firstp) - { - mpmin = p1; - mpmax = p1; - firstp = 0; - } - else - { - mpmin.SetToMin (p1); - mpmax.SetToMax (p1); - } + double hi = Dist (p1, p2); + if (hi > maxh) maxh = hi; + + bbox.Add (p1); } } + + Point3d mpmin = bbox.PMin(); + Point3d mpmax = bbox.PMax(); Point3d mpc = Center (mpmin, mpmax); double d = max3(mpmax.X()-mpmin.X(), mpmax.Y()-mpmin.Y(), @@ -1158,53 +1141,41 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, LocalH loch2 (mpmin, mpmax, 1); + if (mp.maxh < maxh) maxh = mp.maxh; - if (mp.maxh < maxh) - { - maxh = mp.maxh; - //(*testout) << "reducing maxh to " << maxh << " because of mp.maxh" << endl; - } - - int changed; + bool changed; do { mesh.LocalHFunction().ClearFlags(); - for (i = 1; i <= adfront->GetNF(); i++) + for (int i = 1; i <= adfront->GetNF(); i++) { const MiniElement2d & el = adfront->GetFace(i); - Point3d pmin = adfront->GetPoint (el.PNum(1)); - Point3d pmax = pmin; - - for (j = 2; j <= 3; j++) - { - const Point3d & p = adfront->GetPoint (el.PNum(j)); - pmin.SetToMin (p); - pmax.SetToMax (p); - } + Box<3> bbox (adfront->GetPoint (el[0])); + bbox.Add (adfront->GetPoint (el[1])); + bbox.Add (adfront->GetPoint (el[2])); - double filld = filldist * Dist (pmin, pmax); - - pmin = pmin - Vec3d (filld, filld, filld); - pmax = pmax + Vec3d (filld, filld, filld); - // (*testout) << "cut : " << pmin << " - " << pmax << endl; - mesh.LocalHFunction().CutBoundary (pmin, pmax); + + double filld = filldist * bbox.Diam(); + bbox.Increase (filld); + + mesh.LocalHFunction().CutBoundary (bbox); // .PMin(), bbox.PMax()); } - locadfront = adfront; + // locadfront = adfront; mesh.LocalHFunction().FindInnerBoxes (adfront, NULL); npoints.SetSize(0); mesh.LocalHFunction().GetInnerPoints (npoints); - changed = 0; - for (i = 1; i <= npoints.Size(); i++) + changed = false; + for (int i = 1; i <= npoints.Size(); i++) { if (mesh.LocalHFunction().GetH(npoints.Get(i)) > 1.5 * maxh) { mesh.LocalHFunction().SetH (npoints.Get(i), maxh); - changed = 1; + changed = true; } } } @@ -1212,7 +1183,7 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, if (debugparam.slowchecks) (*testout) << "Blockfill with points: " << endl; - for (i = 1; i <= npoints.Size(); i++) + for (int i = 1; i <= npoints.Size(); i++) { if (meshbox.IsIn (npoints.Get(i))) { @@ -1238,13 +1209,13 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, loch2.ClearFlags(); - for (i = 1; i <= adfront->GetNF(); i++) + for (int i = 1; i <= adfront->GetNF(); i++) { const MiniElement2d & el = adfront->GetFace(i); Point3d pmin = adfront->GetPoint (el.PNum(1)); Point3d pmax = pmin; - for (j = 2; j <= 3; j++) + for (int j = 2; j <= 3; j++) { const Point3d & p = adfront->GetPoint (el.PNum(j)); pmin.SetToMin (p); @@ -1254,13 +1225,13 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, loch2.SetH (Center (pmin, pmax), Dist (pmin, pmax)); } - for (i = 1; i <= adfront->GetNF(); i++) + for (int i = 1; i <= adfront->GetNF(); i++) { const MiniElement2d & el = adfront->GetFace(i); Point3d pmin = adfront->GetPoint (el.PNum(1)); Point3d pmax = pmin; - for (j = 2; j <= 3; j++) + for (int j = 2; j <= 3; j++) { const Point3d & p = adfront->GetPoint (el.PNum(j)); pmin.SetToMin (p); @@ -1270,16 +1241,17 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, double filld = filldist * Dist (pmin, pmax); pmin = pmin - Vec3d (filld, filld, filld); pmax = pmax + Vec3d (filld, filld, filld); - loch2.CutBoundary (pmin, pmax); + // loch2.CutBoundary (pmin, pmax); + loch2.CutBoundary (Box<3> (pmin, pmax)); // pmin, pmax); } - locadfront = adfront; + // locadfront = adfront; loch2.FindInnerBoxes (adfront, NULL); npoints.SetSize(0); loch2.GetOuterPoints (npoints); - for (i = 1; i <= npoints.Size(); i++) + for (int i = 1; i <= npoints.Size(); i++) { if (meshbox.IsIn (npoints.Get(i))) { diff --git a/ng/variables.tcl b/ng/variables.tcl index 29028e95..c9ee8c0b 100644 --- a/ng/variables.tcl +++ b/ng/variables.tcl @@ -496,6 +496,8 @@ proc saveoptions { } { puts $datei "stloptions.recalchopt ${stloptions.recalchopt}" puts $datei "visoptions.subdivisions ${visoptions.subdivisions}" + puts $datei "visoptions.autoredraw ${visoptions.autoredraw}" + puts $datei "visoptions.autoredrawtime ${visoptions.autoredrawtime}" # trafo options diff --git a/nglib/ng_stl.cpp b/nglib/ng_stl.cpp index 7fff7242..af8a2ed8 100644 --- a/nglib/ng_stl.cpp +++ b/nglib/ng_stl.cpp @@ -66,7 +66,7 @@ int main (int argc, char ** argv) Ng_Meshing_Parameters mp; mp.maxh = 1.0e+6; mp.fineness = 0.4; - mp.secondorder = 0; + mp.second_order = 0; cout << "Initialise the STL Geometry structure...." << endl; ng_res = Ng_STL_InitSTLGeometry(stl_geom); diff --git a/nglib/ng_vol.cpp b/nglib/ng_vol.cpp index d67b081c..8ccb3f81 100644 --- a/nglib/ng_vol.cpp +++ b/nglib/ng_vol.cpp @@ -63,7 +63,7 @@ int main (int argc, char ** argv) Ng_Meshing_Parameters mp; mp.maxh = 1e6; mp.fineness = 1; - mp.secondorder = 0; + mp.second_order = 0; cout << "start meshing" << endl; Ng_GenerateVolumeMesh (mesh, &mp);