From 15a1e07092ea5383b76c139a08bd7129d1775e62 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 2 Apr 2013 20:29:05 +0000 Subject: [PATCH] thread-save build-searchtree --- libsrc/meshing/meshclass.cpp | 238 ++++++++++++++++------------------- 1 file changed, 107 insertions(+), 131 deletions(-) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index e6acd79c..ca496265 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -150,17 +150,20 @@ namespace netgen PointIndex Mesh :: AddPoint (const Point3d & p, int layer) { + return AddPoint (p, layer, INNERPOINT); + /* NgLock lock(mutex); lock.Lock(); timestamp = NextTimeStamp(); - PointIndex pi = points.Size() + PointIndex::BASE; + PointIndex pi = points.End(); points.Append ( MeshPoint (p, layer, INNERPOINT) ); lock.UnLock(); return pi; + */ } PointIndex Mesh :: AddPoint (const Point3d & p, int layer, POINTTYPE type) @@ -170,7 +173,7 @@ namespace netgen timestamp = NextTimeStamp(); - PointIndex pi = points.Size() + PointIndex::BASE; + PointIndex pi = points.End(); points.Append ( MeshPoint (p, layer, type) ); lock.UnLock(); @@ -676,13 +679,13 @@ namespace netgen */ int cnt_sing = 0; - for (PointIndex pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) + for (PointIndex pi = points.Begin(); pi < points.End(); pi++) if ((*this)[pi].Singularity()>=1.) cnt_sing++; if (cnt_sing) { outfile << "singular_points" << endl << cnt_sing << endl; - for (PointIndex pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) + for (PointIndex pi = points.Begin(); pi < points.End(); pi++) if ((*this)[pi].Singularity()>=1.) outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl; } @@ -1509,9 +1512,6 @@ namespace netgen void Mesh :: CalcSurfacesOfNode () { - int i, j, k; - SurfaceElementIndex sei; - surfacesonnode.SetSize (GetNP()); delete boundaryedges; @@ -1528,18 +1528,18 @@ namespace netgen surfelementht = new INDEX_3_CLOSED_HASHTABLE (3*GetNSE() + 1); segmentht = new INDEX_2_CLOSED_HASHTABLE (3*GetNSeg() + 1); - for (sei = 0; sei < GetNSE(); sei++) + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { const Element2d & sel = surfelements[sei]; if (sel.IsDeleted()) continue; int si = sel.GetIndex(); - for (j = 0; j < sel.GetNP(); j++) + for (int j = 0; j < sel.GetNP(); j++) { PointIndex pi = sel[j]; bool found = 0; - for (k = 0; k < surfacesonnode[pi].Size(); k++) + for (int k = 0; k < surfacesonnode[pi].Size(); k++) if (surfacesonnode[pi][k] == si) { found = 1; @@ -1548,7 +1548,6 @@ namespace netgen if (!found) surfacesonnode.Add (pi, si); - } } /* @@ -1567,7 +1566,7 @@ namespace netgen surfelementht -> AllocateElements(); */ - for (sei = 0; sei < GetNSE(); sei++) + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { const Element2d & sel = surfelements[sei]; if (sel.IsDeleted()) continue; @@ -1584,17 +1583,16 @@ namespace netgen if (dimension == 3) { - for (PointIndex pi = PointIndex::BASE; - pi < np+PointIndex::BASE; pi++) + for (PointIndex pi = points.Begin(); pi < points.End(); pi++) points[pi].SetType (INNERPOINT); if (GetNFD() == 0) { - for (sei = 0; sei < GetNSE(); sei++) + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { const Element2d & sel = surfelements[sei]; if (sel.IsDeleted()) continue; - for (j = 0; j < sel.GetNP(); j++) + for (int j = 0; j < sel.GetNP(); j++) { PointIndex pi = SurfaceElement(sei)[j]; points[pi].SetType(FIXEDPOINT); @@ -1603,11 +1601,11 @@ namespace netgen } else { - for (sei = 0; sei < GetNSE(); sei++) + for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { const Element2d & sel = surfelements[sei]; if (sel.IsDeleted()) continue; - for (j = 0; j < sel.GetNP(); j++) + for (int j = 0; j < sel.GetNP(); j++) { PointIndex pi = sel[j]; int ns = surfacesonnode[pi].Size(); @@ -1621,10 +1619,10 @@ namespace netgen } } - for (i = 0; i < segments.Size(); i++) + for (int i = 0; i < segments.Size(); i++) { const Segment & seg = segments[i]; - for (j = 1; j <= 2; j++) + for (int j = 1; j <= 2; j++) { PointIndex hi = (j == 1) ? seg[0] : seg[1]; @@ -1635,7 +1633,7 @@ namespace netgen } - for (i = 0; i < lockedpoints.Size(); i++) + for (int i = 0; i < lockedpoints.Size(); i++) points[lockedpoints[i]].SetType(FIXEDPOINT); } @@ -1660,7 +1658,7 @@ namespace netgen // eltyps.SetSize (GetNE()); // eltyps = FREEELEMENT; - for (i = 0; i < GetNSeg(); i++) + for (int i = 0; i < GetNSeg(); i++) { const Segment & seg = segments[i]; INDEX_2 i2(seg[0], seg[1]); @@ -1820,7 +1818,7 @@ namespace netgen INDEX_3_CLOSED_HASHTABLE faceht(100); openelements.SetSize(0); - for (PointIndex pi = PointIndex::BASE; pi < np+PointIndex::BASE; pi++) + for (PointIndex pi = points.Begin(); pi < points.End(); pi++) if (selsonpoint[pi].Size()+elsonpoint[pi].Size()) { faceht.SetSize (2 * selsonpoint[pi].Size() + 4 * elsonpoint[pi].Size()); @@ -1974,8 +1972,8 @@ namespace netgen for (int j = 1; j <= 3; j++) { - int pi = sel.PNum(j); - if (pi < points.Size()+PointIndex::BASE) + PointIndex pi = sel.PNum(j); + if (pi < points.End()) points[pi].SetType (FIXEDPOINT); } } @@ -2985,8 +2983,7 @@ namespace netgen pmin = Point3d (1e10, 1e10, 1e10); pmax = Point3d (-1e10, -1e10, -1e10); - for (PointIndex pi = PointIndex::BASE; - pi < GetNP()+PointIndex::BASE; pi++) + for (PointIndex pi = points.Begin(); pi < points.End(); pi++) { pmin.SetToMin ( (*this) [pi] ); pmax.SetToMax ( (*this) [pi] ); @@ -3035,8 +3032,7 @@ namespace netgen pmin = Point3d (1e10, 1e10, 1e10); pmax = Point3d (-1e10, -1e10, -1e10); - for (PointIndex pi = PointIndex::BASE; - pi < GetNP()+PointIndex::BASE; pi++) + for (PointIndex pi = points.Begin(); pi < points.End(); pi++) if (points[pi].Type() <= ptyp) { pmin.SetToMin ( (*this) [pi] ); @@ -3068,8 +3064,7 @@ namespace netgen void Mesh :: Compress () { - int i, j; - Array op2np(GetNP()); + Array op2np(GetNP()); Array hpoints; BitArrayChar pused(GetNP()); @@ -3084,7 +3079,7 @@ namespace netgen (*testout) << "np: " << GetNP() << endl; */ - for (i = 0; i < volelements.Size(); i++) + for (int i = 0; i < volelements.Size(); i++) if (volelements[i][0] <= PointIndex::BASE-1 || volelements[i].IsDeleted()) { @@ -3093,14 +3088,14 @@ namespace netgen } - for (i = 0; i < surfelements.Size(); i++) + for (int i = 0; i < surfelements.Size(); i++) if (surfelements[i].IsDeleted()) { surfelements.Delete(i); i--; } - for (i = 0; i < segments.Size(); i++) + for (int i = 0; i < segments.Size(); i++) if (segments[i][0] <= PointIndex::BASE-1) { segments.Delete(i); @@ -3108,35 +3103,35 @@ namespace netgen } pused.Clear(); - for (i = 0; i < volelements.Size(); i++) + for (int i = 0; i < volelements.Size(); i++) { const Element & el = volelements[i]; - for (j = 0; j < el.GetNP(); j++) + for (int j = 0; j < el.GetNP(); j++) pused.Set (el[j]); } - for (i = 0; i < surfelements.Size(); i++) + for (int i = 0; i < surfelements.Size(); i++) { const Element2d & el = surfelements[i]; - for (j = 0; j < el.GetNP(); j++) + for (int j = 0; j < el.GetNP(); j++) pused.Set (el[j]); } - for (i = 0; i < segments.Size(); i++) + for (int i = 0; i < segments.Size(); i++) { const Segment & seg = segments[i]; pused.Set (seg[0]); pused.Set (seg[1]); } - for (i = 0; i < openelements.Size(); i++) + for (int i = 0; i < openelements.Size(); i++) { const Element2d & el = openelements[i]; - for (j = 0; j < el.GetNP(); j++) + for (int j = 0; j < el.GetNP(); j++) pused.Set(el[j]); } - for (i = 0; i < lockedpoints.Size(); i++) + for (int i = 0; i < lockedpoints.Size(); i++) pused.Set (lockedpoints[i]); @@ -3157,54 +3152,52 @@ namespace netgen int npi = PointIndex::BASE-1; - for (i = PointIndex::BASE; - i < points.Size()+PointIndex::BASE; i++) - if (pused.Test(i)) + for (PointIndex pi = points.Begin(); pi < points.End(); pi++) + if (pused.Test(pi)) { npi++; - op2np[i] = npi; - hpoints.Append (points[i]); + op2np[pi] = npi; + hpoints.Append (points[pi]); } else - op2np[i] = -1; + op2np[pi] = -1; points.SetSize(0); - for (i = 0; i < hpoints.Size(); i++) + for (int i = 0; i < hpoints.Size(); i++) points.Append (hpoints[i]); - - for (i = 1; i <= volelements.Size(); i++) + for (int i = 1; i <= volelements.Size(); i++) { Element & el = VolumeElement(i); - for (j = 0; j < el.GetNP(); j++) + for (int j = 0; j < el.GetNP(); j++) el[j] = op2np[el[j]]; } - for (i = 1; i <= surfelements.Size(); i++) + for (int i = 1; i <= surfelements.Size(); i++) { Element2d & el = SurfaceElement(i); - for (j = 0; j < el.GetNP(); j++) + for (int j = 0; j < el.GetNP(); j++) el[j] = op2np[el[j]]; } - for (i = 0; i < segments.Size(); i++) + for (int i = 0; i < segments.Size(); i++) { Segment & seg = segments[i]; seg[0] = op2np[seg[0]]; seg[1] = op2np[seg[1]]; } - for (i = 1; i <= openelements.Size(); i++) + for (int i = 1; i <= openelements.Size(); i++) { Element2d & el = openelements.Elem(i); - for (j = 0; j < el.GetNP(); j++) + for (int j = 0; j < el.GetNP(); j++) el[j] = op2np[el[j]]; } - for (i = 0; i < lockedpoints.Size(); i++) + for (int i = 0; i < lockedpoints.Size(); i++) lockedpoints[i] = op2np[lockedpoints[i]]; for (int i = 0; i < facedecoding.Size(); i++) @@ -3216,7 +3209,6 @@ namespace netgen facedecoding[ind-1].firstelement = i; } - CalcSurfacesOfNode(); @@ -4089,75 +4081,59 @@ namespace netgen void Mesh :: BuildElementSearchTree () { - if (elementsearchtreets == GetTimeStamp()) - return; + if (elementsearchtreets == GetTimeStamp()) return; - NgLock lock(mutex); - lock.Lock(); - - PrintMessage (4, "Rebuild element searchtree"); - - delete elementsearchtree; - elementsearchtree = NULL; - - Box3d box; - int ne = (dimension == 2) ? GetNSE() : GetNE(); - if (!ne) - { - lock.UnLock(); - return; - } - - if (dimension == 2) - { - box.SetPoint (Point (SurfaceElement(1).PNum(1))); - for (int i = 1; i <= ne; i++) - { - const Element2d & el = SurfaceElement(i); - for (int j = 1; j <= el.GetNP(); j++) - box.AddPoint (Point (el.PNum(j))); - } - - box.Increase (1.01 * box.CalcDiam()); - elementsearchtree = new Box3dTree (box.PMin(), box.PMax()); - - for (int i = 1; i <= ne; i++) - { - const Element2d & el = SurfaceElement(i); - box.SetPoint (Point (el.PNum(1))); - for (int j = 1; j <= el.GetNP(); j++) - box.AddPoint (Point (el.PNum(j))); - - elementsearchtree -> Insert (box.PMin(), box.PMax(), i); - } - } - else - { - box.SetPoint (Point (VolumeElement(1).PNum(1))); - for (int i = 1; i <= ne; i++) - { - const Element & el = VolumeElement(i); - for (int j = 1; j <= el.GetNP(); j++) - box.AddPoint (Point (el.PNum(j))); - } - - box.Increase (1.01 * box.CalcDiam()); - elementsearchtree = new Box3dTree (box.PMin(), box.PMax()); - - for (int i = 1; i <= ne; i++) - { - const Element & el = VolumeElement(i); - box.SetPoint (Point (el.PNum(1))); - for (int j = 1; j <= el.GetNP(); j++) - box.AddPoint (Point (el.PNum(j))); - - elementsearchtree -> Insert (box.PMin(), box.PMax(), i); - } - } - - elementsearchtreets = GetTimeStamp(); - - lock.UnLock(); +#pragma omp critical (buildsearchtree) + { + if (elementsearchtreets != GetTimeStamp()) + { + NgLock lock(mutex); + lock.Lock(); + + PrintMessage (4, "Rebuild element searchtree"); + + delete elementsearchtree; + elementsearchtree = NULL; + + int ne = (dimension == 2) ? GetNSE() : GetNE(); + + if (ne) + { + if (dimension == 2) + { + Box<3> box (Box<3>::EMPTY_BOX); + for (SurfaceElementIndex sei = 0; sei < ne; sei++) + box.Add (points[surfelements[sei].PNums()]); + + box.Increase (1.01 * box.Diam()); + elementsearchtree = new Box3dTree (box); + + for (SurfaceElementIndex sei = 0; sei < ne; sei++) + { + box.Set (points[surfelements[sei].PNums()]); + elementsearchtree -> Insert (box, sei+1); + } + } + else + { + Box<3> box (Box<3>::EMPTY_BOX); + for (ElementIndex ei = 0; ei < ne; ei++) + box.Add (points[volelements[ei].PNums()]); + + box.Increase (1.01 * box.Diam()); + elementsearchtree = new Box3dTree (box); + + for (ElementIndex ei = 0; ei < ne; ei++) + { + box.Set (points[volelements[ei].PNums()]); + elementsearchtree -> Insert (box, ei+1); + } + } + + elementsearchtreets = GetTimeStamp(); + } + } + } } @@ -4507,7 +4483,7 @@ namespace netgen } - int Mesh :: GetElementOfPoint (const Point3d & p, + int Mesh :: GetElementOfPoint (const netgen::Point<3> & p, double lami[3], bool build_searchtree, const int index, @@ -4526,7 +4502,7 @@ namespace netgen - int Mesh :: GetElementOfPoint (const Point3d & p, + int Mesh :: GetElementOfPoint (const netgen::Point<3> & p, double lami[3], const Array * const indices, bool build_searchtree, @@ -4649,7 +4625,7 @@ namespace netgen - int Mesh :: GetSurfaceElementOfPoint (const Point3d & p, + int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p, double lami[3], bool build_searchtree, const int index, @@ -4668,7 +4644,7 @@ namespace netgen - int Mesh :: GetSurfaceElementOfPoint (const Point3d & p, + int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p, double lami[3], const Array * const indices, bool build_searchtree,