diff --git a/.gitignore b/.gitignore
index 9eddc43d82aa89c9d1c0a3b3f0adef19ce0dd100..1fdf329b1c2dcfac1f58a64d11fdecaeb9cb5c5e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -85,6 +85,9 @@ target/
 
 # svn
 .svn
+.cache/
+.pytest_cache/
+
 # coverage
 .coverage
 
@@ -102,4 +105,6 @@ doc/*.ipynb
 *.h5
 *.xdmf
 doc/notebook.rst
+doc/README.rst
+doc/tutorials.rst
 *jitfailure-dolfin*/
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 41c18e7ff07ee1fe04812488f0910a98d6220de2..da4a78247391c57861ef9cf2655fd109aee73164 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -36,8 +36,7 @@ coverage:
   - conda env create -f conda/env.yaml
   - source activate bvpy
   - conda install conda-bld/linux-64/bvpy-*
-  - nosetests -v -x
-  - coverage report
+  - pytest --cov=bvpy -v -x
   retry:
     max: 2
     when: runner_system_failure
@@ -109,7 +108,6 @@ pages:
       - if: '$CI_COMMIT_BRANCH == "master"'
         when: on_success
 
-
 docker_build_deploy:
   image: docker:latest
   stage: deploy
diff --git a/.pkglts/pkg_cfg.json b/.pkglts/pkg_cfg.json
index df471d780cb52b996ae74ef7196a2d1280f8f02d..bfeac7c26925f7c511198dd6dea09d069a3fbbbb 100644
--- a/.pkglts/pkg_cfg.json
+++ b/.pkglts/pkg_cfg.json
@@ -47,11 +47,11 @@
     "namespace_method": "pkg_util"
   },
   "test": {
-    "suite_name": "nose"
+    "suite_name": "pytest"
   },
   "version": {
     "major": 1,
     "minor": 0,
-    "post": 0
+    "post": 1
   }
 }
diff --git a/.pkglts/pkg_hash.json b/.pkglts/pkg_hash.json
index 6db8e484ecf0e3cacec949584e1675c4c963d95e..0ed8526eeec39070a47502edf15ca8fc3edc4f90 100644
--- a/.pkglts/pkg_hash.json
+++ b/.pkglts/pkg_hash.json
@@ -3,11 +3,11 @@
     "coverage": "Wdw0A46bthKCXu+bEDTYhUqOiC6oO1z6vr7BJzujMEDHrqYW6nwpebmtffbKEeQtJ+CllE6GyQ/S+E7WsGRN7g=="
   },
   ".gitignore": {
-    "git": "PQc3hvWK2Bkdaz48B41hV3wb8vasE0g/tt4SgjNxfsaUQQhe7V17qjwtDkE6q9ySjgQDSww10SmmibMnn8gT/A=="
+    "git": "77pwFZlnA/vkpAcyacCHkz6TTwwEUrNQspX9U4c8JSQ39A5NazBHfQFpmcWz/3L46HRH+X+RJ+4jFW9RXvs4MA=="
   },
   "AUTHORS.rst": {
     "doc.authors": "8m7WJgCsHuCrAK8j/orSy/0v36/TEHawg8zn8Fu98kw0G5doIGOSaOLizrkp3e/TCOaxTeCqFHaFmHMByRfgFA==",
-    "doc.contributors": "4ogNKmPWl19ent0JBpc2GefSpdxSyIp17Yq2XIzy5qF7fu6/rkRWPH+NvC1qURedbsgGOO10huCWl+opq16MfA=="
+    "doc.contributors": "pi2IMJvG9beHzzq+o4FlhAA4ncMolpvZrv9pDljrVTdaiwfXeiTySLkG7Rr1N57HcZAb5lJ9dRe1w5l+qOXtTg=="
   },
   "CONTRIBUTING.rst": {},
   "HISTORY.rst": {},
@@ -28,20 +28,20 @@
     "sphinx": "E9sFsnXAm9uPYBn7f0Qb6edct0UKizsXRIJyjZISRCjHbT2HldXK4GIj1yO69R/OVSfnmdAX5lzXQ4wONQNFgA=="
   },
   "requirements.txt": {
-    "reqs": "UY5N7ID42/r9fOX60LaDbFpZbuuklWjYPbn0oq7n3f2WOZ5MnrYZmCatYavLkr8alDbpcD5nN12JBNa52txGzA=="
+    "reqs": "T02TvsmYi9qrLOKFmlKzraUSIwfedls+Z0iG0JzFrz4SlI/wb3QrTsPofbdJV2KE2GA/Q3/j5YeZaVHnC5nN6g=="
   },
   "requirements_minimal.txt": {
     "reqs": "AWuoxM/eZa+Zy1+ouKN+Lrc/SBs640mRZm3y4E/rbAOGZuvR7CtvYjlndWAzxwLd5fQj99R6tu0YJ/9TeDcx9w=="
   },
   "setup.cfg": {
     "sphinx": "mYLeHXWNqXwTk9eWxY4wTrltnDTrAHQvl0l+kosuIfAgrbAuhWMgySq4WR2Y855Q5HDAtGZY7Pnt+BAU6niXIg==",
-    "test.nose": "8BEO3bAAzFeAdtXqubagljVsKlDTTlPhxqeM7FqH6eWNEjAPsfMimlW5LBo35uzE2WbJh7gL4SpkKuNlMU/LsA==",
-    "test.pytest": "vmiIOMqGhuXJBom/KrWFzvETfJmbSMcLkvZ6XDTcFWl7XRHJgu1tcb4eHn97TgcziEqpfD96M5qO0DV3z3S+CQ==",
+    "test.nose": "vmiIOMqGhuXJBom/KrWFzvETfJmbSMcLkvZ6XDTcFWl7XRHJgu1tcb4eHn97TgcziEqpfD96M5qO0DV3z3S+CQ==",
+    "test.pytest": "R5P5XMaQxSXsNQJ3hXYY45iLvnl8Mbec6XuIaFTum52gym+wTM2TfEbW/N2doHNaCkpdAWg07oblbYi15AwQkA==",
     "wheel": "vUJZkCvbK4LQuEURneKnX/SNIwhJ9/9EaT5T6Uw1K1fauKFFiHnq1DTDES7cgBD0OANelyvwCCZO/rV4Iok+ag=="
   },
   "setup.py": {
     "pysetup.call": "t30iIuCMd/sWxStR2XNYb8bCBgwanhKt8YYn2VKACyAqL5cKFRUWpVRIpAR+uxXUZ9PYYy2csMmECcONnRZE8w==",
-    "pysetup.kwds": "f4aACD9E05zMLYsGLGCmuXYWoJYRxvBWIPnYpSODKQihBS2+4TDmMjPDlKAUsZrUN88JKQhdIOxsqmZMdV0+OQ=="
+    "pysetup.kwds": "cAzFCKINVCtdrW5WhuidyviQrjIRX8YUv/rk3SF+rIsw8qmZw/sHOk10U61j0JP5yPh3jTQcLE3CRoC59uJl7A=="
   },
   "src/bvpy/__init__.py": {
     "src": "L2dZgUBRC5qv3o09+1CbMcgYfQtR1ZNSvP8KlP8REQazZiAQAktALjpCC2ijdxvKPv/tN9KsZrzNW+lZMszSfg==",
diff --git a/.pkglts/pkg_version.json b/.pkglts/pkg_version.json
index 27f630add8f9d7353374a750ee170bee20753094..0e6ee391060c18a526902aaadd1ed3c36073ebeb 100644
--- a/.pkglts/pkg_version.json
+++ b/.pkglts/pkg_version.json
@@ -1,13 +1,13 @@
 {
-  "base": "5.1.9",
-  "coverage": "5.1.9",
-  "doc": "5.1.9",
-  "git": "5.1.9",
-  "license": "5.1.9",
-  "pysetup": "5.1.9",
-  "reqs": "5.1.9",
-  "sphinx": "5.1.9",
-  "src": "5.1.9",
-  "test": "5.1.9",
-  "version": "5.1.9"
+  "base": "5.2.0",
+  "coverage": "5.2.0",
+  "doc": "5.2.0",
+  "git": "5.2.0",
+  "license": "5.2.0",
+  "pysetup": "5.2.0",
+  "reqs": "5.2.0",
+  "sphinx": "5.2.0",
+  "src": "5.2.0",
+  "test": "5.2.0",
+  "version": "5.2.0"
 }
\ No newline at end of file
diff --git a/AUTHORS.rst b/AUTHORS.rst
index f84e0d75a081691dfa67fc45ae96d5229315033e..55cf4cab5901be7f642a44f09d8c3163f7eeab6f 100644
--- a/AUTHORS.rst
+++ b/AUTHORS.rst
@@ -16,11 +16,6 @@ Contributors
 
 .. {# pkglts, doc.contributors
 
-* Florian <florian.gacon@inria.fr>
-* Olivier Ali <olivier.ali@inria.fr>
-* Florian <florian.gacon@gmail.com>
-* GACON Florian <florian.gacon@inria.fr>
-* Christophe Godin <christophe.godin@inria.fr>
 * ALI Olivier <olivier.ali@inria.fr>
 
 .. #}
diff --git a/LICENSE b/LICENSE
index a2d60b55ffa79ae0d2c73f846855aea57a8a279f..650a3f2c7442de4804864065ca70f3d0b8951a2e 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,520 +1,158 @@
 {# pkglts, license
+GNU LESSER GENERAL PUBLIC LICENSE
+Version 3, 29 June 2007
 
-CeCILL-C FREE SOFTWARE LICENSE AGREEMENT
+Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
 
+Everyone is permitted to copy and distribute verbatim copies of this license
+document, but changing it is not allowed.
 
-    Notice
+This version of the GNU Lesser General Public License incorporates the terms
+and conditions of version 3 of the GNU General Public License, supplemented
+by the additional permissions listed below.
 
-This Agreement is a Free Software license agreement that is the result
-of discussions between its authors in order to ensure compliance with
-the two main principles guiding its drafting:
+0. Additional Definitions.
 
-    * firstly, compliance with the principles governing the distribution
-      of Free Software: access to source code, broad rights granted to
-      users,
-    * secondly, the election of a governing law, French law, with which
-      it is conformant, both as regards the law of torts and
-      intellectual property law, and the protection that it offers to
-      both authors and holders of the economic rights over software.
+As used herein, “this License” refers to version 3 of the GNU Lesser General
+Public License, and the “GNU GPL” refers to version 3 of the
+GNU General Public License.
 
-The authors of the CeCILL-C (for Ce[a] C[nrs] I[nria] L[ogiciel] L[ibre])
-license are:
+“The Library” refers to a covered work governed by this License, other than
+an Application or a Combined Work as defined below.
 
-Commissariat a l'Energie Atomique - CEA, a public scientific, technical
-and industrial research establishment, having its principal place of
-business at 25 rue Leblanc, immeuble Le Ponant D, 75015 Paris, France.
+An “Application” is any work that makes use of an interface provided by the
+Library, but which is not otherwise based on the Library. Defining a subclass
+of a class defined by the Library is deemed a mode of using an interface
+provided by the Library.
 
-Centre National de la Recherche Scientifique - CNRS, a public scientific
-and technological establishment, having its principal place of business
-at 3 rue Michel-Ange, 75794 Paris cedex 16, France.
+A “Combined Work” is a work produced by combining or linking an Application
+with the Library. The particular version of the Library with which the
+Combined Work was made is also called the “Linked Version”.
 
-Institut National de Recherche en Informatique et en Automatique -
-INRIA, a public scientific and technological establishment, having its
-principal place of business at Domaine de Voluceau, Rocquencourt, BP
-105, 78153 Le Chesnay cedex, France.
+The “Minimal Corresponding Source” for a Combined Work means the Corresponding
+Source for the Combined Work, excluding any source code for portions of the
+Combined Work that, considered in isolation, are based on the Application,
+and not on the Linked Version.
 
+The “Corresponding Application Code” for a Combined Work means the object code
+and/or source code for the Application, including any data and utility programs
+needed for reproducing the Combined Work from the Application, but excluding
+the System Libraries of the Combined Work.
 
-    Preamble
+1. Exception to Section 3 of the GNU GPL.
 
-The purpose of this Free Software license agreement is to grant users
-the right to modify and re-use the software governed by this license.
+You may convey a covered work under sections 3 and 4 of this License without
+being bound by section 3 of the GNU GPL.
 
-The exercising of this right is conditional upon the obligation to make
-available to the community the modifications made to the source code of
-the software so as to contribute to its evolution.
+2. Conveying Modified Versions.
 
-In consideration of access to the source code and the rights to copy,
-modify and redistribute granted by the license, users are provided only
-with a limited warranty and the software's author, the holder of the
-economic rights, and the successive licensors only have limited liability.
-
-In this respect, the risks associated with loading, using, modifying
-and/or developing or reproducing the software by the user are brought to
-the user's attention, given its Free Software status, which may make it
-complicated to use, with the result that its use is reserved for
-developers and experienced professionals having in-depth computer
-knowledge. Users are therefore encouraged to load and test the
-suitability of the software as regards their requirements in conditions
-enabling the security of their systems and/or data to be ensured and,
-more generally, to use and operate it in the same conditions of
-security. This Agreement may be freely reproduced and published,
-provided it is not altered, and that no provisions are either added or
-removed herefrom.
-
-This Agreement may apply to any or all software for which the holder of
-the economic rights decides to submit the use thereof to its provisions.
-
-
-    Article 1 - DEFINITIONS
-
-For the purpose of this Agreement, when the following expressions
-commence with a capital letter, they shall have the following meaning:
-
-Agreement: means this license agreement, and its possible subsequent
-versions and annexes.
-
-Software: means the software in its Object Code and/or Source Code form
-and, where applicable, its documentation, "as is" when the Licensee
-accepts the Agreement.
-
-Initial Software: means the Software in its Source Code and possibly its
-Object Code form and, where applicable, its documentation, "as is" when
-it is first distributed under the terms and conditions of the Agreement.
-
-Modified Software: means the Software modified by at least one
-Integrated Contribution.
-
-Source Code: means all the Software's instructions and program lines to
-which access is required so as to modify the Software.
-
-Object Code: means the binary files originating from the compilation of
-the Source Code.
-
-Holder: means the holder(s) of the economic rights over the Initial
-Software.
-
-Licensee: means the Software user(s) having accepted the Agreement.
-
-Contributor: means a Licensee having made at least one Integrated
-Contribution.
-
-Licensor: means the Holder, or any other individual or legal entity, who
-distributes the Software under the Agreement.
-
-Integrated Contribution: means any or all modifications, corrections,
-translations, adaptations and/or new functions integrated into the
-Source Code by any or all Contributors.
-
-Related Module: means a set of sources files including their
-documentation that, without modification to the Source Code, enables
-supplementary functions or services in addition to those offered by the
-Software.
-
-Derivative Software: means any combination of the Software, modified or
-not, and of a Related Module.
-
-Parties: mean both the Licensee and the Licensor.
-
-These expressions may be used both in singular and plural form.
-
-
-    Article 2 - PURPOSE
-
-The purpose of the Agreement is the grant by the Licensor to the
-Licensee of a non-exclusive, transferable and worldwide license for the
-Software as set forth in Article 5 hereinafter for the whole term of the
-protection granted by the rights over said Software. 
-
-
-    Article 3 - ACCEPTANCE
-
-3.1 The Licensee shall be deemed as having accepted the terms and
-conditions of this Agreement upon the occurrence of the first of the
-following events:
-
-    * (i) loading the Software by any or all means, notably, by
-      downloading from a remote server, or by loading from a physical
-      medium;
-    * (ii) the first time the Licensee exercises any of the rights
-      granted hereunder.
-
-3.2 One copy of the Agreement, containing a notice relating to the
-characteristics of the Software, to the limited warranty, and to the
-fact that its use is restricted to experienced users has been provided
-to the Licensee prior to its acceptance as set forth in Article 3.1
-hereinabove, and the Licensee hereby acknowledges that it has read and
-understood it.
-
-
-    Article 4 - EFFECTIVE DATE AND TERM
-
-
-      4.1 EFFECTIVE DATE
-
-The Agreement shall become effective on the date when it is accepted by
-the Licensee as set forth in Article 3.1.
-
-
-      4.2 TERM
-
-The Agreement shall remain in force for the entire legal term of
-protection of the economic rights over the Software.
-
-
-    Article 5 - SCOPE OF RIGHTS GRANTED
-
-The Licensor hereby grants to the Licensee, who accepts, the following
-rights over the Software for any or all use, and for the term of the
-Agreement, on the basis of the terms and conditions set forth hereinafter.
-
-Besides, if the Licensor owns or comes to own one or more patents
-protecting all or part of the functions of the Software or of its
-components, the Licensor undertakes not to enforce the rights granted by
-these patents against successive Licensees using, exploiting or
-modifying the Software. If these patents are transferred, the Licensor
-undertakes to have the transferees subscribe to the obligations set
-forth in this paragraph.
-
-
-      5.1 RIGHT OF USE
-
-The Licensee is authorized to use the Software, without any limitation
-as to its fields of application, with it being hereinafter specified
-that this comprises:
-
-   1. permanent or temporary reproduction of all or part of the Software
-      by any or all means and in any or all form.
-
-   2. loading, displaying, running, or storing the Software on any or
-      all medium.
-
-   3. entitlement to observe, study or test its operation so as to
-      determine the ideas and principles behind any or all constituent
-      elements of said Software. This shall apply when the Licensee
-      carries out any or all loading, displaying, running, transmission
-      or storage operation as regards the Software, that it is entitled
-      to carry out hereunder.
-
-
-      5.2 RIGHT OF MODIFICATION
-
-The right of modification includes the right to translate, adapt,
-arrange, or make any or all modifications to the Software, and the right
-to reproduce the resulting software. It includes, in particular, the
-right to create a Derivative Software.
-
-The Licensee is authorized to make any or all modification to the
-Software provided that it includes an explicit notice that it is the
-author of said modification and indicates the date of the creation thereof.
-
-
-      5.3 RIGHT OF DISTRIBUTION
-
-In particular, the right of distribution includes the right to publish,
-transmit and communicate the Software to the general public on any or
-all medium, and by any or all means, and the right to market, either in
-consideration of a fee, or free of charge, one or more copies of the
-Software by any means.
-
-The Licensee is further authorized to distribute copies of the modified
-or unmodified Software to third parties according to the terms and
-conditions set forth hereinafter.
-
-
-        5.3.1 DISTRIBUTION OF SOFTWARE WITHOUT MODIFICATION
-
-The Licensee is authorized to distribute true copies of the Software in
-Source Code or Object Code form, provided that said distribution
-complies with all the provisions of the Agreement and is accompanied by:
-
-   1. a copy of the Agreement,
-
-   2. a notice relating to the limitation of both the Licensor's
-      warranty and liability as set forth in Articles 8 and 9,
-
-and that, in the event that only the Object Code of the Software is
-redistributed, the Licensee allows effective access to the full Source
-Code of the Software at a minimum during the entire period of its
-distribution of the Software, it being understood that the additional
-cost of acquiring the Source Code shall not exceed the cost of
-transferring the data.
-
-
-        5.3.2 DISTRIBUTION OF MODIFIED SOFTWARE
-
-When the Licensee makes an Integrated Contribution to the Software, the
-terms and conditions for the distribution of the resulting Modified
-Software become subject to all the provisions of this Agreement.
-
-The Licensee is authorized to distribute the Modified Software, in
-source code or object code form, provided that said distribution
-complies with all the provisions of the Agreement and is accompanied by:
-
-   1. a copy of the Agreement,
-
-   2. a notice relating to the limitation of both the Licensor's
-      warranty and liability as set forth in Articles 8 and 9,
-
-and that, in the event that only the object code of the Modified
-Software is redistributed, the Licensee allows effective access to the
-full source code of the Modified Software at a minimum during the entire
-period of its distribution of the Modified Software, it being understood
-that the additional cost of acquiring the source code shall not exceed
-the cost of transferring the data.
-
-
-        5.3.3 DISTRIBUTION OF DERIVATIVE SOFTWARE
-
-When the Licensee creates Derivative Software, this Derivative Software
-may be distributed under a license agreement other than this Agreement,
-subject to compliance with the requirement to include a notice
-concerning the rights over the Software as defined in Article 6.4.
-In the event the creation of the Derivative Software required modification 
-of the Source Code, the Licensee undertakes that:
-
-   1. the resulting Modified Software will be governed by this Agreement,
-   2. the Integrated Contributions in the resulting Modified Software
-      will be clearly identified and documented,
-   3. the Licensee will allow effective access to the source code of the
-      Modified Software, at a minimum during the entire period of
-      distribution of the Derivative Software, such that such
-      modifications may be carried over in a subsequent version of the
-      Software; it being understood that the additional cost of
-      purchasing the source code of the Modified Software shall not
-      exceed the cost of transferring the data.
-
-
-        5.3.4 COMPATIBILITY WITH THE CeCILL LICENSE
-
-When a Modified Software contains an Integrated Contribution subject to
-the CeCILL license agreement, or when a Derivative Software contains a
-Related Module subject to the CeCILL license agreement, the provisions
-set forth in the third item of Article 6.4 are optional.
-
-
-    Article 6 - INTELLECTUAL PROPERTY
-
-
-      6.1 OVER THE INITIAL SOFTWARE
-
-The Holder owns the economic rights over the Initial Software. Any or
-all use of the Initial Software is subject to compliance with the terms
-and conditions under which the Holder has elected to distribute its work
-and no one shall be entitled to modify the terms and conditions for the
-distribution of said Initial Software.
-
-The Holder undertakes that the Initial Software will remain ruled at
-least by this Agreement, for the duration set forth in Article 4.2.
-
-
-      6.2 OVER THE INTEGRATED CONTRIBUTIONS
-
-The Licensee who develops an Integrated Contribution is the owner of the
-intellectual property rights over this Contribution as defined by
-applicable law.
-
-
-      6.3 OVER THE RELATED MODULES
-
-The Licensee who develops a Related Module is the owner of the
-intellectual property rights over this Related Module as defined by
-applicable law and is free to choose the type of agreement that shall
-govern its distribution under the conditions defined in Article 5.3.3.
-
-
-      6.4 NOTICE OF RIGHTS
-
-The Licensee expressly undertakes:
-
-   1. not to remove, or modify, in any manner, the intellectual property
-      notices attached to the Software;
-
-   2. to reproduce said notices, in an identical manner, in the copies
-      of the Software modified or not;
-
-   3. to ensure that use of the Software, its intellectual property
-      notices and the fact that it is governed by the Agreement is
-      indicated in a text that is easily accessible, specifically from
-      the interface of any Derivative Software.
-
-The Licensee undertakes not to directly or indirectly infringe the
-intellectual property rights of the Holder and/or Contributors on the
-Software and to take, where applicable, vis-a-vis its staff, any and all
-measures required to ensure respect of said intellectual property rights
-of the Holder and/or Contributors.
-
-
-    Article 7 - RELATED SERVICES
-
-7.1 Under no circumstances shall the Agreement oblige the Licensor to
-provide technical assistance or maintenance services for the Software.
-
-However, the Licensor is entitled to offer this type of services. The
-terms and conditions of such technical assistance, and/or such
-maintenance, shall be set forth in a separate instrument. Only the
-Licensor offering said maintenance and/or technical assistance services
-shall incur liability therefor.
-
-7.2 Similarly, any Licensor is entitled to offer to its licensees, under
-its sole responsibility, a warranty, that shall only be binding upon
-itself, for the redistribution of the Software and/or the Modified
-Software, under terms and conditions that it is free to decide. Said
-warranty, and the financial terms and conditions of its application,
-shall be subject of a separate instrument executed between the Licensor
-and the Licensee.
-
-
-    Article 8 - LIABILITY
-
-8.1 Subject to the provisions of Article 8.2, the Licensee shall be
-entitled to claim compensation for any direct loss it may have suffered
-from the Software as a result of a fault on the part of the relevant
-Licensor, subject to providing evidence thereof.
-
-8.2 The Licensor's liability is limited to the commitments made under
-this Agreement and shall not be incurred as a result of in particular:
-(i) loss due the Licensee's total or partial failure to fulfill its
-obligations, (ii) direct or consequential loss that is suffered by the
-Licensee due to the use or performance of the Software, and (iii) more
-generally, any consequential loss. In particular the Parties expressly
-agree that any or all pecuniary or business loss (i.e. loss of data,
-loss of profits, operating loss, loss of customers or orders,
-opportunity cost, any disturbance to business activities) or any or all
-legal proceedings instituted against the Licensee by a third party,
-shall constitute consequential loss and shall not provide entitlement to
-any or all compensation from the Licensor.
-
-
-    Article 9 - WARRANTY
-
-9.1 The Licensee acknowledges that the scientific and technical
-state-of-the-art when the Software was distributed did not enable all
-possible uses to be tested and verified, nor for the presence of
-possible defects to be detected. In this respect, the Licensee's
-attention has been drawn to the risks associated with loading, using,
-modifying and/or developing and reproducing the Software which are
-reserved for experienced users.
-
-The Licensee shall be responsible for verifying, by any or all means,
-the suitability of the product for its requirements, its good working
-order, and for ensuring that it shall not cause damage to either persons
-or properties.
-
-9.2 The Licensor hereby represents, in good faith, that it is entitled
-to grant all the rights over the Software (including in particular the
-rights set forth in Article 5).
-
-9.3 The Licensee acknowledges that the Software is supplied "as is" by
-the Licensor without any other express or tacit warranty, other than
-that provided for in Article 9.2 and, in particular, without any warranty
-as to its commercial value, its secured, safe, innovative or relevant
-nature.
-
-Specifically, the Licensor does not warrant that the Software is free
-from any error, that it will operate without interruption, that it will
-be compatible with the Licensee's own equipment and software
-configuration, nor that it will meet the Licensee's requirements.
-
-9.4 The Licensor does not either expressly or tacitly warrant that the
-Software does not infringe any third party intellectual property right
-relating to a patent, software or any other property right. Therefore,
-the Licensor disclaims any and all liability towards the Licensee
-arising out of any or all proceedings for infringement that may be
-instituted in respect of the use, modification and redistribution of the
-Software. Nevertheless, should such proceedings be instituted against
-the Licensee, the Licensor shall provide it with technical and legal
-assistance for its defense. Such technical and legal assistance shall be
-decided on a case-by-case basis between the relevant Licensor and the
-Licensee pursuant to a memorandum of understanding. The Licensor
-disclaims any and all liability as regards the Licensee's use of the
-name of the Software. No warranty is given as regards the existence of
-prior rights over the name of the Software or as regards the existence
-of a trademark.
-
-
-    Article 10 - TERMINATION
-
-10.1 In the event of a breach by the Licensee of its obligations
-hereunder, the Licensor may automatically terminate this Agreement
-thirty (30) days after notice has been sent to the Licensee and has
-remained ineffective.
-
-10.2 A Licensee whose Agreement is terminated shall no longer be
-authorized to use, modify or distribute the Software. However, any
-licenses that it may have granted prior to termination of the Agreement
-shall remain valid subject to their having been granted in compliance
-with the terms and conditions hereof.
-
-
-    Article 11 - MISCELLANEOUS
-
-
-      11.1 EXCUSABLE EVENTS
-
-Neither Party shall be liable for any or all delay, or failure to
-perform the Agreement, that may be attributable to an event of force
-majeure, an act of God or an outside cause, such as defective
-functioning or interruptions of the electricity or telecommunications
-networks, network paralysis following a virus attack, intervention by
-government authorities, natural disasters, water damage, earthquakes,
-fire, explosions, strikes and labor unrest, war, etc.
-
-11.2 Any failure by either Party, on one or more occasions, to invoke
-one or more of the provisions hereof, shall under no circumstances be
-interpreted as being a waiver by the interested Party of its right to
-invoke said provision(s) subsequently.
-
-11.3 The Agreement cancels and replaces any or all previous agreements,
-whether written or oral, between the Parties and having the same
-purpose, and constitutes the entirety of the agreement between said
-Parties concerning said purpose. No supplement or modification to the
-terms and conditions hereof shall be effective as between the Parties
-unless it is made in writing and signed by their duly authorized
-representatives.
-
-11.4 In the event that one or more of the provisions hereof were to
-conflict with a current or future applicable act or legislative text,
-said act or legislative text shall prevail, and the Parties shall make
-the necessary amendments so as to comply with said act or legislative
-text. All other provisions shall remain effective. Similarly, invalidity
-of a provision of the Agreement, for any reason whatsoever, shall not
-cause the Agreement as a whole to be invalid.
-
-
-      11.5 LANGUAGE
-
-The Agreement is drafted in both French and English and both versions
-are deemed authentic.
-
-
-    Article 12 - NEW VERSIONS OF THE AGREEMENT
-
-12.1 Any person is authorized to duplicate and distribute copies of this
-Agreement.
-
-12.2 So as to ensure coherence, the wording of this Agreement is
-protected and may only be modified by the authors of the License, who
-reserve the right to periodically publish updates or new versions of the
-Agreement, each with a separate number. These subsequent versions may
-address new issues encountered by Free Software.
-
-12.3 Any Software distributed under a given version of the Agreement may
-only be subsequently distributed under the same version of the Agreement
-or a subsequent version.
-
-
-    Article 13 - GOVERNING LAW AND JURISDICTION
-
-13.1 The Agreement is governed by French law. The Parties agree to
-endeavor to seek an amicable solution to any disagreements or disputes
-that may arise during the performance of the Agreement.
-
-13.2 Failing an amicable solution within two (2) months as from their
-occurrence, and unless emergency proceedings are necessary, the
-disagreements or disputes shall be referred to the Paris Courts having
-jurisdiction, by the more diligent Party.
-
-
-Version 1.0 dated 2006-09-05.
+If you modify a copy of the Library, and, in your modifications, a facility
+refers to a function or data to be supplied by an Application that uses the
+facility (other than as an argument passed when the facility is invoked),
+then you may convey a copy of the modified version:
+
+    a) under this License, provided that you make a good faith effort to
+    ensure that, in the event an Application does not supply the function or
+    data, the facility still operates, and performs whatever part of its
+    purpose remains meaningful, or
+
+    b) under the GNU GPL, with none of the additional permissions of this
+    License applicable to that copy.
+
+3. Object Code Incorporating Material from Library Header Files.
+
+The object code form of an Application may incorporate material from a header
+file that is part of the Library. You may convey such object code under terms
+of your choice, provided that, if the incorporated material is not limited to
+numerical parameters, data structure layouts and accessors, or small macros,
+inline functions and templates (ten or fewer lines in length),
+you do both of the following:
+
+    a) Give prominent notice with each copy of the object code that the Library
+    is used in it and that the Library and its use are covered by this License.
+
+    b) Accompany the object code with a copy of the GNU GPL
+    and this license document.
+
+4. Combined Works.
+
+You may convey a Combined Work under terms of your choice that, taken together,
+effectively do not restrict modification of the portions of the Library
+contained in the Combined Work and reverse engineering for debugging such
+modifications, if you also do each of the following:
+
+    a) Give prominent notice with each copy of the Combined Work that the
+    Library is used in it and that the Library and its use are covered
+    by this License.
+
+    b) Accompany the Combined Work with a copy of the GNU GPL and
+    this license document.
+
+    c) For a Combined Work that displays copyright notices during execution,
+    include the copyright notice for the Library among these notices, as well
+    as a reference directing the user to the copies of the GNU GPL
+    and this license document.
+
+    d) Do one of the following:
+
+        0) Convey the Minimal Corresponding Source under the terms of this
+        License, and the Corresponding Application Code in a form suitable
+        for, and under terms that permit, the user to recombine or relink
+        the Application with a modified version of the Linked Version to
+        produce a modified Combined Work, in the manner specified by section 6
+        of the GNU GPL for conveying Corresponding Source.
+
+        1) Use a suitable shared library mechanism for linking with the
+        Library. A suitable mechanism is one that (a) uses at run time a
+        copy of the Library already present on the user's computer system,
+        and (b) will operate properly with a modified version of the Library
+        that is interface-compatible with the Linked Version.
+
+    e) Provide Installation Information, but only if you would otherwise be
+    required to provide such information under section 6 of the GNU GPL, and
+    only to the extent that such information is necessary to install and
+    execute a modified version of the Combined Work produced by recombining
+    or relinking the Application with a modified version of the Linked Version.
+    (If you use option 4d0, the Installation Information must accompany the
+    Minimal Corresponding Source and Corresponding Application Code. If you
+    use option 4d1, you must provide the Installation Information in the
+    manner specified by section 6 of the GNU GPL for
+    conveying Corresponding Source.)
+
+5. Combined Libraries.
+
+You may place library facilities that are a work based on the Library side by
+side in a single library together with other library facilities that are not
+Applications and are not covered by this License, and convey such a combined
+library under terms of your choice, if you do both of the following:
+
+    a) Accompany the combined library with a copy of the same work based on
+    the Library, uncombined with any other library facilities, conveyed under
+    the terms of this License.
+
+    b) Give prominent notice with the combined library that part of it is a
+    work based on the Library, and explaining where to find the accompanying
+    uncombined form of the same work.
+
+6. Revised Versions of the GNU Lesser General Public License.
+
+The Free Software Foundation may publish revised and/or new versions of the
+GNU Lesser General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Library as you
+received it specifies that a certain numbered version of the GNU Lesser
+General Public License “or any later version” applies to it, you have the
+option of following the terms and conditions either of that published version
+or of any later version published by the Free Software Foundation. If the
+Library as you received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser General
+Public License ever published by the Free Software Foundation.
+
+If the Library as you received it specifies that a proxy can decide whether
+future versions of the GNU Lesser General Public License shall apply, that
+proxy's public statement of acceptance of any version is permanent
+authorization for you to choose that version for the Library.
 
 #}
