diff --git a/.coveragerc b/.coveragerc
deleted file mode 100644
index a44f260e59a587768f969fb3cce0b14d5af6bad6..0000000000000000000000000000000000000000
--- a/.coveragerc
+++ /dev/null
@@ -1,17 +0,0 @@
-# {# pkglts, coverage
-[html]
-title = bvpy's coverage
-directory = build/htmlcov
-
-[run]
-source = bvpy
-
-[report]
-show_missing = True
-exclude_lines =
-    # Don't complain if tests don't hit defensive assertion code:
-    raise AssertionError
-    raise NotImplemented
-    raise NotImplementedError
-    __main__
-# #}
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 42974364977d8e882270a11ba42eb2bfcb9ea301..1d02cda35e0dd707a7f2326f2de7955feb36b937 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -13,7 +13,7 @@ conda_build:
     - large
   script:
     - mkdir conda-bld
-    - conda build . -c conda-forge --output-folder conda-bld/
+    - conda build conda/ -c conda-forge --output-folder conda-bld/
   retry:
     max: 2
     when: runner_system_failure
@@ -21,9 +21,15 @@ conda_build:
     expire_in: 1 hour
     paths:
       - conda-bld
-  only:
-    - develop
-    - master
+  rules:
+    - if: $CI_COMMIT_BRANCH =~ /^release/
+      when: manual
+      allow_failure: true
+    - if: '$CI_COMMIT_BRANCH == "develop"'
+      when: manual
+      allow_failure: true
+    - if: '$CI_COMMIT_BRANCH == "master"'
+      when: on_success
 
 
 coverage:
@@ -36,13 +42,13 @@ coverage:
     # Install required 'libgl1' system dependency:
     - apt-get update && apt-get install libgl1 -y
     # Install faster solver 'conda-libmamba-solver' & set it as default solver:
-    - conda install -y -n base -c conda-forge conda-libmamba-solver
-    - conda config --set solver libmamba
+    # - conda install -y -n base -c conda-forge conda-libmamba-solver
+    # - conda config --set solver libmamba
     # Create 'bvpy' conda environment & activate it:
-    - conda env create -f conda/env_20240207.yaml -y
+    - conda env create --file conda/env.yaml
     - source activate bvpy
     # Install `bvpy` from sources:
-    - python -m pip install -e .
+    - python -m pip install -e .[test]
     # Test it with `pytest`:
     - pytest --version
     - pytest -v --cov=bvpy test/
@@ -52,6 +58,7 @@ coverage:
   only:
     - develop
     - master
+    - /^release/
 
 anaconda :
   stage: deploy
@@ -70,9 +77,15 @@ anaconda :
   dependencies:
     - conda_build
     - coverage
-  only:
-    - master
-    - develop
+  rules:
+    - if: $CI_COMMIT_BRANCH =~ /^release/
+      when: manual
+      allow_failure: true
+    - if: '$CI_COMMIT_BRANCH == "develop"'
+      when: manual
+      allow_failure: true
+    - if: '$CI_COMMIT_BRANCH == "master"'
+      when: on_success
 
 pages:
   stage: deploy
@@ -82,19 +95,20 @@ pages:
   image: continuumio/miniconda3
   script:
     # Install required 'libgl1' system dependency:
-    - apt-get update && apt-get install libgl1 -y
+    - apt-get update && apt-get install libgl1 wget -y
+    # Download and install pandoc3.1.12.3 (version must be at least (2.14.2) but less than (4.0.0)):
+    - wget https://github.com/jgm/pandoc/releases/download/3.1.12.3/pandoc-3.1.12.3-1-amd64.deb
+    - dpkg -i pandoc-3.1.12.3-1-amd64.deb
     # Install faster solver 'conda-libmamba-solver' & set it as default solver:
-    - conda install -y -n base -c conda-forge conda-libmamba-solver
-    - conda config --set solver libmamba
+    # - conda install -y -n base -c conda-forge conda-libmamba-solver
+    # - conda config --set solver libmamba
     # Create 'bvpy' conda environment & activate it:
-    - conda env create -f conda/env_20240207.yaml -y
-    - conda activate bvpy
+    - conda env create -f conda/env.yaml -n bvpy
+    - source activate bvpy
     # Install `bvpy` package (from sources):
-    - python -m pip install -e .
-    # Build the documentation:
-    - cd doc && make clean && make html
-    # Publish the documentation:
-    - mv build/html ../public
+    - python -m pip install -e .[doc]
+    # Build the documentation with Sphinx:
+    - sphinx-build -b html doc/ public/
   retry:
     max: 2
     when: runner_system_failure
@@ -104,10 +118,14 @@ pages:
       - public
   dependencies: []
   rules:
-      - if: '$CI_COMMIT_BRANCH == "develop"'
-        when: manual
-      - if: '$CI_COMMIT_BRANCH == "master"'
-        when: on_success
+    - if: $CI_COMMIT_BRANCH =~ /^release/
+      when: manual
+      allow_failure: true
+    - if: '$CI_COMMIT_BRANCH == "develop"'
+      when: manual
+      allow_failure: true
+    - if: '$CI_COMMIT_BRANCH == "master"'
+      when: on_success
 
 docker_build_deploy:
   image: docker:latest
diff --git a/README.rst b/README.rst
index f8a4b83650b2f1a39fa5a3867154d1dac2068e22..40f2f672db46f573a402a41a58688abaf990210c 100644
--- a/README.rst
+++ b/README.rst
@@ -45,9 +45,9 @@ Some tutorials explaining basic manipulations are gathered `here <https://mosaic
 Requirements
 ------------
 
-* Python 3.7+
-* FEniCS
-* GMSH
+* Python 3.9
+* FEniCS 2019.1.0
+* GMSH 4.11
 
 
 Installation
@@ -61,19 +61,15 @@ You will need conda in order to install ``bvpy``, you can download it  `here <ht
 
   git clone https://gitlab.inria.fr/mosaic/bvpy.git
   cd bvpy
-  conda env create -f conda/env_OS.yaml -n bvpy-dev
+  conda create -n bvpy-dev -c conda-forge "python=3.9" "fenics=2019.1.0"
   conda activate bvpy-dev
-  python setup.py develop --prefix=$CONDA_PREFIX
-
-**Note:** depending on your operating system (linux or macOs) you will need to use the dedicated environment file (`env_OS.yaml`) on the third line in the code block above:
-* For linux use `env_linux.yaml`
-* For macOs use `env_osx.yaml`
+  python -m pip install -e .
 
 - From anaconda:
 
 .. code-block:: bash
 
-	conda install -c conda-forge bvpy
+	conda install -c mosaic bvpy
 
 
 **Important Note:** Bvpy is currently running on *FEniCS legacy* and not *FEniCSx*. If you are using an ARM-based computer (mac M1s, M2s), there is currently no ARM-Friendly conda package of *FEniCS legacy*. But you can still run x86-based packages. to do so type the following:
@@ -107,17 +103,17 @@ And use the URL on the terminal to open the tutorials.
 Test
 ----
 
-If you install ``bvpy`` from anaconda install ``pytest`` and ``pytest-cov`` inside your environment:
+If you install ``bvpy`` from sources install ``pytest`` and ``pytest-cov`` inside your environment:
 
 .. code-block:: bash
 
-	conda install -c conda-forge pytest pytest-cov
+	python -m pip install -e ".[test]"
 
 And then run at the root of the project:
 
 .. code-block:: bash
 
-	pytest -v
+	pytest -v --cov=bvpy test/
 
 Support
 -------
