2D poisson problem with piecewise function f

1

Hello everyone. I am quite new to coding and this problem has been bugging me for a long time.

I am trying to code a boundary value problem in the unit square mesh, however, I have a problem defining the function properly.

from __future__ import print_function

"""
Created on Mon Dec 20 21:30:34 2021

@author: Shirley
"""

"""
FEniCS tutorial demo program: Poisson equation with Dirichlet conditions.
Test problem is chosen to give an exact solution at all nodes of the mesh.
  -Laplace(u) = f    in the unit square
            u = u_D  on the boundary
  u_D = 0
         ( 1 \in \Omega'  
    f = {  
         ( 0 \in  \Omega - \Omega'
"""

from fenics import *
import matplotlib.pyplot as plt

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

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

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

tol = 1E-14

class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return ( (x[0] >= 0.35 + tol) &  (x[0] <= 0.7 + tol) & (x[1] >= 0.35 + tol) & (x[1] <= 0.7 + tol) )

class Omega_1_bdry_left(SubDomain):
    def inside(self, x, on_boundary):
        return ( (x[0] >= 0.34 + tol) &  (x[0] < 0.35 + tol) & (x[1] >= 0.34 + tol) & (x[1] <= 0.71 + tol) )

class Omega_1_bdry_right(SubDomain):
    def inside(self, x, on_boundary):
        return ( (x[0] >= 0.7 + tol) &  (x[0] < 0.71 + tol) & (x[1] >= 0.34 + tol) & (x[1] <= 0.71 + tol) )

class Omega_1_bdry_bottom(SubDomain):
    def inside(self, x, on_boundary):
        return ( (x[0] >= 0.35 + tol) &  (x[0] < 0.7 + tol) & (x[1] >= 0.34 + tol) & (x[1] <= 0.35 + tol) )

class Omega_1_bdry_top(SubDomain):
    def inside(self, x, on_boundary):
        return ( (x[0] >= 0.35 + tol) &  (x[0] < 0.7 + tol) & (x[1] >= 0.7 + tol) & (x[1] <= 0.71 + tol) )

# class Omega_2(SubDomain):
#     def inside(self, x, on_boundary):
#         return x[1] >= 0.5 - tol

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

p1 = Constant(1.0)
p2 = Constant(0.0)

f = Expression(' ( (x[0] >= 0.35 + DOLFIN_EPS)&&(x[0] <= 0.7 + DOLFIN_EPS)&&(x[1] >= 0.35 + DOLFIN_EPS)&&(x[1] <= 0.7 + DOLFIN_EPS) ) ? p1 : ( ( (x[0] < 0.35 + DOLFIN_EPS)&&(x[0] > 0.7 + DOLFIN_EPS)&&(x[1] < 0.35 + DOLFIN_EPS)&&(x[1] > 0.7 + DOLFIN_EPS) ) ? p2 : 0 ) ) ',
             p1=p1, p2=p2, degree=2 )

# f = Constant(-6.0)

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

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

# Plot solution and mesh
plot(u)
plot(mesh)

# Save solution to file in VTK format
vtkfile = File('poisson/solution.pvd')
vtkfile << u

# Compute error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')

# Compute maximum error at vertices
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 errors
print('error_L2  =', error_L2)
print('error_max =', error_max)

# Hold plot
plt.show()

I get this error, but I don’t know why:

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File “/Users/Shirley/Downloads/PrgrmngMtrs/fiat-master/test/plsworktest.py”, line 71, in
f = Expression(’ ( (x[0] >= 0.35 + DOLFIN_EPS)&&(x[0] <= 0.7 + DOLFIN_EPS)&&(x[1] >= 0.35 + DOLFIN_EPS)&&(x[1] <= 0.7 + DOLFIN_EPS) ) ? p1 : ( ( (x[0] < 0.35 + DOLFIN_EPS)&&(x[0] > 0.7 + DOLFIN_EPS)&&(x[1] < 0.35 + DOLFIN_EPS)&&(x[1] > 0.7 + DOLFIN_EPS) ) ? p2 : 0 ) ) ',
File “/Users/Shirley/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/function/expression.py”, line 400, in init
self._cpp_object = jit.compile_expression(cpp_code, params)
File “/Users/Shirley/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/function/jit.py”, line 158, in compile_expression
expression = compile_class(cpp_data, mpi_comm=mpi_comm)
File “/Users/Shirley/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/jit/jit.py”, line 170, in compile_class
raise RuntimeError(“Unable to compile C++ code with dijitso”)
RuntimeError: Unable to compile C++ code with dijitso

Also, I have defined the function to be 1 at the subdomain Omega1 and it should be zero outside. I also defined a subdomain around Omega1 to set the function from 1 to 0. Is there a better way to do this?

2

The error occurs because you have an extra ‘)’ inside the expression. The correct f should be:

f = Expression(' ( (x[0] >= 0.35 + DOLFIN_EPS)&&(x[0] <= 0.7 + DOLFIN_EPS)&&(x[1] >= 0.35 + DOLFIN_EPS)&&(x[1] <= 0.7 + DOLFIN_EPS) ) ? p1 : ( ( (x[0] < 0.35 + DOLFIN_EPS)&&(x[0] > 0.7 + DOLFIN_EPS)&&(x[1] < 0.35 + DOLFIN_EPS)&&(x[1] > 0.7 + DOLFIN_EPS) ) ? p2 : 0 ) ',
             p1=p1, p2=p2, degree=2 )

Note that the second condition in this expression, i.e.

(x[0] < 0.35 + DOLFIN_EPS)
&& (x[0] > 0.7 + DOLFIN_EPS)
&& (x[1] < 0.35 + DOLFIN_EPS)
&& (x[1] > 0.7 + DOLFIN_EPS)

will always be false. The && operator is the C++ logical “and” operator. x[0] cannot simultaneously be less than 0.35 and greater than 0.7. It looks like you are trying to define a source term that is equal to 1 inside the subdomain \bigl( [0.35,0.7]\times[0.35,0.7]\bigr) and zero outside. If that is the case, the expression

f = Expression(' ( (x[0] >= 0.35 - DOLFIN_EPS)&&(x[0] <= 0.7 + DOLFIN_EPS)&&(x[1] >= 0.35 - DOLFIN_EPS)&&(x[1] <= 0.7 + DOLFIN_EPS) ) ? p1 : p2',
             p1=p1, p2=p2, degree=2 )

is sufficient to achieve what you desire. Note that I have altered some of the conditionals, i.e.

(x[0] >= 0.35 - DOLFIN_EPS)
(x[1] >= 0.35 - DOLFIN_EPS)

by subtracting DOLFIN_EPS instead of adding it. This allows the point 0.35 to be included in \Omega_1 rather than excluded.

The best way to do this would be to create a mesh that conforms with \Omega_1 (e.g., using gmsh), apply a marker to \Omega_1, then create a Measure to integrate the source term only over \Omega_1. Your current approach using an Expression will not create a piecewise function, since it must first be interpolated into a finite element space, and your mesh elements do not conform with the boundary of \Omega_1.

3 Likes