Extracting a csv of the solution and its gradient at the dofs

I/O
1

Hi, I solved a pde and now have the solution variable uh, of the type

dolfinx.fem.function.Function.

I would like to extract the position of the dofs, the value of uh at the dofs and the value of grad(uh) at the dofs.
By looking at other question I came up with this code but I don’t know how to extract the value of grad(uh) (and I’m not sure it works properly for the dofs and uh).

outfile1 = open("dofs.csv", "w")

coords = V.tabulate_dof_coordinates()
for x in coords:
    print(x[0],x[1],file = outfile1)
    
##---------------------
outfile2 = open("solution.csv", "w")

sol = uh.x.array
for x in sol:
    print(x, file = outfile2)
    
##------------------------
outfile3 = open("grad_solution.csv", "w")


???

Does someone knows how to extract those values?

2

You should interpolate the gradient of uh into an appropriate function space, and then tabulate the dof coordinates of that space.

For arguments sake, say you have a P-th order Lagrange space. Then the gradient is in a P-1 discontinuous Lagrange space.

3

Is the interpolation computed in the following way?

W = fem.FunctionSpace(domain, ("DG", 1))
grad = fem.Function(W)
grad.interpolate(lambda x: ufl.grad(uh))

This makes the kernel of JupyterLab crash, I probably did something wrong.
Moreover, is W the right space? For uh I’m using

V = fem.FunctionSpace(domain, ("CG", 1))

Thank you really much!

4

I also found in the documentation of ufl the class ufl.differentiation.Grad, that should be the type of ufl.grad(uh). This class has the method evaluate(x, mapping, component, index_values, derivatives=()) that could be what I’m looking for, if I use it in this way:

outfile1 = open("dofs.csv", "w")
outfile3 = open("grad_solution.csv", "w")
coords = V.tabulate_dof_coordinates()
for x in coords:
    print(x[0],x[1],file = outfile1)
    ufl.grad(uh).evaluate(x, mapping, component, index_values, derivatives=())

Am I right?

I also didn’t understand the inputs of evaluate, and the documentation doesn’t explain further, is there some resource I can look up? :confused:

5

See, for example, the elasticity demo. Also keep in mind the purpose of ufl. It is a computational symbolic algebra library, and its purpose is not for numeric evaluation.

6

I tried to do as in the demo

W = fem.FunctionSpace(domain, ("DG", 0))
grad = ufl.grad(uh)
grad_expression = fem.Expression(grad, W.element.interpolation_points)
grad_h = fem.Function(W)
grad_h.interpolate(grad_expression)

but it gives me this error

RuntimeError: Function value size not equal to Expression value size

Have I used a wrong space? For uh I used CG1, so DG0 should be fine iirc

7

The gradient of a scalar function is a vector-valued function. i.e.

should be W = fem.VectorFunctionSpace(domain, ("DG", 0))

1 Like
8

I still have two questions, thanks in advance for the patience.

  1. I used the following code
outfile1 = open("dofs.csv", "w")

coords = V.tabulate_dof_coordinates().reshape((-1,gdim))
for x in coords:
    print(x[0],",",x[1],file = outfile1)

##---------------------
outfile2 = open("solution.csv", "w")

sol = uh.x.array
for x in sol:
    print(x, file = outfile2)

The dofs.csv file has 3596 rows, while the solution.csv file has 2137 rows. Shouldn’t they be equal? The first file should be the coordinates of the dofs and the second should be the value of uh at said points, right?

  1. I understood that once I have grad_h as above I should use grad_h.eval(…), but this requires to find out in which cell is contained each dof, is there a quick way to do it?

Thanks again to everyone

9

With only snippets of code it is really hard to help you any further. We would need a minimal code that can reproduce your issue.

1.If V is a first order Lagrange function, the number of vertices should be equal to the number of dofs.

  1. I do not see why you need to use eval after the interpolation into a vector DG space. You can tabulate the dof coordinates, where the [gdim*i:gdim*(i+1)] dofs corresponds to the ith entry in the DG_space.tabulate_dof_coordinates.

In general, i dont understand why you write the data to file if you have the mesh, dof coordinates and function space available.

10

Sure, as for the pde I’m trying to simulate, is the same I stated in this post.

My code is the following:

import numpy as np
import pyvista

from dolfinx.fem import (Constant,  Function, FunctionSpace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)

V = FunctionSpace(mesh, ("CG", 3))
import ufl
import dolfinx
from petsc4py.PETSc import ScalarType

####
n_temp = 1.4
reff= -1.44*n_temp**-2+0.71*n_temp**-1+0.668+0.0636*n_temp
temp = (1+reff)/(1-reff)
A = Constant(mesh, ScalarType(temp)) 
####

x = ufl.SpatialCoordinate(mesh)
u_ex = lambda x: x[0]**3*x[1]+x[0]**2*x[1]+2*x[1]**2+10*x[1]