diff --git a/conda/build.sh b/conda/build.sh
deleted file mode 100644
index f505a90aa29c7db198093e72eec28e7af3219ecc..0000000000000000000000000000000000000000
--- a/conda/build.sh
+++ /dev/null
@@ -1,6 +0,0 @@
-#!/bin/bash
-
-export CONDA_BUILD=1
-
-$PYTHON -m pip install .
-rm -rf ${PREFIX}/include
diff --git a/conda/bvpy_x86.yml b/conda/bvpy_x86.yml
deleted file mode 100755
index 2c70ef150e358149ccbab601fb2aa6e3e255c0b9..0000000000000000000000000000000000000000
--- a/conda/bvpy_x86.yml
+++ /dev/null
@@ -1,182 +0,0 @@
-name: bvpy_x86
-channels:
-  - conda-forge
-  - defaults
-dependencies:
-  - appnope=0.1.3=pyhd8ed1ab_0
-  - argon2-cffi=21.3.0=pyhd8ed1ab_0
-  - argon2-cffi-bindings=21.2.0=py37h69ee0a8_2
-  - attrs=21.4.0=pyhd8ed1ab_0
-  - backcall=0.2.0=pyh9f0ad1d_0
-  - backports=1.0=py_2
-  - backports.functools_lru_cache=1.6.4=pyhd8ed1ab_0
-  - beautifulsoup4=4.11.1=pyha770c72_0
-  - bleach=5.0.0=pyhd8ed1ab_0
-  - boost-cpp=1.74.0=h8b082ac_8
-  # - bvpy=1.0.1=py37h694c41f_0
-  - bzip2=1.0.8=h0d85af4_4
-  - c-ares=1.18.1=h0d85af4_0
-  - ca-certificates=2022.5.18.1=h033912b_0
-  - cached-property=1.5.2=hd8ed1ab_1
-  - cached_property=1.5.2=pyha770c72_1
-  - cffi=1.15.0=py37h446072c_0
-  - cftime=1.6.0=py37h49e79e5_1
-  - cmake=3.23.1=h35a7dd9_0
-  - commonmark=0.9.1=py_0
-  - curl=7.83.1=h23f1065_0
-  - dataclasses=0.8=pyhc8e2a94_3
-  - debugpy=1.6.0=py37h0582d14_0
-  - decorator=5.1.1=pyhd8ed1ab_0
-  - defusedxml=0.7.1=pyhd8ed1ab_0
-  - eigen=3.4.0=h940c156_0
-  - entrypoints=0.4=pyhd8ed1ab_0
-  - expat=2.4.8=h96cf925_0
-  - fenics=2019.1.0=py37hf985489_30
-  - fenics-dijitso=2019.1.0=py37hf985489_30
-  - fenics-dolfin=2019.1.0=py37hc1c2bcd_30
-  - fenics-ffc=2019.1.0=py37hf985489_30
-  - fenics-fiat=2019.1.0=py37hf985489_30
-  - fenics-libdolfin=2019.1.0=h5e46fef_30
-  - fenics-ufl=2019.1.0=py37hf985489_30
-  - fftw=3.3.10=mpi_openmpi_hcb0582e_2
-  - flit-core=3.7.1=pyhd8ed1ab_0
-  - fontconfig=2.13.94=h10f422b_0
-  - freetype=2.10.4=h4cff582_1
-  - future=0.18.2=py37hf985489_5
-  - gmp=6.2.1=h2e338ed_0
-  - gmpy2=2.1.2=py37h60f582e_0
-  - gmsh=4.5.6=h3eaa6b3_0
-  - h5py=3.6.0=nompi_py37h0ac0de7_100
-  - hdf4=4.2.15=hefd3b78_3
-  - hdf5=1.12.1=mpi_openmpi_h0d5d0c2_4
-  - hypre=2.24.0=mpi_openmpi_h001df97_0
-  - icu=70.1=h96cf925_0
-  - importlib-metadata=4.11.4=py37hf985489_0
-  - importlib_metadata=4.11.4=hd8ed1ab_0
-  - importlib_resources=5.7.1=pyhd8ed1ab_1
-  - ipykernel=6.13.0=py37h0a7177a_0
-  - ipython=7.33.0=py37hf985489_0
-  - ipython_genutils=0.2.0=py_1
-  - ipywidgets=7.7.0=pyhd8ed1ab_0
-  - jedi=0.18.1=py37hf985489_1
-  - jinja2=3.1.2=pyhd8ed1ab_0
-  - jpeg=9e=h5eb16cf_1
-  - jsonschema=4.5.1=pyhd8ed1ab_0
-  - jupyter_client=7.3.1=pyhd8ed1ab_0
-  - jupyter_core=4.10.0=py37hf985489_0
-  - jupyterlab_pygments=0.2.2=pyhd8ed1ab_0
-  - jupyterlab_widgets=1.1.0=pyhd8ed1ab_0
-  - krb5=1.19.3=hb98e516_0
-  - libblas=3.9.0=14_osx64_openblas
-  - libcblas=3.9.0=14_osx64_openblas
-  - libcurl=7.83.1=h23f1065_0
-  - libcxx=14.0.3=hc203e6f_0
-  - libedit=3.1.20191231=h0678c8f_2
-  - libev=4.33=haf1e3a3_1
-  - libffi=3.4.2=h0d85af4_5
-  - libgfortran=5.0.0=9_3_0_h6c81a4c_23
-  - libgfortran5=9.3.0=h6c81a4c_23
-  - libiconv=1.16=haf1e3a3_0
-  - liblapack=3.9.0=14_osx64_openblas
-  - libnetcdf=4.8.1=mpi_openmpi_h6e52f38_2
-  - libnghttp2=1.47.0=hca56917_0
-  - libopenblas=0.3.20=openmp_hb3cd9ec_0
-  - libpng=1.6.37=h7cec526_2
-  - libsodium=1.0.18=hbcb3906_1
-  - libssh2=1.10.0=hd3787cc_2
-  - libuv=1.43.0=h0d85af4_0
-  - libxml2=2.9.14=h08a9926_0
-  - libzip=1.8.0=h7e5727d_1
-  - libzlib=1.2.11=h6c3fc93_1014
-  - llvm-openmp=14.0.3=ha654fa7_0
-  - lz4-c=1.9.3=he49afe7_1
-  - markupsafe=2.1.1=py37h69ee0a8_1
-  - matplotlib-inline=0.1.3=pyhd8ed1ab_0
-  - meshio=5.3.4=pyhd8ed1ab_0
-  - metis=5.1.0=h2e338ed_1006
-  - mistune=0.8.4=py37h271585c_1005
-  - mpc=1.2.1=hbb51d92_0
-  - mpfr=4.1.0=h0f52abe_1
-  - mpi=1.0=openmpi
-  - mpi4py=3.1.3=py37hadad684_1
-  - mpmath=1.2.1=pyhd8ed1ab_0
-  - mumps-include=5.2.1=h694c41f_11
-  - mumps-mpi=5.2.1=had0ebf5_11
-  - nbclient=0.6.3=pyhd8ed1ab_0
-  - nbconvert=6.5.0=pyhd8ed1ab_0
-  - nbconvert-core=6.5.0=pyhd8ed1ab_0
-  - nbconvert-pandoc=6.5.0=pyhd8ed1ab_0
-  - nbformat=5.4.0=pyhd8ed1ab_0
-  - ncurses=6.3=h96cf925_1
-  - nest-asyncio=1.5.5=pyhd8ed1ab_0
-  - netcdf4=1.5.8=nompi_py37h5a78667_101
-  - nose=1.3.7=py_1006
-  - notebook=6.4.11=pyha770c72_0
-  - numpy=1.21.6=py37h345d48f_0
-  - occt=7.4.0=hb9b6dc7_3
-  - openmpi=4.1.3=hd33e60e_103
-  - openssl=3.0.3=hfe4f2af_0
-  - packaging=21.3=pyhd8ed1ab_0
-  - pandoc=2.18=h694c41f_0
-  - pandocfilters=1.5.0=pyhd8ed1ab_0
-  - parmetis=4.0.3=h7eda126_1005
-  - parso=0.8.3=pyhd8ed1ab_0
-  - petsc=3.17.1=real_hf50cc79_101
-  - petsc4py=3.17.1=real_hb74dba9_100
-  - pexpect=4.8.0=pyh9f0ad1d_2
-  - pickleshare=0.7.5=py_1003
-  - pip=22.1.1=pyhd8ed1ab_0
-  - pkg-config=0.29.2=ha3d46e9_1008
-  - pkgconfig=1.5.5=py37hf985489_2
-  - plotly=5.8.0=pyhd8ed1ab_0
-  - prometheus_client=0.14.1=pyhd8ed1ab_0
-  - prompt-toolkit=3.0.29=pyha770c72_0
-  - psutil=5.9.1=py37h994c40b_0
-  - ptscotch=6.0.9=h6295d7f_2
-  - ptyprocess=0.7.0=pyhd3deb0d_0
-  - pybind11=2.9.2=py37h18621fa_1
-  - pybind11-global=2.9.2=py37h18621fa_1
-  - pycparser=2.21=pyhd8ed1ab_0
-  - pygments=2.12.0=pyhd8ed1ab_0
-  - pyparsing=3.0.9=pyhd8ed1ab_0
-  - pyrsistent=0.18.1=py37h69ee0a8_1
-  - python=3.7.12=hf3644f1_100_cpython
-  - python-dateutil=2.8.2=pyhd8ed1ab_0
-  - python-fastjsonschema=2.15.3=pyhd8ed1ab_0
-  - python-gmsh=4.5.6=0
-  - python_abi=3.7=2_cp37m
-  - pyzmq=23.0.0=py37h573e605_0
-  - readline=8.1=h05e3726_0
-  - rhash=1.4.1=h0d85af4_0
-  - rich=12.4.4=pyhd8ed1ab_0
-  - scalapack=2.2.0=h208a4c8_1
-  - scotch=6.0.9=h3da7401_2
-  - send2trash=1.8.0=pyhd8ed1ab_0
-  - setuptools=62.3.2=py37hf985489_0
-  - six=1.16.0=pyh6c4a22f_0
-  - slepc=3.17.1=real_h0adec3d_100
-  - slepc4py=3.17.1=real_haf36324_100
-  - soupsieve=2.3.1=pyhd8ed1ab_0
-  - sqlite=3.38.5=hd9f0692_0
-  - suitesparse=5.10.1=h68a9093_0
-  - superlu=5.2.2=h621599f_0
-  - superlu_dist=7.2.0=hcde7739_0
-  - sympy=1.10.1=py37hf985489_0
-  - tbb=2020.2=h940c156_4
-  - tenacity=8.0.1=pyhd8ed1ab_0
-  - terminado=0.15.0=py37hf985489_0
-  - tinycss2=1.1.1=pyhd8ed1ab_0
-  - tk=8.6.12=h5dbffcc_0
-  - tornado=6.1=py37h69ee0a8_3
-  - traitlets=5.2.1.post0=pyhd8ed1ab_0
-  - typing_extensions=4.2.0=pyha770c72_1
-  - wcwidth=0.2.5=pyh9f0ad1d_2
-  - webencodings=0.5.1=py_1
-  - wheel=0.37.1=pyhd8ed1ab_0
-  - widgetsnbextension=3.6.0=py37hf985489_0
-  - xz=5.2.5=haf1e3a3_1
-  - yaml=0.2.5=h0d85af4_2
-  - zeromq=4.3.4=he49afe7_1
-  - zipp=3.8.0=pyhd8ed1ab_0
-  - zlib=1.2.11=h6c3fc93_1014
-  - zstd=1.5.2=ha9df2e0_1
diff --git a/conda/conda_build_config.yaml b/conda/conda_build_config.yaml
index 08ce6cfaf133e821f46287ab83844a6d13690711..40d9784da3c647dacace249147fe49c2ea0343ad 100644
--- a/conda/conda_build_config.yaml
+++ b/conda/conda_build_config.yaml
@@ -1,4 +1,4 @@
 python:
