Issue with cross product in variational form

Question

Hello
I am trying to model a wave propagating in a solid like Aluminum. My domain is a rectangle with a circular cutout in the middle and the wave is propagating from a source point located at the left edge. The source point has an equation like: Sin(omega*t). I have done it once in a scalar function space. Now I want to do that in a vectorial function space. When I run my code I get this error:
" Cross product requires arguments of rank 1"
Here you can find my manipulation for deriving the variational forms:

And here you can find my code

from dolfin import *
import numpy as np
import math
import mshr
# Aluminium material properties

Poisson=0.33
E=70*(10**9)
dens=2700
# G and Lambda

G=E/(2*(1+Poisson))
lambd=E*Poisson/((1+Poisson)*(1-2*Poisson))
# Wave speed (Dilational and Shear waves)

Cp=math.sqrt((lambd+2*G)/dens)
Cs=math.sqrt(G/dens)
# Domain

domain = mshr.Rectangle(Point(0,0), Point(6,3))-mshr.Circle(Point(3,1.5), 0.2)
# Mesh

mesh = mshr.generate_mesh(domain, 30)
# Mesh refinement

mesh1 = refine(mesh)


V=VectorFunctionSpace(mesh1, "CG", 1, dim=2)
# Time variables
dt = 0.000004; t = 0; T = 0.0012

# Previous and current solution (Dilational and Shear waves)
W_const = Constant((0.0, 0.0))
u1= interpolate(W_const, V)
u0= interpolate(W_const, V)


# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)

a = inner(u,v)*dens*dx + dt*dt*(lambd+2*G)*inner(grad(u), grad(v))*dx+dt*dt*(lambd+G)*cross(grad(v),curl(u))*dx

L = 2*dens*inner(u1,v)*dx-inner(u0,v)*dx

#Boundary condition

bc = DirichletBC(V, 0, "on_boundary")
A, b = assemble_system(a, L, bc)

u=Function(V)

# Saving results

f = File("gg.pvd")
#Solve

while t <= T:
    A, b = assemble_system(a, L, bc)

    delta = PointSource(V, Point(0., 1.5), sin(Cp * 5* t))
    delta.apply(b)

    solve(A, u.vector(), b)

    u0.assign(u1)
    u1.assign(u)

    t += dt
    plot(a, interactive=False)
    f<<u

I think the problem might be in the derivation of the weak forms but I am not sure. I also changed the product from cross to inner. This time I got this error:
“Shape mismatch”
Does anybody happen to know where the problem might be?

Answer 1

You use a notation here that is different from what is used in FEniCS.
The term related to (lambda + mu) would be grad(div(u)) and its
weak form div(u)div(v)dx. Nabla squared is equal to div(grad(u)) which
corresponds to inner(grad(u), grad(v))*dx. Have a look at the demos
related to linear elasticity.

Comments

Thanks! It solved my problem.

---

Hi,

I tried implementing (div(u),div(v))*dx in your code. But I not able to get the solution. It is giving me this error:

*** Error: Unable to create Dirichlet boundary condition.
*** Reason: Expecting a vector-valued boundary value but given function is scalar.
*** Where: This error was encountered inside DirichletBC.cpp.

Can you please share your working code. Thank you.

Ratnesh

---

Can you show the code?

---

Please see the code below:

from dolfin import *
import numpy as np
import math
import mshr

Poisson=0.33
E=70*(10**9)
dens=2700

G=E/(2*(1+Poisson))
lambd=E*Poisson/((1+Poisson)*(1-2*Poisson))

Cp=math.sqrt((lambd+2*G)/dens)
Cs=math.sqrt(G/dens)

domain = mshr.Rectangle(Point(0,0), Point(6,3))-mshr.Circle(Point(3,1.5), 0.2)
mesh = mshr.generate_mesh(domain, 30)
mesh1 = refine(mesh)

V=VectorFunctionSpace(mesh1, "CG", 1, dim=2)

dt = 0.000004; t = 0; T = 0.0012

W_const = Constant((0.0, 0.0))
u1= interpolate(W_const, V)
u0= interpolate(W_const, V)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(u,v)*dens*dx + dt*dt*(lambd+2*G)*inner(grad(u), grad(v))*dx+ 
dt*dt*inner(div(u),div(v))*dx
L = 2*dens*inner(u1,v)*dx-inner(u0,v)*dx

bc = DirichletBC(V, 0, "on_boundary")
A, b = assemble_system(a, L, bc)

u=Function(V)

f = File("gg.pvd")

while t <= T:
   A, b = assemble_system(a, L, bc)
   delta = PointSource(V, Point(0., 1.5), sin(Cp * 5* t))
   delta.apply(b)

   solve(A, u.vector(), b)

   u0.assign(u1)
   u1.assign(u)

   t += dt
   plot(a, interactive=False)
   f<<u

---

div(u) returns a scalar so you should not use inner product.
instead: div(u)div(v)dx

---

I still getting the same error when I do div(u)div(v)dx

Here is my code:

# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)

a = inner(u,v)*dens*dx + dt*dt*(lambd+2*G)*inner(grad(u), grad(v))*dx+\
dt*dt*div(u)*div(v)*dx

L = 2*dens*inner(u1,v)*dx-inner(u0,v)*dx

#Boundary conditionhttp://wiki.octave.org/Fem-fenics

bc = DirichletBC(V, 0, "on_boundary")
A, b = assemble_system(a, L, bc)

u=Function(V)