Skip to content
GitLab
Commit 05d09d4b authored by David Thompson's avatar David Thompson
Browse files

Implement 21-point wedge; fix others.

parent 0a2d95aa
Hide whitespace changes
Inline Side-by-side
Common/DataModel/Testing/Cxx/LagrangeHexahedron.cxx
View file @ 05d09d4b
......@@ -83,12 +83,12 @@ static const double expectedFacePoints333[96][3] = {
{ 0, 1, 1 },
{ 0, 0.666667, 0 },
{ 0, 0.333333, 0 },
{ 0, 0, 0.666667 },
{ 0, 0, 0.333333 },
{ 0, 0, 0.666667 },
{ 0, 0.666667, 1 },
{ 0, 0.333333, 1 },
{ 0, 1, 0.666667 },
{ 0, 1, 0.333333 },
{ 0, 1, 0.666667 },
{ 0, 0.666667, 0.333333 },
{ 0, 0.333333, 0.333333 },
{ 0, 0.666667, 0.666667 },
......@@ -111,39 +111,39 @@ static const double expectedFacePoints333[96][3] = {
{ 1, 0.333333, 0.666667 },
{ 1, 0.666667, 0.666667 },
{ 1, 0, 0 },
{ 0, 0, 0 },
{ 0, 0, 1 },
{ 1, 0, 0 },
{ 1, 0, 1 },
{ 0.666667, 0, 0 },
{ 0, 0, 1 },
{ 0.333333, 0, 0 },
{ 0, 0, 0.666667 },
{ 0, 0, 0.333333 },
{ 0.666667, 0, 1 },
{ 0.333333, 0, 1 },
{ 1, 0, 0.666667 },
{ 0.666667, 0, 0 },
{ 1, 0, 0.333333 },
{ 0.666667, 0, 0.333333 },
{ 1, 0, 0.666667 },
{ 0.333333, 0, 1 },
{ 0.666667, 0, 1 },
{ 0, 0, 0.333333 },
{ 0, 0, 0.666667 },
{ 0.333333, 0, 0.333333 },
{ 0.666667, 0, 0.666667 },
{ 0.666667, 0, 0.333333 },
{ 0.333333, 0, 0.666667 },
{ 0.666667, 0, 0.666667 },
{ 0, 1, 0 },
{ 1, 1, 0 },
{ 1, 1, 1 },
{ 0, 1, 0 },
{ 0, 1, 1 },
{ 0.333333, 1, 0 },
{ 1, 1, 1 },
{ 0.666667, 1, 0 },
{ 1, 1, 0.333333 },
{ 1, 1, 0.666667 },
{ 0.333333, 1, 1 },
{ 0.666667, 1, 1 },
{ 0.333333, 1, 0 },
{ 0, 1, 0.333333 },
{ 0, 1, 0.666667 },
{ 0.333333, 1, 0.333333 },
{ 0.666667, 1, 1 },
{ 0.333333, 1, 1 },
{ 1, 1, 0.333333 },
{ 1, 1, 0.666667 },
{ 0.666667, 1, 0.333333 },
{ 0.333333, 1, 0.666667 },
{ 0.333333, 1, 0.333333 },
{ 0.666667, 1, 0.666667 },
{ 0.333333, 1, 0.666667 },
{ 1, 0, 0 },
{ 0, 0, 0 },
......@@ -151,12 +151,12 @@ static const double expectedFacePoints333[96][3] = {
{ 1, 1, 0 },
{ 0.666667, 0, 0 },
{ 0.333333, 0, 0 },
{ 0, 0.666667, 0 },
{ 0, 0.333333, 0 },
{ 0, 0.666667, 0 },
{ 0.666667, 1, 0 },
{ 0.333333, 1, 0 },
{ 1, 0.666667, 0 },
{ 1, 0.333333, 0 },
{ 1, 0.666667, 0 },
{ 0.666667, 0.333333, 0 },
{ 0.333333, 0.333333, 0 },
{ 0.666667, 0.666667, 0 },
......
Common/DataModel/vtkLagrangeCurve.cxx
View file @ 05d09d4b
......@@ -229,11 +229,11 @@ int vtkLagrangeCurve::IntersectWithLine(
if (!intersection || (t >= 0 && (t < tFirst || tFirst < 0)))
{
tFirst = t;
subId = i;
for (int ii = 0; ii < 3; ++ii)
{
x[ii] = tmpX[ii];
pcoords[ii] = tmpP[ii]; // Translate this after we're sure it's the closest hit.
subId = i;
}
}
intersection = true;
......@@ -241,7 +241,8 @@ int vtkLagrangeCurve::IntersectWithLine(
}
if (intersection)
{
this->TransformApproxToCellParams(subId, pcoords);
intersection &= this->TransformApproxToCellParams(subId, pcoords);
t = tFirst;
}
return intersection ? 1 : 0;
}
......
Common/DataModel/vtkLagrangeHexahedron.cxx
View file @ 05d09d4b
......@@ -104,12 +104,12 @@ vtkCell* vtkLagrangeHexahedron::GetEdge(int edgeId)
vtkCell* vtkLagrangeHexahedron::GetFace(int faceId)
{
if (faceId < 0 || faceId >= 6)
{
{
return nullptr;
}
}
// Do we need to flip the face to get an outward-pointing normal?
bool flipFace = (faceId % 2 == 0 ? true : false);
bool flipFace = (faceId % 2 == ((faceId / 2) % 2) ? true : false);
vtkLagrangeQuadrilateral* result = this->FaceCell.GetPointer();
const int* order = this->GetOrder();
......@@ -122,112 +122,123 @@ vtkCell* vtkLagrangeHexahedron::GetFace(int faceId)
// Add vertex DOFs to result
int sn = 0;
if (!flipFace)
{
{
for (int ii = 0; ii < 4; ++ii, ++sn)
{
{
result->Points->SetPoint(sn, this->Points->GetPoint(corners[ii]));
result->PointIds->SetId(sn, this->PointIds->GetId(corners[ii]));
}
}
}
else
{
{
for (int ii = 0; ii < 4; ++ii, ++sn)
{
{
result->Points->SetPoint((5 - sn) % 4, this->Points->GetPoint(corners[ii]));
result->PointIds->SetId((5 - sn) % 4, this->PointIds->GetId(corners[ii]));
}
}
}
// Add edge DOFs to result
int offset;
const int* faceEdges = vtkLagrangeInterpolation::GetEdgeIndicesBoundingHexFace(faceId);
for (int ii = 0; ii < 4; ++ii)
{
{
offset = 8;
if (!flipFace)
{
{
int pp = vtkLagrangeInterpolation::GetVaryingParameterOfHexEdge(faceEdges[ii]);
if (pp == 2)
{
{
offset += 4 * (order[0] - 1 + order[1] - 1);
offset += (faceEdges[ii] - 8) * (order[2] - 1);
}
}
else
{
{
for (int ee = 0; ee < faceEdges[ii]; ++ee)
{
{
offset += order[ee % 2 == 0 ? 0 : 1] - 1;
}
}
}
for (int jj = 0; jj < order[pp] - 1; ++jj, ++sn)
{
{
result->Points->SetPoint(sn, this->Points->GetPoint(offset + jj));
result->PointIds->SetId(sn, this->PointIds->GetId(offset + jj));
}
}
}
else
{
{
// Flip both the edge position among edges (ii => (4 - ii) % 4)
// and the edge's node order (jj => order[pp] - jj - 1).
int pp = vtkLagrangeInterpolation::GetVaryingParameterOfHexEdge(faceEdges[(4 - ii) % 4]);
if (pp == 2)
{
{
offset += 4 * (order[0] - 1 + order[1] - 1);
offset += (faceEdges[(4 - ii) % 4] - 8) * (order[2] - 1);
}
}
else
{
{
for (int ee = 0; ee < faceEdges[(4 - ii) % 4]; ++ee)
{
{
offset += order[ee % 2 == 0 ? 0 : 1] - 1;
}
}
for (int jj = 0; jj < order[pp] - 1; ++jj, ++sn)
}
if (ii % 2 == 0)
{
for (int jj = 0; jj < order[pp] - 1; ++jj, ++sn)
{
result->Points->SetPoint(sn, this->Points->GetPoint(offset + order[pp] - jj - 2));
result->PointIds->SetId(sn, this->PointIds->GetId(offset + order[pp] - jj - 2));
}
}
else
{
for (int jj = 0; jj < order[pp] - 1; ++jj, ++sn)
{
result->Points->SetPoint(sn, this->Points->GetPoint(offset + order[pp] - jj - 2));
result->PointIds->SetId(sn, this->PointIds->GetId(offset + order[pp] - jj - 2));
result->Points->SetPoint(sn, this->Points->GetPoint(offset + jj));
result->PointIds->SetId(sn, this->PointIds->GetId(offset + jj));
}
}
}
}
// Now add face DOF
offset = 8 + 4 * (order[0] - 1 + order[1] - 1 + order[2] - 1);
// skip DOF for other faces of hex before this one
for (int ff = 0; ff < faceId; ++ff)
{
{
vtkVector2i tmp = vtkLagrangeInterpolation::GetVaryingParametersOfHexFace(ff);
offset += (order[tmp[0]] - 1) * (order[tmp[1]] - 1);
}
}
if (!flipFace)
{
{
int nfdof = (order[faceParams[0]] - 1) * (order[faceParams[1]] - 1);
for (int ii = 0; ii < nfdof; ++ii, ++sn)
{
{
result->Points->SetPoint(sn, this->Points->GetPoint(offset + ii));
result->PointIds->SetId(sn, this->PointIds->GetId(offset + ii));
}
}
}
else
{
{
int delta = order[faceParams[0]] - 1;
for (int jj = 0; jj < (order[faceParams[1]] - 1); ++jj)
{
{
for (int ii = delta - 1; ii >= 0; --ii, ++sn)
{
{
result->Points->SetPoint(sn, this->Points->GetPoint(offset + ii + jj * delta));
result->PointIds->SetId(sn, this->PointIds->GetId(offset + ii + jj * delta));
}
}
}
}
/*
std::cout << "Hex Face " << faceId << "\n";
for (int yy = 0; yy < npts; ++yy)
{
vtkVector3d xx;
result->Points->GetPoint(yy, xx.GetData());
std::cout << " " << yy << " " << result->PointIds->GetId(yy) << " " << xx << "\n";
}
*/
{
vtkVector3d xx;
result->Points->GetPoint(yy, xx.GetData());
std::cout << " " << yy << " " << result->PointIds->GetId(yy) << " " << xx << "\n";
}
*/
return result;
}
......@@ -449,27 +460,28 @@ int vtkLagrangeHexahedron::IntersectWithLine(
int tmpId;
this->GetOrder(); // Ensure Order is up to date.
for (int ff = 0; ff < this->GetNumberOfFaces(); ++ff)
{
{
vtkCell* bdy = this->GetFace(ff);
if (bdy->IntersectWithLine(p1, p2, tol, t, tmpX.GetData(), tmpP.GetData(), tmpId))
{
{
intersection = true;
if (t < tFirst)
{
{
tFirst = t;
subId = ff;
for (int ii = 0; ii < 3; ++ii)
{
{
x[ii] = tmpX[ii];
pcoords[ii] = tmpP[ii]; // Translate this after we're sure it's the closest hit.
subId = ff;
}
}
}
}
}
if (intersection)
{
this->TransformFaceToCellParams(subId, pcoords);
}
{
intersection &= this->TransformFaceToCellParams(subId, pcoords);
t = tFirst;
}
return intersection ? 1 : 0;
}
......@@ -512,13 +524,12 @@ int vtkLagrangeHexahedron::Triangulate(
void vtkLagrangeHexahedron::Derivatives(
int vtkNotUsed(subId),
double vtkNotUsed(pcoords)[3],
double* vtkNotUsed(values),
int vtkNotUsed(dim),
double* vtkNotUsed(derivs))
double pcoords[3],
double* values,
int dim,
double* derivs)
{
// TODO: Fill me in?
return;
this->Interp->Tensor3EvaluateDerivative(this->Order, pcoords, values, dim, derivs);
}
double* vtkLagrangeHexahedron::GetParametricCoords()
......@@ -820,7 +831,7 @@ bool vtkLagrangeHexahedron::TransformFaceToCellParams(int bdyFace, double* pcoor
{
pcoords[faceParams[pp]] = tmp[pp];
}
if (bdyFace % 2 == 0)
if (bdyFace % 2 == ((bdyFace / 2) % 2))
{
// Flip first parametric axis of "positive" faces to compensate for GetFace,
// which flips odd faces to obtain inward-pointing normals for each boundary.
......
Common/DataModel/vtkLagrangeHexahedron.h
View file @ 05d09d4b
......@@ -30,6 +30,7 @@ class vtkDoubleArray;
class vtkHexahedron;
class vtkIdList;
class vtkLagrangeCurve;
class vtkLagrangeInterpolation;
class vtkLagrangeQuadrilateral;
class vtkPointData;
class vtkPoints;
......@@ -118,6 +119,7 @@ protected:
vtkNew<vtkIdList> TmpIds;
vtkNew<vtkLagrangeQuadrilateral> FaceCell;
vtkNew<vtkLagrangeCurve> EdgeCell;
vtkNew<vtkLagrangeInterpolation> Interp;
private:
vtkLagrangeHexahedron(const vtkLagrangeHexahedron&) = delete;
......
Common/DataModel/vtkLagrangeInterpolation.cxx
View file @ 05d09d4b
......@@ -608,6 +608,35 @@ int vtkLagrangeInterpolation::Tensor3ShapeDerivatives(const int order[2], const
return sn;
}
void vtkLagrangeInterpolation::Tensor3EvaluateDerivative(
const int order[4],
const double* pcoords,
double* fieldVals,
int fieldDim,
double* fieldDerivs)
{
this->PrepareForOrder(order);
this->Tensor3ShapeDerivatives(order, pcoords, &this->DerivSpace[0]);
// Loop over variables we differentiate with respect to:
for (int vv = 0; vv < 3; ++vv)
{
// Loop over components of the field:
for (int cc = 0; cc < fieldDim; ++cc)
{
// Loop over shape functions (per-DOF values of the cell):
fieldDerivs[3 * cc + vv] = 0.;
for (int pp = 0; pp < order[3]; ++pp)
{
// Note the subtle difference between the indexing of this->DerivSpace here and in WedgeEvaluateDerivative.
// TODO: Choose the better approach and be consistent.
fieldDerivs[3 * cc + vv] += this->DerivSpace[order[3] * pp + vv] * fieldVals[fieldDim * pp + cc];
}
}
}
// ((d(vx)/dx),(d(vx)/dy),(d(vx)/dz), (d(vy)/dx),(d(vy)/dy), (d(vy)/dz), (d(vz)/dx),(d(vz)/dy),(d(vz)/dz))
}
/// Wedge shape function computation
void vtkLagrangeInterpolation::WedgeShapeFunctions(const int order[3], const double pcoords[3], double* shape)
{
......@@ -632,6 +661,42 @@ void vtkLagrangeInterpolation::WedgeShapeFunctions(const int order[3], const dou
return;
}
#ifdef VTK_21_POINT_WEDGE
if (order[3] == 21 && order[0] == 2)
{
const double r = pcoords[0];
const double s = pcoords[1];
const double t = 2 * pcoords[2] - 1.;
const double rsm = 1. - r - s;
const double rs = r * s;
const double tp = 1. + t;
const double tm = 1. - t;
shape[ 0] = -0.5 * t * tm * rsm * (1.0 - 2.0 * (r + s) + 3.0 * rs);
shape[ 1] = -0.5 * t * tm * (r - 2.0 * (rsm * r + rs) + 3.0 * rsm * rs);
shape[ 2] = -0.5 * t * tm * (s - 2.0 * (rsm * s + rs) + 3.0 * rsm * rs);
shape[ 3] = 0.5 * t * tp * rsm * (1.0 - 2.0 * (r + s) + 3.0 * rs);
shape[ 4] = 0.5 * t * tp * (r - 2.0 * (rsm * r + rs) + 3.0 * rsm * rs);
shape[ 5] = 0.5 * t * tp * (s - 2.0 * (rsm * s + rs) + 3.0 * rsm * rs);
shape[ 6] = -0.5 * t * tm * rsm * (4.0 * r - 12.0 * rs);
shape[ 7] = -0.5 * t * tm * (4.0 * rs - 12.0 * rsm * rs);
shape[ 8] = -0.5 * t * tm * rsm * (4.0 * s - 12.0 * rs);
shape[ 9] = 0.5 * t * tp * rsm * (4.0 * r - 12.0 * rs);
shape[10] = 0.5 * t * tp * (4.0 * rs - 12.0 * rsm * rs);
shape[11] = 0.5 * t * tp * rsm * (4.0 * s - 12.0 * rs);
shape[12] = tp * tm * rsm * (1.0 - 2.0 * (r + s) + 3.0 * rs);
shape[13] = tp * tm * (r - 2.0 * (rsm * r + rs) + 3.0 * rsm * rs);
shape[14] = tp * tm * (s - 2.0 * (rsm * s + rs) + 3.0 * rsm * rs);
shape[15] = -0.5 * 27.0 * t * tm * rsm * rs;
shape[16] = 0.5 * 27.0 * t * tp * rsm * rs;
shape[17] = tp * tm * rsm * (4.0 * r - 12.0 * rs);
shape[18] = tp * tm * (4.0 * rs - 12.0 * rsm * rs);
shape[19] = tp * tm * rsm * (4.0 * s - 12.0 * rs);
shape[20] = 27.0 * tp * tm * rsm * rs;
return;
}
#endif
// FIXME: Eventually needs to be varying length.
double ll[vtkLagrangeInterpolation::MaxDegree + 1];
double tt[(vtkLagrangeInterpolation::MaxDegree + 1) * (vtkLagrangeInterpolation::MaxDegree + 2) / 2];
......@@ -709,6 +774,89 @@ void vtkLagrangeInterpolation::WedgeShapeDerivatives(const int order[2], const d
int sn;
int numPts = numtripts * (tOrder + 1);
#ifdef VTK_21_POINT_WEDGE
if (order[3] == 21 && order[0] == 2)
{
const double r = pcoords[0];
const double s = pcoords[1];
const double t = 2 * pcoords[2] - 1.;
const double tm = t - 1.;
const double tp = t + 1.;
const double rsm = 1. - r - s;
const double rs = r * s;
// dN/dr
derivs[ 0] = 0.5*t*tm*(-3.0*rs + 2.0*r + 2.0*s + (3.0*s - 2.0)*rsm - 1.0);
derivs[ 1] = -0.5*t*tm*(3.0*rs - 4.0*r - 3.0*s*rsm + 1.0);
derivs[ 2] = -1.5*s*t*tm*(2*r + s - 1);
derivs[ 3] = 0.5*t*tp*(-3.0*rs + 2.0*r + 2.0*s + (3.0*s - 2.0)*rsm - 1.0);
derivs[ 4] = -0.5*t*tp*(3.0*rs - 4.0*r - 3.0*s*rsm + 1.0);
derivs[ 5] = -1.5*s*t*tp*(2*r + s - 1);
derivs[ 6] = 0.5*t*(12.0*s - 4.0)*tm*(2*r + s - 1);
derivs[ 7] = 0.5*s*t*tm*(24.0*r + 12.0*s - 8.0);
derivs[ 8] = s*t*tm*(12.0*r + 6.0*s - 8.0);
derivs[ 9] = 0.5*t*(12.0*s - 4.0)*tp*(2*r + s - 1);
derivs[10] = 0.5*s*t*tp*(24.0*r + 12.0*s - 8.0);
derivs[11] = s*t*tp*(12.0*r + 6.0*s - 8.0);
derivs[12] = tm*tp*(3.0*rs - 2.0*r - 2.0*s - (3.0*s - 2.0)*rsm + 1.0);
derivs[13] = tm*tp*(3.0*rs - 4.0*r - 3.0*s*rsm + 1.0);
derivs[14] = 3.0*s*tm*tp*(2*r + s - 1);
derivs[15] = 13.5*s*t*tm*(-2*r - s + 1);
derivs[16] = 13.5*s*t*tp*(-2*r - s + 1);
derivs[17] = (12.0*s - 4.0)*tm*tp*(-2*r - s + 1);
derivs[18] = -s*tm*tp*(24.0*r + 12.0*s - 8.0);
derivs[19] = s*tm*tp*(-24.0*r - 12.0*s + 16.0);
derivs[20] = 27.0*s*tm*tp*(2*r + s - 1);
// dN/ds
derivs[21] = 0.5*t*tm*(-3.0*rs + 2.0*r + 2.0*s + (3.0*r - 2.0)*rsm - 1.0);
derivs[22] = -1.5*r*t*tm*(r + 2*s - 1);
derivs[23] = -0.5*t*tm*(3.0*rs - 3.0*r*rsm - 4.0*s + 1.0);
derivs[24] = 0.5*t*tp*(-3.0*rs + 2.0*r + 2.0*s + (3.0*r - 2.0)*rsm - 1.0);
derivs[25] = -1.5*r*t*tp*(r + 2*s - 1);
derivs[26] = -0.5*t*tp*(3.0*rs - 3.0*r*rsm - 4.0*s + 1.0);
derivs[27] = r*t*tm*(6.0*r + 12.0*s - 8.0);
derivs[28] = 0.5*r*t*tm*(12.0*r + 24.0*s - 8.0);
derivs[29] = 0.5*t*(12.0*r - 4.0)*tm*(r + 2*s - 1);
derivs[30] = r*t*tp*(6.0*r + 12.0*s - 8.0);
derivs[31] = 0.5*r*t*tp*(12.0*r + 24.0*s - 8.0);
derivs[32] = 0.5*t*(12.0*r - 4.0)*tp*(r + 2*s - 1);
derivs[33] = tm*tp*(3.0*rs - 2.0*r - 2.0*s - (3.0*r - 2.0)*rsm + 1.0);
derivs[34] = 3.0*r*tm*tp*(r + 2*s - 1);
derivs[35] = tm*tp*(3.0*rs - 3.0*r*rsm - 4.0*s + 1.0);
derivs[36] = 13.5*r*t*tm*(-r - 2*s + 1);
derivs[37] = 13.5*r*t*tp*(-r - 2*s + 1);
derivs[38] = r*tm*tp*(-12.0*r - 24.0*s + 16.0);
derivs[39] = -r*tm*tp*(12.0*r + 24.0*s - 8.0);
derivs[40] = (12.0*r - 4.0)*tm*tp*(-r - 2*s + 1);
derivs[41] = 27.0*r*tm*tp*(r + 2*s - 1);
// dN/dt
derivs[42] = -0.5*(-2*t + 1)*rsm*(3.0*rs - 2.0*r - 2.0*s + 1.0);
derivs[43] = 0.5*r*(-2*t + 1)*(-2.0*r - 3.0*s*rsm + 1.0);
derivs[44] = 0.5*s*(-2*t + 1)*(-3.0*r*rsm - 2.0*s + 1.0);
derivs[45] = 0.5*(2*t + 1)*rsm*(3.0*rs - 2.0*r - 2.0*s + 1.0);
derivs[46] = -0.5*r*(2*t + 1)*(-2.0*r - 3.0*s*rsm + 1.0);
derivs[47] = -0.5*s*(2*t + 1)*(-3.0*r*rsm - 2.0*s + 1.0);
derivs[48] = -0.5*r*(12.0*s - 4.0)*(2*t - 1)*rsm;
derivs[49] = 0.5*rs*(2*t - 1)*(12.0*r + 12.0*s - 8.0);
derivs[50] = -0.5*s*(12.0*r - 4.0)*(2*t - 1)*rsm;
derivs[51] = -0.5*r*(12.0*s - 4.0)*(2*t + 1)*rsm;
derivs[52] = 0.5*rs*(2*t + 1)*(12.0*r + 12.0*s - 8.0);
derivs[53] = 0.5*s*(12.0*r - 4.0)*(2*t + 1)*rsm;
derivs[54] = -2*t*rsm*(3.0*rs - 2.0*r - 2.0*s + 1.0);
derivs[55] = 2*r*t*(-2.0*r + 3.0*s*rsm + 1.0);
derivs[56] = 2*s*t*(-3.0*r*rsm - 2.0*s + 1.0);
derivs[57] = -13.5*rs*(-2*t + 1)*rsm;
derivs[58] = 13.5*rs*(2*t + 1)*rsm;
derivs[59] = 2*r*t*(12.0*s - 4.0)*rsm;
derivs[60] = rs*t*(-24.0*r - 24.0*s + 16.0);
derivs[61] = 2*s*t*(12.0*r - 4.0)*rsm;
derivs[62] = -54.0*rs*t*rsm;
return;
}
#endif
vtkIdType ijk[3];
for (int kk = 0; kk <= tOrder; ++kk)
{
......@@ -732,6 +880,55 @@ void vtkLagrangeInterpolation::WedgeShapeDerivatives(const int order[2], const d
}
}
void vtkLagrangeInterpolation::WedgeEvaluate(
const int order[4],
const double* pcoords,
double* fieldVals,
int fieldDim,
double* fieldAtPCoords)
{
this->PrepareForOrder(order);
this->WedgeShapeFunctions(order, pcoords, &this->ShapeSpace[0]);
// Loop over components of the field:
for (int cc = 0; cc < fieldDim; ++cc)
{
fieldAtPCoords[cc] = 0.;
// Loop over shape functions (per-DOF values of the cell):
for (int pp = 0; pp < order[3]; ++pp)
{
fieldAtPCoords[cc] += this->ShapeSpace[pp] * fieldVals[fieldDim * pp + cc];
}
}