ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-26 13:50:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

thread-save build-searchtree

Browse Source
This commit is contained in:
Joachim Schoeberl 2013-04-02 20:29:05 +00:00
parent 641d40639d
commit 15a1e07092
1 changed files with 107 additions and 131 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

194
libsrc/meshing/meshclass.cpp
View File

@ -150,17 +150,20 @@ namespace netgen
PointIndex Mesh :: AddPoint (const Point3d & p, int layer) PointIndex Mesh :: AddPoint (const Point3d & p, int layer)
{ {
return AddPoint (p, layer, INNERPOINT);
/*
NgLock lock(mutex); NgLock lock(mutex);
lock.Lock(); lock.Lock();
timestamp = NextTimeStamp(); timestamp = NextTimeStamp();
PointIndex pi = points.Size() + PointIndex::BASE; PointIndex pi = points.End();
points.Append ( MeshPoint (p, layer, INNERPOINT) ); points.Append ( MeshPoint (p, layer, INNERPOINT) );
lock.UnLock(); lock.UnLock();
return pi; return pi;
*/
} }
PointIndex Mesh :: AddPoint (const Point3d & p, int layer, POINTTYPE type) PointIndex Mesh :: AddPoint (const Point3d & p, int layer, POINTTYPE type)
@ -170,7 +173,7 @@ namespace netgen
timestamp = NextTimeStamp(); timestamp = NextTimeStamp();
PointIndex pi = points.Size() + PointIndex::BASE; PointIndex pi = points.End();
points.Append ( MeshPoint (p, layer, type) ); points.Append ( MeshPoint (p, layer, type) );
lock.UnLock(); lock.UnLock();
@ -676,13 +679,13 @@ namespace netgen
*/ */
int cnt_sing = 0; 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 ((*this)[pi].Singularity()>=1.) cnt_sing++;
if (cnt_sing) if (cnt_sing)
{ {
outfile << "singular_points" << endl << cnt_sing << endl; 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.) if ((*this)[pi].Singularity()>=1.)
outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl; outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl;
} }
@ -1509,9 +1512,6 @@ namespace netgen
void Mesh :: CalcSurfacesOfNode () void Mesh :: CalcSurfacesOfNode ()
{ {
int i, j, k;
SurfaceElementIndex sei;
surfacesonnode.SetSize (GetNP()); surfacesonnode.SetSize (GetNP());
delete boundaryedges; delete boundaryedges;
@ -1528,18 +1528,18 @@ namespace netgen
surfelementht = new INDEX_3_CLOSED_HASHTABLE<int> (3*GetNSE() + 1); surfelementht = new INDEX_3_CLOSED_HASHTABLE<int> (3*GetNSE() + 1);
segmentht = new INDEX_2_CLOSED_HASHTABLE<int> (3*GetNSeg() + 1); segmentht = new INDEX_2_CLOSED_HASHTABLE<int> (3*GetNSeg() + 1);
for (sei = 0; sei < GetNSE(); sei++) for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
{ {
const Element2d & sel = surfelements[sei]; const Element2d & sel = surfelements[sei];
if (sel.IsDeleted()) continue; if (sel.IsDeleted()) continue;
int si = sel.GetIndex(); int si = sel.GetIndex();
for (j = 0; j < sel.GetNP(); j++) for (int j = 0; j < sel.GetNP(); j++)
{ {
PointIndex pi = sel[j]; PointIndex pi = sel[j];
bool found = 0; 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) if (surfacesonnode[pi][k] == si)
{ {
found = 1; found = 1;
@ -1548,7 +1548,6 @@ namespace netgen
if (!found) if (!found)
surfacesonnode.Add (pi, si); surfacesonnode.Add (pi, si);
} }
} }
/* /*
@ -1567,7 +1566,7 @@ namespace netgen
surfelementht -> AllocateElements(); surfelementht -> AllocateElements();
*/ */
for (sei = 0; sei < GetNSE(); sei++) for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
{ {
const Element2d & sel = surfelements[sei]; const Element2d & sel = surfelements[sei];
if (sel.IsDeleted()) continue; if (sel.IsDeleted()) continue;
@ -1584,17 +1583,16 @@ namespace netgen
if (dimension == 3) if (dimension == 3)
{ {
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
pi < np+PointIndex::BASE; pi++)
points[pi].SetType (INNERPOINT); points[pi].SetType (INNERPOINT);
if (GetNFD() == 0) if (GetNFD() == 0)
{ {
for (sei = 0; sei < GetNSE(); sei++) for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
{ {
const Element2d & sel = surfelements[sei]; const Element2d & sel = surfelements[sei];
if (sel.IsDeleted()) continue; if (sel.IsDeleted()) continue;
for (j = 0; j < sel.GetNP(); j++) for (int j = 0; j < sel.GetNP(); j++)
{ {
PointIndex pi = SurfaceElement(sei)[j]; PointIndex pi = SurfaceElement(sei)[j];
points[pi].SetType(FIXEDPOINT); points[pi].SetType(FIXEDPOINT);
@ -1603,11 +1601,11 @@ namespace netgen
} }
else else
{ {
for (sei = 0; sei < GetNSE(); sei++) for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
{ {
const Element2d & sel = surfelements[sei]; const Element2d & sel = surfelements[sei];
if (sel.IsDeleted()) continue; if (sel.IsDeleted()) continue;
for (j = 0; j < sel.GetNP(); j++) for (int j = 0; j < sel.GetNP(); j++)
{ {
PointIndex pi = sel[j]; PointIndex pi = sel[j];
int ns = surfacesonnode[pi].Size(); 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]; 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]; 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); points[lockedpoints[i]].SetType(FIXEDPOINT);
} }
@ -1660,7 +1658,7 @@ namespace netgen
// eltyps.SetSize (GetNE()); // eltyps.SetSize (GetNE());
// eltyps = FREEELEMENT; // eltyps = FREEELEMENT;
for (i = 0; i < GetNSeg(); i++) for (int i = 0; i < GetNSeg(); i++)
{ {
const Segment & seg = segments[i]; const Segment & seg = segments[i];
INDEX_2 i2(seg[0], seg[1]); INDEX_2 i2(seg[0], seg[1]);
@ -1820,7 +1818,7 @@ namespace netgen
INDEX_3_CLOSED_HASHTABLE<INDEX_2> faceht(100); INDEX_3_CLOSED_HASHTABLE<INDEX_2> faceht(100);
openelements.SetSize(0); 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()) if (selsonpoint[pi].Size()+elsonpoint[pi].Size())
{ {
faceht.SetSize (2 * selsonpoint[pi].Size() + 4 * 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++) for (int j = 1; j <= 3; j++)
{ {
int pi = sel.PNum(j); PointIndex pi = sel.PNum(j);
if (pi < points.Size()+PointIndex::BASE) if (pi < points.End())
points[pi].SetType (FIXEDPOINT); points[pi].SetType (FIXEDPOINT);
} }
} }
@ -2985,8 +2983,7 @@ namespace netgen
pmin = Point3d (1e10, 1e10, 1e10); pmin = Point3d (1e10, 1e10, 1e10);
pmax = Point3d (-1e10, -1e10, -1e10); pmax = Point3d (-1e10, -1e10, -1e10);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
pi < GetNP()+PointIndex::BASE; pi++)
{ {
pmin.SetToMin ( (*this) [pi] ); pmin.SetToMin ( (*this) [pi] );
pmax.SetToMax ( (*this) [pi] ); pmax.SetToMax ( (*this) [pi] );
@ -3035,8 +3032,7 @@ namespace netgen
pmin = Point3d (1e10, 1e10, 1e10); pmin = Point3d (1e10, 1e10, 1e10);
pmax = Point3d (-1e10, -1e10, -1e10); pmax = Point3d (-1e10, -1e10, -1e10);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
pi < GetNP()+PointIndex::BASE; pi++)
if (points[pi].Type() <= ptyp) if (points[pi].Type() <= ptyp)
{ {
pmin.SetToMin ( (*this) [pi] ); pmin.SetToMin ( (*this) [pi] );
@ -3068,8 +3064,7 @@ namespace netgen
void Mesh :: Compress () void Mesh :: Compress ()
{ {
int i, j; Array<PointIndex,PointIndex::BASE,PointIndex> op2np(GetNP());
Array<int,PointIndex::BASE> op2np(GetNP());
Array<MeshPoint> hpoints; Array<MeshPoint> hpoints;
BitArrayChar<PointIndex::BASE> pused(GetNP()); BitArrayChar<PointIndex::BASE> pused(GetNP());
@ -3084,7 +3079,7 @@ namespace netgen
(*testout) << "np: " << GetNP() << endl; (*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 || if (volelements[i][0] <= PointIndex::BASE-1 ||
volelements[i].IsDeleted()) 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()) if (surfelements[i].IsDeleted())
{ {
surfelements.Delete(i); surfelements.Delete(i);
i--; i--;
} }
for (i = 0; i < segments.Size(); i++) for (int i = 0; i < segments.Size(); i++)
if (segments[i][0] <= PointIndex::BASE-1) if (segments[i][0] <= PointIndex::BASE-1)
{ {
segments.Delete(i); segments.Delete(i);
@ -3108,35 +3103,35 @@ namespace netgen
} }
pused.Clear(); pused.Clear();
for (i = 0; i < volelements.Size(); i++) for (int i = 0; i < volelements.Size(); i++)
{ {
const Element & el = volelements[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]); pused.Set (el[j]);
} }
for (i = 0; i < surfelements.Size(); i++) for (int i = 0; i < surfelements.Size(); i++)
{ {
const Element2d & el = surfelements[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]); pused.Set (el[j]);
} }
for (i = 0; i < segments.Size(); i++) for (int i = 0; i < segments.Size(); i++)
{ {
const Segment & seg = segments[i]; const Segment & seg = segments[i];
pused.Set (seg[0]); pused.Set (seg[0]);
pused.Set (seg[1]); pused.Set (seg[1]);
} }
for (i = 0; i < openelements.Size(); i++) for (int i = 0; i < openelements.Size(); i++)
{ {
const Element2d & el = openelements[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]); pused.Set(el[j]);
} }
for (i = 0; i < lockedpoints.Size(); i++) for (int i = 0; i < lockedpoints.Size(); i++)
pused.Set (lockedpoints[i]); pused.Set (lockedpoints[i]);
@ -3157,54 +3152,52 @@ namespace netgen
int npi = PointIndex::BASE-1; int npi = PointIndex::BASE-1;
for (i = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
i < points.Size()+PointIndex::BASE; i++) if (pused.Test(pi))
if (pused.Test(i))
{ {
npi++; npi++;
op2np[i] = npi; op2np[pi] = npi;
hpoints.Append (points[i]); hpoints.Append (points[pi]);
} }
else else
op2np[i] = -1; op2np[pi] = -1;
points.SetSize(0); points.SetSize(0);
for (i = 0; i < hpoints.Size(); i++) for (int i = 0; i < hpoints.Size(); i++)
points.Append (hpoints[i]); points.Append (hpoints[i]);
for (int i = 1; i <= volelements.Size(); i++)
for (i = 1; i <= volelements.Size(); i++)
{ {
Element & el = VolumeElement(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]]; el[j] = op2np[el[j]];
} }
for (i = 1; i <= surfelements.Size(); i++) for (int i = 1; i <= surfelements.Size(); i++)
{ {
Element2d & el = SurfaceElement(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]]; el[j] = op2np[el[j]];
} }
for (i = 0; i < segments.Size(); i++) for (int i = 0; i < segments.Size(); i++)
{ {
Segment & seg = segments[i]; Segment & seg = segments[i];
seg[0] = op2np[seg[0]]; seg[0] = op2np[seg[0]];
seg[1] = op2np[seg[1]]; 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); 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]]; 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]]; lockedpoints[i] = op2np[lockedpoints[i]];
for (int i = 0; i < facedecoding.Size(); i++) for (int i = 0; i < facedecoding.Size(); i++)
@ -3216,7 +3209,6 @@ namespace netgen
facedecoding[ind-1].firstelement = i; facedecoding[ind-1].firstelement = i;
} }
CalcSurfacesOfNode(); CalcSurfacesOfNode();
@ -4089,9 +4081,12 @@ namespace netgen
void Mesh :: BuildElementSearchTree () void Mesh :: BuildElementSearchTree ()
{ {
if (elementsearchtreets == GetTimeStamp()) if (elementsearchtreets == GetTimeStamp()) return;
return;
#pragma omp critical (buildsearchtree)
{
if (elementsearchtreets != GetTimeStamp())
{
NgLock lock(mutex); NgLock lock(mutex);
lock.Lock(); lock.Lock();
@ -4100,64 +4095,45 @@ namespace netgen
delete elementsearchtree; delete elementsearchtree;
elementsearchtree = NULL; elementsearchtree = NULL;
Box3d box;
int ne = (dimension == 2) ? GetNSE() : GetNE(); int ne = (dimension == 2) ? GetNSE() : GetNE();
if (!ne)
{
lock.UnLock();
return;
}
if (ne)
{
if (dimension == 2) if (dimension == 2)
{ {
box.SetPoint (Point (SurfaceElement(1).PNum(1))); Box<3> box (Box<3>::EMPTY_BOX);
for (int i = 1; i <= ne; i++) 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++)
{ {
const Element2d & el = SurfaceElement(i); box.Set (points[surfelements[sei].PNums()]);
for (int j = 1; j <= el.GetNP(); j++) elementsearchtree -> Insert (box, sei+1);
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 else
{ {
box.SetPoint (Point (VolumeElement(1).PNum(1))); Box<3> box (Box<3>::EMPTY_BOX);
for (int i = 1; i <= ne; i++) 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++)
{ {
const Element & el = VolumeElement(i); box.Set (points[volelements[ei].PNums()]);
for (int j = 1; j <= el.GetNP(); j++) elementsearchtree -> Insert (box, ei+1);
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(); elementsearchtreets = GetTimeStamp();
}
lock.UnLock(); }
}
} }
@ -4507,7 +4483,7 @@ namespace netgen
} }
int Mesh :: GetElementOfPoint (const Point3d & p, int Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
double lami[3], double lami[3],
bool build_searchtree, bool build_searchtree,
const int index, 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], double lami[3],
const Array<int> * const indices, const Array<int> * const indices,
bool build_searchtree, 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], double lami[3],
bool build_searchtree, bool build_searchtree,
const int index, 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], double lami[3],
const Array<int> * const indices, const Array<int> * const indices,
bool build_searchtree, bool build_searchtree,