Storing arbitrary values on node edges

On 6 September 2012 00:41, jon pry <email address hidden>wrote:

> New question #207800 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/207800
>
> Hi,
>
> I am developing some mechanical code based on the hyperelastic python
> demo. I can successfully calculate the displacement of the beam under load
> but I would also like to know the stresses. It seems that the calculation
> of the stress is intertwined with the calculation of the displacement since
> the displacement simply seeks to minimize the stress. So the easiest thing
> to do would be to save the calculated stress around the vector onto the
> element currently being used. I can't find any way to do that.
> The only path forward I've come up with is to solve the whole thing
> again after deforming the mesh by the results of the first pass, but this
> seems like a lot of work to retrieve numbers already calculated. I've
> attached the code. All I want is to plot the "psi" variable without
> breaking my back. What is the easiest way to do this?
>

You could do:

psi_proj = project(psi, some_space)
plot(psi_proj)

Kent

>
> """ This demo program solves a hyperelastic problem. It is implemented
> in Python by Johan Hake following the C++ demo by Harish Narayanan"""
>
> # Copyright (C) 2008-2010 Johan Hake 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 Harish Narayanan 2009
> # Modified by Anders Logg 2011
> #
> # First added: 2009-10-11
> # Last changed: 2010-08-28
>
> # Begin demo
>
> from dolfin import *
>
> # Optimization options for the form compiler
> parameters["form_compiler"]["cpp_optimize"] = True
> ffc_options = {"optimize": True, \
> "eliminate_zeros": True, \
> "precompute_basis_const": True, \
> "precompute_ip_const": True}
>
> # Create mesh and define function space
> mesh = UnitCube(16, 16, 16)
> Vs = VectorFunctionSpace(mesh, "Lagrange", 1)
> #Fs = FunctionSpace(mesh, "Lagrange", 1)
>
> # Mark boundary subdomians
> left, right = compile_subdomains(["(std::abs(x[0]) < DOLFIN_EPS) &&
> on_boundary",
> "(std::abs(x[0] - 1.0) < DOLFIN_EPS) &&
> on_boundary"])
>
> # Define Dirichlet boundary (x = 0 or x = 1)
> c = Expression(("0.0", "0.0", "0.0"))
> r = Expression(("scale*0.0",
> "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] -
> z0)*sin(theta) - x[1])",
> "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] -
> z0)*cos(theta) - x[2])"),
> scale = 0.1, y0 = 0.5, z0 = 0.5, theta = 0)
> bcl = DirichletBC(Vs, c, left)
> bcr = DirichletBC(Vs, c, right)
> bcs = [bcl, bcr]
>
> def epsilon(v):
> Dv = grad(v)
> return 0.5*(Dv + Dv.T)
>
> # Define functions
> du = TrialFunction(Vs) # Incremental displacement
> v = TestFunction(Vs) # Test function
> u = Function(Vs) # Displacement from previous iteration
> p = Function(Vs)
> B = Constant((0.0, -5.0, 0.0)) # Body force per unit volume
> T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
>
> # Kinematics
> I = Identity(Vs.cell().d) # Identity tensor
> F = I + grad(u) # Deformation gradient
> C = F.T*F # Right Cauchy-Green tensor
> Ek = .5 * (C - I) # Green–Lagrange finite strain
> tensor
>
> # Invariants of deformation tensors
> Ic = tr(C)
> Iic = (Ic**2 - tr(C*C))/2
> J = det(F)
>
> # Elasticity parameters
> E, nu = 10.0, 0.3
> mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
> K = 4000.0# Constant(E/(2*(1-nu)))
> c1 = 10.0#Constant(.1)
> c2 = 20.0#Constant(.2)
>
> # Stored strain energy density (compressible neo-Hookean model)
> psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
>
> # Stored strain energy density (incompressible mooney-rivlin model)
> #psi = (mu/2)*(Ic - 3) - mu*(J-1) + (lmbda/2)*(J-1)**2
> #psi = (mu/2)*(Ic - 3) + mu*(Iic-3) + (K/2)*(ln(J))**2
>
>
>
> #a = inner(epsilon(v), epsilon(u))*dx
> #p.vector()[:] = a
>
> # St Venant Kirchoff model
> #psi = mu*tr(Ek*Ek) + lmbda/2.0*tr(Ek)**2
> # Hooke's law
> #psi = tr(0*lmbda*tr(Ek)*I) - tr(0*2.0*mu*Ek)# +
>
> # Total potential energy
> Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
>
> F = derivative(Pi, u, v)
> # Compute first variation of Pi (directional derivative about u in the
> direction of v)
>
> # Compute Jacobian of F
> J = derivative(F, u, du)
>
> # Solve variational problem
> solve(F == 0, u, bcs, J=J,
> form_compiler_parameters=ffc_options)
>
> # Save solution in VTK format
> file = File("displacement.pvd");
> file << u;
>
> # Plot and hold solution
> f0 = Expression(('x[0]','x[1]','x[2]'))
> f0 = interpolate(f0, Vs)
> plot(tr(grad(u)))
> plot(u, mode = "displacement")
> interactive()
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>

