Table of Contents

*SfePy*can be used in two basic ways:- a black-box partial differential equation (PDE) solver,
- a Python package to build custom applications involving solving PDEs by the finite element (FE) method.

This tutorial focuses on the first way and introduces the basic concepts
and nomenclature used in the following parts of the documentation. Check
also the *Primer* which focuses on a particular problem in detail.

Users not familiar with the finite element method should start with the
*Notes on solving PDEs by the Finite Element Method*.

The simplest way of using *SfePy* is to solve a system of PDEs defined
in a **problem description file**, also referred to as **input
file**. In such a file, the problem is described using several keywords
that allow one to define the equations, variables, finite element
approximations, solvers, solution domain and subdomains etc., see
*Problem Description File* for a full list of those keywords.

The syntax of the problem description file is very simple yet powerful, as the file itself is just a regular Python module that can be normally imported - no special parsing is necessary. The keywords mentioned above are regular Python variables (usually of the dict type) with special names.

Below we show:

- how to solve a problem given by a problem description file, and
- explain the elements of the file on several examples.

But let us begin with a slight detour...

- A top-level script (usually
`simple.py`, as in this tutorial) reads in an input file. - Following the contents of the input file, a
`Problem`instance is created - this is the input file coming to life. Let us call the instance`problem`.- The
`problem`sets up its domain, regions (various sub-domains), fields (the FE approximations), the equations and the solvers. The equations determine the materials and variables in use - only those are fully instantiated, so the input file can safely contain definitions of items that are not used actually.

- The
- Prior to solution,
`problem.time_update()`function has to be called to setup boundary conditions, material parameters and other potentially time-dependent data. This holds also for stationary problems with a single “time step”. - The solution is then obtained by calling
`problem.solve()`function. - Finally, the solution can be stored using
`problem.save_state()`

The above last three steps are essentially repeated for each time
step. So that is it - using the code a black-box PDE solver shields the
user from having to create the `Problem` instance by hand. But note
that this is possible, and often necessary when the flexibility of the
default solvers is not enough. At the end of the tutorial an example
demonstrating the interactive creation of the `problem` is shown, see
*Interactive Example: Linear Elasticity*.

Now let us continue with running a simulation.

The following commands should be run in the top-level directory of the *SfePy*
source tree after compiling the C extension files. See
*Installation* for full installation instructions.

This section introduces the basics of running *SfePy* on the command line. The
`$` indicates the command prompt of your terminal.

The script

`simple.py`is the most basic starting point in*SfePy*. It is invoked as follows:$ ./simple.py examples/diffusion/poisson_short_syntax.py

`examples/diffusion/poisson_short_syntax.py`is the*SfePy**problem description*file, which defines the problem to be solved in terms*SfePy*can understand- Running the above command creates the output file
`cylinder.vtk`in the*SfePy*top-level directory

*SfePy*can also be invoked interactively using IPython with custom imports, as described in*Using IPython*. In the SfePy top-level directory, run:$ ipython --profile=sfepy

See

*Interactive Example: Linear Elasticity*for more information.

The

`postproc.py`script can be used for quick postprocessing and visualization of the*SfePy*output files. It requires mayavi2 installed on your system.As a simple example, try:

$ ./postproc.py cylinder.vtk

The following interactive 3D window should display:

- The left mouse button by itself orbits the 3D view
- Holding shift and the left mouse button pans the view
- Holding control and the left mouse button rotates about the screen normal axis
- The right mouse button controls the zoom

Here we discuss the contents of the
`examples/diffusion/poisson_short_syntax.py` problem description file. For
additional examples, see the problem description files in the
`examples/` directory of *SfePy*.

The problem at hand is the following:

(1)

where is a subset of the domain boundary. For simplicity, we set , but we still work with the material constant even though it has no influence on the solution in this case. We also assume zero fluxes over , i.e. there. The particular boundary conditions used below are on the left side of the cylindrical domain depicted in the previous section and on the right side.

The first step to do is to write (1) in *weak
formulation* (?). The ,
. So only one term in weak form
(?) remains:

