# Useful Code Snippets and FAQ¶

## Miscellaneous¶

1. No module named ‘sfepy.discrete.common.extmods.mappings’.

When installing SfePy from sources or using the git version, its extension modules have to be compiled before using the package, see Compilation of C Extension Modules.

2. Finite element approximation (field) order and numerical quadrature order.

SfePy supports reading only straight-facet (linear approximation) meshes, nevertheless field orders higher than one can be used, because internally, the mesh elements are enriched with the required additional nodes. The calculation then occurs on such an augmented mesh with appropriate higher order elements.

The quadrature order equal to two-times the field order (used in many examples) works well for bilinear forms with constant (on each element) material parameters. For example, a dot product involves integrating ```u * v```, so if the approximation order of `u` and `v` is 1, their product’s order is 2. Of course, there are terms that could use a lower quadrature order, or higher, depending on the data. Increased quadrature order is required e.g. in terms with highly oscillating material coefficients.

Example:

```approx_order = 2
# The finite element approximation order.
fields = {
'displacement': ('real', 3, 'Omega', approx_order),
}
integrals = {
'i' : 2 * approx_order,
}
```
3. Higher order DOF visualization when using an approximation order greater than one.

By default, the additional, higher order DOFs, are not used in the VTK/HDF5 results files (`'strip'` linearization kind). To see the influence of those DOFs, `'adaptive'` linearization has to be used, see diffusion/sinbc.py (declarative API) and diffusion/laplace_refine_interactive.py or multi_physics/biot_parallel_interactive.py (imperative API, search `linearization`).

4. Numbering of DOFs.

Locally (in a connectivity row), the DOFs are stored DOF-by-DOF (`u_0` in all local nodes, `u_1` in all local nodes, …).

Globally (in a state vector), the DOFs are stored node-by-node (```u_0, u_1, ..., u_X``` in node 0, `u_0, u_1, ..., u_X` in node 1, …).

See also `create_adof_conn()`.

5. How to work with solvers/preconditioners?

See multi_physics/biot_short_syntax.py (user-defined preconditioners) or navier_stokes/stokes_slip_bc.py (petsc solver setup).

6. How to get the linear system components: the matrix and the right-hand side?

To get the residual vector `r` (see Implementation of Essential Boundary Conditions) and the tangent matrix `K`, the imperative API can be used as follows:

```# pb is a Problem instance,
pb.time_update()
state = pb.create_state()
state.apply_ebc()
r = pb.equations.eval_residuals(state())
K = pb.equations.eval_tangent_matrices(state(), pb.mtx_a)
```

## Material Parameters¶

1. How to set material parameters per region in the interactive mode (imperative API)?

Example: define `rho`, `D` to have different values in regions `omega1`, `omega2`:

```m = Material('m', values={'rho': {'omega1': 2700, 'omega2': 6000},
'D': {'omega1': D1, 'omega2': D2}})
```
2. How to implement state dependent materials?

Besides writing a custom solver, one can use pseudo-time-stepping for this purpose, as demonstrated in linear_elasticity/material_nonlinearity.py or diffusion/poisson_field_dependent_material.py. Note that the examples are contrived, and in practice care must be taken to ensure convergence.

3. Why are results of a 2D elasticity simulation not consistent with a properly constrained 3D elasticity simulation?

Possible reason: when using the Young’s modulus and Poisson’s ratio as input parameters, and then calling `stiffness_from_youngpoisson()`, note that the default value of the `plane` argument is `'strain'`, corresponding to the plane strain assumption, see also `lame_from_youngpoisson()`. Try setting `plane='stress'`.

4. How to set (time-dependent) material parameters by a function in the interactive mode (imperative API)?

Example (also showing the full material function signature):

```from sfepy.discrete import Material, Function

def get_pars(ts, coors, mode=None,
equations=None, term=None, problem=None, **kwargs):
value1 = a_function(ts.t, coors)
value2 = another_function(ts.step, coors)
if mode == 'qp':
out = {
'value1' : value1.reshape(coors.shape, 1, 1),
'value2' : value2.reshape(coors.shape, 1, 1),
}
return out
m = Material('m', function=Function('get_pars', get_pars))
```
5. How to get cells corresponding to coordinates in a material function?

The full signature of the material function is:

```def get_pars(ts, coors, mode=None,
equations=None, term=None, problem=None, **kwargs)
```

Thus it has access to `term.region.cells`, hence access to the cells that correspond to the coordinates. The length of the `coors` is `n_cell * n_qp`, where `n_qp` is the number of quadrature points per cell, and `n_cell = len(term.region.cells)`, so that `coors.reshape((n_cell, n_qp, -1))` can be used.