SfePy Documentation

Release version: 2024.2

Robert Cimrman and Contributors

Jun 28, 2024

 CONTENTS

1 Documentation 3
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Supported Platforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2 Installing SfePy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.3 Using SfePy Docker Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.4 Installing SfePy from Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.5 Testing Installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.6 Debugging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.7 Using IPython . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.8 Notes on Multi-platform Python Distributions . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Tutorial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1 Basic SfePy Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.2 Basic Notions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.3 Running a Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3.4 Example Problem Description File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.5 Interactive Example: Linear Elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4 User’s Guide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.4.1 Running a Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.4.2 Visualization of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.4.3 Problem Description File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.4.4 Building Equations in SfePy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
1.4.5 Term Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
1.4.6 Solution Postprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
1.4.7 Probing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
1.4.8 Postprocessing filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
1.4.9 Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
1.4.10 Solving Problems in Parallel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
1.4.11 Isogeometric Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
1.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
1.5.1 Primer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
1.5.2 Using Salome with SfePy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
1.5.3 Preprocessing: FreeCAD/OpenSCAD + Gmsh . . . . . . . . . . . . . . . . . . . . . . . . . 77
1.5.4 Material Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
1.5.5 Mesh parametrization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
1.5.6 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
1.5.7 Example Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 708
1.6 Useful Code Snippets and FAQ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 708
1.6.1 Miscellaneous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 709
1.6.2 Mesh-Related Tasks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 711

i
 1.6.3 Regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 711
1.6.4 Material Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 712
1.7 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 713
1.7.1 Notes on solving PDEs by the Finite Element Method . . . . . . . . . . . . . . . . . . . . . 713
1.7.2 Implementation of Essential Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 715
1.8 Term Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 717
1.8.1 Term Syntax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 717
1.8.2 Term Table . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 719

2 Development 747
2.1 General Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 747
2.2 Possible Topics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 747
2.3 Developer Guide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 748
2.3.1 Retrieving the Latest Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 748
2.3.2 SfePy Directory Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 749
2.3.3 Exploring the Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 750
2.3.4 How to Contribute . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 751
2.3.5 How to Regenerate Documentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756
2.3.6 How to Implement a New Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756
2.3.7 Multi-linear Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 762
2.3.8 How To Make a Release . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764
2.3.9 Module Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 767

Bibliography 1213

Python Module Index 1215

Index 1219

ii
 SfePy Documentation, Release version: 2024.2

• genindex
• modindex
• search

CONTENTS 1
SfePy Documentation, Release version: 2024.2

2 CONTENTS
 CHAPTER

ONE

DOCUMENTATION

1.1 Introduction

SfePy (http://sfepy.org) is a software for solving systems of coupled partial differential equations (PDEs) by the finite
element method in 1D, 2D and 3D. It can be viewed both as black-box PDE solver, and as a Python package which can
be used for building custom applications. The word “simple” means that complex FEM problems can be coded very
easily and rapidly.
There is also a preliminary support for the isogeometric analysis, outlined in Isogeometric Analysis.
The code is written almost entirely in Python, with exception of the most time demanding routines - those are written
in C and wrapped by Cython or written directly in Cython.
SfePy is a free software released under the New BSD License. It relies on NumPy and SciPy (an excellent collection
of tools for scientific computations in Python). It is a multi-platform software that should work on Linux, Mac OS X
and Windows.
SfePy was originally developed as a flexible framework to quickly implement and test the mathematical models devel-
oped during our various research projects. It has evolved, however, to a rather full-featured (yet small) finite element
code. Many terms have been implemented that can be used to build the PDEs, see Term Overview. SfePy comes also
with a number of examples that can get you started, check Examples, Gallery and Tutorial. Some more advanced
features are discussed in Primer.

1.2 Installation

1.2.1 Supported Platforms

SfePy is known to work on various flavors of recent Linux, Intel-based MacOS and Windows. The release 2019.4 was
the last with Python 2.7 support. SfePy should work with any recent Python 3.x that is supported by NumPy and SciPy.
Note: Depending on Python installation and OS used, replacing python by python3 might be required in all the
commands below (e.g. in Compilation of C Extension Modules) in order to use Python 3.

3
SfePy Documentation, Release version: 2024.2

Notes on selecting Python Distribution

It is only matter of taste to use either native OS Python installation, pip, or any other suitable distribution. On all
supported platforms we could recommend multi-platform scientific Python distributions Anaconda as easy-to-use,
stable and up-to-date Python distribution with all the required dependencies (including pre-built sfepy package). See
also Notes on Multi-platform Python Distributions for further details.
On different platforms the following options can be recommended:
• Linux: Anaconda, OS native installation, if available, or pip.
• macOS: Anaconda.
• Windows: free versions of commercial scientific Python distributions Anaconda or Enthought Deployment Man-
ager . In addition a completely free open-source portable distribution WinPython can be used.

1.2.2 Installing SfePy

The released versions of SfePy can be installed as follows.

• Using pip:

pip install sfepy

• Using Anaconda: install sfepy from conda-forge channel:

conda install -c conda-forge sfepy

See Notes on Multi-platform Python Distributions for additional notes.

If the installation succeeded, proceed with Testing Installation.

1.2.3 Using SfePy Docker Images

Besides the classical installation we also provide official Docker images with ready-to-run Anaconda and SfePy instal-
lation.
Before you start using SfePy images, you need to first install and configure Docker on your computer. To do this follow
official Docker documentation.
Currently available all-in-one image is:
• sfepy/sfepy-desktop - an Ubuntu based container containing a full desktop environment in officially supported
flavors accessible via any modern web browser.
For available runtime options and further information see sfepy-docker project on Github.

1.2.4 Installing SfePy from Sources

The latest stable release can be obtained from the download page. Otherwise, download the development version of
the code from SfePy git repository:

git clone git://github.com/sfepy/sfepy.git

In case you wish to use a specific release instead of the latest master version, use:

4 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

git tag -l

to see the available releases - the release tags have form release_<year>.<int>.
See the download page for additional download options.

Requirements

Installation prerequisites, required to build SfePy:

• a C compiler suite,
• Python 3.x,
• NumPy,
• Cython.
• Cmake
• scikit-build
• ninja
Python packages required for using SfePy:
• Pyparsing,
• SciPy,
• meshio for reading and writing mesh files,
• scikit-umfpack for enabling UMFPACK solver for SciPy >= 0.14.0,
• Matplotlib for various plots,
• PyTables for storing results in HDF5 files,
• SymPy for some tests and functions,
• igakit for generating IGA domains,
• petsc4py and mpi4py for running parallel examples and using parallel solvers from PETSc,
• slepc4py for eigenvalue problem solvers from SLEPc,
• pymetis for mesh partitioning using Metis,
• Read the Docs Sphinx theme for building documentation,
• psutil for memory requirements checking,
• PyVista for post-processing.
Make sure the dependencies of those packages are also installed (e.g igakit reguires FORTRAN compiler, scikit-
umfpack does not work without UMFPACK, petsc4py without PETSc etc.). All dependencies of meshio need to be
installed for full mesh file format support (when using pip: pip install meshio[all]).
SfePy should work both with bleeding edge (Git) and last released versions of NumPy and SciPy. Please, submit an
issue at Issues page in case this does not hold.
Other dependencies/suggestions:
• To be able to (re)generate the documentation Sphinx, numpydoc and LaTeX are needed (see How to Regenerate
Documentation).

1.2. Installation 5
SfePy Documentation, Release version: 2024.2

• If doxygen is installed, the documentation of data structures and functions can be automatically generated by
running:

python setup.py doxygendocs

• Mesh generation tools use pexpect and gmsh or tetgen.

• IPython is recommended over the regular Python shell to fluently follow some parts of primer/tutorial (see Using
IPython).
• MUMPS library for using MUMPS linear direct solver (real and complex arithmetic, parallel factorization)

Compilation of C Extension Modules

In the SfePy top-level directory:

1. Look at site_cfg_template.py and follow the instructions therein. Usually no changes are necessary.
2. For in-place use, compile the extension modules:

python setup.py build_ext --inplace

After a successful compilation, SfePy can be used in-place. However, the the sfepy-* commands, such as
sfepy-run are only available after installing the package. Their functionality can be accessed by invoking
directly the corresponding scripts in sfepy/scripts/.

Installation

SfePy can be used without any installation by running its main scripts and examples from the top-level directory of the
distribution or can be installed locally or system-wide:
• system-wide (may require root privileges):

pip install .

• local:

pip install --user .

• development (editable install):

pip install -e .

The editable install allows working in-place and at the same time the sfepy-* commands are available.
If all went well, proceed with Testing Installation.

6 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

1.2.5 Testing Installation

After building and/or installing SfePy you should check if all the functions are working properly by running the auto-
mated tests.

Running Automated Test Suite

The test suite is based on pytest. Install it by:

pip install pytest

or:

conda install pytest

when working in Anaconda. If SfePy was installed, it can be tested using the command:

sfepy-test

that accepts all of the pytest options, for example:

sfepy-test -vv --durations=0 -m 'not slow' -k test_assembling.py

The tests output directory can also be specified:

sfepy-test --output-dir=output-tests

In general. the tests can be run using:

python -c "import sfepy; sfepy.test()"

in the SfePy top-level directory in case of the in-place build and anywhere else when testing the installed package.
Additional pytest options can be passed as arguments to sfepy.test(), for example:

python -c "import sfepy; sfepy.test('-vv', '--durations=0', '-m not slow', '-k test_
˓→assembling.py')"

Analogously to sfepy-test, the tests output directory can be specified using:

python -c "import sfepy; sfepy.test('--output-dir=output-tests')"

See pytest usage instructions for other options and usage patterns.
To test an in-place build (e.g. in a cloned git repository), the following simpler command can be used in the sources
top-level directory:

python -m pytest sfepy/tests

python -m pytest -v sfepy/tests/test_assembling.py

which will also add the current directory to sys.path. If the top-level directory is already in sys.path (e.g. using
export PYTHONPATH=.), the simplest way of invoking pytest is:

pytest sfepy/tests
pytest -v sfepy/tests/test_assembling.py

1.2. Installation 7
SfePy Documentation, Release version: 2024.2

1.2.6 Debugging

If something goes wrong, edit the site_cfg.py config file and set debug_flags = '-DDEBUG_FMF' to turn on
bound checks in the low level C functions, and recompile the code:

python setup.py clean

python setup.py build_ext --inplace

Then re-run your code and report the output to the SfePy mailing list.

1.2.7 Using IPython

We generally recommend to use (a customized) IPython interactive shell over the regular Python interpreter when
following Tutorial or Primer (or even for any regular interactive work with SfePy).
Install IPython (as a generic part of your selected distribution) and then customize it to your choice.
Depending on your IPython usage, you can customize your default profile or create a SfePy specific new one as follows:
1. Create a new SfePy profile:

ipython profile create sfepy

2. Open the ~/.ipython/profile_sfepy/ipython_config.py file in a text editor and add/edit after the c =
get_config() line:

exec_lines = [
'import numpy as nm',
'import matplotlib as mpl',
'mpl.use("WXAgg")',
#
# Add your preferred SfePy customization here...
#
]

c.InteractiveShellApp.exec_lines = exec_lines
c.TerminalIPythonApp.gui = 'wx'
c.TerminalInteractiveShell.colors = 'Linux' # NoColor, Linux, or LightBG

Please note, that generally it is not recommended to use star (*) imports here.
3. Run the customized IPython shell:

ipython --profile=sfepy

1.2.8 Notes on Multi-platform Python Distributions

Anaconda

We highly recommend this scientific-oriented Python distribution.

(Currently regularly tested by developers on SfePy releases with Python 3.6 64-bit on Ubuntu 16.04 LTS, Windows
8.1+ and macOS 10.12+.)

8 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Download appropriate Anaconda Python 3.x installer package and follow install instructions. We recommend to choose
user-level install option (no admin privileges required).
Anaconda can be used for:
1. installing the latest release of SfePy directly from the conda-forge channel (see sfepy-feedstock). In this case,
follow the instructions in Installing SfePy.
Installing/upgrading SfePy from the conda-forge channel can also be achieved by adding conda-forge to your
channels with:

conda config --add channels conda-forge

Once the conda-forge channel has been enabled, SfePy can be installed with:

conda install sfepy

It is possible to list all of the versions of SfePy available on your platform with:

conda search sfepy --channel conda-forge

2. installing the SfePy dependencies only - then proceed with the Installing SfePy from Sources instructions.
In this case, install the missing/required packages using built-in conda package manager:

conda install wxpython

See conda help for further information.

Occasionally, you should check for distribution and/or installed packages updates (there is no built-in automatic update
mechanism available):

conda update conda

conda update anaconda
conda update <package>

or try:

conda update --all

Compilation of C Extension Modules on Windows

To build SfePy extension modules, included mingw-w32/64 compiler tools should work fine. If you encounter any
problems, we recommend to install and use Microsoft Visual C++ Build Tools instead (see Anaconda FAQ).

1.3 Tutorial

1.3.1 Basic SfePy Usage

SfePy package can be used in two basic ways as a:

1. Black-box Partial Differential Equation (PDE) solver,
2. Python package to build custom applications involving solving PDEs by the Finite Element Method (FEM).

1.3. Tutorial 9
SfePy Documentation, Release version: 2024.2

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.

Invoking SfePy from the Command Line

This section introduces the basics of running SfePy from the command line, assuming it was properly installed, see
Installation.
The command sfepy-run is the most basic starting point in SfePy. It can be used to run declarative problem description
files (see below) as follows:

sfepy-run <problem_description_file>

Using SfePy Interactively

All functions of SfePy package can be also used interactively (see Interactive Example: Linear Elasticity for instance).
We recommend to use the IPython interactive shell for the best fluent user experience. You can customize your IPython
startup profile as described in Using IPython.

1.3.2 Basic Notions

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 and solution domain and subdomains (see sec-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. . .

Sneak Peek: What is Going on Under the Hood

1. A command (usually sfepy-run as in this tutorial) or a script reads in an input file.

2. 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 instance 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.

10 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

3. The solution is then obtained by calling problem.solve() function, which in turn calls a top-level time-stepping
solver. In each step, problem.time_update() is called to setup boundary conditions, material parameters and
other potentially time-dependent data. The problem.save_state() is called at the end of each time step to save the
results. This holds also for stationary problems with a single “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 instance is shown, see Interactive
Example: Linear Elasticity.
Now let us continue with running a simulation.

1.3.3 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.
• Download sfepy/examples/diffusion/poisson_short_syntax.py. It represents our sample SfePy sec-
problem-description-file, which defines the problem to be solved in terms SfePy can understand.
• Use the downloaded file in place of <problem_description_file.py> and run sfepy-run as described above. The
successful execution of the command creates output file cylinder.vtk in the SfePy top-level directory.

Postprocessing the Results

• The sfepy-view command can be used for quick postprocessing and visualization of the SfePy output files. It
requires pyvista installed on your system.
• As a simple example, try:
sfepy-view cylinder.vtk
• The following interactive 3D window should display:

• You can manipulate displayed image using:

1.3. Tutorial 11
SfePy Documentation, Release version: 2024.2

– 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.

1.3.4 Example Problem Description File

Here we discuss the contents of the sfepy/examples/diffusion/poisson_short_syntax.py problem descrip-

tion file. For additional examples, see the problem description files in the sfepy/examples/ directory of SfePy.
The problem at hand is the following:

𝑐∆𝑇 = 𝑓 in Ω, 𝑇 (𝑡) = 𝑇¯(𝑡) on Γ , (1.1)

where Γ ⊆ Ω is a subset of the domain Ω boundary. For simplicity, we set 𝑓 ≡ 0, 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.
𝜕𝑛 = 0 there. The particular boundary conditions used below are 𝑇 = 2 on the left side of the cylindrical domain
𝜕𝑇

depicted in the previous section and 𝑇 = −2 on the right side.

The first step to do is to write (1.1) in weak formulation (1.15). The 𝑓 = 0, 𝑔 = 𝜕𝑇
𝜕𝑛 = 0. So only one term in weak
form (1.15) remains:
∫︁
𝑐 ∇𝑇 · ∇𝑠 = 0, ∀𝑠 ∈ 𝑉0 . (1.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 sfepy/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'),
}

12 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

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 (1.2), namely: Find 𝑡 ∈ 𝑉ℎ , such that
∫︁
𝑠𝑇 ( 𝑐 𝐺𝑇 𝐺)𝑡 = 0, ∀𝑠 ∈ 𝑉ℎ0 ,
Ωℎ

where ∇𝑢 ≈ 𝐺𝑢.

1.3. Tutorial 13
SfePy Documentation, Release version: 2024.2

The equations block is the heart of the SfePy problem description 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 kinds and options. See sfepy.solvers.ls, sfepy.solvers.nls
and sfepy.solvers.ts_solvers for available solvers and their parameters.. Even linear problems are solved by a
nonlinear solver (KISS rule) – only one iteration is needed and the final residual is obtained for free. Note that we do
not need to define a time-stepping solver here - the problem is stationary and the default 'ts.stationary' solver is
created automatically.

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

The solvers to use are specified in the options block. We can define multiple solvers with different convergence param-
eters.
That’s it! Now it is possible to proceed as described in Invoking SfePy from the Command Line.

1.3.5 Interactive Example: Linear Elasticity

This example shows how to use SfePy interactively, but also how to make a custom simulation script. We will use
IPython interactive shell which allows more flexible and intuitive work (but you can use standard Python shell as well).
We wish to solve the following linear elasticity problem:

𝜕𝜎𝑖𝑗 (𝑢)
− + 𝑓𝑖 = 0 in Ω, 𝑢 = 0 on Γ1 , ¯1 on Γ2 ,
𝑢1 = 𝑢 (1.3)
𝜕𝑥𝑗
𝜕𝑢
where the stress is defined as 𝜎𝑖𝑗 = 2𝜇𝑒𝑖𝑗 +𝜆𝑒𝑘𝑘 𝛿𝑖𝑗 , 𝜆, 𝜇 are the Lamé’s constants, the strain is 𝑒𝑖𝑗 (𝑢) = 21 ( 𝜕𝑥
𝜕𝑢𝑖
𝑗
+ 𝜕𝑥𝑗𝑖 )
and 𝑓 are volume forces. This can be written in general form as 𝜎𝑖𝑗 (𝑢) = 𝐷𝑖𝑗𝑘𝑙 𝑒𝑘𝑙 (𝑢), where in our case 𝐷𝑖𝑗𝑘𝑙 =
𝜇(𝛿𝑖𝑘 𝛿𝑗𝑙 + 𝛿𝑖𝑙 𝛿𝑗𝑘 ) + 𝜆 𝛿𝑖𝑗 𝛿𝑘𝑙 .
In the weak form the equation (1.3) is
∫︁ ∫︁
𝐷𝑖𝑗𝑘𝑙 𝑒𝑘𝑙 (𝑢)𝑒𝑖𝑗 (𝑣) + 𝑓 𝑖 𝑣𝑖 = 0 , (1.4)
Ω Ω

14 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

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 available below in Complete Example as a Script.
In the SfePy top-level directory run

ipython

In [1]: import numpy as nm

In [2]: from sfepy.discrete.fem import Mesh, FEDomain, Field

Read a finite element mesh, that defines the domain Ω.

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

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

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

Define the regions – the whole domain Ω, where the solution is sought, and Γ1 , Γ2 , 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 [5]: min_x, max_x = domain.get_mesh_bounding_box()[:, 0]

In [6]: eps = 1e-8 * (max_x - min_x)
In [7]: omega = domain.create_region('Omega', 'all')
In [8]: gamma1 = domain.create_region('Gamma1',
...: 'vertices in x < %.10f ' % (min_x + eps),
...: 'facet')
In [9]: 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 [10]: field = Field.from_args('fu', nm.float64, 'vector', omega,

....: approx_order=2)

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

In [11]: from sfepy.discrete import (FieldVariable, Material, Integral, Function,

....: Equation, Equations, Problem)

In [12]: u = FieldVariable('u', 'unknown', field)

In [13]: 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 𝜆 = 1, 𝜇 = 1.
The volume forces will be defined also as a material as a constant (column) vector [0.02, 0.01]𝑇 .

In [14]: from sfepy.mechanics.matcoefs import stiffness_from_lame

(continues on next page)

1.3. Tutorial 15
SfePy Documentation, Release version: 2024.2

(continued from previous page)

In [15]: m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))
In [16]: 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 [17]: integral = Integral('i', order=3)

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

In [18]: from sfepy.terms import Term

In [19]: t1 = Term.new('dw_lin_elastic(m.D, v, u)',

....: integral, omega, m=m, v=v, u=u)

In [20]: t2 = Term.new('dw_volume_lvf(f.val, v)',

....: integral, omega, f=f, v=v)
In [21]: eq = Equation('balance', t1 + t2)
In [22]: eqs = Equations([eq])

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

In [23]: from sfepy.discrete.conditions import Conditions, EssentialBC

In [24]: fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})

In [25]: def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0):
....: val = shift * coors[:,1]**2
....: return val
In [26]: bc_fun = Function('shift_u_fun', shift_u_fun,
....: extra_args={'shift' : 0.01})
In [27]: 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 [28]: from sfepy.base.base import IndexedStruct

In [29]: from sfepy.solvers.ls import ScipyDirect
In [30]: from sfepy.solvers.nls import Newton

In [31]: ls = ScipyDirect({})
In [32]: nls_status = IndexedStruct()
In [33]: nls = Newton({}, lin_solver=ls, status=nls_status)

Now we are ready to create a Problem instance.

In [34]: pb = Problem('elasticity', equations=eqs)

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

In [35]: pb.save_regions_as_groups('regions')

And view them

16 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

In [36]: ! sfepy-view regions.vtk -2 --grid-vector1 "[2, 0, 0]"

You should see this:

Finally, we set the boundary conditions and the top-level solver , solve the problem, save and view the results. For
stationary problems, the top-level solver needs not to be a time-stepping solver - when a nonlinear solver is set instead,
the default 'ts.stationary' time-stepping solver is created automatically.

In [39]: pb.set_bcs(ebcs=Conditions([fix_u, shift_u]))

In [40]: pb.set_solver(nls)

In [41]: status = IndexedStruct()

In [42]: variables = pb.solve(status=status)

In [43]: print('Nonlinear solver status:\n', nls_status)

In [44]: print('Stationary solver status:\n', status)

In [45]: pb.save_state('linear_elasticity.vtk', variables)

This

sfepy-view linear_elasticity.vtk -2

is used to produce the resulting image:

1.3. Tutorial 17
SfePy Documentation, Release version: 2024.2

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

sfepy-view linear_elasticity.vtk -2 -f u:wu:p0 1:vw:p0

And the result is:

18 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Complete Example as a Script

The source code: linear_elastic_interactive.py.

This file should be run from the top-level SfePy source directory so it can find the mesh file correctly. Please note that
the provided example script may differ from above tutorial in some minor details.

1 #!/usr/bin/env python
2 """
3 Linear elasticity example using the imperative API.
4 """
5 from argparse import ArgumentParser
6 import numpy as nm
7

8 import sys
9 sys.path.append('.')
10

11 from sfepy.base.base import IndexedStruct

12 from sfepy.discrete import (FieldVariable, Material, Integral, Function,
13 Equation, Equations, Problem)
14 from sfepy.discrete.fem import Mesh, FEDomain, Field
15 from sfepy.terms import Term
16 from sfepy.discrete.conditions import Conditions, EssentialBC
17 from sfepy.solvers.ls import ScipyDirect
18 from sfepy.solvers.nls import Newton
19 from sfepy.mechanics.matcoefs import stiffness_from_lame
20

21

22 def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0):

23 """
24 Define a displacement depending on the y coordinate.
25 """
26 val = shift * coors[:,1]**2
27

28 return val
29

30

31 def main():
32 from sfepy import data_dir
33

34 parser = ArgumentParser()
35 parser.add_argument('--version', action='version', version='%(prog)s')
36 options = parser.parse_args()
37

38 mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh')

39 domain = FEDomain('domain', mesh)
40

41 min_x, max_x = domain.get_mesh_bounding_box()[:,0]

42 eps = 1e-8 * (max_x - min_x)
43 omega = domain.create_region('Omega', 'all')
44 gamma1 = domain.create_region('Gamma1',
45 'vertices in x < %.10f ' % (min_x + eps),
46 'facet')
47 gamma2 = domain.create_region('Gamma2',
(continues on next page)

1.3. Tutorial 19
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

48 'vertices in x > %.10f ' % (max_x - eps),
49 'facet')
50

51 field = Field.from_args('fu', nm.float64, 'vector', omega,

52 approx_order=2)
53

54 u = FieldVariable('u', 'unknown', field)

55 v = FieldVariable('v', 'test', field, primary_var_name='u')
56

57 m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))

58 f = Material('f', val=[[0.02], [0.01]])
59

60 integral = Integral('i', order=3)

61

62 t1 = Term.new('dw_lin_elastic(m.D, v, u)',
63 integral, omega, m=m, v=v, u=u)
64 t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
65 eq = Equation('balance', t1 + t2)
66 eqs = Equations([eq])
67

68 fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})

69

70 bc_fun = Function('shift_u_fun', shift_u_fun,

71 extra_args={'shift' : 0.01})
72 shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})
73

74 ls = ScipyDirect({})
75

76 nls_status = IndexedStruct()
77 nls = Newton({}, lin_solver=ls, status=nls_status)
78

79 pb = Problem('elasticity', equations=eqs)
80 pb.save_regions_as_groups('regions')
81

82 pb.set_bcs(ebcs=Conditions([fix_u, shift_u]))
83

84 pb.set_solver(nls)
85

86 status = IndexedStruct()
87 variables = pb.solve(status=status)
88

89 print('Nonlinear solver status:\n', nls_status)

90 print('Stationary solver status:\n', status)
91

92 pb.save_state('linear_elasticity.vtk', variables)
93

94

95 if __name__ == '__main__':
96 main()

20 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

1.4 User’s Guide

This manual provides reference documentation to SfePy from a user’s perspective.

1.4.1 Running a Simulation

The following should be run in the top-level directory of the SfePy source tree after installing the sfepy package. See
Installation for installation instructions.

Basic Usage

• sfepy-run sfepy/examples/diffusion/poisson_short_syntax.py

– Creates cylinder.vtk

• sfepy-run sfepy/examples/navier_stokes/stokes.py

– Creates channels_symm944t.vtk

Running Tests

The tests are based on pytest and can be run using:

python -c "import sfepy; sfepy.test()"

See Testing Installation for additional information.

Computations and Examples

Most of the example problems in the sfepy/examples directory can be computed by the sfepy-run command.
Additional (stand-alone) examples can be run directly, e.g.:

python sfepy/examples/large_deformation/compare_elastic_materials.py

Common Tasks

• Run a simulation:

sfepy-run sfepy/examples/diffusion/poisson_short_syntax.py
sfepy-run sfepy/examples/diffusion/poisson_short_syntax.py -o some_results # ->␣
˓→produces some_results.vtk

• Print available terms:

sfepy-run --list=terms

• Run a simulation and also save Dirichlet boundary conditions:

sfepy-run --save-ebc sfepy/examples/diffusion/poisson_short_syntax.py # -> produces␣

˓→an additional .vtk file with BC visualization

1.4. User’s Guide 21

SfePy Documentation, Release version: 2024.2

• Use a restart file to continue an interrupted simulation:

– Warning: This feature is preliminary and does not support terms with internal state.
– Run:

sfepy-run sfepy/examples/large_deformation/balloon.py --save-restart=-1

and break the computation after a while (hit Ctrl-C). The mode --save-restart=-1 is currently the only
supported mode. It saves a restart file for each time step, and only the last computed time step restart file is
kept.
– A file named 'unit_ball.restart-??.h5' should be created, where '??' indicates the last stored time
step. Let us assume it is 'unit_ball.restart-04.h5', i.e. the fifth step.
– Restart the simulation by:

sfepy-run sfepy/examples/large_deformation/balloon.py --load-restart=unit_ball.

˓→restart-04.h5

The simulation should continue from the next time step. Verify that by running:

sfepy-run sfepy/examples/large_deformation/balloon.py

and compare the residuals printed in the corresponding time steps.

1.4.2 Visualization of Results

Quick visualisation of the SfePy results can be done by sfepy.scripts.resview.py command, which uses PyVista
visualization toolkit (need to be installed).
For example, results of a stationary Navier-Stokes flow simulation:

sfepy-run sfepy/examples/navier_stokes/navier_stokes.py -o navier_stokes

can be viewed using:

sfepy-view navier_stokes.vtk

produces:

22 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Using:

sfepy-view navier_stokes.vtk -f p:i5:p0 p:e:o.2:p0 u:t1000:p1 u:o.2:p1

the output is split into plots plot:0 and plot:1, where these plots contain:
• plot:0: field p, 5 isosurfaces, mesh edges switched on, the surface opacity set to 0.2;
• plot:1: vector field u streamlines, the surface opacity set to 0.2;

1.4. User’s Guide 23

SfePy Documentation, Release version: 2024.2

The actual camera position can be printed by pressing the ‘c’ key. The above figures use the following values:

--camera-position="-0.236731,-0.352216,0.237044,0.094458,0.0292563,0.0385116,0.266969,0.
˓→252278,0.930099"

The argument -o filename.png takes the screenshot of the produced view:

sfepy-view navier_stokes.vtk -o image.png

1.4.3 Problem Description File

Here we discuss the basic items that users have to specify in their input files. For complete examples, see the problem
description files in the sfepy/examples/ directory of SfePy.

24 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Long Syntax

Besides the short syntax described below there is (due to history) also a long syntax which is explained in prob-
lem_desc_file_long. The short and long syntax can be mixed together in one description file.

FE Mesh

A FE mesh defining a domain geometry can be stored in several formats:

• legacy VTK (.vtk)
• custom HDF5 file (.h5)
• medit mesh file (.mesh)
• tetgen mesh files (.node, .ele)
• comsol text mesh file (.txt)
• abaqus text mesh file (.inp)
• avs-ucd text mesh file (.inp)
• hypermesh text mesh file (.hmascii)
• hermes3d mesh file (.mesh3d)
• nastran text mesh file (.bdf)
• gambit neutral text mesh file (.neu)
• salome/pythonocc med binary mesh file (.med)
Example:

filename_mesh = 'meshes/3d/cylinder.vtk'

The VTK and HDF5 formats can be used for storing the results. The format can be selected in options, see Miscella-
neous.
The following geometry elements are supported:

1.4. User’s Guide 25

SfePy Documentation, Release version: 2024.2

Note the orientation of the vertices matters, the figure displays the correct orientation when interpreted in a right-handed
coordinate system.

