Generate an empty form

variational formulation
1

I want to implement an example of left ventricle in the paper, here is the constitutive law.

image
image1287×860 78.5 KB

I have the displacement of the structure and I want to calculate the feedback force. Here is my code.

from fenics import *
def first_PK_stress(u):                     # input the displacement
    
    # Parameters
    C = Constant(10.0)                      # C = 10kPa
    bf = Constant(1.0)
    bt = Constant(1.0)
    bfs = Constant(1.0)

    I = Identity(len(u))
    F = I + grad(u)                         # nabla_grad is used someplaces. I think grad is correct.
    E = 0.5*(F.T*F - I)

    e1 = as_vector([ 1.0, 0.0, 0.0 ])
    e2 = as_vector([ 0.0, 1.0, 0.0 ])
    e3 = as_vector([ 0.0, 0.0, 1.0 ])

    E11, E12, E13 = inner(E*e1, e1), inner(E*e1, e2), inner(E*e1, e3)
    E21, E22, E23 = inner(E*e2, e1), inner(E*e2, e2), inner(E*e2, e3)
    E31, E32, E33 = inner(E*e3, e1), inner(E*e3, e2), inner(E*e3, e3)

    Q = bf*E11**2 + bt*(E22**2 + E33**2 + E23**2 + E32**2) \
      + bfs*(E12**2 + E21**2 + E13**2 + E31**2)

    Wpassive = C/2.0 * (exp(Q) - 1)

    FF = variable(F)
    # CC = variable(F.T*F)
    # S = 2*diff(Wpassive,CC)         # calculate the second PK stress tensor
    return diff(Wpassive, FF)         # Calculate the first PK stress tensor



mesh = UnitCubeMesh(10,10,10)
V    = VectorFunctionSpace(mesh,"P",2)
u = TrialFunction(V)
v = TestFunction(V)
disp = Function(V)

# Define the constutive law for the left ventricle
F = inner(first_PK_stress(disp), grad(v))*dx + inner(u, v)*dx
a = lhs(F)
L = rhs(F)
A = assemble(a)
b = assemble(L)

The information is :

Traceback (most recent call last):
  File "test_force.py", line 45, in <module>
    b = assemble(L)
  File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 198, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 58, in _create_dolfin_form
    function_spaces=function_spaces)
  File "/usr/lib/python3/dist-packages/dolfin/fem/form.py", line 23, in __init__
    self.subdomains, = list(sd.values())  # Assuming single domain
ValueError: not enough values to unpack (expected 1, got 0)

I think the reason is the variational formulation generate an empty form L : Form([]), but I don’t know why.

2

If I project the first Piola-Kirchhoff to a tensor function space and then calculate the variational formulation, I can get the correct L which is not empty any more. What is the difference?

T = TensorFunctionSpace(mesh, "P", 2)
t = project(first_PK_stress(disp),T)
F = inner(t, grad(v))*dx + inner(u, v)*dx
L = rhs(F)
3

The problem is that FF is not defined until after defining Wpassive, so Wpassive does not depend on FF and diff(W,FF) is zero. You need to first define the Variable that Wpassive depends on, next define Wpassive in terms of that Variable, then differentiate Wpassive. You can fix this by defining F as as a Variable right away, i.e.,

F = variable(I + grad(u))

and returning diff(Wpassive,F).

4

Thank you very much!

5

Hello @kamesky !

I correct this mistake but another mistake arises. I want to implement fluid structure interaction model with IBM (immersed boundary method). In IBM, the displacement of the structure is calculated at first, then we calculate the displacement gradient, then the first Piola-Kirchhoff stress tensor, finally the feedback force. The feed back force consists of body force and surface force. Here are the variational formulations from E. Griffith - 2017.They are mathematically equivalent.

\int \mathbf{G}\cdot\mathbf{V}\;d\mathbf{x} = \int(\nabla\cdot\mathbb{P})\cdot \mathbf{V}\;d\mathbf{x}-\int(\mathbb{P}\mathbf{n})\cdot \mathbf{V}\;d\mathbf{x}
\int \mathbf{G}\cdot\mathbf{V}\;d\mathbf{x} = -\int\nabla\mathbb{P}:\nabla \mathbf{V}\;d\mathbf{x}

where \mathbb{P} is the first Piola-Kirchhoff stress tensor, \mathbf{n} is the out normal. Here is my minimum code for the second formulation, which generate extremely large feed back force at the surface when there is no displacement.

from fenics import *
from mshr import *

