60
edits
No edit summary |
No edit summary |
||
Line 302: | Line 302: | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
<div style="width: 100%;"> | <div style="width: 100%;"> | ||
<div style="float:left; width: 48%"> | <div style="float:left; width: 48%"> | ||
{{Code|Define HyperElasticity problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | {{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; ... | |||
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 = 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%"> | <div style="float:right; width: 48%"> | ||
{{Code|Define HyperElasticity problem with fenics c++|<syntaxhighlight lang="cpp" style="font-size:13px"> | {{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>}} | </syntaxhighlight>}} |
edits