I have a rather simple script to solve the Poisson equation with a source and a sink term on a certain boundary. After calculating the steady-state solution, I would like to compute some derived quantities from it. I understand that for saving these quantities to a file, I need to convert them from their UFL representation by either using the
from fenics import * import numpy as np # Define tolerance tol = 1E-14 # Create mesh and define function space mesh = IntervalMesh(100, -0.0005, 0.001) V = FunctionSpace(mesh, 'P', 1) # Define subdomains for regions c and a subdomain_c = CompiledSubDomain('x <= mesh_boundary + tol', mesh_boundary=0., tol=tol) subdomain_a = CompiledSubDomain('x >= mesh_boundary - tol', mesh_boundary=0., tol=tol) # Define regions to set coefficients on regions = MeshFunction('size_t', mesh, mesh.topology().dim()) subdomain_c.mark(regions, 0) subdomain_a.mark(regions, 1) # Define coefficient class class Variable_Coefficient(UserExpression): def __init__(self, regions, c_0, c_1, **kwargs): super().__init__(**kwargs) self.regions = regions self.c_0 = c_0 self.c_1 = c_1 def eval_cell(self, values, x, cell): if self.regions[cell.index] == 0: values = self.c_0 elif self.regions[cell.index] == 1: values = self.c_1 # Assign coefficient values D = Variable_Coefficient(regions, 1E-6, 1E-7, degree=0) G = Variable_Coefficient(regions, 1, 0., degree=0) R = Variable_Coefficient(regions, 0., 0.5, degree=0) u_max = Constant(2.5) u_sink = Constant(1.) # Define variational problem u = Function(V) v = TestFunction(V) S = conditional(lt(u, u_sink), u, u_sink) F = D * dot(grad(u), grad(v)) * dx - G * u_max * v * dx + G * u * v * dx + R * S * v * dx # Compute solution solve(F == 0, u) # Compute decomposition and redeposition rates decomposition = project(G * u_max - G * u, V) redeposition = project(R * S, V) #this works #decomposition = interpolate(G * u_max - G * u, V) #redeposition = interpolate(R * S, V) #this does not work # Save solution to file in VTK format vtkfile = File('/mnt/c/Users/margenfeld/Desktop/fenics/output/solution.pvd') vtkfile << u vtkfile << decomposition vtkfile << redeposition
This works with
project, however, I get some aphysical oscillations in the saved quantities. Thus, I wanted to check if this is also the case for
interpolate. When I run the program, I get the error message
Traceback (most recent call last): File "/home/user/fenics/fin_decomposition.py", line 76, in <module> dec = interpolate(decomposition, V) File "/home/user/anaconda3/envs/fenics/lib/python3.8/site-packages/dolfin/fem/interpolation.py", line 73, in interpolate Pv.interpolate(v) File "/home/user/anaconda3/envs/fenics/lib/python3.8/site-packages/dolfin/function/function.py", line 365, in interpolate self._cpp_object.interpolate(u) AttributeError: 'Sum' object has no attribute '_cpp_object'
I was of the impression that
interpolate can be used interchangeably in terms of syntax. That does not seem to be the case, however. I am looking forward to a solution to this problem!