Editing Fem-fenics
Jump to navigation
Jump to search
The edit can be undone. Please check the comparison below to verify that this is what you want to do, and then publish the changes below to finish undoing the edit.
Latest revision | Your text | ||
Line 1: | Line 1: | ||
Package for solving Partial Differential Equations based on Fenics. | Package for solving Partial Differential Equations based on Fenics. | ||
== Tutorials == | == Tutorials == | ||
=== Poisson Equation === | === Poisson Equation === | ||
Here is a first example for the solution of the Poisson equation. | Here is a first example for the solution of the Poisson equation. | ||
The equation being solved is | The equation being solved is | ||
<div style="width: 100%;"> | <div style="width: 100%;"> | ||
Line 33: | Line 10: | ||
{{Code|Define Poisson problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | {{Code|Define Poisson problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
pkg load fem-fenics msh | pkg load fem-fenics msh | ||
import_ufl_Problem ('Poisson') | |||
# Create mesh and define function space | # Create mesh and define function space | ||
Line 50: | Line 16: | ||
mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4)); | mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4)); | ||
# File Poisson.ufl | |||
# element = FiniteElement("Lagrange", triangle, 1) | |||
V = FunctionSpace('Poisson', mesh); | V = FunctionSpace('Poisson', mesh); | ||
Line 61: | Line 29: | ||
# Define variational problem | # Define variational problem | ||
# | # File Poisson.ufl | ||
# u = TrialFunction(element) | # u = TrialFunction(element) | ||
# v = TestFunction(element) | # v = TestFunction(element) | ||
Line 71: | Line 39: | ||
g = Expression ('g', @(x,y) sin (5.0 * x)); | g = Expression ('g', @(x,y) sin (5.0 * x)); | ||
# | # File Poisson.ufl | ||
# a = inner(grad(u), grad(v))*dx | # a = inner(grad(u), grad(v))*dx | ||
# L = f*v*dx + g*v*ds | # L = f*v*dx + g*v*ds | ||
a = BilinearForm ('Poisson' | a = BilinearForm ('Poisson', V); | ||
L = LinearForm ('Poisson', V, f, g); | L = LinearForm ('Poisson', V, f, g); | ||
Line 99: | Line 67: | ||
# Create mesh and define function space | |||
mesh = UnitSquareMesh(32, 32) | |||
V = FunctionSpace(mesh, "Lagrange", 1) | V = FunctionSpace(mesh, "Lagrange", 1) | ||
Line 131: | Line 90: | ||
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + | f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)") | ||
g = Expression("sin(5*x[0])") | g = Expression("sin(5*x[0])") | ||
Line 154: | Line 113: | ||
# Plot solution | # Plot solution | ||
plot(u, interactive=True) | plot(u, interactive=True) | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
</div> | </div> | ||
Line 161: | Line 118: | ||
<div style="clear:both"></div> | <div style="clear:both"></div> | ||
=== Mixed Formulation for Poisson Equation === | |||
<div style="width: 100%;"> | <div style="width: 100%;"> | ||
<div style="float:left; width: 48%"> | <div style="float:left; width: 48%"> | ||
{{Code|Define MixedPoisson problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | {{Code|Define MixedPoisson problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
pkg load fem-fenics msh | pkg load fem-fenics msh | ||
import_ufl_Problem ('MixedPoisson') | |||
# Create mesh | # Create mesh | ||
Line 198: | Line 132: | ||
mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4)); | mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4)); | ||
# ufl | # File MixedPoisson.ufl | ||
# BDM = FiniteElement("BDM", triangle, 1) | # BDM = FiniteElement("BDM", triangle, 1) | ||
# DG = FiniteElement("DG", triangle, 0) | # DG = FiniteElement("DG", triangle, 0) | ||
Line 205: | Line 139: | ||
# Define trial and test function | # Define trial and test function | ||
# ufl | # File MixedPoisson.ufl | ||
# (sigma, u) = TrialFunctions(W) | # (sigma, u) = TrialFunctions(W) | ||
# (tau, v) = TestFunctions(W) | # (tau, v) = TestFunctions(W) | ||
Line 214: | Line 148: | ||
# Define variational form | # Define variational form | ||
# ufl | # File MixedPoisson.ufl | ||
# a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx | # a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx | ||
# L = - f*v*dx | # L = - f*v*dx | ||
a = BilinearForm ('MixedPoisson' | a = BilinearForm ('MixedPoisson', V); | ||
L = LinearForm ('MixedPoisson', V, f); | L = LinearForm ('MixedPoisson', V, f); | ||
Line 260: | Line 194: | ||
{{Code|Define MixedPoisson problem with fenics python|<syntaxhighlight lang="python" style="font-size:13px"> | {{Code|Define MixedPoisson problem with fenics python|<syntaxhighlight lang="python" style="font-size:13px"> | ||
from dolfin import * | from dolfin import * | ||
# Create mesh | # Create mesh | ||
mesh = UnitSquareMesh(32, 32) | mesh = UnitSquareMesh(32, 32) | ||
Line 291: | Line 209: | ||
(sigma, u) = TrialFunctions(W) | (sigma, u) = TrialFunctions(W) | ||
(tau, v) = TestFunctions(W) | (tau, v) = TestFunctions(W) | ||
Line 337: | Line 256: | ||
interactive() | interactive() | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
</div> | </div> | ||
Line 343: | Line 261: | ||
<div style="clear:both"></div> | <div style="clear:both"></div> | ||
=== | === Mixed Formulation for Poisson Equation === | ||
This time we compare the code with the | This time we compare the code with the c++ version of DOLFIN. | ||
{{Code|HyperElasticity Problem: the ufl file|<syntaxhighlight lang="octave" style="font-size:13px"> | {{Code|HyperElasticity Problem: the ufl file|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
# Function spaces | # Function spaces | ||
Line 393: | Line 301: | ||
J = derivative(F, u, du) | J = derivative(F, u, du) | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
Line 402: | Line 309: | ||
{{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 | pkg load fem-fenics msh | ||
problem = 'HyperElasticity'; | |||
import_ufl_Problem (problem); | |||
Line 530: | Line 398: | ||
# Define source and boundary traction functions | # Define source and boundary traction functions | ||
B = | B = Expression ('B', @(x,y,z) [0.0; -0.5; 0.0]); | ||
T = | T = Expression ('T', @(x,y,z) [0.1; 0.0; 0.0]); | ||
Line 539: | Line 407: | ||
E = 10.0; | E = 10.0; | ||
nu = 0.3; | nu = 0.3; | ||
mu = | mu = Expression ('mu', @(x,y,z) E./(2*(1+nu))); | ||
lmbda = | lmbda = Expression ('lmbda', @(x,y,z) E*nu./((1+nu)*(1-2*nu))); | ||
u = Expression ('u', @(x,y,z) [0; 0; 0]); | u = Expression ('u', @(x,y,z) [0; 0; 0]); | ||
Line 554: | Line 422: | ||
u0 = assemble (L, bcs{:}); | u0 = assemble (L, bcs{:}); | ||
# Create function for the resolution of the NL problem | # Create function for the resolution of the NL problem | ||
function [y, jac] = f (problem, xx, V, bc1, bc2, B, T, mu, lmbda) | # 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 | # endfunction | ||
fs = @(xx) f (problem, xx, V, bcl, bcr, B, T, mu, lmbda); | fs = @(xx) f (problem, xx, V, bcl, bcr, B, T, mu, lmbda); | ||
[x, fval, info] = fsolve (fs, u0, optimset ("jacobian", "on")); | [x, fval, info] = fsolve (fs, u0, optimset ("jacobian", "on")); | ||
Line 586: | Line 453: | ||
using namespace dolfin; | using namespace dolfin; | ||
// Sub domain for clamp at left end | // Sub domain for clamp at left end | ||
Line 679: | Line 507: | ||
// New coordinates | // New coordinates | ||
double y = y0 + (x[1]-y0)*cos(theta) - (x[2]-z0)*sin(theta); | 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); | double z = z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta); | ||
// Rotate at right end | // Rotate at right end | ||
Line 734: | Line 562: | ||
// Solve nonlinear variational problem F(u; v) = 0 | // Solve nonlinear variational problem F(u; v) = 0 | ||
solve(F == 0, u, bcs, J); | solve(F == 0, u, bcs, J); | ||
Line 760: | Line 587: | ||
return 0; | return 0; | ||
} | } | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
</div> | </div> | ||
</div> | </div> | ||
<div style="clear:both"></div> | <div style="clear:both"></div> | ||
[[Category:OctaveForge]] | |||
[[Category:Packages]] | |||
[[ | |||
[[Category: |