Regions

Regions serve to select a certain part of the computational domain using topological entities of the FE mesh. They are
used to define the boundary conditions, the domains of terms and materials etc.
Let us denote D the maximal dimension of topological entities. For volume meshes it is also the dimension of space
the domain is embedded in. Then the following topological entities can be defined on the mesh (notation follows
[Logg2012]):

topological entity dimension co-dimension

vertex 0 D
edge 1 D-1
face 2 D-2
facet D-1 1
cell D 0

If D = 2, faces are not defined and facets are edges. If D = 3, facets are faces.
Following the above definitions, a region can be of different kind:
• cell, facet, face, edge, vertex - entities of higher dimension are not included.
• cell_only, facet_only, face_only, edge_only, vertex_only - only the specified entities are included,
other entities are empty sets, so that set-like operators still work, see below.

26 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

• The cell kind is the most general and should be used with cell terms. It is also the default if the kind is not
specified in region definition.
• The facet kind (same as edge in 2D and face in 3D) is to be used with boundary (facet integral) terms.
• The vertex (same as vertex_only) kind can be used with point-wise defined terms (e.g. point loads).
The kinds allow a clear distinction between regions of different purpose (cell integration domains, facet domains, etc.)
and could be used to lower memory usage.
A region definition involves topological entity selections combined with set-like operators. The set-like operators can
result in intermediate regions that have the cell kind. The desired kind is set to the final region, removing unneeded
entities. Most entity selectors are defined in terms of vertices and cells - the other entities are computed as needed.

Topological Entity Selection

topological entity selection explanation

all all entities of the mesh
vertices of surface surface of the mesh
vertices of group <integer> vertices of given group
vertices of set <str> vertices of a given named vertex set2
vertices in <expr> vertices given by an expression3
vertices by <function> vertices given by a function of coordinates4
vertex <id>[, <id>, ...] vertices given by their ids
vertex in r.<name of another region> any single vertex in the given region
cells of group <integer> cells of given group
cells by <efunction> cells given by a function of coordinates5
cell <id>[, <id>, ...], cells given by their ids
copy r.<name of another region> a copy of the given region
r.<name of another region> a reference to the given region

2 Only if mesh format supports reading boundary condition vertices as vertex sets.
3 <expr> is a logical expression like (y <= 0.1) & (x < 0.2). In 2D use x, y, in 3D use x, y and z. & stands for logical and, | stands for
logical or.
4 <function> is a function with signature fun(coors, domain=None), where coors are coordinates of mesh vertices.
5 <efunction> is a function with signature fun(coors, domain=None), where coors are coordinates of mesh cell centroids.

1.4. User’s Guide 27

SfePy Documentation, Release version: 2024.2

topological entity selection footnotes

set-like operator explanation

+v vertex union
+e edge union
+f face union
+s facet union
+c cell union
-v vertex difference
-e edge difference
-f face difference
-s facet difference
-c cell difference
*v vertex intersection
*e edge intersection
*f face intersection
*s facet intersection
*c cell intersection

Region Definition Syntax

Regions are defined by the following Python dictionary:

regions = {
<name> : (<selection>, [<kind>], [<parent>], [{<misc. options>}]),
}

or:

regions = {
<name> : <selection>,
}

Example definitions:

regions = {
'Omega' : 'all',
'Right' : ('vertices in (x > 0.99)', 'facet'),
'Gamma1' : ("""(cells of group 1 *v cells of group 2)
+v r.Right""", 'facet', 'Omega'),
'Omega_B' : 'vertices by get_ball',
}

The Omega_B region illustrates the selection by a function (see Topological Entity Selection). In this example, the
function is:

import numpy as nm

def get_ball(coors, domain=None):

x, y, z = coors[:, 0], coors[:, 1], coors[:, 2]

(continues on next page)

28 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

r = nm.sqrt(x**2 + y**2 + z**2)
flag = nm.where((r < 0.1))[0]

return flag

The function needs to be registered in Functions:

functions = {
'get_ball' : (get_ball,),
}

The mirror region can be defined explicitly as:

regions = {
'Top': ('r.Y *v r.Surf1', 'facet', 'Y', {'mirror_region': 'Bottom'}),
'Bottom': ('r.Y *v r.Surf2', 'facet', 'Y', {'mirror_region': 'Top'}),
}

Fields

Fields correspond to FE spaces:

fields = {
<name> : (<data_type>, <shape>, <region_name>, <approx_order>,
[<space>, <poly_space_base>])
}

where
• <data_type> is a numpy type (float64 or complex128) or ‘real’ or ‘complex’
• <shape> is the number of DOFs per node: 1 or (1,) or ‘scalar’, space dimension (2, or (2,) or 3 or (3,)) or
‘vector’; it can be other positive integer than just 1, 2, or 3
• <region_name> is the name of region where the field is defined
• <approx_order> is the FE approximation order, e.g. 0, 1, 2, ‘1B’ (1 with bubble)
• <space> is the function space
• <poly_space_base> is the basis of the FE (usually polynomial) space
Example: scalar P1 elements in 2D on a region Omega:

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

The following approximation orders can be used:

• simplex elements: 1, 2, ‘1B’, ‘2B’
• tensor product elements: 0, 1, ‘1B’
Optional bubble function enrichment is marked by ‘B’.

1.4. User’s Guide 29

SfePy Documentation, Release version: 2024.2

The implemented combinations of spaces and bases are listed below, the space column corresponds to <space>, the
basis column to <poly_space_base> and region type to the field region type.

Table 1: Fields
space basis region kind description
H1 bernstein cell, facet Bernstein basis approximation with positive-only basis function
values.
H1 iga cell Bezier extraction based NURBS approximation for isogeometric
analysis.
H1 lagrange cell, facet Lagrange basis nodal approximation.
H1 lagrange_discontinuous cell The C0 constant-per-cell approximation.
H1 lobatto cell Hierarchical basis approximation with Lobatto polynomials.
H1 sem cell, facet Spectral element method approximation.
H1 serendipity cell, facet Lagrange basis nodal serendipity approximation with order <= 3.
H1 shell10x cell The approximation for the shell10x element.
L2 constant cell, facet The L2 constant-in-a-region approximation.
DG legendre_discontinuous cell Discontinuous Galerkin method approximation with Legendre
basis.

Variables

Variables use the FE approximation given by the specified field:

variables = {
<name> : (<kind>, <field_name>, <spec>, [<history>])
}

where
• <kind> - ‘unknown field’, ‘test field’ or ‘parameter field’
• <spec> - in case of: primary variable - order in the global vector of unknowns, dual variable - name of
primary variable
• <history> - number of time steps to remember prior to current step
Example:

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

Integrals

Define the integral type and quadrature rule. This keyword is optional, as the integration orders can be specified directly
in equations (see below):

integrals = {
<name> : <order>
}

where

30 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

• <name> - the integral name - it has to begin with ‘i’!

• <order> - the order of polynomials to integrate, or ‘custom’ for integrals with explicitly given values and
weights
Example:

import numpy as nm
N = 2
integrals = {
'i1' : 2,
'i2' : ('custom', zip(nm.linspace( 1e-10, 0.5, N ),
nm.linspace( 1e-10, 0.5, N )),
[1./N] * N),
}

Essential Boundary Conditions and Constraints

The essential boundary conditions set values of DOFs in some regions, while the constraints constrain or transform
values of DOFs in some regions.

Dirichlet Boundary Conditions

The Dirichlet, or essential, boundary conditions apply in a given region given by its name, and, optionally, in selected
times. The times can be given either using a list of tuples (t0, t1) making the condition active for t0 <= t < t1, or by a
name of a function taking the time argument and returning True or False depending on whether the condition is active
at the given time or not.
Dirichlet (essential) boundary conditions:

ebcs = {
<name> : (<region_name>, [<times_specification>,]
{<dof_specification> : <value>[,
<dof_specification> : <value>, ...]})
}

Example:

ebcs = {
'u1' : ('Left', {'u.all' : 0.0}),
'u2' : ('Right', [(0.0, 1.0)], {'u.0' : 0.1}),
'phi' : ('Surface', {'phi.all' : 0.0}),
'u_yz' : ('Gamma', {'u.[1,2]' : 'rotate_yz'}),
}

The u_yz condition illustrates calculating the condition value by a function. In this example, it is a function of coordi-
nates coors of region nodes:

import numpy as nm

def rotate_yz(ts, coor, **kwargs):

from sfepy.linalg import rotation_matrix2d

(continues on next page)

1.4. User’s Guide 31

SfePy Documentation, Release version: 2024.2

(continued from previous page)

vec = coor[:,1:3] - centre

angle = 10.0 * ts.step

mtx = rotation_matrix2d(angle)
vec_rotated = nm.dot(vec, mtx)

displacement = vec_rotated - vec

return displacement

The function needs to be registered in Functions:

functions = {
'rotate_yz' : (rotate_yz,),
}

Periodic Boundary Conditions

The periodic boundary conditions tie DOFs of a single variable in two regions that have matching nodes. Can be used
with functions in sfepy.discrete.fem.periodic.
Periodic boundary conditions:

epbcs = {
<name> : ((<region1_name>, <region2_name>), [<times_specification>,]
{<dof_specification> : <dof_specification>[,
<dof_specification> : <dof_specification>, ...]},
<match_function_name>)
}

Example:

epbcs = {
'up1' : (('Left', 'Right'), {'u.all' : 'u.all', 'p.0' : 'p.0'},
'match_y_line'),
}

Linear Combination Boundary Conditions

The linear combination boundary conditions (LCBCs) are more general than the Dirichlet BCs or periodic BCs. They
can be used to substitute one set of DOFs in a region by another set of DOFs, possibly in another region and of another
variable. The LCBCs can be used only in FEM with nodal (Lagrange) basis.
Available LCBC kinds:
• 'rigid' - in linear elasticity problems, a region moves as a rigid body;
• 'no_penetration' - in flow problems, the velocity vector is constrained to the plane tangent to the surface;
• 'normal_direction' - the velocity vector is constrained to the normal direction;
• 'edge_direction' - the velocity vector is constrained to the mesh edge direction;

32 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

• 'integral_mean_value' - all DOFs in a region are summed to a single new DOF;

• 'shifted_periodic' - generalized periodic BCs that work with two different variables and can have a non-zero
mutual shift;
• 'match_dofs' - tie DOFs of two fields.
Only the 'shifted_periodic' LCBC needs the second region and the DOF mapping function, see below.
Linear combination boundary conditions:

lcbcs = {
'shifted' : (('Left', 'Right'),
{'u1.all' : 'u2.all'},
'match_y_line', 'shifted_periodic',
'get_shift'),
'mean' : ('Middle', {'u1.all' : None}, None, 'integral_mean_value'),
}

Initial Conditions

Initial conditions are applied prior to the boundary conditions - no special care must be used for the boundary dofs:

ics = {
<name> : (<region_name>, {<dof_specification> : <value>[,
<dof_specification> : <value>, ...]},...)
}

Example:

ics = {
'ic' : ('Omega', {'T.0' : 5.0}),
}

Materials

Materials are used to define constitutive parameters (e.g. stiffness, permeability, or viscosity), and other non-field
arguments of terms (e.g. known traction or volume forces). Depending on a particular term, the parameters can be
constants, functions defined over FE mesh nodes, functions defined in the elements, etc.
Example:

material = {
'm' : ({'val' : [0.0, -1.0, 0.0]},),
'm2' : 'get_pars',
'm3' : (None, 'get_pars'), # Same as the above line.
}

Example: different material parameters in regions ‘Yc’, ‘Ym’:

from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson

dim = 3
materials = {
'mat' : ({'D' : {
(continues on next page)

1.4. User’s Guide 33

SfePy Documentation, Release version: 2024.2

(continued from previous page)

'Ym': stiffness_from_youngpoisson(dim, 7.0e9, 0.4),
'Yc': stiffness_from_youngpoisson(dim, 70.0e9, 0.2)}
},),
}

Defining Material Parameters by Functions

The functions for defining material parameters can work in two modes, distinguished by the mode argument. The two
modes are ‘qp’ and ‘special’. The first mode is used for usual functions that define parameters in quadrature points
(hence ‘qp’), while the second one can be used for special values like various flags.
The shape and type of data returned in the ‘special’ mode can be arbitrary (depending on the term used). On the other
hand, in the ‘qp’ mode all the data have to be numpy float64 arrays with shape (n_coor, n_row, n_col), where n_coor
is the number of quadrature points given by the coors argument, n_coor = coors.shape[0], and (n_row, n_col) is the
shape of a material parameter in each quadrature point. For example, for scalar parameters, the shape is (n_coor, 1, 1).
The shape is determined by each term.
Example:

def get_pars(ts, coors, mode=None, **kwargs):

if mode == 'qp':
val = coors[:,0]
val.shape = (coors.shape[0], 1, 1)

return {'x_coor' : val}

The function needs to be registered in Functions:

functions = {
'get_pars' : (get_pars,),
}

If a material parameter has the same value in all quadrature points, than it is not necessary to repeat the constant and
the array can be with shape (1, n_row, n_col).

Equations and Terms

Equations can be built by combining terms listed in Term Table.

Examples

• Laplace equation, named integral:

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

• Laplace equation, simplified integral given by order:

34 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

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

• Laplace equation, automatic integration order (not implemented yet!):

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

• Navier-Stokes equations:

equations = {
'balance' :
"""+ dw_div_grad.i2.Omega( fluid.viscosity, v, u )
+ dw_convect.i2.Omega( v, u )
- dw_stokes.i1.Omega( v, p ) = 0""",
'incompressibility' :
"""dw_stokes.i1.Omega( u, q ) = 0""",
}

Configuring Solvers

In SfePy, a non-linear solver has to be specified even when solving a linear problem. The linear problem is/should be
then solved in one iteration of the nonlinear solver.
Linear and nonlinear solver:

solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton',
{'i_max' : 1}),
}

Solver selection:

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

For the case that a chosen linear solver is not available, it is possible to define the fallback option of the chosen solver
which specifies a possible alternative:

solvers = {
'ls': ('ls.mumps', {'fallback': 'ls2'}),
'ls2': ('ls.scipy_umfpack', {}),
'newton': ('nls.newton', {
'i_max' : 1,
'eps_a' : 1e-10,
}),
}

1.4. User’s Guide 35

SfePy Documentation, Release version: 2024.2

Another possibility is to use a “virtual” solver that ensures an automatic selection of an available solver, see Virtual
Linear Solvers with Automatic Selection.

Functions

Functions are a way of customizing SfePy behavior. They make it possible to define material properties, boundary
conditions, parametric sweeps, and other items in an arbitrary manner. Functions are normal Python functions declared
in the Problem Definition file, so they can invoke the full power of Python. In order for SfePy to make use of the
functions, they must be declared using the function keyword. See the examples below, and also the corresponding
sections above.

Examples

See sfepy/examples/diffusion/poisson_functions.py for a complete problem description file demonstrating

how to use different kinds of functions.
• functions for defining regions:
def get_circle(coors, domain=None):
r = nm.sqrt(coors[:,0]**2.0 + coors[:,1]**2.0)
return nm.where(r < 0.2)[0]

functions = {
'get_circle' : (get_circle,),
}

• functions for defining boundary conditions:

def get_p_edge(ts, coors, bc=None, problem=None):
if bc.name == 'p_left':
return nm.sin(nm.pi * coors[:,1])
else:
return nm.cos(nm.pi * coors[:,1])

functions = {
'get_p_edge' : (get_p_edge,),
}

ebcs = {
'p' : ('Gamma', {'p.0' : 'get_p_edge'}),
}

The values can be given by a function of time, coordinates and possibly other data, for example:
ebcs = {
'f1' : ('Gamma1', {'u.0' : 'get_ebc_x'}),
'f2' : ('Gamma2', {'u.all' : 'get_ebc_all'}),
}

def get_ebc_x(coors, amplitude):

z = coors[:, 2]
val = amplitude * nm.sin(z * 2.0 * nm.pi)
return val
(continues on next page)

36 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

def get_ebc_all(ts, coors):

val = ts.step * coors
return val

functions = {
'get_ebc_x' : (lambda ts, coors, bc, problem, **kwargs:
get_ebc_x(coors, 5.0),),
'get_ebc_all' : (lambda ts, coors, bc, problem, **kwargs:
get_ebc_all(ts, coors),),
}

Note that when setting more than one component as in get_ebc_all() above, the function should return either
an array of shape (coors.shape[0], n_components), or the same array flattened to 1D row-by-row (i.e. node-by-
node), where n_components corresponds to the number of components in the boundary condition definition. For
example, with ‘u.[0, 1]’, n_components is 2.
• function for defining usual material parameters:

def get_pars(ts, coors, mode=None, **kwargs):

if mode == 'qp':
val = coors[:,0]
val.shape = (coors.shape[0], 1, 1)

return {'x_coor' : val}

functions = {
'get_pars' : (get_pars,),
}

The keyword arguments contain both additional use-specified arguments, if any, and the following: equations,
term, problem, for cases when the function needs access to the equations, problem, or term instances that
requested the parameters that are being evaluated. The full signature of the function is:

def get_pars(ts, coors, mode=None,

equations=None, term=None, problem=None, **kwargs)

• function for defining special material parameters, with an extra argument:

def get_pars_special(ts, coors, mode=None, extra_arg=None):

if mode == 'special':
if extra_arg == 'hello!':
ic = 0
else:
ic = 1
return {('x_%s' % ic) : coors[:,ic]}

functions = {
'get_pars1' : (lambda ts, coors, mode=None, **kwargs:
get_pars_special(ts, coors, mode,
extra_arg='hello!'),),
}

(continues on next page)

1.4. User’s Guide 37

SfePy Documentation, Release version: 2024.2

(continued from previous page)

# Just another way of adding a function, besides 'functions' keyword.
function_1 = {
'name' : 'get_pars2',
'function' : lambda ts, coors, mode=None, **kwargs:
get_pars_special(ts, coors, mode, extra_arg='hi!'),
}

• function combining both kinds of material parameters:

def get_pars_both(ts, coors, mode=None, **kwargs):

out = {}

if mode == 'special':

out['flag'] = coors.max() > 1.0

elif mode == 'qp':

val = coors[:,1]
val.shape = (coors.shape[0], 1, 1)

out['y_coor'] = val

return out

functions = {
'get_pars_both' : (get_pars_both,),
}

• function for setting values of a parameter variable:

variables = {
'p' : ('parameter field', 'temperature',
{'setter' : 'get_load_variable'}),
}

def get_load_variable(ts, coors, region=None, variable=None, **kwargs):

y = coors[:,1]
val = 5e5 * y
return val

functions = {
'get_load_variable' : (get_load_variable,)
}

38 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Miscellaneous

The options can be used to select solvers, output file format, output directory, to register functions to be called at various
phases of the solution (the hooks), and for other settings.
Additional options (including solver selection):

options = {
# int >= 0, uniform mesh refinement level
'refinement_level : 0',

# bool, default: False, if True, allow selecting empty regions with no

# entities
'allow_empty_regions' : True,

# string, output directory

'output_dir' : 'output/<output_dir>',

# 'vtk' or 'h5', output file (results) format

'output_format' : 'vtk',

# output file format variant compatible with 'output_format'

'file_format' : 'vtk-ascii',

# string, nonlinear solver name

'nls' : 'newton',

# string, linear solver name

'ls' : 'ls',

# string, time stepping solver name

'ts' : 'ts',

# The times at which results should be saved:

# - a sequence of times
# - or 'all' for all time steps (the default value)
# - or an int, number of time steps, spaced regularly from t0 to t1
# - or a function `is_save(ts)`
'save_times' : 'all',

# save a restart file for each time step, only the last computed time
# step restart file is kept.
'save_restart' : -1,

# string, a function to be called after each time step

'step_hook' : '<step_hook_function>',

# string, a function to be called after each time step, used to

# update the results to be saved
'post_process_hook' : '<post_process_hook_function>',

# string, as above, at the end of simulation

'post_process_hook_final' : '<post_process_hook_final_function>',

(continues on next page)

1.4. User’s Guide 39

SfePy Documentation, Release version: 2024.2

(continued from previous page)

# string, a function to generate probe instances
'gen_probes' : '<gen_probes_function>',

# string, a function to probe data

'probe_hook' : '<probe_hook_function>',

# string, a function to modify problem definition parameters

'parametric_hook' : '<parametric_hook_function>',

# float, default: 1e-9. If the distance between two mesh vertices

# is less than this value, they are considered identical.
# This affects:
# - periodic regions matching
# - mirror regions matching
# - fixing of mesh doubled vertices
'mesh_eps': 1e-7,

# bool, default: True. If True, the (tangent) matrices and residual

# vectors (right-hand sides) contain only active DOFs, otherwise all
# DOFs (including the ones fixed by the Dirichlet or periodic boundary
# conditions) are included. Note that the rows/columns corresponding to
# fixed DOFs are modified w.r.t. a problem without the boundary
# conditions.
'active_only' : False,

# bool, default: False. If True, all DOF connectivities are used to

# pre-allocate the matrix graph. If False, only cell region
# connectivities are used.
'any_dof_conn' : False,

# bool, default: False. If True, automatically transform equations to a

# form suitable for the given solver. Implemented for
# ElastodynamicsBaseTS-based solvers
'auto_transform_equations' : True,
}

• post_process_hook enables computing derived quantities, like stress or strain, from the primary unknown
variables. See the examples in sfepy/examples/large_deformation/ directory.
• parametric_hook makes it possible to run parametric studies by modifying the problem description program-
matically. See sfepy/examples/diffusion/poisson_parametric_study.py for an example.
• output_dir redirects output files to specified directory

40 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

1.4.4 Building Equations in SfePy

Equations in SfePy are built using terms, which correspond directly to the integral forms of weak formulation of a
problem to be solved. As an example, let us consider the Laplace equation in time interval 𝑡 ∈ [0, 𝑡final ]:

𝜕𝑇
+ 𝑐∆𝑇 = 0 in Ω, 𝑇 (𝑡) = 𝑇¯(𝑡) on Γ . (1.5)
𝜕𝑡
The weak formulation of (1.5) is: Find 𝑇 ∈ 𝑉 , such that
∫︁ ∫︁
𝜕𝑇
𝑠 + 𝑐 ∇𝑇 : ∇𝑠 = 0, ∀𝑠 ∈ 𝑉0 , (1.6)
Ω 𝜕𝑡 Ω

where we assume no fluxes over 𝜕Ω ∖ Γ. In the syntax used in SfePy input files, this can be written as:

dw_dot.i.Omega( s, dT/dt ) + dw_laplace.i.Omega( coef, s, T) = 0

which directly corresponds to the discrete version of (1.6): Find 𝑇 ∈ 𝑉ℎ , such that
∫︁ ∫︁
𝑇 𝑇 𝜕𝑇 𝑇
𝑠 ( 𝜑 𝜑) +𝑠 ( 𝑐 𝐺𝑇 𝐺)𝑇 = 0, ∀𝑠 ∈ 𝑉ℎ0 ,
Ωℎ 𝜕𝑡 Ωℎ

where 𝑢 ≈ 𝜑𝑢, ∇𝑢 ≈ 𝐺𝑢 for 𝑢 ∈ {𝑠, 𝑇 }. The integrals over the discrete domain Ωℎ are approximated by a numerical
quadrature, that is named i in our case.

Syntax of Terms in Equations

The terms in equations are written in form:

<term_name>.<i>.<r>( <arg1>, <arg2>, ... )

where <i> denotes an integral name (i.e. a name of numerical quadrature to use) and <r> marks a region (domain
of the integral). In the following, <virtual> corresponds to a test function, <state> to a unknown function and
<parameter> to a known function arguments.
When solving, the individual terms in equations are evaluated in the ‘weak’ mode. The evaluation modes are described
in the next section.

1.4.5 Term Evaluation

Terms can be evaluated in two ways:

1. implicitly by using them in equations;
2. explicitly using Problem.evaluate(). This way is mostly used in the postprocessing.
Each term supports one or more evaluation modes:
• ‘weak’ : Assemble (in the finite element sense) either the vector or matrix depending on diff_var argument (the
name of variable to differentiate with respect to) of Term.evaluate(). This mode is usually used implicitly
when building the linear system corresponding to given equations.
• ‘eval’ : The evaluation mode integrates the term (= integral) over a region. The result has the same dimension
as the quantity being integrated. This mode can be used, for example, to compute some global quantities during
postprocessing such as fluxes or total values of extensive quantities (mass, volume, energy, . . . ).
• ‘el_eval’ : The element evaluation mode results in an array of a quantity integrated over each element of a region.

1.4. User’s Guide 41

SfePy Documentation, Release version: 2024.2

• ‘el_avg’ : The element average mode results in an array of a quantity averaged in each element of a region. This
is the mode for postprocessing.
• ‘qp’ : The quadrature points mode results in an array of a quantity interpolated into quadrature points of each
element in a region. This mode is used when further point-wise calculations with the result are needed. The
same element type and number of quadrature points in each element are assumed.
Not all terms support all the modes - consult the documentation of the individual terms. There are, however, certain
naming conventions:
• ‘dw_*’ terms support ‘weak’ mode
• ‘dq_*’ terms support ‘qp’ mode
• ‘d_*’, ‘di_*’ terms support ‘eval’ and ‘el_eval’ modes
• ‘ev_*’ terms support ‘eval’, ‘el_eval’, ‘el_avg’ and ‘qp’ modes
Note that the naming prefixes are due to history when the mode argument to Problem.evaluate() and Term.
evaluate() was not available. Now they are often redundant, but are kept around to indicate the evaluation purpose
of each term.
Several examples of using the Problem.evaluate() function are shown below.

1.4.6 Solution Postprocessing

A solution to equations given in a problem description file is given by the variables of the ‘unknown field’ kind, that
are set in the solution procedure. By default, those are the only values that are stored into a results file. The solution
postprocessing allows computing additional, derived, quantities, based on the primary variables values, as well as any
other quantities to be stored in the results.
Let us illustrate this using several typical examples. Let us assume that the postprocessing function is
called ‘post_process()’, and is added to options as discussed in Miscellaneous, see ‘post_process_hook’ and
‘post_process_hook_final’. Then:
• compute stress and strain given the displacements (variable u):

def post_process(out, problem, variables, extend=False):

"""
This will be called after the problem is solved.

Parameters
----------
out : dict
The output dictionary, where this function will store additional
data.
problem : Problem instance
The current Problem instance.
variables : Variables instance
The computed state, containing FE coefficients of all the unknown
variables.
extend : bool
The flag indicating whether to extend the output data to the whole
domain. It can be ignored if the problem is solved on the whole
domain already.

Returns
(continues on next page)

42 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

-------
out : dict
The updated output dictionary.
"""
from sfepy.base.base import Struct

# Cauchy strain averaged in elements.

strain = problem.evaluate('ev_cauchy_strain.i.Omega(u)',
mode='el_avg')
out['cauchy_strain'] = Struct(name='output_data',
mode='cell', data=strain,
dofs=None)
# Cauchy stress averaged in elements.
stress = problem.evaluate('ev_cauchy_stress.i.Omega(solid.D, u)',
mode='el_avg')
out['cauchy_stress'] = Struct(name='output_data',
mode='cell', data=stress,
dofs=None)

return out

The full example is linear_elasticity/linear_elastic_probes.py.

• compute diffusion velocity given the pressure:

def post_process(out, pb, state, extend=False):

from sfepy.base.base import Struct

dvel = pb.evaluate('ev_diffusion_velocity.i.Omega(m.K, p)',

mode='el_avg')
out['dvel'] = Struct(name='output_data',
mode='cell', data=dvel, dofs=None)

return out

The full example is biot-biot_npbc.

• store values of a non-homogeneous material parameter:

def post_process(out, pb, state, extend=False):

from sfepy.base.base import Struct

mu = pb.evaluate('ev_integrate_mat.2.Omega(nonlinear.mu, u)',
mode='el_avg', copy_materials=False, verbose=False)
out['mu'] = Struct(name='mu', mode='cell', data=mu, dofs=None)

return out

The full example is linear_elasticity/material_nonlinearity.py.

• compute volume of a region (u is any variable defined in the region Omega):

volume = problem.evaluate('ev_volume.2.Omega(u)')

1.4. User’s Guide 43

SfePy Documentation, Release version: 2024.2

1.4.7 Probing

Probing applies interpolation to output the solution along specified paths. There are two ways of probing:
• VTK probes: It is the simple way of probing using the ‘post_process_hook’. It generates matplotlib figures with
the probing results and previews of the mesh with the probe paths. See Primer or linear_elasticity/its2D_5.py
example.
• SfePy probes: As mentioned in Miscellaneous, it relies on defining two additional functions, namely the
‘gen_probes’ function, that should create the required probes (see sfepy.discrete.probes), and the
‘probe_hook’ function that performs the actual probing of the results for each of the probes. This function can re-
turn the probing results, as well as a handle to a corresponding matplotlib figure. See linear_elasticity/its2D_4.py
for additional explanation.
Using sfepy.discrete.probes allows correct probing of fields with the approximation order greater than one,
see Interactive Example in Primer or linear_elasticity/its2D_interactive.py for an example of interactive use.

1.4.8 Postprocessing filters

The following postprocessing functions based on the VTK filters are available:
• ‘get_vtk_surface’: extract mesh surface
• ‘get_vtk_edges’: extract mesh edges
• ‘get_vtk_by_group’: extract domain by a material ID
• ‘tetrahedralize_vtk_mesh’: 3D cells are converted to tetrahedral meshes, 2D cells to triangles
The following code demonstrates the use of the postprocessing filters:

mesh = problem.domain.mesh
mesh_name = mesh.name[mesh.name.rfind(osp.sep) + 1:]

vtkdata = get_vtk_from_mesh(mesh, out, 'postproc_')

matrix = get_vtk_by_group(vtkdata, 1, 1)

matrix_surf = get_vtk_surface(matrix)
matrix_surf_tri = tetrahedralize_vtk_mesh(matrix_surf)
write_vtk_to_file('%s_mat1_surface.vtk' % mesh_name, matrix_surf_tri)

matrix_edges = get_vtk_edges(matrix)
write_vtk_to_file('%s_mat1_edges.vtk' % mesh_name, matrix_edges)

1.4.9 Solvers

This section describes the time-stepping, nonlinear, linear, eigenvalue and optimization solvers available in SfePy.
There are many internal and external solvers in the sfepy.solvers package that can be called using a uniform interface.

44 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Time-stepping solvers

