mesh = Mesh('mesh.xml')
For a scalar function (pressure acoustics)
V = FunctionSpace(mesh,'CG',2)
p = TrialFunction(V) #pressure trial function
dp = TestFunction(V) #pressure test function
for elasticity:
V = VectorFunctionSpace(mesh,'CG',2)
u = TrialFunction(V) #Displacement trial function
du = TestFunction(V) #Displacement test function
The inhomogeneous Helmholtz equation for pressure is: $$ \nabla^2p+k^2p =f$$ The weak form of the above equation is: $$ \frac{k^2}{\rho}\int \limits_\Omega p \delta p \, d \Omega - \frac{1}{\rho}\int \limits_\Omega \nabla p \cdot \nabla \delta p \,d \Omega + \frac{1}{\rho}\int \limits_{\Gamma_N} \delta p (\nabla p \cdot n) \,d\Gamma_N =\frac{1}{\rho}\int \limits_\Omega f\delta p \, d \Omega \\ p = 0 \quad \text{on} \quad \Gamma_D$$ This is written is FEniCs as
n = FacetNormal(mesh)
StiffnessMatrix = 1/rho*dot(grad(p),grad(dp))*dx
MassMatrix = 1/rho*p*dp*dx
NeumannBC = 1/rho*dp*dot(p,n)*ds
forcingFcn = 1/rho*dp*f*dx
WeakForm = -Stiffness + k^2*MassMatrix + NeumannBC
The stiffness tensor is $$ C_{ijkl} = K\delta_{ij}\delta_{kl} +\mu(\delta_{ik}\delta_{jl} +\delta_{il}\delta_{jk} -\frac{2}{3}\delta_{ij}\delta_{kl})$$ In FEniCS, this is implemented as:
I = identity(2)# Second order identity tensor
C = as_tensor(K*I[i,j]*I[k,l] + mu*(I[i,k]*I[j,l] + I[i,l]*I[j,k] - 2./3*I[i,j]*I[k,l],(i,j,k,l))
The weak form of linear elasticity is: $$ \int \limits_\Omega [\sigma(u)_{ij}\varepsilon(\delta u)_{ij} - \omega^2\rho \delta u_i u_i] \,d\Omega -\int \limits_\Gamma \sigma_{ij}n_j\delta u_i \,d \Gamma=0$$ This is implemented in FEniCS as:
StrainU = 0.5*(grad(u)+transpose(grad(u))) = sym(grad(u))
StrainDu = 0.5*(grad(du)+transpose(grad(du)))
StressU = as_tensor(C[i,j,k,l]*StrainU[k,l],(i,j))
WeakForm = inner(StressU,StrainDu)*dx - omega**2*rho*dot(u,du)*dx - p*dot(du,n)*ds
F = grad(u)+I
C = F.T*F #Right Cauchy-Green tensor
E = 0.5*(C-I) #Finite Strain tensor
W = lmbda/2*tr(E)**2+mu*tr(E*E) # St. Venant-Kirchhoff Strain-energy density
E = variable(E)
S = diff(W,E) #Second Piola-Stress tensor
class GammaN(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], 0, tol)
class GammaD(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], 0, tol)
NeumannBoundary = GammaN()
DirichletBoundary = GammaD()
boundary_markers = FacetFunction('size_t', mesh)
boundary_markers.set_all(0)
DirichletBoundary.mark(boundary_markers,1)
NeumannBoundary.mark(boundary_markers,2)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
bc = DirichletBC(V,Constant(0.0),DirichletBoundary)
solve(WeakForm==forcingFcn,u,bc)
file = File("u.pvd")
file << u
plot(u)