diff --git a/docs/examples.rst b/docs/examples.rst
new file mode 100644
index 00000000..23f82b26
--- /dev/null
+++ b/docs/examples.rst
@@ -0,0 +1,11 @@
+Examples
+========
+
+.. toctree::
+   :maxdepth: 1
+   :caption: Tutorials:
+   
+   examples/beam
+   examples/numcont
+   examples/platewithhole
+   examples/shear
diff --git a/docs/tutorial/beam.rst b/docs/examples/beam.rst
similarity index 81%
rename from docs/tutorial/beam.rst
rename to docs/examples/beam.rst
index b322ac67..9d6ffe00 100644
--- a/docs/tutorial/beam.rst
+++ b/docs/examples/beam.rst
@@ -18,27 +18,27 @@ First, let's create a meshed cube out of hexahedron cells with ``n=(181, 9, 9)``
 ..  code-block:: python
 
     import numpy as np
-    import felupe as fe
+    import felupe as fem
 
-    cube = fe.Cube(a=(0, 0, 0), b=(2000, 100, 100), n=(181, 9, 9))
-    region = fe.RegionHexahedron(cube)
-    displacement = fe.Field(region, dim=3)
-    field = fe.FieldContainer([displacement])
+    cube = fem.Cube(a=(0, 0, 0), b=(2000, 100, 100), n=(181, 9, 9))
+    region = fem.RegionHexahedron(cube)
+    displacement = fem.Field(region, dim=3)
+    field = fem.FieldContainer([displacement])
 
 
 A fixed boundary condition is applied on the left end of the beam. The degrees of freedom are partitioned into active (``dof1``) and fixed or inactive (``dof0``) degrees of freedom.
 
 ..  code-block:: python
 
-    bounds = {"fixed": fe.dof.Boundary(displacement, fx=lambda x: x==0)}
-    dof0, dof1 = fe.dof.partition(field, bounds)
+    bounds = {"fixed": fem.dof.Boundary(displacement, fx=lambda x: x==0)}
+    dof0, dof1 = fem.dof.partition(field, bounds)
 
 
 The material behavior is defined through a built-in isotropic linear-elastic material formulation.
 
 ..  code-block:: python
 
-    umat = fe.LinearElastic(E=206000, nu=0.3)
+    umat = fem.LinearElastic(E=206000, nu=0.3)
     density = 7850 * 1e-12
 
 
@@ -53,7 +53,7 @@ The body force is now assembled. Note that the gravity vector has to be reshaped
 
     gravity = np.array([0, 0, 9.81]) * 1e3
 
