22
edits
Fem-fenics: Difference between revisions
Jump to navigation
Jump to search
→Hyperelasticity: Add ufl snippet to the example m-file
(→Hyperelasticity: Add ufl snippet to the example m-file) |
|||
| Line 372: | Line 372: | ||
{{Code|Define HyperElasticity problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | {{Code|Define HyperElasticity problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
pkg load fem-fenics msh | pkg load fem-fenics msh | ||
ufl start HyperElasticity | |||
ufl # Function spaces | |||
ufl element = VectorElement("Lagrange", tetrahedron, 1) | |||
ufl | |||
ufl # Trial and test functions | |||
ufl du = TrialFunction(element) # Incremental displacement | |||
ufl v = TestFunction(element) # Test function | |||
ufl | |||
ufl # Functions | |||
ufl u = Coefficient(element) # Displacement from previous iteration | |||
ufl B = Coefficient(element) # Body force per unit volume | |||
ufl T = Coefficient(element) # Traction force on the boundary | |||
ufl | |||
ufl # Kinematics | |||
ufl I = Identity(element.cell().d) # Identity tensor | |||
ufl F = I + grad(u) # Deformation gradient | |||
ufl C = F.T*F # Right Cauchy-Green tensor | |||
ufl | |||
ufl # Invariants of deformation tensors | |||
ufl Ic = tr(C) | |||
ufl J = det(F) | |||
ufl | |||
ufl # Elasticity parameters | |||
ufl mu = Constant(tetrahedron) | |||
ufl lmbda = Constant(tetrahedron) | |||
ufl | |||
ufl # Stored strain energy density (compressible neo-Hookean model) | |||
ufl psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 | |||
ufl | |||
ufl # Total potential energy | |||
ufl Pi = psi*dx - inner(B, u)*dx - inner(T, u)*ds | |||
ufl | |||
ufl # First variation of Pi (directional derivative about u in the direction of v) | |||
ufl F = derivative(Pi, u, v) | |||
ufl | |||
ufl # Compute Jacobian of F | |||
ufl J = derivative(F, u, du) | |||
delete HyperElasticity.ufl | |||
problem = 'HyperElasticity'; | |||
| Line 517: | Line 557: | ||
using namespace dolfin; | using namespace dolfin; | ||
// Sub domain for clamp at left end | // Sub domain for clamp at left end | ||
| Line 626: | Line 706: | ||
// Solve nonlinear variational problem F(u; v) = 0 | // Solve nonlinear variational problem F(u; v) = 0 | ||
solve(F == 0, u, bcs, J); | solve(F == 0, u, bcs, J); | ||