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 | ||
<math>-\mathrm{div}\ ( \nabla u(x, y) ) ) = 1 \qquad \mbox{ in } \Omega</math> | <math> -\mathrm{div}\ ( \nabla u(x, y) ) ) = 1 \qquad \mbox{ in } \Omega</math> | ||
<math>u(x, y) = 0 \qquad \mbox{ on } \Gamma_D</math> | <math> u(x, y) = 0 \qquad \mbox{ on } \Gamma_D </math> | ||
<math>(\nabla u(x, y) ) \cdot \mathbf{n} = 0 \qquad \mbox{ on } \Gamma_N</math> | <math> ( \nabla u(x, y) ) \cdot \mathbf{n} = 0 \qquad \mbox{ on } \Gamma_N</math> | ||
<div style="width: 100%;"> | <div style="width: 100%;"> | ||
Line 33: | Line 16: | ||
{{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 22: | ||
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 35: | ||
# Define variational problem | # Define variational problem | ||
# | # File Poisson.ufl | ||
# u = TrialFunction(element) | # u = TrialFunction(element) | ||
# v = TestFunction(element) | # v = TestFunction(element) | ||
Line 71: | Line 45: | ||
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 73: | ||
# Create mesh and define function space | |||
mesh = UnitSquareMesh(32, 32) | |||
V = FunctionSpace(mesh, "Lagrange", 1) | V = FunctionSpace(mesh, "Lagrange", 1) | ||
Line 131: | Line 96: | ||
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 119: | ||
# Plot solution | # Plot solution | ||
plot(u, interactive=True) | plot(u, interactive=True) | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
</div> | </div> | ||
Line 161: | Line 124: | ||
<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 138: | ||
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 145: | ||
# 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 154: | ||
# 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 200: | ||
{{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 215: | ||
(sigma, u) = TrialFunctions(W) | (sigma, u) = TrialFunctions(W) | ||
(tau, v) = TestFunctions(W) | (tau, v) = TestFunctions(W) | ||
Line 337: | Line 262: | ||
interactive() | interactive() | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
</div> | </div> | ||
Line 344: | Line 268: | ||
=== Hyperelasticity === | === Hyperelasticity === | ||
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 307: | ||
J = derivative(F, u, du) | J = derivative(F, u, du) | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
Line 402: | Line 315: | ||
{{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 554: | Line 428: | ||
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 459: | ||
using namespace dolfin; | using namespace dolfin; | ||
// Sub domain for clamp at left end | // Sub domain for clamp at left end | ||
Line 679: | Line 513: | ||
// 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 568: | ||
// 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 593: | ||
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: |