# Basic queries from the fenics book, chapter 1, Set #4

1.) Is this the way to read the ChangeLog?

For example, in the following website

Deprecate subclassing of Expression in Python; new Python class UserExpression introduced for user overloads

Hence, if I use the fenics version 2019.2.0.dev0

I replace Expression with UserExpression

Am I correct to how to read the changelog?

2.a) Let’s say I a matrix A; I write in matlab syntax.
A = [2, 3, 4; 1, 3, 4, 0, 2, 4]

What is the difference between dolfin matrix A and the numpy converted matrix A?

In page 34. of fenics book, it is mentioned

We may convert the matrix and vector data to numpy arrays by calling the array() method as shown before.

2.b) How to print the A matrix and b vector in AU = b. What I read so far, the .vector() is deprecated, but somehow it shows in the matrix in a terminal. The vector b is not showing at all. How do I look for the proper syntax in a change log?

If I change b to b.get_local() and then print it, I can see the vector. I am confused by the difference among b, b.vector().get_local(), and b.get_local()

Here is a full code

“”"
FEniCS tutorial demo program:
Poisson equation with Dirichlet and Neumann conditions.
As dn2_p2D.py, but the linear system is explicitly formed and solved.

-Laplace(u) = f on the unit square.
u = 1 + 2y^2 on x=0.
u = 2 + 2y^2 on x=1.
-du/dn = g on y=0 and y=1.
u = 1 + x^2 + 2y^2, f = -6, g = -4y.
“”"

from dolfin import *
import numpy

# Create mesh and define function space

mesh = UnitSquareMesh(2, 1)
V = FunctionSpace(mesh, ‘Lagrange’, 1)

# Define Dirichlet conditions for x=0 boundary

u_L = Expression(‘1 + 2*x[1]*x[1]’, degree = 2)

class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0]) < tol

Gamma_0 = DirichletBC(V, u_L, LeftBoundary())

# Define Dirichlet conditions for x=1 boundary

u_R = Expression(‘2 + 2*x[1]*x[1]’, degree = 2)

class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0] - 1) < tol

Gamma_1 = DirichletBC(V, u_R, RightBoundary())

bcs = [Gamma_0, Gamma_1]

# Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
g = Expression('-4x[1]', degree = 2)
L = f
v
dx - gvds

# Assemble and solve linear system

A = assemble(a)
b = assemble(L)

if mesh.num_cells() < 16:
print(‘A = assemble(a); b = assemble(L)’)
print((‘A before incorporation of essential BC:\n’, A.array()))
print((‘b before incorporation of essential BC:\n’, b.array()))

for bc in bcs:
bc.apply(A, b)

if mesh.num_cells() < 16:
print((‘A after incorporation of essential BC:\n’, A.array()))
print((‘b after incorporation of essential BC:\n’, b.array()))

# (symmetric modification of boundary conditions)

A, b = assemble_system(a, L, bcs)

if mesh.num_cells() < 16:
print(‘\nA, b = assemble_system(a, L, bcs)’)
print((‘A after incorporation of essential BC:\n’, A.array()))
print((‘b after incorporation of essential BC:\n’, b.array()))

# Compute solution

u = Function(V)
U = u.vector()
solve(A, U, b)

#plot(u)