-    - 3.7
+    # - 3.7
     # - 3.8
-    # - 3.9
+    - 3.9
diff --git a/conda/env.yaml b/conda/env.yaml
index 0f1a3e48219b3b7868f2deecf2ed8ef8d9ed6a73..d4e06a84d3f0326e8a6866b2dd22ac60a3c47b6d 100644
--- a/conda/env.yaml
+++ b/conda/env.yaml
@@ -3,22 +3,6 @@ channels:
   - conda-forge
   - defaults
 dependencies:
-  - python=3.7
-  - ipython
-  - numpy
-  - python-gmsh=4.5.6
-  - pytest
-  - pytest-cov
-  - meshio
-  - fenics=2019.1
-  - importlib_metadata
-  - notebook
-  - plotly
-  - matplotlib
-  - ipywidgets
-  - sphinx
-  - sphinx_rtd_theme
-  - make
-  - nbsphinx
-  - recommonmark
-  - typing_extensions
+  - python=3.9
+  - python-gmsh=4.11.1
+  - fenics=2019.1.0
\ No newline at end of file
diff --git a/conda/env_20240207.yaml b/conda/env_20240207.yaml
index e77ecd18a9090c0d521448c48c76ff220d8117d7..554b135747cd06666e6386da53bd0c805f16f288 100644
--- a/conda/env_20240207.yaml
+++ b/conda/env_20240207.yaml
@@ -3,26 +3,26 @@ channels:
   - conda-forge
   - defaults
 dependencies:
-  - python=3.9
+  - python >=3.9,<4.0
   - ipython
   - numpy
-  - python-gmsh>=4.5
-  - pytest
-  - pytest-cov
+  - pandas
+  - matplotlib
+  - python-gmsh=4.11
   - meshio
   - fenics=2019.1
   - importlib_metadata
   - notebook
+  - vtk
   - plotly
-  - matplotlib
+  - pyvista
+  - trame
   - ipywidgets
+  - pytest
+  - pytest-cov
   - sphinx
   - sphinx_rtd_theme
   - make
   - nbsphinx
   - recommonmark
