Mercurial > fem-fenics-eugenio
changeset 96:38b31ad6ab4e
New example: resolution of Hyperelastic problem using a minimization technique
* Hyperelasticity.m: script for the resolution of the problem. The problem is
one of the demos presented with dolfin.
* f.m: the problem is non-linear, and we thus use fsolve to solve it.
* Hyperelasticity.ufl: definition of the problem
author | gedeone-octave <marcovass89@hotmail.it> |
---|---|
date | Sun, 11 Aug 2013 22:10:47 +0200 |
parents | 77aabcc087dc |
children | dc9987325fea |
files | example/Elasticity/HyperElasticity.m example/Elasticity/HyperElasticity.ufl example/Elasticity/f.m |
diffstat | 3 files changed, 120 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/Elasticity/HyperElasticity.m Sun Aug 11 22:10:47 2013 +0200 @@ -0,0 +1,35 @@ +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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/Elasticity/HyperElasticity.ufl Sun Aug 11 22:10:47 2013 +0200 @@ -0,0 +1,63 @@ +# Copyright (C) 2009-2010 Harish Narayanan and Garth N. Wells +# +# This file is part of DOLFIN. +# +# DOLFIN is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# DOLFIN is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>. +# +# Modified by Anders Logg 2011 +# +# First added: 2009-09-29 +# Last changed: 2011-06-28 +# +# The bilinear form a(du, v) and linear form L(v) for +# a hyperelastic model. +# +# Compile this form with FFC: ffc -l dolfin -feliminate_zeros -fprecompute_basis_const -fprecompute_ip_const HyperElasticity.ufl + +# Function spaces +element = VectorElement("Lagrange", tetrahedron, 1) + +# Trial and test functions +du = TrialFunction(element) # Incremental displacement +v = TestFunction(element) # Test function + +# Functions +u = Coefficient(element) # Displacement from previous iteration +B = Coefficient(element) # Body force per unit volume +T = Coefficient(element) # Traction force on the boundary + +# Kinematics +I = Identity(element.cell().d) # Identity tensor +F = I + grad(u) # Deformation gradient +C = F.T*F # Right Cauchy-Green tensor + +# Invariants of deformation tensors +Ic = tr(C) +J = det(F) + +# Elasticity parameters +mu = Constant(tetrahedron) +lmbda = Constant(tetrahedron) + +# Stored strain energy density (compressible neo-Hookean model) +psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 + +# Total potential energy +Pi = psi*dx - inner(B, u)*dx - inner(T, u)*ds + +# First variation of Pi (directional derivative about u in the direction of v) +F = derivative(Pi, u, v) + +# Compute Jacobian of F +J = derivative(F, u, du)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/Elasticity/f.m Sun Aug 11 22:10:47 2013 +0200 @@ -0,0 +1,22 @@ +function [y, jac] = f (xx) + + global V; + global mu; + global lmbda; + global T; + global B; + global bc1; + global bc2; + + u = Function ('u', V, xx); + + a = HyperElasticity_BilinearForm (V, mu, lmbda, u); + L = HyperElasticity_LinearForm (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