(2)

Comparing the above integral term with the
long table in *Term Overview*, we can see that *SfePy* contains
this term under name dw_laplace. We are now ready to proceed to the
actual problem definition.

Open the `examples/diffusion/poisson_short_syntax.py` file in your favorite text
editor. Note that the file is a regular python source code.

```
from sfepy import data_dir
filename_mesh = data_dir + '/meshes/3d/cylinder.mesh'
```

The `filename_mesh` variable points to the file containing the mesh for the
particular problem. *SfePy* supports a variety of mesh formats.

```
materials = {
'coef': ({'val' : 1.0},),
}
```

Here we define just a constant coefficient of the Poisson
equation, using the `'values'` attribute. Other possible attribute is
`'function'`, for material coefficients computed/obtained at runtime.

Many finite element problems require the definition of material parameters.
These can be handled in *SfePy* with material variables which associate the
material parameters with the corresponding region of the mesh.

```
regions = {
'Omega' : 'all', # or 'cells of group 6'
'Gamma_Left' : ('vertices in (x < 0.00001)', 'facet'),
'Gamma_Right' : ('vertices in (x > 0.099999)', 'facet'),
}
```

Regions assign names to various parts of the finite element mesh. The region names can later be referred to, for example when specifying portions of the mesh to apply boundary conditions to. Regions can be specified in a variety of ways, including by element or by node. Here, Omega is the elemental domain over which the PDE is solved and Gamma_Left and Gamma_Right define surfaces upon which the boundary conditions will be applied.

```
fields = {
'temperature': ('real', 1, 'Omega', 1)
}
```

A field is used mainly to define the approximation on a (sub)domain, i.e. to define the discrete spaces , where we seek the solution.

The Poisson equation can be used to compute e.g. a temperature distribution,
so let us call our field `'temperature'`. On the region `'Omega'`
it will be approximated using linear finite elements.

A field in a given region defines the finite element approximation. Several variables can use the same field, see below.

```
variables = {
't': ('unknown field', 'temperature', 0),
's': ('test field', 'temperature', 't'),
}
```

One field can be used to generate discrete degrees of freedom (DOFs) of
several variables. Here the unknown variable (the temperature) is called
`'t'`, it’s associated DOF name is `'t.0'` — this will be referred
to in the Dirichlet boundary section (`ebc`). The corresponding test
variable of the weak formulation is called `'s'`. Notice that the
`'dual'` item of a test variable must specify the unknown it
corresponds to.

For each unknown (or state) variable there has to be a test (or virtual) variable defined, as usual in weak formulation of PDEs.

```
ebcs = {
't1': ('Gamma_Left', {'t.0' : 2.0}),
't2', ('Gamma_Right', {'t.0' : -2.0}),
}
```

Essential (Dirichlet) boundary conditions can be specified as above.

Boundary conditions place restrictions on the finite element formulation and create a unique solution to the problem. Here, we specify that a temperature of +2 is applied to the left surface of the mesh and a temperature of -2 is applied to the right surface.

```
integrals = {
'i': 2,
}
```

Integrals specify which numerical scheme to use. Here we are using a 2nd order quadrature over a 3 dimensional space.

```
equations = {
'Temperature' : """dw_laplace.i.Omega( coef.val, s, t ) = 0"""
}
```

The equation above directly corresponds to the discrete version of (2), namely: Find , such that

where .

The equations block is the heart of the *SfePy* problem definition file. Here,
we are specifying that the Laplacian of the temperature (in the weak
formulation) is 0, where `coef.val` is a material constant. We are using the
`i` integral defined previously, over the domain specified by the region
Omega.

The above syntax is useful for defining *custom integrals* with
user-defined quadrature points and weights, see *Integrals*. The
above uniform integration can be more easily achieved by:

```
equations = {
'Temperature' : """dw_laplace.2.Omega( coef.val, s, t ) = 0"""
}
```

The integration order is specified directly in place of the integral name. The integral definition is superfluous in this case.

```
solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton',
{'i_max' : 1,
'eps_a' : 1e-10,
}),
}
```

