# How to define a function defined on boundary?

Hello anyone,

I once again do a numerical experiment based on My previous question.

The problem is:
find u = u(\boldsymbol{x}, t) obeying

\begin{alignat}{2} \frac{\partial u}{\partial t} - \nabla \cdot (\kappa \nabla u) &= f, && (\boldsymbol{x}, t) \in \Omega \times [0,T],\\ -\kappa \frac{\partial u}{\partial \boldsymbol{n}} &= g, &\quad& (\boldsymbol{x}, t) \in \partial \Omega \times [0,T],\\ u(0) &= u_0, && \boldsymbol{x} \in \Omega, \end{alignat}

where \Omega is unit square, T =2, and the diffusion coefficient

$$\kappa = \kappa (t).$$

The sources

\begin{align} f &= 2t - \kappa (6x + 12y),\\ g &= \begin{cases} 0, & y=0,\\ -6 \kappa, & y=1,\\ 0, & x=0,\\ -3 \kappa. & x=1, \end{cases}\\ u_0 &= 1 + x^3 + 2y^3 \end{align}

are derived from exact solution

$$u_{ex} = 1 + x^3 + 2y^3 + t^2.$$

The diffusion coefficient

$$\kappa = \sqrt{t}$$

in this demo.

I follow this tutorial to code as:

from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_matrix, create_vector, set_bc
import numpy as np
import time

start = time.time()

t = 0.0
T = 2.0
num_steps = 20
dt = (T - t) / num_steps

nx, ny = 5, 5
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))

class exact_solution():
def __init__(self, t):
self.t = t

def __call__(self, x):
return 1 + x[0] ** 3 + 2 * x[1] ** 3 + self.t ** 2

u_exact = exact_solution(t)

u_n = fem.Function(V)
u_n.interpolate(u_exact)

kappa_expression = np.sqrt(t)
kappa = fem.Constant(domain, kappa_expression)

class source_term():
def __init__(self, t, kappa):
self.t = t
self.kappa = kappa

def __call__(self, x):
return 2 * self.t - self.kappa * (6 * x[0] + 12 * x[1])

f_expression = source_term(t, kappa_expression)

f = fem.Function(V)
f.interpolate(f_expression)

g = fem.Function(V)

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

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
if locator == 1:
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
elif locator == 2:
g.x.array[facets] = np.full_like(facets, -6.0 * kappa_expression, dtype = default_scalar_type)
elif locator == 3:
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
elif locator == 4:
g.x.array[facets] = np.full_like(facets, -3.0 * kappa_expression, dtype = default_scalar_type)
else:
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))

facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.Measure("ds", domain = domain, subdomain_data = facet_tag)

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = u * v * ufl.dx + dt * kappa * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - u_n * v * ufl.dx - dt * f * v * ufl.dx + dt * g * v * ds(1) + dt * g * v * ds(2) + dt * g * v * ds(3) + dt * g * v * ds(4)
a = fem.form(ufl.lhs(F))
L = fem.form(ufl.rhs(F))

A = create_matrix(a)
b = create_vector(L)
uh = fem.Function(V)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

errors_L2 = []
for n in range(num_steps):

t += dt
kappa_expression = np.sqrt(t)
kappa.value = kappa_expression
u_exact.t = t
f_expression.t = t
f_expression.kappa = kappa_expression
f.interpolate(f_expression)

A.zeroEntries()
assemble_matrix(A, a)
A.assemble()

with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [[]])
set_bc(b, [])

solver.solve(b, uh.vector)
uh.x.scatter_forward()

u_n.x.array[:] = uh.x.array

V_ex = fem.functionspace(domain, ("Lagrange", 2))
u_ex = fem.Function(V_ex)
u_ex.interpolate(u_exact)
error_L2 = domain.comm.allreduce(fem.assemble_scalar(fem.form((uh - u_ex)**2 * ufl.dx)), op = MPI.SUM)
errors_L2.append(error_L2)

L2_error = np.sqrt(np.trapz(errors_L2, dx = dt))
print(f"L2-error: {L2_error:.2e}")

end = time.time()
print(f"Execution time: {(end - start) * 10 ** 3:.6f} ms")


then the system tells me

Traceback (most recent call last):
File "/media/math/个人文件空间/WangRan/work-fenics/SIR/demo4/ex.py", line 69, in <module>
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
IndexError: index 37 is out of bounds for axis 0 with size 36


I guess the error results from the definition of g since it is defined on the square, but how to define a function only on facet?
Moreover, is there any simple way to define the function g?
Thanks very much.

g.x.array is indexed over the degrees of freedom, not the facets. Once you have facets, find the corresponding dofs using dolfinx.fem — DOLFINx 0.9.0.0 documentation

Thanks. Now I define g as

g = fem.Function(V)
domain.topology.create_connectivity(fdim, domain.topology.dim)
bottom_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(1))
top_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(2))
left_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(3))
right_dofs = fem.locate_dofs_topological(V, fdim, facet_tag.find(4))
g.x.array[bottom_dofs] = np.full_like(bottom_dofs, 0.0, dtype = default_scalar_type)
g.x.array[top_dofs] = np.full_like(top_dofs, -6.0 * kappa_expression, dtype = default_scalar_type)
g.x.array[left_dofs] = np.full_like(left_dofs, 0.0, dtype = default_scalar_type)
g.x.array[right_dofs] = np.full_like(right_dofs, -3.0 * kappa_expression, dtype = default_scalar_type)


Also, I update g.x.array in the iteration, and I obtain the error

L2-error: 6.86e-01


I don’t know why the error is so large.

You should carry out a convergence analysis, i.e. progressively double nx and ny and see if the FE solution converges at the expected rate.

Thank you for your advice. I change the nx and ny from 2^1 to 2^7, then obtain the errors:

L2-error: 1.81e+00
L2-error: 8.83e-01
L2-error: 3.83e-01
L2-error: 1.28e-01
L2-error: 5.67e-02
L2-error: 1.01e-01
L2-error: 1.32e-01


I think these results show that my code is wrong…
However, I haven’t found the error.

The suggestion would be to start from/go back to a stationary problem, and not progress to an unsteady one until your steady code has the right convergence rates.