-  - typing_extensions
-  - pyvista
-  - trame
-  - vtk=9.2.6
-  - pandas
\ No newline at end of file
+  - typing_extensions
\ No newline at end of file
diff --git a/conda/env_linux.yaml b/conda/env_linux.yaml
deleted file mode 100644
index bb2199942f013707bbc6980c7e437d226b499065..0000000000000000000000000000000000000000
--- a/conda/env_linux.yaml
+++ /dev/null
@@ -1,24 +0,0 @@
-name: bvpy
-channels:
-  - conda-forge
-  - defaults
-dependencies:
-  - python=3.7
-  - ipython
-  - numpy
-  - python-gmsh=4.5.6
-  - pytest
-  - pytest-cov
-  - meshio
-  - fenics =2019.1.0=py37h89c1867_30
-  - importlib_metadata
-  - notebook
-  - plotly
-  - matplotlib
-  - ipywidgets
-  - sphinx
-  - sphinx_rtd_theme
-  - make
-  - nbsphinx
-  - recommonmark
-  - typing_extensions
diff --git a/conda/env_osx.yaml b/conda/env_osx.yaml
deleted file mode 100644
index 04769b9307d4d2c2f901379cefd5445b98e6d89c..0000000000000000000000000000000000000000
--- a/conda/env_osx.yaml
+++ /dev/null
@@ -1,24 +0,0 @@
-name: bvpy
-channels:
-  - conda-forge
-  - defaults
-dependencies:
-  - python=3.7
-  - ipython
-  - numpy
-  - python-gmsh=4.5.6
-  - pytest
-  - pytest-cov
-  - meshio
-  - fenics =2019.1.0=py37hf985489_30
-  - importlib_metadata
-  - notebook
-  - plotly
-  - matplotlib
-  - ipywidgets
-  - sphinx
-  - sphinx_rtd_theme
-  - make
-  - nbsphinx
-  - recommonmark
-  - typing_extensions
diff --git a/conda/meta.yaml b/conda/meta.yaml
index 84b636a2fc13bd716cc768a1e9dd0784f1bb382f..f9c5622a646f687a35748ea573a3a17d0d721f3b 100644
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@@ -1,47 +1,39 @@
+{% set pyproject = load_file_data('../pyproject.toml', 'toml', from_recipe_dir=True) %}
+{% set project = pyproject.get('project', {}) %}
+{% set deps = project.get('dependencies', {}) %}
+{% set urls = project.get('urls', {}) %}
+# https://docs.conda.io/projects/conda-build/en/latest/resources/define-metadata.html#loading-data-from-other-files
+# Source code for Jinja context: https://github.com/conda/conda-build/blob/main/conda_build/jinja_context.py
+# Then search for the `load_file_data` function.
+
 package:
-  name: bvpy
-  version: {{ environ.get('GIT_DESCRIBE_TAG', 'default') }}
+  name: {{ project.get('name') }}
+  version: {{ project.get('version') }}
 
 source:
-  path: ../
+  path: ../  # to build from local sources
 
 build:
-  preserve_egg_dir: True
-  number: {{ environ.get('GIT_DESCRIBE_NUMBER', 0)}}
+  string: {{ 'py'+ PY_VER.replace('.','') + '_' + environ.get('GIT_DESCRIBE_HASH', 'latest') }}
+  script: python -m pip install .
 
 requirements:
   build:
     - setuptools
-    - python {{ python }}
-  host:
-    - python {{ python }}
+    - python  {{ python }}
   run:
-    - python {{ python }}
-    - ipython
-    - numpy
-    - python-gmsh >=4.5.6
-    - meshio
-    - fenics =2019.1.0=py37h89c1867_30 # [linux]
-    - fenics =2019.1.0=py37hf985489_30 # [osx]
-    - importlib_metadata
-    - notebook
-    - plotly
-    - ipywidgets
-
-# test:
-#   source_files:
-#     - test
-#     - data
-#   requires:
-#     - nose
-#     - coverage
-#   commands:
-#     - nosetests -v test --with-coverage --cover-erase --cover-package=bvpy
-#     - coverage report
+    - python  {{ python }}
+    - python-gmsh=4.11
+    - fenics=2019.1.0
+    {% for dep in deps %}
+    - {{ dep }}
+    {% endfor %}
 
 about:
-  license: Cecill-C
+  home: {{ urls.get('Homepage') }}
+  license: LGPL-3.0-or-later
   license_file: LICENSE
-  summary: Package providing tools to solve BVP (Boundary Value Problem) and IBVP (Initial Boundary Value Problem).
-  dev_url: https://gitlab.inria.fr/mosaic/bvpy
-  doc_source_url: https://mosaic.gitlabpages.inria.fr/bvpy/
+  summary: {{ project.get('description') }}
+  description: |
+    Package providing tools to solve BVP (Boundary Value Problem) and IBVP (Initial Boundary Value Problem).
+  dev_url: {{ urls.get('Repository') }}
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000000000000000000000000000000000000..c5153038bcb278f8c94c8aac61d0c5ae71d1fc64
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,64 @@
+[project]
+name = "bvpy"
+version = "1.1.0"
+description = "A package implementing Boundary Value Problem."
+authors = [
+    {name="Florian Gacon", email="florian.gacon@inria.fr"},
+    {name="Olivier Ali", email="olivier.ali@inria.fr"}
+]
+maintainers= [
+    {name="Manuel Petit", email="manuel.petit@inria.fr"},
+]
+readme = "README.rst"
+requires-python = ">=3.9, <4"
+dependencies = [
+  "numpy<2",
+  "meshio",
+  "importlib_metadata",
+  "plotly",
+  "pyvista[all]",
+  "typing_extensions",
+  "h5py",
+  "nbformat>=4.2.0",
+]
+license = {file = "LICENSE"}
+
+[project.urls]
+Homepage = "https://mosaic.gitlabpages.inria.fr/bvpy/"
+Documentation = "https://mosaic.gitlabpages.inria.fr/bvpy/"
+Repository = "https://gitlab.inria.fr/mosaic/bvpy"
+Issues = "https://gitlab.inria.fr/mosaic/bvpy/issues"
+
+[project.optional-dependencies]
+doc = [
+  "matplotlib",
+  "sphinx",
+  "sphinx_rtd_theme",
+  "nbsphinx",
+  "notebook",
+  "chardet",
+  "recommonmark",
+]
+test = [
+  "pytest",
+  "pytest-cov",
+]
+
+[build-system]
+requires = ["setuptools"]
+build-backend = "setuptools.build_meta"
+
+[tool.coverage.run]
+source = ["src"]
+omit = ["*__init__.py", "test/*", "version.py"]
+
+[tool.coverage.report]
+omit = ["*__init__.py", "test/*"]
+show_missing = true
+
+[tool.pytest.ini_options]
+minversion = "6.0"
+addopts = "-v --cov=bvpy"
+testpaths = [
+  "test",
+]
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
deleted file mode 100644
index e528e57cc1c3f44143b31054b9f02a16389c99f5..0000000000000000000000000000000000000000
--- a/setup.cfg
+++ /dev/null
@@ -1,14 +0,0 @@
-[global]
-# verbose=0
-
-# {# pkglts, wheel
-[bdist_wheel]
-
-# #}
-# {# pkglts, sphinx
-[build_sphinx]
-build-dir=build/sphinx
-# #}
-# {# pkglts, test.nose
-
-# #}
diff --git a/setup.py b/setup.py
deleted file mode 100644
index e7c15882eadea509ead39edb53438e3f8cef1997..0000000000000000000000000000000000000000
--- a/setup.py
+++ /dev/null
@@ -1,42 +0,0 @@
-# {# pkglts, pysetup.kwds
-# format setup arguments
-
-from setuptools import setup, find_packages
-
-
-short_descr = "A package implementing Boundary Value Problem."
-readme = open('README.rst').read()
-history = open('HISTORY.rst').read()
-
-# find packages
-pkgs = find_packages('src')
-
-
-setup_kwds = dict(
-    name='bvpy',
-    version="1.0.1",
-    description=short_descr,
-    long_description=readme + '\n\n' + history,
-    author="Florian Gacon",
-    author_email="florian.gacon@inria.fr",
-    url='',
-    license='cecill-c',
-    zip_safe=False,
-
-    packages=pkgs,
-
-    package_dir={'': 'src'},
-    setup_requires=[],
-    install_requires=[
-        ],
-    tests_require=[],
-    entry_points={},
-    keywords='',
-    )
-# #}
-# change setup_kwds below before the next pkglts tag
-
-# do not change things below
-# {# pkglts, pysetup.call
-setup(**setup_kwds)
-# #}
diff --git a/src/bvpy/__init__.py b/src/bvpy/__init__.py
index d2bfa5f2cf73c5ef9b9c3c6c75c186ee9a1f84c2..d6d4d31a625bc1d69999d1b512c204045b5f1b30 100644
--- a/src/bvpy/__init__.py
+++ b/src/bvpy/__init__.py
@@ -1,15 +1,5 @@
-"""
-belle petite description
-"""
-# {# pkglts, src
-# FYEO
-# #}
-# {# pkglts, version, after src
-from . import version
-
-__version__ = version.__version__
-# #}
-
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
 
 import logging
 from fenics import MPI as _MPI
