diff --git a/CHANGELOG.md b/CHANGELOG.md
index a71dde14..7bcf5e41 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -28,6 +28,7 @@ All notable changes to this project will be documented in this file. The format
 - Enforce contiguous arrays for the region shape-function and -gradient arrays `h` and `dhdX`. This recovers the integral-form assembly performance from v8.6.0.
 - Make the private basis classes public (`assembly.expression.Basis`, `assembly.expression.BasisField` and `assembly.expression.BasisArray`) as especially their docstrings are useful to understand how a *Basis* is created on a field.
 - Remove material parameter keyword-arguments in the `function()`-, `gradient()`- and `hessian()`-methods of the core constitutive material formulations. This removes the ability for temporary material parameters in `LinearElastic(E=None, nu=0.3).gradient(x, E=1.0)`, use `LinearElastic(E=1.0, nu=0.3).gradient(x)` instead. This affects the material formulations `LinearElastic`, `LinearElasticPlaneStrain`, `LinearElasticPlaneStress`, `LinearElasticLargeStrain`, `NeoHooke` and `NeoHookeCompressible`. `None` is still supported for the material parameters of `NeoHooke` and `NeoHookeCompressible`.
+- Remove `LinearElasticPlaneStrain` from the top-level package namespace because this should only be used with `Field(region, dim=2)`. The preferred method is to use `FieldPlaneStrain(region, dim=2)` and the default `LinearElastic`. `LinearElasticPlaneStrain` remains available in `constitution.LinearElasticPlaneStrain`.
 
 ### Deprecated
 - Deprecate `SolidBodyGravity`, `SolidBodyForce` should be used instead.
diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py
index 388705a0..4c679413 100644
--- a/src/felupe/__init__.py
+++ b/src/felupe/__init__.py
@@ -24,7 +24,6 @@
     LinearElastic,
     LinearElasticLargeStrain,
     LinearElasticOrthotropic,
-    LinearElasticPlaneStrain,
     LinearElasticPlaneStress,
     LinearElasticPlasticIsotropicHardening,
     LineChange,
@@ -174,7 +173,6 @@
     "LinearElastic",
     "LinearElasticLargeStrain",
     "LinearElasticOrthotropic",
-    "LinearElasticPlaneStrain",
     "LinearElasticPlaneStress",
     "LinearElasticPlasticIsotropicHardening",
     "LineChange",
diff --git a/src/felupe/constitution/linear_elasticity/_linear_elastic.py b/src/felupe/constitution/linear_elasticity/_linear_elastic.py
index 2fcfbab2..a671d75b 100644
--- a/src/felupe/constitution/linear_elasticity/_linear_elastic.py
+++ b/src/felupe/constitution/linear_elasticity/_linear_elastic.py
@@ -336,6 +336,13 @@ class LinearElasticPlaneStrain:
     nu : float
         Poisson ratio.
 
+    Notes
+    -----
+    .. warning::
+
+       This class must not be used with :class:`~felupe.FieldPlaneStrain` but with
+       :class:`~felupe.Field` instead!
+
     """
 
     def __init__(self, E, nu):
diff --git a/src/felupe/mechanics/_solidbody_cauchy_stress.py b/src/felupe/mechanics/_solidbody_cauchy_stress.py
index 06b0dcde..f754a26c 100644
--- a/src/felupe/mechanics/_solidbody_cauchy_stress.py
+++ b/src/felupe/mechanics/_solidbody_cauchy_stress.py
@@ -48,7 +48,9 @@ class SolidBodyCauchyStress:
     Examples
     --------
     ..  pyvista-plot::
+        :force_static:
 
+        >>> import numpy as np
         >>> import felupe as fem
         >>>
         >>> mesh = fem.Rectangle(n=6)
diff --git a/tests/test_planestrain.py b/tests/test_planestrain.py
index 576108bf..fca188a0 100644
--- a/tests/test_planestrain.py
+++ b/tests/test_planestrain.py
@@ -75,7 +75,7 @@ def test_planestrain():
     g = fem.FieldsMixed(r, n=1, planestrain=False)
 
     u = fem.LinearElastic(E=1, nu=0.3)
-    v = fem.LinearElasticPlaneStrain(E=1, nu=0.3)
+    v = fem.constitution.LinearElasticPlaneStrain(E=1, nu=0.3)
 
     b = fem.SolidBody(u, f)
     c = fem.SolidBody(v, g)
diff --git a/tests/test_tools.py b/tests/test_tools.py
index 6ae1cb8b..ff49a51f 100644
--- a/tests/test_tools.py
+++ b/tests/test_tools.py
@@ -225,7 +225,7 @@ def test_newton_plane():
     )
 
     # define the constitutive material behavior
-    umat = fem.LinearElasticPlaneStrain(E=1.0, nu=0.3)
+    umat = fem.constitution.LinearElasticPlaneStrain(E=1.0, nu=0.3)
 
     # newton-rhapson procedure
     res = fem.newtonrhapson(