diff --git a/examples/python/solver_callback/CMakeLists.txt b/examples/python/solver_callback/CMakeLists.txt
new file mode 100644
index 000000000..4f7e8017e
--- /dev/null
+++ b/examples/python/solver_callback/CMakeLists.txt
@@ -0,0 +1,35 @@
+#===============================================================================
+# @file CMakeLists.txt
+#
+# @author Nicolas Richart <nicolas.richart@epfl.ch>
+#
+# @brief CMakeLists for solid mechanics dynamics example
+#
+#
+# @section LICENSE
+#
+# Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+#
+# Akantu is free software: you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License along
+# with Akantu. If not, see <http://www.gnu.org/licenses/>.
+#
+#===============================================================================
+
+
+add_mesh(solver_callback_mesh bar.geo DIM 2)
+register_example(solver_callback
+ SCRIPT solver_callback.py
+ PYTHON
+ FILES_TO_COPY material.dat
+ DEPENDS solver_callback_mesh)
diff --git a/examples/python/solver_callback/bar.geo b/examples/python/solver_callback/bar.geo
new file mode 100644
index 000000000..cf7a819d3
--- /dev/null
+++ b/examples/python/solver_callback/bar.geo
@@ -0,0 +1,35 @@
+// Mesh size
+h = 0.05;
+
+h1 = h;
+h2 = h;
+
+// Dimensions of the bar
+Lx = 10;
+Ly = 1;
+
+// ------------------------------------------
+// Geometry
+// ------------------------------------------
+
+Point(101) = { 0.0, -Ly/2, 0.0, h1};
+Point(102) = { Lx, -Ly/2, 0.0, h2};
+
+Point(103) = { Lx, Ly/2., 0.0, h2};
+Point(104) = { 0.0, Ly/2., 0.0, h1};
+
+Line(201) = {101, 102};
+Line(202) = {102, 103};
+Line(203) = {103, 104};
+Line(204) = {104, 101};
+
+Line Loop(301) = {201, 202, 203, 204};
+Plane Surface(301) = {301};
+
+Transfinite Surface "*";
+Recombine Surface "*";
+Physical Surface("Body") = {301};
+
+Physical Line("YBlocked") = {201, 203};
+Physical Line("Right") = {202};
+Physical Line("Left") = {204};
diff --git a/examples/python/solver_callback/material.dat b/examples/python/solver_callback/material.dat
new file mode 100644
index 000000000..42b9b8ae7
--- /dev/null
+++ b/examples/python/solver_callback/material.dat
@@ -0,0 +1,6 @@
+material elastic [
+ name = fictive
+ rho = 1. # density
+ E = 1. # young modulus
+ nu = .3 # poisson ratio
+]
diff --git a/examples/python/solver_callback/solver_callback.py b/examples/python/solver_callback/solver_callback.py
new file mode 100644
index 000000000..5ef00ca96
--- /dev/null
+++ b/examples/python/solver_callback/solver_callback.py
@@ -0,0 +1,178 @@
+#!/usr/bin/env python3
+""" solver_callback.py: solver_callback overload example"""
+
+__author__ = "Nicolas Richart"
+__credits__ = [
+ "Guillaume Anciaux <guillaume.anciaux@epfl.ch>",
+ "Nicolas Richart <nicolas.richart@epfl.ch>",
+]
+__copyright__ = "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" \
+ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
+ " en Mécanique des Solides)"
+__license__ = "LGPLv3"
+
+import numpy as np
+import akantu as aka
+
+
+class SolverCallback(aka.SolverCallback):
+ def __init__(self, model):
+ super().__init__(model.getDOFManager())
+ self.model = model
+
+ mesh = model.getMesh()
+ left = mesh.getElementGroup("Left").getNodeGroup().getNodes()
+ right = mesh.getElementGroup("Right").getNodeGroup().getNodes()
+ position = mesh.getNodes()
+
+ self.pair = []
+ for node_l in left:
+ node_l = int(node_l)
+ for node_r in right:
+ node_r = int(node_r)
+ if abs(position[node_r, 1] - position[node_l, 1]) < 1e-6:
+ self.pair.append([node_l, node_r])
+
+ blocked_dofs = model.getBlockedDOFs()
+ self.periodic_K_modif = aka.TermsToAssemble()
+
+ matrix_type = self.model.getMatrixType("K")
+ for p in self.pair:
+ #blocked_dofs[p[1]] = True
+ self.periodic_K_modif(p[0]*2, p[0]*2, 1)
+ self.periodic_K_modif(p[0]*2, p[1]*2, -1)
+ if matrix_type == aka._unsymmetric:
+ self.periodic_K_modif(p[1]*2, p[0]*2, -1)
+ self.periodic_K_modif(p[1]*2, p[1]*2, 1)
+
+ self.first = True
+ self.k_release = -1
+
+ def getMatrixType(self, matrix_id):
+ return self.model.getMatrixType(matrix_id)
+
+ def assembleMatrix(self, matrix_id):
+ self.model.assembleMatrix(matrix_id)
+ if matrix_id == "K":
+ release = self.model.getDOFManager().getMatrix("K").getRelease()
+ if release == self.k_release:
+ return
+
+ if self.first:
+ self.model.getDOFManager().getMatrix("K").saveMatrix("K0.mtx")
+
+ self.model.getDOFManager().assemblePreassembledMatrix(
+ 'displacement', 'displacement', "K", self.periodic_K_modif)
+ if self.first:
+ self.model.getDOFManager().getMatrix("K").saveMatrix("K1.mtx")
+
+ self.k_release = self.model.getDOFManager().getMatrix("K").getRelease()
+ self.first = False
+
+ def assembleLumpedMatrix(self, matrix_id):
+ self.model.assembleLumpedMatrix(matrix_id)
+
+ def assembleResidual(self):
+ displacement = self.model.getDisplacement()
+ force = np.zeros(displacement.shape)
+ for p in self.pair:
+ force[p[0], 0] += displacement[p[0], 0] - displacement[p[1], 0]
+ force[p[1], 0] += displacement[p[1], 0] - displacement[p[0], 0]
+ self.model.getDOFManager().assembleToResidual('displacement', force, -1.);
+ self.model.assembleResidual();
+
+# -----------------------------------------------------------------------------
+def main():
+ spatial_dimension = 2
+ mesh_file = 'bar.msh'
+ max_steps = 250
+ time_step = 1e-3
+
+ aka.parseInput('material.dat')
+
+ # -------------------------------------------------------------------------
+ # Initialization
+ # -------------------------------------------------------------------------
+ mesh = aka.Mesh(spatial_dimension)
+ mesh.read(mesh_file)
+
+ model = aka.SolidMechanicsModel(mesh)
+
+ model.initFull(_analysis_method=aka._implicit_dynamic)
+
+ model.setBaseName("solver_callback")
+ model.addDumpFieldVector("displacement")
+ model.addDumpFieldVector("acceleration")
+ model.addDumpFieldVector("velocity")
+ model.addDumpFieldVector("internal_force")
+ model.addDumpFieldVector("external_force")
+ model.addDumpField("strain")
+ model.addDumpField("stress")
+ model.addDumpField("blocked_dofs")
+
+ # -------------------------------------------------------------------------
+ # boundary conditions
+ # -------------------------------------------------------------------------
+ model.applyBC(aka.FixedValue(0, aka._y), "YBlocked")
+
+ # -------------------------------------------------------------------------
+ # initial conditions
+ # -------------------------------------------------------------------------
+ displacement = model.getDisplacement()
+ velocity = model.getVelocity()
+ nb_nodes = mesh.getNbNodes()
+ position = mesh.getNodes()
+
+ L = 1 # pulse_width
+ A = 0.01
+ v = np.sqrt(model.getMaterial(0).getReal('E') /
+ model.getMaterial(0).getReal('rho'))
+ k = 0.1 * 2 * np.pi * 3 / L
+ t = 0.
+ velocity[:, 0] = k * v * A * np.sin(k * ((position[:, 0] - 5.) - v * t))
+ displacement[:, 0] = A * np.cos(k * ((position[:, 0] - 5.) - v * t))
+
+ # -------------------------------------------------------------------------
+ # timestep value computation
+ # -------------------------------------------------------------------------
+ time_factor = 0.8
+ stable_time_step = model.getStableTimeStep() * time_factor
+
+ print("Stable Time Step = {0}".format(stable_time_step))
+ print("Required Time Step = {0}".format(time_step))
+
+ time_step = stable_time_step * time_factor
+
+ model.setTimeStep(time_step)
+ solver_callback = SolverCallback(model)
+
+ solver = model.getNonLinearSolver()
+ solver.set("max_iterations", 100)
+ solver.set("threshold", 1e-7)
+
+ # -------------------------------------------------------------------------
+ # loop for evolution of motion dynamics
+ # -------------------------------------------------------------------------
+ print("step,step * time_step,epot,ekin,epot + ekin")
+ for step in range(0, max_steps + 1):
+
+ model.solveStep(solver_callback)
+ #model.solveStep()
+
+ if step % 10 == 0:
+ model.dump()
+
+ epot = model.getEnergy('potential')
+ ekin = model.getEnergy('kinetic')
+
+ # output energy calculation to screen
+ print("{0},{1},{2},{3},{4}".format(step, step * time_step,
+ epot, ekin,
+ (epot + ekin)))
+
+ return
+
+
+# -----------------------------------------------------------------------------
+if __name__ == "__main__":
+ main()