diff --git a/README.rst b/README.rst
index bd8d4b518aa95c4f8a94c88f73bc5ea9eedac21f..2062c3879090acbd53408b3ba92d82bda6d11cab 100644
--- a/README.rst
+++ b/README.rst
@@ -1,32 +1,50 @@
-========================
-``bvpy``
-========================
+Welcome
+-------
 
-.. {# pkglts, doc
+**Bvpy** is a python library, based on `FEniCS <https://fenicsproject.org/>`_ and `GMSH <http://gmsh.info/>`_, to easily implement and study numerically Boundary Value Problems (BVPs) as well as Initial Boundary Value Problems (IBVPs) through the Finite Element Method (FEM).
 
-.. #}
+Initially built up in the context of developmental biology and morphomechanics, **Bvpy** proposes an intuitive API as close as possible to the formal definition of BVPs. Its purpose is to enable users to quickly and efficiently estimate the behavior of a wide variety of fields (scalars, vectors, tensors) on biologically relevant structures (Riemannian and non-Rieamannian manifolds), inspired by biophysical and biochemical processes (morphogene patterning, active matter mechanics, diffusion-reaction processes, active transports...).
 
-Package providing tools to solve BVP (Boundary Value Problem) and IBVP (Initial Boundary Value Problem).
+Despite this biological motivation, the **Bvpy** library has been implemented in an *agnostic* manner that makes it suitable for many other scientific context.
 
-It rely on FEniCS finite element library (https://fenicsproject.org/) and GMSH (http://gmsh.info/).
 
-The doc can be found  `here <https://mosaic.gitlabpages.inria.fr/bvpy>`_.
+.. sidebar:: **Bvpy**
 
-:Lead Development: Florian Gacon
+  :Lead Development: Florian Gacon
 
-:Coordination: Olivier Ali
+  :Coordination: Olivier Ali
 
-:Active team: Inria project team `Mosaic <https://team.inria.fr/mosaic/>`_
+  :Active team: Inria project team `Mosaic <https://team.inria.fr/mosaic/>`_
 
-:Institutes: `Inria <http://www.inria.fr>`_
+  :Institutes: `Inria <http://www.inria.fr>`_
 
-:Language: Python
+  :Language: Python
 
-:Supported OS: Linux, MacOS
+  :Supported OS: Linux, MacOS
 
-:Licence: `Cecill-C`
+  :Licence: `LGPL`
 
-:Funding: Inria ADT Gnomon (Christophe Godin)
+  :Funding: Inria ADT Gnomon (Christophe Godin)
+
+
+Documentation
+-------------
+
+The documentation can be found `here <https://mosaic.gitlabpages.inria.fr/bvpy/>`_.
+
+A quick introduction to boundary-value problems and Initial-Boundary-value problems, as well as a the general philosophy behind the **Bvpy** library development can be found `here <https://mosaic.gitlabpages.inria.fr/bvpy/general_introduction.html>`_.
+
+A detailed description of the main classes and models of the library can be found `here <https://mosaic.gitlabpages.inria.fr/bvpy/library_description.html>`_.
+
+Some tutorials explaining basic manipulations are gathered `here <https://mosaic.gitlabpages.inria.fr/bvpy/tutorials.html>`_.
+
+
+Requirements
+------------
+
+* Python 3.7+
+* FEniCS
+* GMSH
 
 
 Installation
@@ -49,7 +67,7 @@ You will need conda in order to install ``bvpy``, you can download it  `here <ht
 
 .. code-block:: bash
 
-	conda install -c mosaic -c conda-forge bvpy
+	conda install -c conda-forge bvpy
 
 
 - From docker
@@ -72,14 +90,46 @@ And use the URL on the terminal to open the tutorials.
 Test
 ----
 
-If you install ``bvpy`` from anaconda install ``nose`` and ``coverage`` inside your environment:
+If you install ``bvpy`` from anaconda install ``pytest`` and ``pytest-cov`` inside your environment:
 
 .. code-block:: bash
 
-	conda install -c conda-forge nose coverage
+	conda install -c conda-forge pytest pytest-cov
 
 And then run at the root of the project:
 
 .. code-block:: bash
 
-	nosetests -v
+	pytest -v
+
+Support
+-------
+
+If you encounter an error or need help, please `raise an issue <https://gitlab.com/oali/bvpy/-/issues>`_.
+
+Contributing
+------------
+
+Bvpy is a is a collaborative project and contributions are welcome. If you want to contribute, please contact the `coordinator <mailto:olivier.ali@inria.fr>`_ prior to any `merge request <https://gitlab.com/oali/bvpy/-/merge_requests>`_ and
+check the `gitlab merge request guidelines <https://docs.gitlab.com/ee/development/contributing/merge_request_workflow.html>`_ if needed.
+
+Citation
+--------
+
+If you are using Bvpy in a published work, please use this bibtex entry as reference:
+
+.. code-block::
+
+  @article{gacon_2021,
+  doi = {xxx},
+  url = {xxx},
+  year = {2021},
+  publisher = {xxx},
+  volume = {x},
+  number = {xx},
+  pages = {xxxx},
+  author = {Florian Gacon and Christophe Godin and Olivier Ali},
+  title = {BVPy: A FEniCS-based Python package to ease the expression and study
+           of boundary value problems in Biology.},
+  journal = {Journal of Open Source Software}
+  }
diff --git a/conda/env.yaml b/conda/env.yaml
index 5e09a91862b4454fd59d3237da93147f286efb4f..07e741f7c36eb0a523c7f43db9e4d8888ceb1ee5 100644
--- a/conda/env.yaml
+++ b/conda/env.yaml
@@ -7,8 +7,8 @@ dependencies:
   - ipython
   - numpy
   - python-gmsh=4.5.6
-  - nose
-  - coverage
+  - pytest
+  - pytest-cov
   - meshio
   - fenics
   - importlib_metadata
@@ -21,3 +21,4 @@ dependencies:
   - make
   - nbsphinx
   - recommonmark
+  - typing_extensions
diff --git a/conda/meta.yaml b/conda/meta.yaml
index 2a376a77f3fda8957fdcc195f95738ad913daadf..e01a6c42b29e4178fb66988f6b0f1870df3a32b0 100644
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@@ -20,7 +20,6 @@ requirements:
   - ipython
   - numpy
   - python-gmsh=4.5.6
-  - nose
   - meshio
   - fenics
   - importlib_metadata
diff --git a/data/tutorial_results/tuto_4_result_1_paraview_snapshot.png b/data/tutorial_results/tuto_6_result_1_paraview_snapshot.png
similarity index 100%
rename from data/tutorial_results/tuto_4_result_1_paraview_snapshot.png
rename to data/tutorial_results/tuto_6_result_1_paraview_snapshot.png
diff --git a/data/tutorial_results/tuto_4_result_3_paraview_snapshot.png b/data/tutorial_results/tuto_6_result_3_paraview_snapshot.png
similarity index 100%
rename from data/tutorial_results/tuto_4_result_3_paraview_snapshot.png
rename to data/tutorial_results/tuto_6_result_3_paraview_snapshot.png
diff --git a/data/tutorial_results/tuto_4_result_4_paraview_snapshot.png b/data/tutorial_results/tuto_6_result_4_paraview_snapshot.png
similarity index 100%
rename from data/tutorial_results/tuto_4_result_4_paraview_snapshot.png
rename to data/tutorial_results/tuto_6_result_4_paraview_snapshot.png
diff --git a/data/tutorial_results/tuto_4_result_5_paraview_snapshot.png b/data/tutorial_results/tuto_6_result_5_paraview_snapshot.png
similarity index 100%
rename from data/tutorial_results/tuto_4_result_5_paraview_snapshot.png
rename to data/tutorial_results/tuto_6_result_5_paraview_snapshot.png
diff --git a/data/tutorial_results/tuto_7_result_3_paraview_snapshot.gif.png b/data/tutorial_results/tuto_7_result_3_paraview_snapshot.gif.png
new file mode 100644
index 0000000000000000000000000000000000000000..3587eaaf4dc5593faae6246fb4caa0b74b2df6e5
Binary files /dev/null and b/data/tutorial_results/tuto_7_result_3_paraview_snapshot.gif.png differ
diff --git a/doc/_static/logo_repo.png b/doc/_static/logo_repo.png
new file mode 100644
index 0000000000000000000000000000000000000000..f4a1fdd2c2ed54d7a2f4921408824810d0e6cb8d
Binary files /dev/null and b/doc/_static/logo_repo.png differ
diff --git a/doc/conf.py b/doc/conf.py
index 24f77e2ec61c7ebb9b8062adf6328c3017675be6..32e2bdf9cb750718490debc73f8355254de6d813 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -22,6 +22,12 @@ import sys
 def setup(app):
     app.add_css_file('css/custom.css')
 
+# -------------------- Import Readme --------------------
+
+
+os.system('cp ../README.rst .')
+print('import file: README.rst')
+
 # -------------------- Import Notebooks --------------------
 
 for file in os.listdir('.'):
@@ -36,11 +42,12 @@ to_import = ['bvpy_tutorial_domain.ipynb',
              'bvpy_tutorial_3_vforms.ipynb',
              'bvpy_tutorial_4_bnd_cond.ipynb',
              'bvpy_tutorial_5_reaction_diffusion.ipynb',
-             'bvpy_tutorial_6_linear_elasticity.ipynb']
+             'bvpy_tutorial_6_linear_elasticity.ipynb',
+             'bvpy_tutorial_7_hyper_elasticity.ipynb']
 
 for i in notebooks:
     if '.' != i[0] and i in to_import:
-        print('import file:'+i)
+        print('import file: '+i)
         os.system('cp ../tutorials/'+i+' '+i)
 
 
@@ -150,9 +157,9 @@ copyright = u"2020, bvpy"
 #
 
 # The short X.Y version.
-version = "1.0.0"
+version = "1.0.1"
 # The full version, including alpha/beta/rc tags.
-release = "1.0.0"
+release = "1.0.1"
 
 
 # The language for content autogenerated by Sphinx. Refer to documentation
diff --git a/doc/index.rst b/doc/index.rst
index 15a84b57555bad75935eb9037368f6725e07566a..ba09098b29e67eccc110674599039675a916c1a3 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -1,29 +1,23 @@
-.. include:: substitutions.txt
+.. image:: https://gitlab.inria.fr/mosaic/bvpy/badges/develop/coverage.svg
+  :target: https://gitlab.inria.fr/mosaic/bvpy
 
-========================
-Bvpy
-========================
-
-.. image:: https://gitlab.inria.fr/mosaic/work-in-progress/bvpy/badges/develop/coverage.svg
-  :target: https://gitlab.inria.fr/mosaic/work-in-progress/bvpy
-
-.. image:: https://gitlab.inria.fr/mosaic/work-in-progress/bvpy/badges/develop/pipeline.svg
-  :target: https://gitlab.inria.fr/mosaic/work-in-progress/bvpy/pipelines
+.. image:: https://gitlab.inria.fr/mosaic/bvpy/badges/master/pipeline.svg
+  :target: https://gitlab.inria.fr/mosaic/bvpy/pipelines
 
 .. image:: https://anaconda.org/mosaic/bvpy/badges/platforms.svg
-  :target: https://anaconda.org/MOSAIC/bvpy
+  :target: https://anaconda.org/conda-forge/bvpy
 
-.. image:: https://anaconda.org/mosaic/bvpy/badges/version.svg
-  :target: https://anaconda.org/MOSAIC/bvpy
+.. image:: https://img.shields.io/conda/vn/conda-forge/bvpy.svg
+  :target: https://anaconda.org/conda-forge/bvpy
 
 .. image:: https://anaconda.org/mosaic/bvpy/badges/license.svg
-  :target: https://cecill.info/licences/Licence_CeCILL-C_V1-en.txt
+  :target: https://www.gnu.org/licenses/lgpl-3.0.txts
 
-.. image:: https://anaconda.org/mosaic/bvpy/badges/latest_release_date.svg
-  :target: https://anaconda.org/MOSAIC/bvpy
+.. image:: https://anaconda.org/conda-forge/bvpy/badges/latest_release_date.svg
+  :target: https://anaconda.org/conda-forge/bvpy
 
-.. image:: https://anaconda.org/mosaic/bvpy/badges/downloads.svg
-  :target: https://anaconda.org/MOSAIC/bvpy
+.. image:: https://anaconda.org/conda-forge/bvpy/badges/downloads.svg
+  :target: https://anaconda.org/conda-forge/bvpy
 
 .. toctree::
    :caption: Table of contents
@@ -37,86 +31,4 @@ Bvpy
    Sources <_dvlpt/modules>
 
 
-Welcome
--------
-
-**Bvpy** is a python library, based on `FEniCS <https://fenicsproject.org/>`_ and `GMSH <http://gmsh.info/>`_, to easily implement and study numerically Boundary Value Problems (BVPs) as well as Initial Boundary Value Problems (IBVPs) through the Finite Element Method (FEM).
-
-Initially built up in the context of developmental biology and morphomechanics, **Bvpy** proposes an intuitive API as close as possible to the formal definition of BVPs given in eq.$(\ref{eq:bvp})$. Its purpose is to enable users to quickly and efficiently estimate the behavior of a wide variety of fields (scalars, vectors, tensors) on biologically relevant structures (Riemannian and non-Rieamannian manifolds), inspired by biophysical and biochemical processes (morphogene patterning, active matter mechanics, diffusion-reaction processes, active transports...).
-
-Despite this biological motivation, the **Bvpy** library has been implemented in an *agnostic* manner that makes it suitable for many other scientific context.
-
-
-.. sidebar:: |name|
-
-  :Lead Development: Florian Gacon
-
-  :Coordination: Olivier Ali
-
-  :Active team: Inria project team `Mosaic <https://team.inria.fr/mosaic/>`_
-
-  :Institutes: `Inria <http://www.inria.fr>`_
-
-  :Language: Python
-
-  :Supported OS: Linux, MacOS
-
-  :Licence: `Cecill-C`
-
-  :Funding: Inria ADT Gnomon (Christophe Godin)
-
-
-Requirements
-------------
-
-* Python 3.6+
-* FEniCS
-* GMSH
-
-
-Installation
-------------
-
-In order to install ``bvpy`` you will need `conda <https://docs.conda.io/projects/conda/en/latest/user-guide/install/>`_.
-
-* From sources:
-
-.. code-block:: bash
-
-  git clone https://gitlab.inria.fr/mosaic/bvpy.git
-  cd bvpy
-  conda env create -f conda/env.yaml -n bvpy-dev
-  conda activate bvpy-dev
-  python setup.py develop --prefix=$CONDA_PREFIX
-
-
-* On `Anaconda Cloud <https://anaconda.org/mosaic/bvpy>`_:
-
-.. code-block:: bash
-
-	conda install -c mosaic -c conda-forge bvpy
-
-
-Documentation
--------------
-
-A quick introduction to boundary-value problems and Initial-Boundary-value problems, as well as a the general philosophy behind the **Bvpy** library development can be found :ref:`here<general_introduction>`.
-
-A detailed description of the main classes and moduels of the library can be found :ref:`here<library_description>`.
-
-Some tutorials explaining basic manipulations are gathered :ref:`here<tutorials>`.
-
-.. Finally, the technical description (docblocks) of all the classes and modules can be found :ref:`here<_src>`.
-
-
-Contents
---------
-
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
-
-Licence
--------
-
-You can distribute and/or modify **Bvpy** under the terms of the Inria licence. If **Bvpy** contributes to a project that leads to a scientific publication, please acknowledge this work by Citing the project.z
+.. include:: README.rst
diff --git a/doc/library_description.rst b/doc/library_description.rst
index 118fa1423d621e2605ef81c0df119e1c39620a5e..2d7014e202aaaf4bee72bcd4946e95ba1eaee9b5 100644
--- a/doc/library_description.rst
+++ b/doc/library_description.rst
@@ -1,5 +1,14 @@
 .. _library_description:
 
+.. math::
+  \newcommand{\euclid}[1]{\mathbb{R}^{#1}}
+  \newcommand{\pos}{\boldsymbol{x}}
+  \newcommand{\field}{u}
+  \newcommand{\source}{g}
+  \newcommand{\body}{\Omega}
+  \newcommand{\surf}{\partial\body}
+  \newcommand{\grad}[1]{\boldsymbol{\nabla}(#1)}
+  \newcommand{\tline}{[t_0, t_{\infty}]}
 
 Library description
 ###################
@@ -42,7 +51,7 @@ The purpose of the first four modules listed below is to provide the library wit
 
 Domains
 -------
-The ``bvpy.domains`` module contains the classes that can be used to implement the notion of integration domain, $\body$ in eq.($\ref{eq:bvp}$). It is built around the (abstract) mother class ``AbstractDomain``. All *domain classes* used as arguments for *BVPs*/*IBVPs* inherits from it. The  ``AbstractDomain`` encompasses two complementary description of the integration domain:
+The ``bvpy.domains`` module contains the classes that can be used to implement the notion of integration domain, :math:`\body` in :ref:`eq.(1)<eq_bvp>`. It is built around the (abstract) mother class ``AbstractDomain``. All *domain classes* used as arguments for *BVPs*/*IBVPs* inherits from it. The  ``AbstractDomain`` encompasses two complementary description of the integration domain:
 
 * A ``gmsh.model`` called the *geometry* of the domain and set by the ``geometry()`` abstract method.
 * A ``fenics.mesh`` called the *mesh* of the domain and produced by the ``discretize()`` method. Note that this method is not automatically called at instanciation but either by high-end functions such as ``plot()`` (from the ``bvpy.utils.visu`` module) or `bvp.solve()`.
@@ -119,7 +128,7 @@ The second class of sub-module  specifically tackles the problem differential sy
   \strain = \text{sym}\big(\grad{\mathbf{u}} \big)
   \end{cases}
 
-where $\mathbf{u}$ stands for the displacement field. The tensor $\stiffness$ accounts for the linear stiffness tensor of the considered material.
+where :math:`\mathbf{u}` stands for the displacement field. The tensor $\stiffness$ accounts for the linear stiffness tensor of the considered material.
 
 The latter encompasses hyper-elastic equations, *i.e.* valid beyond the linear deformation regime. They corresponds to the same general mathematical definition given above but with another strain definition:
 
diff --git a/doc/tutorials.rst b/doc/tutorials.rst
deleted file mode 100644
index 18cb5cc476a5f2f2873d52d9550ae2af4bfda44e..0000000000000000000000000000000000000000
--- a/doc/tutorials.rst
+++ /dev/null
@@ -1,16 +0,0 @@
-.. _tutorials:
-
-Tutorials
-*********
-
-.. toctree::
-	:titlesonly:
-	:maxdepth: 1
-
-	bvpy_tutorial_domain
-	bvpy_tutorial_1_hello_world
-	bvpy_tutorial_2_domains
-	bvpy_tutorial_3_vforms
-	bvpy_tutorial_4_bnd_cond
-	bvpy_tutorial_5_reaction_diffusion
-	bvpy_tutorial_6_linear_elasticity
diff --git a/paper/paper.bib b/paper/paper.bib
index 21b94572881d4933129c2aa51b12b81d00fd8ed4..d190997eb45c91c89d1b7fb9e794952a961b826a 100644
--- a/paper/paper.bib
+++ b/paper/paper.bib
@@ -7,6 +7,93 @@ publisher = {Springer Publishing Company, Incorporated},
 abstract = {This book is a tutorial written by researchers and developers behind the FEniCS Project and explores an advanced, expressive approach to the development of mathematical software. The presentation spans mathematical background, software design and the use of FEniCS in applications. Theoretical aspects are complemented with computer code which is available as free/open source software. The book begins with a special introductory tutorial for beginners. Followingare chapters in Part I addressing fundamental aspects of the approach to automating the creation of finite element solvers. Chapters in Part II address the design and implementation of the FEnicS software. Chapters in Part III present the application of FEniCS to a wide range of applications, including fluid flow, solid mechanics, electromagnetics and geophysics.}
 }
 
+@article{Maas_2012,
+title = {{FEBio: Finite Elements for Biomechanics}},
+author = {Maas, Steve A. and Ellis, Benjamin J. and Ateshian, Gerard A. and Weiss, Jeffrey A.},
+journal = {Journal of Biomechanical Engineering},
+issn = {0148-0731},
+doi = {10.1115/1.4005694},
+pmid = {22482660},
+pmcid = {PMC3705975},
+abstract = {{In the field of computational biomechanics, investigators have primarily used commercial software that is neither geared toward biological applications nor sufficiently flexible to follow the latest developments in the field. This lack of a tailored software environment has hampered research progress, as well as dissemination of models and results. To address these issues, we developed the FEBio software suite (http://mrl.sci.utah.edu/software/febio), a nonlinear implicit finite element (FE) framework, designed specifically for analysis in computational solid biomechanics. This paper provides an overview of the theoretical basis of FEBio and its main features. FEBio offers modeling scenarios, constitutive models, and boundary conditions, which are relevant to numerous applications in biomechanics. The open-source FEBio software is written in C++, with particular attention to scalar and parallel performance on modern computer architectures. Software verification is a large part of the development and maintenance of FEBio, and to demonstrate the general approach, the description and results of several problems from the FEBio Verification Suite are presented and compared to analytical solutions or results from other established and verified FE codes. An additional simulation is described that illustrates the application of FEBio to a research problem in biomechanics. Together with the pre- and postprocessing software PRE VIEW and POST VIEW , FEBio provides a tailored solution for research and development in computational biomechanics.}},
+pages = {011005},
+number = {1},
+volume = {134},
+local-url = {file://localhost/Users/oali/Documents/Work/Biblio/Library/Papers%20Library/maas_jbme_2012.pdf},
+year = {2012}
+}
+
+@article{Badia2020,
+  doi = {10.21105/joss.02520},
+  url = {https://doi.org/10.21105/joss.02520},
+  year = {2020},
+  publisher = {The Open Journal},
+  volume = {5},
+  number = {52},
+  pages = {2520},
+  author = {Santiago Badia and Francesc Verdugo},
+  title = {Gridap: An extensible Finite Element toolbox in Julia},
+  journal = {Journal of Open Source Software}
+}
+
+
+@article{hecht_2012,
+  AUTHOR = {Hecht, F.},
+  TITLE = {New development in FreeFem++},
+  JOURNAL = {J. Numer. Math.},
+  FJOURNAL = {Journal of Numerical Mathematics},
+  VOLUME = {20}, YEAR = {2012},
+  NUMBER = {3-4}, PAGES = {251--265},
+  ISSN = {1570-2820},
+  MRCLASS = {65Y15},
+  MRNUMBER = {3043640},
+  URL = {https://freefem.org/}
+}
+
+@article{Sapala:2018gl,
+title = {{Why plants make puzzle cells, and how their shape emerges.}},
+author = {Sapala, Aleksandra and Runions, Adam and Routier-Kierzkowska, Anne-Lise and Gupta, Mainak Das and Hong, Lilan and Hofhuis, Hugo and Verger, Stéphane and Mosca, Gabriella and Li, Chun-Biu and Hay, Angela and Hamant, Olivier and Roeder, Adrienne Hk and Tsiantis, Miltos and Prusinkiewicz, Przemyslaw and Smith, Richard S},
+journal = {eLife},
+doi = {10.7554/elife.32794},
+abstract = {{The shape and function of plant cells are often highly interdependent. The puzzle-shaped cells that appear in the epidermis of many plants are a striking example of a complex cell shape, however their functional benefit has remained elusive. We propose that these intricate forms provide an effective strategy to reduce mechanical stress in the cell wall of the epidermis. When tissue-level growth is isotropic, we hypothesize that lobes emerge at the cellular level to prevent formation of large isodiametric cells that would bulge under the stress produced by turgor pressure. Data from various plant organs and species support the relationship between lobes and growth isotropy, which we test with mutants where growth direction is perturbed. Using simulation models we show that a mechanism actively regulating cellular stress plausibly reproduces the development of epidermal cell shape. Together, our results suggest that mechanical stress is a key driver of cell-shape morphogenesis.}},
+pages = {2061},
+volume = {7},
+language = {English},
+local-url = {file://localhost/Users/oali/Documents/Work/Biblio/Library/Papers%20Library/Smith/Sapala_2018.pdf},
+year = {2018}
+}
+
+@article{deReuille:2015gu,
+title = {{MorphoGraphX: A platform for quantifying morphogenesis in 4D}},
+author = {Reuille, Pierre Barbier de and Routier-Kierzkowska, Anne-Lise and Kierzkowski, Daniel and Bassel, George W and Schuepbach, Thierry and Tauriello, Gerardo and Bajpai, Namrata and Strauss, Soeren and Weber, Alain and Kiss, Annamaria and Burian, Agata and Hofhuis, Hugo and Sapala, Aleksandra and Lipowczan, Marcin and Heimlicher, Maria B and Robinson, Sarah and Bayer, Emmanuelle M and Basler, Konrad and Koumoutsakos, Petros and Roeder, Adrienne H K and Aegerter-Wilmsen, Tinri and Nakayama, Naomi and Tsiantis, Miltos and Hay, Angela and Kwiatkowska, Dorota and Xenarios, Ioannis and Kuhlemeier, Cris and Smith, Richard S},
+journal = {eLife},
+doi = {10.7554/elife.05864},
+abstract = {{Morphogenesis emerges from complex multiscale interactions between genetic and mechanical processes. To understand these processes, the evolution of cell shape, proliferation and gene expression must be quantified. This quantification is usually performed either in full 3D, which is computationally expensive and technically challenging, or on 2D planar projections, which introduces geometrical artifacts on highly curved organs. Here we present MorphoGraphX ( www.MorphoGraphX.org), a software that bridges this gap by working directly with curved surface images extracted from 3D data. In addition to traditional 3D image analysis, we have developed algorithms to operate on curved surfaces, such as cell segmentation, lineage tracking and fluorescence signal quantification. The software's modular design makes it easy to include existing libraries, or to implement new algorithms. Cell geometries extracted with MorphoGraphX can be exported and used as templates for simulation models, providing a powerful platform to investigate the interactions between shape, genes and growth.}},
+volume = {4},
+language = {English},
+local-url = {file://localhost/Users/oali/Documents/Work/Biblio/Library/Papers%20Library/Smith/Reuille_2015.pdf},
+year = {2015}
+}
+
+@incollection{faure:hal-00681539,
+  TITLE = {{SOFA: A Multi-Model Framework for Interactive Physical Simulation}},
+  AUTHOR = {Faure, Fran{\c c}ois and Duriez, Christian and Delingette, Herv{\'e} and Allard, J{\'e}r{\'e}mie and Gilles, Benjamin and Marchesseau, St{\'e}phanie and Talbot, Hugo and Courtecuisse, Hadrien and Bousquet, Guillaume and Peterlik, Igor and Cotin, St{\'e}phane},
+  URL = {https://hal.inria.fr/hal-00681539},
+  BOOKTITLE = {{Soft Tissue Biomechanical Modeling for Computer Assisted Surgery}},
+  EDITOR = {Yohan Payan},
+  PUBLISHER = {{Springer}},
+  SERIES = {Studies in Mechanobiology, Tissue Engineering and Biomaterials},
+  VOLUME = {11},
+  PAGES = {283-321},
+  YEAR = {2012},
+  MONTH = Jun,
+  DOI = {10.1007/8415\_2012\_125},
+  KEYWORDS = {Motion},
+  PDF = {https://hal.inria.fr/hal-00681539/file/main.pdf},
+  HAL_ID = {hal-00681539},
+  HAL_VERSION = {v1},
+}
+
 @article{Geuzaine:2009,
 author = {Geuzaine, Christophe and Remacle, Jean‐François},
 title = {{Gmsh: A 3‐D finite element mesh generator with built‐in pre‐ and post‐processing facilities}},
diff --git a/paper/paper.md b/paper/paper.md
index b4b7566ac24ac4d0c3aecab1b4b3689dd0949df2..cd8af926e6b2eaf584cfe822172f5a287f023a8f 100644
--- a/paper/paper.md
+++ b/paper/paper.md
@@ -1,5 +1,5 @@
 ---
-title: 'BVPy: A Python package based on FEniCS and Gmsh for solving boundary value problems in cell and tissue mechanics.'
+title: 'BVPy: A FEniCS-based Python package to ease the expression and study of boundary value problems in Biology.'
 tags:
   - Python
   - Biomechanics
@@ -18,7 +18,8 @@ authors:
     orcid: 0000-0001-7671-7225
     affiliation: 1
 affiliations:
-  - name: Laboratoire de Reproduction et Développement des plantes, Univ Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRA, Inria, 69342 Lyon, France
+  - name: Laboratoire Reproduction et Développement des Plantes, Université de Lyon, Ecole Normale Supérieure de Lyon, UCB Lyon 1, CNRS, INRAE, INRIA;  46, allée d'Italie,
+69364 LYON Cedex 07, France.
     index: 1
 date: 22 October 2020
 bibliography: paper.bib
@@ -26,18 +27,30 @@ bibliography: paper.bib
 
 # Summary
 
-Bvpy is a python library, based on FEniCS [@Logg:2012], Gmsh [@Geuzaine:2009] and Meshio [@Schlomer:2020], to easily implement and study numerically Boundary Value Problems (*BVPs*) and Initial Boundary Value Problems (*IBVPs*) through the Finite Element Method (*FEM*). Initially built up in the context of developmental biology and morphomechanics, Bvpy proposes an intuitive Application Programming Interface (*API*) as close as possible to the formal definition of *BVPs* given in \autoref{eq:bvp}. Its purpose is to enable users to quickly and efficiently estimate the behavior of a wide variety of fields (scalars, vectors, tensors) on biologically relevant structures, inspired by biophysical and biochemical processes (morphogene patterning, active matter mechanics, active transports...). Despite this biological motivation, the Bvpy library has been implemented in an *agnostic* manner that makes it suitable for many other scientific context.
+BVPy is a python library to easily implement and study numerically Boundary Value Problems (*BVPs*) and Initial Boundary Value Problems (*IBVPs*) through the Finite Element Method (*FEM*). BVPy proposes an intuitive Application Programming Interface (*API*) to harness and combine the core functionalities of three powerful libraries: FEniCS [@Logg:2012] provides the core data structures and solving algorithms; Gmsh [@Geuzaine:2009] defines the domains and their meshing; and Meshio [@Schlomer:2020] handles data reading and writing. Initially built in the context of developmental biology and morphomechanics, its purpose is to enable all users, even with little to none experience in *FEM*, to quickly and efficiently estimate the behavior of a wide variety of fields (scalars, vectors, tensors) on biologically relevant structures, inspired by biophysical and biochemical processes (morphogene patterning, active matter mechanics, active transports...). Despite this biological motivation, the BVPy library has been implemented in an *agnostic* manner that makes it suitable for many other scientific context.
 
 # Statement of need
 
-*FEMs* are becoming an ubiquitous tool in many research areas and the need for a platform accessible to a large audience is growing. Such platforms should: *(i)* enable people with various level of expertise to contribute; and *(ii)* serve as a communication and education mean between experienced users and novices. Bvpy aims at fulfilling these needs through a twofold strategy:
+*FEMs* are becoming an ubiquitous tool in many research areas and the need for a platform accessible to a large audience of non-specialists is growing. Such platforms should: *(i)* provide an easy access to top-tier, open source, existing libraries, usually restricted to expert users; *(ii)* serve as a means for experienced researchers to communicate with and educate novices. BVPy aims to fulfil these needs through a twofold strategy:
 
-* Provide a high-level *API* to soften the learning curve of *FEMs*. So users with little-to-none knowledge can parametrize, run and monitor simulations based on built-in templates.
+* Provide a high-level *API* to soften the learning curve of *FEMs*. So users with little to none knowledge can parametrize, run and monitor simulations based on built-in templates.
 * Enable users experienced in *FEM*-based modeling to develop and fully customize *de novo* templates that could, in turn, be used by non-specialists.
 
+#  State of the field
+
+Although an exhaustive inventory of the existing *FEM* solutions is far beyond the reach of this paper, it is worth mentioning a few existing alternatives; especially in the context of computational morphomechanics, in which BVPy has been developed.
+
+Some *FEM*-frameworks provide *GUI* such as Sofa [@faure:hal-00681539], FEBio [@Maas_2012] or MorphoMechanX (see [@Sapala:2018gl] for example of use) the yet-to-be published add-on to the MorphoGraphX software [@deReuille:2015gu]. While such userfriendly solutions are accessible to a wide range of users; their "monolithic" construction reduces versatility and prevents more experienced users to tune them at will. Moreover, their code sources are not always freely available.
+
+Besides these *GUI*-based solutions, some "lighter" and more open frameworks, such as FREEFEM [@hecht_2012] or FEniCS [@Logg:2012] are also available. Such solutions appear more versatile and opensource. But their use is usually restricted to people already familiar with the theory of Finite Elements.
+
+A gap currently exists between "closed", userfriendly, softwares on one hand and "open", technical frameworks on the other. BVPy attempts to fill this gap. Our strategy was to harness FEniCS richness into an intuitive and evolutive *API*. To that end, the library architecture maps the conceptual mathematical components of *BVPs* to abstract classes built on FEniCS and Gmsh core components (see the library general description below). From these abstract classes, experienced users can implement *concrete* classes, dedicated to their very specific need. Once implemented, these *concrete* classes can easily be instantiated and combined in a manner that do not require specific *FEM*-related skills (see the Example of use below). Although the library comes with a few of these *concrete* classes, covering *classic* equations (*e.g.* Poisson's and Helmholtz's equations, linear elasticity and hyperelasticity, reaction-diffusion equations...); the full potential of the library resides in its ability to integrate *de novo* classes addressing genuine equations and problems.
+
+
+
 #  Library general description
 
-By definition, a *BVP* corresponds to a (set) differential equation(s) together with a set of constrains defined on the boundaries of the integration domain, see \autoref{eq:bvp}.
+By definition, a *BVP* corresponds to a (set) differential equation(s) together with a set of constraints defined on the boundaries of the integration domain, see \autoref{eq:bvp}.
 
 \begin{equation}\label{eq:bvp}
 \begin{cases}
@@ -58,33 +71,33 @@ u(t,\boldsymbol{x}) = u_N(t,\boldsymbol{x}) & \text{on} \quad \partial\Omega_N\\
 \end{cases}
 \end{equation}
 
-the integration domain and its boundaries can be split into one-dimension *time-*line and orthogonal *spatial* hyperplanes: $\Omega = [t_0, t_{\infty}]\times\Omega_{t_0}$, $\partial\Omega_N=[t_0, t_{\infty}]\times\partial\Omega_{t_0N}$ and $\partial\Omega_D=[t_0, t_{\infty}]\times\partial\Omega_{t_0D}$. Assuming that the spatial hyperplane $\Omega_{t_0}$ does not to evolve with time, such *IBVPs* can be implemented by coupling the resolution of time-dependent *BVPs* with time-integration schemes.
+the integration domain and its boundaries can be split into one-dimension *time-*line and orthogonal *spatial* hyperplanes: $\Omega = [t_0, t_{\infty}]\times\Omega_{t_0}$, $\partial\Omega_N=[t_0, t_{\infty}]\times\partial\Omega_{t_0N}$ and $\partial\Omega_D=[t_0, t_{\infty}]\times\partial\Omega_{t_0D}$. Assuming that the spatial hyperplane $\Omega_{t_0}$ does not evolve with time, such *IBVPs* can be implemented by coupling the resolution of time-dependent *BVPs* with time-integration schemes.
 
 
 ## Scope and range
 
-**In terms of equations:** In its current version (1.0.0), Bvpy can handle non-homogeneous second order Partial Differential Equations (*PDEs*) with potentially non-linear first order terms but linear higher order ones. The unknown  within these *PDEs* can be scalar fields, vector fields or second-order tensor fields. Systems of coupled *PDEs* can be implemented. We extended the initial notion of *BVP* to encompass *IBVPs* featuring first-order time derivatives, as described by \autoref{eq:ibvp}.
+**In terms of equations:** In its current version (1.0.0), BVPy can handle non-homogeneous second order Partial Differential Equations (*PDEs*) with potentially non-linear first order terms but linear higher order ones. The unknowns within these *PDEs* can be scalar fields, vector fields or second-order tensor fields. Systems of coupled *PDEs* can be implemented. We extended the initial notion of *BVP* to encompass *IBVPs* featuring first-order time derivatives, as described by \autoref{eq:ibvp}.
 
-**In terms of domains:** Bvpy provides geometrical primitives to generate simple 2D and 3D domains (rectangles, cubes, ellipoids, torus, ...) as well as basic Constructive Solid Geometry (*CSG*) functionalities (addition, substraction and intersection) to combine these primitives into more complex geometries. The library can also handle triangulated, as well as piecewise-polygonal, meshes generated elsewhere. In the current version, `.txt` , `.ply` and `.obj` files are accepted as inputs.
+**In terms of domains:** BVPy provides geometrical primitives to generate simple 2D and 3D domains (rectangles, cubes, ellipoids, torus, ...) as well as basic Constructive Solid Geometry (*CSG*) functionalities (addition, substraction and intersection) to combine these primitives into more complex geometries. The library can also handle triangulated, as well as piecewise-polygonal, meshes generated elsewhere. In the current version, `.txt` , `.ply` and `.obj` files are accepted as inputs. This versatility comes however with a drawback, for the domain generation procedure is not yet compatible with parallelization processes.
 
-**In terms of boundary conditions:** Classic *Dirichlet* and *Neumann* boundary conditions can be implemented as well as combinaisons of these along the domain boundaries. Periodic boundary conditions can also be defined.
+**In terms of boundary conditions:** Classic *Dirichlet* and *Neumann* boundary conditions can be implemented as well as combinations of these along the domain boundaries. Periodic boundary conditions can also be defined.
 
 ## Organization
 
-Bvpy has been developed and organized as close as possible to the *BVP*  and *IBVP* mathematical formulations given in \autoref{eq:bvp} and \autoref{eq:ibvp}. The library is built around the following key components, see \autoref{fig:schema}:
+BVPy has been developed and organized as close as possible to the *BVP*  and *IBVP* mathematical formulations given in \autoref{eq:bvp} and \autoref{eq:ibvp}. The library is built around the following key components, see \autoref{fig:schema}:
 
-*  Two main classes, named `BVP` and `IBVP`, that respectively encapsulate FEniCS implementations of \autoref{eq:bvp}\& \autoref{eq:ibvp}.
+*  Two main classes, named `BVP` and `IBVP`, that respectively encapsulate FEniCS implementations of \autoref{eq:bvp} \& \autoref{eq:ibvp}.
 *  Three main modules, `bvpy.domains`, `bvpy.vforms` and `bvpy.boundary_conditions`  emulating the corresponding mathematical components of a *bvp*.
 *  A `bvpy.solvers` module where a collection of FEniCS-based linear and non-linear solvers are available.
 
 Besides these main components, the `bvpy.utils` module gathers some useful "housekeeping" functions regrouped thematically into sub-modules, *e.g.* `visu`, `io`... See online documentation.
 
-![Representation of the main modules and submodules forming the Bvpy library. \label{fig:schema}](./fig1.png)
+![Representation of the main modules and submodules forming the BVPy library. \label{fig:schema}](./fig1.png)
 
 
 ## Example of use
 
-The following example is loosely inspired by [@Zhao:2020]. The idea is to estimate mechanical stress distribution within a 2D cross section of a pressurized plant tissue with heterogeneous rigidity. The goal here is to illustrate how parcimonuous, intuitive and yet insightful, Bvpy simulations can be.
+The following example is loosely inspired by [@Zhao:2020]. The idea is to estimate the mechanical stress distribution within a 2D cross section of a pressurized plant tissue with heterogeneous rigidity. The goal here is to illustrate how parsimonious, intuitive and yet insightful, BVPy simulations can be.
 
 ```python
 from bvpy import *
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 79c826e4091e66ed7aa23decf1b36510bd397e3d..0000000000000000000000000000000000000000
--- a/requirements.txt
+++ /dev/null
@@ -1,18 +0,0 @@
-# {# pkglts, reqs
-# requirements are managed by pkglts, do not edit this file at all
-# edit .pkglts/pkg_cfg instead
-# section reqs
-
-# doc
-sphinx
-
-# dvlpt
-
-# install
-
-# test
-coverage
-mock
-nose
-
-# #}
diff --git a/requirements_minimal.txt b/requirements_minimal.txt
deleted file mode 100644
index ef47c7c2ba6bde30dba8b10ace3b7d0b464bc0f2..0000000000000000000000000000000000000000
--- a/requirements_minimal.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-# {# pkglts, reqs
-
-
-# #}
diff --git a/setup.cfg b/setup.cfg
index 88246440bf9738d15e9df823e51dd310ad34c014..6a0849feefa5506f2e6d037838f6526d595d7ec1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -11,21 +11,13 @@ build-dir=build/sphinx
 # #}
 # {# pkglts, test.nose
 
-[nosetests]
-verbosity=1
-detailed-errors=1
-
-with-coverage=1
-cover-erase=1
-# cover-inclusive=1
-cover-package=bvpy
-
-# debug=nose.loader
-# pdb=1
-# pdb-failures=1
 # #}
 
 
 # {# pkglts, test.pytest
+[aliases]
+test=pytest
 
+[tool:pytest]
+addopts = --maxfail=2 -rf
 # #}
diff --git a/setup.py b/setup.py
index d88c094b76bbba8ddba21edc8b7f5f2c3e9988bc..e7c15882eadea509ead39edb53438e3f8cef1997 100644
--- a/setup.py
+++ b/setup.py
@@ -12,10 +12,9 @@ history = open('HISTORY.rst').read()
 pkgs = find_packages('src')
 
 
-
 setup_kwds = dict(
     name='bvpy',
-    version="1.0.0",
+    version="1.0.1",
     description=short_descr,
     long_description=readme + '\n\n' + history,
     author="Florian Gacon",
@@ -25,21 +24,14 @@ setup_kwds = dict(
     zip_safe=False,
 
     packages=pkgs,
-    
+
     package_dir={'': 'src'},
-    setup_requires=[
-        ],
+    setup_requires=[],
     install_requires=[
         ],
-    tests_require=[
-        "coverage",
-        "mock",
-        "nose",
-        ],
+    tests_require=[],
     entry_points={},
     keywords='',
-    
-    test_suite='nose.collector',
     )
 # #}
 # change setup_kwds below before the next pkglts tag
diff --git a/src/bvpy/boundary_conditions/boundary.py b/src/bvpy/boundary_conditions/boundary.py
index 2ddb5fa4dc69a2517957ea5d7f26195fd1c2d3d9..945f3c64791baa678d54946df0ce053e77e02542 100644
--- a/src/bvpy/boundary_conditions/boundary.py
+++ b/src/bvpy/boundary_conditions/boundary.py
@@ -11,11 +11,12 @@
 #           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/boundary_conditions/dirichlet.py b/src/bvpy/boundary_conditions/dirichlet.py
index e4fb6786c519fa1407b74449fbb024072efafb74..fa95b9b8dd3548f2fc290c9e97e5a0b9b7e4a9ce 100644
--- a/src/bvpy/boundary_conditions/dirichlet.py
+++ b/src/bvpy/boundary_conditions/dirichlet.py
@@ -11,18 +11,18 @@
 #           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 import fenics as fe
 from bvpy import logger
 
 import numpy as np
-import numpy.linalg as lng
 from sympy import Symbol
 from sympy.parsing.sympy_parser import parse_expr
 
diff --git a/src/bvpy/boundary_conditions/neumann.py b/src/bvpy/boundary_conditions/neumann.py
index d351a864569233f1c51c65c4b32a696c2a19cf86..522a71e8e865382bb6abacbbe6dfd717cefab7ae 100644
--- a/src/bvpy/boundary_conditions/neumann.py
+++ b/src/bvpy/boundary_conditions/neumann.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 import fenics as fe
diff --git a/src/bvpy/boundary_conditions/periodic.py b/src/bvpy/boundary_conditions/periodic.py
index e0f9ad4ee61246f9399989ea9c1dd44fd0a0818c..08f768a573b9c8382223745ba270a6f754e1baf3 100644
--- a/src/bvpy/boundary_conditions/periodic.py
+++ b/src/bvpy/boundary_conditions/periodic.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#          Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/bvp.py b/src/bvpy/bvp.py
index 0fa496bdac2e57e314e45481769f9905e4b8ec05..606b13bdeb8f1377ea804b592208693c9ae8f5b1 100644
--- a/src/bvpy/bvp.py
+++ b/src/bvpy/bvp.py
@@ -8,13 +8,15 @@
 #
 #       File contributor(s):
 #           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 import fenics as fe
@@ -201,7 +203,7 @@ class BVP(object):
             self.domain.discretize()
 
     def set_vform(self, vform: AbstractVform, constrain=None):
-        """Short summary.
+        """Instanciates the function space and the trial and test functions.
 
         Parameters
         ----------
diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py
index f598cc3ece34538334626a329730a981dc70c448..4db6b4b385bea6371e8460563e346f0cbcf9086a 100644
--- a/src/bvpy/domains/abstract.py
+++ b/src/bvpy/domains/abstract.py
@@ -11,11 +11,12 @@
 #           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/domains/custom_domain.py b/src/bvpy/domains/custom_domain.py
index c1c54b9090fb66ce1246dfde1856e52a652abc28..82e2906bb9eb9f859428f594ce2e591d108a401b 100644
--- a/src/bvpy/domains/custom_domain.py
+++ b/src/bvpy/domains/custom_domain.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 import gmsh
diff --git a/src/bvpy/domains/custom_polygonal.py b/src/bvpy/domains/custom_polygonal.py
index dadfa7508e2eaca7a732546aabd4d7a3c33ea0e0..258d49b5270751bafa6d7db6da814537a83584f2 100644
--- a/src/bvpy/domains/custom_polygonal.py
+++ b/src/bvpy/domains/custom_polygonal.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 import numpy as np
diff --git a/src/bvpy/domains/geometry.py b/src/bvpy/domains/geometry.py
index 1a8d80926004c6799e1715f98547f5e96dd77897..ac779377df8d1e0a3caa3f73ca7f06f6ff38111c 100644
--- a/src/bvpy/domains/geometry.py
+++ b/src/bvpy/domains/geometry.py
@@ -11,12 +11,12 @@
 #           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
 #           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/domains/primitives.py b/src/bvpy/domains/primitives.py
index b37f8b4c8f1dd877dced8203a0a93d8db7ff899e..1444f46a82bac9718bf49e611ee64eae08311ef3 100644
--- a/src/bvpy/domains/primitives.py
+++ b/src/bvpy/domains/primitives.py
@@ -11,12 +11,12 @@
 #           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
 #           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/ibvp.py b/src/bvpy/ibvp.py
index 8437aef27f31ac77ae322af9bdc251301a6566d5..bc5d58f584cf99a19700d8b819ffe34c2e0f6fdd 100644
--- a/src/bvpy/ibvp.py
+++ b/src/bvpy/ibvp.py
@@ -8,13 +8,15 @@
 #
 #       File contributor(s):
 #           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/solvers/linear.py b/src/bvpy/solvers/linear.py
index 831a4987c7c413f35826164ae0e75f75d1d447fd..82d5102494d35ad2cac5c83481518b516377bd80 100644
--- a/src/bvpy/solvers/linear.py
+++ b/src/bvpy/solvers/linear.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 import fenics as fe
diff --git a/src/bvpy/solvers/nonlinear.py b/src/bvpy/solvers/nonlinear.py
index 7aaf17ba0b5ee899a20ce5c9ad8ebac4764ba38a..980255a75f51f439c6ba3d231626b4a7d2f7b12e 100644
--- a/src/bvpy/solvers/nonlinear.py
+++ b/src/bvpy/solvers/nonlinear.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/solvers/parameter.py b/src/bvpy/solvers/parameter.py
index 237257cf28b1c0da5fd2b9583649290a08cbde10..150bca4cdbf5b6904a614e5b03d938995ec01458 100644
--- a/src/bvpy/solvers/parameter.py
+++ b/src/bvpy/solvers/parameter.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
@@ -80,7 +81,4 @@ class SolverParameter():
         None
 
         """
-        try:
-            self._param[item] = value
-        except KeyError:
-            raise AttributeError(item)
+        self._param[item] = value
diff --git a/src/bvpy/solvers/time_integration.py b/src/bvpy/solvers/time_integration.py
index 9de6e1f650eef929b00450703ea2446ae073281c..0bbb8aaa7f6a2aa2c456a487933a592bb598b5fa 100644
--- a/src/bvpy/solvers/time_integration.py
+++ b/src/bvpy/solvers/time_integration.py
@@ -10,11 +10,12 @@
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/utils/io.py b/src/bvpy/utils/io.py
index 501de51f5e23cd8f8aed5e9b27a7a21fe516a8a6..39fb924474fe2cdf5ca24c2a8aa00caa8622d480 100644
--- a/src/bvpy/utils/io.py
+++ b/src/bvpy/utils/io.py
@@ -1,3 +1,24 @@
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+#       bvpy.utils.interface
+#
+#       File author(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#
+#       File contributor(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       File maintainer(s):
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
+#       See accompanying file LICENSE.txt or copy at
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
+#
+# -----------------------------------------------------------------------
 import numpy as np
 import fenics as fe
 from bvpy.domains.abstract import AbstractDomain
diff --git a/src/bvpy/utils/post_processing.py b/src/bvpy/utils/post_processing.py
index a4e54d37c286dc77a943cda0b1d5dd5faf715ea4..85c75e723066969bfa73ec3b508101aa2dcbf63f 100644
--- a/src/bvpy/utils/post_processing.py
+++ b/src/bvpy/utils/post_processing.py
@@ -1,9 +1,29 @@
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+#       bvpy.utils.interface
+#
+#       File author(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#
+#       File contributor(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       File maintainer(s):
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
+#       See accompanying file LICENSE.txt or copy at
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
+#
+# -----------------------------------------------------------------------
 import fenics as fe
 import numpy as np
 import operator
 
-from bvpy.utils.pre_processing import (create_expression,
-                                       _expression_from_function)
+from bvpy.utils.pre_processing import create_expression
 
 
 class SolutionExplorer(object):
@@ -145,7 +165,8 @@ class SolutionExplorer(object):
         comparator : string
             The comparison to apply (could be '<', '<=', '>', '>=' or '==')
         sort : string
-            Sort the coordinates and the Function on the specify axe (could be 'x', 'y' or 'z').
+            Sort the coordinates and the Function on the specify axe
+            (could be 'x', 'y' or 'z').
             Do not sort if None.
 
         Returns
@@ -155,10 +176,11 @@ class SolutionExplorer(object):
 
         """
         d = dict(x=0, y=1, z=2)
-        c = {'>': operator.gt, '>=':operator.ge, '<':operator.lt, '<=':operator.le, '==': operator.eq}
+        c = {'>': operator.gt, '>=': operator.ge, '<': operator.lt,
+             '<=': operator.le, '==': operator.eq}
         assert comparator in c.keys()
         assert axe in d.keys()
-        mask = c[comparator](self.get_vertex_position()[:,d[axe]],threshold)
+        mask = c[comparator](self.get_vertex_position()[:, d[axe]], threshold)
 
         pos = self.get_vertex_position()[mask]
         val = self.to_vertex_values()[mask]
diff --git a/src/bvpy/utils/pre_processing.py b/src/bvpy/utils/pre_processing.py
index d94eddad352852735211fb11934567eb141b9919..90bd0b8ea56211f5bd8ab7a05e0a8a2c4ab58da7 100644
--- a/src/bvpy/utils/pre_processing.py
+++ b/src/bvpy/utils/pre_processing.py
@@ -8,13 +8,15 @@
 #
 #       File contributor(s):
 #           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
@@ -57,7 +59,7 @@ def _expression_from_string(source, functionspace, degree, **kwargs):
 def _expression_from_list(source, functionspace, degree):
     source = np.asarray(source)
     dim = source[0].shape
-    if len(source.shape)==1:
+    if len(source.shape) == 1:
         source = np.array([[i] for i in source])
 
     class UserExpr(fe.UserExpression):
@@ -81,7 +83,6 @@ def _expression_from_list(source, functionspace, degree):
     else:
         return expr
 
-
 def _expression_from_function(source, functionspace, degree):
     x = [0, 0, 0]
     dim = np.array(source(x)).shape
@@ -97,7 +98,7 @@ def _expression_from_function(source, functionspace, degree):
             def value_shape(self):
                 return dim
 
-    elif len(dim)==1:
+    elif len(dim) == 1:
 
         class UserExpr(fe.UserExpression):
             def eval(self, value, x):
@@ -108,7 +109,7 @@ def _expression_from_function(source, functionspace, degree):
             def value_shape(self):
                 return dim
 
-    elif len(dim)==2:
+    elif len(dim) == 2:
 
         class UserExpr(fe.UserExpression):
             def eval(self, value, x):
diff --git a/src/bvpy/utils/progress_bar.py b/src/bvpy/utils/progress_bar.py
index f23ff73c82b7ce43f71cd41d90ba1878e661516f..e165d35b52c217cbc1bb5b7c682db508d5610efa 100644
--- a/src/bvpy/utils/progress_bar.py
+++ b/src/bvpy/utils/progress_bar.py
@@ -1,3 +1,23 @@
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+#       bvpy.utils.interface
+#
+#       File author(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#
+#       File contributor(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#
+#       File maintainer(s):
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
+#       See accompanying file LICENSE.txt or copy at
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
+#
+# -----------------------------------------------------------------------
 from math import floor
 
 
diff --git a/src/bvpy/utils/tensor.py b/src/bvpy/utils/tensor.py
new file mode 100644
index 0000000000000000000000000000000000000000..5f73635f789058c68f20b320b970a30a19fdba63
--- /dev/null
+++ b/src/bvpy/utils/tensor.py
@@ -0,0 +1,163 @@
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+#       bvpy.utils.interface
+#
+#       File author(s):
+#           Olivier Ali <olviier.ali@inria.fr>
+#
+#       File contributor(s):
+#           Olivier Ali <olviier.ali@inria.fr>
+#           Florian Gacon <florian.gacon@inria.fr>
+#
+#       File maintainer(s):
+#           Olivier Ali <olviier.ali@inria.fr>
+#
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
+#       See accompanying file LICENSE.txt or copy at
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
+#
+# -----------------------------------------------------------------------
+
+import numpy as np
+import numpy.linalg as lng
+import fenics as fe
+
+
+def applyElementwise(f, T):
+    """Apply a function to each component of tensor T.
+
+    Parameters
+    ----------
+    f : function applied to tensor
+
+    T : UFL Tensor
+
+    Returns
+    -------
+    UFL Tensor
+
+    """
+
+    shape = T.ufl_shape
+    print('shape tensor:', shape)
+
+    t_array = []
+    for i in range(shape[0]):
+        t_j = []
+        for j in range(shape[1]):
+            t_j.append(f(T[i, j]))
+        t_array.append(t_j)
+
+    return fe.as_tensor(t_array)
+
+
+def amplitude(tensor):
+    """Measures the amplitude of a tensor.
+
+    We choose to define the amplitude of a tensor as its Frobenius norm divided
+    by the square root of its dimension.
+    So the identity tensor has always a unit amplitude.
+
+    Parameters
+    ----------
+    tensor : numpy.nDarray
+
+    Returns
+    -------
+    float
+
+    """
+    dim, _ = tensor.shape
+
+    assert(dim == _)
+
+    return lng.norm(tensor) / np.sqrt(dim)
+
+
+def sph(tensor):
+    """Extract the spherical part of a symmetric tensor.
+
+    Parameters
+    ----------
+    tensor : numpy.nDarray
+
+    Returns
+    -------
+    numpy.nDarray
+    """
+
+    dim, _ = tensor.shape
+
+    assert(dim == _)
+
+    return np.trace(tensor) * np.eye(dim) / dim
+
+
+def dev(tensor):
+    """Extracts the deviatoric part of a symmetric tensor.
+
+    Parameters
+    ----------
+    tensor : numpy.nDarray
+
+    Returns
+    -------
+    numpy.nDarray
+    """
+    dim, _ = tensor.shape
+
+    assert(dim == _)
+
+    return tensor - sph(tensor)
+
+
+def intensity(tensor):
+    """Measures the intesity of a tensor.
+
+    We chose to define the intensity of a symmetric tensor as the amplitude of
+    its spherical part.
+
+    Parameters
+    ----------
+    tensor : numpy.nDarray
+
+    Returns
+    -------
+    float
+
+    """
+    dim, _ = tensor.shape
+
+    assert(dim == _)
+
+    return amplitude(sph(tensor))
+
+
+def anisotropy(tensor):
+    """Measures the anisotropy degree of a tensor.
+
+    We chose to define the anistropy of a symmetric tensor as the ratio
+    between the amplitudes of its deviatoric and spherical parts.
+
+    Parameters
+    ----------
+    tensor : numpy.nDarray
+
+    Returns
+    -------
+    float
+
+    """
+    dim, _ = tensor.shape
+
+    assert(dim == _)
+
+    try:
+        return float(amplitude(dev(tensor))) / float(amplitude(sph(tensor)))
+
+    except ZeroDivisionError:
+        print("WARNING: tensor has a null amplitude.")
+
+        return False
diff --git a/src/bvpy/utils/visu.py b/src/bvpy/utils/visu.py
index 1386be84a8eaabd868b49ab3989738bd6f491a9e..b2fbcdcc3e973a4689be249493db929e0390beee 100644
--- a/src/bvpy/utils/visu.py
+++ b/src/bvpy/utils/visu.py
@@ -1,3 +1,24 @@
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+#       bvpy.utils.interface
+#
+#       File author(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#
+#       File contributor(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       File maintainer(s):
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
+#       See accompanying file LICENSE.txt or copy at
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
+#
+# -----------------------------------------------------------------------
 import fenics as fe
 import numpy as np
 import plotly.io as pio
diff --git a/src/bvpy/vforms/abstract.py b/src/bvpy/vforms/abstract.py
index 9d1ee5c2ed65f8fafaad12ae87c6199b393344b0..d2a07f42f8b901f9c15f31f85447b676be6f758f 100644
--- a/src/bvpy/vforms/abstract.py
+++ b/src/bvpy/vforms/abstract.py
@@ -1,20 +1,22 @@
 # -*- python -*-
 # -*- coding: utf-8 -*-
 #
-#       bvpy.vforms.abstract
+#       bvpy.utils.interface
 #
 #       File author(s):
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File contributor(s):
 #           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 import fenics as fe
diff --git a/src/bvpy/vforms/elasticity.py b/src/bvpy/vforms/elasticity.py
index c992ee7ff700d309a56b04429bd8855d198426d4..05b4bdae65d15ceb52c4c95f9cb663eb94864e5c 100644
--- a/src/bvpy/vforms/elasticity.py
+++ b/src/bvpy/vforms/elasticity.py
@@ -13,19 +13,19 @@
 #       File maintainer(s):
 #           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
 import fenics as fe
-from abc import abstractmethod
 
-from ufl import nabla_grad, nabla_div
+from ufl import nabla_grad, nabla_div, ln
 
 from bvpy.vforms.abstract import AbstractVform, Element
-from bvpy.utils.pre_processing import create_expression
+from bvpy.utils.tensor import applyElementwise
 
 __all__ = ['ElasticForm', 'LinearElasticForm', 'HyperElasticForm']
 
@@ -54,7 +54,9 @@ class ElasticForm(AbstractVform):
         If considering plane stress or not.
 
     """
-    def __init__(self, source=[0, 0, 0], young=1, poisson=.5, plane_stress=False, thickness=1):
+
+    def __init__(self, source=[0, 0, 0], young=1, poisson=.5,
+                 plane_stress=False, thickness=1):
         super().__init__()
         self._elements = [Element('P', 1, 'vector')]
 
@@ -70,14 +72,18 @@ class ElasticForm(AbstractVform):
         self._update_mecha_parameters()
 
     def _update_mecha_parameters(self):
-        self._parameters['mu'] = .5 * self._parameters['young'] / (1 + self._parameters['poisson'])
+        self._parameters['mu'] = .5 * self._parameters['young']\
+                                / (1 + self._parameters['poisson'])
 
         if self.plane_stress:
-            self._parameters['lmbda'] = self._parameters['young'] * self._parameters['poisson']\
-                          / (1 - self._parameters['poisson'] ** 2)
+            self._parameters['lmbda'] = self._parameters['young']\
+                                      * self._parameters['poisson']\
+                                      / (1 - self._parameters['poisson'] ** 2)
         else:
-            self._parameters['lmbda'] = self._parameters['young'] * self._parameters['poisson']\
-                          / ((1 + self._parameters['poisson'])*(1-2*self._parameters['poisson']))
+            self._parameters['lmbda'] = self._parameters['young']\
+                                      * self._parameters['poisson']\
+                                      / ((1 + self._parameters['poisson'])
+                                         * (1-2*self._parameters['poisson']))
 
     @property
     def plane_stress(self):
@@ -213,6 +219,7 @@ class HyperElasticForm(ElasticForm):
         Description of attribute `model`.
 
     """
+
     def __init__(self, source=None, young=1, poisson=.5,
                  material_model='st_venant', thickness=1, plane_stress=False):
         super().__init__(source, young, poisson, plane_stress, thickness)
@@ -220,14 +227,75 @@ class HyperElasticForm(ElasticForm):
         self.elastic_potential = None
         self.model = material_model
 
-    def construct_form(self, u, v, sol):
-        d = sol.geometric_dimension()
+    def _strain(self, u):
+        """Computes the Green Lagrange strain tensor field.
+
+        Parameters
+        ----------
+        u : :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
+            Displacement field.
+
+        type : str
+            Optional. Name of the strain measure we want to compute.
+            The default is 'GL' (Green-Lagrange).
+            Available measures: 'GL', 'biot', 'natural' and 'almansi'.
+
+        Returns
+        -------
+        :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
+            The computed strain.
+
+        """
+        d = u.geometric_dimension()
         Id = fe.Identity(d)
-        F = Id + fe.nabla_grad(sol)
+        F = Id + fe.nabla_grad(u)
         C = F.T * F
 
-        E = 0.5 * (nabla_grad(sol) + nabla_grad(sol).T
-                   + nabla_grad(sol).T * nabla_grad(sol))
+        if type == 'GL':
+            return 0.5 * (C - Id)
+
+        elif type == 'biot':
+            return applyElementwise(lambda x: x**0.5, C) - Id
+
+        elif type == 'natural':
+            return 0.5 * applyElementwise(lambda x: ln(x), C)
+
+        elif type == 'almansi':
+            return 0.5 * (Id - fe.inv(C))
+
+        return 0.5 * (nabla_grad(u) + nabla_grad(u).T
+                      + nabla_grad(u).T * nabla_grad(u))
+
+    def _stress(self, u=None, model='st_venant'):
+        """Computes the second Piola-Kirchhoff stress tensor field.
+
+        Parameters
+        ----------
+        u : :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
+            Displacement field.
+
+        model : str.
+            The name of the considered hyperelastic model.
+            Default value is 'st_venant'.
+
+        Returns
+        -------
+        :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
+            The computed stress.
+
+        """
+        d = u.geometric_dimension()
+
+        if model == 'st_venant':
+            return self._parameters['lmbda'] * fe.tr(self._strain(u))\
+                   * fe.Identity(d) + 2 * self._parameters['mu']\
+                   * self._strain(u)
+        else:
+            print(f"WARNING: stress derivation from the {model} model is"
+                  + "not implemented yet in the stress method.")
+
+    def construct_form(self, u, v, sol):
+        E = self._strain(sol)
 
         # -- Constitutive models
         if self.model == 'st_venant':
@@ -236,6 +304,11 @@ class HyperElasticForm(ElasticForm):
             W += self._parameters['mu'] * fe.tr(E * E)
 
         elif self.model == 'neo_hookean':
+            d = u.geometric_dimension()
+            Id = fe.Identity(d)
+            F = Id + fe.nabla_grad(sol)
+            C = F.T * F
+
             Ic = fe.tr(C)
             J = fe.det(F)
             W = (self._parameters['mu'] / 2) * (Ic - 3)\
diff --git a/src/bvpy/vforms/helmholtz.py b/src/bvpy/vforms/helmholtz.py
index 16c52fcc5601f5a83620ed6f5af6deb6ec68c7b9..7a431768b8a0ffa08fe4455788adbb22223bbde9 100644
--- a/src/bvpy/vforms/helmholtz.py
+++ b/src/bvpy/vforms/helmholtz.py
@@ -11,12 +11,12 @@
 #           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
 #           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 from math import pi
diff --git a/src/bvpy/vforms/plasticity.py b/src/bvpy/vforms/plasticity.py
index 4b4bba665177e432c33d03834c775a997fd38ea9..0a626f741e2abc53d1fe9e984ccc63288f7ff9cb 100644
--- a/src/bvpy/vforms/plasticity.py
+++ b/src/bvpy/vforms/plasticity.py
@@ -1,3 +1,24 @@
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+#       bvpy.bvp
+#
+#       File author(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#
+#       File contributor(s):
+#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       File maintainer(s):
+#           Olivier Ali <olivier.ali@inria.fr>
+#
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
+#       See accompanying file LICENSE.txt or copy at
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
+#
+# -----------------------------------------------------------------------
 from .elasticity import ElasticForm
 import numpy as np
 
@@ -13,6 +34,7 @@ class PlasticForm(ElasticForm):
         The plastic part of the deformation gradient.
 
     """
+
     def __init__(self, source=None, young=1, poisson=.3, plane_stress=False,
                  F_plastic=np.diag([1, 1, 1]), thickness=1):
         super().__init__(source, young, poisson, plane_stress, thickness)
diff --git a/src/bvpy/vforms/poisson.py b/src/bvpy/vforms/poisson.py
index 00a3904170b5d0f3cbcda3a6a465a477a21c8967..e13ff0baa651ae2c7ce878b7c3f9cacbbf3bde42 100644
--- a/src/bvpy/vforms/poisson.py
+++ b/src/bvpy/vforms/poisson.py
@@ -11,12 +11,12 @@
 #           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
 #           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/src/bvpy/vforms/transport.py b/src/bvpy/vforms/transport.py
index cfb4337cea065b8935bb377218a3b182d0e9498c..926099fa0551d8bdecaefdd42fe77681459084e9 100644
--- a/src/bvpy/vforms/transport.py
+++ b/src/bvpy/vforms/transport.py
@@ -1,20 +1,22 @@
 # -*- python -*-
 # -*- coding: utf-8 -*-
 #
-#       bvpy.templates.vforms.poisson
+#       bvpy.bvp
 #
 #       File author(s):
 #           Florian Gacon <florian.gacon@inria.fr>
 #
 #       File contributor(s):
 #           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
 #       File maintainer(s):
-#           Florian Gacon <florian.gacon@inria.fr>
+#           Olivier Ali <olivier.ali@inria.fr>
 #
-#       Distributed under the Cecill-C License.
+#       Copyright © by Inria
+#       Distributed under the LGPL License..
 #       See accompanying file LICENSE.txt or copy at
-#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
 
diff --git a/test/test_boundary.py b/test/test_boundary.py
index ba0531108d4509a8eac3c0e28da7d25e27d931e6..9cab4bd8d213e59dfaf58d862e0ee66a7be50164 100644
--- a/test/test_boundary.py
+++ b/test/test_boundary.py
@@ -1,7 +1,8 @@
 import unittest
-from bvpy.boundary_conditions.boundary import *
+from bvpy.boundary_conditions.boundary import Boundary
 import fenics as fe
 
+
 class TestBoundary(unittest.TestCase):
     def setUp(self):
         """Initiates the test.
@@ -25,7 +26,7 @@ class TestBoundary(unittest.TestCase):
 
     def test_mark(self):
         b1 = Boundary("near(x, 0)")
-        mesh = fe.UnitSquareMesh(2,2)
+        mesh = fe.UnitSquareMesh(2, 2)
         mf = fe.MeshFunction("size_t", mesh, 0, 0)
         b1.markMeshFunction(1, meshfunc=mf)
         self.assertEqual(mf[0], 1)
diff --git a/test/test_bvp.py b/test/test_bvp.py
index 58cd4ecfd27e65a5b55082f986593e0ae7fc28b1..343bf0a27b0259428d907beb4bd978e84d2dd878 100644
--- a/test/test_bvp.py
+++ b/test/test_bvp.py
@@ -27,8 +27,9 @@ class TestPoisson(unittest.TestCase):
         self.assertIsNotNone(bvp.domain.mesh)
         bvp.info()
         print(bvp)
-        
+
         bvp = BVP()
+        self.assertEqual("BVP", bvp.__repr__()[1:4])
         self.assertIsNone(bvp.domain)
         self.assertIsNone(bvp.vform)
 
diff --git a/test/test_dirichlet.py b/test/test_dirichlet.py
index 5507f14a13bd81dc3da943cad74146e85f223baf..d5e2ae9f0312d2ba73e79fbaa6a47848c8541054 100644
--- a/test/test_dirichlet.py
+++ b/test/test_dirichlet.py
@@ -1,10 +1,11 @@
 import unittest
-from bvpy.boundary_conditions.dirichlet import *
+from bvpy.boundary_conditions.dirichlet import NormalDirichlet, ZeroDirichlet
 from bvpy.boundary_conditions import Boundary
 import fenics as fe
 
 from numpy.testing import assert_array_equal
 
+
 class TestDirichlet(unittest.TestCase):
     def setUp(self):
         """Initiates the test.
@@ -15,11 +16,12 @@ class TestDirichlet(unittest.TestCase):
         """
 
     def test_normal(self):
-        mesh = fe.UnitSquareMesh(2,2)
+        mesh = fe.UnitSquareMesh(2, 2)
         V = fe.VectorFunctionSpace(mesh, 'P', 1)
         diri = NormalDirichlet(boundary='all')
         diri = NormalDirichlet(boundary=Boundary('all'))
-        self.assertEqual("<NormalDirichlet object, location: all, value: 1>", diri.__repr__())
+        self.assertEqual("<NormalDirichlet object, location: all, value: 1>",
+                         diri.__repr__())
 
         d = diri.apply(V)
         self.assertEqual(d.get_boundary_values()[5], 1)
@@ -27,5 +29,4 @@ class TestDirichlet(unittest.TestCase):
     def test_zero(self):
         b = Boundary("all")
         diri = ZeroDirichlet(b, shape=3)
-
-        assert_array_equal(diri._val, [0,0,0])
+        assert_array_equal(diri._val, [0, 0, 0])
diff --git a/test/test_elasticity.py b/test/test_elasticity.py
index 6672525362544f70470e90db3b98422189a24b06..8f2e8ca019b48090803e9e056f41db1a82de485c 100644
--- a/test/test_elasticity.py
+++ b/test/test_elasticity.py
@@ -1,11 +1,9 @@
 import unittest
 
 from bvpy import BVP, logger
-from bvpy.domains import Rectangle, CustomDomain
+from bvpy.domains import CustomDomain
 from bvpy.vforms import LinearElasticForm, HyperElasticForm
 from bvpy.boundary_conditions import dirichlet
-from bvpy.utils.pre_processing import create_expression
-import tempfile
 import fenics as fe
 
 
@@ -50,9 +48,6 @@ class TestLinearElasticityProblem(unittest.TestCase):
         err_rel = abs(self.theo_def_body + comp_def) / self.theo_def_body
         self.assertLess(err_rel, 1e-3)
 
-        strain = vform.strain(bvp.solution)
-        stress = vform.stress(bvp.solution)
-
     def test_hyper(self):
         vform = HyperElasticForm(source=[0, -self.rho_g],
                                  young=1e5,
diff --git a/test/test_io.py b/test/test_io.py
index 0b66113e51bfd7af4d64538a11911a709881347b..5a8fa37b72d0ee8341b8ae502df392df279d6031 100644
--- a/test/test_io.py
+++ b/test/test_io.py
@@ -75,6 +75,7 @@ class TestIO(unittest.TestCase):
     def test_numerize(self):
         self.assertEqual(numerize('3'), 3)
         self.assertEqual(numerize('3.4'), 3.4)
+        self.assertEqual(numerize('three'), 'three')
 
     def test_save_domain(self):
         domain = Rectangle()
diff --git a/test/test_pre_process.py b/test/test_pre_process.py
index e1a6957bed346e64cb242c07f799cbee8eccdda1..3b2b314001f5ef8b44dc10e3ed7d8d5d3040df7c 100644
--- a/test/test_pre_process.py
+++ b/test/test_pre_process.py
@@ -10,9 +10,9 @@ class TestPreProcess(unittest.TestCase):
     def setUp(self):
         """Initiates the test.
         """
-        self.mesh = UnitSquareMesh(2,2)
-        self.V = FunctionSpace(self.mesh ,'P', 1)
-        self.VV = VectorFunctionSpace(self.mesh ,'P', 1)
+        self.mesh = UnitSquareMesh(2, 2)
+        self.V = FunctionSpace(self.mesh,'P', 1)
+        self.VV = VectorFunctionSpace(self.mesh,'P', 1)
 
     def tearDown(self):
         """Concludes and closes the test.
diff --git a/test/test_solvers.py b/test/test_solvers.py
index 9f17ee68d4128b09ea90fc366786b69a5de10936..7e05a06cd7d047acc359c648d2d5ce057e758c7a 100644
--- a/test/test_solvers.py
+++ b/test/test_solvers.py
@@ -13,10 +13,10 @@
 # Code:
 
 import unittest
-
 from bvpy import logger
 from bvpy.solvers.linear import LinearSolver
 from bvpy.solvers.nonlinear import NonLinearSolver
+# from bvpy.solvers.parameter import SolverParameter
 import fenics as fe
 
 # -- test class
@@ -62,7 +62,10 @@ class TestLinearSolver(unittest.TestCase):
         solver.solve()
         error_L2 = fe.errornorm(self.u_D, self.u, 'L2')
         logger.debug(error_L2)
-        assert(error_L2<1e-2)
+        assert(error_L2 < 1e-2)
+
+        with self.assertRaises(AttributeError):
+            solver.parameters['fake_parameter']
 
     def test_krylov(self):
         solver = LinearSolver()
@@ -74,7 +77,7 @@ class TestLinearSolver(unittest.TestCase):
         solver.solve()
         error_L2 = fe.errornorm(self.u_D, self.u, 'L2')
         logger.debug(error_L2)
-        assert(error_L2<1e-2)
+        assert(error_L2 < 1e-2)
 
 
 class TestNonLinearSolver(unittest.TestCase):
@@ -90,7 +93,8 @@ class TestNonLinearSolver(unittest.TestCase):
         import sympy as sym
         x, y = sym.symbols('x[0], x[1]')
         u = 1 + x + 2*y
-        f = - sym.diff(q(u)*sym.diff(u, x), x) - sym.diff(q(u)*sym.diff(u, y), y)
+        f = - sym.diff(q(u)*sym.diff(u, x), x)\
+            - sym.diff(q(u)*sym.diff(u, y), y)
         f = sym.simplify(f)
         u_code = sym.printing.ccode(u)
         f_code = sym.printing.ccode(f)
@@ -129,7 +133,6 @@ class TestNonLinearSolver(unittest.TestCase):
         """
         pass
 
-
     def test_GMRES(self):
         solver = NonLinearSolver()
         solver.parameters['krylov_solver']['monitor_convergence'] = False
@@ -137,4 +140,4 @@ class TestNonLinearSolver(unittest.TestCase):
         solver.set_problem(self.u, self.F, self.J, [self.bc], [self.ps])
         solver.solve()
         error_L2 = fe.errornorm(self.u_e, self.u, 'L2')
-        assert(error_L2<1e-10), error_L2
+        assert(error_L2 < 1e-10), error_L2
diff --git a/test/test_tensor.py b/test/test_tensor.py
new file mode 100644
index 0000000000000000000000000000000000000000..16950ecbef1395477218c781f14ee0f170640242
--- /dev/null
+++ b/test/test_tensor.py
@@ -0,0 +1,87 @@
+import unittest
+
+import fenics as fe
+import numpy as np
+from numpy.testing import assert_array_equal
+from bvpy import logger
+from bvpy.utils.tensor import (applyElementwise, amplitude, sph,
+                               dev, intensity, anisotropy)
+
+
+class TestTensor(unittest.TestCase):
+
+    def setUp(self):
+        """Initiates the test.
+        """
+        logger.setLevel(50)
+        self.Zr2 = np.zeros((2, 2))
+        self.Zr3 = np.zeros((3, 3))
+        self.Id2 = np.eye(2)
+        self.Id3 = np.eye(3)
+        self.ev1 = 2
+        self.ev2 = 3
+        self.mtrx = np.array([[self.ev1, 0],
+                              [0, self.ev2]])
+
+    def tearDown(self):
+        """Concludes and closes the test.
+        """
+        pass
+
+    def test_applyEmentwise(self):
+        def square(x):
+            return x**2
+
+        tnsr_fspace = fe.TensorFunctionSpace(fe.UnitSquareMesh(2, 2), 'DG', 0)
+        fe_tensor = fe.Constant([[1, 2], [3, 4]])
+
+        tensor = fe.project(fe_tensor, tnsr_fspace)
+        squared_tensor = applyElementwise(square, tensor)
+
+        for mtrx in fe.project(squared_tensor,
+                               tnsr_fspace).vector().get_local().reshape(-1,
+                                                                         2,
+                                                                         2):
+            assert_array_equal(mtrx, [[1, 4], [9, 16]])
+
+    def test_amplitude(self):
+        self.assertEqual(amplitude(self.Id2), 1)
+        self.assertEqual(amplitude(self.Id3), 1)
+        self.assertEqual(amplitude(self.Zr2), 0)
+        self.assertEqual(amplitude(self.Zr3), 0)
+        self.assertAlmostEqual(amplitude(self.mtrx),
+                               np.sqrt((self.ev1**2 + self.ev2**2) / 2),
+                               delta=1e-7)
+
+    def test_sph(self):
+        assert_array_equal(sph(self.Id2), self.Id2)
+        assert_array_equal(sph(self.Id3), self.Id3)
+        assert_array_equal(sph(self.Zr2), self.Zr2)
+        assert_array_equal(sph(self.Zr3), self.Zr3)
+
+        sph_mtrx = .5 * (self.ev1 + self.ev2) * self.Id2
+        assert_array_equal(sph(self.mtrx), sph_mtrx)
+
+    def test_dev(self):
+        mtrx_dev = .5 * (self.ev1 - self.ev2) * np.array([[1, 0],
+                                                          [0, -1]])
+
+        assert_array_equal(dev(self.Id2), self.Zr2)
+        assert_array_equal(dev(self.Id3), self.Zr3)
+        assert_array_equal(dev(self.mtrx), mtrx_dev)
+
+    def test_intensity(self):
+        self.assertEqual(intensity(self.Id2), 1)
+        self.assertEqual(intensity(self.Id3), 1)
+        self.assertEqual(intensity(self.Zr2), 0)
+        self.assertEqual(intensity(self.Zr3), 0)
+        self.assertEqual(intensity(self.mtrx), .5 * (self.ev1 + self.ev2))
+
+    def test_anisotropy(self):
+        self.assertEqual(anisotropy(self.Id2), 0)
+        self.assertEqual(anisotropy(self.Id3), 0)
+
+        self.assertEqual(anisotropy(self.mtrx),
+                         np.abs((self.ev1 - self.ev2) / (self.ev1 + self.ev2)))
+
+        self.assertFalse(anisotropy(self.Zr2))
diff --git a/test/test_transport.py b/test/test_transport.py
index 7542790ed59bf0cc00976405b03a9e42d09b88b4..5f730ca8cc1da3e5d097178ebba8f82bdb0682be 100644
--- a/test/test_transport.py
+++ b/test/test_transport.py
@@ -88,9 +88,16 @@ class TestTransport(unittest.TestCase):
         state = [[1+0.04*k1**2+0.1*r, 0.2*k2+0.1*r] for r in random]
         init = create_expression(state)
 
+        periodic_bc = SquarePeriodicBoundary(length=10, width=10)
+        x = [10, 10]
+        y = [10, 0]
+        z = [10, 0]
+        periodic_bc.map(x, z)
+        periodic_bc.map(y, z)
+
         ibvp = IBVP(domain=dom,
                     vform=vform,
-                    bc=SquarePeriodicBoundary(length=10, width=10),
+                    bc=periodic_bc,
                     scheme=ImplicitEuler,
                     initial_state=init)
 
diff --git a/test/test_visu.py b/test/test_visu.py
index 1e59cf5412268ef96c95fa71a592e3bb578be7e6..66c8409da6d3a92becf4291585896cd336a1ea19 100644
--- a/test/test_visu.py
+++ b/test/test_visu.py
@@ -11,11 +11,12 @@ import fenics as fe
 
 set_renderer('iframe')
 
+
 class TestVisu(unittest.TestCase):
     def setUp(self):
         """Initiates the test.
         """
-        self.mesh = fe.UnitSquareMesh(2,2)
+        self.mesh = fe.UnitSquareMesh(2, 2)
         self.V = fe.FunctionSpace(self.mesh, 'P', 1)
 
     def tearDown(self):
@@ -24,11 +25,14 @@ class TestVisu(unittest.TestCase):
         pass
 
     def test_cone(self):
-        pos = [[0,0,0],
-               [1,0,0]]
-        val = [[0,0,0.1],
-               [0,0,0.1]]
+        pos = [[0, 0, 0],
+               [1, 0, 0]]
+        val = [[0, 0, 0.1],
+               [0, 0, 0.1]]
+        val_2D = [[0, 0.1],
+                  [0, 0.1]]
         self.assertIsNotNone(_cone_plot(pos, val))
+        self.assertIsNotNone(_cone_plot(pos, val_2D))
 
     def test_wireframe_mesh(self):
         self.assertIsNotNone(_wireframe_plot_mesh(self.mesh))
diff --git a/tutorials/bvpy_tutorial_1_hello_world.ipynb b/tutorials/bvpy_tutorial_1_hello_world.ipynb
index 3cda1ff181ffa2c9304c1619613851432120a36a..6dfff63154185bdee9a12ebf2cebee669237e62a 100644
--- a/tutorials/bvpy_tutorial_1_hello_world.ipynb
+++ b/tutorials/bvpy_tutorial_1_hello_world.ipynb
@@ -23,16 +23,23 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## General introduction\n",
+    "## General introduction"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "Mathematically the Poisson problem reads:\n",
     "\n",
-    "$$\n",
-    "\\begin{cases}\n",
+    "\\begin{equation*}\n",
+    "\\left\\{\n",
+    "\\begin{array}{rl}\n",
     "-\\nabla^2 u(\\mathbf{x}) = f(\\mathbf{x}) &\\text{in } \\Omega \\\\\n",
     "u(\\mathbf{x}) =  u_D(\\mathbf{x}) &\\text{on } \\partial\\Omega,\n",
-    "\\end{cases}\n",
-    "\\label{eq:poisson}\n",
-    "$$\n",
+    "\\end{array}\n",
+    "\\right.\n",
+    "\\end{equation*}\n",
     "\n",
     "where $\\Omega$ corresponds to the integration domain and $\\partial\\Omega$ to its border. $f(\\mathbf{x})$ is called the *source* function and needs to be specified prior to the cell_size. Similarly $u_D(\\mathbf{x})$ is the function that will prescribe the constrained values $u(\\mathbf{x})$ must verify on $\\partial\\Omega$.\n",
     "\n",
@@ -119,7 +126,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "**N.B.:** As mentionned in the FEniCS tutorial, the chosen function $u_D$ here is an exact solution of the considered Poisson equation. This will be useful to estimate the precision of our numerical estimation downrange."
+    "> **Note:** As mentionned in the FEniCS tutorial, the chosen function $u_D$ here is an exact solution of the considered Poisson equation. This will be useful to estimate the precision of our numerical estimation downrange."
    ]
   },
   {
@@ -180,7 +187,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "**Note:** the `set_renderer('notebook')` is mandatory only when working within jupyter notebooks.\n",
+    "> **Note:** the `set_renderer('notebook')` is mandatory only when working within jupyter notebooks.\n",
     "\n",
     "In its current state, **Bvpy** provides minimalist visualization and analysis tools. One may want to turn to other frameworks, such as **Paraview** for instance, to visualize and analyze simulation results. To do so, **Bvpy** enable the recording of simulation in the `.xdmf` format."
    ]
@@ -236,7 +243,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## First attempt to play a bit with the simulation parameters\n",
+    "## First attempt to play a bit with the simulation parameters"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "Let's now look at the simplest modifications we can perform.\n",
     "We can first check all the properties and parameter values of our simulation, as follow:"
    ]
@@ -316,6 +329,158 @@
    "source": [
     "And this time, indeed, we get a very accurate result."
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Inter-operability with FEniCS"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**BVPy** goal is not only to provide a clear and simple *API* on top of the **FEniCS** library but also to make this *API* as transparent as possible. Users already experienced in **FEniCS** can therefore easily access familiar objects they are used to.\n",
+    "\n",
+    "To that end, many attributes from the main **BVPy** classes are **FEniCS** objects. Let's review the main ones."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### FEniCS mesh within the Domain class"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's consider the `rectangle` domain previously defined and check the kind of attributes it contains:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(f\"rectangle is of type {type(rectangle)} \\n\")\n",
+    "\n",
+    "for name, attribute in rectangle.__dict__.items():\n",
+    "    print(f\"{name} is of type {type(attribute)}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As you can see from above, some attributes of the `Domain` class are **Dolfin** object. In particular, one can easily access the **FEniCS** mesh within the domain:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import fenics as fe\n",
+    "\n",
+    "fe_mesh = rectangle.mesh\n",
+    "isinstance(fe_mesh, fe.Mesh)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Developpers used to the **FEniCS** *API* can benefit from it. For instance if one wants to get the positions of the vertices of the mesh defined on the domain, one can use the *classic* functions `fenics.vertices()` and `mesh.coordinates()`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "points = np.array([fe_mesh.coordinates()[vtx.index()][:2] for vtx in fe.vertices(fe_mesh)]).transpose()\n",
+    "\n",
+    "plt.scatter(points[0], points[1])\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Note:** The `cdata` and `bdata` attributes are `MeshFunctionSizet`, *i.e.* **FEniCS** functions defined on a given mesh returning integer values. They can be useful for parametrization, as we will see in tutorial #6."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### FEniCS function space within the BVP class "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Another important concept in *FEM* is the one of *function space*. This concept is implemented within the **FEniCS** library within the `fenics.FunctionSpace` class, inherited from the **Dolfin** class `dolfin.functions.functionspace.FunctionSpace`. \n",
+    "\n",
+    "Within **BVPy** the instanciation of a function space is done in two steps:\n",
+    "* First, when the `vform` is instanced, one can define the family, degree and type of elements to consider, *e.g.* `someForm._elements = [Element('P', 1, 'scalar')]`.\n",
+    "* Then, once the corresponding `BVP` is instanciated, the `set_vform` methods instanciates the `fenics.FunctionSpace` to consider according the the various elements defined in the previous step."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fe_func_space = prblm_1.functionSpace\n",
+    "if isinstance(fe_func_space, fe.FunctionSpace):\n",
+    "    print(\"There is a FEniCS function space deep inside the BVP class!\")\n",
+    "    print(f\"it is of dimension {fe_func_space.dim()}\")\n",
+    "    print(f\"and contains elements of type {fe_func_space.ufl_element()}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "One can also get this function space, once the simulation is done by considering the solution of the considered `BVPy`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol = prblm_1.solution\n",
+    "isinstance(sol.function_space(), fe.FunctionSpace)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Note:** For users not at ease with the **FEniCS** *API*, the `SolutionExplorer` class from the `utils.post_processing` module is designed to handle the conversion between **FEniCS** objects and *more usual* ones such as `numpy.nDarray`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
diff --git a/tutorials/bvpy_tutorial_2_domains.ipynb b/tutorials/bvpy_tutorial_2_domains.ipynb
index 242f4b23e587fef4839cdffe2c4ac9a3799fd312..a21e058a4353330bd7d13c2d2b7ba148a5e4b9c4 100644
--- a/tutorials/bvpy_tutorial_2_domains.ipynb
+++ b/tutorials/bvpy_tutorial_2_domains.ipynb
@@ -282,7 +282,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "dsk1 = Disk(center=[0, -.25, 0])\n",
+    "dsk1 = Disk(center=[0, -.25, 0], clear=True, cell_size=.1)\n",
     "dsk2 = Disk(center=[0, .25, 0])\n",
     "dsk3 = Disk(radius=.3)\n",
     "rec1 = Rectangle(y=-.05, x=-.5, length=1, width=.1)\n",
@@ -291,6 +291,13 @@
     "plot(cmplx_dom)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Note:** You can see in the first line above we added the argument `clear=True` at the first domain instanciation. When multiple domains are instanciated, they share the same **Gmsh** factory and are, so to speak, “super-imposed\" on one another. This can create conflicts latter on when domains are processed. For instance, when the `plot` function probes the factory for a specific domain, It might not return the desired one but all the previously instanced ones contained within the factory. to avoid this, we added a call to the `gmsh.clear()` (triggered with the argument `clear` is set to `True`) method within our `AbstractDomain` class. This **Gmsh** method resets the domain factory and avoids visualization conflits."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -466,7 +473,7 @@
    "toc_cell": false,
    "toc_position": {},
    "toc_section_display": true,
-   "toc_window_display": true
+   "toc_window_display": false
   },
   "varInspector": {
    "cols": {
diff --git a/tutorials/bvpy_tutorial_3_vforms.ipynb b/tutorials/bvpy_tutorial_3_vforms.ipynb
index 6b2fa36a066426f63790e8541e827b48c6da264d..459e89c16fe6b788d35d016cc0367ff5dfa50488 100644
--- a/tutorials/bvpy_tutorial_3_vforms.ipynb
+++ b/tutorials/bvpy_tutorial_3_vforms.ipynb
@@ -101,6 +101,7 @@
     "The [*Helmholtz equation*](https://en.wikipedia.org/wiki/Helmholtz_equation):\n",
     "$$\n",
     " \\Delta u + k(\\mathbf{x})u +f(\\mathbf{x}) = 0 \\quad\\quad \\forall \\mathbf{x} \\in \\Omega\n",
+    " \\tag{1}\n",
     " \\label{eq:hlmtz}\n",
     "$$\n",
     "> **Note:** In general, the coefficient of the first order term (second term in the *lhs* of eq.($\\ref{eq:hlmtz}$)) is constant. Within our framework, a spacial dependency comes for free."
@@ -112,22 +113,25 @@
    "source": [
     "### Defining the variational form to implement\n",
     "\n",
-    "The first thing to do is to compute the variational expression of eq.($\\ref{eq:hlmtz}$). This is done by multiplying eq.($\\ref{eq:hlmtz}$) by a trial function (noted $\\psi$ hereafter) and to integrate the resulting expression on the considered domain $\\Omega$:\n",
+    "The first thing to do is to compute the variational expression of eq.(1). This is done by multiplying eq.(1) by a trial function (noted $\\psi$ hereafter) and to integrate the resulting expression on the considered domain $\\Omega$:\n",
     "$$\n",
     "\\int_{\\Omega}d\\mathbf{x}\\Delta(u) v + \\int_{\\Omega}d\\mathbf{x}k(\\mathbf{x})uv + \\int_{\\Omega}d\\mathbf{x}f(\\mathbf{x})v = 0\n",
+    "\\tag{2}\n",
     "\\label{eq:start_eq}\n",
     "$$\n",
     "\n",
-    "Integration by part of the first term of the *rhs* of eq.($\\ref{eq:start_eq}$) yields:\n",
+    "Integration by part of the first term of the *rhs* of eq.(2) yields:\n",
     "$$\n",
     "\\int_{\\Omega}d\\mathbf{x}\\Delta(u)v =  \\Big[ \\nabla(u)v \\Big]_{\\partial\\Omega}-\\int_{\\Omega}d\\mathbf{x}\\nabla(u)\\nabla(v)\n",
+    "\\tag{3}\n",
     "\\label{eq:ipp}\n",
     "$$\n",
     "\n",
-    "Assuming the integrated term of the *lhs* of eq.($\\ref{eq:ipp}$) vanishes; combining eqs.($\\ref{eq:start_eq} $\\&$\\ref{eq:ipp}$) yields:\n",
+    "Assuming the integrated term of the *lhs* of eq.(3) vanishes; combining eqs.(2 \\& 3) yields:\n",
     "\n",
     "$$\n",
     "\\int_{\\Omega}d\\mathbf{x}\\Big(\\nabla(u) \\nabla(v) - k(\\mathbf{x})uv\\Big) = \\int_{\\Omega}d\\mathbf{x}f(\\mathbf{x})v\n",
+    "\\tag{4}\n",
     "\\label{eq:end_eq}\n",
     "$$\n",
     "\n",
@@ -140,6 +144,7 @@
     "B(u,v) = \\nabla(u) \\nabla(v) - kuv\\\\\n",
     "L(v) = f v\n",
     "\\end{cases}\n",
+    "\\tag{5}\n",
     "\\label{eq:linbilin_terms}\n",
     "$$"
    ]
@@ -172,7 +177,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "* `set_expression(self)`: This method has a purely informative role. It provides the class with a string description of the variational form, close to eq.($\\ref{eq:end_eq}$). This description is printed when the `.info()` method is called. In the present case, it should read:\n",
+    "* `set_expression(self)`: This method has a purely informative role. It provides the class with a string description of the variational form, close to eq.(4). This description is printed when the `.info()` method is called. In the present case, it should read:\n",
     "```python\n",
     "def set_expression(self):\n",
     "    self._expression = \"grad(u)*grad(v)*dx - linear_coef*u*v*dx\"\n",
@@ -336,14 +341,14 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.7.9"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,
    "autoclose": false,
    "autocomplete": true,
    "bibliofile": "biblio.bib",
-   "cite_by": "apalike",
+   "cite_by": "key",
    "current_citInitial": 1,
    "eqLabelWithNumbers": true,
    "eqNumInitial": 1,
diff --git a/tutorials/bvpy_tutorial_4_bnd_cond.ipynb b/tutorials/bvpy_tutorial_4_bnd_cond.ipynb
index d3684da1e62d2edc1cfeb25aad192563295ac47a..a0d9981a3dfa3769e224baf1843e230607f49a15 100644
--- a/tutorials/bvpy_tutorial_4_bnd_cond.ipynb
+++ b/tutorials/bvpy_tutorial_4_bnd_cond.ipynb
@@ -236,7 +236,7 @@
    "toc_cell": false,
    "toc_position": {},
    "toc_section_display": true,
-   "toc_window_display": true
+   "toc_window_display": false
   }
  },
  "nbformat": 4,
diff --git a/tutorials/bvpy_tutorial_5_reaction_diffusion.ipynb b/tutorials/bvpy_tutorial_5_reaction_diffusion.ipynb
index 5c15d1a0137c713764b32bcfb0bd616d089d9e51..96974f1e59c583966e1b364f503b515ec1e294ca 100644
--- a/tutorials/bvpy_tutorial_5_reaction_diffusion.ipynb
+++ b/tutorials/bvpy_tutorial_5_reaction_diffusion.ipynb
@@ -551,7 +551,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.7.9"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,
diff --git a/tutorials/bvpy_tutorial_6_linear_elasticity.ipynb b/tutorials/bvpy_tutorial_6_linear_elasticity.ipynb
index 49014e520de4b003a474e37d1d644684678d8d19..8b7b9eeb8b7de5054a441e78b28433d2b5617b9b 100644
--- a/tutorials/bvpy_tutorial_6_linear_elasticity.ipynb
+++ b/tutorials/bvpy_tutorial_6_linear_elasticity.ipynb
@@ -24,7 +24,7 @@
     "\n",
     "\n",
     "\n",
-    "Download the notebook: [bvpy_tutorial_4](https://gitlab.inria.fr/mosaic/publications/bvpy/-/raw/develop/tutorials/bvpy_tutorial_4_linear_elasticity.ipynb?inline=false)."
+    "Download the notebook: [bvpy_tutorial_6](https://gitlab.inria.fr/mosaic/publications/bvpy/-/raw/develop/tutorials/bvpy_tutorial_6_linear_elasticity.ipynb?inline=false)."
    ]
   },
   {
@@ -262,7 +262,7 @@
    "source": [
     "from bvpy.utils.io import save\n",
     "\n",
-    "path = '../data/tutorial_results/tuto_4_result_1.xdmf'\n",
+    "path = '../data/tutorial_results/tuto_6_result_1.xdmf'\n",
     "save(stretching, path)"
    ]
   },
@@ -274,7 +274,7 @@
    "source": [
     "from IPython.display import Image\n",
     "\n",
-    "Image(filename='../data/tutorial_results/tuto_4_result_1_paraview_snapshot.png')"
+    "Image(filename='../data/tutorial_results/tuto_6_result_1_paraview_snapshot.png')"
    ]
   },
   {
@@ -390,8 +390,8 @@
    "source": [
     "wound_displacement = wound_healing.solution\n",
     "plot(wound_displacement, size=1)\n",
-    "path = '../data/tutorial_results/tuto_4_result_2.xdmf'\n",
-    "save(wound, path)"
+    "path = '../data/tutorial_results/tuto_6_result_2.xdmf'\n",
+    "save(wound_healing, path)"
    ]
   },
   {
@@ -445,7 +445,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "save(stress, '../data/tutorial_results/tuto_4_result_3.xdmf')"
+    "save(stress, '../data/tutorial_results/tuto_6_result_3.xdmf')"
    ]
   },
   {
@@ -454,7 +454,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Image(filename=\"../data/tutorial_results/tuto_4_result_3_paraview_snapshot.png\")"
+    "Image(filename=\"../data/tutorial_results/tuto_6_result_3_paraview_snapshot.png\")"
    ]
   },
   {
@@ -553,7 +553,7 @@
     "displacement = prblm.solution\n",
     "#plot(displacement, size=1, norm=True)\n",
     "\n",
-    "save(displacement, '../data/tutorial_results/tuto_4_result_4.xdmf')"
+    "save(displacement, '../data/tutorial_results/tuto_6_result_4.xdmf')"
    ]
   },
   {
@@ -580,8 +580,8 @@
    "outputs": [],
    "source": [
     "stress = heterogeneous_elastic_response.stress(displacement)\n",
-    "save(stress, \"../data/tutorial_results/tuto_4_result_5.xdmf\")\n",
-    "Image(filename=\"../data/tutorial_results/tuto_4_result_5_paraview_snapshot.png\")"
+    "save(stress, \"../data/tutorial_results/tuto_6_result_5.xdmf\")\n",
+    "Image(filename=\"../data/tutorial_results/tuto_6_result_5_paraview_snapshot.png\")"
    ]
   },
   {
@@ -630,84 +630,15 @@
     "inner_stress = stress_function.cell_label_filter(tissue.cdata, 0)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "code_folding": [
-     3,
-     15,
-     25,
-     34,
-     43
-    ]
-   },
-   "outputs": [],
-   "source": [
-    "# -- Then for each zone, let's compute the stress intensities\n",
-    "import numpy as np\n",
-    "import numpy.linalg as lng\n",
-    "\n",
-    "def amplitude(tensor):\n",
-    "    \"\"\"Measures the amplitude of a tensor.\n",
-    "    \n",
-    "    We choose to define the amplitude of a tensor as its Frobenius norm divided\n",
-    "    by the square root of its dimension. So the identity tensor has always a unit amplitude.\n",
-    "    \"\"\"\n",
-    "    dim, _ = tensor.shape\n",
-    "    \n",
-    "    assert(dim == _)\n",
-    "    \n",
-    "    return lng.norm(tensor) / np.sqrt(dim)\n",
-    "\n",
-    "def sph(tensor):\n",
-    "    \"\"\"Extract the spherical part of a symmetric tensor.\n",
-    "    \"\"\"\n",
-    "    \n",
-    "    dim, _ = tensor.shape\n",
-    "    \n",
-    "    assert(dim == _)\n",
-    "    \n",
-    "    return np.trace(tensor) * np.eye(dim) / dim\n",
-    "\n",
-    "def dev(tensor):\n",
-    "    \"\"\"Extracts the deviatoric part of a symmetric tensor.\n",
-    "    \"\"\"\n",
-    "    dim, _ = tensor.shape\n",
-    "    \n",
-    "    assert(dim == _)\n",
-    "    \n",
-    "    return tensor - sph(tensor)\n",
-    "\n",
-    "def intensity(tensor):\n",
-    "    \"\"\"Measures the intesity of a tensor.\n",
-    "    \"\"\"\n",
-    "    dim, _ = tensor.shape\n",
-    "    \n",
-    "    assert(dim == _)\n",
-    "    \n",
-    "    return amplitude(sph(tensor))\n",
-    "\n",
-    "def anisotropy(tensor):\n",
-    "    \"\"\"Measures the anisotropy degree of a tensor.\n",
-    "    \"\"\"\n",
-    "    dim, _ = tensor.shape\n",
-    "    \n",
-    "    assert(dim == _)\n",
-    "    \n",
-    "    try:\n",
-    "        return amplitude(dev(tensor)) / amplitude(sph(tensor))\n",
-    "\n",
-    "    except ZeroDivisionError:\n",
-    "        print(f\"WARNING: tensor {tensor} has a null amplitude.\")"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
+    "# -- Let's compute the intensity and the anisotropy of the stress tensor field.\n",
+    "from bvpy.utils.tensor import intensity, anisotropy\n",
+    "\n",
     "epiderm_strs_int = list(map(lambda x: intensity(x), epiderm_stress))\n",
     "inner_strs_int = list(map(lambda x: intensity(x), inner_stress))\n",
     "\n",
@@ -878,7 +809,7 @@
    "source": [
     "displacement = pressurized_shell.solution\n",
     "\n",
-    "save(displacement, '../data/tutorial_results/tuto_4_result_6.xdmf')\n",
+    "save(displacement, '../data/tutorial_results/tuto_6_result_6.xdmf')\n",
     "\n",
     "plot(displacement, size=.25, norm=True)"
    ]
@@ -953,6 +884,15 @@
     "#### Problem instanciation and resolution"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prblm.info()"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -971,7 +911,7 @@
    "source": [
     "displacement = prblm.solution\n",
     "\n",
-    "save(displacement, '../data/tutorial_results/tuto_4_result_7.xdmf')\n",
+    "save(displacement, '../data/tutorial_results/tuto_6_result_7.xdmf')\n",
     "\n",
     "plot(displacement, size=1000, norm=True)"
    ]
diff --git a/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb b/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..a4c93ce25742b1c3699831e16ebb1a1650cb3f88
--- /dev/null
+++ b/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb
@@ -0,0 +1,906 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# More advanced examples of use\n",
+    "\n",
+    "In this tutorial we are going to show how one can use of **BVPy** to tackle less obvious problems than Poisson's or linear elasticity. To that end, we propose two examples:\n",
+    "* A quick comparison between linear-elastic and hyperelastic behaviors\n",
+    "* Modelization of the coupling between a reacto-diffusive scalar field and the elastic modulus of the considered medium.\n",
+    "\n",
+    "\n",
+    "\n",
+    "**Covered topics:**\n",
+    "\n",
+    "- Hyper elasticity variational form manipulation: `HyperElasticVform()`, `set_parameters()`, `info()`, `stress()`.\n",
+    "\n",
+    "- ...\n",
+    "\n",
+    "\n",
+    "\n",
+    "Download the notebook: [bvpy_tutorial_7](https://gitlab.inria.fr/mosaic/publications/bvpy/-/raw/develop/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb?inline=false)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Comparison between Linear elastic and hyper-elastic models"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Basic notions about hyperelasticity"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Note:** The following notes are by no means an exhaustive description of the hyperelastic theory. Those are just quick remarks and comments gathered to provide some context to the following developments."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The theory of hyperelasticity can be seen as the generalization of linear elasticity. It relies on two key ingredients: \n",
+    "* The formalization of the elastic behavior of a medium as a potential energy functional. \n",
+    "* The quantification of the deformation (between a resting and a loaded configuration) through a strain tensor fields.\n",
+    "\n",
+    "The general idea is to build from these ingredients a mathematical expression of the potential energy as a function of invariants of the chosen strain tensor field. The combinaison of a particular expression of the energy functional and a strain measure corresponds to a specific hyperelastic model.\n",
+    "\n",
+    "For instance, two common hyperelastic models frequently encounter within the litterature are the **Saint-Venant Kirchhoff** and the **Neo-Hookean** models, which correspond respectively to the following elastic potentials:\n",
+    "$$\n",
+    "\\mathcal{E}_{\\text{SVK}} = \\frac{\\lambda}{2} \\text{tr}(\\boldsymbol{E})^2 + \\mu\\text{tr}(\\boldsymbol{E}^2),\n",
+    "$$\n",
+    "and\n",
+    "$$\n",
+    "\\mathcal{E}_{\\text{NH}} = C_1 \\big(\\text{det}(\\boldsymbol{F})  -3\\big)\n",
+    "$$\n",
+    "where $\\lambda, \\mu, C_1$ are coefficients, $\\boldsymbol{E}$ stands for the **Euler-Lagrange** strain tensor and $\\boldsymbol{F}$ the deformation gradient.\n",
+    "\n",
+    "> **Note:** Both of these models are implemented within the `HyperElasticForm` class from the `bvpy.vforms.elasticity` sub-module."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### The Euler-Lagrange strain measure"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's consider a deformation field $\\varphi: \\boldsymbol{X} \\to \\boldsymbol{x}$ between two manifolds, describing a continuum in its resting and loaded configurations. \n",
+    "\n",
+    "A small vector $\\boldsymbol{dX}$ from the resting configuration is streched into $\\boldsymbol{dx} = \\boldsymbol{F}\\cdot \\boldsymbol{dX}$ in the loaded configuration. $\\boldsymbol{F}$ is a tensor field, called the *deformation gradient*: $\\boldsymbol{F} = \\boldsymbol{\\nabla}(\\varphi)$.\n",
+    "\n",
+    "From here, one can compute the (squared) length of any streched incremental element:\n",
+    "$\n",
+    "dx^2 = ||\\boldsymbol{dx}||^2 = \\boldsymbol{dX}^T\\cdot\\boldsymbol{F}^T\\cdot \\boldsymbol{F}\\cdot \\boldsymbol{dX}.\n",
+    "$\n",
+    "The streching, *i.e.* the difference between the size of the same incremental vector in both configurations, is then quantified by a symmetric second order tensor field:\n",
+    "$$\n",
+    "dl^2 = dx^2 - dX^2 = 2 \\boldsymbol{dX}^T\\cdot\\boldsymbol{E}\\cdot \\boldsymbol{dX}.\n",
+    "$$\n",
+    "\n",
+    "The tensor within the right hand-side term: \n",
+    "$$\n",
+    "\\boldsymbol{E} = \\frac{1}{2}(\\boldsymbol{F}^T\\cdot \\boldsymbol{F} - \\boldsymbol{I}_{\\text{d}}),\n",
+    "$$\n",
+    "is called the **Euler-Lagrange** strain tensor and is a central strain measure in hyperelasticity."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### How is it related to the linear elasticity theory ?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we make the assumption that the deformation between the resting and loaded configuations is small, the gradient $\\boldsymbol{F}$ can be taylor-expanded as follow:\n",
+    "$$\n",
+    "\\boldsymbol{F} \\approx \\boldsymbol{I}_{\\text{d}} + \\boldsymbol{\\nabla}(\\boldsymbol{u}),\n",
+    "$$\n",
+    "where $\\boldsymbol{u}$ is the displacement field presented in tutorial #6. Injecting this decomposition into the **Euler-Lagrange** strain expression above yields:\n",
+    "$$\n",
+    "\\boldsymbol{E} = \\boldsymbol{\\varepsilon} + \\frac{1}{2}\\boldsymbol{\\nabla}(\\boldsymbol{u})^T\\cdot\\boldsymbol{\\nabla}(\\boldsymbol{u}),\n",
+    "$$\n",
+    "with:\n",
+    "$$\n",
+    "\\boldsymbol{\\varepsilon} = \\frac{1}{2} \\big(\\boldsymbol{\\nabla}(\\boldsymbol{u})^T + \\boldsymbol{\\nabla}(\\boldsymbol{u}) \\big)\n",
+    "$$\n",
+    "the linear strain measure.\n",
+    "\n",
+    "**Morality:** From the quick developments above we see that in linear elasticity, the quadradic term $\\boldsymbol{\\nabla}(\\boldsymbol{u})^T\\cdot\\boldsymbol{\\nabla}(\\boldsymbol{u})$ is missing. In the following excercice, we are going to assess the importance of this term by comparing the response of a linear elastic and an hyperelastic model to the same loading."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "heading_collapsed": true
+   },
+   "source": [
+    "### Domain definition\n",
+    "Let's consider a simple circular patch."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "hidden": true
+   },
+   "outputs": [],
+   "source": [
+    "from bvpy.domains import Disk\n",
+    "from bvpy.utils.visu import plot, set_renderer\n",
+    "\n",
+    "patch = Disk(radius=2, cell_size= .1, clear=True, dim=2)\n",
+    "\n",
+    "set_renderer('notebook')\n",
+    "plot(patch)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "heading_collapsed": true
+   },
+   "source": [
+    "### Boundary conditions\n",
+    "Let's assume that this domain is deformed by an imposed displacement field appliedto its border."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "hidden": true
+   },
+   "outputs": [],
+   "source": [
+    "from bvpy.boundary_conditions import NormalDirichlet\n",
+    "\n",
+    "prescribed_boundary_displacement = 1\n",
+    "\n",
+    "traction = NormalDirichlet(val=prescribed_boundary_displacement, boundary='all')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Two different elastic responses\n",
+    "Let's know consider to types of elastic responses:\n",
+    "* A *classic* linear response as the one depicted in the previous tutorial.\n",
+    "* A non-linear response, often called hyperelasticity.\n",
+    "\n",
+    "While the former response will be implemented with the `LinearElasticForm` already presented in tutorial 6; the latter will be rely on a dedicated class `HyperElasticForm` from the `bvpy.vforms.elasticity` sub-module.\n",
+    "\n",
+    "Both will use the same elastic rheological parameters: Young's modulus and Poisson's ratio, with the same values."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Young = 1e3\n",
+    "Poisson = .5"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Linear elastic response"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from bvpy.vforms import LinearElasticForm\n",
+    "\n",
+    "linear_response = LinearElasticForm(young=Young, poisson=Poisson, source=[0., 0.], plane_stress=True)\n",
+    "\n",
+    "linear_response.info()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Hyper elastic response"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from bvpy.vforms import HyperElasticForm\n",
+    "\n",
+    "non_linear_response = HyperElasticForm(young=Young, poisson=Poisson, source=[0., 0.], plane_stress=True)\n",
+    "\n",
+    "non_linear_response.info()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Note:** Hyperelasticiy is a general term that encompasses various rheological models. In the current version of the library, we implemented *only* two of them: The **Saint-Venant Kirchhoff** model and the **Neo-Hookean** model. Selection between is controled by the `material_model` argument that sets the `self.model` attribute of the class. By default the chosen model is the Saint-Venant Kirchhoff one. Other models could of course be either added to the classes or could even be implemented in dedicated classes deriving from this one."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Models comparison"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The idea now is to compare the strain values computed by these two models for various displacements prescribed on the boundary of the structure.\n",
+    "\n",
+    "To that end, we are going to interate over prescribed displacement values, compute the corresponding strain fields for both models and compare them.\n",
+    "\n",
+    "> **Note:** We are going to use the domain `patch` previously defined as well as the two elastic *vform* classes `linear_response` and `non_linear_response`, we therefore do not need to re-implement them."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Running the simulations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "from bvpy import BVP\n",
+    "\n",
+    "from bvpy.utils.tensor import intensity\n",
+    "from bvpy.utils.post_processing import SolutionExplorer\n",
+    "\n",
+    "boundary_displacements = np.arange(.05, 1.05, .05)\n",
+    "\n",
+    "strn_int = {'ln': [], 'nl': []}\n",
+    "    \n",
+    "for bnd_displ in boundary_displacements:\n",
+    "    \n",
+    "    # -- Defining the boundary loading conditions\n",
+    "    traction = NormalDirichlet(val=bnd_displ, boundary='near(x*x + y*y, 4, 1e-2)')\n",
+    "    \n",
+    "    for name, response in {'ln': linear_response, 'nl': non_linear_response}.items():\n",
+    "        # -- Setting the problem\n",
+    "        stretching = BVP(patch, response, traction)\n",
+    "    \n",
+    "        # -- Solving the problems\n",
+    "        stretching.solve()\n",
+    "        \n",
+    "        # -- Computing the corresponding strain field\n",
+    "        fe_strain = response.strain(stretching.solution)\n",
+    "        np_strain = SolutionExplorer(fe_strain).to_vertex_values()\n",
+    "        \n",
+    "        strain_intensity = np.array(list(map(lambda x: intensity(x), np_strain)))\n",
+    "        \n",
+    "        strn_int[name].append((strain_intensity.mean(), strain_intensity.std()))     "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Note:** We used the `SolutionExplorer` class to convert the strain field initially defined as a `Fenics.Function` (`fe_strain` in the code above) into a `numpy.nDarray` version (`np_strain` above). The idea behind this move is to facilitate the analysis of the simulation outputs by providing them in the very common format of `numpy.nDarray`."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Quick visualization of the results"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Since the models are completely isotropic (and homogeneous) is shape, loading and properties, we expect the resulting strain fields to be also isotropic and homogeneous. We aim therefore to compare their intensity values averaged over the whole domain.\n",
+    "\n",
+    "To that end, for each value of the prescribed boundary displacement, we are going to compute the average value of the *linear* and the *non-linear* strain field intensity. We will plot one againts the other and see how the resulting curve shifts away from the first bisector."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.lines as mlines\n",
+    "\n",
+    "threshold = .05 # Comparison threshold\n",
+    "\n",
+    "ln_strain = np.array(strn_int['ln'])\n",
+    "nl_strain = np.array(strn_int['nl'])\n",
+    "\n",
+    "patch_radius = patch._characteristic_size['Radius']\n",
+    "relat_prescr_displ = []\n",
+    "\n",
+    "fig = plt.figure(figsize=(15,10))\n",
+    "plt.plot([0, 60], [0, 60], ls=':',c='k', alpha=.5)\n",
+    "\n",
+    "for bnd_displ, ln_strn, nl_strn in zip(boundary_displacements, ln_strain, nl_strain):\n",
+    "\n",
+    "    if nl_strn[0]/ln_strn[0]-1 < threshold:\n",
+    "        marker, color = 'o', 'g'\n",
+    "    else:\n",
+    "        marker, color = '^','r'\n",
+    "        relat_prescr_displ.append([ln_strn[0], bnd_displ/patch_radius])\n",
+    "\n",
+    "    plt.scatter(ln_strn[0]*100, nl_strn[0]*100, c=color, marker=marker)\n",
+    "\n",
+    "plt.plot([relat_prescr_displ[0][0]*100, relat_prescr_displ[0][0]*100], [0, 20], c='k', alpha=.5)\n",
+    "\n",
+    "# -- legend formatting\n",
+    "red_triangles = mlines.Line2D([], [], ls='', color='r', marker='^',\n",
+    "                              label=f\"Above {threshold:.0%} discrepency\")\n",
+    "green_dots = mlines.Line2D([], [], ls='', color='g', marker='o',\n",
+    "                           label=f\"Below {threshold:.0%} discrepency\")\n",
+    "gray_line = mlines.Line2D([], [], color='k', alpha=.5,\n",
+    "                          label=f\"Boundary displacement = {relat_prescr_displ[0][1]:.0%}\\n of initial radius\")\n",
+    "plt.legend(handles=[red_triangles, green_dots, gray_line], loc=\"lower right\", fontsize=16)\n",
+    "\n",
+    "plt.xlabel('Linear elasticity strain intensity (in %)', fontsize=24)\n",
+    "plt.ylabel('Hyper elasticity strain intensity (in %)', fontsize=24)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Morality:** Hyperelastic models and non-linear strain measures (such as the Green-Lagrange one we used here) are usually considered when the studied displacement field is large. A rule of thumb says that when the strain values exceed 10%, one would better use an hyperelastic model to describe its system. The figure above illustrate this by showing that above 10% the strain intensities of both models deviate for more than 5%. The quadratic term ($\\boldsymbol{\\nabla}(\\boldsymbol{u})^T\\cdot\\boldsymbol{\\nabla}(\\boldsymbol{u})$) discharded in the linear approximation is not negligeable anymore.\n",
+    "\n",
+    "> **Note:** Beside this illustration of basic continuum mechanics facts, let's remark that the implementation of the simulation loop generating the data, within the **BVPy** formalism, is smaller than the code snippet required to visualize the data (resp. 30 and 36 lines). ;)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Toward Multiphysics modeling: Combining BVPs and IBVPs into \"mixted models\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The goal now is to illustrate how we can use the **BVPy** library to combine *simple* *atomic* models into more complex ones. To that end, we are going to consider two models presented in previous tutorials:\n",
+    "* The linear elastic model of a 2D curved domain loaded with a normal force field, introduced in tutorial #5.\n",
+    "* The *Turing-like* Lengyel-Epstein model, introduced in tutorial #6.\n",
+    "\n",
+    "The resulting *toy system* we want to build up will feature the following properties:\n",
+    "* The domain is the surface of a sphere.\n",
+    "* The Young modulus quantifying the elasticity of the domain is function of the concentration of one of the two chemical species forming the reaction-diffusion model.\n",
+    "* From a mechanical perspective, the system is assumed to be in a quasi-static evolution; *i.e.* chemical reacitons and diffusion processes are far slower than mechanical relaxation.\n",
+    "\n",
+    "The latter point means that, in terms of modeling strategy we are going to discretize time into steps upon which chemicals will evolve according to an *IBVP*. In parallel, the mechanical state of the structure will be computed as a *BVP* depending on the concentrations computed at the previous time."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Domain definition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's consider a sphere."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from bvpy.domains import Sphere\n",
+    "\n",
+    "globe = Sphere(radius=3, cell_size=.1, clear=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Definition of the mechanics of the system"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The particularity of the considered system is that the Young modulus of the material is a function of a chemical compound concentration. The first thing to do is then to write the corresponding function that takes a concentration as an inputs and returns the corresponding Young's modulus value."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "code_folding": [
+     2
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import fenics as fe\n",
+    "\n",
+    "def compute_young(concentrations, Young_iso=1e3, domain=globe, component=0):\n",
+    "    \"\"\"Computes Young's modulus given a concentration field.\n",
+    "    \n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    concentrations : Fenics.UserExpression or Fenics.Function.\n",
+    "        The two concentrations of the interacting compounds from the Turing model.\n",
+    "    Young_iso : float.\n",
+    "        The basal value of the Young modulus without influence of the chemical field. \n",
+    "        Optional, the default value is 1e3.\n",
+    "    domain : Bvpy.Domains.\n",
+    "        The domain to consider.\n",
+    "        Optional, the default is globe (defined above).\n",
+    "    component : int.\n",
+    "        The index of the chemical compound to consider (either 0 or 1).\n",
+    "        Optional, the default is 0.\n",
+    "    \n",
+    "    Returns\n",
+    "    -------\n",
+    "    Young : Fenics.Function\n",
+    "        The scalar field of the Young modulus to consider in the linear elastic model.\n",
+    "    \n",
+    "    \"\"\"\n",
+    "    \n",
+    "    vct_func_space = fe.VectorFunctionSpace(domain.mesh, \"P\", 1, dim=2)\n",
+    "    \n",
+    "    if isinstance(concentrations, fe.UserExpression):\n",
+    "        c_func = fe.interpolate(concentrations, vct_func_space).sub(component)\n",
+    "    \n",
+    "    elif isinstance(concentrations, fe.Function):\n",
+    "        c_func = fe.project(concentrations, vct_func_space).sub(component)\n",
+    "    \n",
+    "    else:\n",
+    "        print(\"WARNING: the type of the input is not supported\")\n",
+    "    \n",
+    "    c_max = c_func.vector().max()\n",
+    "    c_min = c_func.vector().min()\n",
+    "\n",
+    "    Young = fe.Constant(1)\n",
+    "    Young += fe.Constant(100 / (c_max - c_min)) * (fe.Constant(c_max) - c_func)\n",
+    "    Young *= fe.Constant(Young_iso)\n",
+    "    \n",
+    "    return Young"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Notes:** \n",
+    "> * We are considering two types of input for the function `compute_young`, either `Fenics.UserExpression` or `Fencis.Function`. This comes from the fact that initial conditions for the *Turing-like* model are defined with *user expression* while the model itself ouputs *functions*.\n",
+    "> * The `concentrations` variable contains the concentrations of the two chemical species. One can choose one or the other as the *stiffening agent* by setting the `component` variable either to 0 or 1.\n",
+    "\n",
+    "**Inter-operability with FEniCS:**\n",
+    "The function above is a good example of how **FEniCS** objects and methods can easily be used within **BVPy**. Such functions can be implemented by users already familiar with **FEniCS** and used by everyone."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's now instanciate the mechanical model.\n",
+    "\n",
+    "To that end, we are going to define a function that takes concentrations as inputs and return a stress field as outputs. This stress field will correspond to the distribution of mechanical constrains within the considered domain, when mechanical equilibrium between pressure forces and elastic response is reached. The elastic response, parametrized by the Young modulus of the structure is dependent upon the concentration field via the function defined above."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "code_folding": [
+     3
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from bvpy.domains.geometry import normal\n",
+    "from bvpy.boundary_conditions import dirichlet\n",
+    "\n",
+    "def compute_mecha_equi(concentrations, **kwargs):\n",
+    "    \"\"\"Simulates the loading of an linear elastic domain by a normal force field.\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    concentrations : Fenics.UserExpression or Fenics.Function.\n",
+    "        The two concentrations of the interacting compounds from the Turing model.\n",
+    "    \n",
+    "    Other parameters\n",
+    "    ----------------\n",
+    "    domain : bvpy.Domain\n",
+    "        Optional. The default is globe.\n",
+    "    pressure : float\n",
+    "        The value of the pressure loading the structure.\n",
+    "        Optional. The default is 1e2.\n",
+    "    Young_iso : float\n",
+    "        Basal value of the structure Young's modulus.\n",
+    "        Optional. The default is 1e3.\n",
+    "    Poisson : float\n",
+    "        Poisson's coefficient value for the structure.\n",
+    "        Optional. The default is .5\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    Fenics.Function\n",
+    "        The stress tensor field at mechanical equilibrium.\n",
+    "\n",
+    "    \"\"\"\n",
+    "    # -- Parameters\n",
+    "    domain = kwargs.get('domain', globe)\n",
+    "    pressure = kwargs.get('pressure', 1e3)\n",
+    "    Young_iso = kwargs.get('Young_iso', 1e2)\n",
+    "    Poisson = kwargs.get('Poisson', .5)\n",
+    "\n",
+    "    # -- Setting boundary conditions\n",
+    "    fixed_center = [dirichlet(0, boundary='near(x, 0, 1e-2)', subspace=0),\n",
+    "                    dirichlet(0, boundary='near(y, 0, 1e-2)', subspace=1),\n",
+    "                    dirichlet(0, boundary='near(z, 0, 1e-2)', subspace=2)]\n",
+    "\n",
+    "    # -- Setting the loading force field as a source term within the vform\n",
+    "    pressure_force = normal(domain, scale=pressure)\n",
+    "    \n",
+    "    # -- Setting the elastic response\n",
+    "    Young = compute_young(concentrations, Young_iso=Young_iso, domain=domain)\n",
+    "    \n",
+    "    elastic_response = LinearElasticForm(young=Young, poisson=Poisson, source=pressure_force,\n",
+    "                                         plane_stress=True)\n",
+    "\n",
+    "    # -- Setting the problem\n",
+    "    mecha_equi = BVP(domain, elastic_response, fixed_center)\n",
+    "\n",
+    "    # -- Solving the problem\n",
+    "    mecha_equi.solve()\n",
+    "    \n",
+    "    displacement = mecha_equi.solution\n",
+    "    \n",
+    "    return elastic_response.stress(displacement)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Definition of the chemistry of  the system"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's now define the same kind of function but to handle the chemical side of the system."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "code_folding": [
+     5
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from bvpy.ibvp import IBVP\n",
+    "from bvpy.vforms import CoupledTransportForm\n",
+    "\n",
+    "def compute_turing_pattern(initial_concentrations, **kwargs):\n",
+    "    \"\"\"Computes the time evolution of two chemicals interacting with one another in a Turing-like manner.\n",
+    "    \n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    concentrations : Fenics.UserExpression or Fenics.Function.\n",
+    "        The two concentrations of the interacting compounds from the Turing model.\n",
+    "    \n",
+    "    Other parameters\n",
+    "    ----------------\n",
+    "    domain : bvpy.Domain\n",
+    "        Optional. The default is globe.\n",
+    "    Da : float\n",
+    "        The value of diffusion coefficient of the first chemical compound (named \"a\").\n",
+    "        Optional. The default is 1.\n",
+    "    ka : float\n",
+    "        The value of production rate of the first chemical compound (named \"a\").\n",
+    "        Optional. The default is 9.\n",
+    "    Db : float\n",
+    "        The value of diffusion coefficient of the second chemical compound (named \"b\").\n",
+    "        Optional. The default is 2e-2.\n",
+    "    ka : float\n",
+    "        The value of production rate of the second chemical compound (named \"b\").\n",
+    "        Optional. The default is 11.\n",
+    "    t_min : float\n",
+    "        Initial time of the simulation.\n",
+    "        Optional. The default is 0.\n",
+    "    t_max : float\n",
+    "        Final time of the simulation.\n",
+    "        Optional. The default is 2.5.\n",
+    "    dt : float\n",
+    "        Time step increment.\n",
+    "        Optional. The default is .1.\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    Fenics.Function\n",
+    "        The concentrations (a, b) vector field after a t_max - t_min time.\n",
+    "\n",
+    "    \"\"\"\n",
+    "    # -- Parameters\n",
+    "    domain = kwargs.get('domain', globe)\n",
+    "    Da = kwargs.get('Da', 1)\n",
+    "    Db = kwargs.get('Db', 2e-2)\n",
+    "    ka = kwargs.get('ka', 9)\n",
+    "    kb = kwargs.get('kb', 11)\n",
+    "    \n",
+    "    t_min = kwargs.get('t_min', 0)\n",
+    "    t_max = kwargs.get('t_max',2.5)\n",
+    "    dt = kwargs.get('dt',.1)\n",
+    "    \n",
+    "    save_path = kwargs.get('saving_path', '../data/tutorial_results/tuto_7_result_1.xdmf')\n",
+    "    \n",
+    "    # -- Instancing Turing model\n",
+    "    LE_model = CoupledTransportForm()\n",
+    "\n",
+    "    LE_model.add_species('a', Da, lambda a, b: ka*(b-(a*b)/(1+b*b)))\n",
+    "    LE_model.add_species('b', Db, lambda a, b: kb-b-(4*a*b)/(1+b*b))\n",
+    "\n",
+    "    # -- Building the problem\n",
+    "    turing_system = IBVP(domain, LE_model, initial_state=initial_concentrations)\n",
+    "\n",
+    "    # -- Solving the problem\n",
+    "    turing_system.integrate(t_min, t_max, dt, save_path)\n",
+    "    \n",
+    "    return turing_system.solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Let's combine both models into a big loop"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First, we need to define an initial state in term of chemical concentrations."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from bvpy.utils.pre_processing import create_expression\n",
+    "\n",
+    "# -- Creating the initial concentration distribution\n",
+    "ka = 9\n",
+    "kb = 11\n",
+    "\n",
+    "np.random.seed(1093)\n",
+    "globe.discretize()\n",
+    "random = np.random.uniform(low=-1, high=1, size=globe.mesh.num_entities(2))\n",
+    "\n",
+    "random_pattern = [[1+0.04*ka**2+0.1*r, 0.2*kb+0.1*r] for r in random]\n",
+    "\n",
+    "concentrations = create_expression(random_pattern)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Notes:** \n",
+    "> * The `create_expression` function, from the `bvpy.utils.pre_processing` sub-module is used to *translate* a *classical* `list` into a `Fenics.UserExpression`. \n",
+    "> * Similarly, the `SolutionExplorer` class from the `bvpy.utils.post_processing` sub-module enables *translation* of `Fenics.Function` into more common objects, such as `numpy.nDarray`.\n",
+    "> * The general idea is that the `pre_processing` and `post_processing` sub-modules aim at interfacing the ***FEniCS***-*based* core with more commun data representations."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we can combine both models together. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from bvpy.utils.io import save\n",
+    "\n",
+    "# -- Running the simulations\n",
+    "step_nbr = 1\n",
+    "for step in np.arange(1, 1+step_nbr):\n",
+    "    \n",
+    "    print(f\"\\n step #{step} over {step_nbr}\")\n",
+    "        \n",
+    "    stress = compute_mecha_equi(concentrations)\n",
+    "    concentrations = compute_turing_pattern(concentrations)\n",
+    "\n",
+    "    save(stress, f\"../data/tutorial_results/tuto_7_result_{1+step}.xdmf\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> **Note:** For the sake of computational time (especially in the CI context), we only considered 1 time step in the above loop. We encourage the interested user to download the notebook and run the simulation on more steps. Below we show an animation of the simulation outputs computed over 10 steps."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's have a look at the resulting stress field and its quasi-static evolution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from IPython.display import Image\n",
+    "\n",
+    "Image(filename=\"../data/tutorial_results/tuto_7_result_3_paraview_snapshot.gif.png\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> The animation above, generated from the **paraview** software, shows the stress field quasi-static time evolution across the ten steps of the simulation. The colors code for its intensity (blue = low and red = high). The ellopsoid depicts its local orientation and anisotropy (their radii are aligned in the eigenvectors directions and proportional to the corresponding eigenvalues)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Concluding remark\n",
+    "With this simple *artificial* example we wanted to show how the main `bvpy.bvp` and `bvpy.ibvp` classes can be combined into complex models.\n",
+    "\n",
+    "Models inspired by real life systems can sometimes feature interplay between fields of various natures, *e.g.* scalar, vector or tensor; each featuring its own dynamics happening putatively at different time scales.\n",
+    "This example illustrates how one can build such a complex system in a rather intuitive manner within the **BVPy** framework.\n",
+    "Moreover, the fact that **FEniCS** objects and methods are directly accessible enables fine tuning and manipulations of these models."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.9"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": true,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  },
+  "varInspector": {
+   "cols": {
+    "lenName": 16,
+    "lenType": 16,
+    "lenVar": 40
+   },
+   "kernels_config": {
+    "python": {
+     "delete_cmd_postfix": "",
+     "delete_cmd_prefix": "del ",
+     "library": "var_list.py",
+     "varRefreshCmd": "print(var_dic_list())"
+    },
+    "r": {
+     "delete_cmd_postfix": ") ",
+     "delete_cmd_prefix": "rm(",
+     "library": "var_list.r",
+     "varRefreshCmd": "cat(var_dic_list()) "
+    }
+   },
+   "types_to_exclude": [
+    "module",
+    "function",
+    "builtin_function_or_method",
+    "instance",
+    "_Feature"
+   ],
+   "window_display": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/tutorials/bvpy_tutorial_plasticity.ipynb b/tutorials/bvpy_tutorial_plasticity.ipynb
index de68833309315452106e4c19f5ec7d098e647419..64883d5d91d50586fb67046a61a847c3244fbddc 100644
--- a/tutorials/bvpy_tutorial_plasticity.ipynb
+++ b/tutorials/bvpy_tutorial_plasticity.ipynb
@@ -164,7 +164,38 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.7.9"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
   },
   "varInspector": {
    "cols": {