Fem-fenics: Difference between revisions

Jump to navigation Jump to search

Fem-fenics (edit)

Revision as of 11:12, 10 June 2019

9,077 bytes added ,  10 June 2019
m
Remove redundant Category:Packages.
m (Remove redundant Category:Packages.)
 
(23 intermediate revisions by 6 users not shown)
Line 8: Line 8:
== Tutorials ==
== Tutorials ==


A generic problem has to be solved in two steps:
The solution of a problem can be logically divided in two steps. According to convenience or personal preference, they can be addressed with different files or just in one Octave script:


* 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]
* the description of the '''abstract problem''': this should be done via the 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] As mentioned before, the problem can be defined in a separate .ufl file or handled directly in an m-file using ufl blocks.
* 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.
* the implementation of a '''specific problem''', an instance of the abstract one: this is done in a script file ('''.m''') where the fem-fenics functions are used and the problem is solved. 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')
 
ufl start Poisson
ufl element = FiniteElement("Lagrange", triangle, 1)
ufl
ufl u = TrialFunction(element)
ufl v = TestFunction(element)
ufl f = Coefficient(element)
ufl g = Coefficient(element)
ufl
ufl a = inner(grad(u), grad(v))*dx
ufl L = f*v*dx + g*v*ds
ufl end


# Create mesh and define function space
# Create mesh and define function space
Line 39: Line 50:
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 52: Line 61:


# Define variational problem
# Define variational problem
# File Poisson.ufl
# Operation performed in the ufl snippet above
# u = TrialFunction(element)
# u = TrialFunction(element)
# v = TestFunction(element)
# v = TestFunction(element)
Line 62: Line 71:
g = Expression ('g', @(x,y) sin (5.0 * x));
g = Expression ('g', @(x,y) sin (5.0 * x));


