### How to build functions using subspaces from a mixed element

298
views
0
6 months ago by
I am working on an elasticity PDE, and since I have complex variables, I need to use "MixedElement" in Fenics. For simplicity, here I just copy the code for a problem with real (not complex) values. Here is how I defined the function space:

mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)

Vu_ele = VectorElement('P', mesh.ufl_cell(), 1, dim=3)
V      = FunctionSpace(mesh, MixedElement([Vu_ele, Vu_ele]))​
Vmu    = FunctionSpace(mesh, "P", 1)

Here is how I defined my variables and BC:

# Define boundary condition
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < DOLFIN_EPS

bc     = DirichletBC(V.sub(0), Constant((0., 0., 0.)), clamped_boundary)
bc_adj = DirichletBC(V.sub(1), Constant((0., 0., 0.)), clamped_boundary)

# Define variables
w = Function(V)

u_trial, p_trial  = TrialFunctions(V)
u_test , p_test   = TestFunctions(V)

mu_true  = interpolate(Expression('1', degree=1),Vmu)

u  = w.split(deepcopy=True)[0]
p  = w.split(deepcopy=True)[1]​

d = u.geometric_dimension()  # space dimension
f = Constant((0., 0., 0.))

I get an error when trying to solve this equation:

# Define state equation
def state(mu):
u  = w.split(deepcopy=True)[0]
a  = inner(grad(u_test), mu*(grad(u_trial) + grad(u_trial).T) + lambda_*div(u_trial)*Identity(d))*dx
L  = inner(f, u_test)*dx
# Compute solution
solve(a == L, u, bc)
return u

um = state(mu_true)​

Here is the error:

*** Error: Unable to define linear variational problem a(u, v) = L(v) for all v.
*** Reason: Expecting the solution variable u to be a member of the trial space.
*** Where: This error was encountered inside LinearVariationalProblem.cpp.

Why is u not a member of the trial space? How do I solve this and still use the mixed element and NOT use separate "FunctionSpace" and  "VectorFunctionSpace"?

Here is the whole code (simple version):

from __future__ import print_function
from dolfin import *
import numpy as np

# Scaled variables
L = 1; W = 0.2
#mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)

Vu_ele = VectorElement('P', mesh.ufl_cell(), 1, dim=3)
V      = FunctionSpace(mesh, MixedElement([Vu_ele, Vu_ele]))
Vmu    = FunctionSpace(mesh, "P", 1)

# Define boundary condition
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < DOLFIN_EPS

bc     = DirichletBC(V.sub(0), Constant((0., 0., 0.)), clamped_boundary)
bc_adj = DirichletBC(V.sub(1), Constant((0., 0., 0.)), clamped_boundary)

# Define variables
w = Function(V)

u_trial, p_trial  = TrialFunctions(V)
u_test , p_test   = TestFunctions(V)

mu_true  = interpolate(Expression('1', degree=1),Vmu)

u  = w.split(deepcopy=True)[0]
p  = w.split(deepcopy=True)[1]

d = u.geometric_dimension()  # space dimension
f = Constant((0., 0., 0.))

# Define state equation
def state(mu):
u  = w.split(deepcopy=True)[0]
a  = inner(grad(u_test), mu*(grad(u_trial) + grad(u_trial).T) + lambda_*div(u_trial)*Identity(d))*dx
L  = inner(f, u_test)*dx
# Compute solution
solve(a == L, u, bc)
return u

um = state(mu_true)​
Community: FEniCS Project

### 1 Answer

4
6 months ago by
The "splitting" of mixed functions in FEniCS is always source of confusion and errors. There is known issue related to this, see https://bitbucket.org/fenics-project/dolfin/issues/194/split-u-versus-usplit

First of all, when you call solve(a == L, u, bc) the object u must be mixed function, i.e. of class Function which is a member of mixed function space, that is e.g. your w. Your u is a Function, that is correct, but it is a function on vector element space of degree 1.
Just replace the solve line with solve(a == L, w, bc). The solution will be mixed function and u = w.split(deepcopy=True)[0] will get you the u-part of it.

Btw. don't know what equations are you trying to compute, but you dont have equations for p, system is not complete.
EDIT: Now I see the note about complex variable -- do not forget to add it to the system.

Thank you Michal for the quick response. Your suggestion helps this simple version which only has real variables (as you noted) run. I still have a question though. You mentioned that u is a Function on vector element space of degree 1. Aren't u_trial and u_test also trial function and test function on the same vector element space? If yes, then what is the source of the error? If not, what space do they belong to then?

written 6 months ago by Navid Nazari
I am definitely not an UFL (and FFC) expert. But very briefly, dolfin.split(w) returns you a symbolic UFL-way split function w. You can use them inside UFL form and assemble matrices/vectors/scalars, to set up system of equations, couple the equations, ... . This symbolic split doesn't produce a fully-fledged dolfin Function, you cannot save it = you cannot see the actual expansion coefficients wrt basis.

On the other hand, w.split(True) gives you functions as true members of individual finite element spaces, you can save them for visualization. But they are  "withdrawn" from the mixed space context, they do not refer to any "parent" space, you can't tell they are actually siblings (at least I do not know how).

EDIT: You call TrialFunctions and TestFunctions methods. They work essentially the same as dolfin.split() described above.
written 6 months ago by Michal Habera
I just looked at the results generated from running the code after implementing your suggestion, and all the values reported by print(u.vector()[i]) for any i is [nan]. This happens even though the code finishes its run without any errors. Do you have any idea what might be going on?

Edit: I know that the code body, other than definition of the function spaces, works because it gives me correct answers if I use V = VectorFunctionSpace(...) and make the necessary changes.
written 6 months ago by Navid Nazari
Your system is incomplete, the matrix is singular, what do you expect to get? Complete the system by adding equations for p...

If you do not want to take care of p, then add d.inner(p_test, p_trial)*dx
written 6 months ago by Michal Habera
Oh, right. I forgot that I was solving for w now. Thanks again!
written 6 months ago by Navid Nazari
Please login to add an answer/comment or follow this question.

Similar posts:
Search »