pkg load fem-fenics msh
problem = 'HyperElasticity';
import_ufl_Problem (problem);

Rotation = @(x,y,z) ...
   [0; ...
   0.5*(0.5 + (y - 0.5)*cos(pi/3) - (z-0.5)*sin(pi/3) - y);...
   0.5*(0.5 + (y - 0.5)*sin(pi/3) + (z-0.5)*cos(pi/3) - z)];

#Create mesh and define function space
x = y = z = linspace (0, 1, 17);
mshd = Mesh (msh3m_structured_mesh (x, y, z, 1, 1:6));
V  = FunctionSpace (problem, mshd);

# Create Dirichlet boundary conditions
bcl = DirichletBC (V, @(x,y,z) [0; 0; 0], 1);
bcr = DirichletBC (V, Rotation, 2);
bcs = {bcl, bcr};
 
# Define source and boundary traction functions
B = Constant ('B', [0.0; -0.5; 0.0]);
T = Constant ('T', [0.1; 0.0; 0.0]);

# Set material parameters
E = 10.0;
nu = 0.3;
mu = Constant ('mu', E./(2*(1+nu)));
lmbda = Constant ('lmbda', E*nu./((1+nu)*(1-2*nu)));
u = Expression ('u', @(x,y,z) [0; 0; 0]);

# Create (linear) form defining (nonlinear) variational problem
L = ResidualForm (problem, V, mu, lmbda, B, T, u);

# Solve nonlinear variational problem F(u; v) = 0

# Create function for the resolution of the NL problem
  function [y, jac] = f (problem, xx, V, bc1, bc2, B, T, mu, lmbda)
    u = Function ('u', V, xx);
    a = JacobianForm (problem, V, mu, lmbda, u);
    L = ResidualForm (problem, V, mu, lmbda, B, T, u);
    if (nargout == 1)
      [y, xx] = assemble (L, xx, bc1, bc2);
    elseif (nargout == 2)
      [jac, y, xx] = assemble_system (a, L, xx, bc1, bc2);
    endif
  endfunction

u0 = assemble (L, bcs{:});
fs = @(xx) f (problem, xx, V, bcl, bcr, B, T, mu, lmbda);
[x, fval, info] = fsolve (fs, u0, optimset ("jacobian", "on"));
func = Function ('u', V, x);
 
# Save solution in VTK format
save (func, 'displacement');

# Plot solution
plot (func);