All PDEs that can be described in a problem description file are solved internally by a time-stepping solver. This holds
even for stationary problems, where the default single-step solver ('ts.stationary') is created automatically. In this
way, all problems are treated in a uniform way. The same holds when building a problem interactively, or when writing
a script, whenever the Problem.solve() function is used for a problem solution.
The following solvers are available:
• ts.adaptive: Implicit time stepping solver with an adaptive time step.
• ts.bathe: Solve elastodynamics problems by the Bathe method.
• ts.central_difference: Solve elastodynamics problems by the explicit central difference method.
• ts.euler: Simple forward euler method
• ts.generalized_alpha: Solve elastodynamics problems by the generalized 𝛼 method.
• ts.multistaged: Explicit time stepping solver with multistage solve_step method
• ts.newmark: Solve elastodynamics problems by the Newmark method.
• ts.runge_kutta_4: Classical 4th order Runge-Kutta method,
• ts.simple: Implicit time stepping solver with a fixed time step.
• ts.stationary: Solver for stationary problems without time stepping.
• ts.tvd_runge_kutta_3: 3rd order Total Variation Diminishing Runge-Kutta method
• ts.velocity_verlet: Solve elastodynamics problems by the explicit velocity-Verlet method.
See sfepy.solvers.ts_solvers for available time-stepping solvers and their options.
The following time step controllers are available:
• tsc.ed_basic: Adaptive time step I-controller for elastodynamics.
• tsc.ed_linear: Adaptive time step controller for elastodynamics and linear problems.
• tsc.ed_pid: Adaptive time step PID controller for elastodynamics.
• tsc.fixed: Fixed (do-nothing) time step controller.
• tsc.time_sequence: Given times sequence time step controller.
See sfepy.solvers.ts_controllers for available time step controllers and their options.

Nonlinear Solvers

Almost every problem, even linear, is solved in SfePy using a nonlinear solver that calls a linear solver in each iteration.
This approach unifies treatment of linear and non-linear problems, and simplifies application of Dirichlet (essential)
boundary conditions, as the linear system computes not a solution, but a solution increment, i.e., it always has zero
boundary conditions.
The following solvers are available:
• nls.newton: Solves a nonlinear system 𝑓 (𝑥) = 0 using the Newton method.
• nls.oseen: The Oseen solver for Navier-Stokes equations.
• nls.petsc: Interface to PETSc SNES (Scalable Nonlinear Equations Solvers).
• nls.scipy_broyden_like: Interface to Broyden and Anderson solvers from scipy.optimize.
• nls.semismooth_newton: The semi-smooth Newton method.

1.4. User’s Guide 45

SfePy Documentation, Release version: 2024.2

See sfepy.solvers.nls, sfepy.solvers.oseen and sfepy.solvers.semismooth_newton for all available

nonlinear solvers and their options.

Linear Solvers

Choosing a suitable linear solver is key to solving efficiently stationary as well as transient PDEs. SfePy allows using
a number of external solvers with a unified interface.
The following solvers are available:
• ls.cholesky: Interface to scikit-sparse.cholesky solver.
• ls.cm_pb: Conjugate multiple problems.
• ls.mumps: Interface to MUMPS solver.
• ls.mumps_par: Interface to MUMPS parallel solver.
• ls.petsc: PETSc Krylov subspace solver.
• ls.pyamg: Interface to PyAMG solvers.
• ls.pyamg_krylov: Interface to PyAMG Krylov solvers.
• ls.rmm: Special solver for explicit transient elastodynamics.
• ls.schur_mumps: Mumps Schur complement solver.
• ls.scipy_direct: Direct sparse solver from SciPy.
• ls.scipy_iterative: Interface to SciPy iterative solvers.
• ls.scipy_superlu: SuperLU - direct sparse solver from SciPy.
• ls.scipy_umfpack: UMFPACK - direct sparse solver from SciPy.
See sfepy.solvers.ls for all available linear solvers and their options.

Virtual Linear Solvers with Automatic Selection

A “virtual” solver can be used in case it is not clear which external linear solvers are available. Each “virtual” solver
selects the first available solver from a pre-defined list.
The following solvers are available:
• ls.auto_direct: The automatically selected linear direct solver.
• ls.auto_iterative: The automatically selected linear iterative solver.
See sfepy.solvers.auto_fallback for all available virtual solvers.

Eigenvalue Problem Solvers

The following eigenvalue problem solvers are available:

• eig.matlab: Matlab eigenvalue problem solver.
• eig.octave: Octave eigenvalue problem solver.
• eig.primme: PRIMME eigenvalue problem solver.
• eig.scipy: SciPy-based solver for both dense and sparse problems.
• eig.scipy_lobpcg: SciPy-based LOBPCG solver for sparse symmetric problems.

46 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

• eig.sgscipy: SciPy-based solver for dense symmetric problems.

• eig.slepc: General SLEPc eigenvalue problem solver.
See sfepy.solvers.eigen for available eigenvalue problem solvers and their options.

Quadratic Eigenvalue Problem Solvers

The following quadratic eigenvalue problem solvers are available:

• eig.qevp: Quadratic eigenvalue problem solver based on the problem linearization.
See sfepy.solvers.qeigen for available quadratic eigenvalue problem solvers and their options.

Optimization Solvers

The following optimization solvers are available:

• nls.scipy_fmin_like: Interface to SciPy optimization solvers scipy.optimize.fmin_*.
• opt.fmin_sd: Steepest descent optimization solver.
See sfepy.solvers.optimize for available optimization solvers and their options.

1.4.10 Solving Problems in Parallel

The PETSc-based nonlinear equations solver 'nls.petsc' and linear system solver 'ls.petsc' can be used for
parallel computations, together with the modules in sfepy.parallel package. This feature is very preliminary, and can
be used only with the commands for interactive use - problem description files are not supported (yet). The key module
is sfepy.parallel.parallel that takes care of the domain and field DOFs distribution among parallel tasks, as well
as parallel assembling to PETSc vectors and matrices.

Current Implementation Drawbacks

• The partitioning of the domain and fields DOFs is not done in parallel and all tasks need to load the whole mesh
and define the global fields - those must fit into memory available to each task.
• While all KSP and SNES solver are supported, in principle, most of their options have to be passed using the
command-line parameters of PETSc - they are not supported yet in the SfePy solver parameters.
• There are no performance statistics yet. The code was tested on a single multi-cpu machine only.
• The global solution is gathered to task 0 and saved to disk serially.
• The vertices of surface region selector does not work in parallel, because the region definition is applied
to a task-local domain.

1.4. User’s Guide 47

SfePy Documentation, Release version: 2024.2

Examples

The examples demonstrating the use parallel problem solving in SfePy are:
• diffusion/poisson_parallel_interactive.py
• multi_physics/biot_parallel_interactive.py
See their help messages for further information.

1.4.11 Isogeometric Analysis

Isogeometric analysis (IGA) is a recently developed computational approach that allows using the NURBS-based do-
main description from CAD design tools also for approximation purposes similar to the finite element method.
The implementation is SfePy is based on Bezier extraction of NURBS as developed in1 . This approach allows reusing
the existing finite element assembling routines, as still the evaluation of weak forms occurs locally in “elements” and
the local contributions are then assembled to the global system.

Current Implementation

The IGA code is still very preliminary and some crucial components are missing. The current implementation is also
very slow, as it is in pure Python.
The following already works:
• single patch tensor product domain support in 2D and 3D
• region selection based on topological Bezier mesh, see below
• Dirichlet boundary conditions using projections for non-constant values
• evaluation in arbitrary point in the physical domain
• both scalar and vector cell terms work
• term integration over the whole domain as well as a cell subdomain
• simple linearization (output file generation) based on sampling the results with uniform parametric vectors
• basic domain generation with sfepy/scripts/gen_iga_patch.py based on igakit
The following is not implemented yet:
• tests
• theoretical convergence rate verification
• facet terms
• other boundary conditions
• proper (adaptive) linearization for post-processing
• support for multiple NURBS patches
1 Michael J. Borden, Michael A. Scott, John A. Evans, Thomas J. R. Hughes: Isogeometric finite element data structures based on Bezier

extraction of NURBS, Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, Texas, March 2010.

48 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Domain Description

The domain description is in custom HDF5-based files with .iga extension. Such a file contains:
• NURBS patch data (knots, degrees, control points and weights). Those can either be generated using igakit,
created manually or imported from other tools.
• Bezier extraction operators and corresponding DOF connectivity (computed by SfePy).
• Bezier mesh control points, weights and connectivity (computed by SfePy).
The Bezier mesh is used to create a topological Bezier mesh - a subset of the Bezier mesh containing the Bezier
element corner vertices only. Those vertices are interpolatory (are on the exact geometry) and so can be used for region
selections.

Region Selection

The domain description files contain vertex sets for regions corresponding to the patch sides, named 'xiIJ', where I
is the parametric axis (0, 1, or 2) and J is 0 or 1 for the beginning and end of the axis knot span. Other regions can be
defined in the usual way, using the topological Bezier mesh entities.

Examples

The examples demonstrating the use of IGA in SfePy are:

• diffusion/poisson_iga.py
• linear_elasticity/linear_elastic_iga.py
• navier_stokes/navier_stokes2d_iga.py
Their problem description files are almost the same as their FEM equivalents, with the following differences:
• There is filename_domain instead of filename_mesh.
• Fields are defined as follows:

fields = {
't1' : ('real', 1, 'Omega', None, 'H1', 'iga'),
't2' : ('real', 1, 'Omega', 'iga', 'H1', 'iga'),
't3' : ('real', 1, 'Omega', 'iga+%d', 'H1', 'iga'),
}

The approximation order in the first definition is None as it is given by the NURBS degrees in the domain
description. The second definition is equivalent to the first one. The third definition, where %d should be a non-
negative integer, illustrates how to increase the field’s NURBS degrees (while keeping the continuity) w.r.t. the
domain NURBS description. It is applied in the navier_stokes/navier_stokes2d_iga.py example to the velocity
field.

1.4. User’s Guide 49

SfePy Documentation, Release version: 2024.2

1.5 Examples

This section contains domain-specific tutorials as well as the automatically generated list of the standard examples that
come with SfePy.

1.5.1 Primer

A beginner’s tutorial highlighting the basics of SfePy.

Introduction

This primer presents a step-by-step walk-through of the process to solve a simple mechanics problem. The typical
process to solve a problem using SfePy is followed: a model is meshed, a problem definition file is drafted, SfePy is run
to solve the problem and finally the results of the analysis are visualised.

Problem statement

A popular test to measure the tensile strength of concrete or asphalt materials is the indirect tensile strength (ITS) test
pictured below. In this test a cylindrical specimen is loaded across its diameter to failure. The test is usually run by
loading the specimen at a constant deformation rate of 50 mm/minute (say) and measuring the load response. When
the tensile stress that develops in the specimen under loading exceeds its tensile strength then the specimen will fail.
To model this problem using finite elements the indirect tensile test can be simplified to represent a diametrically point
loaded disk as shown in the schematic.

The tensile and compressive stresses that develop in the specimen as a result of the point loads P are a function of the
diameter 𝐷 and thickness 𝑡 of the cylindrical specimen. At the centre of the specimen, the compressive stress is 3 times
the tensile stress and the analytical formulation for these are, respectively:

2𝑃
𝜎𝑡 = (1.7)
𝜋𝑡𝐷

50 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

6𝑃
𝜎𝑐 = (1.8)
𝜋𝑡𝐷
These solutions may be approximated using finite element methods. To solve this problem using SfePy the first step is
meshing a suitable model.

Meshing

Assuming plane strain conditions, the indirect tensile test may be modelled using a 2D finite element mesh. Further-
more, the geometry of the model is symmetrical about the x- and y-axes passing through the centre of the circle. To
take advantage of this symmetry only one quarter of the 2D model will be meshed and boundary conditions will be
established to indicate this symmetry. The meshing program Gmsh is used here to very quickly mesh the model. Follow
these steps to model the ITS:
1. The ITS specimen has a diameter of 150 mm. Using Gmsh add three new points (geometry elementary entities)
at the following coordinates: (0.00.0), (75.0, 0.0) and (0.0, 75.0).
2. Next add two straight lines connecting the points.
3. Next add a Circle arc connecting two of the points to form the quarter circle segment.
4. Still under Geometry add a ruled surface.
5. With the geometry of the model defined, add a mesh by clicking on the 2D button under the Mesh functions.
The figures that follow show the various stages in the model process.

1.5. Examples 51
SfePy Documentation, Release version: 2024.2

That’s the meshing done. Save the mesh in a format that SfePy recognizes. For now use the medit .mesh format e.g.
its2D.mesh.
Hint: Check the drop down in the Save As dialog for the different formats that Gmsh can save to.
If you open the its2D.mesh file using a text editor you’ll notice that Gmsh saves the mesh in a 3D format and includes
some extra geometry items that should be deleted. Reformatted the mesh file to a 2D format and delete the Edges block.
Note that when you do this the file cannot be reopened by Gmsh so it is always a good idea to also save your meshes in
Gmsh’s native format as well (Shift-Ctrl-S). Click here to download the reformatted mesh file that will be used in the
tutorial.

You’ll notice that the mesh contains 55 vertices (nodes) and 83 triangle elements. The mesh file provides the coordinates
of the nodes and the element connectivity. It is important to note that node and element numbering in SfePy start at 0
and not 1 as is the case in Gmsh and some other meshing programs.
To view .mesh files you can use a demo of medit. After loading your mesh file with medit you can see the node and
element numbering by pressing P and F respectively. The numbering in medit starts at 1 as shown. Thus the node
at the center of the model in SfePy numbering is 0, and elements 76 and 77 are connected to this node. Node and

52 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

element numbers can also be viewed in Gmsh – under the mesh option under the Visibility tab enable the node and
surface labels. Note that the surface labels as numbered in Gmsh follow on from the line numbering. So to get the
corresponding element number in SfePy you’ll need to subtract the number of lines in the Gmsh file + 1. Confused yet?
Luckily, SfePy provides some useful mesh functions to indicate which elements are connected to which nodes. Nodes
and elements can also be identified by defining regions, which is addressed later.
Another open source python option to view .mesh files is the appropriately named Python Mesh Viewer.
The next step in the process is coding the SfePy problem definition file.

Problem description

The programming of the problem description file is well documented in the SfePy User’s Guide. The problem descrip-
tion file used in the tutorial follows:

r"""
Diametrically point loaded 2-D disk. See :ref:`sec-primer`.

Find :math:`\ul{u}` such that:

.. math::
\int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= 0
\;, \quad \forall \ul{v} \;,

where

.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl}
\;.
"""
from __future__ import absolute_import
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
from sfepy.discrete.fem.utils import refine_mesh
from sfepy import data_dir

# Fix the mesh file name if you run this file outside the SfePy directory.
filename_mesh = data_dir + '/meshes/2d/its2D.mesh'

refinement_level = 0
filename_mesh = refine_mesh(filename_mesh, refinement_level)

output_dir = '.' # set this to a valid directory you have write access to

young = 2000.0 # Young's modulus [MPa]

poisson = 0.4 # Poisson's ratio

options = {
'output_dir' : output_dir,
}

regions = {
'Omega' : 'all',
(continues on next page)

1.5. Examples 53
SfePy Documentation, Release version: 2024.2

(continued from previous page)

'Left' : ('vertices in (x < 0.001)', 'facet'),
'Bottom' : ('vertices in (y < 0.001)', 'facet'),
'Top' : ('vertex 2', 'vertex'),
}

materials = {
'Asphalt' : ({'D': stiffness_from_youngpoisson(2, young, poisson)},),
'Load' : ({'.val' : [0.0, -1000.0]},),
}

fields = {
'displacement': ('real', 'vector', 'Omega', 1),
}

equations = {
'balance_of_forces' :
"""dw_lin_elastic.2.Omega(Asphalt.D, v, u)
= dw_point_load.0.Top(Load.val, v)""",
}

variables = {
'u' : ('unknown field', 'displacement', 0),
'v' : ('test field', 'displacement', 'u'),
}

ebcs = {
'XSym' : ('Bottom', {'u.1' : 0.0}),
'YSym' : ('Left', {'u.0' : 0.0}),
}

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

Download the Problem description file and open it in your favourite Python editor. Note that you may wish to
change the location of the output directory to somewhere on your drive. You may also need to edit the mesh file name.
For the analysis we will assume that the material of the test specimen is linear elastic and isotropic. We define two
material constants i.e. Young’s modulus and Poisson’s ratio. The material is assumed to be asphalt concrete having a
Young’s modulus of 2,000 MPa and a Poisson’s ration of 0.4.
Note: Be consistent in your choice and use of units. In the tutorial we are using Newton (N), millimeters (mm) and
megaPascal (MPa). The sfepy.mechanics.units module might help you in determining which derived units correspond
to given basic units.
The following block of code defines regions on your mesh:

regions = {
'Omega' : 'all',
'Left' : ('vertices in (x < 0.001)', 'facet'),
(continues on next page)

54 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

'Bottom' : ('vertices in (y < 0.001)', 'facet'),
'Top' : ('vertex 2', 'vertex'),
}

Four regions are defined:

1. ‘Omega’: all the elements in the mesh,
2. ‘Left’: the y-axis,
3. ‘Bottom’: the x-axis,
4. ‘Top’: the topmost node. This is where the load is applied.
Having defined the regions these can be used in other parts of your code. For example, in the definition of the boundary
conditions:

ebcs = {
'XSym' : ('Bottom', {'u.1' : 0.0}),
'YSym' : ('Left', {'u.0' : 0.0}),
}

Now the power of the regions entity becomes apparent. To ensure symmetry about the x-axis, the vertical or y-
displacement of the nodes in the ‘Bottom’ region are prevented or set to zero. Similarly, for symmetry about the
y-axis, any horizontal or displacement in the x-direction of the nodes in the ‘Left’ region or y-axis is prevented.
The load is specified in terms of the ‘Load’ material as follows:

materials = {
'Asphalt' : ({
'lam' : lame_from_youngpoisson(young, poisson)[0],
'mu' : lame_from_youngpoisson(young, poisson)[1],
},),
'Load' : ({'.val' : [0.0, -1000.0]},),
}

Note the dot in ‘.val’ – this denotes a special material value, i.e., a value that is not to be evaluated in quadrature points.
The load is then applied in equations using the ‘dw_point_load.0.Top(Load.val, v)’ term in the topmost node (region
‘Top’).
We provided the material constants in terms of Young’s modulus and Poisson’s ratio, but the linear elastic isotropic
equation used requires as input Lamé’s parameters. The lame_from_youngpoisson() function is thus used for conver-
sion. Note that to use this function it was necessary to import the function into the code, which was done up front:

from sfepy.mechanics.matcoefs import lame_from_youngpoisson

Hint: Check out the sfepy.mechanics.matcoefs module for other useful material related functions.
That’s it – we are now ready to solve the problem.

1.5. Examples 55
SfePy Documentation, Release version: 2024.2

Running SfePy