# File Poisson.ufl
# As in the ufl snippet above
# 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 88: Line 97:
{{Code|Define Poisson problem with fenics python|<syntaxhighlight lang="python" style="font-size:13px">  
{{Code|Define Poisson problem with fenics python|<syntaxhighlight lang="python" style="font-size:13px">  
from dolfin import *
from dolfin import *




Line 93: Line 113:


mesh = UnitSquareMesh(32, 32)
mesh = UnitSquareMesh(32, 32)


V = FunctionSpace(mesh, "Lagrange", 1)
V = FunctionSpace(mesh, "Lagrange", 1)
Line 144: Line 162:


=== 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 Dicontinuos element of order 0 is used.
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 Discontinuos elements 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 154: Line 172:
<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 avilable on the [http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/mixed-poisson/python/documentation.html Fenics website.]
A complete description of the problem is available 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')
 
ufl start MixedPoisson
ufl
ufl BDM = FiniteElement("BDM", triangle, 1)
ufl DG  = FiniteElement("DG", triangle, 0)
ufl W = BDM * DG
ufl
ufl "(sigma, u)" = TrialFunctions(W)
ufl "(tau, v)" = TestFunctions(W)
ufl
ufl CG = FiniteElement("CG", triangle, 1)
ufl f = Coefficient(CG)
ufl
ufl a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
ufl L = - f*v*dx
ufl end


# Create mesh
# Create mesh
Line 165: Line 198:
mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));


# File MixedPoisson.ufl
# ufl snippet above
#  BDM = FiniteElement("BDM", triangle, 1)
#  BDM = FiniteElement("BDM", triangle, 1)
#  DG  = FiniteElement("DG", triangle, 0)
#  DG  = FiniteElement("DG", triangle, 0)
Line 172: Line 205:


# Define trial and test function
# Define trial and test function
# File MixedPoisson.ufl
# ufl snippet above
#  (sigma, u) = TrialFunctions(W)
#  (sigma, u) = TrialFunctions(W)
#  (tau, v)  = TestFunctions(W)
#  (tau, v)  = TestFunctions(W)
Line 181: Line 214:


# Define variational form
# Define variational form
# File MixedPoisson.ufl
# ufl snippet above
#  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', V);
a = BilinearForm ('MixedPoisson', V, V);
L = LinearForm ('MixedPoisson', V, f);
L = LinearForm ('MixedPoisson', V, f);


Line 227: Line 260:
{{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 242: Line 291:
(sigma, u) = TrialFunctions(W)
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
(tau, v) = TestFunctions(W)




Line 296: Line 344:


=== Hyperelasticity ===
=== Hyperelasticity ===
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
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 354: Line 402:
{{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);




ufl start HyperElasticity
ufl # Function spaces
ufl element = VectorElement("Lagrange", tetrahedron, 1)
ufl
ufl # Trial and test functions
ufl du = TrialFunction(element)    # Incremental displacement
ufl v  = TestFunction(element)      # Test function
ufl
ufl # Functions
ufl u = Coefficient(element)        # Displacement from previous iteration
ufl B = Coefficient(element)        # Body force per unit volume
ufl T = Coefficient(element)        # Traction force on the boundary
ufl
ufl # Kinematics
ufl I = Identity(element.cell().d)  # Identity tensor
ufl F = I + grad(u)                # Deformation gradient
ufl C = F.T*F                      # Right Cauchy-Green tensor
ufl
ufl # Invariants of deformation tensors
ufl Ic = tr(C)
ufl J  = det(F)
ufl
ufl # Elasticity parameters
ufl mu    = Constant(tetrahedron)
ufl lmbda = Constant(tetrahedron)
ufl
ufl # Stored strain energy density (compressible neo-Hookean model)
ufl psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
ufl
ufl # Total potential energy
ufl Pi = psi*dx - inner(B, u)*dx - inner(T, u)*ds
ufl
ufl # First variation of Pi (directional derivative about u in the direction of v)
ufl F = derivative(Pi, u, v)
ufl
ufl # Compute Jacobian of F
ufl J = derivative(F, u, du)
ufl end
problem = 'HyperElasticity';




Line 499: Line 586:


using namespace dolfin;
using namespace dolfin;


// Sub domain for clamp at left end
// Sub domain for clamp at left end
Line 608: Line 734:
   // 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 719: Line 846:
#  eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx \
#  eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx \
#      + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
#      + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
a1 = BilinearForm ('TentativeVelocity', V, k);
a1 = BilinearForm ('TentativeVelocity', V, V, k);




Line 725: Line 852:
#  a = inner(grad(p), grad(q))*dx
#  a = inner(grad(p), grad(q))*dx
#  L = -(1/k)*div(u1)*q*dx
#  L = -(1/k)*div(u1)*q*dx
a2 = BilinearForm ('PressureUpdate', Q);
a2 = BilinearForm ('PressureUpdate', Q, Q);


# Velocity update
# Velocity update
#  a = inner(u, v)*dx
#  a = inner(u, v)*dx
#  L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
#  L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
a3 = BilinearForm ('VelocityUpdate', V);
a3 = BilinearForm ('VelocityUpdate', V, V);


# Assemble matrices
# Assemble matrices
Line 932: Line 1,059:
# Hold plot
# Hold plot
interactive()
interactive()
© Copyright 2011, The FEniCS Project
</syntaxhighlight>}}
  </div>
</div>
<div style="clear:both"></div>
=== Obstacles in the Domain ===
[[File:Fem_fenics_Subdomains.png|right|500px]]
<math> - \mathrm {div} (a \nabla u) = 1 \quad \mbox { in } \Omega, </math>
<math> u = 5 \quad \mbox { on } \Gamma_T, </math>
<math> u = 0 \quad \mbox { on } \Gamma_B, </math>
<math> \nabla u \cdot \mathbf {n} = - 10 e^{- (y - 0.5)^2} \quad \mbox { on } \Gamma_L, </math>
<math> \nabla u \cdot \mathbf {n} = 1 \quad \mbox { on } \Gamma_R </math>
The above example is a weighted Poisson problem on the unit square. The diffusion coefficient <math> a </math> assumes value 0.01 on the obstacle <math> \Omega_1 = [0.5, 0.7] \times [0.2, 1.0] </math>, whilst 1 outside on <math> \Omega_0 = \Omega \setminus \Omega_1 </math>. You can find a detailed explanation on the [http://fenicsproject.org/documentation/dolfin/1.4.0/python/demo/documented/subdomains-poisson/python/documentation.html FEniCS website]. On the side, you can see the solution.
<div style="clear: both"></div>
<div style="width: 100%;">
  <div style="float:left; width: 48%">
{{Code|Define Poisson problem with obstacle with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px">
pkg load fem-fenics msh
# Create mesh
x = y = linspace (0, 1, 65);
[msh, facets] = Mesh (msh2m_structured_mesh (x, y, 0, 4:-1:1));
ufl start Subdomains
ufl fe = FiniteElement "(""CG"", triangle, 2)"
ufl u = TrialFunction (fe)
ufl v = TestFunction (fe)
ufl
ufl a0 = Coefficient (fe)
ufl a1 = Coefficient (fe)
ufl g_L = Coefficient (fe)
ufl g_R = Coefficient (fe)
ufl f = Coefficient (fe)
ufl
ufl a= "inner(a0*grad(u),grad(v))*dx(0) + inner(a1*grad(u),grad(v))*dx(1)"
ufl L = g_L*v*ds(1) + g_R*v*ds(3) + f*v*dx(0) + f*v*dx(1)
ufl end
V = FunctionSpace ("Subdomains", msh);
# Define problem coefficients
a0 = Constant ("a0", 1.0);
a1 = Constant ("a1", 0.01);
g_L = Expression ("g_L", @(x, y) - 10*exp(- (y - 0.5) ^ 2));
g_R = Constant ("g_R", 1.0);
f = Constant ("f", 1.0);
# Define subdomains
domains = MeshFunction ("dx", msh, 2, 0);
obstacle = SubDomain (@(x,y) (y >= 0.5) && (y <= 0.7) && ...
                      (x >= 0.2) && (x <= 1.0), false);
domains = mark (obstacle, domains, 1);
# Define boundary conditions
bc1 = DirichletBC (V, @(x, y) 5.0, facets, 2);
bc2 = DirichletBC (V, @(x, y) 0.0, facets, 4);
# Define variational form
a = BilinearForm ("Subdomains", V, V, a0, a1, domains);
L = LinearForm ("Subdomains", V, g_L, g_R, f, facets, domains);
# Solve problem
[A, b] = assemble_system (a, L, bc1, bc2);
sol = A \ b;
u = Function ("u", V, sol);
# Plot solution
[X, Y] = meshgrid (x, y);
U = u (X, Y);
surf (X, Y, U);
</syntaxhighlight>}}
  </div>
  <div style="float:right; width: 48%">
{{Code|Define Poisson problem with obstacle with FEniCS python |<syntaxhighlight lang="python" style="font-size:13px">
from dolfin import *
# Define mesh
mesh = UnitSquareMesh(64, 64)
# Create classes for defining parts of the boundaries
class Left(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 0.0)
class Right(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 1.0)
class Bottom(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 0.0)
class Top(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 1.0)
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)
# Define input data
a0 = Constant(1.0)
a1 = Constant(0.01)
g_L = Expression("- 10*exp(- pow(x[1] - 0.5, 2))")
g_R = Constant("1.0")
f = Constant(1.0)
# Initialize mesh function for interior domains
class Obstacle(SubDomain):
  def inside(self, x, on_boundary):
    return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))
obstacle = Obstacle()
domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)
# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, 5.0, boundaries, 2),
      DirichletBC(V, 0.0, boundaries, 4)]
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]
# Define variational form
F = (inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)
    - g_L*v*ds(1) - g_R*v*ds(3)
    - f*v*dx(0) - f*v*dx(1))
# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)
# Solve problem
u = Function(V)
solve(a == L, u, bcs)
# Plot solution
plot(u, title="u")
interactive()
© Copyright 2011, The FEniCS Project
© Copyright 2011, The FEniCS Project
</syntaxhighlight>}}
</syntaxhighlight>}}
Line 942: Line 1,303:
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.
* 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 ''Expression'' 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.


*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.
* Unfortunately the feature used in previous versions of the package to handle ''DirichletBC'' is not implemented yet in the FEniCS library for distributed execution. This means that prior to running code via MPI it should be adapted to use the newly introduced ''MeshFunction'' to mark border subsets, so that essential boundary conditions are correctly applied.


*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.
== Known issues ==
* Fem-fenics needs both in the installation phase and in normal usage some includes that are not in standard system include directories, namely those related to MPI and the Eigen template library. If you experience an installation failure or errors when importing your UFL files, then you should properly set your CPPFLAGS environment variable before launching Octave, as per [http://octave.1599824.n4.nabble.com/fem-fenics-0-0-1-released-tp4662182p4662195.html]. This can be done in {{Path|sh}} compatible shells with the command below:
<syntaxhighlight lang="bash" style="font-size:13px">
export CPPFLAGS="-I/usr/include/eigen3 $(mpicxx -showme:compile)"
</syntaxhighlight>
Notice that {{Path|/usr/include/eigen3}} has to be changed accordingly to the path in which Eigen are to be found on your system.
  It should readily work in Ubuntu 13.10 and 14.04.
 
* There is a known bug with Openmpi, discussed on the maintainers list [http://octave.1599824.n4.nabble.com/fem-fenics-0-0-1-released-tp4662182p4662234.html], that needs the preloading of a MPI shared library with the following command:
<syntaxhighlight lang="bash" style="font-size:13px">
export LD_PRELOAD=/usr/lib/openmpi/lib/libmpi.so
</syntaxhighlight>
Again, notice that {{Path|/usr/lib/openmpi/lib/libmpi.so}} has to be changed to fit to your system.


== Additional functionality / TODOS ==
== Additional functionality / TODOS ==
Line 953: Line 1,330:
* '''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.
* <strike>  @function/'''feval''': the function should accept as input also an array of values. </strike> Show how it can be used in an example with odepkg.
 
== External Links ==
* User manual [https://drive.google.com/file/d/0ByWLfuWVSWHbaXN3T3diaXEwU0k/edit?usp=sharing].
* Repository [http://sourceforge.net/p/octave/fem-fenics/]
* Functions reference [http://octave.sourceforge.net/fem-fenics/overview.html]
* Presentation at MOX [https://drive.google.com/file/d/0ByWLfuWVSWHbZWZzRzY2em5PU28/edit?usp=sharing]


[[Category:OctaveForge]]
[[Category:Octave Forge]]
[[Category:Packages]]
Retrieved from "https://wiki.octave.org/Special:MobileDiff/3197...11972"

Navigation menu