To get the stress as a Function (so you can plot and examine the values), follow these steps:

1. Compute the solution as you do know.

2. (re)define the stress in terms of you solution Function u.

3. Project the stress expression onto a matrix-valued space by calling the project function.

--
Anders

6 sep 2012 kl. 00:41 skrev jon pry <email address hidden>:

> New question #207800 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/207800
>
> Hi,
>
> I am developing some mechanical code based on the hyperelastic python demo. I can successfully calculate the displacement of the beam under load but I would also like to know the stresses. It seems that the calculation of the stress is intertwined with the calculation of the displacement since the displacement simply seeks to minimize the stress. So the easiest thing to do would be to save the calculated stress around the vector onto the element currently being used. I can't find any way to do that.
> The only path forward I've come up with is to solve the whole thing again after deforming the mesh by the results of the first pass, but this seems like a lot of work to retrieve numbers already calculated. I've attached the code. All I want is to plot the "psi" variable without breaking my back. What is the easiest way to do this?
>
> """ This demo program solves a hyperelastic problem. It is implemented
> in Python by Johan Hake following the C++ demo by Harish Narayanan"""
>
> # Copyright (C) 2008-2010 Johan Hake 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 Harish Narayanan 2009
> # Modified by Anders Logg 2011
> #
> # First added: 2009-10-11
> # Last changed: 2010-08-28
>
> # Begin demo
>
> from dolfin import *
>
> # Optimization options for the form compiler
> parameters["form_compiler"]["cpp_optimize"] = True
> ffc_options = {"optimize": True, \
> "eliminate_zeros": True, \
> "precompute_basis_const": True, \
> "precompute_ip_const": True}
>
> # Create mesh and define function space
> mesh = UnitCube(16, 16, 16)
> Vs = VectorFunctionSpace(mesh, "Lagrange", 1)
> #Fs = FunctionSpace(mesh, "Lagrange", 1)
>
> # Mark boundary subdomians
> left, right = compile_subdomains(["(std::abs(x[0]) < DOLFIN_EPS) && on_boundary",
> "(std::abs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary"])
>
> # Define Dirichlet boundary (x = 0 or x = 1)
> c = Expression(("0.0", "0.0", "0.0"))
> r = Expression(("scale*0.0",
> "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
> "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
> scale = 0.1, y0 = 0.5, z0 = 0.5, theta = 0)
> bcl = DirichletBC(Vs, c, left)
> bcr = DirichletBC(Vs, c, right)
> bcs = [bcl, bcr]
>
> def epsilon(v):
> Dv = grad(v)
> return 0.5*(Dv + Dv.T)
>
> # Define functions
> du = TrialFunction(Vs) # Incremental displacement
> v = TestFunction(Vs) # Test function
> u = Function(Vs) # Displacement from previous iteration
> p = Function(Vs)
> B = Constant((0.0, -5.0, 0.0)) # Body force per unit volume
> T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
>
> # Kinematics
> I = Identity(Vs.cell().d) # Identity tensor
> F = I + grad(u) # Deformation gradient
> C = F.T*F # Right Cauchy-Green tensor
> Ek = .5 * (C - I) # Green–Lagrange finite strain tensor
>
> # Invariants of deformation tensors
> Ic = tr(C)
> Iic = (Ic**2 - tr(C*C))/2
> J = det(F)
>
> # Elasticity parameters
> E, nu = 10.0, 0.3
> mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
> K = 4000.0# Constant(E/(2*(1-nu)))
> c1 = 10.0#Constant(.1)
> c2 = 20.0#Constant(.2)
>
> # Stored strain energy density (compressible neo-Hookean model)
> psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
>
> # Stored strain energy density (incompressible mooney-rivlin model)
> #psi = (mu/2)*(Ic - 3) - mu*(J-1) + (lmbda/2)*(J-1)**2
> #psi = (mu/2)*(Ic - 3) + mu*(Iic-3) + (K/2)*(ln(J))**2
>
>
>
> #a = inner(epsilon(v), epsilon(u))*dx
> #p.vector()[:] = a
>
> # St Venant Kirchoff model
> #psi = mu*tr(Ek*Ek) + lmbda/2.0*tr(Ek)**2
> # Hooke's law
> #psi = tr(0*lmbda*tr(Ek)*I) - tr(0*2.0*mu*Ek)# +
>
> # Total potential energy
> Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
>
> F = derivative(Pi, u, v)
> # Compute first variation of Pi (directional derivative about u in the direction of v)
>
> # Compute Jacobian of F
> J = derivative(F, u, du)
>
> # Solve variational problem
> solve(F == 0, u, bcs, J=J,
> form_compiler_parameters=ffc_options)
>
> # Save solution in VTK format
> file = File("displacement.pvd");
> file << u;
>
> # Plot and hold solution
> f0 = Expression(('x[0]','x[1]','x[2]'))
> f0 = interpolate(f0, Vs)
> plot(tr(grad(u)))
> plot(u, mode = "displacement")
> interactive()
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.