diff --git a/src/bvpy/boundary_conditions/boundary.py b/src/bvpy/boundary_conditions/boundary.py
index 17666cb4182b3cad5c75483e8fb76ff19dbfbb1f..b35ff5bedadf9f07eb9537ae8b3e898e7ea41755 100644
--- a/src/bvpy/boundary_conditions/boundary.py
+++ b/src/bvpy/boundary_conditions/boundary.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.boundary_conditions.boundary
diff --git a/src/bvpy/boundary_conditions/dirichlet.py b/src/bvpy/boundary_conditions/dirichlet.py
index 7d5d90f4f631d209495a023b0c93c5a14c887516..7255301d6c23f571b5e24183ac6ed39252af2cd7 100644
--- a/src/bvpy/boundary_conditions/dirichlet.py
+++ b/src/bvpy/boundary_conditions/dirichlet.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.boundary_conditions.dirichlet
diff --git a/src/bvpy/boundary_conditions/neumann.py b/src/bvpy/boundary_conditions/neumann.py
index 0d46e5efbc183c31df3b00b231e0e7599f55b8da..a6b4c3c7e9e887e0f05161ee007baeb3addd45a5 100644
--- a/src/bvpy/boundary_conditions/neumann.py
+++ b/src/bvpy/boundary_conditions/neumann.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       femtk.boundary_conditions.neumann
diff --git a/src/bvpy/boundary_conditions/periodic.py b/src/bvpy/boundary_conditions/periodic.py
index 8d59d977be436233b5999ab3946989b634d96bb5..8291265ca9819b4c773065df473ac23357e3b959 100644
--- a/src/bvpy/boundary_conditions/periodic.py
+++ b/src/bvpy/boundary_conditions/periodic.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.boundary_conditions.periodic
diff --git a/src/bvpy/bvp.py b/src/bvpy/bvp.py
index b587a8299f00639e8a59991230cc1886ddab0e4a..4cc6fcf4e37ee7e2552f9dc607a2a51c8135d953 100644
--- a/src/bvpy/bvp.py
+++ b/src/bvpy/bvp.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.bvp
@@ -387,7 +387,15 @@ class BVP(object):
         parameters = self.solver.parameters._param
         krylov_parameter = self.solver.parameters._param['krylov_solver']
         for key in kwargs.keys():
-            if key in parameters:
+            if (key in krylov_parameter) and (key in parameters) and isinstance(self.solver, LinearSolver):
+                krylov_parameter[key] = kwargs[key]
+            elif key == 'krylov_solver':
+                for subkey in kwargs[key].keys():
+                    if subkey in krylov_parameter:
+                        krylov_parameter[subkey] = kwargs[key][subkey]
+                    else:
+                        logger.warn(str(subkey) + ' is not in the list of parameter for the krylov solver')
+            elif key in parameters:
                 parameters[key] = kwargs[key]
             elif key in krylov_parameter:
                 krylov_parameter[key] = kwargs[key]
diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py
index 40f689c9021cd975e2c4d107ac102838997ee43d..21836c3e5b6fecdb63a3b6e2f96a8157360fdf9f 100644
--- a/src/bvpy/domains/abstract.py
+++ b/src/bvpy/domains/abstract.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.domains.abstract
diff --git a/src/bvpy/domains/custom_domain.py b/src/bvpy/domains/custom_domain.py
index 82e2906bb9eb9f859428f594ce2e591d108a401b..d4726f64b7d9e523475b78d156ced987602f34ab 100644
--- a/src/bvpy/domains/custom_domain.py
+++ b/src/bvpy/domains/custom_domain.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.templates.custom_polygonal
diff --git a/src/bvpy/domains/custom_polygonal.py b/src/bvpy/domains/custom_polygonal.py
index e8b18bc2deefa5b0da0a7f1c8f881a4b1c446577..dba0dcecef88a9876c8db94c7bec275c22bbc797 100644
--- a/src/bvpy/domains/custom_polygonal.py
+++ b/src/bvpy/domains/custom_polygonal.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.templates.custom_polygonal
diff --git a/src/bvpy/domains/geometry.py b/src/bvpy/domains/geometry.py
index 98bc8719bfeafd4a14fdd25a57a58adfec9d9d33..2b21b35e3b0b333e9e708546355c0f5fef8462d7 100644
--- a/src/bvpy/domains/geometry.py
+++ b/src/bvpy/domains/geometry.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       femtk.utils.geometry
diff --git a/src/bvpy/domains/primitives.py b/src/bvpy/domains/primitives.py
index 8c759f0ced54089499fdba466f86e24671ded081..717f06525ac8978f1bae3e2d44552db2583b3bc5 100644
--- a/src/bvpy/domains/primitives.py
+++ b/src/bvpy/domains/primitives.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.domains.primitives
diff --git a/src/bvpy/ibvp.py b/src/bvpy/ibvp.py
index bc5d58f584cf99a19700d8b819ffe34c2e0f6fdd..1408e8f85e29f3da25712971c7203b55d74c37ca 100644
--- a/src/bvpy/ibvp.py
+++ b/src/bvpy/ibvp.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.ibvp
@@ -68,10 +68,14 @@ class IBVP(BVP):
     scheme
 
     """
+
     def __init__(self, domain=None, vform=None, bc=None,
                  initial_state=None, scheme=ImplicitEuler):
         super().__init__(domain=domain, vform=vform, bc=bc)
 
+        self._relative_tolerance = None
+        self._absolute_tolerance = None
+        self._linear_solver = None
         self.time = Time()
         self.scheme = scheme(vform)
         self.previous_solution = self.solution.copy(deepcopy=True)
@@ -110,43 +114,38 @@ class IBVP(BVP):
 
         if store_steps:
             self.solution_steps = OrderedDict()
-            self.solution_steps[self.time.current] =\
+            self.solution_steps[self.time.current] = \
                 self.solution.copy(deepcopy=True)
 
         if filename is not None:
-            file = fe.XDMFFile(MPI_COMM_WORLD, filename)
-            file.parameters['functions_share_mesh'] = True
-            file.write(self.solution, self.time.current)
+            with fe.XDMFFile(MPI_COMM_WORLD, filename) as file:
+                file.parameters['functions_share_mesh'] = True  # functions share the same mesh (reduce file size)
+                file.write(self.solution, self.time.current)
             logger.setLevel(logging.WARNING)
             self._progressbar.set_range(t_i, t_f)
 
-            while t_f-self.time.current > dt:
+            while t_f - self.time.current >= dt:
                 self.step(dt, **kwargs)
-                file.write(self.solution, self.time.current)
+                with fe.XDMFFile(MPI_COMM_WORLD, filename) as file:
+                    file.write(self.solution, self.time.current)
                 if store_steps:
-                    self.solution_steps[self.time.current] =\
+                    self.solution_steps[self.time.current] = \
                         self.solution.copy(deepcopy=True)
                 self._progressbar.update(self.time.current)
                 self._progressbar.show()
-
-            self.step(t_f-self.time.current, **kwargs)
-            file.write(self.solution, self.time.current)
-            self._progressbar.update(self.time.current)
-            self._progressbar.show()
-            file.close()
-
         else:
             logger.setLevel(logging.WARNING)
             self._progressbar.set_range(t_i, t_f)
 
-            while t_f-self.time.current > dt:
+            while t_f - self.time.current >= dt:
                 self.step(dt, **kwargs)
                 if store_steps:
-                    self.solution_steps[self.time.current]
+                    self.solution_steps[self.time.current] = \
+                        self.solution.copy(deepcopy=True)
                 self._progressbar.update(self.time.current)
                 self._progressbar.show()
 
-            self.step(t_f-self.time.current, **kwargs)
+            self.step(t_f - self.time.current, **kwargs)
             self._progressbar.update(self.time.current)
             self._progressbar.show()
 
@@ -160,9 +159,10 @@ class IBVP(BVP):
     def solve(self, lhs, rhs, **kwargs):
         if isinstance(rhs, Zero):
             for neumann in self.neumannCondition['constant']:
-                lhs -= fe.dot(neumann.apply(self.domain.bdata),self.testFunction)*self.domain.ds(neumann._id)
+                lhs -= fe.dot(neumann.apply(self.domain.bdata), self.testFunction) * self.domain.ds(neumann._id)
             for neumann in self.neumannCondition['variable']:
-                lhs -= fe.dot(neumann.apply(self.functionSpace, self.domain.bdata),self.testFunction)*self.domain.ds(neumann._id)
+                lhs -= fe.dot(neumann.apply(self.functionSpace, self.domain.bdata), self.testFunction) * self.domain.ds(
+                    neumann._id)
 
             self.solver = NonLinearSolver()
             self._set_solver_parameter(**kwargs)
@@ -173,9 +173,10 @@ class IBVP(BVP):
         else:
 
             for neumann in self.neumannCondition['constant']:
-                rhs += fe.dot(neumann.apply(self.domain.bdata),self.testFunction)*self.domain.ds(neumann._id)
+                rhs += fe.dot(neumann.apply(self.domain.bdata), self.testFunction) * self.domain.ds(neumann._id)
             for neumann in self.neumannCondition['variable']:
-                rhs += fe.dot(neumann.apply(self.functionSpace, self.domain.bdata),self.testFunction)*self.domain.ds(neumann._id)
+                rhs += fe.dot(neumann.apply(self.functionSpace, self.domain.bdata), self.testFunction) * self.domain.ds(
+                    neumann._id)
 
             self.solver = LinearSolver()
             self._set_solver_parameter(**kwargs)
diff --git a/src/bvpy/solvers/linear.py b/src/bvpy/solvers/linear.py
index 7a4b6dee90af601f87df4f1ab4dbf82d029efcc8..e3ac7b0641c95f3351694c682c9c2cd718a68b57 100644
--- a/src/bvpy/solvers/linear.py
+++ b/src/bvpy/solvers/linear.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.solvers.linear
diff --git a/src/bvpy/solvers/nonlinear.py b/src/bvpy/solvers/nonlinear.py
index a4028bd4c4dd22a98627a3eb4051338c5d94294e..b136475409c446211c4080961cb21d9253413389 100644
--- a/src/bvpy/solvers/nonlinear.py
+++ b/src/bvpy/solvers/nonlinear.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.solvers.nonlinear
@@ -18,6 +18,7 @@
 #           https://www.gnu.org/licenses/lgpl-3.0.en.html
 #
 # -----------------------------------------------------------------------
+import logging
 
 import fenics as fe
 import time
@@ -180,16 +181,33 @@ class NonLinearSolver():
                 self._F, self._Jac, self._bnd, self._ps)
             self._solver = fe.PETScSNESSolver(MPI_COMM_WORLD)
 
-            paramters_dict = {}
-            paramters_dict.update(self.parameters['krylov_solver'])
-            paramters_dict['linear_solver'] = self.parameters['linear_solver']
+            for key, val in self.parameters._param.items():
+                if isinstance(val, dict):
+                    for subkey, subval in val.items():
+                        try:
+                            self._solver.parameters[key][subkey] = subval
+                        except:
+                            logging.warning(f'The parameter {key} {subkey} was not set!')
+                else:
+                    try:
+                        self._solver.parameters[key] = val
+                    except:
+                        logging.warning(f'The parameter {key} was not set!')
+
+            #paramters_dict = {}
+            #paramters_dict.update(self.parameters['krylov_solver'])
+            #paramters_dict['linear_solver'] = self.parameters['linear_solver']
             # TODO :: Check with Flo : is this normal that the NL solver has a 'linear_solver' param ?
-            paramters_dict['preconditioner'] = self.parameters['preconditioner']
-            paramters_dict['line_search'] = self.parameters['line_search']
-            paramters_dict['maximum_iterations'] = self.parameters['maximum_iterations']
-            monitor = paramters_dict.pop('monitor_convergence')
-            self._solver.parameters.update(paramters_dict)
-            self._solver.parameters['report'] = monitor
+            #paramters_dict['preconditioner'] = self.parameters['preconditioner']
+            #paramters_dict['line_search'] = self.parameters['line_search']
+            #paramters_dict['maximum_iterations'] = self.parameters['maximum_iterations']
+
+
+            #self._solver.parameters.update(paramters_dict)
+            #self._solver.parameters['report'] = monitor
+            #self._solver.parameters['krylov_solver'].update(self.parameters['krylov_solver']) # set the krylov parameters
+            #self._solver.parameters['absolute_tolerance']
+
             self._snes = self._solver.snes()
             self._snes.setConvergenceHistory()
 
@@ -203,17 +221,16 @@ class NonLinearSolver():
             of the performed computation.
         """
         if self._type == 'petsc':
-            if MPI_COMM_WORLD.Get_rank() == 0: # print only process of rank 0
-                logger.info(' Backend       : petsc')
-                logger.info(' Linear Solver : '
-                            + str(self._solver.parameters['linear_solver']))
-                logger.info(' Preconditioner: '
-                            + str(self._solver.parameters['preconditioner']))
-                logger.info(f" atol: {self._solver.parameters['absolute_tolerance']}")
-                logger.info(f" rtol: {self._solver.parameters['relative_tolerance']}")
-                logger.info(' Size          : '
-                            + str(self._u.function_space().dim()))
-                logger.info(' Solving NonLinearProblem ...')
+            logger.info(' Backend       : petsc')
+            logger.info(' Linear Solver : '
+                        + str(self._solver.parameters['linear_solver']))
+            logger.info(' Preconditioner: '
+                        + str(self._solver.parameters['preconditioner']))
+            logger.info(f" atol: {self._solver.parameters['absolute_tolerance']}")
+            logger.info(f" rtol: {self._solver.parameters['relative_tolerance']}")
+            logger.info(' Size          : '
+                        + str(self._u.function_space().dim()))
+            logger.info(' Solving NonLinearProblem ...')
 
             start = time.time()
             self._solver.solve(self._problem, self._u.vector())
@@ -221,10 +238,9 @@ class NonLinearSolver():
 
             residuals = self._snes.getConvergenceHistory()[0]
 