Here, we specify the linear and nonlinear solver kind and options. The convergence parameters can be adjusted if necessary, otherwise leave the default.

Even linear problems are solved by a nonlinear solver (KISS rule) - only one iteration is needed and the final rezidual is obtained for free.

```
options = {
'nls' : 'newton',
'ls' : 'ls',
}
```

The solvers to use are specified in the options block. We can define multiple solvers with different convergence parameters if necessary.

That’s it! Now it is possible to proceed as described in
*Invoking SfePy from the command line*.

This example shows how to use *SfePy* interactively, but also how to make a
custom simulation script. We will use IPython with custom imports, as
described in *Using IPython*, for the explanation, but regular Python
shell would do as well, provided the proper modules are imported.

We wish to solve the following linear elasticity problem:

(3)

where the stress is defined as , , are the Lamé’s constants, the strain is and are volume forces. This can be written in general form as , where in our case .

In the weak form the equation (3) is

(4)

where is the test function, and both , belong to a suitable function space.

**Hint:** Whenever you create a new object (e.g. a Mesh instance, see
below), try to print it using the print statement - it will give you
insight about the object internals.

The whole example summarized in a script is below in
*Complete Example as a Script*.

In the SfePy top-level directory, run:

```
$ ipython --profile=sfepy
```

Read a finite element mesh, that defines the domain .

```
In [1]: mesh = Mesh.from_file('meshes/2d/rectangle_tri.mesh')
```

Create a domain. The domain allows defining regions or subdomains.

```
In [2]: domain = FEDomain('domain', mesh)
```

Define the regions - the whole domain , where the solution is sought, and , , where the boundary conditions will be applied. As the domain is rectangular, we first get a bounding box to get correct bounds for selecting the boundary edges.

```
In [3]: min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
In [4]: eps = 1e-8 * (max_x - min_x)
In [5]: omega = domain.create_region('Omega', 'all')
In [6]: gamma1 = domain.create_region('Gamma1',
...: 'vertices in x < %.10f' % (min_x + eps),
...: 'facet')
In [7]: gamma2 = domain.create_region('Gamma2',
...: 'vertices in x > %.10f' % (max_x - eps),
...: 'facet')
```

Next we define the actual finite element approximation using the
`Field` class.

```
In [8]: field = Field.from_args('fu', nm.float64, 'vector', omega,
...: space='H1', poly_space_base='lagrange',
...: approx_order=2)
```

Using the field fu, we can define both the unknown variable and the test variable .

```
In [9]: u = FieldVariable('u', 'unknown', field)
In [10]: v = FieldVariable('v', 'test', field, primary_var_name='u')
```

Before we can define the terms to build the equation of linear elasticity, we have to create also the materials, i.e. define the (constitutive) parameters. The linear elastic material m will be defined using the two Lamé constants , . The volume forces will be defined also as a material, as a constant (column) vector .

```
In [11]: m = Material('m', lam=1.0, mu=1.0)
In [12]: f = Material('f', val=[[0.02], [0.01]])
```

One more thing needs to be defined - the numerical quadrature that will be used to integrate each term over its domain.

```
In [14]: integral = Integral('i', order=3)
```

Now we are ready to define the two terms and build the equations.

```
In [15]: from sfepy.terms import Term
In [16]: t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u)',
integral, omega, m=m, v=v, u=u)
In [17]: t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
In [18]: eq = Equation('balance', t1 + t2)
In [19]: eqs = Equations([eq])
```

The equations have to be completed by boundary conditions. Let us clamp the left edge , and shift the right edge in the direction a bit, depending on the coordinate.

```
In [20]: from sfepy.discrete.conditions import Conditions, EssentialBC
In [21]: fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})
In [22]: def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0):
....: val = shift * coors[:,1]**2
....: return val
In [23]: bc_fun = Function('shift_u_fun', shift_u_fun,
....: extra_args={'shift' : 0.01})
In [24]: shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})
```

The last thing to define before building the problem are the solvers. Here we just use a sparse direct SciPy solver and the SfePy Newton solver with default parameters. We also wish to store the convergence statistics of the Newton solver. As the problem is linear, it should converge in one iteration.

