Handle pyramids in smoothing

This commit is contained in:
Matthias Hochsteger 2020-07-23 12:19:12 +02:00
parent e17de17385
commit bb46dd6561

View File

@ -339,9 +339,9 @@ namespace netgen
{ {
static Timer tim("PointFunction - build elementsonpoint table"); RegionTimer reg(tim); static Timer tim("PointFunction - build elementsonpoint table"); RegionTimer reg(tim);
for (int i = 0; i < elements.Size(); i++) for (int i = 0; i < elements.Size(); i++)
if (elements[i].NP() == 4) if (!elements[i].IsDeleted())
for (int j = 0; j < elements[i].NP(); j++) for (int j = 0; j < elements[i].GetNP(); j++)
elementsonpoint.Add (elements[i][j], i); elementsonpoint.Add (elements[i][j], i);
} }
void PointFunction :: SetPointIndex (PointIndex aactpind) void PointFunction :: SetPointIndex (PointIndex aactpind)
@ -362,8 +362,7 @@ namespace netgen
for (int j = 0; j < elementsonpoint[actpind].Size(); j++) for (int j = 0; j < elementsonpoint[actpind].Size(); j++)
{ {
const Element & el = elements[elementsonpoint[actpind][j]]; const Element & el = elements[elementsonpoint[actpind][j]];
badness += CalcTetBadness (points[el[0]], points[el[1]], badness += CalcBad(points, el, -1, mp);
points[el[2]], points[el[3]], -1, mp);
} }
points[actpind] = Point<3> (hp); points[actpind] = Point<3> (hp);
@ -391,6 +390,19 @@ namespace netgen
vgrad += vgradi; vgrad += vgradi;
} }
if(el.GetType()==PYRAMID)
{
f += CalcTetBadnessGrad (points[el[0]],
points[el[1]],
points[el[2]],
points[el[4]], -1, 4, vgradi, mp);
vgrad += vgradi;
f += CalcTetBadnessGrad (points[el[2]],
points[el[3]],
points[el[0]],
points[el[4]], -1, 4, vgradi, mp);
vgrad += vgradi;
}
} }
points[actpind] = Point<3> (hp); points[actpind] = Point<3> (hp);
@ -423,6 +435,20 @@ namespace netgen
vgrad += vgradi; vgrad += vgradi;
} }
if(el.GetType()==PYRAMID)
{
f += CalcTetBadnessGrad (points[el[0]],
points[el[1]],
points[el[2]],
points[el[4]], -1, 4, vgradi, mp);
vgrad += vgradi;
f += CalcTetBadnessGrad (points[el[2]],
points[el[3]],
points[el[0]],
points[el[4]], -1, 4, vgradi, mp);
vgrad += vgradi;
}
} }
points[actpind] = Point<3> (hp); points[actpind] = Point<3> (hp);