One option to solve the problem is to run the sfepy-run` command: from the command line:

sfepy-run its2D_1.py

Note: For the purpose of this tutorial it is assumed that the problem description file (its2D_1.py) is in the current
directory. If you have the its2D_1.py file in another directory then make sure you include the path to this file as well.
SfePy solves the problem and outputs the solution to the output path (output_dir) provided in the script. The output file
will be in the VTK format by default if this is not explicitly specified and the name of the output file will be the same
as that used for the mesh file except with the ‘.vtk’ extension i.e. its2D.vtk.
The VTK format is an ASCII format. Open the file using a text editor. You’ll notice that the output file includes separate
sections:
• POINTS (these are the model nodes),
• CELLS (the model element connectivity),
• VECTORS (the node displacements in the x-, y- and z- directions).
SfePy provides a script to quickly view the solution. To run this script you need to have pyvista installed. From the
command line issue the following (assuming the correct paths):

sfepy-view its2D.vtk -2 -e

The sfepy-view command generates the image shown below, which shows by default the displacements in the model
as arrows and their magnitude as color scale. Cool, but we are more interested in the stresses. To get these we need to
modify the problem description file and do some post-processing.

Post-processing

SfePy provides functions to calculate stresses and strains. We’ll include a function to calculate these and update the
problem material definition and options to call this function as a post_process_hook(). Save this file as its2D_2.py.

r"""
Diametrically point loaded 2-D disk with postprocessing. See
:ref:`sec-primer`.

Find :math:`\ul{u}` such that:

.. math::
(continues on next page)

56 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

\int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= 0
\;, \quad \forall \ul{v} \;,

where

.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl}
\;.
"""

from __future__ import absolute_import

from sfepy.examples.linear_elasticity.its2D_1 import *

from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson

def stress_strain(out, pb, state, extend=False):

"""
Calculate and output strain and stress for given displacements.
"""
from sfepy.base.base import Struct

ev = pb.evaluate
strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.2.Omega(Asphalt.D, u)', mode='el_avg',
copy_materials=False)

out['cauchy_strain'] = Struct(name='output_data', mode='cell',

data=strain, dofs=None)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress, dofs=None)

return out

asphalt = materials['Asphalt'][0]
asphalt.update({'D' : stiffness_from_youngpoisson(2, young, poisson)})
options.update({'post_process_hook' : 'stress_strain',})

The updated file imports all of the previous definitions in its2D_1.py. The stress function (de_cauchy_stress())
requires as input the stiffness tensor – thus it was necessary to update the materials accordingly. The problem options
were also updated to call the stress_strain() function as a post_process_hook().
Run SfePy to solve the updated problem and view the solution (assuming the correct paths):

sfepy-run its2D_2.py
sfepy-view its2D.vtk -2 --max-plots 2

In addition to the node displacements, the VTK output shown below now also includes the stresses and strains averaged
in the elements:

1.5. Examples 57
SfePy Documentation, Release version: 2024.2

Remember the objective was to determine the stresses at the centre of the specimen under a load 𝑃 . The solution as
currently derived is expressed in terms of a global displacement vector 𝑢. The global (residual) force vector 𝑓 is a
function of the global displacement vector and the global stiffness matrix 𝐾 as: 𝑓 = 𝐾𝑢. Let’s determine the force
vector interactively.

Running SfePy in interactive mode

In addition to solving problems using the sfepy-run command you can also run SfePy interactively (we will use IPython
interactive shell in following examples).
In the SfePy top-level directory run

ipython

issue the following commands:

In [1]: from sfepy.applications import solve_pde

In [2]: pb, variables = solve_pde('its2D_2.py')

The problem is solved and the problem definition and solution are provided in the pb and variables variables respec-
tively. The solution, or in this case, the global displacement vector 𝑢, contains the x- and y-displacements at the nodes
in the 2D model:

In [3]: u = variables()

In [4]: u
Out[4]:
array([ 0. , 0. , 0.37376671, ..., -0.19923848,
0.08820237, -0.11201528])

In [5]: u.shape
Out[5]: (110,)

In [6]: u.shape = (55, 2)

In [7]: u
Out[7]:
array([[ 0. , 0. ],
[ 0.37376671, 0. ],
[ 0. , -1.65318152],
(continues on next page)

58 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

...,
[ 0.08716448, -0.23069047],
[ 0.27741356, -0.19923848],
[ 0.08820237, -0.11201528]])

Note: We have used the fact, that the state vector contains only one variable (u). In general, the following can be used:

In [8]: u = variables.get_state_parts()['u']

In [9]: u
Out[9]:
array([[ 0. , 0. ],
[ 0.37376671, 0. ],
[ 0. , -1.65318152],
...,
[ 0.08716448, -0.23069047],
[ 0.27741356, -0.19923848],
[ 0.08820237, -0.11201528]])

Both variables() and variables.get_state_parts() return a view of the DOF vector, that is why in Out[8] the vector is
reshaped according to Out[6]. It is thus possible to set the values of state variables by manipulating the state vector,
but shape changes such as the one above are not advised (see In [15] below) - work on a copy instead.
From the above it can be seen that u holds the displacements at the 55 nodes in the model and that the displacement
at node 2 (on which the load is applied) is (0, −1.65318152). The global stiffness matrix is saved in pb as a sparse
matrix:

In [10]: K = pb.mtx_a

In [11]: K
Out[11]:
<94x94 sparse matrix of type '<type 'numpy.float64'>'
with 1070 stored elements in Compressed Sparse Row format>

In [12]: print(K)
(0, 0) 2443.95959851
(0, 7) -2110.99917491
(0, 14) -332.960423597
(0, 15) 1428.57142857
(1, 1) 2443.95959852
(1, 13) -2110.99917492
(1, 32) 1428.57142857
(1, 33) -332.960423596
(2, 2) 4048.78343529
(2, 3) -1354.87004384
(2, 52) -609.367453538
(2, 53) -1869.0018791
(2, 92) -357.41672785
(2, 93) 1510.24654193
(3, 2) -1354.87004384
(3, 3) 4121.03202907
(3, 4) -1696.54911732
(3, 48) 76.2400806561
(continues on next page)

1.5. Examples 59
SfePy Documentation, Release version: 2024.2

(continued from previous page)

(3, 49) -1669.59247304
(3, 52) -1145.85294856
(3, 53) 2062.13955556
(4, 3) -1696.54911732
(4, 4) 4410.17902905
(4, 5) -1872.87344838
(4, 42) -130.515009576
: :
(91, 81) -1610.0550578
(91, 86) -199.343680224
(91, 87) -2330.41406097
(91, 90) -575.80373408
(91, 91) 7853.23899229
(92, 2) -357.41672785
(92, 8) 1735.59411191
(92, 50) -464.976034459
(92, 51) -1761.31189004
(92, 52) -3300.45367361
(92, 53) 1574.59387937
(92, 88) -250.325600254
(92, 89) 1334.11823335
(92, 92) 9219.18643706
(92, 93) -2607.52659081
(93, 2) 1510.24654193
(93, 8) -657.361661955
(93, 50) -1761.31189004
(93, 51) 54.1134516246
(93, 52) 1574.59387937
(93, 53) -315.793227627
(93, 88) 1334.11823335
(93, 89) -4348.13351285
(93, 92) -2607.52659081
(93, 93) 9821.16012014

In [13]: K.shape
Out[13]: (94, 94)

One would expect the shape of the global stiffness matrix 𝐾 to be (110, 110) i.e. to have the same number of rows and
columns as u. This matrix has been reduced by the fixed degrees of freedom imposed by the boundary conditions set at
the nodes on symmetry axes. To restore the shape of the matrix, temporarily remove the imposed boundary conditions:

In [14]: pb.remove_bcs()

Note that the matrix is reallocated, so it contains only zeros at this moment.
Now we can calculate the force vector 𝑓 :

In [15]: f = pb.evaluator.eval_residual(u)

This leads to:

ValueError: shape mismatch: value array of shape (55,2) could not be broadcast to␣
˓→indexing result of shape (110,)

60 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

– the original shape of the DOF vector needs to be restored:

In [16]: u.shape = variables.vec.shape = (110,)

In [17]: f = pb.evaluator.eval_residual(u)

In [18]: f.shape
Out[18]: (110,)

In [19]: f
Out[19]:
array([ -4.73618436e+01, 1.42752386e+02, 1.56921124e-13, ...,
-2.06057393e-13, 2.13162821e-14, -2.84217094e-14])

Remember to restore the original boundary conditions previously removed in step [14]:

In [20]: pb.time_update()

To view the residual force vector, we can save it to a VTK file. This requires setting f to (a copy of) the variables as
follows:

In [21]: fvars = variables.copy() # Shallow copy, .vec is shared!

In [22]: fvars.init_state(f) # Initialize new .vec with f.
In [23]: out = fvars.create_output()
In [24]: pb.save_state('file.vtk', out=out)

Running the sfepy-view command on file.vtk, i.e.

In [25]: ! sfepy-view file.vtk

displays the average nodal forces as shown below:

The forces in the x- and y-directions at node 2 are:

In [26]: f.shape = (55, 2)

In [27]: f[2]
Out[27]: array([ 6.20373272e+02, -1.13686838e-13])

Great, we have an almost zero residual vertical load or force apparent at node 2 i.e. -1.13686838e-13 Newton. Note
that these relatively very small numbers can be slightly different, depending on the linear solver kind and version.
Let us now check the stress at node 0, the centre of the specimen.

1.5. Examples 61
SfePy Documentation, Release version: 2024.2

Generating output at element nodes

Previously we had calculated the stresses in the model but these were averaged from those calculated at Gauss quadrature
points within the elements. It is possible to provide custom integrals to allow the calculation of stresses with the Gauss
quadrature points at the element nodes. This will provide us a more accurate estimate of the stress at the centre of the
specimen located at node 0. The code below outlines one way to achieve this.

r"""
Diametrically point loaded 2-D disk with nodal stress calculation. See
:ref:`sec-primer`.

Find :math:`\ul{u}` such that:

.. math::
\int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= 0
\;, \quad \forall \ul{v} \;,

where

.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl}
\;.
"""
from __future__ import print_function
from __future__ import absolute_import
from sfepy.examples.linear_elasticity.its2D_1 import *

from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson

from sfepy.discrete.fem.geometry_element import geometry_data
from sfepy.discrete import FieldVariable
from sfepy.discrete.fem import Field
import numpy as nm

gdata = geometry_data['2_3']
nc = len(gdata.coors)

def nodal_stress(out, pb, state, extend=False, integrals=None):

'''
Calculate stresses at nodal points.
'''

# Point load.
mat = pb.get_materials()['Load']
P = 2.0 * mat.get_data('special', 'val')[1]

# Calculate nodal stress.

pb.time_update()

if integrals is None: integrals = pb.get_integrals()

stress = pb.evaluate('ev_cauchy_stress.ivn.Omega(Asphalt.D, u)', mode='qp',

(continues on next page)

62 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

integrals=integrals, copy_materials=False)
sfield = Field.from_args('stress', nm.float64, (3,),
pb.domain.regions['Omega'])
svar = FieldVariable('sigma', 'parameter', sfield,
primary_var_name='(set-to-None)')
svar.set_from_qp(stress, integrals['ivn'])

print('\n==================================================================')
print('Given load = %.2f N' % -P)
print('\nAnalytical solution')
print('===================')
print('Horizontal tensile stress = %.5e MPa/mm' % (-2.*P/(nm.pi*150.)))
print('Vertical compressive stress = %.5e MPa/mm' % (-6.*P/(nm.pi*150.)))
print('\nFEM solution')
print('============')
print('Horizontal tensile stress = %.5e MPa/mm' % (svar()[0]))
print('Vertical compressive stress = %.5e MPa/mm' % (-svar()[1]))
print('==================================================================')
return out

asphalt = materials['Asphalt'][0]
asphalt.update({'D' : stiffness_from_youngpoisson(2, young, poisson)})
options.update({'post_process_hook' : 'nodal_stress',})

integrals = {
'ivn' : ('custom', gdata.coors, [gdata.volume / nc] * nc),
}

The output:

==================================================================
Given load = 2000.00 N

Analytical solution
===================
Horizontal tensile stress = 8.48826e+00 MPa/mm
Vertical compressive stress = 2.54648e+01 MPa/mm

FEM solution
============
Horizontal tensile stress = 7.57220e+00 MPa/mm
Vertical compressive stress = 2.58660e+01 MPa/mm
==================================================================

Not bad for such a coarse mesh! Re-running the problem using a finer mesh provides a more accurate solution:

==================================================================
Given load = 2000.00 N

Analytical solution
===================
Horizontal tensile stress = 8.48826e+00 MPa/mm
(continues on next page)

1.5. Examples 63
SfePy Documentation, Release version: 2024.2

(continued from previous page)

Vertical compressive stress = 2.54648e+01 MPa/mm

FEM solution
============
Horizontal tensile stress = 8.50042e+00 MPa/mm
Vertical compressive stress = 2.54300e+01 MPa/mm

To see how the FEM solution approaches the analytical one, try to play with the uniform mesh refinement level in the
Problem description file, namely lines 25, 26:

refinement_level = 0
filename_mesh = refine_mesh(filename_mesh, refinement_level)

The above computation could also be done in the IPython shell:

In [23]: from sfepy.applications import solve_pde

In [24]: from sfepy.discrete import (Field, FieldVariable, Material,
Integral, Integrals)
In [25]: from sfepy.discrete.fem.geometry_element import geometry_data

In [26]: gdata = geometry_data['2_3']

In [27]: nc = len(gdata.coors)
In [28]: ivn = Integral('ivn', order=-1,
....: coors=gdata.coors, weights=[gdata.volume / nc] * nc)

In [29]: pb, variables = solve_pde('sfepy/examples/linear_elasticity/its2D_2.py')

In [30]: stress = pb.evaluate('ev_cauchy_stress.ivn.Omega(Asphalt.D,u)',

....: mode='qp', integrals=Integrals([ivn]))
In [31]: sfield = Field.from_args('stress', nm.float64, (3,), pb.domain.regions['Omega'])
In [32]: svar = FieldVariable('sigma', 'parameter', sfield,
....: primary_var_name='(set-to-None)')
In [33]: svar.set_from_qp(stress, ivn)

In [34]: print('Horizontal tensile stress = %.5e MPa/mm' % (svar()[0]))

Horizontal tensile stress = 7.57220e+00 MPa/mm
In [35]: print('Vertical compressive stress = %.5e MPa/mm' % (-svar()[1]))
Vertical compressive stress = 2.58660e+01 MPa/mm

In [36]: mat = pb.get_materials()['Load']

In [37]: P = 2.0 * mat.get_data('special', 'val')[1]
In [38]: P
Out[38]: -2000.0

In [39]: print('Horizontal tensile stress = %.5e MPa/mm' % (-2.*P/(nm.pi*150.)))

Horizontal tensile stress = 8.48826e+00 MPa/mm
In [40]: print('Vertical compressive stress = %.5e MPa/mm' % (-6.*P/(nm.pi*150.)))
Vertical compressive stress = 2.54648e+01 MPa/mm

To wrap this tutorial up let’s explore SfePy’s probing functions.

64 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Probing

As a bonus for sticking to the end of this tutorial see the following Problem description file that provides SfePy
functions to quickly and neatly probe the solution.

r"""
Diametrically point loaded 2-D disk with postprocessing and probes. See
:ref:`sec-primer`.

Find :math:`\ul{u}` such that:

.. math::
\int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= 0
\;, \quad \forall \ul{v} \;,

where

.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl}
\;.
"""
from __future__ import absolute_import
from sfepy.examples.linear_elasticity.its2D_1 import *

from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson

from sfepy.postprocess.probes_vtk import Probe

import os
from six.moves import range

def stress_strain(out, pb, state, extend=False):

"""
Calculate and output strain and stress for given displacements.
"""
from sfepy.base.base import Struct
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm

ev = pb.evaluate
strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.2.Omega(Asphalt.D, u)', mode='el_avg')

out['cauchy_strain'] = Struct(name='output_data', mode='cell',

data=strain, dofs=None)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress, dofs=None)

probe = Probe(out, pb.domain.mesh, probe_view=True)

ps0 = [[0.0, 0.0, 0.0], [ 0.0, 0.0, 0.0]]

ps1 = [[75.0, 0.0, 0.0], [ 0.0, 75.0, 0.0]]
(continues on next page)

1.5. Examples 65
SfePy Documentation, Release version: 2024.2

(continued from previous page)

n_point = 10

labels = ['%s -> %s' % (p0, p1) for p0, p1 in zip(ps0, ps1)]
probes = []
for ip in range(len(ps0)):
p0, p1 = ps0[ip], ps1[ip]
probes.append('line%d' % ip)
probe.add_line_probe('line%d' % ip, p0, p1, n_point)

for ip, label in zip(probes, labels):

fig = plt.figure()
plt.clf()
fig.subplots_adjust(hspace=0.4)
plt.subplot(311)
pars, vals = probe(ip, 'u')
for ic in range(vals.shape[1] - 1):
plt.plot(pars, vals[:,ic], label=r'$u_{%d}$' % (ic + 1),
lw=1, ls='-', marker='+', ms=3)
plt.ylabel('displacements')
plt.xlabel('probe %s' % label, fontsize=8)
plt.legend(loc='best', prop=fm.FontProperties(size=10))

sym_labels = ['11', '22', '12']

plt.subplot(312)
pars, vals = probe(ip, 'cauchy_strain')
for ii in range(vals.shape[1]):
plt.plot(pars, vals[:, ii], label=r'$e_{%s}$' % sym_labels[ii],
lw=1, ls='-', marker='+', ms=3)
plt.ylabel('Cauchy strain')
plt.xlabel('probe %s' % label, fontsize=8)
plt.legend(loc='best', prop=fm.FontProperties(size=8))

plt.subplot(313)
pars, vals = probe(ip, 'cauchy_stress')
for ii in range(vals.shape[1]):
plt.plot(pars, vals[:, ii], label=r'$\sigma_{%s}$' % sym_labels[ii],
lw=1, ls='-', marker='+', ms=3)
plt.ylabel('Cauchy stress')
plt.xlabel('probe %s' % label, fontsize=8)
plt.legend(loc='best', prop=fm.FontProperties(size=8))

opts = pb.conf.options
filename_results = os.path.join(opts.get('output_dir'),
'its2D_probe_%s.png' % ip)

fig.savefig(filename_results)

return out

materials['Asphalt'][0].update({'D' : stiffness_from_youngpoisson(2, young, poisson)})

(continues on next page)

66 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

options.update({
'post_process_hook' : 'stress_strain',
})

Probing applies interpolation to output the solution along specified paths. For the tutorial, line probing is done along
the x- and y-axes of the model.
Run SfePy to solve the problem and apply the probes:

sfepy-run its2D_5.py

The probing function will generate the following figures that show the displacements, normal stresses and strains as
well as shear stresses and strains along the probe paths. Note that you need matplotlib installed to run this example.

The probing function also generates previews of the mesh with the probe paths.

1.5. Examples 67
SfePy Documentation, Release version: 2024.2

Interactive Example

SfePy can be used also interactively by constructing directly the classes that corresponds to the keywords in the problem
description files. The following listing shows a script with the same (and more) functionality as the above examples:

#!/usr/bin/env python
"""
Diametrically point loaded 2-D disk, using commands for interactive use. See
:ref:`sec-primer`.

The script combines the functionality of all the ``its2D_?.py`` examples and
allows setting various simulation parameters, namely:

- material parameters
- displacement field approximation order
- uniform mesh refinement level

The example shows also how to probe the results as in

:ref:`linear_elasticity-its2D_4`. Using :mod:`sfepy.discrete.probes` allows
correct probing of fields with the approximation order greater than one.

In the SfePy top-level directory the following command can be used to get usage
information::

python sfepy/examples/linear_elasticity/its2D_interactive.py -h
"""
from __future__ import absolute_import
(continues on next page)

68 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

import sys
from six.moves import range
sys.path.append('.')
from argparse import ArgumentParser, RawDescriptionHelpFormatter

import numpy as nm
import matplotlib.pyplot as plt

from sfepy.base.base import assert_, output, ordered_iteritems, IndexedStruct

from sfepy.discrete import (FieldVariable, Material, Integral, Integrals,
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.mechanics.matcoefs import stiffness_from_youngpoisson
from sfepy.solvers.auto_fallback import AutoDirect
from sfepy.solvers.nls import Newton
from sfepy.discrete.fem.geometry_element import geometry_data
from sfepy.discrete.probes import LineProbe
from sfepy.discrete.projections import project_by_component

from sfepy.examples.linear_elasticity.its2D_2 import stress_strain

from sfepy.examples.linear_elasticity.its2D_3 import nodal_stress

def gen_lines(problem):
"""
Define two line probes.

Additional probes can be added by appending to `ps0` (start points) and

`ps1` (end points) lists.
"""
ps0 = [[0.0, 0.0], [0.0, 0.0]]
ps1 = [[75.0, 0.0], [0.0, 75.0]]

# Use enough points for higher order approximations.

n_point = 1000

labels = ['%s -> %s' % (p0, p1) for p0, p1 in zip(ps0, ps1)]
probes = []
for ip in range(len(ps0)):
p0, p1 = ps0[ip], ps1[ip]
probes.append(LineProbe(p0, p1, n_point))

return probes, labels

def probe_results(u, strain, stress, probe, label):

"""
Probe the results using the given probe and plot the probed values.
"""
results = {}

pars, vals = probe(u)

(continues on next page)

1.5. Examples 69
SfePy Documentation, Release version: 2024.2

(continued from previous page)

results['u'] = (pars, vals)
pars, vals = probe(strain)
results['cauchy_strain'] = (pars, vals)
pars, vals = probe(stress)
results['cauchy_stress'] = (pars, vals)

fig = plt.figure()
plt.clf()
fig.subplots_adjust(hspace=0.4)
plt.subplot(311)
pars, vals = results['u']
for ic in range(vals.shape[1]):
plt.plot(pars, vals[:,ic], label=r'$u_{%d}$' % (ic + 1),
lw=1, ls='-', marker='+', ms=3)
plt.ylabel('displacements')
plt.xlabel('probe %s' % label, fontsize=8)
plt.legend(loc='best', fontsize=10)

sym_indices = ['11', '22', '12']

plt.subplot(312)
pars, vals = results['cauchy_strain']
for ic in range(vals.shape[1]):
plt.plot(pars, vals[:,ic], label=r'$e_{%s}$' % sym_indices[ic],
lw=1, ls='-', marker='+', ms=3)
plt.ylabel('Cauchy strain')
plt.xlabel('probe %s' % label, fontsize=8)
plt.legend(loc='best', fontsize=10)

plt.subplot(313)
pars, vals = results['cauchy_stress']
for ic in range(vals.shape[1]):
plt.plot(pars, vals[:,ic], label=r'$\sigma_{%s}$' % sym_indices[ic],
lw=1, ls='-', marker='+', ms=3)
plt.ylabel('Cauchy stress')
plt.xlabel('probe %s' % label, fontsize=8)
plt.legend(loc='best', fontsize=10)

return fig, results

helps = {
'young' : "the Young's modulus [default: %(default)s]",
'poisson' : "the Poisson's ratio [default: %(default)s]",
'load' : "the vertical load value (negative means compression)"
" [default: %(default)s]",
'order' : 'displacement field approximation order [default: %(default)s]',
'refine' : 'uniform mesh refinement level [default: %(default)s]',
'probe' : 'probe the results',
}

def main():
from sfepy import data_dir
(continues on next page)

70 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

parser = ArgumentParser(description=__doc__,
formatter_class=RawDescriptionHelpFormatter)
parser.add_argument('--version', action='version', version='%(prog)s')
parser.add_argument('--young', metavar='float', type=float,
action='store', dest='young',
default=2000.0, help=helps['young'])
parser.add_argument('--poisson', metavar='float', type=float,
action='store', dest='poisson',
default=0.4, help=helps['poisson'])
parser.add_argument('--load', metavar='float', type=float,
action='store', dest='load',
default=-1000.0, help=helps['load'])
parser.add_argument('--order', metavar='int', type=int,
action='store', dest='order',
default=1, help=helps['order'])
parser.add_argument('-r', '--refine', metavar='int', type=int,
action='store', dest='refine',
default=0, help=helps['refine'])
parser.add_argument('-p', '--probe',
action="store_true", dest='probe',
default=False, help=helps['probe'])
options = parser.parse_args()

assert_((0.0 < options.poisson < 0.5),

"Poisson's ratio must be in ]0, 0.5[!")
assert_((0 < options.order),
'displacement approximation order must be at least 1!')

output('using values:')
output(" Young's modulus:", options.young)
output(" Poisson's ratio:", options.poisson)
output(' vertical load:', options.load)
output('uniform mesh refinement level:', options.refine)

# Build the problem definition.

mesh = Mesh.from_file(data_dir + '/meshes/2d/its2D.mesh')
domain = FEDomain('domain', mesh)

if options.refine > 0:
for ii in range(options.refine):
output('refine %d...' % ii)
domain = domain.refine()
output('... %d nodes %d elements'
% (domain.shape.n_nod, domain.shape.n_el))

omega = domain.create_region('Omega', 'all')

left = domain.create_region('Left',
'vertices in x < 0.001', 'facet')
bottom = domain.create_region('Bottom',
'vertices in y < 0.001', 'facet')
top = domain.create_region('Top', 'vertex 2', 'vertex')
(continues on next page)

1.5. Examples 71
SfePy Documentation, Release version: 2024.2

(continued from previous page)

field = Field.from_args('fu', nm.float64, 'vector', omega,

approx_order=options.order)

u = FieldVariable('u', 'unknown', field)

v = FieldVariable('v', 'test', field, primary_var_name='u')

D = stiffness_from_youngpoisson(2, options.young, options.poisson)

asphalt = Material('Asphalt', D=D)

load = Material('Load', values={'.val' : [0.0, options.load]})

integral = Integral('i', order=2*options.order)

integral0 = Integral('i', order=0)

t1 = Term.new('dw_lin_elastic(Asphalt.D, v, u)',
integral, omega, Asphalt=asphalt, v=v, u=u)
t2 = Term.new('dw_point_load(Load.val, v)',
integral0, top, Load=load, v=v)
eq = Equation('balance', t1 - t2)
eqs = Equations([eq])

xsym = EssentialBC('XSym', bottom, {'u.1' : 0.0})

ysym = EssentialBC('YSym', left, {'u.0' : 0.0})

ls = AutoDirect({})

nls_status = IndexedStruct()
nls = Newton({}, lin_solver=ls, status=nls_status)

pb = Problem('elasticity', equations=eqs)

pb.set_bcs(ebcs=Conditions([xsym, ysym]))

pb.set_solver(nls)

# Solve the problem.

variables = pb.solve()
output(nls_status)

# Postprocess the solution.

out = variables.create_output()
out = stress_strain(out, pb, variables, extend=True)
pb.save_state('its2D_interactive.vtk', out=out)

gdata = geometry_data['2_3']
nc = len(gdata.coors)

integral_vn = Integral('ivn', coors=gdata.coors,

weights=[gdata.volume / nc] * nc)

nodal_stress(out, pb, variables, integrals=Integrals([integral_vn]))

(continues on next page)

72 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

if options.probe:
# Probe the solution.
probes, labels = gen_lines(pb)

sfield = Field.from_args('sym_tensor', nm.float64, 3, omega,

approx_order=options.order - 1)
stress = FieldVariable('stress', 'parameter', sfield,
primary_var_name='(set-to-None)')
strain = FieldVariable('strain', 'parameter', sfield,
primary_var_name='(set-to-None)')

cfield = Field.from_args('component', nm.float64, 1, omega,

approx_order=options.order - 1)
component = FieldVariable('component', 'parameter', cfield,
primary_var_name='(set-to-None)')

ev = pb.evaluate
order = 2 * (options.order - 1)
strain_qp = ev('ev_cauchy_strain.%d.Omega(u)' % order, mode='qp')
stress_qp = ev('ev_cauchy_stress.%d.Omega(Asphalt.D, u)' % order,
mode='qp', copy_materials=False)

project_by_component(strain, strain_qp, component, order)

project_by_component(stress, stress_qp, component, order)

all_results = []
for ii, probe in enumerate(probes):
fig, results = probe_results(u, strain, stress, probe, labels[ii])

fig.savefig('its2D_interactive_probe_%d.png' % ii)
all_results.append(results)

for ii, results in enumerate(all_results):

output('probe %d:' % ii)
output.level += 2
for key, res in ordered_iteritems(results):
output(key + ':')
val = res[1]
output(' min: %+.2e, mean: %+.2e, max: %+.2e'
% (val.min(), val.mean(), val.max()))
output.level -= 2

if __name__ == '__main__':
main()

The script can be run from the SfePy top-level directory, assuming the in-place build, as follows:

python sfepy/examples/linear_elasticity/its2D_interactive.py

The script allows setting several parameters that influence the solution, see:

1.5. Examples 73
SfePy Documentation, Release version: 2024.2

python sfepy/examples/linear_elasticity/its2D_interactive.py -h

for the complete list. Besides the material parameters, a uniform mesh refinement level and the displacement field
approximation order can be specified. The script demonstrates how to
• project a derived quantity, that is evaluated in quadrature points (e.g. a strain or stress), into a field variable;
• probe the solution defined in the field variables.
Using sfepy.discrete.probes allows correct probing of fields with the approximation order greater than one.
The end.

1.5.2 Using Salome with SfePy

NOTE This tutorial was created in 2014, so it may be obsolete.

Introduction

Salome is a powerful open-source tool for generating meshes for numerical simulation and post processing the results.
This is a short tutorial on using Salome as a preprocessor for preparing meshes for use with SfePy.

Tutorial prerequisites

This tutorial assumes that you have a working copy of Salome. It is possible to build Salome from source code.
Fortunately, for the less brave, many pre-compiled binaries for different platforms are available at the Salome download
page. Registration for a free account may be required to download from the preceding site.
In addition, this tutorial assumes you have a working copy of SfePy with MED read support. See the Installation for
help. Note that it is not actually necessary to “install” SfePy; one may run the code from the source directory (see
notation below) after compilation of the C extension modules (again, see the installation notes if you are confused).

Note on notation used in this tutorial

We are using the following notations:

• <sfepy_root>: the root directory of the SfePy source code
• <work_dir>: the working directory where you plan to save your files

Step 1: Using Salome

Salome has its own set of tutorials and community resources. It is suggested you look around on Salome web site to
familiarize yourself with the available resources.
This tutorial follows the EDF Exercise 1 available from the Salome Tutorial Site. Go ahead and complete this tutorial
now. We will use the result from there in the following.
This is the mesh you should end up with:

74 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Step 2: Exporting mesh from Salome

In the Salome MESH module, right click on the mesh object Mesh_Partition_Hexa you created in the Salome EDF
Exercise 1 Tutorial and click Export to MED file. Save the file as Mesh_Partition_Hexa.med in your working
directory <work_dir>.

Step 3: Copy SfePy project description files

In this tutorial, we will assume that we need to solve a linear elasticity problem on the mesh generated by Salome. Since
the Salome mesh looks a bit like a fish, we will try to simulate the fish waving its tail.
Copy the file <sfepy_root>/sfepy/examples/linear_elasticity/linear_elastic.py to <work_dir>. Use
your favorite python editor to load this file. We will customize this file for our purposes.

Step 4: Modify linear_elastic.py

Mesh specification

The first thing we have to do is tell SfePy to use our new mesh. Change the line

filename_mesh = data_dir + '/meshes/3d/cylinder.mesh'

to

filename_mesh = 'Mesh_Partition_Hexa.med'

1.5. Examples 75
SfePy Documentation, Release version: 2024.2

Region specification

Next, we have to define sensible Regions for the mesh. We will apply a displacement to the Tail and keep the Top and
Bottom of the fish fixed. Change the lines

regions = {
'Omega' : 'all',
'Left' : ('vertices in (x < 0.001)', 'facet'),
'Right' : ('vertices in (x > 0.099)', 'facet'),
'SomewhereTop' : ('vertices in (z > 0.017) & (x > 0.03) & (x < 0.07)', 'vertex'),
}

to

regions = {
'Omega' : 'all',
'Tail' : ('vertices in (x < -94)', 'facet'),
'TopFixed' : ('vertices in (z > 9.999) & (x > 54)', 'facet'),
'BotFixed' : ('vertices in (z < 0.001) & (x > 54)', 'facet'),
}

Field specification

The Salome mesh uses hexahedral linear order elements; in SfePy notation these are called 3_8, see User’s Guide.
Just keep the lines

fields = {
'displacement': ('real', 'vector', 'Omega', 1),
}

Boundary condition specifications

In this section, we tell SfePy to fix the top and bottom parts of the “head” of the fish and move the tail 10 units to the
side (z direction).
Change the lines

ebcs = {
'Fixed' : ('Left', {'u.all' : 0.0}),
'Displaced' : ('Right', {'u.0' : 0.01, 'u.[1,2]' : 0.0}),
'PerturbedSurface' : ('SomewhereTop', {'u.2' : 0.005}),
}

to

ebcs = {
'TopFixed' : ('TopFixed', {'u.all' : 0.0}),
'BotFixed' : ('BotFixed', {'u.all' : 0.0}),
'Displaced' : ('Tail', {'u.2' : 10, 'u.[0,1]' : 0.0}),
}

76 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Step 5: Run SfePy

Save your changes to linear_elastic.py. Now it’s time to run the SfePy calculation. In your <work_dir> in your
terminal type:

sfepy-run linear_elastic.py

This will run the SfePy calculation. Some progress information is printed to your screen and the residual (a measure
of the convergence of the solution) is printed for each iteration of the solver. The solver terminates when this residual
is less than a certain value. It should only take 1 iteration since we are solving a linear problem. The results will be
saved to Mesh_Partition_Hexa.vtk.
Now we can view the results of our work. In your terminal, type:

sfepy-view Mesh_Partition_Hexa.vtk -f u:wu:f2.0:p0 0:vw:p0

You should get the plot with the deformed and undeformed meshs. Notice how the fish is bending its tail in response
to the applied displacement.
Now you should be able to use meshes created in Salome with SfePy!

1.5.3 Preprocessing: FreeCAD/OpenSCAD + Gmsh

Introduction

There are several open source tools for preparing 2D and 3D finite element meshes like Salome, FreeCAD, Gmsh,
Netgen, etc. Most of them are GUI based geometrical modeling and meshing environments/tools but they also usually
allow using their libraries in user scripts. Some of the above mentioned tools are handy for solid modeling, some
of them are great for meshing. This tutorial shows how to combine solid geometry modeling functions provided by
FreeCAD or OpenSCAD with meshing functions of Gmsh.
The collaboration of modeling, meshing and conversion tools and the workflow are illustrated in the following scheme.

1.5. Examples 77
 SfePy Documentation, Release version: 2024.2

Creating geometry using FreeCAD

Functionalities of FreeCAD are accessible to Python and can be used to define geometrical models in simple Python
scripts. There is a tutorial related to Python scripting in FreeCAD.
The first step in creating a Python script is to set up a path to the FreeCAD libraries and import all required modules:

1 import sys
2 FREECADPATH = '/usr/lib/freecad/lib/'
3 sys.path.append(FREECADPATH)
4

5 from FreeCAD import Base, newDocument

6 import Part
7 import Draft
8 import ProfileLib.RegularPolygon as Poly

Now, a new empty FreeCAD document can be defined as:

doc = newDocument()

All new objects describing the geometry will be added to this document.
In the following lines a geometrical model of a screwdriver handle will be created. Let’s start by defining a sphere and
a cylinder and join these objects into the one called uni:

1 radius = 0.01
2 height = 0.1
3

4 cyl = doc.addObject("Part::Cylinder", "cyl")

5 cyl.Radius = radius
6 cyl.Height = height
7

8 sph = doc.addObject("Part::Sphere", "sph")

9 sph.Radius = radius
10

11 uni = doc.addObject("Part::MultiFuse", "uni")

12 uni.Shapes = [cyl, sph]

Create a polygon, revolve it around the z-axis to create a solid and use the result as the cutting tool applied to uni
object:

1 ske = doc.addObject('Sketcher::SketchObject', 'Sketch')

2 ske.Placement = Base.Placement(Base.Vector(0, 0, 0),
3 Base.Rotation(-0.707107, 0, 0, -0.707107))
4 Poly.makeRegularPolygon('Sketch', 5,
5 Base.Vector(-1.2 * radius, 0.9 * height, 0),
6 Base.Vector(-0.8 * radius, 0.9 * height, 0))
7

8 cut = doc.addObject("PartDesign::Revolution", "Revolution")

9 cut.Sketch = ske
10 cut.ReferenceAxis = (ske, ['V_Axis'])
11 cut.Angle = 360.0
12

13 dif = doc.addObject("Part::Cut", "dif")

(continues on next page)

78 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

14 dif.Base = uni
15 dif.Tool = cut

Create a cylinder, make a polar array of the cylinder objects and subtract it from the previous result:

1 cyl1 = doc.addObject("Part::Cylinder", "cyl1")

2 cyl1.Radius = 0.2 * radius
3 cyl1.Height = 1.1 * height
4 cyl1.Placement = Base.Placement(Base.Vector(-1.1 * radius, 0, -0.2 * height),
5 Base.Rotation(0, 0, 0, 1))
6

7 arr = Draft.makeArray(cyl1, Base.Vector(1, 0, 0), Base.Vector(0, 1, 0), 2, 2)

8 arr.ArrayType = "polar"
9 arr.NumberPolar = 6
10

11 dif2 = doc.addObject("Part::Cut", "dif2")

12 dif2.Base = dif
13 dif2.Tool = arr

Create a middle hole for the screwdriver metal part:

1 cyl2 = doc.addObject("Part::Cylinder", "cyl2")

2 cyl2.Radius = 0.3 * radius
3 cyl2.Height = height
4

5 dif3 = doc.addObject("Part::Cut", "dif3")

6 dif3.Base = dif2
7 dif3.Tool = cyl2

Finally, recompute the geometry, export the part to the STEP file and save the document in FreeCAD format (not really
needed for subsequent mesh generation, but may be useful for visualization and geometry check):

1 doc.recompute()
2

3 Part.export([dif3], 'screwdriver_handle.step')
4

5 doc.saveAs('screwdriver_handle.FCStd')

A finite element mesh can be generated directly in FreeCAD using MeshPart module:

1 import MeshPart
2

3 mesh = doc.addObject("Mesh::Feature", "Mesh")

4 mesh.Mesh = MeshPart.meshFromShape(Shape=dif3.Shape, MaxLength=0.002)
5 mesh.Mesh.write("./screwdriver_handle.bdf", "NAS", "mesh")

The meshing function of MeshPart module is limited to triangular grids so it is better to use Gmsh mesh generator
which can provide triangular and quadrilateral meshes in 2D or tetrahedral and hexahedral meshes in 3D. Gmsh allows
to control the meshing process through a wide range of parameters. Meshing by Gmsh will be described in section
Gmsh - generating finite element mesh.

1.5. Examples 79
 SfePy Documentation, Release version: 2024.2

The example of screwdriver handle: screwdriver_handle.py.

There are two simple ways how to discover Python calls of FreeCAD functions. You can enable “show script
commands in python console” in Edit->Preferences->General->Macro and the Python console by selecting
View->Views->Python Console and all subsequent operations will be printed in the console as the Python code.
The second way is to switch on the macro recording function (Macro->Macro recording ...) which generates a
Python script (FCMacro file) containing all the code related to actions in the FreeCAD graphical interface.

Creating geometry using OpenSCAD

The alternative tool for solid geometrical modeling is OpenSCAD - “The Programmers Solid 3D CAD Modeller”.
It has its own description language based on functional programming that is used to construct solid models using
geometrical primitives similar to FreeCAD. Solid geometries can be exported to several file formats including STL and
CSG. OpenSCAD allows solid modeling based on Constructive Solid Geometry (CSG) principles and extrusion of 2D
objects into 3D. The model of a screwdriver handle presented in the previous section can be defined in OpenSCAD by
the following code (screwdriver_handle.scad):
1 radius = 0.01;
2 height = 0.1;
3 $fn = 50;
4

5 difference() {
6 difference() {
7 difference() {
8 union() {
9 cylinder(center=false, h=height, r=radius);
10 sphere(radius);
11 };
(continues on next page)

80 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

12 translate([0, 0, 0.9*height])
13 rotate_extrude()
14 polygon([[0.8*radius, 0], [1.8*radius, -0.577*radius], [1.8*radius, 0.
˓→577*radius]]);

15 }
16 cylinder(center=false, h=1.1*height, r=0.3*radius);
17 }
18 for (i = [1:6]) {
19 rotate([0, 0, 360/6*i])
20 translate([-1.1*radius, 0.0, -0.2*height])
21 cylinder(center=false, h=1.1*height, r=0.2*radius);
22 }
23 }

To generate a finite element mesh of the solid geometry the model must be exported to a suitable file format. OpenSCAD
has limited export options, but by using FreeCAD import/export functions, it is possible to find a workaround. The
OpenSCAD model can be exported to the CSG file format and FreeCAD can be used as a mesh converter to the STEP
format:

1 import sys
2 sys.path.append('/usr/lib/freecad/lib/')
3 sys.path.append('/usr/lib/freecad/Mod/OpenSCAD/')
4

5 import FreeCAD
6 import Part
7 import importCSG
8

9 importCSG.open('screwdriver_handle.csg')
10 Part.export([FreeCAD.ActiveDocument.Objects[-1]], 'screwdriver_handle.step')

1.5. Examples 81
 SfePy Documentation, Release version: 2024.2

Gmsh - generating finite element mesh

Gmsh can create finite element meshes using geometrical models imported from STEP, IGES and BRep files (has to
be compiled with OpenCASCADE support).
The following GEO file imports screwdriver_handle.step file and defines a field controlling the mesh size
(screwdriver_handle.geo):

1 Merge "screwdriver_handle.step";
2

3 Field[1] = MathEval;
4 Field[1].F = "0.002";
5 Background Field = 1;

Now, run Gmsh generator and export the mesh into the MSH format in which all surface and volumetric elements are
stored:

gmsh -3 -format msh -o screwdriver_handle.msh screwdriver_handle.geo

By converting the MSH file into the VTK format using sfepy-convert:

sfepy-convert -d 3 screwdriver_handle.msh screwdriver_handle.vtk

the surface elements are discarded and only the volumetric mesh is preserved.

Note: planar 2D meshes

To create a planar 2D mesh, such as

82 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

that can be described by this Gmsh code, the mesh generator can be called as follows:

gmsh -2 -format msh -o circle_in_square.msh circle_in_square.geo

This, however is not enough to create a truly 2D mesh - the created mesh vertices still have the third, 𝑧, component
which is equal to zero. In order to remove the third component, use:

sfepy-convert --2d circle_in_square.msh circle_in_square.h5

Now, in the resulting circle_in_square.h5, each vertex has only two coordinates. Another way of generating the
2D mesh is to use the legacy VTK format as follows:

gmsh -2 -format vtk -o circle_in_square.vtk circle_in_square.geo

sfepy-convert circle_in_square.vtk circle_in_square.h5

This is due to the fact that the legacy VTK does not support 2D vertices and so the VTKMeshIO reader tries to detect
the planar geometry by comparing the 𝑧 components to zero - the --2d option of sfepy-convert is not needed in
this case.

Multipart models

Meshing models composed of parts with different material groups is a little bit tricky task. But there are some more or
less general ways of doing that. Here, the method using functions of Gmsh for periodic meshes will be shown.
The screwdriver handle example is extended by adding a screwdriver shank. The new part is composed of a cylinder
trimmed at one end:

1 cyl3 = doc.addObject("Part::Cylinder", "cyl3")

2 cyl3.Radius = 0.3 * radius
3 cyl3.Height = height
4 cyl3.Placement = Base.Placement(Base.Vector(0, 0, height),
5 Base.Rotation(0, 0, 0, 1))
6

7 tip1 = doc.addObject("Part::Box", "tip1")

8 tip1.Length = radius
9 tip1.Width = 2 * radius
10 tip1.Height = 3 * radius
11 tip1.Placement = Base.Placement(Base.Vector(0, -radius, 1.71 * height),
12 Base.Rotation(Base.Vector(0, 1, 0), -10),
13 Base.Vector(0, 0, 3 * radius))
14

15 tip2 = doc.addObject("Part::Mirroring", "tip2")

16 tip2.Source = tip1
17 tip2.Normal = (1, 0, 0)
18

19 tip3 = doc.addObject("Part::MultiFuse", "tip3")

20 tip3.Shapes = [tip1, tip2]
21

22 dif4 = doc.addObject("Part::Cut", "dif4")

23 dif4.Base = cyl3
24 dif4.Tool = tip3
25

26 uni2 = doc.addObject("Part::MultiFuse", "uni2")

27 uni2.Shapes = [cyl2, dif4]

1.5. Examples 83
 SfePy Documentation, Release version: 2024.2

The handle and shank are exported to the STEP file as two separated parts:

1 doc.recompute()
2

3 Part.export([dif3, uni2], 'screwdriver_full.step')

4 doc.saveAs('screwdriver_full.FCStd')

The full screwdriver example (handle + shank): screwdriver_full.py.

To create a coincidence mesh on the handle and shank interface, it is necessary to identify the interface surfaces and
declare them to be periodic in the GEO file. The identification has to be done manually in the Gmsh graphical interface.

The input file for Gmsh is than as follows (screwdriver_full.geo):

1 Merge "screwdriver_full.step";
2

3 Periodic Surface 5 {7} = 26 {67};

4 Periodic Surface 3 {6, 2, -6, 7} = 27 {68, 69, -68, 67};
5

6 Physical Volume(1) = {1};

7 Physical Volume(2) = {2};
(continues on next page)

84 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

8

9 Field[1] = MathEval;
10 Field[1].F = "0.0015";
11 Background Field = 1;

where the first pair of periodic surfaces corresponds to the common circle faces (bottom of the shank) and the second
pair to the common cylindrical surfaces. See Gmsh Reference manual for details on periodic meshing.
Using the above stated GEO file, Gmsh creates a mesh containing duplicate vertices on the handle/shank interface.
These duplicate vertices can be removed during the conversion to the VTK format by giving --merge (or just -m)
argument to convert_mesh.py script:

sfepy-convert -m screwdriver_full.msh screwdriver_full.vtk

In order to extract the cells by the physical groups use the conversion script with --save-per-mat argument:
sfepy-convert --save-per-mat screwdriver_full.vtk screwdriver.vtk

It produces screwdriver.vtk contaning the original mesh and screwdriver_matid_1.vtk, screwdriver_matid_2.vtk files
containing only the cells of a given physical group and all vertices of the original mesh.

When using OpenSCAD, define the full screwdriver geometry as (screwdriver_full.scad):

1 radius = 0.01;
2 height = 0.1;
3 $fn = 50;
4

5 module tip() {
6 rotate([0, -10, 0])
7 translate([0, -radius, -3*radius])
8 cube([radius, 2*radius, 3*radius], center=false);
9 }
10

11 difference() {
12 difference() {
13 difference() {
14 union() {
(continues on next page)

1.5. Examples 85
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

15 cylinder(center=false, h=height, r=radius);
16 sphere(radius);
17 };
18 translate([0, 0, 0.9*height])
19 rotate_extrude()
20 polygon([[0.8*radius, 0], [1.8*radius, -0.577*radius], [1.8*radius, 0.
˓→577*radius]]);

21 }
22 cylinder(center=false, h=height, r=0.3*radius);
23 }
24 for (i = [1:6]) {
25 rotate([0, 0, 360/6*i])
26 translate([-1.1*radius, 0.0, -0.2*height])
27 cylinder(center=false, h=1.1*height, r=0.2*radius);
28 }
29 }
30

31 union() {
32 difference() {
33 translate([0, 0, height])
34 cylinder(center=false, h=height, r=0.3*radius);
35 translate([0, 0, 1.71*height + 3*radius])
36 union() {
37 tip();
38 mirror ([1, 0, 0]) tip();
39 }
40 }
41 cylinder(center=false, h=height, r=0.3*radius);
42 }

and convert the CSG file to the STEP file by:

1 importCSG.open('screwdriver_full.csg')
2 top_group = FreeCAD.ActiveDocument.Objects[-1]
3 Part.export(top_group.OutList, 'screwdriver_full.step')

Since the different tools for geometry definition have been used, the numbering of geometric objects may differ and the
surface and edge numbers have to be changed in the GEO file:

Periodic Surface 5 {6} = 26 {66};

Periodic Surface 3 {5, 2, -5, 6} = 27 {67, 68, -67, 66};

Note: The numbering of objects may vary between FreeCAD, OpenSCAD and Gmsh versions.

86 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

1.5.4 Material Identification

Introduction

This tutorial shows identification of material parameters of a composite structure using data (force-displacement curves)
obtained by a standard tensile test.

Composite structure

The unidirectional long fiber carbon-epoxy composite is considered. Its microstructure was analysed by the scanning
electron microscopy and the data, volume fractions and fibers cross-sections, were used to generate a periodic finite
element mesh (representative volume element - RVE) representing the random composite structure at the microscopic
level (the random structure generation algorithm is described in1 ):

This RVE is used in the micromechanical FE analysis which is based on the two-scale homogenization method.

Material testing

Several carbon-expoxy specimens with different fiber orientations (0, 30, 60 and 90 degrees) were subjected to the
tensile test in order to obtain force-elongation dependencies, see2 . The slopes of the linearized dependencies were used
in an objective function of the identification process.
1 Lubachevsky B. D., How to Simulate Billiards and Similar Systems, Journal of Computational Physics, 94(2), 1991. http://arxiv.org/PS_cache/

cond-mat/pdf/0503/0503627v2.pdf
2 Srbová H., Kroupa T., Zemčík R., Identification of the Material Parameters of a Unidirectional Fiber Composite Using a Micromodel, Materiali

in Tehnologije, 46(5), 2012, 431-434.

1.5. Examples 87
SfePy Documentation, Release version: 2024.2

Numerical simulation

The linear isotropic material model is used for both components (fiber and matrix) of the composite so only four material
parameters (Young’s modulus and Poisson’s ratio for each component) are necessary to fully describe the mechanical
behavior of the structure.
The numerical simulations of the tensile tests are based on the homogenization method applied to the linear elastic prob-
lem3 . The homogenization procedure results in the microscopic problem solved within the RVE and the macroscopic
problem that involves the homogenized elastic coefficients.

3 Pinho-da-Cruz L., Oliveira J. A. and Teixeira-Dias F., Asymptotic homogenization in linear elasticity. Part I: Mathematical formulation and

finite element modeling, Computational Materials Science, 45(4), 2009, 1073–1080.

88 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

Homogenized coefficients

The problem at the microscopic level is formulated in terms of characteristic response functions and its solution is used
to evaluate the homogenized elasticity tensor. The microscopic problem has to be solved with the periodic boundary
conditions.
The following SfePy description file is used for definition of the microscopic problem: homogenization_opt_src.
In the case of the identification process function get_mat() obtains the material parameters (Young’s modules, Poisson’s
ratios) from the outer identification loop. Otherwise these parameters are given by values.
Notice the use of parametric_hook (Miscellaneous) to pass around the optimization parameters.

Macroscopic simulation

The homogenized elasticity problem is solved for the unknown macroscopic displacements and the elongation of the
composite specimen is evaluated for a given loading. These values are used to determine the slopes of the calculated
force-elongation dependencies which are required by the objective function.
The SfePy description file for the macroscopic analysis: linear_elasticity_opt_src.

Identification procedure

The identification of material parameters, i.e. the Young’s modulus and Poisson’s ratio, of the epoxy matrix (𝐸𝑚 , 𝜈𝑚 )
and carbon fibers (𝐸𝑓 , 𝜈𝑓 ) can be formulated as a minimization of the following objective function:
(︃ )︃2
𝑖
∑︁ 𝑘𝑐𝑜𝑚𝑝 (x)
Φ(x) = 1− 𝑖
, (1.9)
𝑘𝑒𝑥𝑝
𝑖∈{0,30,60,90}

where 𝑘𝑐𝑜𝑚𝑝
𝑖
and 𝑘𝑒𝑥𝑝
𝑖
are the computed and measured slopes of the force-elongation tangent lines for a given fiber ori-
entation. This function is minimized using scipy.optimize.fmin_tnc(), considering bounds of the identified parameters.
Tho following steps are performed in each iteration of the optimization loop:
1. Solution of the microscopic problem, evaluation of the homogenized elasticity tensor.
2. Solution of the macroscopic problems for different fiber orientations (0, 30, 60, 90), this is incorporated by
appropriate rotation of the elasticity tensor.
3. Evaluation of the objective function.
Python script for material identification: material_opt_src.

Running identification script

Run the script from the command shell as (from the top-level directory of SfePy):

$ python sfepy/examples/homogenization/material_opt.py

The iteration process is monitored using graphs where the values of the objective function and material parameters are
plotted.

1.5. Examples 89
SfePy Documentation, Release version: 2024.2

The resulting values of 𝐸𝑓 , 𝜈𝑓 , 𝐸𝑚 , 𝜈𝑚 can be found at the end of the script output:

>>> material optimization FINISHED <<<

material_opt_micro: terminated
optimized parameters: [1.71129526e+11 3.20844131e-01 2.33507829e+09 2.00000000e-01]

So that:
𝐸𝑓 = 171.13 GPa
𝜈𝑓 = 3.21
𝐸𝑚 = 2.34 GPa
𝜈𝑚 = 0.20
Note: The results may vary across SciPy versions and related libraries.

1.5.5 Mesh parametrization

Introduction

When dealing with shape optimization we usually need to modify a FE mesh using a few optimization parameters
describing the mesh geometry. The B-spline parametrization offers an efficient way to do that. A mesh region (2D or
3D) that is to be parametrized is enclosed in the so called spline-box and the positions of all vertices inside the box can
be changed by moving the control points of the B-spline curves.
There are two different classes for the B-spline parametrization implemented in SfePy (module sfepy.mesh.
splinebox): SplineBox and SplineRegion2D. The first one defines a rectangular parametrization box in 2D or
3D while the second one allows to set up an arbitrary shaped region of parametrization in 2D.

90 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

SplineBox

The rectangular B-spline parametrization is created as follows:

1 from sfepy.mesh.splinebox import SplineBox

2

3 spb = SplineBox(<bbox>, <coors>, <nsg>)

the first parameter defines the range of the box in each dimension, the second parameter is the array of coordinates
(vertices) to be parametrized and the last one (optional) determines the number of control points in each dimension.
The number of the control points (𝑛𝑐𝑝) is calculated as:

𝑛𝑐𝑝𝑖 = 𝑛𝑠𝑔𝑖 + 𝑑𝑒𝑔𝑟𝑒𝑒, 𝑖 = 1, 2(, 3) (1.10)

where 𝑑𝑒𝑔𝑟𝑒𝑒 is the degree of the B-spline curve (default value: 3 = cubic spline) and 𝑛𝑠𝑔 is the number of the spline
segments (default value: [1,1(,1)] = 4 control points for all dimensions).
The position of the vertices can be modified by moving the control points:

spb.move_control_point(<cpoint>, <val>)

where <cpoint> is the index or position of the control point, for explanation see the following figure.

The displacement is given by <val>. The modified coordinates of the vertices are evaluated by:

new_coors = spb.evaluate()

Example

• Create a new 2D SplineBox with the left bottom corner at [-1,-1] and the right top corner at [1, 0.6] which has
5 control points in x-direction and 4 control points in y-direction:

1 from sfepy.mesh.splinebox import SplineBox

2 from sfepy.discrete.fem import Mesh
3

4 mesh = Mesh.from_file('meshes/2d/square_tri1.mesh')
5 spb = SplineBox([[-1, 1], [-1, 0.6]], mesh.coors, nsg=[2,1])

• Modify the position of mesh coordinates by moving three control points (with indices 1,2 and 3):

1.5. Examples 91
 SfePy Documentation, Release version: 2024.2

1 spb.move_control_point(1, [0.1, -0.2])

2 spb.move_control_point(2, [0.2, -0.3])
3 spb.move_control_point(3, [0.0, -0.1])

• Evaluate the new coordinates:

mesh.cmesh.coors[:] = spb.evaluate()

• Write the deformed mesh and the spline control net (the net of control points) into vtk files:

spb.write_control_net('square_tri1_spbox.vtk')
mesh.write('square_tri1_deform.vtk')

The following figures show the undeformed (left) and deformed (right) mesh and the control net.

SplineRegion2D

In this case, the region (only in 2D) of parametrization is defined by four B-spline curves:

1 from sfepy.mesh.splinebox import SplineRegion2D

2

3 spb = SplineRegion2D([<bspl1>, <bspl2>, <bspl3>, <bspl4>], <coors>)

The curves must form a closed loop, must be oriented counterclockwise and the opposite curves (<bspl1>, <bspl3>
and <bspl2>, <bspl4>) must have the same number of control points and the same knot vectors, see the figure below,
on the left.

92 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

The position of the selected vertices, depicted in the figure on the right, are driven by the control points in the same
way as explained above for SplineBox.
Note: Initializing SplineRegion2D may be time consuming due to the fact that for all vertex coordinates the spline
parameters have to be found using an optimization method in which the B-spline basis is repeatedly evaluated.

Example

• First of all, define four B-spline curves (the default degree of the spline curve is 3) representing the boundary of
a parametrization area:

1 from sfepy.mesh.bspline import BSpline

2

3 # left / right boundary

4 line_l = nm.array([[-1, 1], [-1, .5], [-1, 0], [-1, -.5]])
5 line_r = nm.array([[0, -.2], [.1, .2], [.3, .6], [.4, 1]])
6

7 sp_l = BSpline()
8 sp_l.approximate(line_l, ncp=4)
9 kn_lr = sp_l.get_knot_vector()
10

11 sp_r = BSpline()
12 sp_r.approximate(line_r, knots=kn_lr)
13

14 # bottom / top boundary

15 line_b = nm.array([[-1, -.5], [-.8, -.6], [-.5, -.4], [-.2, -.2], [0, -.2]])
16 line_t = nm.array([[.4, 1], [0, 1], [-.2, 1], [-.6, 1], [-1, 1]])
17

18 sp_b = BSpline()
19 sp_b.approximate(line_b, ncp=5)
20 kn_bt = sp_b.get_knot_vector()
21

22 sp_t = BSpline()
23 sp_t.approximate(line_t, knots=kn_bt)

• Create a new 2D SplineRegion2D object:

1.5. Examples 93
SfePy Documentation, Release version: 2024.2

1 from sfepy.mesh.splinebox import SplineRegion2D

2

3 spb = SplineRegion2D([sp_b, sp_r, sp_t, sp_l], mesh.coors)

• Move the control points:

1 spb.move_control_point(5, [-.2, .1])

2 spb.move_control_point(10, [-.3, .2])
3 spb.move_control_point(15, [-.1, .2])

• Evaluate the new coordinates:

mesh.cmesh.coors[:] = spb.evaluate()

The figures below show the undeformed (left) and deformed (right) mesh and the control net.

1.5.6 Examples

acoustics

acoustics/acoustics.py

Description
Acoustic pressure distribution.
This example shows how to solve a problem in complex numbers, note the ‘accoustic_pressure’ field definition.
Find 𝑝 such that:
∫︁ ∫︁ ∫︁ ∫︁
𝑐2 ∇𝑞 · ∇𝑝 − 𝑤2 𝑞𝑝 − 𝑖𝑤𝑐 𝑞𝑝 = 𝑖𝑤𝑐2 𝜌𝑣𝑛 𝑞, ∀𝑞 .
Ω Ω Γ𝑜𝑢𝑡 Γ𝑖𝑛

94 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

source code

r"""
Acoustic pressure distribution.

This example shows how to solve a problem in complex numbers, note the
'accoustic_pressure' field definition.

Find :math:`p` such that:

.. math::
c^2 \int_{\Omega} \nabla q \cdot \nabla p
- w^2 \int_{\Omega} q p
- i w c \int_{\Gamma_{out}} q p
= i w c^2 \rho v_n \int_{\Gamma_{in}} q
\;, \quad \forall q \;.
"""
from __future__ import absolute_import
from sfepy import data_dir

filename_mesh = data_dir + '/meshes/2d/special/two_rectangles.mesh'

v_n = 1.0 # m/s

(continues on next page)

1.5. Examples 95
SfePy Documentation, Release version: 2024.2

(continued from previous page)

w = 1000.0
c = 343.0 # m/s
rho = 1.55 # kg/m^3

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

materials = {
'one' : ({'one' : 1.0},),
}

regions = {
'Omega' : 'all',
'Gamma_in' : ('vertices in (x < 0.01)', 'facet'),
'Gamma_out' : ('vertices in (x > 0.99)', 'facet'),
}

fields = {
'accoustic_pressure' : ('complex', 1, 'Omega', 1),
}

variables = {
'p' : ('unknown field', 'accoustic_pressure', 0),
'q' : ('test field', 'accoustic_pressure', 'p'),
}

ebcs = {
}

integrals = {
'i' : 2,
}

equations = {
'Acoustic pressure' :
"""%s * dw_laplace.i.Omega( one.one, q, p )
- %s * dw_dot.i.Omega( q, p )
- %s * dw_dot.i.Gamma_out( q, p )
= %s * dw_integrate.i.Gamma_in( q )"""
% (c*c, w*w, 1j*w*c, 1j*w*c*c*rho*v_n)
}

solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton', {
'i_max' : 1,
'eps_a' : 1e-1,
'eps_r' : 1.0,
'macheps' : 1e-16,
'lin_red' : 1e-1, # Linear system error < (eps_a * lin_red).
(continues on next page)

96 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

(continued from previous page)

'ls_red' : 0.1,
'ls_red_warp' : 0.001,
'ls_on' : 1.1,
'ls_min' : 1e-5,
'check' : 0,
'delta' : 1e-6,
})
}

acoustics/acoustics3d.py

Description
Acoustic pressure distribution in 3D.
Two Laplace equations, one in Ω1 , other in Ω2 , connected on the interface region Γ12 using traces of variables.
Find two complex acoustic pressures 𝑝1 , 𝑝2 such that:
∫︁ ∫︁
𝑘 2 𝑞𝑝 − ∇𝑞 · ∇𝑝
Ω
∫︁ ∫︁ ∫︁ Ω
−𝑖𝑤/𝑐 𝑞𝑝 + 𝑖𝑤𝜌/𝑍 𝑞(𝑝2 − 𝑝1 ) + 𝑖𝑤𝜌/𝑍 𝑞(𝑝1 − 𝑝2 )
Γ𝑜𝑢𝑡 Γ2 Γ
∫︁ 1
= 𝑖𝑤𝜌 𝑣𝑛 𝑞 , ∀𝑞 .
Γ𝑖𝑛

1.5. Examples 97
SfePy Documentation, Release version: 2024.2

98 Chapter 1. Documentation
 SfePy Documentation, Release version: 2024.2

source code

r"""
Acoustic pressure distribution in 3D.

Two Laplace equations, one in :math:`\Omega_1`, other in

:math:`\Omega_2`, connected on the interface region :math:`\Gamma_{12}`
using traces of variables.

Find two complex acoustic pressures :math:`p_1`, :math:`p_2` such that:

.. math::
\int_{\Omega} k^2 q p - \int_{\Omega} \nabla q \cdot \nabla p \\
- i w/c \int_{\Gamma_{out}} q p
+ i w \rho/Z \int_{\Gamma_2} q (p_2 - p_1)
+ i w \rho/Z \int_{\Gamma_1} q (p_1 - p_2) \\
= i w \rho \int_{\Gamma_{in}} v_n q
\;, \quad \forall q \;.
"""

from __future__ import absolute_import

from sfepy import data_dir

(continues on next page)

1.5. Examples 99
SfePy Documentation, Release version: 2024.2

(continued from previous page)

filename_mesh = data_dir + '/meshes/3d/acoustics_mesh3d.mesh'

freq = 1200
v_n = 1.0 # m/s
c = 343.0 # m/s
rho = 1.55 # kg/m^3
R = 1000
w = 2.0 * freq

k1 = w / c
rhoc1 = rho * c

coef_k = ((1.0 + 0.1472 * (freq / R)**(-0.577))

+ 1j * (-0.1734 * (freq / R)**(-0.595)))
coef_r = ((1.0 + 0.0855 * (freq / R)**(-0.754))
+ 1j * (-0.0765 * (freq / R)**(-0.732)))

k2 = k1 * coef_k
rhoc2 = rhoc1 * coef_r

# perforation geometry parameters

tw = 0.9e-3
dh = 2.49e-3
por = 0.08

# acoustic impedance
Z = rho * c / por * (0.006 + 1j * k1 * (tw + 0.375 * dh
* (1 + rhoc2/rhoc1 * k2/k1)))

regions = {
'Omega' : 'all',
'Omega_1' : 'cells of group 1',
'Omega_2' : 'cells of group 2',
'Gamma_12' : ('r.Omega_1 *v r.Omega_2', 'facet'),
'Gamma_12_1' : ('copy r.Gamma_12', 'facet', 'Omega_1'),
'Gamma_12_2' : ('copy r.Gamma_12', 'facet', 'Omega_2'),
'Gamma_in' : ('vertices in (z < 0.001)', 'facet'),
'Gamma_out' : ('vertices in (z > 0.157)', 'facet'),
}

materials = {
}

fields = {
'accoustic_pressure_1' : ('complex', 'scalar', 'Omega_1', 1),
'accoustic_pressure_2' : ('complex', 'scalar', 'Omega_2', 1),
}

variables = {
'p_1' : ('unknown field', 'accoustic_pressure_1'),
'q_1' : ('test field', 'accoustic_pressure_1', 'p_1'),
'p_2' : ('unknown field', 'accoustic_pressure_2'),
(continues on next page)

100 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

'q_2' : ('test field', 'accoustic_pressure_2', 'p_2'),
}

ebcs = {
}

integrals = {
'i' : 2,
}

equations = {
'Acoustic pressure' :
"""%s * dw_dot.i.Omega_1(q_1, p_1)
+ %s * dw_dot.i.Omega_2(q_2, p_2)
- dw_laplace.i.Omega_1(q_1, p_1)
- dw_laplace.i.Omega_2(q_2, p_2)
- %s * dw_dot.i.Gamma_out(q_1, p_1)
+ %s * dw_jump.i.Gamma_12_1(q_1, p_1, tr(p_2))
+ %s * dw_jump.i.Gamma_12_2(q_2, p_2, tr(p_1))
= %s * dw_integrate.i.Gamma_in(q_1)"""
% (k1*k1, k2*k2,
1j*k1,
1j*k1*rhoc1 / Z, 1j*k2*rhoc2 / Z,
1j*k1*rhoc1 * v_n)
}

options = {
'nls': 'newton',
'ls': 'ls',
'split_results_by': 'region',
'output_dir': 'output',
}

solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton', {
'i_max' : 1,
'eps_a' : 1e-10,
'eps_r' : 1.0,
'macheps' : 1e-16,
'lin_red' : 1e-1,
'ls_red' : 0.1,
'ls_red_warp' : 0.001,
'ls_on' : 1.1,
'ls_min' : 1e-5,
'check' : 0,
'delta' : 1e-6,
})
}

1.5. Examples 101

SfePy Documentation, Release version: 2024.2

acoustics/helmholtz_apartment.py

Description
A script demonstrating the solution of the scalar Helmholtz equation for a situation inspired by the physical problem
of WiFi propagation in an apartment.
The example is an adaptation of the project found here: https://bthierry.pages.math.cnrs.fr/course-fem/projet/
2017-2018/
The boundary conditions are defined as perfect radiation conditions, meaning that on the boundary waves will radiate
outside.
The PDE for this physical process implies to find 𝐸(𝑥, 𝑡) for 𝑥 ∈ Ω such that:
{︃
∆𝐸 + 𝑘 2 𝑛(𝑥)2 𝐸 = 𝑓𝑠 ∀𝑥 ∈ Ω
𝜕𝑛 𝐸(𝑥) − 𝑖𝑘𝑛(𝑥)𝐸(𝑥) = 0 ∀𝑥 on 𝜕Ω

Usage

The mesh of the appartement and the different material ids can be visualized with the following:

sfepy-view meshes/2d/helmholtz_apartment.vtk -e -2

The example is run from the top level directory as:

sfepy-run sfepy/examples/acoustics/helmholtz_apartment.py

The result of the computation can be visualized as follows:

sfepy-view helmholtz_apartment.vtk --color-map=seismic -f imag.E:wimag.E:f10%:p0 --

˓→camera-position="-5.14968,-7.27948,7.08783,-0.463265,-0.23358,-0.350532,0.160127,0.

˓→664287,0.730124"

102 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

source code

r"""
A script demonstrating the solution of the scalar Helmholtz equation for a
situation inspired by the physical problem of WiFi propagation in an apartment.

The example is an adaptation of the project found here:

https://bthierry.pages.math.cnrs.fr/course-fem/projet/2017-2018/

The boundary conditions are defined as perfect radiation conditions, meaning

that on the boundary waves will radiate outside.

The PDE for this physical process implies to find :math:`E(x, t)` for :math:`x
\in \Omega` such that:

.. math::
\left\lbrace
\begin{aligned}
\Delta E+ k^2 n(x)^2 E = f_s && \forall x \in \Omega \\
\partial_n E(x)-ikn(x)E(x)=0 && \forall x \text{ on } \partial \Omega
\end{aligned}
\right.

(continues on next page)

1.5. Examples 103

SfePy Documentation, Release version: 2024.2

(continued from previous page)

Usage
-----

The mesh of the appartement and the different material ids can be visualized
with the following::

sfepy-view meshes/2d/helmholtz_apartment.vtk -e -2

The example is run from the top level directory as::

sfepy-run sfepy/examples/acoustics/helmholtz_apartment.py

The result of the computation can be visualized as follows::

sfepy-view helmholtz_apartment.vtk --color-map=seismic -f imag.E:wimag.E:f10%:p0 --

˓→ camera-position="-5.14968,-7.27948,7.08783,-0.463265,-0.23358,-0.350532,0.160127,0.
˓→664287,0.730124"

"""
import numpy as nm
from sfepy import data_dir

f = 2.4 * 1e9 # change this to 2.4 or 5 if you want to simulate frequencies

# typically used by your Wifi
c0 = 3e8 # m/s
k = 2 * nm.pi * f / c0

filename_mesh = data_dir + "/meshes/2d/helmholtz_apartment.vtk"

regions = {
'Omega': 'all',
'Walls': ('cells of group 1'),
'Air': ('cells of group 2'),
'Source': ('cells of group 3'),
'Gamma': ('vertices of surface', 'facet'),
}

# air and walls have different material parameters, hence the 1. and 2.4 factors
materials = {
'material': ({'kn_square': {
'Air': (k * 1.) ** 2,
'Walls': (k * 2.4) ** 2}, },),
'boundary': ({'kn': 1j * k * 2.4},),
'source_func': 'source_func',
}

def source_func(ts, coors, mode=None, **kwargs):

"""The source function for the antenna."""
epsilon = 0.1
x_center = nm.array([-3.23, -1.58])
if mode == 'qp':
(continues on next page)

104 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

dists = nm.abs(nm.sum(nm.square(coors - x_center), axis=1))
val = 3 / nm.pi / epsilon ** 2 * (1 - dists / epsilon)
return {'val': val[:, nm.newaxis, nm.newaxis]}

functions = {
'source_func': (source_func,),
}

fields = {
'electric_field': ('complex', 1, 'Omega', 1),
}

variables = {
'E': ('unknown field', 'electric_field', 1),
'q': ('test field', 'electric_field', 'E'),
}

ebcs = {
}

integrals = {
'i': 2,
}

equations = {
'Helmholtz equation':
"""- dw_laplace.i.Omega(q, E)
+ dw_dot.i.Omega(material.kn_square, q, E)
+ dw_dot.i.Gamma(boundary.kn, q, E)
= dw_integrate.i.Source(source_func.val, q)
"""
}

solvers = {
'ls': ('ls.auto_direct', {}),
'newton': ('nls.newton', {
'i_max': 1,
}),
}

options = {
'refinement_level': 3,
}

1.5. Examples 105

SfePy Documentation, Release version: 2024.2

acoustics/vibro_acoustic3d.py

Description
Vibro-acoustic problem
3D acoustic domain with 2D perforated deforming interface.
Problem definition - find 𝑝 (acoustic pressure), 𝑔 (transversal acoustic velocity), 𝑤 (plate deflection) and 𝜃 (rotation)
such that:
∫︁ ∫︁ ∫︁ ∫︁ ∫︁ ∫︁
2 2 2 + −
𝑐 ∇𝑞 · ∇𝑝 − 𝜔 𝑞𝑝 + 𝑖𝜔𝑐 𝑞𝑝 + 𝑖𝜔𝑐 𝑞𝑝 − 𝑖𝜔𝑐 (𝑞 − 𝑞 )𝑔 = 2𝑖𝜔𝑐 𝑞 𝑝¯ , ∀𝑞 ,
Ω Ω Γ𝑖𝑛
∫︁ Γ𝑜𝑢𝑡 Γ0
∫︁ ∫︁ Γ𝑖𝑛

−𝑖𝜔 𝑓 (𝑝+ − 𝑝− ) − 𝜔 2 𝐹 𝑓 𝑔 + 𝜔2 𝐶𝑓 𝑤 = 0 , ∀𝑓 ,
Γ0 Γ0 Γ0
∫︁ ∫︁ ∫︁ ∫︁
𝜔2 𝐶𝑧𝑔 − 𝜔 2 𝑆𝑧𝑤 + ∇𝑧 · 𝐺 · ∇𝑤 − 𝜃 · 𝐺 · ∇𝑧 = 0 , ∀𝑧 ,
Γ0 Γ0 Γ0 Γ0
∫︁ ∫︁ ∫︁ ∫︁
−𝜔 2 𝑅𝜈 ·𝜃 + 𝐷𝑖𝑗𝑘𝑙 𝑒𝑖𝑗 (𝜈)𝑒𝑘𝑙 (𝜃) − 𝜈 · 𝐺 · ∇𝑤 + 𝜈 · 𝐺 · 𝜃 = 0 , ∀𝜈 ,
Γ0 Γ0 Γ0 Γ0

106 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

1.5. Examples 107

SfePy Documentation, Release version: 2024.2

source code

r"""
Vibro-acoustic problem

3D acoustic domain with 2D perforated deforming interface.

Problem definition - find :math:`p` (acoustic pressure),

:math:`g` (transversal acoustic velocity),
:math:`w` (plate deflection) and :math:`\ul{\theta}` (rotation) such that:

.. math::
c^2 \int_{\Omega} \nabla q \cdot \nabla p
- \omega^2 \int_{\Omega} q p
+ i \omega c \int_{\Gamma_{in}} q p
+ i \omega c \int_{\Gamma_{out}} q p
- i \omega c^2 \int_{\Gamma_0} (q^+ - q^-) g
= 2i \omega c \int_{\Gamma_{in}} q \bar{p}
\;, \quad \forall q \;,

- i \omega \int_{\Gamma_0} f (p^+ - p^-)

- \omega^2 \int_{\Gamma_0} F f g
+ \omega^2 \int_{\Gamma_0} C f w
(continues on next page)

108 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

= 0
\;, \quad \forall f \;,

\omega^2 \int_{\Gamma_0} C z g
- \omega^2 \int_{\Gamma_0} S z w
+ \int_{\Gamma_0} \nabla z \cdot \ull{G} \cdot \nabla w
- \int_{\Gamma_0} \ul{\theta} \cdot \ull{G} \cdot \nabla z
= 0
\;, \quad \forall z \;,

- \omega^2 \int_{\Gamma_0} R\, \ul{\nu} \cdot \ul{\theta}

+ \int_{\Gamma_0} D_{ijkl} e_{ij}(\ul{\nu}) e_{kl}(\ul{\theta})
- \int_{\Gamma_0} \ul{\nu} \cdot \ull{G} \cdot \nabla w
+ \int_{\Gamma_0} \ul{\nu} \cdot \ull{G} \cdot \ul{\theta}
= 0
\;, \quad \forall \ul{\nu} \;,
"""
import numpy as nm
from sfepy import data_dir
from sfepy.mechanics.matcoefs import stiffness_from_lame

def define(sound_speed=343.0,
wave_num=5.5,
p_inc=300,
thickness=0.01,
filename_mesh='/meshes/3d/acoustic_wg.vtk'):

filename_mesh = data_dir + filename_mesh

c = sound_speed
c2 = c**2
w = wave_num * c
w2 = w**2
wc = w * c
wc2 = w * c2

regions = {
'Omega1': 'cells of group 1',
'Omega2': 'cells of group 2',
'GammaIn': ('vertices of group 1', 'face'),
'GammaOut': ('vertices of group 2', 'face'),
'Gamma_aux': ('r.Omega1 *v r.Omega2', 'face'),
'Gamma0_1': ('copy r.Gamma_aux', 'face', 'Omega1'),
'Gamma0_2': ('copy r.Gamma_aux', 'face', 'Omega2'),
'Gamma0': ('copy r.Gamma_aux', 'cell', None, {'mesh_dim': 2}),
'Left_': ('vertices in (x < 0.001)', 'edge'),
'Right_': ('vertices in (x > 0.299)', 'edge'),
'Gamma0_Left': ('r.Gamma_aux *v r.Left_', 'edge'),
'Gamma0_Right': ('r.Gamma_aux *v r.Right_', 'edge'),
}

(continues on next page)

1.5. Examples 109

SfePy Documentation, Release version: 2024.2

(continued from previous page)

fields = {
'pressure1': ('complex', 'scalar', 'Omega1', 1),
'pressure2': ('complex', 'scalar', 'Omega2', 1),
'tvelocity': ('complex', 'scalar', 'Gamma0', 1),
'deflection': ('complex', 'scalar', 'Gamma0', 1),
'rotation': ('complex', 'vector', 'Gamma0', 1),
}

variables = {
'p1': ('unknown field', 'pressure1', 0),
'q1': ('test field', 'pressure1', 'p1'),
'p2': ('unknown field', 'pressure2', 1),
'q2': ('test field', 'pressure2', 'p2'),
'g0': ('unknown field', 'tvelocity', 2),
'f0': ('test field', 'tvelocity', 'g0'),
'w': ('unknown field', 'deflection', 3),
'z': ('test field', 'deflection', 'w'),
'theta': ('unknown field', 'rotation', 4),
'nu': ('test field', 'rotation', 'theta'),
}

ebcs = {
'fixed_l': ('Gamma0_Left', {'w.0': 0.0, 'theta.all': 0.0}),
'fixed_r': ('Gamma0_Right', {'w.0': 0.0, 'theta.all': 0.0}),
}

options = {
'split_results_by': 'region',
'output_dir': 'output',
}

functions = {
}

materials = {
'ac': ({
'F': -2.064e+00,
'c': -1.064e+00,
'T': 9.202e-01,
'hG': thickness * 4.5e10 * nm.eye(2),
'hR': thickness * 0.71,
'h3R': thickness**3 / 3.0 * 0.71,
'h3C': thickness**3 / 3.0 * stiffness_from_lame(2, 1e1, 1e0)}, ),
}

equations = {
'eq_p1': """
%e * dw_laplace.5.Omega1(q1, p1)
- %e * dw_dot.5.Omega1(q1, p1)
+ %s * dw_dot.5.GammaIn(q1, p1)
- %s * dw_dot.5.Gamma0_1(q1, tr(Gamma0, g0))
= %s * dw_integrate.5.GammaIn(q1)""" % (c2, w2, 1j * wc,
(continues on next page)

110 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

1j * wc2, 2j * wc * p_inc),
'eq_p2': """
+ %e * dw_laplace.5.Omega2(q2, p2)
- %e * dw_dot.5.Omega2(q2, p2)
+ %s * dw_dot.5.GammaOut(q2, p2)
+ %s * dw_dot.5.Gamma0_2(q2, tr(Gamma0, g0))
= 0""" % (c2, w2, 1j * wc, 1j * wc2),
'eq_g0': """
- %s * dw_dot.5.Gamma0(f0, tr(Gamma0_1, p1))
+ %s * dw_dot.5.Gamma0(f0, tr(Gamma0_2, p2))
- %e * dw_dot.5.Gamma0(ac.F, f0, g0)
+ %e * dw_dot.5.Gamma0(ac.c, f0, w)
= 0""" % (1j * w, 1j * w, w2, w2),
'eq_w': """
%e * dw_dot.5.Gamma0(ac.c, z, g0)
- %e * dw_dot.5.Gamma0(ac.T, z, w)
- %e * dw_dot.5.Gamma0(ac.hR, z, w)
+ dw_diffusion.5.Gamma0(ac.hG, z, w)
- dw_v_dot_grad_s.5.Gamma0(ac.hG, theta, z)
= 0""" % (w2, w2, w2),
'eq_theta': """
- %e * dw_dot.5.Gamma0(ac.h3R, nu, theta)
+ dw_lin_elastic.5.Gamma0(ac.h3C, nu, theta)
- dw_v_dot_grad_s.5.Gamma0(ac.hG, nu, w)
+ dw_dot.5.Gamma0(ac.hG, nu, theta)
= 0""" % (w2, ),
}

solvers = {
'ls': ('ls.auto_direct', {}),
'nls': ('nls.newton', {
'i_max': 1,
'eps_a': 1e-4,
'eps_r': 1e-6,
})
}

return locals()

dg

dg/advection_1D.py

Description
Transient advection equation in 1D solved using discontinous galerkin method.

𝑑𝑝
+ 𝑎 · 𝑑𝑝/𝑑𝑥 = 0
𝑑𝑡
𝑝(𝑡, 0) = 𝑝(𝑡, 1)

1.5. Examples 111

SfePy Documentation, Release version: 2024.2

Usage Examples

Run:

sfepy-run sfepy/examples/dg/advection_1D.py

To view animated results use sfepy/examples/dg/dg_plot_1D.py specifing name of the output in output/ folder,
default is dg/advection_1D:

python3 sfepy/examples/dg/dg_plot_1D.py dg/advection_1D

dg_plot_1D.py also accepts full and relative paths:

python3 sfepy/examples/dg/dg_plot_1D.py output/dg/advection_1D

source code

r"""
Transient advection equation in 1D solved using discontinous galerkin method.

.. math:: \frac{dp}{dt} + a \cdot dp/dx = 0

p(t,0) = p(t,1)

(continues on next page)

112 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

Usage Examples
--------------
Run::

sfepy-run sfepy/examples/dg/advection_1D.py

To view animated results use ``sfepy/examples/dg/dg_plot_1D.py`` specifing name

of the output in ``output/`` folder, default is ``dg/advection_1D``::

python3 sfepy/examples/dg/dg_plot_1D.py dg/advection_1D

``dg_plot_1D.py`` also accepts full and relative paths::

python3 sfepy/examples/dg/dg_plot_1D.py output/dg/advection_1D

"""
from sfepy.examples.dg.example_dg_common import *
from sfepy.discrete.dg.limiters import MomentLimiter1D

dim = 1

def define(filename_mesh=None,
approx_order=2,

adflux=0.0,
limit=True,

cw=None,
diffcoef=None,
diffscheme="symmetric",

cfl=0.4,
dt=None,
t1=0.1
):

t0 = 0
transient = True

mstart = 0
mend = 1

diffcoef = None
cw = None

example_name = "advection_1D"
dim = 1

if filename_mesh is None:
filename_mesh = get_gen_1D_mesh_hook(0, 1, 100)

materials = {
(continues on next page)

1.5. Examples 113

SfePy Documentation, Release version: 2024.2

(continued from previous page)

'a': ({'val': [1.0], '.flux': adflux},),

regions = {
'Omega': 'all',
'Gamma': ('vertices of surface', 'facet'),
'left': ('vertices in x == 0', 'vertex'),
'right': ('vertices in x == 1', 'vertex')
}

fields = {
'f': ('real', 'scalar', 'Omega', str(approx_order) + 'd', 'DG', 'legendre')
}

variables = {
'p': ('unknown field', 'f', 0, 1),
'v': ('test field', 'f', 'p'),
}

dgebcs = {
'u_left': ('left', {'p.all': 0}),
'u_righ': ('right', {'p.all': 0}),
}

dgepbc_1 = {
'name' : 'u_rl',
'region': ['right', 'left'],
'dofs': {'p.all': 'p.all'},
'match': 'match_y_line',
}

integrals = {
'i': 2 * approx_order,
}

equations = {
'Advection': """
dw_dot.i.Omega(v, p)
- dw_s_dot_mgrad_s.i.Omega(a.val, p[-1], v)
+ dw_dg_advect_laxfrie_flux.i.Omega(a.flux, a.val, v, p[-1]) = 0
"""
}

solvers = {
"tss": ('ts.tvd_runge_kutta_3',
{"t0" : t0,
"t1" : t1,
'limiters': {"f": MomentLimiter1D} if limit else {}
}),
'nls': ('nls.newton', {}),
'ls' : ('ls.scipy_direct', {})
(continues on next page)

114 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

}

options = {
'ts' : 'tss',
'nls' : 'newton',
'ls' : 'ls',
'save_times' : 100,
'active_only' : False,
'pre_process_hook': get_cfl_setup(cfl)
if dt is None else
get_cfl_setup(dt=dt),
'output_dir' : 'output/dg/' + example_name,
'output_format' : "vtk",
}

functions = {}

def local_register_function(fun):
try:
functions.update({fun.__name__: (fun,)})

except AttributeError: # Already a sfepy Function.

fun = fun.function
functions.update({fun.__name__: (fun,)})

return fun

def four_step_p(x):
"""
piecewise constant (-inf, 1.8],(1.8, a + 4](a+4, a + 5](a + 5, inf)
"""
return nm.piecewise(x,
[x <= mstart,
x <= mstart + .4,
mstart + .4 < x,
mstart + .5 <= x],
[0, 0, .5, 0])

@local_register_function
def get_ic(x, ic=None):
return four_step_p(x)

def analytic_sol(coors, t=None, uset=False):

x = coors[..., 0]
if uset:
res = get_ic(x[..., None] - t[None, ...])
return res # for animating transient problem

res = get_ic(x[..., None])

return res[..., 0]

(continues on next page)

1.5. Examples 115

SfePy Documentation, Release version: 2024.2

(continued from previous page)

@local_register_function
def sol_fun(ts, coors, mode="qp", **kwargs):
t = ts.time
if mode == "qp":
return {"p": analytic_sol(coors, t)[..., None, None]}

ics = {
'ic': ('Omega', {'p.0': 'get_ic'}),
}

return locals()

dg/advection_2D.py

Description
Transient advection equation in 2D solved by discontinous Galerkin method.
𝑑𝑝
+ 𝑎 · 𝑔𝑟𝑎𝑑 𝑝 = 0
𝑑𝑡

Usage Examples

Run:

sfepy-run sfepy/examples/dg/advection_2D.py

Results are saved to output/dg/advection_2D folder by default as .msh files, the best way to view them is through GMSH
(http://gmsh.info/) version 4.6 or newer. Start GMSH and use File | Open menu or Crtl + O shortcut, navigate to the
output folder, select all .msh files and hit Open, all files should load as one item in Post-processing named p_cell_nodes.
GMSH is capable of rendering high order approximations in individual elements, to modify fidelity of rendering, double
click the displayed mesh, quick options menu should pop up, click on All view options.... This brings up the
Options window with View [0] selected in left column. Under the tab General ensure that Adapt visualization
grid is ticked, then you can adjust Maximum recursion depth and `Target visualization error to tune the
visualization. To see visualization elements (as opposed to mesh elements) go to Visibility tab and tick Draw
element outlines, this option is also available from quick options menu as View element outlines or under
shortcut Alt+E. In the quick options menu, you can also modify normal raise by clicking View Normal Raise to see
solution rendered as surface above the mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference elements in GMSH and Sfepy and might get
patched in the future.
source code