print((“”"
Solution of the Poisson problem -Laplace(u) = f,
with u = u0 on x=0,1 and -du/dn = g at y=0,1.
%s
“”" % mesh))

# Dump solution to the screen

u_nodal_values = u.vector()
u_array = u_nodal_values.array()
coor = mesh.coordinates()
for i in range(len(u_array)):
print((‘u(%8g,%8g) = %g’ % (coor[i][0], coor[i][1], u_array[i])))

# Exact solution:

u_exact = Expression(‘1 + x[0]x[0] + 2x[1]*x[1]’, degree = 2)

# Verification

u_e = interpolate(u_exact, V)
u_e_array = u_e.vector().array()
print((‘Max error:’, numpy.abs(u_e_array - u_array).max()))

# Compare numerical and exact solution

center = (0.5, 0.5)
print((‘numerical u at the center point:’, u(center)))
print((‘exact u at the center point:’, u_exact(center)))

#interactive()

2c.) In principle, I should be able to print and see the mesh with the followign command in the following section of the above code

print(“”"
Solution of the Poisson problem -Laplace(u) = f,
with u = u0 on x=0,1 and -du/dn = g at y=0,1.
%s
“”" % mesh)

but I get

<dolfin.cpp.generation.UnitSquareMesh object at 0x7f230c0dde50>

I am not sure what is going on?

A minimal working code is

``````from dolfin import *

mesh = IntervalMesh(5, -2, 0)
plot(mesh)
plt.show()
print('%.2f' % mesh)
``````

3.) How to plot the unknown in 1D problem? Like the familiary y = sin(x), problem, how do I plot here? I input the dimension of the problem via command line.

For 3D, I type (works great):
python3 d1_p2D_GKC.py 2 10 3 4

For 2D I type (works great):
python3 d1_p2D_GKC.py 2 10 3

For 1D I type (Does not plot):
python3 d1_p2D_GKC.py 2 10

Here is the minimal working code.

``````from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt

import sys

degree = int(sys.argv[1])
print(degree)
divisions = [int(arg) for arg in sys.argv[2:]]
d = len(divisions)
domain_type = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh]
mesh = domain_type[d-1](*divisions)
V = FunctionSpace(mesh, 'Lagrange', degree)

# Create mesh and define function space
#mesh = UnitCubeMesh(6, 4, 5)
#V = FunctionSpace(mesh, 'Lagrange', 1)

# Define boundary conditions
u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree = 2)

tol = 1E-14 # tolerance for coordiante comparisons
def Dirichlet_boundary0(x, on_boundary):
return on_boundary and abs(x[0]) < tol

def Dirichlet_boundary1(x, on_boundary):
return on_boundary and abs(x[0] - 1) < tol

bc0 = DirichletBC(V, Constant(0), Dirichlet_boundary0)
bc1 = DirichletBC(V, Constant(1), Dirichlet_boundary1)
bcs = [bc0, bc1]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-2.0)
L = f*v*dx

# Compute solution
#u = Function(V)
#solve(a == L, u, bcs)

u = Function(V)
solve(a == L, u, bcs, solver_parameters = {'linear_solver': 'cg', 'preconditioner': 'ilu'})

#File('poisson/parameters.xml') << parameters

from vedo.dolfin import plot as plot

plot(u)
#plt.show()
#plot(mesh)
#plt.show()

# Dump solution to file in VTK format
#file = File("poisson/solution.pvd")
#file << u
``````

4.) What does the ; mean in F(u;v)?

5.) As per previous discussion, I think I have correctly set the iterative solver, however, I get this error

File “vp2_np.py”, line 78, in
prm[“linear_solver”] = ‘gmres’

here is the minimal working code. I typed

python3 vp2_np.py a g 1 3 4

``````"""
FEniCS tutorial demo program:
Nonlinear Poisson equation with Dirichlet conditions
in x-direction and homogeneous Neumann (symmetry) conditions
in all other directions. The domain is the unit hypercube in
of a given dimension.

u = 0 at x=0, u=1 at x=1, du/dn=0 at all other boundaries.
q(u) = (1+u)^m

Solution method: automatic, i.e., by a NonlinearVariationalProblem/Solver
(Newton method), with automatic UFL computation of the derivative.
"""
from dolfin import *
import numpy, sys

# Usage:   ./vp2_np.py m|a |g|l degree nx ny nz
# Example: ./vp2_np.py m    l   1      3  4
J_comp = sys.argv[1]  # m (manual) or a (automatic) computation of J
answer = sys.argv[2]  # g (GMRES) or l (sparse LU) solver
iterative_solver = True if answer == 'g' else False

# Create mesh and define function space
degree = int(sys.argv[3])
divisions = [int(arg) for arg in sys.argv[4:]]
d = len(divisions)
domain_type = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh]
mesh = domain_type[d-1](*divisions)
V = FunctionSpace(mesh, 'Lagrange', degree)

# Define boundary conditions
tol = 1E-14
def left_boundary(x, on_boundary):
return on_boundary and abs(x[0]) < tol

def right_boundary(x, on_boundary):
return on_boundary and abs(x[0]-1) < tol

