diff --git a/examples/visco_elasticity/spring_kelvin_model.py b/examples/visco_elasticity/spring_kelvin_model.py
new file mode 100644
index 0000000..b02dcb3
--- /dev/null
+++ b/examples/visco_elasticity/spring_kelvin_model.py
@@ -0,0 +1,98 @@
+from __future__ import annotations
+
+import numpy as np
+
+from fenics_constitutive import Constraint, IncrSmallStrainModel, strain_from_grad_u
+
+
+class SpringKelvinModel(IncrSmallStrainModel):
+    ''' viscoelastic model based on 1D Three Parameter Model with spring and Kelvin body in row
+
+                               |--- E_1: spring          ---|
+    --- E_0: spring modulus ---|                            |
+                               |--- tau: relaxation time ---|
+
+    with deviatoric assumptions for 3D generalization (volumetric part of visco strain == 0 damper just working on deviatoric part)
+    time integration: backward Euler
+
+    '''
+    def __init__(self, parameters: dict[str, float], constraint: Constraint):
+        self._constraint = constraint
+        self.E0 = parameters["E0"]
+        self.E1 = parameters["E1"]
+        self.tau = parameters["tau"]
+        self.nu = parameters["nu"]
+        self.mu0 = self.E0 / (2.0 * (1.0 + self.nu))
+        self.lam0 = self.E0 * self.nu / ((1.0 + self.nu) * (1.0 - 2.0 * self.nu))
+        self.mu1 = self.E1 / (2.0 * (1.0 + self.nu))
+        self.lam1 = self.E1 * self.nu / ((1.0 + self.nu) * (1.0 - 2.0 * self.nu))
+
+        match constraint:
+            case Constraint.FULL:
+                self.D = None
+
+            case Constraint.PLANE_STRAIN:
+                self.D = None
+
+            case Constraint.PLANE_STRESS:
+                self.D = None
+
+            case Constraint.UNIAXIAL_STRESS:
+                self.D_E0 = np.array([[self.E0]])
+            case _:
+                msg = "Constraint not implemented"
+                raise NotImplementedError(msg)
+
+    def evaluate(
+        self,
+        del_t: float,
+        grad_del_u: np.ndarray,
+        mandel_stress: np.ndarray,
+        tangent: np.ndarray,
+        history: np.ndarray | dict[str, np.ndarray] | None,
+    ) -> None:
+        assert (
+            grad_del_u.size // (self.geometric_dim**2)
+            == mandel_stress.size // self.stress_strain_dim
+            == tangent.size // (self.stress_strain_dim**2)
+        )
+        # n_gauss = grad_del_u.size // (self.geometric_dim**2)
+        mandel_view = mandel_stress.reshape(-1, self.stress_strain_dim)
+        strain_increment = strain_from_grad_u(grad_del_u, self.constraint).reshape(-1, self.stress_strain_dim)
+        strain_visco_n = history['strain_visko'].reshape(-1, self.stress_strain_dim)
+
+        # loop over gauss points
+        for n, eps in enumerate(strain_increment):
+
+            print('eps', eps)
+            print('n', n)
+            stress_n = mandel_view[n]
+            strain_visco_n = strain_visco_n[n]
+
+            # compute visko strain 1D case only !!!
+            factor = (1 / del_t + 1 / self.tau + self.E_0 / (self.tau * self.E_1))
+            deps_visko = (stress_n / (self.tau * self.E_1) + self.E_0 / (self.tau * self.E_1) * eps -
+                          strain_visco_n / self.tau) / factor
+            dstress = self.E_0 * (eps - deps_visko)
+            D = self.D_E0 * (1 - self.E0**2/(self.E1*self.tau *factor))
+
+            print('D', D)
+
+            # update values
+            mandel_view[n] += dstress
+            history['strain_visco'][n] += deps_visko
+            history['strain'][n] += eps
+            tangent[n] = D.flatten()
+
+
+    @property
+    def constraint(self) -> Constraint:
+        return self._constraint
+
+    @property
+    def history_dim(self) -> None:
+        return {'strain_visco': self.stress_strain_dim, 'strain': self.stress_strain_dim}
+
+    def update(self) -> None:
+        pass
+
diff --git a/examples/visco_elasticity/test_viscoelasticity.py b/examples/visco_elasticity/test_viscoelasticity.py
new file mode 100644
index 0000000..1e69229
--- /dev/null
+++ b/examples/visco_elasticity/test_viscoelasticity.py
@@ -0,0 +1,71 @@
+from __future__ import annotations
+
+import dolfinx as df
+import numpy as np
+import pytest
+import ufl
+from dolfinx.nls.petsc import NewtonSolver
+from spring_kelvin_model import SpringKelvinModel
+from linear_elasticity.linear_elasticity_model import LinearElasticityModel
+from mpi4py import MPI
+
+from fenics_constitutive import Constraint, IncrSmallStrainProblem
+
+youngs_modulus = 42.0
+poissons_ratio = 0.3
+visco_modulus = 10.0
+relaxation_time = 10.0
+
+
+def test_creep_uniaxial_stress():
+    mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2)
+    V = df.fem.FunctionSpace(mesh, ("CG", 1))
+    u = df.fem.Function(V)
+    law = SpringKelvinModel(
+        parameters={"E0": youngs_modulus, "E1": visco_modulus, "tau": relaxation_time, "nu": poissons_ratio},
+        constraint=Constraint.UNIAXIAL_STRESS,
+    )
+    # law = LinearElasticityModel(
+    #     parameters={"E": youngs_modulus, "nu": poissons_ratio},
+    #     constraint=Constraint.UNIAXIAL_STRESS,
+    # )
+
+    def left_boundary(x):
+        return np.isclose(x[0], 0.0)
+
+    def right_boundary(x):
+        return np.isclose(x[0], 1.0)
+
+    displacement = df.fem.Constant(mesh, 0.01)
+    dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary)
+    dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary)
+    bc_left = df.fem.dirichletbc(df.fem.Constant(mesh, 0.0), dofs_left, V)
+    bc_right = df.fem.dirichletbc(displacement, dofs_right, V)
+
+    problem = IncrSmallStrainProblem(
+        law,
+        u,
+        [bc_left, bc_right],
+    )
+
+    solver = NewtonSolver(MPI.COMM_WORLD, problem)
+    n, converged = solver.solve(u)
+    # assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (
+    #     youngs_modulus * 0.01
+    # )
+    #
+    # problem.update()
+    # assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (
+    #     youngs_modulus * 0.01
+    # )
+    # assert np.max(problem._u0.x.array) == displacement.value
+    #
+    # displacement.value = 0.02
+    # n, converged = solver.solve(u)
+    # assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / (
+    #     youngs_modulus * 0.02
+    # )
+
+
+if __name__ == "__main__":
+    test_creep_uniaxial_stress()
\ No newline at end of file