r"""
Transient advection equation in 2D solved by discontinous Galerkin method.

.. math:: \frac{dp}{dt} + a\cdot grad\,p = 0

Usage Examples
--------------

(continues on next page)

116 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

Run::

sfepy-run sfepy/examples/dg/advection_2D.py

Results are saved to output/dg/advection_2D folder by default as ``.msh`` files,

the best way to view them is through GMSH (http://gmsh.info/) version 4.6 or
newer. Start GMSH and use ``File | Open`` menu or Crtl + O shortcut, navigate to
the output folder, select all ``.msh`` files and hit Open, all files should load
as one item in Post-processing named p_cell_nodes.

GMSH is capable of rendering high order approximations in individual elements,

to modify fidelity of rendering, double click the displayed mesh, quick options
menu should pop up, click on ``All view options...``. This brings up the Options
window with ``View [0]`` selected in left column. Under the tab ``General``
ensure that ``Adapt visualization grid`` is ticked, then you can adjust
``Maximum recursion depth`` and ```Target visualization error`` to tune
the visualization. To see visualization elements (as opposed to mesh elements)
go to ``Visibility`` tab and tick ``Draw element outlines``, this option is also
available from quick options menu as ``View element outlines`` or under shortcut
``Alt+E``. In the quick options menu, you can also modify normal raise by
clicking ``View Normal Raise`` to see solution rendered as surface above the
mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference
elements in GMSH and Sfepy and might get patched in the future.
"""
from sfepy.examples.dg.example_dg_common import *
from sfepy.discrete.dg.limiters import MomentLimiter2D

