Interpolate from P0 to P1
In dune-tectonic
, in a branch that I haven't pushed yet, I have the following few lines of code:
ScalarVector surroundingArea(vertexBasis.size());
surroundingArea = 0.0;
Vector nodalTractionAverage(vertexBasis.size());
nodalTractionAverage = 0.0;
auto const &indexSet = gridView.indexSet();
for (auto const &intersection : frictionalBoundary) {
auto const &inside = intersection.inside();
auto const &refElement =
Dune::ReferenceElements<double, dimension>::general(inside.type());
auto const globalElementIndex = indexSet.index(inside);
auto const localFaceIndex = intersection.indexInInside();
auto const &face = inside.template subEntity<1>(localFaceIndex);
auto const area = face.geometry().volume();
int const numVertices = refElement.size(localFaceIndex, 1, dimension);
for (int i = 0; i < numVertices; i++) {
auto const localVertexIndex =
refElement.subEntity(localFaceIndex, 1, i, dimension);
auto const globalVertexIndex =
indexSet.subIndex(inside, localVertexIndex, dimension);
surroundingArea[globalVertexIndex] += area;
nodalTractionAverage[globalVertexIndex] += traction[globalElementIndex];
}
}
for (size_t i = 0; i < vertexBasis.size(); ++i)
if (frictionalBoundary.containsVertex(i))
nodalTractionAverage[i] /= surroundingArea[i];
What this does is: It takes P0 data from a vector traction
, which lives on the BoundaryPatch frictionalBoundary
only, and turns it into P1 data (the vector nodalTractionAverage
). To do so, for each vertex, it computes the average over all neighbouring faces (that belong to the BoundaryPatch).
On the one hand, I think that this kind of code is more widely applicable and could be moved to e.g. dune-fufem (or another module). On the other hand, it seems quite specific to me in the sense that
- it works only with the (ordered!) pair (P0,P1) and
- it works for a BoundaryPatch only.
I'm interested in learning all about the ifs, wheres, and hows.