Sampling the Solution over a Line
Let us revisit our simple 3d example from the introduction
[1]:
from matplotlib import pyplot
import numpy
try:
import pygmsh
with pygmsh.geo.Geometry() as geom:
poly = geom.add_polygon([
[ 0.0, 0.5, 0.0], [-0.1, 0.1, 0.0], [-0.5, 0.0, 0.0],
[-0.1, -0.1, 0.0], [ 0.0, -0.5, 0.0], [ 0.1, -0.1, 0.0],
[ 0.5, 0.0, 0.0], [ 0.1, 0.1, 0.0] ], mesh_size=0.05)
geom.twist(
poly,
translation_axis=[0, 0, 1], rotation_axis=[0, 0, 1], point_on_axis=[0, 0, 0],
angle=numpy.pi / 3,
)
mesh = geom.generate_mesh(verbose=False)
points, cells = mesh.points, mesh.cells_dict
domain3d = {"vertices":points.astype("float"), "simplices":cells["tetra"]}
except ImportError: # pygmsh not installed - use a simple cartesian domain
print("pygmsh module not found using a simple Cartesian domain - ignored")
from dune.grid import cartesianDomain
domain3d = cartesianDomain([-0.25,-0.25,0],[0.25,0.25,1],[30,30,60])
from dune.alugrid import aluSimplexGrid as leafGridView3d
gridView3d = leafGridView3d(domain3d)
As before we solve a simple Laplace problem
[2]:
from dune.fem.space import lagrange as solutionSpace
from dune.fem.scheme import galerkin as solutionScheme
from ufl import TrialFunction, TestFunction, SpatialCoordinate, dot, grad, dx, conditional, sqrt
space3d = solutionSpace(gridView3d, order=1)
u = TrialFunction(space3d)
v = TestFunction(space3d)
x = SpatialCoordinate(space3d)
scheme3d = solutionScheme((dot(grad(u),grad(v))+u*v)*dx ==
conditional(dot(x,x)<.01,100,0)*v*dx,
solver='cg')
uh3d = space3d.interpolate(0,name="solution")
info = scheme3d.solve(target=uh3d)
Instead of plotting this using paraview we want to only study the solution along a single line. This requires findings points \(x_i = x_0+\frac{i}{N}(x1-x0)\) for \(i=0,\dots,N\) within the unstructured grid. This would be expensive to compute on the Python so we implement this algorithm in C++ using the LineSegmentSampler
class available in Dune-Fem
. The resulting algorithm
returns a pair of two lists with coordinates \(x_i\) and the values of the grid function at these
points:
#include <vector>
#include <utility>
#include <dune/fem/misc/linesegmentsampler.hh>
template <class GF, class DT>
std::pair<std::vector<DT>, std::vector<typename GF::RangeType>>
sample(const GF &gf, DT &start, DT &end, int n)
{
Dune::Fem::LineSegmentSampler<typename GF::GridPartType> sampler(gf.gridPart(),start,end);
std::vector<DT> coords(n);
std::vector<typename GF::RangeType> values(n);
sampler(gf,values);
sampler.samplePoints(coords);
return std::make_pair(coords,values);
}
[3]:
import dune.generator.algorithm as algorithm
from dune.common import FieldVector
x0, x1 = FieldVector([0,0,0]), FieldVector([0,0,1])
p,v = algorithm.run('sample', 'utility.hh', uh3d, x0, x1, 100)
x,y = numpy.zeros(len(p)), numpy.zeros(len(p))
length = (x1-x0).two_norm
for i in range(len(x)):
x[i] = (p[i]-x0).two_norm / length
y[i] = v[i][0]
pyplot.plot(x,y)
[3]:
[<matplotlib.lines.Line2D at 0x732b5fb84980>]
Note: the coordinates returned are always in the interval \([0,1]\) so if physical coordinates are required, they need to be rescaled. Also, function values returned by the sample
function are always of a FieldVector
type, so that even for a scalar example a v[i]
is a vector of dimension one, so that y[i]=v[i][0]
has to be used.
A mentioned above any grid function can be passed in as argument to the sample
function. So for example plotting \(|\nabla u_h|\) is straight forward using the corresponding ufl expression. Since in this case automatic conversion from the ufl expression (available for example in the plotting function) to a grid function, we need to do this explicitly:
[4]:
from dune.fem.function import gridFunction
absGrad = gridFunction(sqrt(dot(grad(uh3d),grad(uh3d))))
p,v = algorithm.run('sample', 'utility.hh', absGrad, x0, x1, 100)
for i in range(len(x)):
y[i] = v[i][0]
pyplot.plot(x,y)
[4]:
[<matplotlib.lines.Line2D at 0x732b32729520>]
Similar we can plot both partial derivatives of the solution over the given line:
[5]:
absGrad = gridFunction(grad(uh3d))
p,v = algorithm.run('sample', 'utility.hh', absGrad, x0, x1, 100)
dx,dy = numpy.zeros(len(p)), numpy.zeros(len(p))
for i in range(len(x)):
dx[i] = v[i][0]
dy[i] = v[i][1]
pyplot.plot(x,dx)
pyplot.plot(x,dy)
[5]:
[<matplotlib.lines.Line2D at 0x732b30901bb0>]
Remark: plotting over a line has been included as a utility function so can easily be called as shown below
[6]:
from dune.fem.utility import lineSample
x,y = lineSample(uh3d,[0,0,0],[0,0,1],100)
pyplot.plot(x,y)
[6]:
[<matplotlib.lines.Line2D at 0x732b30961760>]
The functions takes the two endpoints of the line and the number of (equidistant) sample points as arguments. By default the method throws an exception if part of the line is outside of the computational domain. This behavior can be disabled by setting the additional optional :code:bool
parameter :code:allowOutside
to :code:True
.
Similar functions are available to obtain a :code:pointSample
and to sample a grid function along the boundary (:code:boundarySample
).