clear all
pkg load msh
pkg load fem-fenics

fem_init_env ();
problem = 'HyperElasticity';
fem_create_all (problem);

x = y = z = linspace (0, 1, 16);
msho = msh3m_structured_mesh (x, y, z, 1, 1:6);
mshd = Mesh (msho);

global V  = HyperElasticity_FunctionSpace (mshd);

global bc2 = DirichletBC (V, @(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)], 2);

global bc1 = DirichletBC (V, @(x,y,z) [0; 0; 0], 1);

global B = Expression ('B', @(x,y,z) [0.0; -0.5; 0.0]);
global T = Expression ('T', @(x,y,z) [0.1; 0.0; 0.0]);
global E = 10.0;
global nu = 0.3;
global mu = Expression ('mu', @(x,y,z) E./(2*(1+nu)));
global lmbda = Expression ('lmbda', @(x,y,z) E*nu./((1+nu)*(1-2*nu)));
u = Expression ('u', @(x,y,z) [0; 0; 0]);

L = HyperElasticity_LinearForm (V, mu, lmbda, B, T, u);
u0 = assemble (L, bc1, bc2);

[x, fval, info] = fsolve (@f, u0, optimset ("jacobian", "on"));

func = Function ('u', V, x);
fem_plot (func);