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 8: | Line 8: | ||
== Tutorials == | == Tutorials == | ||
A generic problem has to be solved in two steps: | |||
* | * a file where the '''abstract problem''' is described: this file has to be written in Unified Form Language ('''UFL'''), which is ''a domain specific language for defining discrete variational forms and functionals in a notation '''close to pen-and-paper formulation'''.'' UFL is easy to learn, and in any case the User manual provides explanations and examples. [http://fenicsproject.org/documentation/ufl/1.2.0/user/user_manual.html#ufl-user-manual] | ||
* | * a script file ('''.m''') where the abstract problem is imported and a '''specific problem''' is implemented and solved: this is the script file where the fem-fenics functions are used. Their '''syntax is as close as possible to the python interface''', so that Fenics users should be comfortable with it, but it is also quite intuitive for beginners. The examples below show the equivalence between the different programming languages. | ||
Line 33: | Line 33: | ||
{{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 39: | ||
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 52: | ||
# Define variational problem | # Define variational problem | ||
# | # File Poisson.ufl | ||
# u = TrialFunction(element) | # u = TrialFunction(element) | ||
# v = TestFunction(element) | # v = TestFunction(element) | ||
Line 71: | Line 62: | ||
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 | ||
Line 99: | Line 90: | ||
# Create mesh and define function space | |||
mesh = UnitSquareMesh(32, 32) | |||
V = FunctionSpace(mesh, "Lagrange", 1) | V = FunctionSpace(mesh, "Lagrange", 1) | ||
Line 162: | Line 144: | ||
=== Mixed Formulation for Poisson Equation === | === Mixed Formulation for Poisson Equation === | ||
In this example the Poisson equation is solved with a '''mixed approach''': the stable FE space obtained using Brezzi-Douglas-Marini polynomial of order 1 and | In this example the Poisson equation is solved with a '''mixed approach''': the stable FE space obtained using Brezzi-Douglas-Marini polynomial of order 1 and Dicontinuos element of order 0 is used. | ||
<math>-\mathrm{div}\ ( \mathbf{\sigma} (x, y) ) ) = f (x, y) \qquad \mbox{ in } \Omega</math> | <math>-\mathrm{div}\ ( \mathbf{\sigma} (x, y) ) ) = f (x, y) \qquad \mbox{ in } \Omega</math> | ||
Line 172: | Line 154: | ||
<math>(\sigma (x, y) ) \cdot \mathbf{n} = \sin (5x) \qquad \mbox{ on } \Gamma_N</math> | <math>(\sigma (x, y) ) \cdot \mathbf{n} = \sin (5x) \qquad \mbox{ on } \Gamma_N</math> | ||
A complete description of the problem is | A complete description of the problem is avilable on the [http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/mixed-poisson/python/documentation.html Fenics website.] | ||
<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 165: | ||
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 172: | ||
# 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 181: | ||
# 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 | ||
Line 260: | Line 227: | ||
{{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 242: | ||
(sigma, u) = TrialFunctions(W) | (sigma, u) = TrialFunctions(W) | ||
(tau, v) = TestFunctions(W) | (tau, v) = TestFunctions(W) | ||
Line 344: | Line 296: | ||
=== Hyperelasticity === | === Hyperelasticity === | ||
This time we compare the code with the | This time we compare the code with the c++ version of DOLFIN. The problem for an elastic material can be expressed as a minimization problem | ||
<math> \min_{u \in V} \Pi</math> | <math> \min_{u \in V} \Pi</math> | ||
Line 402: | Line 354: | ||
{{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 586: | Line 499: | ||
using namespace dolfin; | using namespace dolfin; | ||
// Sub domain for clamp at left end | // Sub domain for clamp at left end | ||
Line 734: | Line 608: | ||
// 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 1,059: | Line 932: | ||
# Hold plot | # Hold plot | ||
interactive() | interactive() | ||
© Copyright 2011, The FEniCS Project | © Copyright 2011, The FEniCS Project | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
Line 1,303: | Line 942: | ||
The relevant implementation details which the user should know are: | The relevant implementation details which the user should know are: | ||
* | *all the objects are managed using boost::shared_ptr <>. It means that '''the same resource can be shared by more objects''' and useless copies should be avoided. For example, if we have two different functional spaces in the same problem, like with Navier-Stokes for the velocity and the pressure, the mesh is shared between them and no one has its own copy. | ||
*The '''essential BC are imposed directly to the matrix''' with the command '''assemble()''', which sets to zero all the off diagonal elements in the corresponding line, sets to 1 the diagonal element and sets to the exact value the rhs. This means that we could loose the symmetry of the matrix, if any. To avoid this problem and preserve the symmetry of the system it is available the '''assemble_system()''' command which builds at once the lhs and the rhs. | |||
* | |||
* | *The '''coefficient of the variational problem''' can be specified using either an ''Expression()'', a ''Constant()'' or a ''Function()''. They are different objects which behave in different ways: an ''Expession'' or a ''Constant'' object overloads the eval() method of the dolfin::Expression class and it is evaluated at run time using the octave function feval (). A ''Function'' instead doesn't need to be evaluated because it is assembled copying element-by-element the values contained in the input vector. | ||
== Additional functionality / TODOS == | == Additional functionality / TODOS == | ||
Line 1,330: | Line 953: | ||
* '''Norma'''l: add the possibility to use a reserved keyword (normal ?) to be used with the DirichletBC. <syntaxhighlight lang="octave" style="font-size:13px"> bc = DirichletBC (V, @(x, y, normal) [sin(x)*normal; 0], [3, 4]);</syntaxhighlight> | * '''Norma'''l: add the possibility to use a reserved keyword (normal ?) to be used with the DirichletBC. <syntaxhighlight lang="octave" style="font-size:13px"> bc = DirichletBC (V, @(x, y, normal) [sin(x)*normal; 0], [3, 4]);</syntaxhighlight> | ||
* | * @function/'''feval''': the function should accept as input also an array of values. Show how it can be used in an example with odepkg. | ||
== External Links == | == External Links == | ||
Line 1,338: | Line 961: | ||
* Presentation at MOX [https://drive.google.com/file/d/0ByWLfuWVSWHbZWZzRzY2em5PU28/edit?usp=sharing] | * Presentation at MOX [https://drive.google.com/file/d/0ByWLfuWVSWHbZWZzRzY2em5PU28/edit?usp=sharing] | ||
[[Category:Octave Forge]] | [[Category:Octave-Forge]] |