def first_PK_stress(u):                       # input the displacement

    C = Constant(10.0)                          # C = 10kPa
    bf = Constant(1.0)
    bt = Constant(1.0)
    bfs = Constant(1.0)

    I = Identity(len(u))
    F = variable(I + grad(u))
    E = 0.5*(F.T*F - I)
    
    e1 = as_vector([ 1.0, 0.0, 0.0 ])
    e2 = as_vector([ 0.0, 1.0, 0.0 ])
    e3 = as_vector([ 0.0, 0.0, 1.0 ])

    E11, E12, E13 = inner(E*e1, e1), inner(E*e1, e2), inner(E*e1, e3)
    E21, E22, E23 = inner(E*e2, e1), inner(E*e2, e2), inner(E*e2, e3)
    E31, E32, E33 = inner(E*e3, e1), inner(E*e3, e2), inner(E*e3, e3)

    Q = bf*E11**2 + bt*(E22**2 + E33**2 + E23**2 + E32**2) \
      + bfs*(E12**2 + E21**2 + E13**2 + E31**2)

    Wpassive = C/2.0 * (exp(Q) - 1)

    return diff(Wpassive, F)       # Calculate the first PK stress tensor


# generate the mesh
mesh_size = 20
big = Ellipsoid(Point(0.0,0.0,0.0),10.0,10.0,20.0)
small = Ellipsoid(Point(0.0,0.0,0.0),7.0,7.0,17.0)
gaizi = Box(Point(-11.0,-11.0,5.0),Point(11.0,11.0,21.0))
domain = big-small-gaizi
mesh = generate_mesh(domain, mesh_size)

# define the variables
V = VectorFunctionSpace(mesh,"P",2)
u = TrialFunction(V)
v = TestFunction(V)
disp = interpolate(Expression(("x[0]","x[1]","x[2]"),degree=1), V)
P = first_PK_stress(disp)

# define the variatinal form
F = inner(u, v)*dx + inner(P,grad(v))*dx
a = lhs(F)
L = rhs(F)
A = assemble(a)
b = assemble(L)
f = Function(V)

# solve the equation.
solve(A, f.vector(), b, "cg", "sor")
File("force_1_0.pvd") << f

Here is the result:

image
image857×733 50.6 KB

6

I think there might just be a mix-up between position and displacement here. In the text, you say there is “no displacement”, but that would correspond to disp being zero, not Expression(("x[0]","x[1]","x[2]"),degree=1).

7

Thank you very much, it yields correct results. But I still have a question. Can I separated the formulation:

\int \mathbf{G}\cdot\mathbf{V}\;d\mathbf{x} = \int(\nabla\cdot\mathbb{P})\cdot \mathbf{V}\;d\mathbf{x}-\int(\mathbb{P}\mathbf{n})\cdot \mathbf{V}\;d\mathbf{s}

as two formulations:

\int \mathbf{G_1}\cdot\mathbf{V}\;d\mathbf{x} = \int(\nabla\cdot\mathbb{P})\cdot \mathbf{V}\;d\mathbf{x}
\int \mathbf{G_2}\cdot\mathbf{V}\;d\mathbf{s} = -\int(\mathbb{P}\mathbf{n})\cdot \mathbf{V}\;d\mathbf{s}

I found G\neq G_1+G_2, the results of separated formulations are much smoother and smaller than that of unified formulation. For example, let the position be (x,y,z)\rightarrow(2x,2y,2z), the same as the former post, the results of G_1 and G_2 are

image
image845×747 46.7 KB

image
image820×703 44.7 KB

The reason I use d \mathbf{s} instead of d\mathbf{x} is because it yields more strange results.

\int \mathbf{G_2}\cdot\mathbf{V}\;d\mathbf{x} = -\int(\mathbb{P}\mathbf{n})\cdot \mathbf{V}\;d\mathbf{s}

image
image866×737 57.5 KB

# define the variables
V = VectorFunctionSpace(mesh,"P",2)
u = TrialFunction(V)
v = TestFunction(V)
disp = interpolate(Expression(("x[0]","x[1]","x[2]"),degree=1), V)
P = first_PK_stress(disp)

# define the variatinal form
F = inner(u, v)*dx - inner(div(P),v)*dx
a = lhs(F)
L = rhs(F)
A = assemble(a)
b = assemble(L)
f = Function(V)

# solve the equation.
solve(A, f.vector(), b, "cg", "sor")
File("G_1.pvd") << f

# define the variatinal form
n = FacetNormal(mesh)
F = inner(u, v)*ds + inner(P*n,v)*ds
a = lhs(F)
L = rhs(F)
A = assemble(a, keep_diagonal=True) 
A.ident_zeros()
b = assemble(L)
f = Function(V)

# solve the equation.
solve(A, f.vector(), b, "cg", "sor")
File("G_2.pvd") << f

# define the variatinal form
n = FacetNormal(mesh)
F = inner(u, v)*dx + inner(P*n,v)*ds
a = lhs(F)
L = rhs(F)
A = assemble(a)
b = assemble(L)
f = Function(V)

# solve the equation.
solve(A, f.vector(), b, "cg", "sor")
File("G_2_.pvd") << f