Output solution across surface of mesh with element coordinates to a .txt file

I/O
1

If I have a very basic simulation like this:

from __future__ import print_function
from dolfin import *
from ufl import as_tensor
from ufl import Index

from mpi4py import MPI

mesh = UnitCubeMesh(16,16,16)
Space = FunctionSpace(mesh,'P', 1)
T = project(Expression("x[0]+0.1*x[1] -0.4*x[2]", degree=3), Space)

How can I output the values of T on each element on the surface, along with the relative coordinates of the element to a .txt file or similar. For example, something like:

x    y    z    T
0.0 0.0 0.0 100
0.3 0.3 0.2 102
.
.
.
.

Is this even possible? Thanks in advance

2

You can use the tabulate_dof_coordinates method of FunctionSpace, as in the following example:

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(1,1,1)
Space = FunctionSpace(mesh,'P', 1)
T = project(Expression("x[0]+0.1*x[1] -0.4*x[2]", degree=3), Space)

# Coordinates of DoFs in Space, as a matrix with a row for each point's
# x, y, and z coordinates:
x_Space = Space.tabulate_dof_coordinates().reshape((-1,mesh.geometry().dim()))

# Values of DoFs in Function T, as a single-column matrix:
T_values = T.vector().get_local().reshape((-1,1))

# Put together into a matrix of the desired form:
print(np.append(x_Space,T_values,1))
3 Likes
3

Ok great! Thank you! How would I extend this to output the coordinates of the FacetNormal for each element on the surface of the mesh?

4

When you say the coordinates of the FacetNormal, do you mean:

  • A tuple of coordinates describing the direction of the normal vector for a given facet?

If so, you should note a few things:
The FacetNormal is constant per facet, i.e. it would be equivalent to having an element with three vector DG functions, one located at each edge.

To get the physical facet normal (and not an approximation through projection or similar methods), you need to use a normalized covariant Piola transform
n_phys = J^{-T} n_ref / ||J^{-T} n_ref||
See for instance DOI: 10.1137/08073901X, where J is the jacobian of the cell, n_re$ the facet normal of the reference element.

6

Yes that is correct, my hope is to try and extract a physical normal vector for every element on the surface of the mesh, so I would need a tuple of coordinates describing the direction of the normal vector for a given facet on the surface of the mesh. My intent is to use this for a further emission simulation in Rust where the direction of emission is given by the normal vector at that specific element.

How would I then implement getting the normal vector at each point, and then doing this transformation that you described and then outputting this to a .txt file like what was given above for the previous answer?

7

It is quite complicated to compute the Jacobian in dolfin. It is quite alot easier in dolfinx

Here is the complicated dolfin version

from dolfin import *
import numpy as np

from FIAT.reference_element import UFCTetrahedron, UFCTriangle
from FIAT import ufc_simplex, Lagrange
mesh = UnitSquareMesh(10, 10)

ufc_element = UFCTetrahedron() if mesh.ufl_cell() == tetrahedron else UFCTriangle()

tdim = mesh.topology().dim()
fdim = tdim - 1
top = mesh.topology()
mesh.init(tdim, fdim)
mesh.init(fdim, tdim)
f_to_c = top(fdim, tdim)

bndry_facets = []
for facet in facets(mesh):
    f_index = facet.index()
    cells = f_to_c(f_index)
    if len(cells) == 1:
        bndry_facets.append(f_index)

# As we are using simplices, we can compute the Jacobian at an arbitrary point
fiat_el = Lagrange(ufc_simplex(tdim), 1)  
dphi_ = fiat_el.tabulate(1, np.array([0,]*tdim))
dphi = np.zeros((tdim, Cell(mesh, 0).num_entities(0)), dtype=np.float64)
dphi[0] = dphi_[(1,0)]
dphi[1] = dphi_[(0,1)]


