# Output value of solution at midpoint of facet

I currently have a UnitCubeMesh and I wish to output the values of my solution Temperature at the midpoints of each facet on a boundary, rather than on the nodes. Currently I have a method to project Temperature over my mesh and output the midpoints of the facets on my boundary, but I am unsure how to extend this to include values of Temperature at the midpoints of these facets:

How can this be done?

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(10,10,10)

class marker(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 1) and on_boundary

Space = FunctionSpace(mesh,'P', 1)
Temperature = project(Expression("x[0]+0.1*x[1]+ x[2]", degree=2), Space) # Temperature solution across mesh
V = FunctionSpace(mesh,'P', 1)

tdim = mesh.topology().dim()
fdim = tdim - 1
mf = MeshFunction("size_t", mesh, fdim, 0)
value = 1
marker().mark(mf, value)
facet_indices = np.flatnonzero(mf.array() == value)
surface_facets = np.array(list(facets(mesh)))[facet_indices]

mesh.init(fdim, tdim)
f_to_c = mesh.topology()(fdim, tdim)
c_to_f = mesh.topology()(tdim, fdim)

for facet in surface_facets:
midpoint = facet.midpoint().array()

You can use function evaluation, as shown in: Bitbucket

For the suggested function:

def test_near_evaluations(R, mesh):
# Test that we allow point evaluation that are slightly outside
parameters["allow_extrapolation"] = False

u0 = Function(R)
u0.vector()[:] = 1.0
a = Vertex(mesh, 0).point()
offset = 0.99*DOLFIN_EPS

a_shift_x = Point(a[0] - offset, a[1], a[2])
assert round(u0(a) - u0(a_shift_x), 7) == 0

a_shift_xyz = Point(a[0] - offset / sqrt(3),
a[1] - offset / sqrt(3),
a[2] - offset / sqrt(3))
assert round(u0(a) - u0(a_shift_xyz), 7) == 0

I am unsure why allow_extrapolation must be set to False, I would have thought this entire code relies on extrapolation to some degree?

What I meant is that you can use u(point). But you should allow extrwpolation as the computed midpoint might be outside the cell by Machine precision.

Currently I have this:

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(10,10,10)

class marker(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 1) and on_boundary

Space = FunctionSpace(mesh,'P', 1)
T = project(Expression("x[0]+0.1*x[1]+ x[2]", degree=2), Space) # Temperature solution across mesh
V = FunctionSpace(mesh,'P', 1)

tdim = mesh.topology().dim()
fdim = tdim - 1
mf = MeshFunction("size_t", mesh, fdim, 0)
value = 1
marker().mark(mf, value)
facet_indices = np.flatnonzero(mf.array() == value)
surface_facets = np.array(list(facets(mesh)))[facet_indices]

mesh.init(fdim, tdim)
f_to_c = mesh.topology()(fdim, tdim)
c_to_f = mesh.topology()(tdim, fdim)

coordinates = []
for facet in surface_facets:
midpoint = facet.midpoint().array()
coordinates.append(midpoint)

T.vector()[:] = T
T_values_mid = []
for i in np.arange(0,len(coordinates)+1, 1):
T_values_mid.append(T(coordinates[i].point()))

print(T_values_mid)

Following the u(point) method. Instead I currently get this error:

Traceback (most recent call last):
File "tempoutput.py", line 31, in <module>
T.vector()[:] = T
TypeError: __setitem__(): incompatible function arguments. The following argument types are supported:
1. (self: dolfin.cpp.la.GenericVector, arg0: slice, arg1: float) -> None
2. (self: dolfin.cpp.la.GenericVector, arg0: slice, arg1: dolfin.cpp.la.GenericVector) -> None
3. (self: dolfin.cpp.la.GenericVector, arg0: slice, arg1: numpy.ndarray[float64]) -> None
4. (self: dolfin.cpp.la.GenericVector, arg0: numpy.ndarray[int32], arg1: float) -> None
5. (self: dolfin.cpp.la.GenericVector, arg0: numpy.ndarray[int32], arg1: numpy.ndarray[float64]) -> None

Invoked with: <dolfin.cpp.la.PETScVector object at 0x7fdeb2155bf8>, slice(None, None, None), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Lagrange', tetrahedron, 1)), 12)

Im not sure what you are trying to do here, as you aretrying to assign the dolfin.Function to the array of the same function.
I think you should just remove this line.