@@ Line 302: / Line 302: @@
 </syntaxhighlight>}}
 <div style="width: 100%;">
     <div style="float:left; width: 48%">
 {{Code|Define HyperElasticity problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px">
+pkg load fem-fenics msh
+problem = 'HyperElasticity';
+import_ufl_Problem (problem);
+Rotation = @(x,y,z) ...
+ [0; ...
+.5*(0.5 + (y - 0.5)*cos(pi/3) - (z-0.5)*sin(pi/3) - y);...
+.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 = Expression ('B', @(x,y,z) [0.0; -0.5; 0.0]);
+T = Expression ('T', @(x,y,z) [0.1; 0.0; 0.0]);
+# Set material parameters
+E = 10.0;
+nu = 0.3;
+mu = Expression ('mu', @(x,y,z) E./(2*(1+nu)));
+lmbda = Expression ('lmbda', @(x,y,z) 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
+u0 = assemble (L, bcs{:});
+# 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
+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);
@@ Line 314: / Line 449: @@
     <div style="float:right; width: 48%">
 {{Code|Define HyperElasticity problem with fenics c++|<syntaxhighlight lang="cpp" style="font-size:13px">
+#include <dolfin.h>
+#include "HyperElasticity.h"
+using namespace dolfin;
+// Sub domain for clamp at left end
+class Left : public SubDomain
+{
+  bool inside(const Array<double>& x, bool on_boundary) const
+  {
+    return (std::abs(x[0]) < DOLFIN_EPS) && on_boundary;
+  }
+};
+// Sub domain for rotation at right end
+class Right : public SubDomain
+{
+  bool inside(const Array<double>& x, bool on_boundary) const
+  {
+    return (std::abs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary;
+  }
+};
+// Dirichlet boundary condition for clamp at left end
+class Clamp : public Expression
+{
+public:
+  Clamp() : Expression(3) {}
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    values[0] = 0.0;
+    values[1] = 0.0;
+    values[2] = 0.0;
+  }
+};
+// Dirichlet boundary condition for rotation at right end
+class Rotation : public Expression
+{
+public:
+  Rotation() : Expression(3) {}
+  void eval(Array<double>& values, const Array<double>& x) const
+  {
+    const double scale = 0.5;
+    // Center of rotation
+    const double y0 = 0.5;
+    const double z0 = 0.5;
+    // Large angle of rotation (60 degrees)
+    double theta = 1.04719755;
+    // New coordinates
+    double y = y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta);
+    double z = z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta);
+    // Rotate at right end
+    values[0] = 0.0;
+    values[1] = scale*(y - x[1]);
+    values[2] = scale*(z - x[2]);
+  }
+};
+int main()
+{
+  // Create mesh and define function space
+  UnitCubeMesh mesh (16, 16, 16);
+  HyperElasticity::FunctionSpace V(mesh);
+  // Define Dirichlet boundaries
+  Left left;
+  Right right;
+  // Define Dirichlet boundary functions
+  Clamp c;
+  Rotation r;
+  // Create Dirichlet boundary conditions
+  DirichletBC bcl(V, c, left);
+  DirichletBC bcr(V, r, right);
+  std::vector<const BoundaryCondition*> bcs;
+  bcs.push_back(&bcl); bcs.push_back(&bcr);
+  // Define source and boundary traction functions
+  Constant B(0.0, -0.5, 0.0);
+  Constant T(0.1,  0.0, 0.0);
+  // Define solution function
+  Function u(V);
+  // Set material parameters
+  const double E  = 10.0;
+  const double nu = 0.3;
+  Constant mu(E/(2*(1 + nu)));
+  Constant lambda(E*nu/((1 + nu)*(1 - 2*nu)));
+  // Create (linear) form defining (nonlinear) variational problem
+  HyperElasticity::ResidualForm F(V);
+  F.mu = mu; F.lmbda = lambda; F.B = B; F.T = T; F.u = u;
+  // Create jacobian dF = F' (for use in nonlinear solver).
+  HyperElasticity::JacobianForm J(V, V);
+  J.mu = mu; J.lmbda = lambda; J.u = u;
+  // Solve nonlinear variational problem F(u; v) = 0
+  solve(F == 0, u, bcs, J);
+  // Save solution in VTK format
+  File file("displacement.pvd");
+  file << u;
+  // Plot solution
+  plot(u);
+  interactive();
+  return 0;
+}
 </syntaxhighlight>}}

Fem-fenics: Difference between revisions

Fem-fenics (edit)

Revision as of 16:32, 6 September 2013