c_to_f = top(tdim, fdim)
for facet in bndry_facets:
    cell = f_to_c(facet)[0]
    facets = c_to_f(cell)
    local_index = np.flatnonzero(facets==facet)[0]
    normal = ufc_element.compute_normal(local_index)
    normal_n = normal/np.linalg.norm(normal)
    coord_dofs = mesh.cells()[cell]
    coords = mesh.coordinates()[coord_dofs]
    J = np.dot(dphi, coords).T
    Jinv = np.linalg.inv(J)
    n_glob = np.dot(Jinv.T, normal_n)
    n_glob = n_glob/np.linalg.norm(n_glob)
    print(facet, n_glob)
8

OK great! I will try and expand this to 3D as the problem I am working with is in 3D. By the looks of things, I will have to expand the dphi terms more? Also out of curiosity, what would a version in dolphinx for a UnitCubeMesh look like? Could this be implemented effectively with the dolphin code for the output to a .txt file that I quoted in the original question?

9

This should automatically work for 3D meshes. I don not have time to sketch it for dolfinx now, i might have time tomorrow

10

After thinking about this, some minor modifications are needed in this line, and previous line to work in 3D, but it should be straightforward to change.

11

ok, I see. However I am still attempting to understand the code you have shown, since it seems quit alien to me as someone who is not a computer scientist nor particular gifted with Fenics.

fiat_el = Lagrange(ufc_simplex(tdim), 1)  
dphi_ = fiat_el.tabulate(1, np.array([0,]*tdim))
dphi = np.zeros((tdim, Cell(mesh, 0).num_entities(0)), dtype=np.float64)
dphi[0] = dphi_[(1,0)]
dphi[1] = dphi_[(0,1)]

For this section of code, what exactly are you doing here? What does each line represent? I am at a loss how I can expand this to a 3D mesh.

12

Here is a generalized code:

from dolfin import *
import numpy as np

from FIAT.reference_element import UFCTetrahedron, UFCTriangle
from FIAT import ufc_simplex, Lagrange
mesh = UnitCubeMesh(3, 3, 3)
# mesh = UnitSquareMesh(10, 10)

ufc_element = UFCTetrahedron() if mesh.ufl_cell() == tetrahedron else UFCTriangle()

tdim = mesh.topology().dim()
fdim = tdim - 1
top = mesh.topology()
mesh.init(tdim, fdim)
mesh.init(fdim, tdim)
f_to_c = top(fdim, tdim)

bndry_facets = []
for facet in facets(mesh):
    f_index = facet.index()
    cells = f_to_c(f_index)
    if len(cells) == 1:
        bndry_facets.append(f_index)

# As we are using simplices, we can compute the Jacobian at an arbitrary point
fiat_el = Lagrange(ufc_simplex(tdim), 1)  
dphi_ = fiat_el.tabulate(1, np.array([0,]*tdim))
dphi = np.zeros((tdim, Cell(mesh, 0).num_entities(0)), dtype=np.float64)
for i in range(tdim):
    t = tuple()
    for j in range(tdim):
        if i == j:
            t += (1, )
        else:
            t += (0, )
    dphi[i] = dphi_[t]


c_to_f = top(tdim, fdim)
for facet in bndry_facets:
    cell = f_to_c(facet)[0]
    facets = c_to_f(cell)
    local_index = np.flatnonzero(facets==facet)[0]
    normal = ufc_element.compute_normal(local_index)
    normal_n = normal/np.linalg.norm(normal)
    coord_dofs = mesh.cells()[cell]
    coords = mesh.coordinates()[coord_dofs]
    J = np.dot(dphi, coords).T
    Jinv = np.linalg.inv(J)
    n_glob = np.dot(Jinv.T, normal_n)
    n_glob = n_glob/np.linalg.norm(n_glob)
    print(facet, n_glob)

What the snipped that you are asking about does it creates a first order Lagrange element, representing the mesh. The next steps is to compute \frac{\partial \phi}{\partial x_i} such that we in turn can compute the Jacobian at a point in the simplex (for the covariant Piola transform).

13

Ok great! Thank you so much! I am noticing one issue however, if you have the mesh be UnitCubeMesh(1,1,1), I would expect six normal vectors:

[0,0,-1]
[0,0,1]
[0,-1,0]
[0,1,0]
[1,0,0]
[-1,0,0]