mesh_center = (0.5, 0.25)

mesh_size = (1.0, 0.5)

def define(filename_mesh=None,
approx_order=2,

adflux=0,
limit=True,

cw=None,
diffcoef=None,
diffscheme="symmetric",

cfl=0.4,
dt=None,
t1=0.01

):

example_name = "advection_2D"
dim = 2

diffcoef = None
cw = None
(continues on next page)

1.5. Examples 117

SfePy Documentation, Release version: 2024.2

(continued from previous page)

if filename_mesh is None:
filename_mesh = get_gen_block_mesh_hook((1., 1.), (20, 20), (.5, .5))

t0 = 0.

angle = 0
# get_common(approx_order, cfl, t0, t1, None, get_ic)
rotm = nm.array([[nm.cos(angle), -nm.sin(angle)],
[nm.sin(angle), nm.cos(angle)]])
velo = nm.sum(rotm.T * nm.array([1., 0.]), axis=-1)[:, None]
materials = {
'a': ({'val': [velo], '.flux': adflux},),
}

regions = {
'Omega' : 'all',
'left': ('vertices in x == 0', 'edge'),
'right': ('vertices in x == 1', 'edge'),
'top': ('vertices in y == 1', 'edge'),
'bottom': ('vertices in y == 0', 'edge')
}

fields = {
'f': ('real', 'scalar', 'Omega', str(approx_order) + 'd', 'DG', 'legendre') #
}

variables = {
'p': ('unknown field', 'f', 0, 1),
'v': ('test field', 'f', 'p'),
}

def gsmooth(x):
r"""
.. :math: C_0^{\inf}
"""
return .3 * nm.piecewise(x, [x <= 0.1, x >= 0.1, .3 < x],
[0, lambda x:
nm.exp(1 / ((10 * (x - .2)) ** 2 - 1) + 1),
0])

def analytic_sol(coors, t):

x_1 = coors[..., 0]
x_2 = coors[..., 1]
sin = nm.sin
pi = nm.pi
exp = nm.exp
# res = four_step_u(x_1) * four_step_u(x_2)
res = gsmooth(x_1) * gsmooth(x_2)
return res

@local_register_function
(continues on next page)

118 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

def sol_fun(ts, coors, mode="qp", **kwargs):
t = ts.time
if mode == "qp":
return {"p": analytic_sol(coors, t)[..., None, None]}

def get_ic(x, ic=None):

return gsmooth(x[..., 0:1]) * gsmooth(x[..., 1:])

functions = {
'get_ic': (get_ic,)
}

ics = {
'ic': ('Omega', {'p.0': 'get_ic'}),
}

dgepbc_1 = {
'name': 'u_rl',
'region': ['right', 'left'],
'dofs': {'p.all': 'p.all'},
'match': 'match_y_line',
}

integrals = {
'i': 3 * approx_order,
}

equations = {
'Advection': """
dw_dot.i.Omega(v, p)
- dw_s_dot_mgrad_s.i.Omega(a.val, p[-1], v)
+ dw_dg_advect_laxfrie_flux.i.Omega(a.flux, a.val, v, p[-1]) = 0
"""
}

solvers = {
"tss": ('ts.tvd_runge_kutta_3',
{"t0" : t0,
"t1" : t1,
'limiters': {"f": MomentLimiter2D} if limit else {}}),
'nls': ('nls.newton',{}),
'ls' : ('ls.scipy_direct', {})
}

options = {
'ts' : 'tss',
'nls' : 'newton',
'ls' : 'ls',
'save_times' : 100,
'active_only' : False,
'output_dir' : 'output/dg/' + example_name,
(continues on next page)

1.5. Examples 119

SfePy Documentation, Release version: 2024.2

(continued from previous page)

'output_format' : 'msh',
'file_format' : 'gmsh-dg',
'pre_process_hook': get_cfl_setup(cfl) if dt is None else get_cfl_setup(dt=dt)
}

return locals()

dg/advection_diffusion_2D.py

Description
Static advection-diffusion equation in 2D solved by discontinous Galerkin method.

𝑎 · 𝑔𝑟𝑎𝑑 𝑝 − 𝑑𝑖𝑣(𝑔𝑟𝑎𝑑 𝑝) = 0

Based on
Antonietti, P., & Quarteroni, A. (2013). Numerical performance of discontinuous
and stabilized continuous Galerkin methods for convection-diffusion problems.

Usage Examples

Run:

sfepy-run sfepy/examples/dg/advection_diffusion_2D.py