mu_s = Constant(mesh, ScalarType(1.22))         
D = Constant(mesh, ScalarType(1/(3*1.22)))   
mu_a = fem.Function(V)
mu_a.interpolate(lambda x: 0.036 + 0.082*(np.square(x[0]-0.5)+np.square(x[1]-0.3)<=np.square(0.2)))

n = ufl.FacetNormal(mesh)
g = -ufl.dot(D*ufl.grad(u_ex(x)),n)
   

S = fem.Constant(domain, ScalarType(0.0))

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

F = D * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + mu_a * ufl.inner(u,v) * ufl.dx - ufl.inner(S, v) * ufl.dx

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = np.array(facet_tag.indices[facet_tag.values == marker])
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type


boundary_conditions = [BoundaryCondition("Dirichlet", 1, lambda x: np.zeros(x[0].size)),
                       BoundaryCondition("Dirichlet", 2, lambda x: np.zeros(x[0].size)),
                       BoundaryCondition("Robin", 3, (1/(2*A), Constant(mesh, 0.0))),]
                       #BoundaryCondition("Neumann", 4, g)]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.name = 'uh'

#Data extraction

W = fem.VectorFunctionSpace(domain, ("DG", 0))
grad = ufl.grad(uh)
grad_expression = fem.Expression(grad, W.element.interpolation_points)
grad_h = fem.Function(W)
grad_h.interpolate(grad_expression)

outfile1 = open("dofs.csv", "w")

coords = V.tabulate_dof_coordinates().reshape((-1,gdim))
for x in coords:
    print(x[0],",",x[1],file = outfile1)

##---------------------
outfile2 = open("solution.csv", "w")

sol = uh.x.array
for x in sol:
    print(x, file = outfile2)
    
##------------------------
outfile3 = open("grad_solution.csv", "w")

???
    
##------------------------
outfile4 = open("mu_a.csv", "w")

data_mua = mu_a.x.array

for x in data_mua:
    print(x, file = outfile4)

So what I’m trying to do is extract the dofs, the solution at the dofs, the gradient of the solution at the dofs and the constant mu_a at the dofs

11

So here V is a 3rd order space, which means that you have more degrees of freedom than you have vertices in the mesh.
Please note that for the different spaces V and W has a different number of degrees of freedom ( and you should) use W = fem.VectorFunctionSpace(domain, ("DG", degree_of_V -1))

I would strongly suggest making this a minimal working example.
The code you just posted is littered with typos, undefined variables etc.
I’ve cleaned up some of the code:

import ufl
import numpy as np

from dolfinx.fem import (Constant,  Function, FunctionSpace, VectorFunctionSpace,
                         dirichletbc, Expression, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (Measure, grad, inner, lhs, rhs)


mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
gdim = mesh.geometry.dim
V = FunctionSpace(mesh, ("CG", 3))

####
n_temp = 1.4
reff = -1.44*n_temp**-2+0.71*n_temp**-1+0.668+0.0636*n_temp
temp = (1+reff)/(1-reff)
A = Constant(mesh, ScalarType(temp))
####

x = ufl.SpatialCoordinate(mesh)
def u_ex(x): return x[0]**3*x[1]+x[0]**2*x[1]+2*x[1]**2+10*x[1]


mu_s = Constant(mesh, ScalarType(1.22))
D = Constant(mesh, ScalarType(1/(3*1.22)))
mu_a = Function(V)
mu_a.interpolate(lambda x: 0.036 + 0.082 *
                 (np.square(x[0]-0.5)+np.square(x[1]-0.3) <= np.square(0.2)))

n = ufl.FacetNormal(mesh)
g = -ufl.dot(D*ufl.grad(u_ex(x)), n)


S = Constant(mesh, ScalarType(0.0))

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

F = D * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + mu_a * \
    ufl.inner(u, v) * ufl.dx - ufl.inner(S, v) * ufl.dx

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(
    mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)


class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = np.array(facet_tag.indices[facet_tag.values == marker])
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
            self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v) * ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))

    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type


boundary_conditions = [BoundaryCondition("Dirichlet", 1, lambda x: np.zeros(x[0].size)),
                       BoundaryCondition(
                           "Dirichlet", 2, lambda x: np.zeros(x[0].size)),
                       BoundaryCondition("Robin", 3, (1/(2*A), Constant(mesh, 0.0))), ]
# BoundaryCondition("Neumann", 4, g)]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={
                        "ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.name = 'uh'

# Data extraction

W = VectorFunctionSpace(mesh, ("DG", 2))
grad = ufl.grad(uh)
grad_expression = Expression(grad, W.element.interpolation_points)
grad_h = Function(W)
grad_h.interpolate(grad_expression)

outfile1 = open("dofs.csv", "w")