Corresponding to each face of the mesh. However, this code give me:

0 [ 0.  0. -1.]
1 [ 0. -1.  0.]
3 [ 0.  0. -1.]
4 [-1.  0.  0.]
7 [ 0. -1.  0.]
8 [-1.  0.  0.]
12 [-1.00000000e+00  4.14340874e-16  0.00000000e+00]
13 [-1.00000000e+00  0.00000000e+00  4.14340874e-16]
14 [ 4.14340874e-16 -1.00000000e+00  0.00000000e+00]
15 [ 0.00000000e+00 -1.00000000e+00  4.14340874e-16]
16 [ 4.14340874e-16  0.00000000e+00 -1.00000000e+00]
17 [ 0.00000000e+00  4.14340874e-16 -1.00000000e+00]
14

Some minor modifications to the code, replacing compute_normal with compute_reference_normal gives you the following:

from dolfin import *
import numpy as np

from FIAT.reference_element import UFCTetrahedron, UFCTriangle
from FIAT import ufc_simplex, Lagrange
mesh = UnitCubeMesh(1, 1, 1)
# mesh = UnitSquareMesh(10, 10)

ufc_element = UFCTetrahedron() if mesh.ufl_cell() == tetrahedron else UFCTriangle()

tdim = mesh.topology().dim()
fdim = tdim - 1
top = mesh.topology()
mesh.init(tdim, fdim)
mesh.init(fdim, tdim)
f_to_c = top(fdim, tdim)

bndry_facets = []
for facet in facets(mesh):
    f_index = facet.index()
    cells = f_to_c(f_index)
    if len(cells) == 1:
        bndry_facets.append(f_index)

# As we are using simplices, we can compute the Jacobian at an arbitrary point
fiat_el = Lagrange(ufc_simplex(tdim), 1)  
dphi_ = fiat_el.tabulate(1, np.array([0,]*tdim))
dphi = np.zeros((tdim, Cell(mesh, 0).num_entities(0)), dtype=np.float64)
for i in range(tdim):
    t = tuple()
    for j in range(tdim):
        if i == j:
            t += (1, )
        else:
            t += (0, )
    dphi[i] = dphi_[t]

c_to_f = top(tdim, fdim)
for facet in bndry_facets:
    cell = f_to_c(facet)[0]
    facets = c_to_f(cell)
    local_index = np.flatnonzero(facets==facet)[0]
    normal = ufc_element.compute_reference_normal(fdim, local_index)
    coord_dofs = mesh.cells()[cell]
    coords = mesh.coordinates()[coord_dofs]
    J = np.dot(dphi, coords).T
    Jinv = np.linalg.inv(J)
    n_glob = np.dot(Jinv.T, normal)
    n_glob = n_glob/np.linalg.norm(n_glob)
    print(facet, n_glob)

with the result:

0 [ 0.  0. -1.]
1 [ 0. -1.  0.]
3 [ 0.  0. -1.]
4 [-1.  0.  0.]
7 [ 0. -1.  0.]
8 [-1.  0.  0.]
12 [  1.00000000e+00  -4.44089210e-16   0.00000000e+00]
13 [  1.00000000e+00   0.00000000e+00  -4.44089210e-16]
14 [ -4.44089210e-16   1.00000000e+00   0.00000000e+00]
15 [  0.00000000e+00   1.00000000e+00  -4.44089210e-16]
16 [ -4.44089210e-16   0.00000000e+00   1.00000000e+00]
17 [  0.00000000e+00  -4.44089210e-16   1.00000000e+00]
1 Like
15

However this still gives you double the expected amount of unit vectors since I assume this code finds the normal vector for both faces of each facet, whereas I would like to get the normal vectors pointing outwards, away from the mesh. Can this be done?

16

I dont understand what you are meaning. For a 1x1x1 unit cube of tetrahedral elements, there are 12 exterior facets (two on each side of the cube), thus you get two normal vectors pointing in the outwards direction for each of the sides of the cube.

17

Ah I see, my mistake, thank you very much!

18