Gamma_0 = DirichletBC(V, Constant(0.0), left_boundary)
Gamma_1 = DirichletBC(V, Constant(1.0), right_boundary)
bcs = [Gamma_0, Gamma_1]

# Choice of nonlinear coefficient
m = 2

def q(u):
return (1+u)**m

def Dq(u):
return m*(1+u)**(m-1)

# Define variational problem
v  = TestFunction(V)
du = TrialFunction(V)
u_ = Function(V)  # most recently computed solution

# J must be a Jacobian (Gateaux derivative in direction of du)
if J_comp == 'm':
else:
J = derivative(F, u_, du)

# Compute solution
problem = NonlinearVariationalProblem(F, u_, bcs, J)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters
#info(prm, True)
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['maximum_iterations'] = 25
prm['newton_solver']['relaxation_parameter'] = 1.0
if iterative_solver:
prm['linear_solver'] = 'gmres'
prm['preconditioner'] = 'ilu'
prm['krylov_solver']['absolute_tolerance'] = 1E-9
prm['krylov_solver']['relative_tolerance'] = 1E-7
prm['krylov_solver']['maximum_iterations'] = 1000
prm['krylov_solver']['monitor_convergence'] = True
prm['krylov_solver']['nonzero_initial_guess'] = False
prm['krylov_solver']['gmres']['restart'] = 40
prm['krylov_solver']['preconditioner']['same_nonzero_pattern'] = True
prm['krylov_solver']['preconditioner']['ilu']['fill_level'] = 0
PROGRESS = 16
set_log_level(LogLevel.PROGRESS)

solver.solve()

print ("""
Solution of the nonlinear Poisson problem div(q(u)*grad(u)) = f,
with f=0, q(u) = (1+u)^m, u=0 at x=0 and u=1 at x=1.
%s
""" % mesh)

# Find max error
u_exact = Expression('pow((pow(2, m+1)-1)*x[0] + 1, 1.0/(m+1)) - 1', m=m, degree = 2)
u_e = interpolate(u_exact, V)
import numpy
diff = numpy.abs(u_e.vector().get_local() - u_.vector().get_local()).max()
print ('Max error:', diff)

from vedo.dolfin import plot
plot(u_)

# Define variational problem
#v  = TestFunction(V)
#du = TrialFunction(V)
#u = Function(V)  # the unknown
#J = derivative(F, u, du)
``````

6.) I came across this, while doing the problem from the fencis book. I have a minimum working example at the end of the problem. I don’t know how to install scitools following instructions here

https://scitools.org.uk/iris/docs/v1.9.2/installing.html

In brief, I

apt-get install python3-numpy
python3 setup.py install

ModuleNotFoundError: No module named ‘pyke’
error: command ‘/usr/bin/python3’ failed with exit status 1

sudo apt-get install pyke

Building dependency tree
E: Unable to locate package pyke

Minimum working example is

``````"""Temperature variations in the ground."""

from dolfin import *
import sys, numpy, time

# Usage:   sin_daD.py degree D   nx ny nz
# Example: sin_daD.py   1   1.5  4  40

# Create mesh and define function space
degree = int(sys.argv[1])
D = float(sys.argv[2])
W = D/2.0
divisions = [int(arg) for arg in sys.argv[3:]]
print('degree=%d, D=%g, W=%g, %s cells' % \
(degree, D, W, 'x'.join(sys.argv[3:])))

d = len(divisions)  # no of space dimensions
if d == 1:
mesh = IntervalMesh(divisions[0], -D, 0)
elif d == 2:
mesh = RectangleMesh(Point(-W/2, -D), Point(W/2, 0), divisions[0], divisions[1])
elif d == 3:
mesh = BoxMesh(Point(-W/2, -W/2), Point(-D, W/2), Point(W/2, 0),
divisions[0], divisions[1], divisions[2])
V = FunctionSpace(mesh, 'Lagrange', degree)

# Define boundary conditions

T_R = 0
T_A = 1.0
omega = 2*pi
# soil:
T_R = 10
T_A = 10
omega = 7.27E-5

T_0 = Expression('T_R + T_A*sin(omega*t)',
T_R=T_R, T_A=T_A, omega=omega, t=0.0, degree = 2)