coords = V.tabulate_dof_coordinates().reshape((-1, 3))
for x in coords:
    print(x[0], ",", x[1], file=outfile1)

# ---------------------
outfile2 = open("solution.csv", "w")

sol = uh.x.array
print(len(sol), coords.shape)

for x in sol:
    print(x, file=outfile2)

which now returns the same shape for both the function and the dof coordinates.

12

Sorry, I replied late at night and both missed a question and butchered my code trying to cut useless lines…
The question was:

This is because I need the data for a neural network.

You’re right about the fact that the shape is the same for the dofs and uh, but the csv files have different lenghts.

Using this code

import ufl
import numpy as np

from dolfinx.fem import (Constant,  Function, FunctionSpace, VectorFunctionSpace,
                         dirichletbc, Expression, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (Measure, grad, inner, lhs, rhs)


mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
gdim = mesh.geometry.dim
V = FunctionSpace(mesh, ("CG", 3))

####
n_temp = 1.4
reff = -1.44*n_temp**-2+0.71*n_temp**-1+0.668+0.0636*n_temp
temp = (1+reff)/(1-reff)
A = Constant(mesh, ScalarType(temp))
####

x = ufl.SpatialCoordinate(mesh)
def u_ex(x): return x[0]**3*x[1]+x[0]**2*x[1]+2*x[1]**2+10*x[1]


mu_s = Constant(mesh, ScalarType(1.22))
D = Constant(mesh, ScalarType(1/(3*1.22)))
mu_a = Function(V)
mu_a.interpolate(lambda x: 0.036 + 0.082 *
                 (np.square(x[0]-0.5)+np.square(x[1]-0.3) <= np.square(0.2)))

n = ufl.FacetNormal(mesh)
g = -ufl.dot(D*ufl.grad(u_ex(x)), n)


S = -ufl.div(D*ufl.grad(u_ex(x)))+mu_a*u_ex(x)

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

F = D * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + mu_a * \
    ufl.inner(u, v) * ufl.dx - ufl.inner(S, v) * ufl.dx

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(
    mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)


class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = np.array(facet_tag.indices[facet_tag.values == marker])
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
            self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v) * ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))

    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type


boundary_conditions = [BoundaryCondition("Dirichlet", 1, lambda x: np.zeros(x[0].size)),
                       BoundaryCondition(
                           "Dirichlet", 2, lambda x: np.zeros(x[0].size)),
                       BoundaryCondition("Robin", 3, (1/(2*A), Constant(mesh, 0.0))), ]
# BoundaryCondition("Neumann", 4, g)]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={
                        "ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.name = 'uh'

# Data extraction

W = VectorFunctionSpace(mesh, ("DG", 2))
grad = ufl.grad(uh)
grad_expression = Expression(grad, W.element.interpolation_points)
grad_h = Function(W)
grad_h.interpolate(grad_expression)

outfile1 = open("dofs.csv", "w")

coords = V.tabulate_dof_coordinates().reshape((-1, 3))
for x in coords:
    print(x[0], ",", x[1], file=outfile1)

# ---------------------
outfile2 = open("solution.csv", "w")

sol = uh.x.array
print(len(sol), coords.shape)

for x in sol:
    print(x, file=outfile2)

Both outputs should be of 961 lines, but dofs.csv contains 901 lines and solution.csv contains 848 lines. Am I missing something?

Also, if I print all three components of the dofs

outfile1 = open("dofs2.csv", "w")

coords = V.tabulate_dof_coordinates().reshape((-1, 3))
for x in coords:
    print(x[0], ",", x[1], ",", x[2], file=outfile1)

I get 771 lines, that’s odd right?

Lastly, I didn’t understand this:

In this way I would have grad(uh) at the DG_space dofs, but I would like to find the value of grad(uh) at the same points as uh, so at the CG_space dofs

13

I cannot reproduce this by running your code. I get 961 lines in each of the files.

The gradient is not well defined at the vertices of the mesh (or at the degrees of freedom of the uh-space), as the quantity is discontinuous (each element will have an entry for each vertex). Even if you interpolate into a DG space of same order, dofs on vertices, edges and facets will have unique entries for each cell.

14

That’s strange, maybe there’s some issue with docker or jupyter, I will investigate that.

So if we want the values of uh and grad(uh) on a grid would the best choice be the dofs of the DG-space or would a custom grid in the domain be better?
In any case I should use eval to evaluate uh and grad(uh) in points not contained in the dofs of, repectively, CG-space and DG-space, right?

15

The best thing would be to interpolate both functions into a function space containing both of them. For instance for a Q-th order Lagrange space and a Q-1 order Discontinuous Lagrange space you could interpolate them into a Q-th order Discontinuous Lagrange space.

If you interpolate them into a common space, there is no need for using eval for anything.

16

Ok, thank you really much!