Using functions in expressions

Hi everyone!

i calculated an Eigenfunction, which i want to use, to generate a load, to my model, i have all the parameters i need, but my console thorws an exception, for the eigenfunction, which is technically a vector with two entries.

mesh = UnitSquareMesh(10,10)
X = FiniteElement('CG', mesh.ufl_cell(), 2)
Y = FiniteElement('CG', mesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(mesh, Z)

u = TrialFunctions(S)
phi = TestFunctions(S)

a = inner(grad(u[1]), grad(phi[0]))*dx + inner(grad(u[0]), grad(phi[1]))*dx

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(S, Constant((0,0)), boundary)

L_dummy = Constant(0)*phi[0]*dx
m = u[0]*phi[0]*dx + u[1]*phi[1]*dx

A, _ = assemble_system(a,L_dummy,bc)
B = assemble(m)
eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B))
eigensolver.parameters['spectrum'] = 'target real'
eigensolver.parameters['tolerance'] = 1e-6
r, c, rx, cx = eigensolver.get_eigenpair(0)
print("Eigenvalue: ", r)
ef = Function(S)
ef.vector()[:] = rx

# %% Computation of time dependent problem

T = 0.1      # final time
num_steps = 500    # number of time steps
dt = T / num_steps # time step size
t = 0

u_m = Function(S)
u_n = Function(S)
u = TrialFunctions(S)
phi = TestFunctions(S)

f = Expression ("-2*r*u_lambda*-2*r*sin(r*t)" , r = r, t = 0, u_lambda = ef[0])

# Create VTK file for saving solution
vtkfile = File('solution/ResonatingPlate.pvd')

#define F and assemble unknown and known variables
a_evo = -dt*dt*(inner(grad(u[0]),grad(phi[0])) + inner(grad(u[1]),grad(phi[1])) + u[1]*phi[0])*dx + (u[0]*phi[1])*dx
L = (2*u_n[0]*phi[1] - u_m[0]*phi[1])*dx + dt*dt*(f*phi[1])*dx


for n in range(num_steps):   
    # Compute solution
    solve(a_evo == L, l, bc)
    # Save to file 
    vtkfile << (l, t)
    #vtkfile2 << (w_temp,t)
    # Update previous solution
    # Update current time
    t += dt
    f.t = t

I have the feeling that there is a simple thing i am missing to use the prior calculated function in my expression. Very hint is very much appreciated.

You could try using a UserExpression, as follows:

import numpy as np
class expr(UserExpression):
    def __init__(self, t, u, r, **kwargs):
        self.t = t
        self.u_lambda = u
        self.r = r
    def eval(self, values, x):
        values[0] = -2*self.r*self.u_lambda[0](x)*(-2)*self.r*np.sin(self.r*self.t)
    def value_shape(self):
        return ()

f = expr(0, ef, r)

Hi, i tried using your UserExpression, it does lead to now further erros in regards to computing, however is still can’t get the desired results. I figured, that the function i want to use for the expression consists of two functions, so the assumption that only using ef[0] proved to be wrong, as i also can’t produce a nice sine wave in the following time dependent computation.
My next idea would be to use both entries. How would i change the UserExpression to realize this?
Or are there more manual like ressources to better build expressions and UserExpression?

You can change the shape of the UserExpression to contain two entries per coordinate, as illustrated below. You could also have a look a this documented demo.

from dolfin import *
mesh = UnitSquareMesh(10,10)
S = VectorFunctionSpace(mesh, "CG", 2)
x = SpatialCoordinate(mesh)
s = project(as_vector((x[0], x[1])), S)
import numpy as np
class expr(UserExpression):
    def __init__(self, t, u, r, **kwargs):
        self.t = t
        self.u_lambda = u
        self.r = r
    def eval(self, values, x):
        values[0] = -2*self.r*self.u_lambda[0](x)*(-2)*self.r*np.sin(self.r*self.t)
        values[1] = -2*self.r*self.u_lambda[1](x)*(-2)*self.r*np.sin(self.r*self.t)
    def value_shape(self):
        return (2, )

r = 1
t = 0.1
user_expr = expr(t,s, r)
print(assemble(inner(user_expr, user_expr)*dx(domain=mesh)))
1 Like

HI, i found a workaround which works really fine and is super simple. It is actually able to split to expression and the function, meaning that the function i calculated doesn’t need to apear in the expression i can just build the expression and multiply it with the function in the following formulation.
Like this:

    f = Constant(2*r*sin(r*t))
    L = (2*u_n[0]*phi[1] - u_m[0]*phi[1])*dx + dt*dt*(inner(f*ef[0],phi[0]))*dx 

this actually produced the desired results.
Thanks for the help!