Using functions in expressions

1

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
eigensolver.solve()
r, c, rx, cx = eigensolver.get_eigenpair(0)
print("Eigenvalue: ", r)
ef = Function(S)
ef.vector()[:] = rx
plot(ef[0])
plt.show()

#
# %% 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


l=Function(S)

for n in range(num_steps):   
    # Compute solution
    solve(a_evo == L, l, bc)
    # Save to file 
    
    vtkfile << (l, t)
    #vtkfile2 << (w_temp,t)
    u_m.assign(u_n)
    # Update previous solution
    u_n.assign(l)
    # 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.

2

You could try using a UserExpression, as follows:

import numpy as np
class expr(UserExpression):
    def __init__(self, t, u, r, **kwargs):
        super().__init__(**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)
3

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?

4

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):
        super().__init__(**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
5

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!