def surface(x, on_boundary):
return on_boundary and abs(x[d-1]) < 1E-14

bc = DirichletBC(V, T_0, surface)

period = 2*pi/omega
t_stop = 5*period
#t_stop = period/4
dt = period/14
theta = 1

kappa_0 = 2.3  # KN/s (K=Kelvin, temp. unit), for soil
kappa_0 = 12.3  # KN/s (K=Kelvin, temp. unit)
kappa_1 = 1E+4
kappa_1 = kappa_0

kappa_str = {}
kappa_str[1] = 'x[0] > -D/2 && x[0] < -D/2 + D/4 ? kappa_1 : kappa_0'
kappa_str[2] = 'x[0] > -W/4 && x[0] < W/4 '\
'&& x[1] > -D/2 && x[1] < -D/2 + D/4 ? '\
'kappa_1 : kappa_0'
kappa_str[3] = 'x[0] > -W/4 && x[0] < W/4 '\
'x[1] > -W/4 && x[1] < W/4 '\
'&& x[2] > -D/2 && x[2] < -D/2 + D/4 ?'\
'kappa_1 : kappa_0'

# Alternative way of defining the kappa function
class Kappa(Expression):
def eval(self, value, x):
"""x: spatial point, value[0]: function value."""
d = len(x)  # no of space dimensions
material = 0  # 0: outside, 1: inside
if d == 1:
if -D/2. < x[d-1] < -D/2. + D/4.:
material = 1
elif d == 2:
if -D/2. < x[d-1] < -D/2. + D/4. and \
-W/4. < x[0] < W/4.:
material = 1
elif d == 3:
if -D/2. < x[d-1] < -D/2. + D/4. and \
-W/4. < x[0] < W/4. and -W/4. < x[1] < W/4.:
material = 1
value[0] = kappa_0 if material == 0 else kappa_1

# Physical parameters
kappa = Expression(kappa_str[d],
D=D, W=W, kappa_0=kappa_0, kappa_1=kappa_1, degree = 2)
#kappa = Kappa(V)

# soil:
rho = 1500
c = 1600

print('Thermal diffusivity:', kappa_0/rho/c)

# Define initial condition
T_1 = interpolate(Constant(T_R), V)

# Define variational problem
T = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = rho*c*T*v*dx + theta*dt*kappa*\
L = (rho*c*T_1*v + dt*f*v -

A = assemble(a)
b = None  # variable used for memory savings in assemble calls

# Compute solution
T = Function(V)   # unknown at the new time level

# hack for plotting (first surface must span all values (?))
dummy = Expression('T_R - T_A/2.0 + 2*T_A*(x[%g]+D)' % (d-1),
T_R=T_R, T_A=T_A, D=D, degree = 2)

# Make all plot commands inctive
import scitools.misc
plot = scitools.misc.DoNothing(silent=True)
# Need initial dummy plot
viz = plot(dummy, axes=True,
title='Temperature', wireframe=False)
viz.elevate(-65)
#time.sleep(1)
viz.update(T_1)

import scitools.BoxField
start_pt = [0]*d; start_pt[-1] = -D  # start pt for line plot
import scitools.easyviz as ev

def line_plot():
"""Make a line plot of T along the vertical direction."""
if T.ufl_element().degree() != 1:
T2 = interpolate(T, FunctionSpace(mesh, 'Lagrange', 1))
else:
T2 = T
T_box = scitools.BoxField.dolfin_function2BoxField(
T2, mesh, divisions, uniform_mesh=True)
#T_box = scitools.BoxField.update_from_dolfin_array(
#        T.vector().array(), T_box)
coor, Tval, fixed, snapped = \
T_box.gridline(start_pt, direction=d-1)

# Use just one ev.plot command, not hold('on') and two ev.plot
# etc for smooth movie on the screen
if kappa_0 == kappa_1:  # analytical solution?
ev.plot(coor, Tval, 'r-',
coor, T_exact(coor), 'b-',
axis=[-D, 0, T_R-T_A, T_R+T_A],
xlabel='depth', ylabel='temperature',
legend=['numerical', 'exact, const kappa=%g' % kappa_0],
legend_loc='upper left',
title='t=%.4f' % t)
else:
ev.plot(coor, Tval, 'r-',
axis=[-D, 0, T_R-T_A, T_R+T_A],
xlabel='depth', ylabel='temperature',
title='t=%.4f' % t)

ev.savefig('tmp_%04d.png' % counter)
time.sleep(0.1)

def T_exact(x):
a = sqrt(omega*rho*c/(2*kappa_0))
return T_R + T_A*exp(a*x)*numpy.sin(omega*t + a*x)

n = FacetNormal(mesh)  # for flux computation at the top boundary

t = dt
counter = 0
while t <= t_stop:
print('time =', t)
b = assemble(L, tensor=b)
T_0.t = t
help = interpolate(T_0, V)
print('T_0:', help.vector().array()[0])
bc.apply(A, b)
solve(A, T.vector(), b)
viz.update(T)
#viz.write_ps('diff2_%04d' % counter)
line_plot()

total_flux = assemble(flux)
print('Total flux:', total_flux)

t += dt
T_1.assign(T)
counter += 1

viz = plot(T, title='Final solution')
viz.elevate(-65)
viz.azimuth(-40)
viz.update(T)

#interactive()
``````

Hi,
I will try to answer as precise as I can.

2.a) Let’s say I a matrix A; I write in matlab syntax.
A = [2, 3, 4; 1, 3, 4, 0, 2, 4]