-            if MPI_COMM_WORLD.Get_rank() == 0:  # print only process of rank 0
-                logger.info(' ... Done [' + str(end - start) + ' s]')
-                conv = self._snes.getConvergedReason()
-                logger.info(' Iterations    : ' + str(len(residuals)))
-                logger.info(' Resiudal      : ' + str(residuals[-1]))
+            logger.info(' ... Done [' + str(end - start) + ' s]')
+            conv = self._snes.getConvergedReason()
+            logger.info(' Iterations    : ' + str(len(residuals)))
+            logger.info(' Resiudal      : ' + str(residuals[-1]))
 
             return residuals
\ No newline at end of file
diff --git a/src/bvpy/solvers/parameter.py b/src/bvpy/solvers/parameter.py
index 3e50ca9043cd0572ef9030d97d2a2c1f3996e7c9..8607b9b9dfcc50711817eed18ec1550a7517797f 100644
--- a/src/bvpy/solvers/parameter.py
+++ b/src/bvpy/solvers/parameter.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.solvers.parameter
@@ -63,7 +63,8 @@ class SolverParameter():
                         'maximum_iterations': 50,
                         'maximum_residual_evaluations': 2000,
                         'method': 'default',
-                        'preconditioner': 'default',
+                        'report': False,
+                        'preconditioner': 'none',
                         'relative_tolerance': 1e-9,
                         'sign': 'default',
                         'solution_tolerance': 0.000000000000}
diff --git a/src/bvpy/solvers/time_integration.py b/src/bvpy/solvers/time_integration.py
index 0bbb8aaa7f6a2aa2c456a487933a592bb598b5fa..864aabcaf07e7d847d38c7bc3a61bb47cd0cf5ba 100644
--- a/src/bvpy/solvers/time_integration.py
+++ b/src/bvpy/solvers/time_integration.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.solvers.time_integration
diff --git a/src/bvpy/utils/io.py b/src/bvpy/utils/io.py
index 791b76575f95d4cf2e1b9f318d7619bbff8cbed5..6a30e947cb37c28521c6996f676521387edc3a75 100644
--- a/src/bvpy/utils/io.py
+++ b/src/bvpy/utils/io.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
diff --git a/src/bvpy/utils/post_processing.py b/src/bvpy/utils/post_processing.py
index d8d24ff4e3f2ce91c59d0d0934dad593024f7613..0c041449f86d6d0cb5e1fb85d9283f5a2c39a8dc 100644
--- a/src/bvpy/utils/post_processing.py
+++ b/src/bvpy/utils/post_processing.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
diff --git a/src/bvpy/utils/pre_processing.py b/src/bvpy/utils/pre_processing.py
index fb574bc107f997c4f46cad8b81d31e7da5829db7..0185b9cdf0599dfbd735821d1e52ce8006cbba64 100644
--- a/src/bvpy/utils/pre_processing.py
+++ b/src/bvpy/utils/pre_processing.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
diff --git a/src/bvpy/utils/progress_bar.py b/src/bvpy/utils/progress_bar.py
index e165d35b52c217cbc1bb5b7c682db508d5610efa..fa9cf770c96a3ff71b53d8d4299c89a791a4b207 100644
--- a/src/bvpy/utils/progress_bar.py
+++ b/src/bvpy/utils/progress_bar.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
diff --git a/src/bvpy/utils/tensor.py b/src/bvpy/utils/tensor.py
index 6beb7bd28669d672cd6d5af9b2b712fcc67a2254..9c75581383c00730966cc75686718c059fffc602 100644
--- a/src/bvpy/utils/tensor.py
+++ b/src/bvpy/utils/tensor.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
diff --git a/src/bvpy/utils/visu.py b/src/bvpy/utils/visu.py
index 7fd7c0f796fdc70ae5c9e776c61a99ffa874a465..fed1357b06a14154e5a207929aefef70f4b3a4d1 100644
--- a/src/bvpy/utils/visu.py
+++ b/src/bvpy/utils/visu.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
diff --git a/src/bvpy/utils/visu_pyvista.py b/src/bvpy/utils/visu_pyvista.py
index 3a21cc752e1f7c066d5a72ea266b4763c94dc309..c5775ba68895b2ce809fcb4e2ae058ac3c3f388a 100644
--- a/src/bvpy/utils/visu_pyvista.py
+++ b/src/bvpy/utils/visu_pyvista.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
@@ -340,7 +340,7 @@ def _visu_feDiriclet(diri, val_range = 'auto', cmap: str = 'inferno', plotter: p
             coords_diri = np.array([dofs[i] for i in values.keys()])
             scalars = np.array(list(values.values()))
 
-            return pv.PointSet(coords_diri), scalars
+            return pv.PointSet(_check_is3D(coords_diri)), scalars
 
         elif (L > 1) and (L < 4):  # vector values
             # - reorganize the values
diff --git a/src/bvpy/version.py b/src/bvpy/version.py
deleted file mode 100644
index a4a549a23a2f42613729527d1252f80f455df8d0..0000000000000000000000000000000000000000
--- a/src/bvpy/version.py
+++ /dev/null
@@ -1,18 +0,0 @@
-"""
-Maintain version for this package.
-Do not edit this file, use 'version' section of config.
-"""
-# {# pkglts, version
-#  -*- coding: utf-8 -*-
-
-MAJOR = 1
-"""(int) Version major component."""
-
-MINOR = 0
-"""(int) Version minor component."""
-
-POST = 0
-"""(int) Version post or bugfix component."""
-
-__version__ = f"{MAJOR:d}.{MINOR:d}.{POST:d}"
-# #}
diff --git a/src/bvpy/vforms/abstract.py b/src/bvpy/vforms/abstract.py
index d2a07f42f8b901f9c15f31f85447b676be6f758f..47318ae59fdc44d792edce2313f38ce5aa9ee92c 100644
--- a/src/bvpy/vforms/abstract.py
+++ b/src/bvpy/vforms/abstract.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.utils.interface
diff --git a/src/bvpy/vforms/elasticity.py b/src/bvpy/vforms/elasticity.py
index 3feae0b9f94a20fcef18dbf8cb55c3d86a7d1463..c1e95312abece7d36f9e9a9254639944b0d72ee7 100644
--- a/src/bvpy/vforms/elasticity.py
+++ b/src/bvpy/vforms/elasticity.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.templates.vforms.linear_elasticity
@@ -310,7 +310,7 @@ class HyperElasticForm(ElasticForm):
         elif self.model == 'neo_hookean':
             d = u.geometric_dimension()
             Id = fe.Identity(d)
-            F = Id + fe.nabla_grad(sol)
+            F = Id + fe.grad(sol)
             C = F.T * F
 
             Ic = fe.tr(C)
diff --git a/src/bvpy/vforms/helmholtz.py b/src/bvpy/vforms/helmholtz.py
index 7a431768b8a0ffa08fe4455788adbb22223bbde9..67f295459e8939eca1af8f80077e08001081057e 100644
--- a/src/bvpy/vforms/helmholtz.py
+++ b/src/bvpy/vforms/helmholtz.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.templates.vforms.helmholtz
diff --git a/src/bvpy/vforms/plasticity.py b/src/bvpy/vforms/plasticity.py
index a2c2ca9e3c049b1c4bac25b621666156bdfbea9a..794bb25ed02b0419395cf3965a4a485b972c1ce3 100644
--- a/src/bvpy/vforms/plasticity.py
+++ b/src/bvpy/vforms/plasticity.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.bvp
diff --git a/src/bvpy/vforms/poisson.py b/src/bvpy/vforms/poisson.py
index e13ff0baa651ae2c7ce878b7c3f9cacbbf3bde42..245faa4f853a41d85ed5a7218e2f0e5ed8662762 100644
--- a/src/bvpy/vforms/poisson.py
+++ b/src/bvpy/vforms/poisson.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.templates.vforms.poisson
diff --git a/src/bvpy/vforms/transport.py b/src/bvpy/vforms/transport.py
index 926099fa0551d8bdecaefdd42fe77681459084e9..7245d20a922071c3c50e050902aed399e71ce383 100644
--- a/src/bvpy/vforms/transport.py
+++ b/src/bvpy/vforms/transport.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       bvpy.bvp
diff --git a/test/test_elasticity.py b/test/test_elasticity.py
index 8f2e8ca019b48090803e9e056f41db1a82de485c..20baecdf92903a14033ea25e87d05c1dff5828be 100644
--- a/test/test_elasticity.py
+++ b/test/test_elasticity.py
@@ -42,7 +42,8 @@ class TestLinearElasticityProblem(unittest.TestCase):
 
         diri = dirichlet([0, 0], boundary='near(x, 0)')
         bvp = BVP(self.domain, vform, diri)
-        bvp.solve()
+        bvp.solve(krylov_solver={'absolute_tolerance':1e-10}) # assert that we can set absolute_tolerance for the krylov_solver
+        bvp.solve(absolute_tolerance=1e-10) # assert that the simplest way is also valid (autoset krylov parameter)
 
         comp_def = bvp.solution(self.L, self.W / 2)[1]
         err_rel = abs(self.theo_def_body + comp_def) / self.theo_def_body
@@ -56,7 +57,7 @@ class TestLinearElasticityProblem(unittest.TestCase):
         print(vform._parameters)
         diri = dirichlet([0, 0], boundary='near(x, 0)')
         bvp = BVP(self.domain, vform, diri)
-        bvp.solve(linear_solver='cg')
+        bvp.solve(linear_solver='superlu') # try another linear solver
 
         comp_def = bvp.solution(self.L, self.W / 2)[1]
         err_rel = abs(self.theo_def_body + comp_def) / self.theo_def_body
@@ -70,7 +71,7 @@ class TestLinearElasticityProblem(unittest.TestCase):
                                  plane_stress=True)
         diri = dirichlet([0, 0], boundary='near(x, 0)')
         bvp = BVP(self.domain, vform, diri)
