The following 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 info. The
$ indicates the command prompt of your terminal.

The example problems in the examples directory can be computed by the script
simple.py which is in the top-level directory of the SfePy distribution.
If it is run without arguments, a help message is printed:

$ ./simple.py
Usage: simple.py [options] filename_in
Solve partial differential equations given in a SfePy problem definition file.
Example problem definition files can be found in ``examples/`` directory of the
SfePy top-level directory. This script works with all the examples except those
in ``examples/standalone/``.
Both normal and parametric study runs are supported. A parametric study allows
repeated runs for varying some of the simulation parameters - see
``examples/diffusion/poisson_parametric_study.py`` file.
Options:
--version show program's version number and exit
-h, --help show this help message and exit
-c "key : value, ...", --conf="key : value, ..."
override problem description file items, written as
python dictionary without surrounding braces
-O "key : value, ...", --options="key : value, ..."
override options item of problem description, written
as python dictionary without surrounding braces
-d "key : value, ...", --define="key : value, ..."
pass given arguments written as python dictionary
without surrounding braces to define() function of
problem description file
-o filename basename of output file(s) [default: <basename of
input file>]
--format=format output file format, one of: {vtk, h5} [default: vtk]
--save-restart=mode if given, save restart files according to the given
mode.
--load-restart=filename
if given, load the given restart file
--log=file log all messages to specified file (existing file will
be overwritten!)
-q, --quiet do not print any messages to screen
--save-ebc save a zero solution with applied EBCs (Dirichlet
boundary conditions)
--save-ebc-nodes save a zero solution with added non-zeros in EBC
(Dirichlet boundary conditions) nodes - scalar
variables are shown using colors, vector variables
using arrows with non-zero components corresponding to
constrained components
--save-regions save problem regions as meshes
--save-regions-as-groups
save problem regions in a single mesh but mark them by
using different element/node group numbers
--save-field-meshes save meshes of problem fields (with extra DOF nodes)
--solve-not do not solve (use in connection with --save-*)
--list=what list data, what can be one of: {terms, solvers}

Additional (stand-alone) examples are in the examples/ directory, e.g.:

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

Run:

./simple.pyexamples/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.

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

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]):

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.

The cell kind is the most general and should be used with volume
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 (surface 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
(volume integration domains, surface domains, etc.) and could be uses 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.

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.

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;

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

Only the 'shifted_periodic' LCBC needs the second region and the DOF
mapping function, see below.

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

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.

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

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.

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:

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

defget_pars_special(ts,coors,mode=None,extra_arg=None):ifmode=='special':ifextra_arg=='hello!':ic=0else:ic=1return{('x_%s'%ic):coors[:,ic]}functions={'get_pars1':(lambdats,coors,mode=None,**kwargs:get_pars_special(ts,coors,mode,extra_arg='hello!'),),}# Just another way of adding a function, besides 'functions' keyword.function_1={'name':'get_pars2','function':lambdats,coors,mode=None,**kwargs:get_pars_special(ts,coors,mode,extra_arg='hi!'),}

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':'h5',# 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>',# 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>',}

post_process_hook enables computing derived quantities, like
stress or strain, from the primary unknown variables. See the
examples in examples/large_deformation/ directory.

parametric_hook makes it possible to run parametric studies by
modifying the problem description programmatically. See
examples/diffusion/poisson_parametric_study.py for an example.

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 :

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.

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

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

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):

defpost_process(out,problem,state,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. state : State 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 ------- out : dict The updated output dictionary. """fromsfepy.base.baseimportStruct# 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)returnout

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 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 return the probing results, as well as a
handle to a corresponding matplotlib figure. See
linear_elasticity/its2D_4.py for additional explanation.

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.

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.

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

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 verticesofsurface region selector does not work in parallel,
because the region definition is applied to a task-local domain.

Isogeometric analysis (IGA) is a recently developed computational approach
that allows using the NURBS-based domain 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
in [1]. 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.

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.

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.

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.