What is the difference between dolfin matrix A and the numpy converted matrix A?

First Matlab and numpy arrays are discrete arrays, while fenics arrays are continuous (can be evaluated outside of the vertices of the mesh).

2.b) How to print the A matrix and b vector in AU = b. What I read so far, the .vector() is deprecated, but somehow it shows in the matrix in a terminal. The vector b is not showing at all. How do I look for the proper syntax in a change log?

For matrix I don’t know how and why you would print it.
For Function you can do plot(u), or if you want the vector print(u.vector()[:]). Vectors already assembled they don’t have the same meaning as function, so I don’t see why you would print it…

2c.) In principle, I should be able to print and see the mesh with the followign command in the following section of the above code

To evaluate a function on a point you can do:

``````center = Point((0.5, 0.5))
print(u(center))
``````

When you do ` print(mesh)` since mesh is an object it will return its memory emplacement and not a good looking representation.
To represent the mesh just do `plot(mesh)`

3.) How to plot the unknown in 1D problem? Like the familiary y = sin(x), problem, how do I plot here? I input the dimension of the problem via command line.

First there is a mistake in your code. To define u0 you call x[1] which doesn’t exist in 1D space.

4.) What does the ; mean in F(u;v)?

The ; just separate the two variables of the bi-linear form (linear in u and linear in v)

File “vp2_np.py”, line 78, in
prm[“linear_solver”] = ‘gmres’

I am not sure what you are trying to do here… Please remember that the examples in the fenics book are not up-to-date. So perhaps this parameter “linear_solver” is not available anymore. But I find weird to try linear method on non-linear problem.

1 Like

In 5.) of my post, although the problem is non-linear, I don’t think the problem remains nonlinear after manipulation to Eq.(1.89) of fenics book, page 45.

I came to learn that the fenics book codes are outdated. But this is what I / we have at hand, and I am trying to work through it, while reading the book. I am trying my best to actively read, do numerical hands on computation.

The problem is nonlinear. It’s just that you are using Newton’s method to solver the resulting system of nonlinear equations. Newton’s method solves the system of equations incrementally, and therefore solving a linear system at each (newton) iteration. That is where you need your linear solver, and for more advanced problems solving this resulting linear system is highly nontrivial and needs fine tuning of the solver parameters.

The reference you point out isn’t accurate or maybe we are looking at different versions of the same tutorial. book and tutorial; you could add links for these for people to pull up instead of mentioning page and equation numbers.

As suggested before; you can use `solver.parameters.keys()` to start looking at the list of available options and suboptions.

IMHO the best way would be to take a look at the set of demos in the latest release of FEniCS and try to implement these yourself, rather than going over the complete FEniCS book. Another very useful resource for more advanced problems is @bleyerj’s COMET

2 Likes