-        bvp.solve(linear_solver='cg', absolute_tolerance=1e-7)
+        bvp.solve(linear_solver='cg', absolute_tolerance=1e-7) # try another linear solver and update absolute_tol
 
         comp_def = bvp.solution(self.L, self.W / 2)[1]
         err_rel = abs(self.theo_def_body + comp_def) / self.theo_def_body
diff --git a/test/test_transport.py b/test/test_transport.py
index 5f730ca8cc1da3e5d097178ebba8f82bdb0682be..c932ed64c22457cdd171a6f0bb45224d150d853f 100644
--- a/test/test_transport.py
+++ b/test/test_transport.py
@@ -73,7 +73,7 @@ class TestTransport(unittest.TestCase):
 
         k1 = 9
         k2 = 11
-        dt = 0.1
+        dt = 0.05
 
         vform = CoupledTransportForm()
         vform.add_species('A', 1, lambda a, b: k1*(b-(a*b)/(1+b*b)))
@@ -104,4 +104,4 @@ class TestTransport(unittest.TestCase):
         ibvp.integrate(0, 1, dt,
                        linear_solver='mumps',
                        absolute_tolerance=1e-10,
-                       relative_tolerance=1e-10)
+                       relative_tolerance=1e-10, report=True)
diff --git a/test/test_visu_pyvista.py b/test/test_visu_pyvista.py
new file mode 100644
index 0000000000000000000000000000000000000000..de0c2239313605fa936685182334bc9e77b3d6d3
--- /dev/null
+++ b/test/test_visu_pyvista.py
@@ -0,0 +1,93 @@
+import unittest
+import shutil
+# -- test class
+
+from bvpy.utils.visu_pyvista import (_check_is3D, _check_scalar_bar_title, _visu_feDiriclet, _visu_feFunction,
+                                     _visu_feFunctionSizet, _visu_feMesh, visualize)
+
+import pyvista as pv
+import fenics as fe
+import numpy as np
+class TestVisu(unittest.TestCase):
+    def setUp(self):
+        """Initiates the test.
+        """
+        self.mesh = fe.UnitSquareMesh(2, 2)
+        self.V = fe.FunctionSpace(self.mesh, 'P', 1)
+        self.f = fe.project(fe.Constant(0), self.V)
+        self.bc = fe.DirichletBC(self.V, fe.Constant(0), "on_boundary")
+        self.mf = fe.MeshFunction("size_t", self.mesh, 2)
+
+    def tearDown(self):
+        """Concludes and closes the test.
+        """
+        pass
+
+    def test_is3D(self):
+        self.assertEqual(_check_is3D(np.ones((10, 1))).shape[1], 3)  # check if 1D data are converted in 3D data
+        self.assertEqual(_check_is3D(np.ones((10, 2))).shape[1], 3) # check if 2D data are converted in 3D data
+        self.assertEqual(_check_is3D(np.ones((10, 3))).shape[1], 3)  # check if 3D data re unchanged
+        self.assertRaises(AssertionError, _check_is3D, np.ones((10, 4)))  # wrong input if 4D
+
+    def test_check_scalar_bar_title(self):
+        pl = pv.Plotter(off_screen=True)
+        pl = _visu_feFunction(self.f, plotter=pl, show_plot=False)
+        scalar_bar_name = list(pl.scalar_bars.keys())[0]
+        self.assertNotEqual(scalar_bar_name, _check_scalar_bar_title(pl, scalar_bar_name)) # generate unique name ?
+        pl.close()
+
+    def test_visu_mesh(self):
+        # - use dedicated function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(_visu_feMesh(self.mesh, plotter=pl, show_plot=False))
+        pl.close()
+
+        # - use generic function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(visualize(self.mesh, plotter=pl, show_plot=False))
+        pl.close()
+
+    def test_visu_function(self):
+        # - use dedicated function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(_visu_feFunction(self.f, plotter=pl, scale=1, show_plot=False))
+        pl.close()
+
+        # - use generic function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(visualize(self.f, plotter=pl, scale=1, show_plot=False))
+        pl.close()
+    def test_visu_dirichlet(self):
+        # - use dedicated function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(_visu_feDiriclet(self.bc, plotter=pl, show_plot=False))
+        pl.close()
+
+        # - use generic function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(visualize(self.bc, plotter=pl, show_plot=False))
+        pl.close()
+    def test_visu_meshfunc(self):
+        # - use dedicated function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(_visu_feFunctionSizet(self.mf, plotter=pl, show_plot=False, cmap="inferno"))
+        pl.close()
+
+        # - use generic function
+        pl = pv.Plotter(off_screen=True)
+        self.assertIsNotNone(visualize(self.mf, plotter=pl, show_plot=False, cmap="inferno"))
+        pl.close()
+    def test_multiplot(self):
+        # - use generic function
+        pl = pv.Plotter(off_screen=True, shape=(1, 4))
+        pl.subplot(0, 0)
+        visualize(self.mesh, plotter=pl, show_plot=False)
+        pl.subplot(0, 1)
+        visualize(self.bc, plotter=pl, show_plot=False)
+        pl.subplot(0, 2)
+        visualize(self.mf, plotter=pl, show_plot=False, cmap="inferno")
+        pl.subplot(0, 3)
+        visualize(self.f, plotter=pl, scale=1, show_plot=False)
+        self.assertIsNotNone(pl)
+        pl.close()
+