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?