Output value of solution at midpoint of facet

I/O
1

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()
2

You can use function evaluation, as shown in: Bitbucket

3

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?

4

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.

5

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)
6

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.