```
In [25]: from sfepy.solvers.ls import ScipyDirect
In [26]: from sfepy.solvers.nls import Newton
In [27]: ls = ScipyDirect({})
In [28]: nls_status = IndexedStruct()
In [29]: nls = Newton({}, lin_solver=ls, status=nls_status)
```

Now we are ready to create a `Problem` instance. Note
that the step above is not really necessary - the above solvers are
constructed by default. We did them to get the nls_status.

```
In [30]: pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls)
```

The `Problem` has several handy methods for
debugging. Let us try saving the regions into a VTK file.

```
In [31]: pb.save_regions_as_groups('regions')
```

And view them.

```
In [32]: view = Viewer('regions.vtk')
In [33]: view()
```

You should see this:

Finally, we apply the boundary conditions, solve the problem, save and view the results.

```
In [34]: pb.time_update(ebcs=Conditions([fix_u, shift_u]))
In [35]: vec = pb.solve()
In [36]: print nls_status
In [37]: pb.save_state('linear_elasticity.vtk', vec)
In [38]: view = Viewer('linear_elasticity.vtk')
In [39]: view()
```

This is the resulting image:

The default view is not very fancy. Let us show the displacements by shifting the mesh. Close the previous window and do:

```
In [56]: view(vector_mode='warp_norm', rel_scaling=2,
....: is_scalar_bar=True, is_wireframe=True)
```

And the result is:

See the docstring of view() and play with its options.

The source code: `linear_elasticity.py`. It should be
run from the *SfePy* source directory so that it finds the mesh file.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 | ```
#!/usr/bin/env python
from __future__ import print_function
from __future__ import absolute_import
from argparse import ArgumentParser
import numpy as nm
import sys
sys.path.append('.')
from sfepy.base.base import IndexedStruct
from sfepy.discrete import (FieldVariable, Material, Integral, Function,
Equation, Equations, Problem)
from sfepy.discrete.fem import Mesh, FEDomain, Field
from sfepy.terms import Term
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.postprocess.viewer import Viewer
from sfepy.mechanics.matcoefs import stiffness_from_lame
def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0):
"""
Define a displacement depending on the y coordinate.
"""
val = shift * coors[:,1]**2
return val
helps = {
'show' : 'show the results figure',
}
def main():
from sfepy import data_dir
parser = ArgumentParser()
parser.add_argument('--version', action='version', version='%(prog)s')
parser.add_argument('-s', '--show',
action="store_true", dest='show',
default=False, help=helps['show'])
options = parser.parse_args()
mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh')
domain = FEDomain('domain', mesh)
min_x, max_x = domain.get_mesh_bounding_box()[:,0]
eps = 1e-8 * (max_x - min_x)
omega = domain.create_region('Omega', 'all')
gamma1 = domain.create_region('Gamma1',
'vertices in x < %.10f' % (min_x + eps),
'facet')
gamma2 = domain.create_region('Gamma2',
'vertices in x > %.10f' % (max_x - eps),
'facet')
field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2)
u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')
m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))
f = Material('f', val=[[0.02], [0.01]])
integral = Integral('i', order=3)
t1 = Term.new('dw_lin_elastic(m.D, v, u)',
integral, omega, m=m, v=v, u=u)
t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
eq = Equation('balance', t1 + t2)
eqs = Equations([eq])
fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})
bc_fun = Function('shift_u_fun', shift_u_fun, extra_args={'shift' : 0.01})
shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})
ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({}, lin_solver=ls, status=nls_status)
pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls)
pb.save_regions_as_groups('regions')
pb.time_update(ebcs=Conditions([fix_u, shift_u]))
vec = pb.solve()
print(nls_status)
pb.save_state('linear_elasticity.vtk', vec)
if options.show:
view = Viewer('linear_elasticity.vtk')
view(vector_mode='warp_norm', rel_scaling=2,
is_scalar_bar=True, is_wireframe=True)
if __name__ == '__main__':
main()
``` |