Actually, one last thing. I have tried stitching together your provided code and the code from @kamensky such that I can output the spatial coordinates of the facet, value of T at that facet and the normal vector components to a txt file, but it does not seem to be working and I am unsure if the spatial coordinate values, T values and normal vector components actually correspond to each other:

from dolfin import *
import numpy as np

from FIAT.reference_element import UFCTetrahedron, UFCTriangle
from FIAT import ufc_simplex, Lagrange
mesh = UnitCubeMesh(1, 1, 1)
#mesh = UnitSquareMesh(10, 10)
Space = FunctionSpace(mesh,'P', 1)
Temperature = project(Expression("x[0]+0.1*x[1] -0.4*x[2]", degree=3), Space)
x_Space = Space.tabulate_dof_coordinates().reshape((-1,mesh.geometry().dim()))
T_values = abs(Temperature.vector().get_local().reshape((-1,1)))

ufc_element = UFCTetrahedron() if mesh.ufl_cell() == tetrahedron else UFCTriangle()

tdim = mesh.topology().dim()
fdim = tdim - 1
top = mesh.topology()
mesh.init(tdim, fdim)
mesh.init(fdim, tdim)
f_to_c = top(fdim, tdim)

bndry_facets = []
for facet in facets(mesh):
    f_index = facet.index()
    cells = f_to_c(f_index)
    if len(cells) == 1:
        bndry_facets.append(f_index)

# As we are using simplices, we can compute the Jacobian at an arbitrary point
fiat_el = Lagrange(ufc_simplex(tdim), 1)  
dphi_ = fiat_el.tabulate(1, np.array([0,]*tdim))
dphi = np.zeros((tdim, Cell(mesh, 0).num_entities(0)), dtype=np.float64)
for i in range(tdim):
    t = tuple()
    for j in range(tdim):
        if i == j:
            t += (1, )
        else:
            t += (0, )
    dphi[i] = dphi_[t]


c_to_f = top(tdim, fdim)
n_values = []
for facet in bndry_facets:
    cell = f_to_c(facet)[0]
    facets = c_to_f(cell)
    local_index = np.flatnonzero(facets==facet)[0]
    normal = ufc_element.compute_reference_normal(fdim, local_index)
    normal_n = normal/np.linalg.norm(normal)
    coord_dofs = mesh.cells()[cell]
    coords = mesh.coordinates()[coord_dofs]
    J = np.dot(dphi, coords).T
    Jinv = np.linalg.inv(J)
    n_glob = np.dot(Jinv.T, normal_n)
    n_glob = n_glob/np.linalg.norm(n_glob)
    n_values.append(n_glob)
#    print(facet,coords, n_glob)    

filename = "output.txt"

data = [x_Space, T_values, n_values]
with open(filename, 'w') as f:
  np.savetxt(filename, data, header = 'x,y,z,T(K),nx,ny,nz')
19

Please note that a P-1 function space has its degrees of freedom at the vertices. Therefore there is not a clear one to one relationship between a degree of freedom at a vertex, and the facet (imagine the degree of freedom at a corner of the cube).
Each facet on the boundary has three vertices, which each are associated with a degree of freedom each. This means that there is no constant value T over a single facet (except when all dofs at the vertices have the same value).
Therefore, it is unclear to me what you want to achieve with stitching together these codes.

20

So my initial idea was to try and output a list of spatial coordinates on the mesh along with the corresponding normal vectors at that point, since I would like to use this txt file that I would generate in a separate piece of code, which is specifically looking at emission of particles from the 3D mesh.

The emission itself would follow a cosine distribution relative to the normal vector at that point and the overall strength of emission would be proportional to the value of T at this point. I would then try and integrate this distribution across the whole mesh (in 3D) and compute the emission flux distribution at a distance x away from my mesh.

So that was my reason for hoping this could be accomplished this way. Is there no way in which this desired data can be output?

21

As I said, there is no unique normal vector at the dof coordinates of a CG 1 space (which is at the mesh vertices).
A simple way of getting an approximation to such a normal at the vertices is by projecting the FacetNormal into a vector CG 1 function space.