Results are saved to output/dg/advection_diffusion_2D folder by default as ` .msh` files, the best way to view them is
through GMSH (http://gmsh.info/) version 4.6 or newer. Start GMSH and use File | Open menu or Crtl + O shortcut,
navigate to the output folder, select all .msh files and hit Open, all files should load as one item in Post-processing named
p_cell_nodes.
GMSH is capable of rendering high order approximations in individual elements, to modify fidelity of rendering, double
click the displayed mesh, quick options menu should pop up, click on All view options.... This brings up the
Options window with View [0] selected in left column. Under the tab General ensure that Adapt visualization
grid is ticked, then you can adjust Maximum recursion depth and `Target visualization error to tune the
visualization. To see visualization elements (as opposed to mesh elements) go to Visibility tab and tick Draw
element outlines, this option is also available from quick options menu as View element outlines or under
shortcut Alt+E. In the quick options menu, you can also modify normal raise by clicking View Normal Raise to see
solution rendered as surface above the mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference elements in GMSH and Sfepy and might get
patched in the future.
source code

r"""
Static advection-diffusion equation in 2D solved by discontinous Galerkin method.

.. math:: a \cdot grad\, p - div(grad\,p) = 0

Based on

Antonietti, P., & Quarteroni, A. (2013). Numerical performance of discontinuous

(continues on next page)

120 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

and stabilized continuous Galerkin methods for convection-diffusion problems.

Usage Examples
--------------

Run::

sfepy-run sfepy/examples/dg/advection_diffusion_2D.py

Results are saved to output/dg/advection_diffusion_2D folder by default as `

`.msh`` files, the best way to view them is through GMSH (http://gmsh.info/)
version 4.6 or newer. Start GMSH and use ``File | Open`` menu or Crtl + O
shortcut, navigate to the output folder, select all ``.msh`` files and hit Open,
all files should load as one item in Post-processing named p_cell_nodes.

GMSH is capable of rendering high order approximations in individual elements,

to modify fidelity of rendering, double click the displayed mesh, quick options
menu should pop up, click on ``All view options...``. This brings up the Options
window with ``View [0]`` selected in left column. Under the tab ``General``
ensure that ``Adapt visualization grid`` is ticked, then you can adjust
``Maximum recursion depth`` and ```Target visualization error`` to tune
the visualization. To see visualization elements (as opposed to mesh elements)
go to ``Visibility`` tab and tick ``Draw element outlines``, this option is also
available from quick options menu as ``View element outlines`` or under shortcut
``Alt+E``. In the quick options menu, you can also modify normal raise by
clicking ``View Normal Raise`` to see solution rendered as surface above the
mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference
elements in GMSH and Sfepy and might get patched in the future.
"""

from sfepy.examples.dg.example_dg_common import *

def define(filename_mesh=None,
approx_order=3,

adflux=0,
limit=False,

cw=1000,
diffcoef=1,
diffscheme="symmetric",

cfl=None,
dt=None,
):

cfl = None
dt = None

(continues on next page)

1.5. Examples 121

SfePy Documentation, Release version: 2024.2

(continued from previous page)

functions = {}
def local_register_function(fun):
try:
functions.update({fun.__name__: (fun,)})

except AttributeError: # Already a sfepy Function.

fun = fun.function
functions.update({fun.__name__: (fun,)})

return fun

example_name = "advection_diffusion_2D"
dim = 2

if filename_mesh is None:
filename_mesh = get_gen_block_mesh_hook((1., 1.), (20, 20), (.5, .5))

velo = [1., 1.]

angle = 0.0 # - nm.pi / 5

rotm = nm.array([[nm.cos(angle), -nm.sin(angle)],
[nm.sin(angle), nm.cos(angle)]])
velo = nm.sum(rotm.T * nm.array(velo), axis=-1)[:, None]

regions = {
'Omega' : 'all',
'left' : ('vertices in x == 0', 'edge'),
'right': ('vertices in x == 1', 'edge'),
'top' : ('vertices in y == 1', 'edge'),
'bottom': ('vertices in y == 0', 'edge')
}

fields = {
'f': ('real', 'scalar', 'Omega', str(approx_order) + 'd', 'DG', 'legendre')
}

variables = {
'p': ('unknown field', 'f', 0),
'v': ('test field', 'f', 'p'),
}

integrals = {
'i': 2 * approx_order,
}

@local_register_function
def bc_funs(ts, coors, bc, problem):
# return 2*coors[..., 1]
t = ts.dt*ts.step
x_1 = coors[..., 0]
x_2 = coors[..., 1]
(continues on next page)

122 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

res = nm.zeros(nm.shape(x_1))

sin = nm.sin
cos = nm.cos
exp = nm.exp
pi = nm.pi

if bc.diff == 0:
if "left" in bc.name:
res[:] = 0
elif "right" in bc.name:
res[:] = 0
elif "bottom" in bc.name:
res[:] = 0 #-2*sin(2*pi*x_1)
elif "top" in bc.name:
res[:] = 0

elif bc.diff == 1:
if "left" in bc.name:
res = nm.stack((-2*pi*(x_2**2 - x_2),
res),
axis=-2)
elif "right" in bc.name:
res = nm.stack((-2*pi*(x_2**2 - x_2), res,),
axis=-2)
elif "bot" in bc.name:
res = nm.stack((res,
sin(2*pi*x_1)),
axis=-2)
elif "top" in bc.name:
res = nm.stack((res,
-sin(2*pi*x_1)),
axis=-2)

return res

@local_register_function
def source_fun(ts, coors, mode="qp", **kwargs):
# t = ts.dt * ts.step
eps = diffcoef
sin = nm.sin
cos = nm.cos
exp = nm.exp
sqrt = nm.sqrt
pi = nm.pi
if mode == "qp":
x_1 = coors[..., 0]
x_2 = coors[..., 1]
res = -2*pi*(x_2**2 - x_2)*cos(2*pi*x_1)\
- 2*(2*pi**2*(x_2**2 - x_2)*sin(2*pi*x_1) - sin(2*pi*x_1))*eps\
- (2*x_2 - 1)*sin(2*pi*x_1)
(continues on next page)

1.5. Examples 123

SfePy Documentation, Release version: 2024.2

(continued from previous page)

return {"val": res[..., None, None]}

def analytic_sol(coors, t):

x_1 = coors[..., 0]
x_2 = coors[..., 1]
sin = nm.sin
pi = nm.pi
res = -(x_2 ** 2 - x_2) * sin(2 * pi * x_1)
return res

@local_register_function
def sol_fun(ts, coors, mode="qp", **kwargs):
t = ts.time
if mode == "qp":
return {"p": analytic_sol(coors, t)[..., None, None]}

dgebcs = {
'u_left' : ('left', {'p.all': "bc_funs", 'grad.p.all' : "bc_funs"}),
'u_top' : ('top', {'p.all': "bc_funs", 'grad.p.all' : "bc_funs"}),
'u_bot' : ('bottom', {'p.all': "bc_funs", 'grad.p.all' : "bc_funs"}),
'u_right': ('right', {'p.all': "bc_funs", 'grad.p.all' : "bc_funs"}),
}

materials = {
'a' : ({'val': [velo], '.flux': adflux},),
'D' : ({'val': [diffcoef], '.cw': cw},),
'g' : 'source_fun'
}

equations = {
'balance': """
- dw_s_dot_mgrad_s.i.Omega(a.val, p, v)
+ dw_dg_advect_laxfrie_flux.i.Omega(a.flux, a.val, v, p)
"""
+
" + dw_laplace.i.Omega(D.val, v, p) " +
diffusion_schemes_implicit[diffscheme] +
" + dw_dg_interior_penalty.i.Omega(D.val, D.cw, v, p)" +
" - dw_volume_lvf.i.Omega(g.val, v)" +
"= 0"
}

solver_0 = {
'name' : 'ls',
'kind' : 'ls.scipy_direct',
}

solver_1 = {
'name' : 'newton',
(continues on next page)

124 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

'kind' : 'nls.newton',

'i_max' : 5,
'eps_a' : 1e-8,
'eps_r' : 1.0,
'macheps' : 1e-16,
'lin_red' : 1e-2, # Linear system error < (eps_a * lin_red).
'ls_red' : 0.1,
'ls_red_warp' : 0.001,
'ls_on' : 0.99999,
'ls_min' : 1e-5,
'check' : 0,
'delta' : 1e-6,
}

options = {
'nls' : 'newton',
'ls' : 'ls',
'output_dir' : 'output/dg/' + example_name,
'output_format' : 'msh',
'file_format' : 'gmsh-dg'
}
return locals()

dg/burgers_2D.py

Description
Burgers equation in 2D solved using discontinous Galerkin method

𝑑𝑝
+ 𝑑𝑖𝑣 𝑓 (𝑝) − 𝑑𝑖𝑣(𝑔𝑟𝑎𝑑 𝑝) = 0
𝑑𝑡
Based on
Kučera, V. (n.d.). Higher order methods for the solution of compressible flows. Charles University. p. 21 eq. (1.39)

Usage Examples

Run:

sfepy-run sfepy/examples/dg/burgers_2D.py

Results are saved to output/dg/burgers_2D folder by default as .msh files, the best way to view them is through GMSH
(http://gmsh.info/) version 4.6 or newer. Start GMSH and use File | Open menu or Crtl + O shortcut, navigate to the
output folder, select all .msh files and hit Open, all files should load as one item in Post-processing named p_cell_nodes.
GMSH is capable of rendering high order approximations in individual elements, to modify fidelity of rendering, double
click the displayed mesh, quick options menu should pop up, click on All view options.... This brings up the
Options window with View [0] selected in left column. Under the tab General ensure that Adapt visualization
grid is ticked, then you can adjust Maximum recursion depth and `Target visualization error to tune the
visualization. To see visualization elements (as opposed to mesh elements) go to Visibility tab and tick Draw
element outlines, this option is also available from quick options menu as View element outlines or under

1.5. Examples 125

SfePy Documentation, Release version: 2024.2

shortcut Alt+E. In the quick options menu, you can also modify normal raise by clicking View Normal Raise to see
solution rendered as surface above the mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference elements in GMSH and Sfepy and might get
patched in the future.
source code

r"""
Burgers equation in 2D solved using discontinous Galerkin method

.. math:: \frac{dp}{dt} + div\,f(p) - div(grad\,p) = 0

Based on

Kučera, V. (n.d.). Higher order methods for the solution of compressible flows.
Charles University. p. 21 eq. (1.39)

Usage Examples
--------------

Run::

sfepy-run sfepy/examples/dg/burgers_2D.py

Results are saved to output/dg/burgers_2D folder by default as ``.msh`` files,

the best way to view them is through GMSH (http://gmsh.info/) version 4.6 or
newer. Start GMSH and use ``File | Open`` menu or Crtl + O shortcut, navigate to
the output folder, select all ``.msh`` files and hit Open, all files should load
as one item in Post-processing named p_cell_nodes.

GMSH is capable of rendering high order approximations in individual elements,

to modify fidelity of rendering, double click the displayed mesh, quick options
menu should pop up, click on ``All view options...``. This brings up the Options
window with ``View [0]`` selected in left column. Under the tab ``General``
ensure that ``Adapt visualization grid`` is ticked, then you can adjust
``Maximum recursion depth`` and ```Target visualization error`` to tune
the visualization. To see visualization elements (as opposed to mesh elements)
go to ``Visibility`` tab and tick ``Draw element outlines``, this option is also
available from quick options menu as ``View element outlines`` or under shortcut
``Alt+E``. In the quick options menu, you can also modify normal raise by
clicking ``View Normal Raise`` to see solution rendered as surface above the
mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference
elements in GMSH and Sfepy and might get patched in the future.
"""

from sfepy.examples.dg.example_dg_common import *

from sfepy import data_dir

from sfepy.discrete.dg.limiters import MomentLimiter2D, IdentityLimiter

mesh_center = (0, 0)
(continues on next page)

126 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

mesh_size = (2, 2)

def define(filename_mesh=None,
approx_order=2,

adflux=0,
limit=False,

cw=10,
diffcoef=0.002,
diffscheme="symmetric",

cfl=None,
dt=1e-5,
t1=0.01
):

functions = {}
def local_register_function(fun):
try:
functions.update({fun.__name__: (fun,)})

except AttributeError: # Already a sfepy Function.

fun = fun.function
functions.update({fun.__name__: (fun,)})

return fun

example_name = "burgers_2D"
dim = 2

if filename_mesh is None:
filename_mesh = data_dir + "/meshes/2d/square_tri2.mesh"

t0 = 0.
if dt is None and cfl is None:
dt = 1e-5

velo = [1., 1.]

angle = 0 # - nm.pi / 5
rotm = nm.array([[nm.cos(angle), -nm.sin(angle)],
[nm.sin(angle), nm.cos(angle)]])
velo = nm.sum(rotm.T * nm.array(velo), axis=-1)[:, None]
burg_velo = velo.T / nm.linalg.norm(velo)

regions = {
'Omega': 'all',
'left' : ('vertices in x == -1', 'edge'),
'right': ('vertices in x == 1', 'edge'),
'top' : ('vertices in y == 1', 'edge'),
(continues on next page)

1.5. Examples 127

SfePy Documentation, Release version: 2024.2

(continued from previous page)

'bottom': ('vertices in y == -1', 'edge')
}

fields = {
'f': ('real', 'scalar', 'Omega',
str(approx_order) + 'd', 'DG', 'legendre')
}

variables = {
'p': ('unknown field', 'f', 0, 1),
'v': ('test field', 'f', 'p'),
}

integrals = {
'i': 5,
}

def analytic_sol(coors, t):

x_1 = coors[..., 0]
x_2 = coors[..., 1]
sin = nm.sin
pi = nm.pi
exp = nm.exp
res = -(exp(-t) - 1)*(sin(5*x_1*x_2) + sin(-4*x_1*x_2 + 4*x_1 + 4*x_2))
return res

@local_register_function
def sol_fun(ts, coors, mode="qp", **kwargs):
t = ts.time
if mode == "qp":
return {"p": analytic_sol(coors, t)[..., None, None]}

@local_register_function
def bc_funs(ts, coors, bc, problem):
# return 2*coors[..., 1]
t = ts.dt*ts.step
x_1 = coors[..., 0]
x_2 = coors[..., 1]
sin = nm.sin
cos = nm.cos
exp = nm.exp
if bc.diff == 0:
if "left" in bc.name:
res = -(exp(-t) - 1)*(sin(-5*x_2) + sin(8*x_2 - 4))
elif "bottom" in bc.name:
res = -(exp(-t) - 1) * (sin(-5 * x_1) + sin(8 * x_1 - 4))
elif "right" in bc.name:
res = -(exp(-t) - 1)*(sin(4) + sin(5*x_2))
elif "top" in bc.name:
res = -(exp(-t) - 1)*(sin(4) + sin(5*x_1))

elif bc.diff == 1:
(continues on next page)

128 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

if "left" in bc.name:
res = nm.stack(((4*(x_2 - 1)*cos(4) - 5*x_2*cos(5*x_2))*
(exp(-t) - 1),
-5*(exp(-t) - 1)*cos(5*x_2)),
axis=-2)
elif "bottom" in bc.name:
res = nm.stack(((5*cos(-5*x_1) - 8*cos(8*x_1 - 4))*(exp(-t) - 1),
-(5*x_1*cos(-5*x_1) - 4*(x_1 - 1)*cos(8*x_1 - 4))*
(exp(-t) - 1)),
axis=-2)

elif "right" in bc.name:

res = nm.stack(((4*(x_2 - 1)*cos(4) - 5*x_2*cos(5*x_2))*
(exp(-t) - 1),
-5*(exp(-t) - 1)*cos(5*x_2)),
axis=-2)
elif "top" in bc.name:
res = nm.stack((-5*(exp(-t) - 1)*cos(5*x_1),
(4*(x_1 - 1)*cos(4) - 5*x_1*cos(5*x_1))*
(exp(-t) - 1)),
axis=-2)

return res

@local_register_function
def source_fun(ts, coors, mode="qp", **kwargs):
if mode == "qp":
t = ts.dt * ts.step
x_1 = coors[..., 0]
x_2 = coors[..., 1]
sin = nm.sin
cos = nm.cos
exp = nm.exp
res = (
+ (5 * x_1 * cos(5 * x_1 * x_2)
- 4 * (x_1 - 1) * cos(4 * x_1 * x_2 - 4 * x_1 - 4 * x_2)) *
(exp(-t) - 1) ** 2 * (sin(5 * x_1 * x_2)
- sin(4 * x_1 * x_2 - 4 * x_1 - 4 * x_2))
+ (5 * x_2 * cos(5 * x_1 * x_2)
- 4 * (x_2 - 1) * cos(4 * x_1 * x_2 - 4 * x_1 - 4 * x_2)) *
(exp(-t) - 1) ** 2 * (sin(5 * x_1 * x_2)
- sin(4 * x_1 * x_2 - 4 * x_1 - 4 * x_2))
- diffcoef *
((25 * x_1 ** 2 * sin(5 * x_1 * x_2) - 16 * (x_1 - 1) ** 2 *
sin(4 * x_1 * x_2 - 4 * x_1 - 4 * x_2)) * (exp(-t) - 1)
+ (25 * x_2 ** 2 * sin(5 * x_1 * x_2) - 16 * (x_2 - 1) ** 2 *
sin(4 * x_1 * x_2 - 4 * x_1 - 4 * x_2)) * (exp(-t) - 1))
+ (sin(5 * x_1 * x_2) - sin(4 * x_1 * x_2 - 4 * x_1 - 4 * x_2))*
exp(-t)
)
return {"val": res[..., None, None]}

(continues on next page)

1.5. Examples 129

SfePy Documentation, Release version: 2024.2

(continued from previous page)

def adv_fun(p):
vu = velo.T * p[..., None]
return vu

def adv_fun_d(p):
v1 = velo.T * nm.ones(p.shape + (1,))
return v1

def burg_fun(p):
vu = .5*burg_velo * p[..., None] ** 2
return vu

def burg_fun_d(p):
v1 = burg_velo * p[..., None]
return v1

materials = {
'a' : ({'val': [velo], '.flux':adflux},),
'D' : ({'val': [diffcoef], '.Cw': cw},),
'g' : 'source_fun'
}

ics = {
'ic': ('Omega', {'p.0': 0}),
}

dgebcs = {
'u_left' : ('left', {'p.all': 'bc_funs', 'grad.p.all': 'bc_funs'}),
'u_right' : ('right', {'p.all': 'bc_funs', 'grad.p.all': 'bc_funs'}),
'u_bottom' : ('bottom', {'p.all': 'bc_funs', 'grad.p.all': 'bc_funs'}),
'u_top' : ('top', {'p.all': 'bc_funs', 'grad.p.all': 'bc_funs'}),
}

equations = {
'balance':
"dw_dot.i.Omega(v, p)" +
# non-linear hyperbolic terms
" - dw_ns_dot_grad_s.i.Omega(burg_fun, burg_fun_d, p[-1], v)" +
" + dw_dg_nonlinear_laxfrie_flux.i.Omega(a.flux, burg_fun, burg_fun_d, v, p[-1])
˓→" +

# diffusion
" + dw_laplace.i.Omega(D.val, v, p[-1])" +
diffusion_schemes_explicit[diffscheme] +
" - dw_dg_interior_penalty.i.Omega(D.val, D.Cw, v, p[-1])"
# source
+ " - dw_volume_lvf.i.Omega(g.val, v)"
" = 0"
}

solvers = {
"tss.tvd_runge_kutta_3": ('ts.tvd_runge_kutta_3',
(continues on next page)

130 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

{"t0": t0,
"t1": t1,
'limiters': {
"f": MomentLimiter2D} if limit else {}}),
"tss.euler": ('ts.euler',
{"t0" : t0,
"t1" : t1,
'limiters': {"f": MomentLimiter2D} if limit else {}}),
'nls': ('nls.newton', {}),
'ls' : ('ls.scipy_direct', {})
}

options = {
'ts' : 'tss.euler',
'nls' : 'nls.newton',
'ls' : 'ls.mumps',
'save_times' : 100,
'output_dir' : 'output/dg/' + example_name,
'output_format' : 'msh',
'file_format' : 'gmsh-dg',
'pre_process_hook': get_cfl_setup(CFL=cfl, dt=dt)
}

return locals()

dg/dg_plot_1D.py

Description
Script for plotting 1D DG FEM data stored in VTK files
source code
#!/usr/bin/env python
# -*- coding:utf-8 -*-
"""
Script for plotting 1D DG FEM data stored in VTK files
"""

import glob

import matplotlib
matplotlib.use('Qt5Agg')
import matplotlib.pyplot as plt
import argparse

from os.path import join as pjoin

import os
import sys
from glob import glob

from sfepy.discrete.dg.dg_1D_vizualizer import \

(continues on next page)

1.5. Examples 131

SfePy Documentation, Release version: 2024.2

(continued from previous page)

(load_1D_vtks, animate_1D_DG_sol, load_state_1D_vtk, plot1D_legendre_dofs,
reconstruct_legendre_dofs)

def load_and_plot_fun(folder, filename, t0, t1, tn,

ic_fun=None, exact=None,
compare=False, polar=False):
"""
Parameters
----------
folder : str
folder where to look for files
filename : str
used in {name}.i.vtk, i = 0,1, ... tns - 1
t0 : float
starting time
t1 : int
final time
tn : int
number of time steps
ic_fun : callable
initital condition
exact : callable
exact solution, for transient problems function of space ant time
compare : bool
polar : bool
"""

# load time data

lmesh, u = load_1D_vtks(folder, filename)
animate_1D_DG_sol(lmesh, t0, t1, u, tn=tn, ic=ic_fun, exact=exact,
delay=100, polar=polar)

if compare:
files = glob(pjoin(folder, "*.vtk"))

n_first = sorted(files)[0].split(".")[-2]
n_last = sorted(files)[-1].split(".")[-2]

print("Plotting files {} and {} to compare first and last".format(n_first, n_

˓→ last))

coors, u_end = load_state_1D_vtk(pjoin(folder, filename + "." + n_last + ".vtk"))

coors, u_start = load_state_1D_vtk(pjoin(folder, filename + "." + n_first + ".vtk
˓→ "))

plot1D_legendre_dofs(coors, [u_start.swapaxes(0, 1)[:, :, 0], u_end.swapaxes(0,␣

˓→ 1)[:, :, 0]])

plt.figure("Reconstructed solution")
ww_s, xx = reconstruct_legendre_dofs(coors, None, u_end)
ww_e, _ = reconstruct_legendre_dofs(coors, None, u_start)
(continues on next page)

132 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

plt.plot(xx, ww_s[:, 0], label="IC")
plt.plot(xx, ww_e[:, 0], label="Last")
plt.legend()
plt.show()

def main(argv):
parser = argparse.ArgumentParser(description='Plotting of 1D DG data in VTK files',
epilog='(c) 2019 T. Zitka , Man-machine'
+' Interaction at NTC UWB')
parser.add_argument("input_name", help="Folder or name of the example in "
+ "output folder with VTK data, file names "
+ "in format <name>.[0-9]*.vtk. If not "
+ "provided asks for the name of the example"
, nargs="?")

parser.add_argument("-s", "--search", help="Search for different orders in "

"provided folder",
action="store_true", default=False)
parser.add_argument("-t0", "--start_time", type=float, default=0,
help="Start time of the simulation")
parser.add_argument("-t1", "--end_time", type=float, default=1.,
help="End time of the simulation")
parser.add_argument("-o", "--order", type=int, default=None,
help="Order of the approximation, when example folder"
+ " cantains more orders this chooses the one")
parser.add_argument("-cf", "--compare-final", action="store_true",
help="To compare starting and final time - " +
"usefull for periodic boundary problems")
parser.add_argument("-p", "--polar", help="Plot in polar projection",
action="store_true")

if argv is None:
argv = sys.argv[1:]
args = parser.parse_args(argv)

t0 = args.start_time
t1 = args.end_time
cf = args.compare_final
pol = args.polar
if args.input_name is None:
input_name = str(input("Please provide name of the example in output/" +
" folder or full path to data: "))
else:
input_name = args.input_name

if os.path.isdir(input_name):
full_infolder_path = os.path.abspath(input_name)
else:
input_name = pjoin("output", input_name)
if os.path.isdir(input_name):
full_infolder_path = os.path.abspath(input_name)
(continues on next page)

1.5. Examples 133

SfePy Documentation, Release version: 2024.2

(continued from previous page)

else:
print("Example {} not found in {}".format(input_name,
os.path.abspath("output")))
return

print("Input folder found: {}".format(full_infolder_path))

if args.search:
print("Input folder contains results for orders {}"
.format([os.path.basename(fol)
for fol in
glob(pjoin(full_infolder_path, "[0-9]*"))]))

if args.order is None:
order = str(input("Please provide order of approximation, default is 1: "))
if len(order) == 0:
order = 1
else:
try:
order = int(order)
except ValueError:
print("Value {} for order not understood!".format(order))
return
else:
order = args.order
print("Looking for results of order {}".format(order))
full_infolder_path = pjoin(full_infolder_path, str(order))
if not os.path.isdir(full_infolder_path):
print("Input folder with order {} not found".format(full_infolder_path))
return
else:
order = args.order

contents = glob(pjoin(full_infolder_path, "*.vtk"))

print()
tn = len(contents) # we assume the contents are time step data files
if tn == 0:
print("Input folder {} is empty!".format(full_infolder_path))
return
base_name = os.path.basename(contents[0]).split(".")[0]
print("Found {} files, basename is {}".format(tn, base_name))
print("Plotting ...")
load_and_plot_fun(full_infolder_path, base_name, t0, t1, tn,
compare=cf, polar=pol)

if __name__ == '__main__':
main(sys.argv[1:])

134 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

dg/example_dg_common.py

Description
Functions common to DG examples
source code

"""
Functions common to DG examples
"""
import os
from glob import glob

import numpy as nm

from sfepy.base.base import output

from sfepy.discrete.fem import Mesh
from sfepy.discrete.fem.meshio import UserMeshIO
from sfepy.mesh.mesh_generators import gen_block_mesh

diffusion_schemes_implicit = {
"symmetric":
" + dw_dg_diffusion_flux.i.Omega(D.val, p, v)"
+ " + dw_dg_diffusion_flux.i.Omega(D.val, v, p)",
"non-symmetric":
" + dw_dg_diffusion_flux.i.Omega(D.val, p, v)"
+ " - dw_dg_diffusion_flux.i.Omega(D.val, v, p)",
"incomplete":
" + dw_dg_diffusion_flux.i.Omega(D.val, p, v)"}

diffusion_schemes_explicit = {
"symmetric":
" - dw_dg_diffusion_flux.i.Omega(D.val, p[-1], v)"
+ " - dw_dg_diffusion_flux.i.Omega(D.val, v, p[-1])",
"non-symmetric":
" - dw_dg_diffusion_flux.i.Omega(D.val, p[-1], v)"
+ " + dw_dg_diffusion_flux.i.Omega(D.val, v, p[-1])",
"incomplete":
" - dw_dg_diffusion_flux.i.Omega(D.val, p[-1], v)"}

functions = {}
def local_register_function(fun):
try:
functions.update({fun.__name__: (fun,)})

except AttributeError: # Already a sfepy Function.

fun = fun.function
functions.update({fun.__name__: (fun,)})

return fun

(continues on next page)

1.5. Examples 135

SfePy Documentation, Release version: 2024.2

(continued from previous page)

def get_cfl_setup(CFL=None, dt=None):
"""
Provide either CFL or dt to create preprocess hook that sets up
Courant-Friedrichs-Levi stability condition for either advection or
diffusion.

Parameters
----------
CFL : float, optional
dt: float, optional

Returns
-------
setup_cfl_condition : callable
expects sfepy.discrete.problem as argument

"""

if CFL is None and dt is None:

raise ValueError("Specifiy either CFL or dt in CFL setup")

def setup_cfl_condition(problem):
"""
Sets up CFL condition for problem ts_conf in problem

Parameters
----------
problem : discrete.problem.Problem
"""
ts_conf = problem.ts_conf
mesh = problem.domain.mesh
dim = mesh.dim
first_field = list(problem.fields.values())[0]
first_field_name = list(problem.fields.keys())[0]
approx_order = first_field.approx_order
mats = problem.create_materials(['a', 'D'])
try:
# make this more general?
# maybe require material name in parameter
velo = problem.conf_materials['material_a__0'].values["val"]
max_velo = nm.max(nm.linalg.norm(velo))
except KeyError:
max_velo = 1

try:
# make this more general?
# maybe require material name in parameter
diffusion = problem.conf_materials['material_D__0'].values["val"]
max_diffusion = nm.max(nm.linalg.norm(diffusion))
except KeyError:
max_diffusion = None

(continues on next page)

136 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

dx = nm.min(problem.domain.mesh.cmesh.get_volumes(dim))

output("Preprocess hook - setup_cfl_condition:...")

output("Approximation order of field {}({}) is {}"
.format(first_field_name, first_field.family_name, approx_order))
output("Space divided into {0} cells, {1} steps, step size {2}"
.format(mesh.n_el, len(mesh.coors), dx))

if dt is None:
adv_dt = get_cfl_advection(max_velo, dx, approx_order, CFL)
diff_dt = get_cfl_diffusion(max_diffusion, dx, approx_order, CFL)
_dt = min(adv_dt, diff_dt)
else:
output("CFL coefficient {0} ignored, dt specified directly"
.format(CFL))
_dt = dt

tn = int(nm.ceil((ts_conf.t1 - ts_conf.t0) / _dt))

dtdx = _dt / dx

ts_conf.dt = _dt
ts_conf.n_step = tn
ts_conf.cour = max_velo * dtdx

output("Time divided into {0} nodes, {1} steps, step size is {2}"
.format(tn - 1, tn, _dt))
output("Courant number c = max(norm(a)) * dt/dx = {0}"
.format(ts_conf.cour))
output("Time stepping solver is {}".format(ts_conf.kind))
output("... CFL setup done.")

return setup_cfl_condition

def get_cfl_advection(max_velo, dx, approx_order, CFL):

"""

Parameters
----------
max_velo : float
dx : float
approx_order : int
CFL : CFL

Returns
-------
dt : float
"""
order_corr = 1. / (2 * approx_order + 1)

dt = dx / max_velo * CFL * order_corr

(continues on next page)

1.5. Examples 137

SfePy Documentation, Release version: 2024.2

(continued from previous page)

if not (nm.isfinite(dt)):
dt = 1
output(("CFL advection: CFL coefficient was {0} " +
"and order correction 1/{1} = {2}")
.format(CFL, (2 * approx_order + 1), order_corr))
output("CFL advection: resulting dt={}".format((dt)))
return dt

def get_cfl_diffusion(max_diffusion, dx, approx_order, CFL,

do_order_corr=False):
"""

Parameters
----------
max_diffusion : float
dx : float
approx_order : int
CFL : float
do_order_corr : bool

Returns
-------
dt : float
"""
if max_diffusion is None:
return 1

if do_order_corr:
order_corr = 1. / (2 * approx_order + 1)
else:
order_corr = 1

dt = dx**2 / max_diffusion * CFL * order_corr

if not (nm.isfinite(dt)):
dt = 1
output(("CFL diffusion: CFL coefficient was {0} " +
"and order correction 1/{1} = {2}")
.format(CFL, (2 * approx_order + 1), order_corr))
output("CFL diffusion: resulting dt={}".format(dt))
return dt

def get_gen_1D_mesh_hook(XS, XE, n_nod):

"""

Parameters
----------
XS : float
(continues on next page)

138 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

leftmost coordinate
XE : float
rightmost coordinate
n_nod : int
number of nodes, number of cells is then n_nod - 1

Returns
-------
mio : UserMeshIO instance
"""
def mesh_hook(mesh, mode):
"""
Generate the 1D mesh.
"""
if mode == 'read':

coors = nm.linspace(XS, XE, n_nod).reshape((n_nod, 1))

conn = nm.arange(n_nod, dtype=nm.int32) \
.repeat(2)[1:-1].reshape((-1, 2))
mat_ids = nm.zeros(n_nod - 1, dtype=nm.int32)
descs = ['1_2']

mesh = Mesh.from_data('uniform_1D{}'.format(n_nod), coors, None,

[conn], [mat_ids], descs)
return mesh

elif mode == 'write':

pass

mio = UserMeshIO(mesh_hook)
return mio

def get_gen_block_mesh_hook(dims, shape, centre, mat_id=0, name='block',

coors=None, verbose=True):
"""

Parameters
----------
dims : array of 2 or 3 floats
Dimensions of the block.
shape : array of 2 or 3 ints
Shape (counts of nodes in x, y, z) of the block mesh.
centre : array of 2 or 3 floats
Centre of the block.
mat_id : int, optional
The material id of all elements.
name : string
Mesh name.
verbose : bool
If True, show progress of the mesh generation.

Returns
(continues on next page)

1.5. Examples 139

SfePy Documentation, Release version: 2024.2

(continued from previous page)

-------
mio : UserMeshIO instance
"""
def mesh_hook(mesh, mode):
"""
Generate the 1D mesh.
"""
if mode == 'read':

mesh = gen_block_mesh(dims, shape, centre, mat_id=mat_id, name=name,

coors=coors, verbose=verbose)
return mesh

elif mode == 'write':

pass

mio = UserMeshIO(mesh_hook)
return mio

def clear_folder(clear_format, confirm=False, doit=True):

"""
Deletes files matching the format

Parameters
----------
clear_format : str
confirm : bool
doit : bool
if False do not delete anything no matter the confirmation

Returns
-------
deleted_anything :
True if there was something to delete
"""
files = glob(clear_format)
if confirm:
for file in files:
output("Will delete file {}".format(file))
doit = input("--------------\nDelete files [Y/n]? ").strip() == "Y"

if doit:
for file in files:
os.remove(file)
return bool(files)

140 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

dg/imperative_burgers_1D.py

Description
Burgers equation in 1D solved using discontinous Galerkin method
source code

#!/usr/bin/env python
"""
Burgers equation in 1D solved using discontinous Galerkin method
"""
import argparse
import sys
sys.path.append('.')
from os.path import join as pjoin

import numpy as nm

from sfepy.examples.dg.example_dg_common import \

clear_folder, get_gen_1D_mesh_hook
from sfepy.examples.dg.dg_plot_1D import load_and_plot_fun

# sfepy imports
from sfepy.base.base import IndexedStruct
from sfepy.base.base import Struct, configure_output, output
from sfepy.discrete import (FieldVariable, Material, Integral, Function,
Equation, Equations, Problem)
from sfepy.discrete.conditions import InitialCondition, EssentialBC, Conditions
from sfepy.discrete.dg.fields import DGField
from sfepy.discrete.dg.limiters import MomentLimiter1D
from sfepy.discrete.fem import FEDomain
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.solvers.ts_dg_solvers import TVDRK3StepSolver
from sfepy.terms.terms_dg import Term

def parse_args(argv=None):
if argv is None:
argv = sys.argv

parser = argparse.ArgumentParser(
description='Solve Burgers equation and display animated results, '
'change script code to modify the problem.',
epilog='(c) 2019 T. Zitka , Man-machine Interaction at NTC UWB')
parser.add_argument('-o', '--output-dir', default='.',
help='output directory')
parser.add_argument('-p', '--plot',
action='store_true', dest='plot',
default=False, help='plot animated results')
options = parser.parse_args(argv[1:])
return options

def main(argv=None):
(continues on next page)

1.5. Examples 141

SfePy Documentation, Release version: 2024.2

(continued from previous page)

options = parse_args(argv=argv)

# vvvvvvvvvvvvvvvv #
approx_order = 2
# ^^^^^^^^^^^^^^^^ #

# Setup output names

outputs_folder = options.output_dir

domain_name = "domain_1D"
problem_name = "iburgers_1D"
output_folder = pjoin(outputs_folder, problem_name, str(approx_order))
output_format = "vtk"
save_timestn = 100
clear_folder(pjoin(output_folder, "*." + output_format))
configure_output({'output_screen': True,
'output_log_name':
pjoin(output_folder,
f"last_run_{problem_name}_{approx_order}.txt")})

# ------------
# | Get mesh |
# ------------
X1 = 0.
XN = 1.
n_nod = 100
n_el = n_nod - 1
mesh = get_gen_1D_mesh_hook(X1, XN, n_nod).read(None)

# -----------------------------
# | Create problem components |
# -----------------------------

integral = Integral('i', order=approx_order * 2)

domain = FEDomain(domain_name, mesh)
omega = domain.create_region('Omega', 'all')
left = domain.create_region('Gamma1',
'vertices in x == %.10f ' % X1,
'vertex')
right = domain.create_region('Gamma2',
'vertices in x == %.10f ' % XN,
'vertex')
field = DGField('dgfu', nm.float64, 'scalar', omega,
approx_order=approx_order)

u = FieldVariable('u', 'unknown', field, history=1)

v = FieldVariable('v', 'test', field, primary_var_name='u')

MassT = Term.new('dw_dot(v, u)', integral, omega, u=u, v=v)

velo = nm.array(1.0)

(continues on next page)

142 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

def adv_fun(u):
vu = velo.T * u[..., None]
return vu

def adv_fun_d(u):
v1 = velo.T * nm.ones(u.shape + (1,))
return v1

burg_velo = velo.T / nm.linalg.norm(velo)

def burg_fun(u):
vu = burg_velo * u[..., None] ** 2
return vu

def burg_fun_d(u):
v1 = 2 * burg_velo * u[..., None]
return v1

StiffT = Term.new('dw_ns_dot_grad_s(fun, fun_d, u[-1], v)',

integral, omega,
u=u, v=v,
fun=burg_fun, fun_d=burg_fun_d)

# alpha = Material('alpha', val=[.0])

# FluxT = AdvectDGFluxTerm("adv_lf_flux(a.val, v, u)", "a.val, v, u[-1]",
# integral, omega, u=u, v=v, a=a, alpha=alpha)

FluxT = Term.new('dw_dg_nonlinear_laxfrie_flux(fun, fun_d, v, u[-1])',

integral, omega,
u=u, v=v,
fun=burg_fun, fun_d=burg_fun_d)

eq = Equation('balance', MassT - StiffT + FluxT)

eqs = Equations([eq])

# ------------------------------
# | Create boundary conditions |
# ------------------------------
left_fix_u = EssentialBC('left_fix_u', left, {'u.all': 1.0})
right_fix_u = EssentialBC('right_fix_u', right, {'u.all': 0.0})

# ----------------------------
# | Create initial condition |
# ----------------------------
def ghump(x):
(continues on next page)

1.5. Examples 143

SfePy Documentation, Release version: 2024.2

(continued from previous page)

"""
Nice gaussian.
"""
return nm.exp(-200 * x ** 2)

def ic_wrap(x, ic=None):

return ghump(x - .3)

ic_fun = Function('ic_fun', ic_wrap)

ics = InitialCondition('ic', omega, {'u.0': ic_fun})

# ------------------
# | Create problem |
# ------------------
pb = Problem(problem_name,
equations=eqs,
conf=Struct(options={"save_times": save_timestn},
ics={}, ebcs={}, epbcs={}, lcbcs={},
materials={}),
active_only=False)
pb.setup_output(output_dir=output_folder, output_format=output_format)
pb.set_ics(Conditions([ics]))

# ------------------
# | Create limiter |
# ------------------
limiter = MomentLimiter1D

# ---------------------------
# | Set time discretization |
# ---------------------------
CFL = .2
max_velo = nm.max(nm.abs(velo))
t0 = 0
t1 = .2
dx = nm.min(mesh.cmesh.get_volumes(1))
dt = dx / max_velo * CFL / (2 * approx_order + 1)
tn = int(nm.ceil((t1 - t0) / dt))
dtdx = dt / dx

# ------------------
# | Create solver |
# ------------------
ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({'is_linear': True}, lin_solver=ls, status=nls_status)

tss_conf = {'t0' : t0,

't1' : t1,
'n_step' : tn,
'limiters': {"dgfu": limiter}}
(continues on next page)

144 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

tss = TVDRK3StepSolver(tss_conf,
nls=nls, context=pb, verbose=True)

# ---------
# | Solve |
# ---------
pb.set_solver(tss)
state_end = pb.solve()

output("Solved equation \n\n\t\t u_t - div(f(u))) = 0\n")

output(f"With IC: {ic_fun.name}")
# output("and EBCs: {}".format(pb.ebcs.names))
# output("and EPBCS: {}".format(pb.epbcs.names))
output("-------------------------------------")
output(f"Approximation order is {approx_order}")
output(f"Space divided into {mesh.n_el} cells, " +
f"{len(mesh.coors)} steps, step size is {dx}")
output(f"Time divided into {tn - 1} nodes, {tn} steps, step size is {dt}")
output(f"CFL coefficient was {CFL} and " +
f"order correction {1 / (2 * approx_order + 1)}")
output(f"Courant number c = max(abs(u)) * dt/dx = {max_velo * dtdx}")
output("------------------------------------------")
output(f"Time stepping solver is {tss.name}")
output(f"Limiter used: {limiter.name}")
output("======================================")

# ----------
# | Plot 1D|
# ----------
if options.plot:
load_and_plot_fun(output_folder, domain_name,
t0, t1, min(tn, save_timestn),
ic_fun)

if __name__ == '__main__':
main()

dg/laplace_2D.py

Description
Laplace equation solved in 2d by discontinous Galerkin method

−𝑑𝑖𝑣(𝑔𝑟𝑎𝑑 𝑝) = 0

on rectangle
p = 0 p_y = 0
[0,b]—————————–[a, b]
|

1.5. Examples 145

SfePy Documentation, Release version: 2024.2

|
p_x = -a | p(x,y) | p_x = 0 p = 0 | | p = 0
|
[0,0]—————————–[a, 0]
p_y = b p = 0
solution to this is

𝑝(𝑥, 𝑦) = 1/2 * 𝑥 * *2 − 1/2 * 𝑦 * *2 − 𝑎 * 𝑥 + 𝑏 * 𝑦

Usage Examples

Run:

sfepy-run sfepy/examples/dg/laplace_2D.py

Results are saved to output/dg/laplace_2D folder by default as .msh files, the best way to view them is through GMSH
(http://gmsh.info/) version 4.6 or newer. Start GMSH and use File | Open menu or Crtl + O shortcut, navigate to the
output folder, select all .msh files and hit Open, all files should load as one item in Post-processing named p_cell_nodes.
GMSH is capable of rendering high order approximations in individual elements, to modify fidelity of rendering, double
click the displayed mesh, quick options menu should pop up, click on All view options.... This brings up the
Options window with View [0] selected in left column. Under the tab General ensure that Adapt visualization
grid is ticked, then you can adjust Maximum recursion depth and `Target visualization error to tune the
visualization. To see visualization elements (as opposed to mesh elements) go to Visibility tab and tick Draw
element outlines, this option is also available from quick options menu as View element outlines or under
shortcut Alt+E. In the quick options menu, you can also modify normal raise by clicking View Normal Raise to see
solution rendered as surface above the mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference elements in GMSH and Sfepy and might get
patched in the future.
source code

r"""
Laplace equation solved in 2d by discontinous Galerkin method

.. math:: - div(grad\,p) = 0

on rectangle
p = 0
p_y = 0
[0,b]-----------------------------[a, b]
| |
| |
p_x = -a | p(x,y) | p_x = 0
p = 0 | | p = 0
| |
[0,0]-----------------------------[a, 0]
p_y = b
p = 0

solution to this is
.. math:: p(x,y) = 1/2*x**2 - 1/2*y**2 - a*x + b*y
(continues on next page)

146 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

Usage Examples
--------------

Run::

sfepy-run sfepy/examples/dg/laplace_2D.py

Results are saved to output/dg/laplace_2D folder by default as ``.msh`` files,

the best way to view them is through GMSH (http://gmsh.info/) version 4.6 or
newer. Start GMSH and use ``File | Open`` menu or Crtl + O shortcut, navigate to
the output folder, select all ``.msh`` files and hit Open, all files should load
as one item in Post-processing named p_cell_nodes.

GMSH is capable of rendering high order approximations in individual elements,

to modify fidelity of rendering, double click the displayed mesh, quick options
menu should pop up, click on ``All view options...``. This brings up the Options
window with ``View [0]`` selected in left column. Under the tab ``General``
ensure that ``Adapt visualization grid`` is ticked, then you can adjust
``Maximum recursion depth`` and ```Target visualization error`` to tune
the visualization. To see visualization elements (as opposed to mesh elements)
go to ``Visibility`` tab and tick ``Draw element outlines``, this option is also
available from quick options menu as ``View element outlines`` or under shortcut
``Alt+E``. In the quick options menu, you can also modify normal raise by
clicking ``View Normal Raise`` to see solution rendered as surface above the
mesh. Note that for triangular meshes normal raise -1 produces expected raise
above the mesh. This is due to the opposite orientation of the reference
elements in GMSH and Sfepy and might get patched in the future.
"""

from sfepy.examples.dg.example_dg_common import *

def define(filename_mesh=None,
approx_order=2,

adflux=None,
limit=False,

cw=100,
diffcoef=1,
diffscheme="symmetric",

cfl=None,
dt=None,
):

cfl = None
dt = None

functions = {}
def local_register_function(fun):
(continues on next page)

1.5. Examples 147

SfePy Documentation, Release version: 2024.2

(continued from previous page)

try:
functions.update({fun.__name__: (fun,)})

except AttributeError: # Already a sfepy Function.

fun = fun.function
functions.update({fun.__name__: (fun,)})

return fun

example_name = "laplace_2D"
dim = 2

if filename_mesh is None:
filename_mesh = get_gen_block_mesh_hook((1., 1.), (16, 16), (.5, .5))

a = 1
b = 1
c = 0

regions = {
'Omega' : 'all',
'left' : ('vertices in x == 0', 'edge'),
'right': ('vertices in x == 1', 'edge'),
'top' : ('vertices in y == 1', 'edge'),
'bottom': ('vertices in y == 0', 'edge')
}
fields = {
'f': ('real', 'scalar', 'Omega', str(approx_order) + 'd', 'DG', 'legendre') #
}

variables = {
'p': ('unknown field', 'f', 0, 1),
'v': ('test field', 'f', 'p'),
}

def analytic_sol(coors, t):

x_1, x_2 = coors[..., 0], coors[..., 1]
res = 1/2*x_1**2 - 1/2*x_2**2 - a*x_1 + b*x_2 + c
return res

@local_register_function
def sol_fun(ts, coors, mode="qp", **kwargs):
t = ts.time
if mode == "qp":
return {"p": analytic_sol(coors, t)[..., None, None]}

@local_register_function
def bc_funs(ts, coors, bc, problem):
t = ts.time
x_1, x_2 = coors[..., 0], coors[..., 1]
res = nm.zeros(x_1.shape)
(continues on next page)

148 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

if bc.diff == 0:
res[:] = analytic_sol(coors, t)

elif bc.diff == 1:
res = nm.stack((x_1 - a, -x_2 + b),
axis=-2)
return res

materials = {
'D' : ({'val': [diffcoef], '.Cw': cw},),
}

dgebcs = {
'u_left' : ('left', {'p.all': "bc_funs", 'grad.p.all': "bc_funs"}),
'u_right' : ('right', {'p.all': "bc_funs", 'grad.p.all': "bc_funs"}),
'u_bottom' : ('bottom', {'p.all': "bc_funs", 'grad.p.all': "bc_funs"}),
'u_top' : ('top', {'p.all': "bc_funs", 'grad.p.all': "bc_funs"}),

integrals = {
'i': 2 * approx_order,
}

equations = {
'laplace': " dw_laplace.i.Omega(D.val, v, p) " +
diffusion_schemes_implicit[diffscheme] +
" + dw_dg_interior_penalty.i.Omega(D.val, D.Cw, v, p)" +
" = 0"
}

solver_0 = {
'name' : 'ls',
'kind' : 'ls.scipy_direct',
}

solver_1 = {
'name' : 'newton',
'kind' : 'nls.newton',

# 'i_max' : 5,
# 'eps_a' : 1e-8,
# 'eps_r' : 1.0,
# 'macheps' : 1e-16,
# 'lin_red' : 1e-2, # Linear system error < (eps_a * lin_red).
# 'ls_red' : 0.1,
# 'ls_red_warp' : 0.001,
# 'ls_on' : 0.99999,
# 'ls_min' : 1e-5,
# 'check' : 0,
# 'delta' : 1e-6,
}
(continues on next page)

1.5. Examples 149

SfePy Documentation, Release version: 2024.2

(continued from previous page)

options = {
'nls' : 'newton',
'ls' : 'ls',
'output_dir' : 'output/dg/' + example_name,
'output_format' : 'msh',
'file_format' : 'gmsh-dg',
# 'pre_process_hook': get_cfl_setup(cfl)
}
return locals()

diffusion

diffusion/cube.py

Description
Laplace equation (e.g. temperature distribution) on a cube geometry with different boundary condition values on the
cube sides. This example was used to create the SfePy logo.
Find 𝑇 such that:
∫︁
𝑐∇𝑠 · ∇𝑇 = 0 , ∀𝑠 .
Ω

150 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

source code

r"""
Laplace equation (e.g. temperature distribution) on a cube geometry with
different boundary condition values on the cube sides. This example was
used to create the SfePy logo.

Find :math:`T` such that:

.. math::
\int_{\Omega} c \nabla s \cdot \nabla T
= 0
\;, \quad \forall s \;.
"""
from __future__ import absolute_import
from sfepy import data_dir

#filename_mesh = data_dir + '/meshes/3d/cube_big_tetra.mesh'

filename_mesh = data_dir + '/meshes/3d/cube_medium_hexa.mesh'

############# Laplace.

material_1 = {
(continues on next page)

1.5. Examples 151

SfePy Documentation, Release version: 2024.2

(continued from previous page)

'name' : 'coef',
'values' : {'val' : 1.0},
}

field_1 = {
'name' : 'temperature',
'dtype' : 'real',
'shape' : (1,),
'region' : 'Omega',
'approx_order' : 1,
}

if filename_mesh.find('cube_medium_hexa.mesh') >= 0:
region_1000 = {
'name' : 'Omega',
'select' : 'cells of group 0',
}
integral_1 = {
'name' : 'i',
'order' : 1,
}
solver_0 = {
'name' : 'ls',
'kind' : 'ls.scipy_direct',
}

elif filename_mesh.find('cube_big_tetra.mesh') >= 0:

region_1000 = {
'name' : 'Omega',
'select' : 'cells of group 6',
}
integral_1 = {
'name' : 'i',
'quadrature' : 'custom',
'vals' : [[1./3., 1./3., 1./3.]],
'weights' : [0.5]
}
solver_0 = {
'name' : 'ls',
'kind' : 'ls.scipy_iterative',

'method' : 'cg',
'i_max' : 1000,
'eps_r' : 1e-12,
}

variable_1 = {
'name' : 'T',
'kind' : 'unknown field',
'field' : 'temperature',
'order' : 0, # order in the global vector of unknowns
}
(continues on next page)

152 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

variable_2 = {
'name' : 's',
'kind' : 'test field',
'field' : 'temperature',
'dual' : 'T',
}

region_0 = {
'name' : 'Surface',
'select' : 'vertices of surface',
'kind' : 'facet',
}
region_1 = {
'name' : 'Bottom',
'select' : 'vertices in (z < -0.4999999)',
'kind' : 'facet',
}
region_2 = {
'name' : 'Top',
'select' : 'vertices in (z > 0.4999999)',
'kind' : 'facet',
}
region_03 = {
'name' : 'Left',
'select' : 'vertices in (x < -0.4999999)',
'kind' : 'facet',
}

ebc_1 = {
'name' : 'T0',
'region' : 'Surface',
'dofs' : {'T.0' : -3.0},
}
ebc_4 = {
'name' : 'T1',
'region' : 'Top',
'dofs' : {'T.0' : 1.0},
}
ebc_3 = {
'name' : 'T2',
'region' : 'Bottom',
'dofs' : {'T.0' : -1.0},
}
ebc_2 = {
'name' : 'T3',
'region' : 'Left',
'dofs' : {'T.0' : 2.0},
}

equations = {
'nice_equation' : """dw_laplace.i.Omega( coef.val, s, T ) = 0""",
(continues on next page)

1.5. Examples 153

SfePy Documentation, Release version: 2024.2

(continued from previous page)

}

solver_1 = {
'name' : 'newton',
'kind' : 'nls.newton',

'i_max' : 1,
'eps_a' : 1e-10,
'eps_r' : 1.0,
'macheps' : 1e-16,
'lin_red' : 1e-2, # Linear system error < (eps_a * lin_red).
'ls_red' : 0.1,
'ls_red_warp' : 0.001,
'ls_on' : 1.1,
'ls_min' : 1e-5,
'check' : 0,
'delta' : 1e-6,
}

diffusion/darcy_flow_multicomp.py

Description
Each of the two equations describes a flow in one compartment of a porous medium. The equations are based on the
Darcy flow and the i-th compartment is defined in Ω𝑖 .
∫︁ ∫︁ ∑︁ ∫︁
𝐾 𝑖 ∇𝑝𝑖 · ∇𝑞 𝑖 + ¯ 𝑘 𝑝𝑖 − 𝑝𝑗 𝑞 𝑖 = 𝑓 𝑖 𝑞𝑖 ,
(︀ )︀
𝐺𝛼
Ω𝑖 Ω𝑖 𝑗 Ω𝑖

∀𝑞 𝑖 ∈ 𝑄𝑖 , 𝑖, 𝑗 = 1, 2 and 𝑖 ̸= 𝑗,
where 𝐾 𝑖 is the local permeability of the i-th compartment, 𝐺𝛼 ¯ 𝑘 = 𝐺𝑖 is the perfusion coefficient related to the
𝑗
compartments 𝑖 and 𝑗, 𝑓 𝑖 are sources or sinks which represent the external flow into the i-th compartment and 𝑝𝑖 is the
pressure in the i-th compartment.

154 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

source code

r"""
Each of the two equations describes a flow in one compartment of a porous
medium. The equations are based on the Darcy flow and the i-th compartment is
defined in :math:`\Omega_{i}`.

.. math::
\int_{\Omega_{i}} K^{i} \nabla p^{i} \cdot \nabla q^{i}+\int_{\Omega_{i}}
\sum_{j} \bar{G}\alpha_{k} \left( p^{i}-p^{j} \right)q^{i}
= \int_{\Omega_{i}} f^{i} q^{i},
.. math::
\forall q^{i} \in Q^{i}, \quad i,j=1,2 \quad \mbox{and} \quad i\neq j,

where :math:`K^{i}` is the local permeability of the i-th compartment,

:math:`\bar{G}\alpha_{k} = G^{i}_{j}` is the perfusion coefficient
related to the compartments :math:`i` and :math:`j`, :math:`f^i` are
sources or sinks which represent the external flow into the i-th
compartment and :math:`p^{i}` is the pressure in the i-th compartment.
"""

from __future__ import absolute_import

from sfepy.base.base import Struct
(continues on next page)

1.5. Examples 155

SfePy Documentation, Release version: 2024.2

(continued from previous page)

import numpy as nm
from sfepy import data_dir

filename_mesh = data_dir + '/meshes/3d/cube_medium_hexa.mesh'

G_bar = 2.0
alpha1 = 1.3
alpha2 = 1.0

materials = {
'mat': ('mat_fun')
}

regions = {
'Omega': 'cells of group 0',
'Sigma_1': ('vertex 0', 'vertex'),
'Omega1': ('copy r.Omega', 'cell', 'Omega'),
'Omega2': ('copy r.Omega', 'cell', 'Omega'),
'Source': 'cell 24',
'Sink': 'cell 1',
}

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

variables = {
'p1': ('unknown field', 'pressure'),
'q1': ('test field', 'pressure', 'p1'),
'p2': ('unknown field', 'pressure'),
'q2': ('test field', 'pressure', 'p2'),
}

ebcs = {
'P1': ('Sigma_1', {'p1.0' : 0.0}),
}

equations = {
'komp1': """dw_diffusion.5.Omega1(mat.K, q1, p1)
+ dw_dot.5.Omega1(mat.G_alfa, q1, p1)
- dw_dot.5.Omega1(mat.G_alfa, q1, p2)
= dw_integrate.5.Source(mat.f_1, q1)""",

'komp2': """dw_diffusion.5.Omega2(mat.K, q2, p2)

+ dw_dot.5.Omega2(mat.G_alfa, q2, p2)
- dw_dot.5.Omega2(mat.G_alfa, q2, p1)
= dw_integrate.5.Sink(mat.f_2, q2)"""
}

solvers = {
'ls': ('ls.scipy_direct', {}),
'newton': ('nls.newton',
{'i_max' : 1,
(continues on next page)

156 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

'eps_a' : 1e-6,
'eps_r' : 1.0,
})
}

def mat_fun(ts, coors, mode=None, **kwargs):

if mode == 'qp':
nqp, dim = coors.shape
alpha = nm.zeros((nqp,1,1), dtype=nm.float64)
alpha[0:nqp // 2,...] = alpha1
alpha[nqp // 2:,...] = alpha2
K = nm.eye(dim, dtype=nm.float64)
K2 = nm.tile(K, (nqp,1,1))
out = {
'K' : K2,
'f_1': 20.0 * nm.ones((nqp,1,1), dtype=nm.float64),
'f_2': -20.0 * nm.ones((nqp,1,1), dtype=nm.float64),
'G_alfa': G_bar * alpha,
}

return out

functions = {
'mat_fun': (mat_fun,),
}

options = {
'post_process_hook': 'postproc',
}

def postproc(out, pb, state, extend=False):

alpha = pb.evaluate('ev_integrate_mat.5.Omega(mat.G_alfa, p1)',
mode='el_avg')
out['alpha'] = Struct(name='output_data',
mode='cell',
data=alpha.reshape(alpha.shape[0], 1, 1, 1),
dofs=None)
return out

diffusion/laplace_1d.py

Description
Laplace equation in 1D with a variable coefficient.
Because the mesh is trivial in 1D, it is generated by mesh_hook(), and registered using UserMeshIO.
Find 𝑡 such that:
∫︁
d𝑠 d𝑡
𝑐(𝑥) =0, ∀𝑠 ,
Ω d𝑥 d𝑥
where the coefficient 𝑐(𝑥) = 0.1 + sin(2𝜋𝑥)2 is computed in get_coef().
View the results using:

1.5. Examples 157

SfePy Documentation, Release version: 2024.2

sfepy-view laplace_1d.vtk -f t:wt 1:vw

source code

r"""
Laplace equation in 1D with a variable coefficient.

Because the mesh is trivial in 1D, it is generated by :func:`mesh_hook()`, and

registered using :class:`UserMeshIO <sfepy.discrete.fem.meshio.UserMeshIO>`.

Find :math:`t` such that:

.. math::
\int_{\Omega} c(x) \tdiff{s}{x} \tdiff{t}{x}
= 0
\;, \quad \forall s \;,

where the coefficient :math:`c(x) = 0.1 + \sin(2 \pi x)^2` is computed in

:func:`get_coef()`.

View the results using::

sfepy-view laplace_1d.vtk -f t:wt 1:vw

(continues on next page)

158 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

"""
from __future__ import absolute_import
import numpy as nm
from sfepy.discrete.fem import Mesh
from sfepy.discrete.fem.meshio import UserMeshIO

def mesh_hook(mesh, mode):

"""
Generate the 1D mesh.
"""
if mode == 'read':
n_nod = 101

coors = nm.linspace(0.0, 1.0, n_nod).reshape((n_nod, 1))

conn = nm.arange(n_nod, dtype=nm.int32).repeat(2)[1:-1].reshape((-1, 2))
mat_ids = nm.zeros(n_nod - 1, dtype=nm.int32)
descs = ['1_2']

mesh = Mesh.from_data('laplace_1d', coors, None,

[conn], [mat_ids], descs)
return mesh

elif mode == 'write':

pass

def get_coef(ts, coors, mode=None, **kwargs):

if mode == 'qp':
x = coors[:, 0]

val = 0.1 + nm.sin(2 * nm.pi * x)**2

val.shape = (coors.shape[0], 1, 1)

return {'val' : val}

filename_mesh = UserMeshIO(mesh_hook)

materials = {
'coef' : 'get_coef',
}

functions = {
'get_coef' : (get_coef,),
}

regions = {
'Omega' : 'all',
'Gamma_Left' : ('vertices in (x < 0.00001)', 'facet'),
'Gamma_Right' : ('vertices in (x > 0.99999)', 'facet'),
}

fields = {
'temperature' : ('real', 1, 'Omega', 1),
(continues on next page)

1.5. Examples 159

SfePy Documentation, Release version: 2024.2

(continued from previous page)

}

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

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

integrals = {
'i' : 2,
}

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

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

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

diffusion/laplace_coupling_lcbcs.py

Description
Two Laplace equations with multiple linear combination constraints.
The two equations are coupled by a periodic-like boundary condition constraint with a shift, given as a non-
homogeneous linear combination boundary condition.

160 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

Find 𝑢 such that:

∫︁
∇𝑣1 · ∇𝑢1 = 0 , ∀𝑣1 ,
Ω1
∫︁
∇𝑣2 · ∇𝑢2 = 0 , ∀𝑣2 ,
Ω2
𝑢1 = 0 on Γ𝑏𝑜𝑡𝑡𝑜𝑚 ,
𝑢2 = 1 on Γ𝑡𝑜𝑝 ,
𝑢1 (𝑥) = 𝑢2 (𝑥) + 𝑎(𝑥) for 𝑥 ∈ Γ = Ω̄1 ∩ Ω̄2
𝑢1 (𝑥) = 𝑢1 (𝑦) + 𝑏(𝑦) for 𝑥 ∈ Γ𝑙𝑒𝑓 𝑡 , 𝑦 ∈ Γ𝑟𝑖𝑔ℎ𝑡 , 𝑦 = 𝑃 (𝑥) ,
𝑢1 = 𝑐11 in Ω𝑚11 ⊂ Ω1 ,
𝑢1 = 𝑐12 in Ω𝑚12 ⊂ Ω1 ,
𝑢2 = 𝑐2 in Ω𝑚2 ⊂ Ω2 ,

where 𝑎(𝑥), 𝑏(𝑦) are given functions (shifts), 𝑃 is the periodic coordinate mapping and 𝑐11 , 𝑐12 and 𝑐2 are unknown
constant values - the unknown DOFs in Ω𝑚11 , Ω𝑚12 and Ω𝑚2 are replaced by the integral mean values.
View the results using:

sfepy-view square_quad.vtk -f u1:wu1:p0 1:vw:p0 u2:wu2:p1 1:vw:p1

source code

1.5. Examples 161

SfePy Documentation, Release version: 2024.2

r"""
Two Laplace equations with multiple linear combination constraints.

The two equations are coupled by a periodic-like boundary condition constraint

with a shift, given as a non-homogeneous linear combination boundary condition.

Find :math:`u` such that:

.. math::
\int_{\Omega_1} \nabla v_1 \cdot \nabla u_1
= 0
\;, \quad \forall v_1 \;,

\int_{\Omega_2} \nabla v_2 \cdot \nabla u_2

= 0
\;, \quad \forall v_2 \;,

u_1 = 0 \mbox{ on } \Gamma_{bottom} \;,

u_2 = 1 \mbox{ on } \Gamma_{top} \;,

u_1(\ul{x}) = u_2(\ul{x}) + a(\ul{x}) \mbox{ for }

\ul{x} \in \Gamma = \bar\Omega_1 \cap \bar\Omega_2

u_1(\ul{x}) = u_1(\ul{y}) + b(\ul{y}) \mbox{ for }

\ul{x} \in \Gamma_{left}, \ul{y} \in \Gamma_{right}, \ul{y} = P(\ul{x}) \;,

u_1 = c_{11} \mbox{ in } \Omega_{m11} \subset \Omega_1 \;,

u_1 = c_{12} \mbox{ in } \Omega_{m12} \subset \Omega_1 \;,

u_2 = c_2 \mbox{ in } \Omega_{m2} \subset \Omega_2 \;,

where :math:`a(\ul{x})`, :math:`b(\ul{y})` are given functions (shifts),

:math:`P` is the periodic coordinate mapping and :math:`c_{11}`, :math:`c_{12}`
and :math:`c_2` are unknown constant values - the unknown DOFs in
:math:`\Omega_{m11}`, :math:`\Omega_{m12}` and :math:`\Omega_{m2}` are replaced
by the integral mean values.

View the results using::

sfepy-view square_quad.vtk -f u1:wu1:p0 1:vw:p0 u2:wu2:p1 1:vw:p1

"""
from __future__ import absolute_import
import numpy as nm

import sfepy.discrete.fem.periodic as per

from sfepy import data_dir

filename_mesh = data_dir + '/meshes/2d/square_quad.mesh'

options = {
'nls' : 'newton',
(continues on next page)

162 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

(continued from previous page)

'ls' : 'ls',
}

def get_shift1(ts, coors, region):

val = 0.1 * coors[:, 0]

return val

def get_shift2(ts, coors, region):

val = nm.empty_like(coors[:, 1])
val.fill(0.3)

return val

functions = {
'get_shift1' : (get_shift1,),
'get_shift2' : (get_shift2,),
'match_y_line' : (per.match_y_line,),
'match_x_line' : (per.match_x_line,),
}

fields = {
'scalar1': ('real', 1, 'Omega1', 1),
'scalar2': ('real', 1, 'Omega2', 1),
}

materials = {
}

variables = {
'u1' : ('unknown field', 'scalar1', 0),
'v1' : ('test field', 'scalar1', 'u1'),
'u2' : ('unknown field', 'scalar2', 1),
'v2' : ('test field', 'scalar2', 'u2'),
}

regions = {
'Omega1' : 'cells of group 1',
'Omega2' : 'cells of group 2',
'Omega_m1' : 'r.Omega1 -v (r.Gamma +s vertices of surface)',
'Omega_m11' : 'r.Omega_m1 *v vertices in (x < 0)',
'Omega_m12' : 'r.Omega_m1 *v vertices in (x > 0)',
'Omega_m2' : 'r.Omega2 -v (r.Gamma +s vertices of surface)',
'Left' : ('vertices in (x < -0.499)', 'facet'),
'Right' : ('vertices in (x > 0.499)', 'facet'),
'Bottom' : ('vertices in ((y < -0.499))', 'facet'),
'Top' : ('vertices in ((y > 0.499))', 'facet'),
'Gamma' : ('r.Omega1 *v r.Omega2 -v vertices of surface', 'facet'),
'Gamma1' : ('copy r.Gamma', 'facet', 'Omega1'),
'Gamma2' : ('copy r.Gamma', 'facet', 'Omega2'),
}

(continues on next page)

1.5. Examples 163

SfePy Documentation, Release version: 2024.2

(continued from previous page)

ebcs = {
'fix1' : ('Top', {'u2.all' : 1.0}),
'fix2' : ('Bottom', {'u1.all' : 0.0}),
}

lcbcs = {
'shifted1' : (('Gamma1', 'Gamma2'),
{'u1.all' : 'u2.all'},
'match_x_line', 'shifted_periodic',
'get_shift1'),
'shifted2' : (('Left', 'Right'),
{'u1.all' : 'u1.all'},
'match_y_line', 'shifted_periodic',
'get_shift2'),
'mean11' : ('Omega_m11', {'u1.all' : None}, None, 'integral_mean_value'),
'mean12' : ('Omega_m12', {'u1.all' : None}, None, 'integral_mean_value'),
'mean2' : ('Omega_m2', {'u2.all' : None}, None, 'integral_mean_value'),
}

integrals = {
'i1' : 2,
}

equations = {
'eq1' : """
dw_laplace.i1.Omega1(v1, u1) = 0
""",
'eq2' : """
dw_laplace.i1.Omega2(v2, u2) = 0
""",
}

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

diffusion/laplace_fluid_2d.py

Description
A Laplace equation that models the flow of “dry water” around an obstacle shaped like a Citroen CX.

164 Chapter 1. Documentation

 SfePy Documentation, Release version: 2024.2

Description

As discussed e.g. in the Feynman lectures Section 12-5 of Volume 2 (https://www.feynmanlectures.caltech.edu/II_12.

html#Ch12-S5), the flow of an irrotational and incompressible fluid can be modelled with a potential 𝑣 = 𝑔𝑟𝑎𝑑(𝜑)
that obeys

∇ · ∇𝜑 = ∆𝜑 = 0

The weak formulation for this problem is to find 𝜑 such that:

∫︁ ∫︁ ∫︁ ∫︁ ∫︁
∇𝜓 · ∇𝜑 = 𝑣0 · 𝑛 𝜓 + 𝑣0 · 𝑛 𝜓 + 𝑣0 · 𝑛 𝜓 + 𝑣0 · 𝑛 𝜓 , ∀𝜓 ,
Ω Γ𝑙𝑒𝑓 𝑡 Γ𝑟𝑖𝑔ℎ𝑡 Γ𝑡𝑜𝑝 Γ𝑏𝑜𝑡𝑡𝑜𝑚

where 𝑣 0 is the 2D vector defining the far field velocity that generates the incompressible flow.
Since the value of the potential is defined up to a constant value, a Dirichlet boundary condition is set at a single vertex
to avoid having a singular matrix.

Usage examples

This example can be run with the following:

sfepy-run sfepy/examples/diffusion/laplace_fluid_2d.py
sfepy-view citroen.vtk -f phi:p0 phi:t50:p0 --2d-view

Generating the mesh

The mesh can be generated with:

gmsh -2 -f msh22 meshes/2d/citroen.geo -o meshes/2d/citroen.msh

sfepy-convert --2d meshes/2d/citroen.msh meshes/2d/citroen.h5

1.5. Examples 165

SfePy Documentation, Release version: 2024.2

source code

r"""
A Laplace equation that models the flow of "dry water" around an obstacle
shaped like a Citroen CX.

Description
-----------

As discussed e.g. in the Feynman lectures Section 12-5 of Volume 2

(https://www.feynmanlectures.caltech.edu/II_12.html#Ch12-S5),
the flow of an irrotational and incompressible fluid can be modelled with a
potential :math:`\ul{v} = \ul{grad}(\phi)` that obeys

.. math::
\nabla \cdot \ul{\nabla}\,\phi = \Delta\,\phi = 0

The weak formulation for this problem is to find :math:`\phi` such that:

.. math::
\int_{\Omega} \nabla \psi \cdot \nabla \phi
= \int_{\Gamma_{left}} \ul{v}_0 \cdot n \, \psi
(continues on next page)

166