-    bodyforce = fe.IntegralForm(
+    bodyforce = fem.IntegralForm(
         fun=[density * gravity.reshape(-1, 1, 1)], 
         v=field, 
         dV=region.dV,
@@ -69,7 +69,7 @@ The weak form of linear elasticity is assembled into the stiffness matrix, where
 
 ..  code-block:: python
     
-    stiffness = fe.IntegralForm(
+    stiffness = fem.IntegralForm(
         fun=umat.elasticity(), 
         v=field, 
         dV=region.dV, 
@@ -80,10 +80,10 @@ The linear equation system may now be solved. First, a partition into active and
 
 ..  code-block:: python
 
-    system = fe.solve.partition(field, stiffness, dof1, dof0, r=-bodyforce)
-    field += fe.solve.solve(*system)
+    system = fem.solve.partition(field, stiffness, dof1, dof0, r=-bodyforce)
+    field += fem.solve.solve(*system)
 
-    fe.save(region, field, filename="bodyforce.vtk")
+    fem.save(region, field, filename="bodyforce.vtk")
 
 
 .. image:: images/beam_bodyforce.png
diff --git a/docs/tutorial/images/animation.ogv b/docs/examples/images/animation.ogv
similarity index 100%
rename from docs/tutorial/images/animation.ogv
rename to docs/examples/images/animation.ogv
diff --git a/docs/tutorial/images/beam.png b/docs/examples/images/beam.png
similarity index 100%
rename from docs/tutorial/images/beam.png
rename to docs/examples/images/beam.png
diff --git a/docs/tutorial/images/beam_bodyforce.png b/docs/examples/images/beam_bodyforce.png
similarity index 100%
rename from docs/tutorial/images/beam_bodyforce.png
rename to docs/examples/images/beam_bodyforce.png
diff --git a/docs/tutorial/images/getting-started.svg b/docs/examples/images/getting-started.svg
similarity index 100%
rename from docs/tutorial/images/getting-started.svg
rename to docs/examples/images/getting-started.svg
diff --git a/docs/tutorial/images/lab.png b/docs/examples/images/lab.png
similarity index 100%
rename from docs/tutorial/images/lab.png
rename to docs/examples/images/lab.png
diff --git a/docs/tutorial/images/platehole.svg b/docs/examples/images/platehole.svg
similarity index 100%
rename from docs/tutorial/images/platehole.svg
rename to docs/examples/images/platehole.svg
diff --git a/docs/tutorial/images/platewithhole_mesh.png b/docs/examples/images/platewithhole_mesh.png
similarity index 100%
rename from docs/tutorial/images/platewithhole_mesh.png
rename to docs/examples/images/platewithhole_mesh.png
diff --git a/docs/tutorial/images/platewithhole_stress.png b/docs/examples/images/platewithhole_stress.png
similarity index 100%
rename from docs/tutorial/images/platewithhole_stress.png
rename to docs/examples/images/platewithhole_stress.png
diff --git a/docs/tutorial/images/platewithhole_stressplot.png b/docs/examples/images/platewithhole_stressplot.png
similarity index 100%
rename from docs/tutorial/images/platewithhole_stressplot.png
rename to docs/examples/images/platewithhole_stressplot.png
diff --git a/docs/tutorial/images/platewithhole_stressplot2.png b/docs/examples/images/platewithhole_stressplot2.png
similarity index 100%
rename from docs/tutorial/images/platewithhole_stressplot2.png
rename to docs/examples/images/platewithhole_stressplot2.png
diff --git a/docs/tutorial/images/plot_force-displacement.png b/docs/examples/images/plot_force-displacement.png
similarity index 100%
rename from docs/tutorial/images/plot_force-displacement.png
rename to docs/examples/images/plot_force-displacement.png
diff --git a/docs/tutorial/images/shear.svg b/docs/examples/images/shear.svg
similarity index 100%
rename from docs/tutorial/images/shear.svg
rename to docs/examples/images/shear.svg
diff --git a/docs/tutorial/images/shear_deformed.png b/docs/examples/images/shear_deformed.png
similarity index 100%
rename from docs/tutorial/images/shear_deformed.png
rename to docs/examples/images/shear_deformed.png
diff --git a/docs/tutorial/images/shear_deformed_old.png b/docs/examples/images/shear_deformed_old.png
similarity index 100%
rename from docs/tutorial/images/shear_deformed_old.png
rename to docs/examples/images/shear_deformed_old.png
diff --git a/docs/tutorial/images/shear_mesh.png b/docs/examples/images/shear_mesh.png
similarity index 100%
rename from docs/tutorial/images/shear_mesh.png
rename to docs/examples/images/shear_mesh.png
diff --git a/docs/tutorial/images/shear_plot.svg b/docs/examples/images/shear_plot.svg
similarity index 100%
rename from docs/tutorial/images/shear_plot.svg
rename to docs/examples/images/shear_plot.svg
diff --git a/docs/tutorial/numcont.rst b/docs/examples/numcont.rst
similarity index 100%
rename from docs/tutorial/numcont.rst
rename to docs/examples/numcont.rst
diff --git a/docs/tutorial/platewithhole.rst b/docs/examples/platewithhole.rst
similarity index 90%
rename from docs/tutorial/platewithhole.rst
rename to docs/examples/platewithhole.rst
index e0089c20..4a1c232b 100644
--- a/docs/tutorial/platewithhole.rst
+++ b/docs/examples/platewithhole.rst
@@ -21,7 +21,7 @@ Let's create a meshed plate with a hole out of quad cells with the help of ``pyg
 
 ..  code-block:: python
 
-    import felupe as fe
+    import felupe as fem
     import pygmsh
 
     h = 1.0
@@ -46,7 +46,7 @@ The points and cells of the above mesh are used to initiate a FElupe mesh.
 
 ..  code-block:: python
 
-    mesh = fe.Mesh(
+    mesh = fem.Mesh(
         points=mesh.points[:, :2], 
         cells=mesh.cells[1].data, 
         cell_type=mesh.cells[1].type
@@ -60,11 +60,11 @@ A numeric quad-region created on the mesh in combination with a vector-valued di
 
 ..  code-block:: python
 
-    region = fe.RegionQuad(mesh)
-    displacement = fe.Field(region, dim=2)
-    field = fe.FieldContainer([displacement])
+    region = fem.RegionQuad(mesh)
+    displacement = fem.Field(region, dim=2)
+    field = fem.FieldContainer([displacement])
 
-    boundaries, loadcase = fe.dof.uniaxial(
+    boundaries, loadcase = fem.dof.uniaxial(
         field, move=0.001, right=L, clamped=False
     )
 
@@ -73,7 +73,7 @@ The material behavior is defined through a built-in isotropic linear-elastic mat
 
 ..  code-block:: python
 
-    umat = fe.LinearElasticPlaneStress(E=210000, nu=0.3)
+    umat = fem.LinearElasticPlaneStress(E=210000, nu=0.3)
     
 
 The weak form of linear elasticity is assembled into the stiffness matrix, where the constitutive elasticity matrix is generated with :func:`umat.hessian` (or the alias :func:`umat.elasticity`).
@@ -85,7 +85,7 @@ The weak form of linear elasticity is assembled into the stiffness matrix, where
 
 ..  code-block:: python
 
-    K = fe.IntegralForm(
+    K = fem.IntegralForm(
         fun=umat.elasticity(), 
         v=field, 
         dV=region.dV, 
@@ -100,8 +100,8 @@ The linear equation system may now be solved. First, a partition into active and
     dof0 = loadcase["dof0"]
     ext0 = loadcase["ext0"]
 
-    system = fe.solve.partition(field, K, dof1, dof0)
-    field += fe.solve.solve(*system, ext0)
+    system = fem.solve.partition(field, K, dof1, dof0)
+    field += fem.solve.solve(*system, ext0)
 
 Let's evaluate the deformation gradient from the displacement field and calculate the stress tensor. This process is also called *stress recovery*.
 
@@ -126,8 +126,8 @@ However, the stress results are still located at the numeric integration points.
         stress[0, 0] * stress[1, 1]
     )
     
-    stress_projected = fe.project(stress, region)
-    vonmises_projected = fe.project(vonmises, region)
+    stress_projected = fem.project(stress, region)
+    vonmises_projected = fem.project(vonmises, region)
 
 
 Results are saved as VTK-files, where additional point-data is passed within the ``point_data`` argument. Stresses are normalized by the mean value of the stress at the right end-face in order to visualize a normalized stress distribution over the plate.
@@ -136,7 +136,7 @@ Results are saved as VTK-files, where additional point-data is passed within the
 
     right = mesh.points[:, 0] == L
 
-    fe.save(
+    fem.save(
         region, 
         field,
         filename="plate_with_hole.vtk",
diff --git a/docs/tutorial/shear.rst b/docs/examples/shear.rst
similarity index 89%
rename from docs/tutorial/shear.rst
rename to docs/examples/shear.rst
index 0153ba6a..d57c9837 100644
--- a/docs/tutorial/shear.rst
+++ b/docs/examples/shear.rst
@@ -34,7 +34,7 @@ degrees of freedom. Hence, we have to drop our MPC-centerpoint from that list.
 ..  code-block:: python
 
     import numpy as np
-    import felupe as fe
+    import felupe as fem
 
     H = 10
     L = 20
@@ -43,7 +43,7 @@ degrees of freedom. Hence, we have to drop our MPC-centerpoint from that list.
     n = 21
     a = min(L / n, H / n)
     
-    mesh = fe.Rectangle((0, 0), (L, H), n=(round(L / a), round(H / a)))
+    mesh = fem.Rectangle((0, 0), (L, H), n=(round(L / a), round(H / a)))
     mesh.points = np.vstack((mesh.points, [0, 2 * H]))
     mesh.update(mesh.cells)
     mesh.points_without_cells = np.array([], dtype=bool)
@@ -56,18 +56,18 @@ displacement field for plane-strain as well as scalar-valued fields for the hydr
 
 ..  code-block:: python
 
-    region = fe.RegionQuad(mesh)
-    fields = fe.FieldsMixed(region, n=3, planestrain=True)
+    region = fem.RegionQuad(mesh)
+    fields = fem.FieldsMixed(region, n=3, planestrain=True)
     
     f0 = lambda y: np.isclose(y, 0)
     f2 = lambda y: np.isclose(y, 2* H)
     
     boundaries = {
-        "fixed": fe.Boundary(fields[0], fy=f0),
-        "control": fe.Boundary(fields[0], fy=f2, skip=(0, 1)),
+        "fixed": fem.Boundary(fields[0], fy=f0),
+        "control": fem.Boundary(fields[0], fy=f2, skip=(0, 1)),
     }
     
-    dof0, dof1 = fe.dof.partition(fields, boundaries)
+    dof0, dof1 = fem.dof.partition(fields, boundaries)
 
 
 The micro-sphere material formulation is used for the rubber. It is defined
@@ -97,7 +97,7 @@ as a hyperelastic material in matADi. The material formulation is finally applie
         bulk=5000.0,
     )
     
-    rubber = fe.SolidBody(umat=mat.ThreeFieldVariation(umat), field=fields)
+    rubber = fem.SolidBody(umat=mat.ThreeFieldVariation(umat), field=fields)
 
 At the centerpoint of a multi-point constraint (MPC) the external shear
 movement is prescribed. It also ensures a force-free top plate in direction 
@@ -105,7 +105,7 @@ movement is prescribed. It also ensures a force-free top plate in direction
 
 ..  code-block:: python
 
-    mpc = fe.MultiPointConstraint(
+    mpc = fem.MultiPointConstraint(
         field=fields,
         points=np.arange(mesh.npoints)[mesh.points[:, 1] == H],
         centerpoint=mesh.npoints - 1,
@@ -126,7 +126,7 @@ job. A job returns a generator object with the results of all substeps.
 
 ..  code-block:: python
 
-    UX = fe.math.linsteps([0, 15], 15)
+    UX = fem.math.linsteps([0, 15], 15)
     UY = []
     FX = []
 
@@ -146,12 +146,12 @@ job. A job returns a generator object with the results of all substeps.
 
 ..  code-block:: python
     
-    step = fe.Step(
+    step = fem.Step(
         items=[rubber, mpc], 
         ramp={boundaries["control"]: UX}, 
         boundaries=boundaries
     )
-    job = fe.Job(steps=[step], callback=callback)
+    job = fem.Job(steps=[step], callback=callback)
     res = job.evaluate()
 
 For the maximum deformed model a VTK-file containing principal stretches
@@ -164,9 +164,9 @@ projected to mesh points is exported.
     F = fields[0].extract()
     C = dot(transpose(F), F)
     
-    stretches = fe.project(np.sqrt(eigh(C)[0]), region)
+    stretches = fem.project(np.sqrt(eigh(C)[0]), region)
     
-    fe.save(region, fields, point_data={
+    fem.save(region, fields, point_data={
         "Maximum-principal-stretch": np.max(stretches, axis=1),
         "Minimum-principal-stretch": np.min(stretches, axis=1),
     })
diff --git a/docs/felupe/global/constitution.rst b/docs/felupe/global/constitution.rst
index cf21042d..a34b960d 100644
--- a/docs/felupe/global/constitution.rst
+++ b/docs/felupe/global/constitution.rst
@@ -2,49 +2,55 @@ Constitution
 ~~~~~~~~~~~~
 
 .. autoclass:: felupe.NeoHooke
-   :noindex:
    :members:
    :undoc-members:
    :inherited-members:
 
 .. autoclass:: felupe.LinearElastic
-   :noindex:
    :members:
    :undoc-members:
    :inherited-members:
 
 .. autoclass:: felupe.LinearElasticPlaneStress
-   :noindex:
    :members:
    :undoc-members:
    :inherited-members:
 
 .. autoclass:: felupe.LinearElasticPlaneStrain
-   :noindex:
+   :members:
+   :undoc-members:
+   :inherited-members:
+
+.. autoclass:: felupe.LinearElasticPlasticIsotropicHardening
    :members:
    :undoc-members:
    :inherited-members:
 
 .. autoclass:: felupe.ThreeFieldVariation
-   :noindex:
    :members:
    :undoc-members:
    :inherited-members:
 
+.. autoclass:: felupe.UserMaterialStrain
+   :members:
+   :undoc-members:
+   :inherited-members:
+
+.. autofunction:: felupe.constitution.linear_elastic
+
+.. autofunction:: felupe.constitution.linear_elastic_plastic_isotropic_hardening
+
 .. autoclass:: felupe.LineChange
-   :noindex:
    :members:
    :undoc-members:
    :inherited-members:
 
 .. autoclass:: felupe.AreaChange
-   :noindex:
    :members:
    :undoc-members:
    :inherited-members:
 
 .. autoclass:: felupe.VolumeChange
-   :noindex:
    :members:
    :undoc-members:
    :inherited-members:
\ No newline at end of file
diff --git a/docs/felupe/global/mechanics.rst b/docs/felupe/global/mechanics.rst
index 90457ca9..f84d8d4f 100644
--- a/docs/felupe/global/mechanics.rst
+++ b/docs/felupe/global/mechanics.rst
@@ -6,7 +6,7 @@ Mechanics
    :undoc-members:
    :show-inheritance:
 
-.. autoclass:: felupe.SolidBodyTensor
+.. autoclass:: felupe.SolidBodyNearlyIncompressible
    :members:
    :undoc-members:
    :show-inheritance:
diff --git a/docs/howto.rst b/docs/howto.rst
index b2d51093..f474b44c 100644
--- a/docs/howto.rst
+++ b/docs/howto.rst
@@ -5,6 +5,7 @@ How-To Guides
    :maxdepth: 1
    :caption: How-To:
    
+   howto/building_blocks
    howto/meshgen
    howto/rbe2
    howto/mixed
@@ -12,4 +13,5 @@ How-To Guides
    howto/solvers
    howto/composite
    howto/forms
-   howto/solid
\ No newline at end of file
+   howto/solid
+   howto/umat
\ No newline at end of file
diff --git a/docs/howto/axi.rst b/docs/howto/axi.rst
index f529dddb..94daa5e8 100644
--- a/docs/howto/axi.rst
+++ b/docs/howto/axi.rst
@@ -5,12 +5,12 @@ For axisymmetric analyses an axisymmetric vector-valued field has to be created
 
 ..  code-block:: python
 
-    import felupe as fe
+    import felupe as fem
 
-    mesh = fe.Rectangle(n=3)
-    region = fe.RegionQuad(mesh)
-    u = fe.FieldAxisymmetric(region, dim=2)
-    field = fe.FieldContainer([u])
+    mesh = fem.Rectangle(n=3)
+    region = fem.RegionQuad(mesh)
+    u = fem.FieldAxisymmetric(region, dim=2)
+    field = fem.FieldContainer([u])
 
 Now it gets important: The 3x3 deformation gradient for an axisymmetric problem is obtained with :meth:`felupe.FieldAxisymmetric.grad` or :meth:`felupe.FieldAxisymmetric.extract` methods. For instances of :class:`felupe.FieldAxisymmetric` the gradient is modified to return a 3x3 gradient as described in :ref:`theory-axi`.
 
@@ -22,7 +22,7 @@ For simplicity, let's assume a (built-in) Neo-Hookean material.
 
 ..  code-block:: python
 
-    umat = fe.NeoHooke(mu=1, bulk=5)
+    umat = fem.NeoHooke(mu=1, bulk=5)
 
 
 FElupe provides an adopted :class:`felupe.IntegralFormAxisymmetric` class for the integration and the sparse matrix assemblage of axisymmetric problems under hood. It uses the additional information (e.g. radial coordinates at integration points) stored in :class:`felupe.FieldAxisymmetric` to provide a consistent interface in comparison to default IntegralForms.
@@ -31,7 +31,7 @@ FElupe provides an adopted :class:`felupe.IntegralFormAxisymmetric` class for th
 
     dA = region.dV
 
-    r = fe.IntegralForm(umat.gradient(F), field, dA).assemble()
-    K = fe.IntegralForm(umat.hessian(F), field, dA, field).assemble()
+    r = fem.IntegralForm(umat.gradient(F), field, dA).assemble()
+    K = fem.IntegralForm(umat.hessian(F), field, dA, field).assemble()
 
 To sum up, for axisymmetric problems use :class:`felupe.FieldAxisymmetric`. Of course, Mixed-field formulations may be applied on axisymmetric scenarios too.
diff --git a/docs/tutorial/gettingstarted.rst b/docs/howto/building_blocks.rst
similarity index 96%
rename from docs/tutorial/gettingstarted.rst
rename to docs/howto/building_blocks.rst
index d985be72..ffefaa68 100644
--- a/docs/tutorial/gettingstarted.rst
+++ b/docs/howto/building_blocks.rst
@@ -1,6 +1,6 @@
-.. _tutorial-getting-started:
+.. _howto-building-blocks:
 
-Getting started
+Building Blocks
 ---------------
 
 .. admonition:: Have a look at the building blocks of FElupe.
@@ -143,7 +143,7 @@ The integral (or weak) forms of equilibrium equations are defined by the :class:
 
 ..  code-block:: python
 
-    linearform   = felupe.IntegralForm(P(F), field, dV, grad_v=[True])
+    linearform = felupe.IntegralForm(P(F)[:-1], field, dV, grad_v=[True])
     bilinearform = felupe.IntegralForm(A(F), field, dV, u=field, grad_v=[True], grad_u=[True])
 
 
@@ -194,7 +194,7 @@ A very simple newton-rhapson code looks like this:
     for iteration in range(8):
         F = field.extract()
 
-        linearform = felupe.IntegralForm(P(F), field, dV)
+        linearform = felupe.IntegralForm(P(F)[:-1], field, dV)
         bilinearform = felupe.IntegralForm(A(F), field, dV, field)
 
         r = linearform.assemble()
@@ -220,11 +220,11 @@ A very simple newton-rhapson code looks like this:
     4 6.017153213511068e-09
     5 5.675484825228616e-16
 
-Alternatively, one may also use the Newton-Rhapson procedure as shown in :ref:`tutorial-hello-felupe`.
+Alternatively, one may also use the Newton-Rhapson function of FElupe.
 
 ..  code-block:: python
 
-    res = fe.newtonrhapson(field, umat=umat, dof1=dof1, dof0=dof0, ext0=ext0)
+    res = fem.newtonrhapson(field, umat=umat, dof1=dof1, dof0=dof0, ext0=ext0)
     field = res.x
 
 
diff --git a/docs/howto/composite.rst b/docs/howto/composite.rst
index 5d7d7531..217e9b1c 100644
--- a/docs/howto/composite.rst
+++ b/docs/howto/composite.rst
@@ -1,107 +1,90 @@
 Composite Regions with Solid Bodies
 -----------------------------------
 
-This section demonstrates how to set up a problem with two regions, each associated to a seperated solid body. First, a (total) region is created.
+This section demonstrates how to set up a problem with two regions, each associated to a seperated solid body.
 
 ..  code-block:: python
 
-    import felupe as fe
+    import felupe as fem
     import numpy as np
 
-    n = 16
-    mesh = fe.Cube(n=n)
-    region = fe.RegionHexahedron(mesh)
+    m = fem.Rectangle(n=21)
 
 
-In a second step, sub-sets for points and cells are created from which two sub-regions and sub-fields are initiated.
+In a second step, sub-meshes are created.
     
 ..  code-block:: python
 
-    points = np.arange(mesh.npoints)[np.logical_or.reduce((
-        mesh.points[:,0] == 0,
-        mesh.points[:,0] == 0 + 1/(n - 1),
-        mesh.points[:,0] == 0.5 - 1/(n - 1) / 2,
-        mesh.points[:,0] == 0.5 + 1/(n - 1) / 2,
-        mesh.points[:,0] == 1 - 1/(n - 1),
-        mesh.points[:,0] == 1,
-    ))]
-    cells = np.isin(mesh.cells, points).sum(1) == mesh.cells.shape[1]
-
-    mesh_rubber = mesh.copy()
-    mesh_rubber.update(mesh_rubber.cells[~cells])
-
-    mesh_steel = mesh.copy()
-    mesh_steel.update(mesh_steel.cells[cells])
+    # take some points from the inside for the fiber-reinforced area
+    eps = 1e-3
+    mask = np.arange(m.npoints)[np.logical_and.reduce([
+        m.points[:, 0] >= 0.3,
+        m.points[:, 0] <= 0.7 + eps,
+        m.points[:, 1] >= 0.3,
+        m.points[:, 1] <= 0.7 + eps,
+    ])]
     
-    region_rubber = fe.RegionHexahedron(mesh_rubber)
-    field_rubber = fe.FieldsMixed(region_rubber, n=3)
-
-    region_steel = fe.RegionHexahedron(mesh_steel)
-    field_steel = fe.FieldsMixed(region_steel, n=1)
+    # copies of the mesh
+    mesh = [m.copy(), m.copy()]
+    
+    # create sub-meshes (fiber, matrix)
+    mesh[0].update(m.cells[ np.all(np.isin(m.cells, mask), axis=1)])
+    mesh[1].update(m.cells[~np.all(np.isin(m.cells, mask), axis=1)])
 
-This is followed by the creation of a global (mixed) field. Note that is approach is only valid for a simulation model with one body formulated on mixed fields.
+This is followed by the creation of a global region/field and two sub-regions/sub-fields.
 
 ..  code-block:: python
     
-    pressure, volumeratio = field_rubber[1:]
-    field = fe.FieldContainer([fe.Field(region, dim=3), pressure, volumeratio])
+    # a global and two sub-regions
+    region = fem.RegionQuad(m)
+    regions = [fem.RegionQuad(me) for me in mesh]
+    
+    # a global and two sub-fields
+    field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)])
+    fields = [
+        fem.FieldContainer([fem.FieldPlaneStrain(region[0], dim=2)]),
+        fem.FieldContainer([fem.FieldPlaneStrain(region[1], dim=2)]),
+    ]
     
 The displacement boundaries are created on the total field.
 
 ..  code-block:: python
 
-    boundaries, loadcase = fe.dof.uniaxial(field, move=-0.1)
+    boundaries = dict(
+        fixed=fem.Boundary(field[0], fx=0),
+        move=fem.Boundary(field[0], fx=1),
+    )
 
 
 The rubber is associated to a Neo-Hookean material formulation whereas the steel is modeled by a linear elastic material formulation. For each material a solid body is created.
 
 ..  code-block:: python
 
-    neohooke = fe.ThreeFieldVariation(fe.NeoHooke(mu=1.0, bulk=5000.0))
-    linearelastic = fe.LinearElastic(E=210000.0, nu=0.3)
-    
-    rubber = fe.SolidBody(neohooke, field_rubber)
-    steel = fe.SolidBody(linearelastic, field_steel)
-
-
-Inside the Newton-Rhapson iterations both the internal force vector and the tangent stiffness matrix are assembled and summed up from contributions of both solid bodies. All field values are linked, that means they share their values array.
-
-..  code-block:: python
-
-    res = fe.newtonrhapson(field, items=[rubber, steel], **loadcase)
-
-..  code-block:: shell
-
-    Newton-Rhapson solver
-    =====================
+    # two material model formulations
+    neo_hooke = fem.NeoHooke(mu=1, bulk=1)
+    linear_elastic = fem.LinearElastic(E=100, nu=0.3)
     
-    | # |  norm(dx) |
-    |---|-----------|
-    | 1 | 9.075e+01 |
-    | 2 | 1.370e+01 |
-    | 3 | 6.998e-01 |
-    | 4 | 1.105e-02 |
-    | 5 | 1.658e-06 |
-    | 6 | 5.446e-13 |
-    
-    Solution converged in 6 iterations within 63.52 seconds.
+    # the solid bodies
+    fiber = fem.SolidBody(umat=linear_elastic, field=fields[0])
+    matrix = fem.SolidBody(umat=neo_hooke, field=fields[1])
 
-Results may be exported either for the total region or with stresses for sub-regions only.
 
-.. image:: images/composite_total.png
-   :width: 600px
+A step is created and further added to a job. The global field must be passed to the ``x0`` argument during the evaluation of the job. Internally, all field values are linked automatically, i.e. they share their ``values`` attribute.
 
 ..  code-block:: python
 
-    s = rubber.evaluate.cauchy_stress()
-    cauchy_stress = fe.project(fe.math.tovoigt(s), region_rubber)
+    # prepare a step with substeps
+    move = fem.math.linsteps([0, 0.5], num=10)
+    step = fem.Step(
+        items=[matrix, fiber],
+        ramp={boundaries["move"]: move}, 
+        boundaries=boundaries
+    )
     
-    fe.save(region, res.x, filename="result.vtk")
+    # take care of the x0-argument
+    job = fem.Job(steps=[step])
+    job.evaluate(x0=field, filename="result.xdmf")
 
-    fe.save(region_rubber, field_rubber,
-        filename="result_rubber.vtk", 
-        point_data={"CauchyStress": cauchy_stress}
-    )
 
-.. image:: images/composite_rubber_cauchy.png
+.. image:: images/composite.png
    :width: 600px
diff --git a/docs/howto/forms.rst b/docs/howto/forms.rst
index af04d37c..75b7c693 100644
--- a/docs/howto/forms.rst
+++ b/docs/howto/forms.rst
@@ -5,12 +5,12 @@ FElupe requires a pre-evaluated array for the definition of a bilinear :class:`f
 
 ..  code-block:: python
 
-    import felupe as fe
+    import felupe as fem
     
-    mesh = fe.Cube(n=11)
-    region = fe.RegionHexahedron(mesh)
-    displacement = fe.Field(region, dim=3)
-    field = fe.FieldContainer([displacement])
+    mesh = fem.Cube(n=11)
+    region = fem.RegionHexahedron(mesh)
+    displacement = fem.Field(region, dim=3)
+    field = fem.FieldContainer([displacement])
 
 The bilinear form of linear elasticity is defined as
 
diff --git a/docs/howto/images/composite.png b/docs/howto/images/composite.png
new file mode 100644
index 00000000..0233c29c
Binary files /dev/null and b/docs/howto/images/composite.png differ
diff --git a/docs/howto/images/composite_rubber_cauchy.png b/docs/howto/images/composite_rubber_cauchy.png
deleted file mode 100644
index 3bb84d79..00000000
Binary files a/docs/howto/images/composite_rubber_cauchy.png and /dev/null differ
diff --git a/docs/howto/images/composite_total.png b/docs/howto/images/composite_total.png
deleted file mode 100644
index 711bf197..00000000
Binary files a/docs/howto/images/composite_total.png and /dev/null differ
diff --git a/docs/tutorial/images/deformed_mesh.png b/docs/howto/images/deformed_mesh.png
similarity index 100%
rename from docs/tutorial/images/deformed_mesh.png
rename to docs/howto/images/deformed_mesh.png
diff --git a/docs/tutorial/images/region.svg b/docs/howto/images/region.svg
similarity index 100%
rename from docs/tutorial/images/region.svg
rename to docs/howto/images/region.svg
diff --git a/docs/tutorial/images/undeformed_mesh.png b/docs/howto/images/undeformed_mesh.png
similarity index 100%
rename from docs/tutorial/images/undeformed_mesh.png
rename to docs/howto/images/undeformed_mesh.png
diff --git a/docs/howto/meshgen.rst b/docs/howto/meshgen.rst
index 0a978659..eee59819 100644
--- a/docs/howto/meshgen.rst
+++ b/docs/howto/meshgen.rst
@@ -6,7 +6,7 @@ FElupe provides a simple mesh generation module :mod:`felupe.mesh`. A Mesh insta
 ..  code-block:: python
 
     import numpy as np
-    import felupe as fe
+    import felupe as fem
 
     points = np.array([
         [ 0, 0], # point 1
@@ -19,7 +19,7 @@ FElupe provides a simple mesh generation module :mod:`felupe.mesh`. A Mesh insta
         [ 0, 1, 3, 2], # point-connectivity of first cell
     ])
 
-    mesh = fe.Mesh(points, cells, cell_type="quad")
+    mesh = fem.Mesh(points, cells, cell_type="quad")
 
 .. image:: images/quad.png
    :width: 400px
@@ -32,9 +32,9 @@ First let's start with the generation of a line from ``x=1`` to ``x=3`` with ``n
 
 ..  code-block:: python
 
-    line = fe.mesh.Line(a=1, b=3, n=7)
-    rect = fe.mesh.expand(line, n=2, z=5)
-    hexa = fe.mesh.expand(rect, n=2, z=3)
+    line = fem.mesh.Line(a=1, b=3, n=7)
+    rect = fem.mesh.expand(line, n=2, z=5)
+    hexa = fem.mesh.expand(rect, n=2, z=3)
 
 .. image:: images/cube.png
    :width: 400px
@@ -43,7 +43,7 @@ First let's start with the generation of a line from ``x=1`` to ``x=3`` with ``n
 Lines, rectangles and cubes
 ***************************
 
-Of course lines, rectangles, cubes and cylinders do not have to be constructed manually each time. Instead, some easier to use classes are povided by FElupe like :class:`felupe.mesh.Line`, :class:`felupe.Rectangle` or :class:`felupe.Cube`.
+Of course lines, rectangles, cubes and cylinders do not have to be constructed manually each time. Instead, some easier to use classes are povided by FElupe like :class:`felupe.mesh.Line`, :class:`felupe.Rectangle` or :class:`felupe.Cube`. For non equi-distant points per axis use :class:`felupe.Grid`.
 
 Triangle and Tetrahedron meshes
 *******************************
@@ -52,11 +52,11 @@ Any quad or tetrahedron mesh may be subdivided (triangulated) to meshes out of T
 
 ..  code-block:: python
 
-    rectangle_quad = fe.Rectangle(n=5)
-    rectangle_tetra = fe.mesh.triangulate(rect_quad)
+    rectangle_quad = fem.Rectangle(n=5)
+    rectangle_tetra = fem.mesh.triangulate(rect_quad)
 
-    cube_hexahedron = fe.Cube(n=5)
-    cube_tetra = fe.mesh.triangulate(cube_hexahedron)
+    cube_hexahedron = fem.Cube(n=5)
+    cube_tetra = fem.mesh.triangulate(cube_hexahedron)
 
 Meshes with midpoints
 *********************
@@ -65,13 +65,13 @@ If a mesh with midpoints is required by a region, functions for edge, face and v
 
 ..  code-block:: python
     
-    rectangle_quad4 = fe.Rectangle(n=6)
-    rectangle_quad8 = fe.mesh.convert(rectangle_quad4, order=2)
-    rectangle_quad9 = fe.mesh.convert(rectangle_quad4, order=2, calc_midfaces=True)
+    rectangle_quad4 = fem.Rectangle(n=6)
+    rectangle_quad8 = fem.mesh.convert(rectangle_quad4, order=2)
+    rectangle_quad9 = fem.mesh.convert(rectangle_quad4, order=2, calc_midfaces=True)
 
 The same also applies on meshes with triangles.
 
 ..  code-block:: python
 
-    rectangle_triangle3 = fe.mesh.triangulate(fe.Rectangle(n=6))
-    rectangle_triangle6 = fe.mesh.add_midpoints_edges(rectangle_triangle3)
\ No newline at end of file
+    rectangle_triangle3 = fem.mesh.triangulate(fem.Rectangle(n=6))
+    rectangle_triangle6 = fem.mesh.add_midpoints_edges(rectangle_triangle3)
\ No newline at end of file
diff --git a/docs/howto/mixed.rst b/docs/howto/mixed.rst
index 925e2a85..0eb6b1b0 100644
--- a/docs/howto/mixed.rst
+++ b/docs/howto/mixed.rst
@@ -1,33 +1,33 @@
 Mixed-Field Problems
 ~~~~~~~~~~~~~~~~~~~~
 
-FElupe supports mixed-field formulations in a similar way it can handle (default) single-field variations. The definition of a mixed-field variation is shown for the hydrostatic-volumetric selective three-field-variation with independend fields for displacements :math:`\boldsymbol{u}`, pressure :math:`p` and volume ratio :math:`J`. The total potential energy for nearly-incompressible hyperelasticity is formulated with a determinant-modified deformation gradient. We take the :ref:`tutorial-getting-started` tutorial and modify it accordingly. The built-in Neo-Hookean material model is used as an argument of :class:`felupe.ThreeFieldVariation` for mixed-field problems.
+FElupe supports mixed-field formulations in a similar way it can handle (default) single-field formulations. The definition of a mixed-field formulation is shown for the hydrostatic-volumetric selective three-field-variation with independend fields for displacements :math:`\boldsymbol{u}`, pressure :math:`p` and volume ratio :math:`J`. The total potential energy for nearly-incompressible hyperelasticity is formulated with a determinant-modified deformation gradient. The built-in Neo-Hookean material model is used as an argument of :class:`felupe.ThreeFieldVariation` for mixed-field problems.
 
 ..  code-block:: python
 
-    import felupe as fe
+    import felupe as fem
 
-    neohooke = fe.constitution.NeoHooke(mu=1.0, bulk=5000.0)
-    umat = fe.constitution.ThreeFieldVariation(neohooke)
+    neohooke = fem.constitution.NeoHooke(mu=1.0, bulk=5000.0)
+    umat = fem.constitution.ThreeFieldVariation(neohooke)
 
 Next, let's create a meshed cube. Two regions, one for the displacements and another one for the pressure and the volume ratio are created.
 
 ..  code-block:: python
 
-    mesh  = fe.Cube(n=6)
+    mesh  = fem.Cube(n=6)
 
-    region  = fe.RegionHexahedron(mesh)
-    region0 = fe.RegionConstantHexahedron(mesh)
+    region  = fem.RegionHexahedron(mesh)
+    region0 = fem.RegionConstantHexahedron(mesh)
 
     dV = region.dV
 
-    displacement = fe.Field(region,  dim=3)
-    pressure     = fe.Field(region0, dim=1)
-    volumeratio  = fe.Field(region0, dim=1, values=1)
+    displacement = fem.Field(region,  dim=3)
+    pressure     = fem.Field(region0, dim=1)
+    volumeratio  = fem.Field(region0, dim=1, values=1)
 
-    field = fe.FieldContainer((displacement, pressure, volumeratio))
+    field = fem.FieldContainer(fields=[displacement, pressure, volumeratio])
 
-Boundary conditions are enforced in the same way as in Getting Started.
+Boundary conditions are enforced on the displacement field.
 
 ..  code-block:: python
 
@@ -35,14 +35,14 @@ Boundary conditions are enforced in the same way as in Getting Started.
 
     f1 = lambda x: np.isclose(x, 1)
 
-    boundaries = fe.dof.symmetry(displacement)
-    boundaries["right"] = fe.Boundary(displacement, fx=f1, skip=(1, 0, 0))
-    boundaries["move" ] = fe.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=-0.4)
+    boundaries = fem.dof.symmetry(displacement)
+    boundaries["right"] = fem.Boundary(displacement, fx=f1, skip=(1, 0, 0))
+    boundaries["move" ] = fem.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=-0.4)
 
-    dof0, dof1 = fe.dof.partition(field, boundaries)
-    ext0 = fe.dof.apply(field, boundaries, dof0)
+    dof0, dof1 = fem.dof.partition(field, boundaries)
+    ext0 = fem.dof.apply(field, boundaries, dof0)
 
-The Newton-Rhapson iterations are coded quite similar to the one used in :ref:`tutorial-getting-started`. FElupe provides a Mixed-field version of it's :class:`felupe.IntegralForm`, called :class:`felupe.IntegralFormMixed`. It assumes that the first field operates on the gradient and all the others don't. The resulting system vector with incremental values of the fields has to be splitted at the field-offsets in order to update the fields.
+The Newton-Rhapson iterations are coded quite similar. For mixed-fields, FElupe assumes that the first field operates on the gradient and all the others don't. The resulting system vector with incremental values of the fields has to be splitted at the field-offsets in order to update the fields.
 
 ..  code-block:: python
 
@@ -50,14 +50,14 @@ The Newton-Rhapson iterations are coded quite similar to the one used in :ref:`t
 
         F, p, J = field.extract()
         
-        linearform   = fe.IntegralForm(umat.gradient([F, p, J]), field, dV)
-        bilinearform = fe.IntegralForm(umat.hessian([F, p, J]), field, dV, field)
+        linearform = fem.IntegralForm(umat.gradient([F, p, J])[:-1], field, dV)
+        bilinearform = fem.IntegralForm(umat.hessian([F, p, J]), field, dV, field)
 
         r = linearform.assemble().toarray()[:, 0]
         K = bilinearform.assemble()
         
-        system = fe.solve.partition(field, K, dof1, dof0, r)
-        dfield = np.split(fe.solve.solve(*system, ext0), field.offsets)
+        system = fem.solve.partition(field, K, dof1, dof0, r)
+        dfield = np.split(fem.solve.solve(*system, ext0), field.offsets)
         
         field += dfield
 
@@ -67,7 +67,7 @@ The Newton-Rhapson iterations are coded quite similar to the one used in :ref:`t
         if norm < 1e-12:
             break
 
-    fe.tools.save(region, field, filename="result.vtk")
+    fem.tools.save(region, field, filename="result.vtk")
 
 The deformed cube is finally visualized by a VTK output file with the help of Paraview.
 
diff --git a/docs/howto/rbe2.rst b/docs/howto/rbe2.rst
index 79bcfd7c..968d370b 100644
--- a/docs/howto/rbe2.rst
+++ b/docs/howto/rbe2.rst
@@ -1,27 +1,27 @@
 Multi-Point Constraints
 -----------------------
 
-This How-To demonstrates the usage of multi-point constraints (also called MPC or RBE2 rigid-body-elements) with an independent centerpoint and one or more dependent points. First, a centerpoint has to be added to the mesh. MPC objects are supported as items of the Newton-Rhapson procedure.
+This How-To demonstrates the usage of multi-point constraints (also called MPC or RBE2 rigid-body-elements) with an independent centerpoint and one or more dependent points. First, a centerpoint has to be added to the mesh. MPC objects are supported as items of a Step and the Newton-Rhapson procedure.
 
 ..  code-block:: python
 
     import numpy as np
-    import felupe as fe
+    import felupe as fem
 
     # mesh with one additional rbe2-control point
-    mesh = fe.Cube(n=11)
+    mesh = fem.Cube(n=11)
     mesh.points = np.vstack((mesh.points, [2.0, 0.0, 0.0]))
     mesh.update(mesh.cells)
     
-    region = felupe.RegionHexahedron(mesh)
-    displacement = felupe.Field(region, dim=3)
-    field = felupe.FieldContainer([displacement])
+    region = fem.RegionHexahedron(mesh)
+    displacement = fem.Field(region, dim=3)
+    field = fem.FieldContainer([displacement])
 
 An instance of :class:`felupe.MultiPointConstraint` defines the multi-point constraint. This instance provides two methods, :meth:`felupe.MultiPointConstraint.assemble.vector` and :meth:`felupe.MultiPointConstraint.assemble.matrix`.
 
 ..  code-block:: python
 
-    MPC = fe.MultiPointConstraint(
+    MPC = fem.MultiPointConstraint(
         field=field, 
         points=np.arange(mesh.npoints)[mesh.points[:, 0] == 1], 
         centerpoint=mesh.npoints - 1, 
@@ -32,12 +32,8 @@ Finally, add the results of these methods to the internal force vector or the st
 
 ..  code-block:: python
 
-    umat = felupe.constitution.NeoHooke(mu=1.0, bulk=2.0)
+    umat = fem.NeoHooke(mu=1.0, bulk=2.0)
+    body = fem.SolidBody(umat=umat, field=field)
 
-    K = fe.IntegralForm(
-        umat.hessian(field.extract()), field, region.dV, field
-    ).assemble() + MPC.assemble.matrix()
-
-    r = fe.IntegralForm(
-        umat.gradient(field.extract()), field, region.dV
-    ).assemble() + MPC.assemble.vector(field)
\ No newline at end of file
+    K = body.assemble.matrix() + MPC.assemble.matrix()
+    r = body.assemble.vector(field) + MPC.assemble.vector(field)
\ No newline at end of file
diff --git a/docs/howto/solid.rst b/docs/howto/solid.rst
index 57c0e850..4524a0ca 100644
--- a/docs/howto/solid.rst
+++ b/docs/howto/solid.rst
@@ -1,29 +1,32 @@
 Solid Mechanics
 ~~~~~~~~~~~~~~~
 
-The mechanics submodule contains classes for the generation of solid bodies. Solid body objects are supported as items of the Newton-Rhapson procedure.
+The mechanics submodule contains classes for the generation of solid bodies. Solid body objects are supported as items of a Step and the Newton-Rhapson procedure.
 
 Solid Body
 ----------
 
-The generation of internal force vectors or stiffness matrices of solid bodies are provided as assembly-methods of a :class:`felupe.SolidBody`. The correct integral form is chosen based on the :class:`felupe.Field` inside the :class:`felupe.FieldContainer`.
+The generation of internal force vectors and stiffness matrices of solid bodies are provided as assembly-methods of a :class:`felupe.SolidBody` or a :class:`felupe.SolidBodyNearlyIncompressible`.
 
 ..  code-block:: python
 
-    import felupe as fe
+    import felupe as fem
 
-    neohooke = fe.NeoHooke(mu=1.0, bulk=5000.0)
-    mesh = fe.Cube(n=6)
-    region = fe.RegionHexahedron(mesh)
-    displacement = fe.Field(region, dim=3)
-    field = fe.FieldContainer([displacement])
+    mesh = fem.Cube(n=6)
+    region = fem.RegionHexahedron(mesh)
+    field = fem.FieldContainer([fem.Field(region, dim=3)])
+    
+    # a solid body
+    body = fem.SolidBody(umat=fem.NeoHooke(mu=1, bulk=2), field=field)
+    
+    # a (nearly) incompressible solid body (to be used with quads and hexahedrons)
+    body = fem.SolidBodyNearlyIncompressible(umat=fem.NeoHooke(mu=1), field=field, bulk=5000)
     
-    body = fe.SolidBody(umat=neohooke, field=field)
     internal_force = body.assemble.vector(field, parallel=False, jit=False)
     stiffness_matrix = body.assemble.matrix(field, parallel=False, jit=False)
 
 
-During assembly, several results are stored, e.g. the gradient of the strain energy density function per unit undeformed volume (first Piola-Kirchhoff stress tensor). Other results are the deformation gradient or the fourth-order elasticity tensor associated to the first Piola-Kirchhoff stress tensor.
+During assembly, several results are stored, e.g. the gradient of the strain energy density function per unit undeformed volume w.r.t. the deformation gradient (first Piola-Kirchhoff stress tensor). Other results are the deformation gradient or the fourth-order elasticity tensor associated to the first Piola-Kirchhoff stress tensor.
 
 ..  code-block:: python
     
@@ -44,43 +47,41 @@ The Cauchy stress tensor, as well as the gradient and the hessian of the strain
 Body Force (Gravity) on a Solid Body
 ------------------------------------
 
-The generation of internal force vectors or stiffness matrices of body forces acting on solid bodies are provided as assembly-methods of a :class:`felupe.SolidBodyGravity`. The correct integral form is chosen based on the :class:`felupe.Field` inside the :class:`felupe.FieldContainer`. If the internal field is a mixed field, the assembled vectors from the gravity contribution have to be resized to the dimensions of the internal force vector.
+The generation of internal force vectors or stiffness matrices of body forces acting on solid bodies are provided as assembly-methods of a :class:`felupe.SolidBodyGravity`.
 
 
 ..  code-block:: python
     
-    body = fe.SolidBodyGravity(field=field, gravity=[9810, 0, 0], density=7.85e-9)
+    body = fem.SolidBodyGravity(field=field, gravity=[9810, 0, 0], density=7.85e-9)
     
-    internal_force_gravity = body.assemble.vector(
-        field, parallel=False, jit=False, resize=internal_force
-    )
+    force_gravity = body.assemble.vector(field, parallel=False, jit=False)
 
 
 Pressure Boundary on a Solid Body
 ---------------------------------
 
-The generation of internal force vectors or stiffness matrices of pressure boundaries on solid bodies are provided as assembly-methods of a :class:`felupe.SolidBodyPressure`. The correct integral form is chosen based on the :class:`felupe.Field` inside the :class:`felupe.FieldContainer`. If the internal field is a mixed field, the assembled vectors and matrices from the pressure contribution have to be resized to the dimensions of the internal force vector and the stiffness matrix.
+The generation of force vectors or stiffness matrices of pressure boundaries on solid bodies are provided as assembly-methods of a :class:`felupe.SolidBodyPressure`.
 
 ..  code-block:: python
     
-    region_pressure = fe.RegionHexahedronBoundary(
+    region_pressure = fem.RegionHexahedronBoundary(
         mesh=mesh,
         only_surface=True, # select only faces on the outline
         mask=mesh.points[:, 0] == 0, # select a subset of faces on the surface
     )
     
-    displacement_boundary = fe.Field(region_pressure, dim=3)
-    field_boundary = fe.FieldContainer([displacement_boundary])
-    displacement_boundary.values = displacement.values # link field values
+    displacement_boundary = 
+    field_boundary = fem.FieldContainer([fem.Field(region_pressure, dim=3)])
+    field_boundary.link(field)
     
-    body_pressure = fe.SolidBodyPressure(field=field_boundary)
+    body_pressure = fem.SolidBodyPressure(field=field_boundary)
     
-    internal_force_pressure = body_pressure.assemble.vector(
-        field=field_boundary, parallel=False, jit=False, resize=internal_force
+    force_pressure = body_pressure.assemble.vector(
+        field=field_boundary, parallel=False, jit=False
     )
     
     stiffness_matrix_pressure = body_pressure.assemble.matrix(
-        field=field_boundary, parallel=False, jit=False, resize=stiffness_matrix
+        field=field_boundary, parallel=False, jit=False
     )
 
 
@@ -88,16 +89,16 @@ For axisymmetric problems the boundary region has to be created with the attribu
 
 ..  code-block:: python
     
-    mesh = fe.Rectangle(a=(0, 30), b=(20, 40), n=(21, 11))
-    region = fe.RegionQuad(mesh)
+    mesh = fem.Rectangle(a=(0, 30), b=(20, 40), n=(21, 11))
+    region = fem.RegionQuad(mesh)
     
-    region_pressure = fe.RegionQuadBoundary(
+    region_pressure = fem.RegionQuadBoundary(
         mesh=mesh,
         only_surface=True, # select only faces on the outline
         mask=mesh.points[:, 0] == 0, # select a subset of faces on the surface
         ensure_3d=True, # flag for axisymmetric boundary region
     )
     
-    displacement = fe.FieldAxisymmetric(region)
-    displacement_boundary = fe.FieldAxisymmetric(region_pressure)
-    displacement_boundary.values = displacement.values # link field values
\ No newline at end of file
+    field = fem.FieldContainer([fem.FieldAxisymmetric(region)])
+    field_boundary = fem.FieldContainer([fem.FieldAxisymmetric(region_pressure)])
+    field_boundary.link(field)
\ No newline at end of file
diff --git a/docs/howto/solvers.rst b/docs/howto/solvers.rst
index a0a32b84..d74fc5c3 100644
--- a/docs/howto/solvers.rst
+++ b/docs/howto/solvers.rst
@@ -9,18 +9,18 @@ Solvers from SciPy Sparse:
 
    .. code-block:: python
       
-      import felupe as fe
+      import felupe as fem
       
       # ...
       
-      system = fe.solve.partition(field, K, dof1, dof0)
-      fe.solve.solve(*system)
+      system = fem.solve.partition(field, K, dof1, dof0)
+      fem.solve.solve(*system)
 
 .. tab:: SciPy Sparse (direct, symmetric)
 
    .. code-block:: python
       
-      import felupe as fe
+      import felupe as fem
       import scipy.sparse.linalg as spla
       from scipy.sparse import tril
       
@@ -29,8 +29,8 @@ Solvers from SciPy Sparse:
       def solver(A, b):
           return spla.spsolve_triangular(tril(A).tocsr(), b).squeeze()
       
-      system = fe.solve.partition(field, K, dof1, dof0)
-      fe.solve.solve(*system)
+      system = fem.solve.partition(field, K, dof1, dof0)
+      fem.solve.solve(*system)
 
 .. tab:: SciPy Sparse (iterative)
 
@@ -38,7 +38,7 @@ Solvers from SciPy Sparse:
 
    ..  code-block:: python
         
-       import felupe as fe
+       import felupe as fem
        import scipy.sparse.linalg as spla
        
        # ...
@@ -48,8 +48,8 @@ Solvers from SciPy Sparse:
            
            return spla.minres(A, b)[0]
        
-       system = fe.solve.partition(field, K, dof1, dof0)
-       fe.solve.solve(*system, solver=solver)
+       system = fem.solve.partition(field, K, dof1, dof0)
+       fem.solve.solve(*system, solver=solver)
 
 Solvers from external packages:
 
@@ -63,7 +63,7 @@ Solvers from external packages:
 
    ..  code-block:: python
       
-       import felupe as fe
+       import felupe as fem
        from pypardiso import spsolve as solver
        
        # undocumented, untested workaround if multiple blas libaries are installed
@@ -72,8 +72,8 @@ Solvers from external packages:
        
        # ...
        
-       system = fe.solve.partition(field, K, dof1, dof0)
-       fe.solve.solve(*system, solver=solver)
+       system = fem.solve.partition(field, K, dof1, dof0)
+       fem.solve.solve(*system, solver=solver)
 
 .. tab:: PyPardiso (direct, symmetric)
 
@@ -85,7 +85,7 @@ Solvers from external packages:
 
    ..  code-block:: python
       
-       import felupe as fe
+       import felupe as fem
        from pypardiso import PyPardisoSolver
        from scipy.sparse import triu
       
@@ -99,8 +99,8 @@ Solvers from external packages:
            # mtype = 6: complex and symmetric
            return PyPardisoSolver(mtype=2).solve(triu(A).tocsr(), b).squeeze()
       
-       system = fe.solve.partition(field, K, dof1, dof0)
-       fe.solve.solve(*system, solver=solver)
+       system = fem.solve.partition(field, K, dof1, dof0)
+       fem.solve.solve(*system, solver=solver)
 
 
     
\ No newline at end of file
diff --git a/docs/howto/umat.rst b/docs/howto/umat.rst
new file mode 100644
index 00000000..1e048633
--- /dev/null
+++ b/docs/howto/umat.rst
@@ -0,0 +1,46 @@
+Small-Strain based User Materials
+---------------------------------
+
+User materials based on the incremental small-strain tensor, suitable for elastic-plastic material formulations, are to be created with :class:`felupe.UserMaterialStrain`. A user-defined function must be created with the arguments
+
++---------------+-----------------------------+
+| **Arguments** | **Description**             |
++---------------+-----------------------------+
+|      dε       | strain increment            |
++---------------+-----------------------------+
+|      εn       | old strain tensor           |
++---------------+-----------------------------+
+|      σn       | old stress tensor           |
++---------------+-----------------------------+
+|      ζn       | list of old state variables |
++---------------+-----------------------------+
+
+and must return:
+
++-------------+-----------------------------+
+| **Returns** | **Description**             |
++-------------+-----------------------------+
+|     dσdε    | tangent modulus             |
++-------------+-----------------------------+
+|      σ      | new stress tensor           |
++-------------+-----------------------------+
+|      ζ      | list of new state variables |
++-------------+-----------------------------+
+
+..  code-block:: python
+
+    def material(dε, εn, σn, ζn, **kwargs):
+        return dσdε, σ, ζ
+
+This function is further added as the ``material`` argument of :class:`felupe.UserMaterialStrain`. Optionally, the shapes of internal state variables may be passed.
+
+..  code-block:: python
+    
+    import felupe as fem
+    
+    umat = fem.UserMaterialStrain(material=material, statevars=(0,), **kwargs)
+
+FElupe contains two reference user materials, one for linear elastic materials and another one for linear elastic-plastic materials with isotropic hardening:
+
+* :func:`felupe.constitution.linear_elastic`
+* :func:`felupe.constitution.linear_elastic_plastic_isotropic_hardening`
\ No newline at end of file
diff --git a/docs/index.rst b/docs/index.rst
index 4ec463c4..069e84a5 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -74,6 +74,7 @@ Install Python, fire up a terminal and run ``pip install felupe[all]``, where ``
    :caption: Contents:
    
    tutorial
+   examples
    howto
    explanation
 
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index de0e8c4b..e56f254d 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -7,9 +7,5 @@ This section is all about learning. Each tutorial focusses on a some lessons to
    :maxdepth: 1
    :caption: Tutorials:
    
-   tutorial/hellofelupe
-   tutorial/gettingstarted
-   tutorial/platewithhole
-   tutorial/beam
-   tutorial/shear
-   tutorial/numcont
+   tutorial/getting_started
+   tutorial/jobs
diff --git a/docs/tutorial/getting_started.rst b/docs/tutorial/getting_started.rst
new file mode 100644
index 00000000..4be1b360
--- /dev/null
+++ b/docs/tutorial/getting_started.rst
@@ -0,0 +1,62 @@
+.. _tutorial-getting-started:
+
+Getting Started
+---------------
+
+..  admonition:: Your very first steps with FElupe.
+    :class: note
+    
+    * create a meshed cube with hexahedron cells
+    
+    * define a numeric region along with a displacement field
+    
+    * load a Neo-Hookean material formulation
+    
+    * apply a uniaxial loadcase
+    
+    * solve the problem
+    
+    * export the displaced mesh
+
+This tutorial covers the essential high-level parts of creating and solving problems with FElupe. As an introductory example, a quarter model of a solid cube with hyperelastic material behaviour is subjected to a uniaxial elongation applied at a clamped end-face. First, let's import FElupe and create a meshed cube out of hexahedron cells with ``n=11`` points per axis. A numeric region, pre-defined for hexahedrons, is created on the mesh. A vector-valued displacement field is initiated on the region. Next, a field container is created on top of the displacement field.
+
+..  code-block:: python
+
+    import felupe as fem
+    
+    mesh = fem.Cube(n=11)
+    region = fem.RegionHexahedron(mesh=mesh)
+    displacement = fem.Field(region=region, dim=3)
+    field = fem.FieldContainer(fields=[displacement])
+
+A uniaxial load case is applied on the displacement field stored inside the field container. This involves setting up symmetry planes as well as the absolute value of the prescribed displacement at the mesh-points on the right-end face of the cube. The right-end face is *clamped*: only displacements in direction x are allowed. The dict of boundary conditions for this pre-defined load case are returned as ``boundaries`` and the partitioned degrees of freedom as well as the external displacements are stored within the returned dict ``loadcase``.
+
+..  code-block:: python
+
+    boundaries, loadcase = fem.dof.uniaxial(field, move=0.2, right=1, clamped=True)
+
+The material behaviour is defined through a built-in Neo-Hookean material formulation. The constitutive isotropic hyperelastic material formulation is applied on the displacement field by the definition of a solid body.
+
+..  code-block:: python
+
+    umat = fem.NeoHooke(mu=1.0, bulk=2.0)
+    body = fem.SolidBody(umat=umat, field=field)
+
+The problem is solved by an iterative `Newton-Rhapson`_ procedure.
+
+..  code-block:: python
+
+    res = fem.newtonrhapson(items=[body], **loadcase)
+
+Results are saved within a VTK-file.
+    
+..  code-block:: python
+
+    fem.save(region, field, filename="result.vtk")
+
+
+..  image:: images/readme.png
+    :width: 600px
+
+
+..  _Newton-Rhapson: https://en.wikipedia.org/wiki/Newton%27s_method
\ No newline at end of file
diff --git a/docs/tutorial/hellofelupe.rst b/docs/tutorial/hellofelupe.rst
deleted file mode 100644
index c53f7fcf..00000000
--- a/docs/tutorial/hellofelupe.rst
+++ /dev/null
@@ -1,69 +0,0 @@
-.. _tutorial-hello-felupe:
-
-Hello, FElupe!
---------------
-
-..  admonition:: Your very first steps with FElupe.
-    :class: note
-    
-    * create a meshed cube with hexahedron cells
-    
-    * define a numeric region along with a displacement field
-    
-    *  load a Neo-Hookean material formulation
-    
-    * apply a uniaxial loadcase
-    
-    * solve the problem
-    
-    * export the displaced mesh
-
-
-A quarter model of a solid cube with hyperelastic material behavior is subjected to a uniaxial elongation applied at a clamped end-face.
-
-..  image:: images/getting-started.svg
-    :width: 400px
-
-
-First, let's import FElupe and create a meshed cube out of hexahedron cells with ``n=11`` points per axis. A numeric region, pre-defined for hexahedrons, is created on the mesh. A vector-valued displacement field is initiated on the region. Next, a field container is created on top of the displacement field.
-
-..  code-block:: python
-
-    import felupe as fe
-    
-    mesh = fe.Cube(n=11)
-    region = fe.RegionHexahedron(mesh=mesh)
-    displacement = fe.Field(region=region, dim=3)
-    field = fe.FieldContainer(fields=[displacement])
-
-A uniaxial load case is applied on the displacement field stored inside the field container. This involves setting up symmetry planes as well as the absolute value of the prescribed displacement at the mesh-points on the right-end face of the cube. The right-end face is clamped: only displacements in direction x are allowed. The dict of boundary conditions for this pre-defined load case are returned as ``boundaries`` and the partitioned degrees of freedom as well as the external displacements are stored within the returned dict ``loadcase``.
-
-..  code-block:: python
-
-    boundaries, loadcase = fe.dof.uniaxial(field, move=0.2, clamped=True)
-
-The material behavior is defined through a built-in Neo-Hookean material formulation. The  constitutive isotropic hyperelastic material formulation is applied on the displacement field by the definition of a solid body.
-
-..  code-block:: python
-
-    umat = fe.NeoHooke(mu=1.0, bulk=2.0)
-    body = fe.SolidBody(umat=umat, field=field)
-
-Inside a `Newton-Rhapson`_ procedure, the internal force vector and the tangent stiffness matrix are generated by assembling both linear and bilinear forms of static equilibrium. Finally, the solution of the incremental displacements is calculated und updated until convergence is reached.
-
-..  code-block:: python
-
-    res = fe.newtonrhapson(items=[body], **loadcase)
-
-The results are saved within a VTK-file.
-    
-..  code-block:: python
-
-    fe.save(region, res.x, filename="result.vtk")
-
-
-..  image:: images/readme.png
-    :width: 600px
-
-
-..  _Newton-Rhapson: https://en.wikipedia.org/wiki/Newton%27s_method
\ No newline at end of file
diff --git a/docs/_static/readme_xdmf_inc10.png b/docs/tutorial/images/job.png
similarity index 100%
rename from docs/_static/readme_xdmf_inc10.png
rename to docs/tutorial/images/job.png
diff --git a/docs/tutorial/jobs.rst b/docs/tutorial/jobs.rst
new file mode 100644
index 00000000..b64ba711
--- /dev/null
+++ b/docs/tutorial/jobs.rst
@@ -0,0 +1,49 @@
+.. _tutorial-jobs:
+
+Run a Job
+---------
+
+..  admonition:: Learn how to apply boundary conditions in a ramped manner within a **Step** and run a **Job**.
+    :class: note
+    
+    * create a **Step** with ramped boundary conditions
+    
+    * run a **Job** and export a XDMF time-series file
+
+
+This tutorial once again covers the essential high-level parts of creating and solving problems with FElupe. This time, however, the external displacements are applied in a ramped manner. The prescribed displacements of a cube under non-homogenous uniaxial loading will be controlled within a step. The Ogden-Roxburgh pseudo-elastic Mullins softening model is combined with an isotropic hyperelastic Neo-Hookean material formulation, which is further applied on a (nearly) incompressible solid body for a realistic analysis of rubber-like materials. Note that the bulk modulus is now an argument of the (nearly) incompressible solid body instead of the constitutive Neo-Hookean material definition.
+
+..  code-block:: python
+
+    import felupe as fem
+    
+    mesh = fem.Cube(n=11)
+    region = fem.RegionHexahedron(mesh=mesh)
+    field = fem.FieldContainer([fem.Field(region=region, dim=3)])
+    
+    boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
+    
+    umat = fem.OgdenRoxburgh(material=fem.NeoHooke(mu=1), r=3, m=1, beta=0)
+    body = fem.SolidBodyNearlyIncompressible(umat=umat, field=field, bulk=5000)
+
+The ramped prescribed displacements for 20 substeps are created with ``linsteps``. A **Step** is created with a list of items to be considered (here, one single solid body) and a dict of ramped boundary conditions along with the prescribed values.
+
+..  code-block:: python
+
+    move = fem.math.linsteps([0, 2, 0], num=10)
+    uniaxial = fem.Step(
+        items=[solid], 
+        ramp={boundaries["move"]: move}, 
+        boundaries=boundaries
+    )
+
+This step is now added to a **Job**. The results are exported after each completed and successful substep as a time-series XDMF-file.
+
+..  code-block:: python
+
+    job = fem.Job(steps=[uniaxial])
+    job.evaluate(filename="result.xdmf")
+
+
+..  image:: images/job.png
+    :width: 600px
\ No newline at end of file
diff --git a/felupe/constitution/_user_materials_models.py b/felupe/constitution/_user_materials_models.py
index c3f2d1f8..3b968622 100644
--- a/felupe/constitution/_user_materials_models.py
+++ b/felupe/constitution/_user_materials_models.py
@@ -37,12 +37,12 @@
 )
 
 
-def linear_elastic(de, εn, σn, ζn, λ, μ, **kwargs):
+def linear_elastic(dε, εn, σn, ζn, λ, μ, **kwargs):
     """3D linear-elastic material formulation.
 
     Arguments
     ---------
-    de : ndarray
+    dε : ndarray
         Strain increment.
     εn : ndarray
         Old strain tensor.
@@ -58,7 +58,7 @@ def linear_elastic(de, εn, σn, ζn, λ, μ, **kwargs):
 
     # change of stress due to change of strain
     I = identity(dim=3, shape=(1, 1))
-    dσ = 2 * μ * de + λ * trace(de) * I
+    dσ = 2 * μ * dε + λ * trace(dε) * I
 
     # update stress
     σ = σn + dσ
@@ -75,13 +75,13 @@ def linear_elastic(de, εn, σn, ζn, λ, μ, **kwargs):
     return dσdε, σ, ζ
 
 
-def linear_elastic_plastic_isotropic_hardening(de, εn, σn, ζn, λ, μ, σy, K, **kwargs):
+def linear_elastic_plastic_isotropic_hardening(dε, εn, σn, ζn, λ, μ, σy, K, **kwargs):
     """Linear-elastic-plastic material formulation with linear isotropic
     hardening (return mapping algorithm).
 
     Arguments
     ---------
-    de : ndarray
+    dε : ndarray
         Strain increment.
     εn : ndarray
         Old strain tensor.
@@ -99,7 +99,7 @@ def linear_elastic_plastic_isotropic_hardening(de, εn, σn, ζn, λ, μ, σy, K
         Isotropic hardening modulus.
     """
 
-    I = identity(de)
+    I = identity(dε)
 
     # elasticity tensor
     if kwargs["tangent"]:
@@ -108,7 +108,7 @@ def linear_elastic_plastic_isotropic_hardening(de, εn, σn, ζn, λ, μ, σy, K
         dσdε = None
 
     # elastic hypothetic (trial) stress and deviatoric stress
-    dσ = 2 * μ * de + λ * trace(de) * I
+    dσ = 2 * μ * dε + λ * trace(dε) * I
     σ = σn + dσ
     s = σ - 1 / 3 * trace(σ) * I
 
diff --git a/felupe/mechanics/_solidbody_incompressible.py b/felupe/mechanics/_solidbody_incompressible.py
index 1db59726..f22d5864 100644
--- a/felupe/mechanics/_solidbody_incompressible.py
+++ b/felupe/mechanics/_solidbody_incompressible.py
@@ -36,15 +36,15 @@
 
 
 class SolidBodyNearlyIncompressible:
-    """A (nearly) incompressible SolidBody with methods for the assembly of
+    r"""A (nearly) incompressible SolidBody with methods for the assembly of
     sparse vectors/matrices based on a ``MaterialTensor`` with state variables.
-
+    
     The volumetric material behaviour is defined by a strain energy function.
-
+    
     ..  math::
-
+        
         U(J) = \frac{K}{2} (J - 1)^2
-
+    
     """
 
     def __init__(self, umat, field, bulk, state=None, statevars=None):