Fem-fenics: Difference between revisions

1,289 bytes added ,  28 May 2014
→‎Hyperelasticity: Add ufl snippet to the example m-file
(Undo revision 4907 by Eg123 (talk))
(→‎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
problem = 'HyperElasticity';
import_ufl_Problem (problem);




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);




22

edits