Update of old Fenics Code for Poisson Problem

variational formulation
1

I have got the following old Python code which uses Fenics to solve the Poisson problem. I used to run it some years ago but it does not run with my current system.

I am basically following “User-defined expression by subclassing” as described on this page:
https://fenicsproject.org/docs/dolfin/1.5.0/python/programmers-reference/functions/expression/Expression.html

The error message that I get is an uncaught exception:
RuntimeError: Must supply C++ code to Expression. You may want to use UserExpression

I wonder whether anybody can confirm that the code does not run and how to fix the issue.

The code follows:

#!/usr/bin/python3

from fenics import *
from dolfin import *
import numpy as np

mesh = UnitSquareMesh(4, 4)

initial_refinements = 0                         # number of initial refinements before we start 

polydegree = 3                                  # polynomial degree of finite element space 

expressiondegree = 10                           # combined polynomial degree 

number_of_iterations = 7                        # number of iterations for main loop 

for i in range( 1, initial_refinements ):
    mesh = refine( mesh )

old_error_L2  = 1
old_error_H1  = 1
old_error_H10 = 1
old_error_max = 1

class Expression_u_D(Expression):
    def eval(self, value, x):
        value[0] = exp( - x[0]*x[0] - x[1]*x[1] )

class Expression_f(Expression):
    def eval(self, value, x):
        value[0] = 4 * exp( - x[0]*x[0] - x[1]*x[1] ) * ( 1 - x[0]*x[0] - x[1]*x[1] )





for i in range(0,number_of_iterations):

    print( "Iteration: %d" % i )
    print( "Number of cells: %d" % mesh.num_cells() )

    V = FunctionSpace(mesh, 'Lagrange', polydegree)

    u_D = Expression_u_D( degree = expressiondegree )

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    u  = TrialFunction(V)
    v  = TestFunction(V)

    f = Expression_f( degree = expressiondegree )

    a = dot( grad(u), grad(v) ) * dx
    L = f * v * dx

    parms = parameters["krylov_solver"]
    parms["relative_tolerance"] = 1.e-10
    parms["absolute_tolerance"] = 1.e-12
        
    u = Function(V)
    solve(a == L, u, bc, solver_parameters={"linear_solver" : "cg", "preconditioner" : "amg", "krylov_solver": parms} )

    plot(u)
    plot(mesh)

    vtkfile = File('poisson/solution.pvd')
    vtkfile << u

    error_L2  = errornorm( u_D, u, 'L2',  degree_rise=7, mesh=mesh )
    error_H1  = errornorm( u_D, u, 'H1',  degree_rise=7, mesh=mesh )
    error_H10 = errornorm( u_D, u, 'H10', degree_rise=7, mesh=mesh )

    vertex_values_u_D = u_D.compute_vertex_values(mesh)
    vertex_values_u = u.compute_vertex_values(mesh)
    import numpy as np
    error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

    print('error_L2  = {:e}\t : {:f}'.format( error_L2,  np.log2(old_error_L2 /error_L2 ) ) )
    print('error_H1  = {:e}\t : {:f}'.format( error_H1,  np.log2(old_error_H1 /error_H1 ) ) )
    print('error_H10 = {:e}\t : {:f}'.format( error_H10, np.log2(old_error_H10/error_H10) ) )
    print('error_max = {:e}\t : {:f}'.format( error_max, np.log2(old_error_max/error_max) ) )
    old_error_L2  = error_L2
    old_error_H1  = error_H1
    old_error_H10 = error_H10
    old_error_max = error_max

    mesh = refine( mesh )

# Hold plot
# interactive()
2

In this case, it’s as simple as just updating your Expression subclasses to instead be subclasses of UserExpression, e.g.,

class Expression_u_D(UserExpression):
    def eval(self, value, x):
        value[0] = exp( - x[0]*x[0] - x[1]*x[1] )
    def value_shape(self):
        return ()

(The value_shape method is not strictly necessary for scalar expressions, because that is the default, but it gets rid of a warning message.)

1 Like