diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index eaec99953ba5d4c8a84ed45811373b238a3472ad..cb1c83736af54c01c899e7c8fa609a9572858d31 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -11,80 +11,66 @@ default:
         - mn2l
 
 variables:
-    GIT_CLONE_PATH: /builds/$CI_PROJECT_PATH/FPM_BUILD/fpm # cannot work in /builds/$CI_PROJECT_PATH
+    GIT_SUBMODULE_STRATEGY: recursive
+    GIT_STRATEGY: clone # workaround full clone for each pipeline (https://gitlab.com/gitlab-org/gitlab-runner/-/issues/26993)
+    GIT_LFS_SKIP_SMUDGE: 1 # do not pull LFS
 
 stages:
     - build
-#    - fmt_dox
+    - test
 
 build:
     <<: *global_tag_def
     stage: build
     script:
-        - printenv | sort
-        - cd ..
-        # get and build waves
-        - rm -rf waves
-        - wget -q https://gitlab.uliege.be/am-dept/waves/-/archive/master/waves-master.tar.bz2
-        - tar xf waves-master.tar.bz2
-        - rm waves-master.tar.bz2
-        - mv waves-master waves
-        - cd waves
+        - git submodule init
+        - git submodule update
+        - rm -rf build workspace
         - mkdir build
         - cd build
-        - cmake -Wno-dev -C ../CMake/disable-trilinos.cmake ..
+        - cmake -Wno-dev ..
         - make -j $(nproc)
-        # check code format
+    artifacts:
+        paths:
+            - build/
+        expire_in: 1 day
+
+format:
+    <<: *global_tag_def
+    stage: build
+    script:
         - clang-format --version # we use clang-format-10 exclusively
-        - cd ../../fpm
-        - ${CI_PROJECT_DIR}/../waves/scripts/format_code.py
+        - ./ext/amfe/scripts/format_code.py
         - mkdir -p patches
         - if git diff --patch --exit-code > patches/clang-format.patch; then echo "Clang format changed nothing"; else echo "Clang format found changes to make!"; false; fi
-        # build fpm
-        - mkdir build
-        - cd build
-        - cmake -Wno-dev .. # -DCMAKE_PREFIX_PATH=${CI_PROJECT_DIR}/../waves (handled by default)
-        - make -j $(nproc)
-        # build documentation
-        - make dox
-        # test
-        - ctest -j $(nproc) --output-on-failure #--verbose
-        - mv ${CI_PROJECT_DIR}/../waves/scripts/format_code.py . # ulgy way to keep a script we need later...
-    #timeout: 10 hours  # will be available in 12.3
     artifacts:
         paths:
-            - build/
             - patches/
         expire_in: 1 day
+        when: on_failure
+    dependencies:
+        - build
+    allow_failure: true
 
-#format:
-#    <<: *global_tag_def
-#    stage: fmt_dox
-#    script:
-#        - clang-format --version # we use clang-format-10 exclusively
-#        - ./build/format_code.py
-#        - mkdir -p patches
-#        - if git diff --patch --exit-code > patches/clang-format.patch; then echo "Clang format changed nothing"; else echo "Clang format found changes to make!"; false; fi
-#    artifacts:
-#        paths:
-#            - patches/
-#        expire_in: 1 day
-#        when: on_failure
-#    dependencies:
-#        - build_test
-#    allow_failure: true
-#
-#            
-#doxygen:
-#    <<: *global_tag_def
-#    stage: fmt_dox
-#    script:
-#        - cd build
-#        - make dox
-#    artifacts:
-#        paths:
-#            - build/doxygen/
-#        expire_in: 1 day
-#    dependencies:
-#        - build_test
+doxygen:
+    <<: *global_tag_def
+    stage: test
+    script:
+        - cd build
+        - make dox
+    artifacts:
+        paths:
+            - build/doxygen/
+        expire_in: 1 week
+    dependencies:
+        - build
 
+ctest:
+    <<: *global_tag_def
+    stage: test
+    script:
+        - cd build
+        - ctest --output-on-failure -j $(nproc)
+    #timeout: 10 hours  # will be available in 12.3
+    dependencies:
+        - build
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000000000000000000000000000000000000..e84a099a335811239805e43e471f8ab027d40310
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "ext/amfe"]
+	path = ext/amfe
+	url = ../amfe.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9df63d97551759d93dd78048f03de57116c0bce5..0258585af917402000de581408f66c6bdd818d4c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -116,9 +116,19 @@ ELSE()
     MESSAGE("Doxygen need to be installed to generate the doxygen documentation")
 ENDIF()
 
+# -- Define (for SWIG to detect definitions)
+INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}) # to find "amfe_def.h"
+
 # -- CTest
 ENABLE_TESTING()
 
 # -- Modules
+ADD_SUBDIRECTORY( ext )
 ADD_SUBDIRECTORY( fpm )
 
+# -- Final
+MESSAGE(STATUS "PROJECT: ${CMAKE_PROJECT_NAME}")
+MESSAGE(STATUS "* SYSTEM NAME=\"${CMAKE_SYSTEM_NAME}\"")
+MESSAGE(STATUS "* CXX COMPILER: ${CMAKE_CXX_COMPILER_ID}")
+MESSAGE(STATUS "* CXX STANDARD: ${CMAKE_CXX_STANDARD}")
+MESSAGE(STATUS "* BUILD TYPE: ${CMAKE_BUILD_TYPE}")
diff --git a/README.md b/README.md
index 2212ff5aedf5dbafa1f92f4efea1cfac81893ea6..1650e532c977f98b427792bdc1dac5b14ed9fa05 100644
--- a/README.md
+++ b/README.md
@@ -10,12 +10,12 @@ Set of python/C++ modules solving the linear potential equation and correcting i
 - [x] [eigen3](http://eigen.tuxfamily.org/) support for linear algebra
 
 ## Build
-fpm is linked against [waves](https://gitlab.uliege.be/am-dept/waves). Once all the required packages and waves are installed:
+fpm depends on [amfe](https://gitlab.uliege.be/am-dept/amfe). Once all the required packages are installed:
 
 ```bash
 git clone git@gitlab.uliege.be:am-dept/fpm.git # get the code
 mkdir build && cd build
-cmake -DCMAKE_PREFIX_PATH=/path/to/waves .. # configure and reference waves package
+cmake .. # configure
 make # compile
 ctest # test
 ```
diff --git a/ext/CMakeLists.txt b/ext/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fddb82f468f6916b39b7b0129565ecad2bb01d05
--- /dev/null
+++ b/ext/CMakeLists.txt
@@ -0,0 +1 @@
+ADD_SUBDIRECTORY( amfe)
\ No newline at end of file
diff --git a/ext/amfe b/ext/amfe
new file mode 160000
index 0000000000000000000000000000000000000000..054ba00c327e4f2745ff9c559832247f90c3cec9
--- /dev/null
+++ b/ext/amfe
@@ -0,0 +1 @@
+Subproject commit 054ba00c327e4f2745ff9c559832247f90c3cec9
diff --git a/fpm/_src/CMakeLists.txt b/fpm/_src/CMakeLists.txt
index a137601db122112709848b228ee3a3d3e3c36695..bd2af2eefd81b35831013f4ad3824436d845b133 100644
--- a/fpm/_src/CMakeLists.txt
+++ b/fpm/_src/CMakeLists.txt
@@ -22,16 +22,12 @@ FILE(GLOB ISRCS *.i)
 SET(CMAKE_SWIG_FLAGS "")
 SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
 
-FOREACH(dir ${WAVES_INCLUDE_DIRS})
-    LIST(APPEND WAVES_SINCLUDE_DIRS "-I${dir}")
-ENDFOREACH()
-FOREACH(dir ${WAVES_SWIG_DIRS})
-    LIST(APPEND WAVES_SINCLUDE_DIRS "-I${dir}")
-ENDFOREACH()
-
-SET(SWINCFLAGS
+SET(SWINCFLAGS 
 -I${PROJECT_SOURCE_DIR}/fpm/src
-${WAVES_SINCLUDE_DIRS}
+-I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/src
+-I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
+-I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/src
+-I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
 )
 SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
 
@@ -43,11 +39,13 @@ endif()
 MACRO_DebugPostfix(_fpmw)
 
 TARGET_INCLUDE_DIRECTORIES(_fpmw PRIVATE ${PROJECT_SOURCE_DIR}/fpm/_src
-                                         ${WAVES_SWIG_DIRS}
-                                         ${PYTHON_INCLUDE_PATH}
+                                          ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
+                                          ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
+                                          ${PYTHON_INCLUDE_PATH}
 )
 
-SWIG_LINK_LIBRARIES(fpmw
-                    fpm ${WAVES_LIBRARIES} ${PYTHON_LIBRARIES}
+SWIG_LINK_LIBRARIES(fpmw 
+                    fpm tbox fwk ${PYTHON_LIBRARIES}
 )
 
+
diff --git a/fpm/src/CMakeLists.txt b/fpm/src/CMakeLists.txt
index 15f283e16e513b3fbea89479fd032b522be76ef2..6664a4f037e3b22139c1d0eb106bfab555d8f453 100644
--- a/fpm/src/CMakeLists.txt
+++ b/fpm/src/CMakeLists.txt
@@ -26,12 +26,6 @@ FIND_PACKAGE(EIGEN 3.3.4 REQUIRED)
 TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${EIGEN_INCLUDE_DIRS})
 TARGET_COMPILE_DEFINITIONS(fpm PUBLIC EIGEN_DONT_PARALLELIZE)
 
-# -- WAVES (if no path is provided, assume that it is located next to fpm)
-IF(NOT DEFINED CMAKE_PREFIX_PATH)
-    SET(CMAKE_PREFIX_PATH "${PROJECT_SOURCE_DIR}/../waves/build")
-ENDIF()
-FIND_PACKAGE(WAVES REQUIRED)
-TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${WAVES_INCLUDE_DIRS})
-TARGET_LINK_LIBRARIES(fpm ${WAVES_LIBRARIES})
+TARGET_LINK_LIBRARIES(fpm tbox)
 
 SOURCE_GROUP(base       REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)")
diff --git a/fpm/src/fBody.cpp b/fpm/src/fBody.cpp
index 2b9bf3e9c68900d383cd0ed1185f746cae531f97..398e3c0172b91968025e28860db168cd3b348593 100644
--- a/fpm/src/fBody.cpp
+++ b/fpm/src/fBody.cpp
@@ -30,6 +30,10 @@ using namespace fpm;
 Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
            std::vector<std::string> const &teIds, double xF) : Group(_msh, id), Cl(0), Cd(0), Cs(0), Cm(0)
 {
+    // Initialize element memory
+    for (auto e : tag->elems)
+        e->initValues(false);
+
     // If wakes are requested, check if they are already available, otherwise create them
     bool hasChanged = false;
     size_t j = 0;
@@ -70,7 +74,7 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
                     err << "fpm::Body: zero or negative length wake requested for " << *tag << "!\n";
                     throw std::runtime_error(err.str());
                 }
-                Node *nN = new Node(msh->nodes.back()->no + 1, Eigen::Vector3d(xF, n->pos(1), n->pos(2)));
+                Node *nN = new Node(msh->nodes.back()->no + 1, Eigen::Vector3d(xF, n->pos(1), n->pos(2)), msh->nodes.back()->row + 1);
                 wkNodes.push_back(nN);
                 msh->nodes.push_back(nN);
             }
@@ -85,7 +89,7 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
             std::map<Node *, Node *> teMap;
             for (auto n : teNodes)
             {
-                Node *nN = new Node(msh->nodes.back()->no + 1, n->pos);
+                Node *nN = new Node(msh->nodes.back()->no + 1, n->pos, msh->nodes.back()->row + 1);
                 teMap[n] = nN;
                 msh->nodes.push_back(nN);
             }
diff --git a/fpm/src/fBuilder.cpp b/fpm/src/fBuilder.cpp
index fe5a573d64007dc76083056a06892a3897b2424e..e349f557b3282752e30a6f1a2f810e3fc860f283 100644
--- a/fpm/src/fBuilder.cpp
+++ b/fpm/src/fBuilder.cpp
@@ -23,7 +23,6 @@
 #include "wTag.h"
 #include "wElement.h"
 #include "wNode.h"
-#include "wMem.h"
 
 #define PI 3.14159
 
@@ -80,7 +79,7 @@ void Builder::run()
             Eigen::MatrixXd trsfCoordSrc(3, nVertices);
 
             // Computations related to source panels
-            glob2loc = Glob2loc(false, ej); // Transformation matrix
+            glob2loc = Glob2loc(false, ej);                   // Transformation matrix
             trsfCoordSrc = TrsfCoordSrc(false, ej, glob2loc); // Local coordinates of the source panel
 
             // body panels
@@ -94,8 +93,8 @@ void Builder::run()
                     Eigen::Vector2d mutau;
 
                     // Computations related to target panels
-                    trsfCoordTgt = TrsfCoordTgt(false, ei, ej, glob2loc); // Local coordinates of the target panel
-                    distance = Distance(false, ei, ej, glob2loc); // Distances between source and target panels
+                    trsfCoordTgt = TrsfCoordTgt(false, ei, ej, glob2loc);               // Local coordinates of the target panel
+                    distance = Distance(false, ei, ej, glob2loc);                       // Distances between source and target panels
                     mutau = MuTau(false, ei, ej, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
 
                     // Assembly of the matrices
@@ -119,7 +118,7 @@ void Builder::run()
                 Eigen::MatrixXd trsfCoordSrc(3, nVertices);
 
                 // Computations related to source panels
-                glob2loc = Glob2loc(true, ej); // Transformation matrix
+                glob2loc = Glob2loc(true, ej);                   // Transformation matrix
                 trsfCoordSrc = TrsfCoordSrc(true, ej, glob2loc); // Local coordinates of the source panel
 
                 // body panels
@@ -133,8 +132,8 @@ void Builder::run()
                         Eigen::Vector2d mutau;
 
                         // Computations related to target panels
-                        trsfCoordTgt = TrsfCoordTgt(true, ei, ej, glob2loc); // Local coordinates of the target panel
-                        distance = Distance(true, ei, ej, glob2loc); // Distances between source and target panels
+                        trsfCoordTgt = TrsfCoordTgt(true, ei, ej, glob2loc);               // Local coordinates of the target panel
+                        distance = Distance(true, ei, ej, glob2loc);                       // Distances between source and target panels
                         mutau = MuTau(true, ei, ej, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
 
                         // Assembly of the matrices
@@ -159,7 +158,7 @@ void Builder::run()
                 Eigen::MatrixXd trsfCoordSrc(3, nVertices);
 
                 // Computations related to source panels
-                glob2loc = Glob2loc(false, ew); // Transformation matrix
+                glob2loc = Glob2loc(false, ew);                   // Transformation matrix
                 trsfCoordSrc = TrsfCoordSrc(false, ew, glob2loc); // Local coordinates of the source panel
 
                 // wake panels
@@ -173,8 +172,8 @@ void Builder::run()
                         Eigen::Vector2d wmutau;
 
                         // Computations related to target panels
-                        trsfCoordTgt = TrsfCoordTgt(false, ei, ew, glob2loc); // Local coordinates of the target panel
-                        distance = Distance(false, ei, ew, glob2loc); // Distances between source and target panels
+                        trsfCoordTgt = TrsfCoordTgt(false, ei, ew, glob2loc);                // Local coordinates of the target panel
+                        distance = Distance(false, ei, ew, glob2loc);                        // Distances between source and target panels
                         wmutau = MuTau(false, ei, ew, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
 
                         // Assembly of the matrices
@@ -201,7 +200,7 @@ void Builder::run()
                     Eigen::MatrixXd trsfCoordSrc(3, nVertices);
 
                     // Computations related to source panels
-                    glob2loc = Glob2loc(true, ew); // Transformation matrix
+                    glob2loc = Glob2loc(true, ew);                   // Transformation matrix
                     trsfCoordSrc = TrsfCoordSrc(true, ew, glob2loc); // Local coordinates of the source panel
 
                     // wake panels
@@ -215,8 +214,8 @@ void Builder::run()
                             Eigen::Vector2d wmutau;
 
                             // Computations related to target panels
-                            trsfCoordTgt = TrsfCoordTgt(true, ei, ew, glob2loc); // Local coordinates of the target panel
-                            distance = Distance(true, ei, ew, glob2loc); // Distances between source and target panels
+                            trsfCoordTgt = TrsfCoordTgt(true, ei, ew, glob2loc);                // Local coordinates of the target panel
+                            distance = Distance(true, ei, ew, glob2loc);                        // Distances between source and target panels
                             wmutau = MuTau(true, ei, ew, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
 
                             // Assembly of the matrices
@@ -249,11 +248,11 @@ void Builder::run()
  */
 Eigen::Matrix3d Builder::Glob2loc(bool symy, Element *ej)
 {
-    Eigen::Matrix3d glob2loc; // Transformation matrix
-    Eigen::Vector3d n = ej->normal(); // Unit normal vector of the source panel
+    Eigen::Matrix3d glob2loc;            // Transformation matrix
+    Eigen::Vector3d n = ej->getNormal(); // Unit normal vector of the source panel
 
     // Computation of the unit longitudinal vector - Source panel
-    Eigen::Vector3d cg = ej->getVMem().getCg();
+    Eigen::Vector3d cg = ej->getCg();
     Eigen::Vector3d x0 = ej->nodes[0]->pos; // Coordinates of node 1
     Eigen::Vector3d l = (x0 - cg).normalized();
 
@@ -265,7 +264,7 @@ Eigen::Matrix3d Builder::Glob2loc(bool symy, Element *ej)
     glob2loc.row(1) = p;
     glob2loc.row(2) = n;
 
-    if(symy)
+    if (symy)
     {
         glob2loc.row(1) = -glob2loc.row(1);
         glob2loc.col(1) = -glob2loc.col(1);
@@ -287,7 +286,7 @@ Eigen::MatrixXd Builder::TrsfCoordSrc(bool symy, Element *ej, Eigen::Matrix3d co
         for (int j = 0; j < nVertices; ++j)
             coords(i, j) = ej->nodes[(nVertices - 1) - j]->pos[i]; // Global coordinates of the nodes of the source panel
 
-    if(symy)
+    if (symy)
         coords.row(1) = -coords.row(1);
 
     for (int i = 0; i < nVertices; ++i)
@@ -302,8 +301,8 @@ Eigen::MatrixXd Builder::TrsfCoordSrc(bool symy, Element *ej, Eigen::Matrix3d co
 Eigen::Vector3d Builder::TrsfCoordTgt(bool symy, Element *ei, Element *ej, Eigen::Matrix3d const &glob2loc)
 {
     Eigen::Vector3d trsfCoordTgt;
-    Eigen::Vector3d cg = ei->getVMem().getCg(); // Global coordinates of the CG of the target panel
-    trsfCoordTgt = glob2loc * cg; // Local coordinates of the target panel
+    Eigen::Vector3d cg = ei->getCg(); // Global coordinates of the CG of the target panel
+    trsfCoordTgt = glob2loc * cg;     // Local coordinates of the target panel
 
     return trsfCoordTgt;
 }
@@ -314,13 +313,13 @@ Eigen::Vector3d Builder::TrsfCoordTgt(bool symy, Element *ei, Element *ej, Eigen
 Eigen::Vector3d Builder::Distance(bool symy, Element *ei, Element *ej, Eigen::Matrix3d const &glob2loc)
 {
     Eigen::Vector3d distance; // Distance between the target and source panels
-    Eigen::Vector3d cgSrc = ej->getVMem().getCg();
-    Eigen::Vector3d cgTgt = ei->getVMem().getCg();
+    Eigen::Vector3d cgSrc = ej->getCg();
+    Eigen::Vector3d cgTgt = ei->getCg();
 
     for (int i = 0; i < 3; ++i)
         distance(i) = cgTgt(i) - cgSrc(i);
 
-    if(symy)
+    if (symy)
         distance(1) = cgTgt(1) + 2 * cgSrc(1);
 
     distance = glob2loc * distance;
@@ -347,24 +346,24 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri
         // Coefficients
         for (int i = 1; i <= nVertices; ++i)
         {
-            e(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + distance(2) * distance(2);
-            h(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1));
-            r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) + distance(2) * distance(2));
+            e(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) + distance(2) * distance(2);
+            h(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1));
+            r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) + distance(2) * distance(2));
         }
         for (int i = 1; i <= nVertices; ++i)
         {
             int j = i % nVertices + 1;
 
-            F(i - 1) = (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * e(i - 1) - (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * h(i - 1);
-            G(i - 1) = (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * e(j - 1) - (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * h(j - 1);
+            F(i - 1) = (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) * e(i - 1) - (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * h(i - 1);
+            G(i - 1) = (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) * e(j - 1) - (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * h(j - 1);
         }
         // Doublet term - Influence of a panel on another panel
         for (int i = 1; i <= nVertices; ++i)
         {
             int j = i % nVertices + 1;
 
-            mu = mu + atan2(distance(2) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (F(i - 1) * r(j - 1) - G(i - 1) * r(i - 1)),
-                            distance(2) * distance(2) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * r(i - 1) * r(j - 1) + F(i - 1) * G(i - 1));
+            mu = mu + atan2(distance(2) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * (F(i - 1) * r(j - 1) - G(i - 1) * r(i - 1)),
+                            distance(2) * distance(2) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * r(i - 1) * r(j - 1) + F(i - 1) * G(i - 1));
         }
         mu = -1 / (4 * PI) * mu;
     }
@@ -375,8 +374,8 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri
     {
         int j = i % nVertices + 1;
 
-        d(i - 1) = sqrt((trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) + (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)));
-        r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) + distance(2) * distance(2));
+        d(i - 1) = sqrt((trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1)) + (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) * (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)));
+        r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) + distance(2) * distance(2));
     }
 
     if (ei->no == ej->no && !symy)
@@ -386,7 +385,7 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri
         {
             int j = i % nVertices + 1;
 
-            tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1)));
+            tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1)));
         }
         tau = -1 / (4 * PI) * tau;
     }
@@ -397,7 +396,7 @@ Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::Matri
         {
             int j = i % nVertices + 1;
 
-            tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1)));
+            tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0, i - 1)) * (trsfCoordSrc(1, j - 1) - trsfCoordSrc(1, i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1, i - 1)) * (trsfCoordSrc(0, j - 1) - trsfCoordSrc(0, i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1)));
         }
 
         tau = -1 / (4 * PI) * tau - distance(2) * mu;
diff --git a/fpm/src/fBuilder.h b/fpm/src/fBuilder.h
index 0e721a828abe1a439bfff6eddad9c61e62de1a86..10069df62cfa39d00a31dbe792b3d788ec7f6540 100644
--- a/fpm/src/fBuilder.h
+++ b/fpm/src/fBuilder.h
@@ -36,8 +36,7 @@ namespace fpm
  */
 class FPM_API Builder : public fwk::wSharedObject
 {
-
-std::map<tbox::Element *, int> rows; ///< map linking an element to a matrix cell
+    std::map<tbox::Element *, int> rows; ///< map linking an element to a matrix cell
 
 public:
     std::shared_ptr<Problem> pbl; ///< problem definition
diff --git a/fpm/src/fField.cpp b/fpm/src/fField.cpp
index 37f76c14f7c00553bb30e9df1bfbb06482f5ce37..b0f06a06ecc2262b4463a138b332d505c056f033 100644
--- a/fpm/src/fField.cpp
+++ b/fpm/src/fField.cpp
@@ -15,11 +15,19 @@
  */
 
 #include "fField.h"
+#include "wElement.h"
 #include "wTag.h"
 
 using namespace tbox;
 using namespace fpm;
 
+Field::Field(std::shared_ptr<tbox::MshData> _msh, std::string const &name) : Group(_msh, name)
+{
+    // Initialize element memory
+    for (auto e : tag->elems)
+        e->initValues(true);
+}
+
 void Field::write(std::ostream &out) const
 {
     out << "fpm::Field " << *tag << std::endl;
diff --git a/fpm/src/fField.h b/fpm/src/fField.h
index 6ef5bb8d99a16907ecf09a473d848a04d23348d1..360c422e0725f4a5e2af133b3b7b6fef4a2755e5 100644
--- a/fpm/src/fField.h
+++ b/fpm/src/fField.h
@@ -30,7 +30,7 @@ namespace fpm
 class FPM_API Field : public tbox::Group
 {
 public:
-    Field(std::shared_ptr<tbox::MshData> _msh, std::string const &name) : Group(_msh, name) {}
+    Field(std::shared_ptr<tbox::MshData> _msh, std::string const &name);
     ~Field() { std::cout << "~Field()\n"; }
 
 #ifndef SWIG
diff --git a/fpm/src/fProblem.cpp b/fpm/src/fProblem.cpp
index 5d1b84063a9542cfc192faf9fa33d75fb2e9084c..f2409059bc6916a0a70c775fdd4f6984bb3865d3 100644
--- a/fpm/src/fProblem.cpp
+++ b/fpm/src/fProblem.cpp
@@ -22,7 +22,6 @@
 #include "wTag.h"
 #include "wElement.h"
 #include "wNode.h"
-#include "wMem.h"
 #include "wCache.h"
 using namespace tbox;
 using namespace fpm;
@@ -57,69 +56,69 @@ void Problem::add(std::shared_ptr<Body> b)
 Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &mu, double const &tau)
 {
     Eigen::Vector3d U(0, 0, 0);
-    size_t nVertices = e.nodes.size();                          // Number of nodes per panel
-    size_t nGP = e.getVCache().getVGauss().getN();              // Number of Gauss points on the element
-    Eigen::MatrixXd coords(3,nVertices);
+    size_t nVertices = e.nodes.size();             // Number of nodes per panel
+    size_t nGP = e.getVCache().getVGauss().getN(); // Number of Gauss points on the element
+    Eigen::MatrixXd coords(3, nVertices);
     for (int i = 0; i < 3; ++i)
         for (int j = 0; j < nVertices; ++j)
-            coords(i,j) = e.nodes[(nVertices - 1) - j]->pos[i]; // Coordinates of the panel
-        
+            coords(i, j) = e.nodes[(nVertices - 1) - j]->pos[i]; // Coordinates of the panel
+
     // Import of the shape functions at the Gauss points
-    Eigen::MatrixXd ff(nVertices,nGP);
-    for(int i = 0; i < nVertices; ++i)
-        for(int j = 0; j < nGP; ++j)
-            ff(i,j) = e.getVCache().getSf(j)[i];
+    Eigen::MatrixXd ff(nVertices, nGP);
+    for (int i = 0; i < nVertices; ++i)
+        for (int j = 0; j < nGP; ++j)
+            ff(i, j) = e.getVCache().getSf(j)[i];
 
     // Import of the derivatives of the shape functions at the Gauss points
-    Eigen::MatrixXd dff(2 * nGP,nVertices);
-    for(int i = 0; i < nGP; ++i)
-        for(int j = 0; j < 2; ++j)
+    Eigen::MatrixXd dff(2 * nGP, nVertices);
+    for (int i = 0; i < nGP; ++i)
+        for (int j = 0; j < 2; ++j)
             dff.row(j + (2 * i)) = e.getVCache().getDsf(i).row(j);
-    
+
     // Computation of the Jacobian matrix - xi-derivative (first row of the Jacobian)
-    Eigen::MatrixXd dXi(3,nGP);
-    for(int i = 0; i < nGP; ++i)
+    Eigen::MatrixXd dXi(3, nGP);
+    for (int i = 0; i < nGP; ++i)
     {
         dXi.col(i) << 0, 0, 0;
 
-        for(int j = 0; j < nVertices; ++j)
+        for (int j = 0; j < nVertices; ++j)
         {
-            dXi(0,i) += dff((2 * i),j) * coords(0,j);
-            dXi(1,i) += dff((2 * i),j) * coords(1,j);
-            dXi(2,i) += dff((2 * i),j) * coords(2,j);
+            dXi(0, i) += dff((2 * i), j) * coords(0, j);
+            dXi(1, i) += dff((2 * i), j) * coords(1, j);
+            dXi(2, i) += dff((2 * i), j) * coords(2, j);
         }
     }
 
     // Computation of the Jacobian matrix - eta-derivative (second row of the Jacobian)
-    Eigen::MatrixXd dEta(3,nGP);
-    for(int i = 0; i < nGP; ++i)
+    Eigen::MatrixXd dEta(3, nGP);
+    for (int i = 0; i < nGP; ++i)
     {
         dEta.col(i) << 0, 0, 0;
 
-        for(int j = 0; j < nVertices; ++j)
+        for (int j = 0; j < nVertices; ++j)
         {
-            dEta(0,i) += dff(1 + (2 * i),j) * coords(0,j);
-            dEta(1,i) += dff(1 + (2 * i),j) * coords(1,j);
-            dEta(2,i) += dff(1 + (2 * i),j) * coords(2,j);
+            dEta(0, i) += dff(1 + (2 * i), j) * coords(0, j);
+            dEta(1, i) += dff(1 + (2 * i), j) * coords(1, j);
+            dEta(2, i) += dff(1 + (2 * i), j) * coords(2, j);
         }
     }
 
     // Computation of the Jacobian matrix - cross-product (third row of the Jacobian)
-    Eigen::MatrixXd dZeta(3,nGP);
+    Eigen::MatrixXd dZeta(3, nGP);
 
-    for(int i = 0; i < nGP; ++i)
+    for (int i = 0; i < nGP; ++i)
     {
         Eigen::Vector3d dxi = dXi.col(i);
         Eigen::Vector3d deta = dEta.col(i);
         Eigen::Vector3d dzeta;
         dzeta = (dxi.cross(deta)).normalized();
-        dZeta.col(i) = dzeta; 
+        dZeta.col(i) = dzeta;
     }
 
     // Assembly of the Jacobian and partial computation of the velocity
     // Done for each Gauss point -> For n Gauss points, there are n values of the velocity -> Average
-    Eigen::MatrixXd v(3,nGP);
-    for(int i = 0; i < nGP; ++i)
+    Eigen::MatrixXd v(3, nGP);
+    for (int i = 0; i < nGP; ++i)
     {
         // Creation of the Jacobian
         Eigen::Matrix3d J0;
@@ -130,9 +129,9 @@ Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &mu, doub
         // Inversion of the Jacobian
         Eigen::Matrix3d Jinv = J0.inverse();
 
-        // Removal of the third column 
-        Eigen::MatrixXd J(3,2);
-        for(int j = 0; j < 2; ++j)
+        // Removal of the third column
+        Eigen::MatrixXd J(3, 2);
+        for (int j = 0; j < 2; ++j)
         {
             J.col(j) = Jinv.col(j);
         }
@@ -145,28 +144,28 @@ Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &mu, doub
         }
 
         // Computation of the velocity
-        Eigen::MatrixXd dff0(2,4);
+        Eigen::MatrixXd dff0(2, 4);
         dff0.row(0) = dff.row(2 * i);
         dff0.row(1) = dff.row(1 + (2 * i));
         v.col(i) = J * dff0 * mu_element;
     }
 
     // Average of the velocity
-    for(int i = 0; i < nGP; ++i)
+    for (int i = 0; i < nGP; ++i)
     {
         U += v.col(i);
     }
-    U = U/nGP;
+    U = U / nGP;
 
     // Import of the unit normal vector - Target panel
-    Eigen::Vector3d n = e.normal(); 
+    Eigen::Vector3d n = e.getNormal();
 
     // Contribution of the normal velocity
-    U += (-tau)*n;
+    U += (-tau) * n;
 
     // Contribution of the freestream velocity
     U += Ui();
-    
+
     return U;
 }
 
@@ -221,17 +220,17 @@ void Problem::updateMem()
     // Update volume Jacobian
     if (field)
         for (auto e : field->tag->elems)
-            e->getVMem().update(true);
+            e->update();
     // Update surface Jacobian
     for (auto s : bodies)
     {
         for (auto e : s->tag->elems)
         {
-            e->getVMem().update(false);
+            e->update();
         }
-        for(auto w : s-> wakes)
+        for (auto w : s->wakes)
             for (auto e : w->tag->elems)
-                e->getVMem().update(false);
+                e->update();
     }
 }
 
diff --git a/fpm/src/fSolver.cpp b/fpm/src/fSolver.cpp
index 6e49acc0a6bb2b2d3416971dbc55529f2f07e8b4..36328b2a77c3d8ea40833d79eb116daa014aba02 100644
--- a/fpm/src/fSolver.cpp
+++ b/fpm/src/fSolver.cpp
@@ -22,7 +22,6 @@
 #include "wMshData.h"
 #include "wNode.h"
 #include "wElement.h"
-#include "wMem.h"
 #include "wTag.h"
 #include "wResults.h"
 #include "wMshExport.h"
@@ -36,20 +35,14 @@ using namespace fpm;
 
 /**
  * @brief Initialize the solver and perform sanity checks
- * @authors Adrien Crovato & Axel Dechamps
  */
 Solver::Solver(std::shared_ptr<Builder> _aic) : aic(_aic), Cl(0), Cd(0), Cs(0), Cm(0)
 {
     // Say hi
     std::cout << std::endl;
     std::cout << "*******************************************************************************" << std::endl;
-    std::cout << "**                                    \\_/                                    **" << std::endl;
-    std::cout << "**                           \\_______O(_)O_______/                           **" << std::endl;
-    std::cout << "**                                                                           **" << std::endl;
-    std::cout << "*******************************************************************************" << std::endl;
-    std::cout << "** Hi! My name is FPM v0.1-2007                                              **" << std::endl;
-    std::cout << "** Adrien Crovato & Romain Boman                                             **" << std::endl;
-    std::cout << "** ULiege 2020-2021                                                          **" << std::endl;
+    std::cout << "** Hi! My name is FPM v1.0-2201                                              **" << std::endl;
+    std::cout << "** ULiege 2020-2022                                                          **" << std::endl;
     std::cout << "*******************************************************************************" << std::endl
               << std::endl;
 
@@ -89,8 +82,8 @@ void Solver::run()
     {
         for (int i = 0; i < body->tag->elems.size(); ++i, ++ig)
         {
-            tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->normal());
-            taue[body->tag->elems[i]->no-1] = aic->pbl->Ui().dot(body->tag->elems[i]->normal());
+            tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->getNormal());
+            taue[body->tag->elems[i]->no - 1] = aic->pbl->Ui().dot(body->tag->elems[i]->getNormal());
         }
     }
 
@@ -148,7 +141,7 @@ void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau,
     Eigen::VectorXd phie = -aic->A * mu - aic->B * tau - aic->C * sigma;
 
     mue.resize(mu.size());
-    for(size_t i = 0 ; i < mue.size() ; ++i)
+    for (size_t i = 0; i < mue.size(); ++i)
         mue[i] = mu(i);
 
     auto pbl = aic->pbl;
@@ -173,10 +166,10 @@ void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau,
         // average at nodes
         for (auto n : body->nodes)
         {
-            for (auto e : body->neMap.at(n)) 
+            for (auto e : body->neMap.at(n))
             {
-                vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); 
-                muNodes[n->no - 1] += mapMu.at(e) / body->neMap.at(n).size(); 
+                vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size();
+                muNodes[n->no - 1] += mapMu.at(e) / body->neMap.at(n).size();
             }
         }
     }
@@ -190,10 +183,10 @@ void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau,
         int ig = 0;
         for (auto e : body->tag->elems)
         {
-            U[e->no - 1] = pbl->U(*e, muNodes, tau(ig));    // velocity
-            rho[e->no - 1] = pbl->Rho(U[e->no - 1]);        // density
-            M[e->no - 1] = pbl->Mach(U[e->no - 1]);         // Mach number
-            Cp[e->no - 1] = pbl->Cp(U[e->no - 1]);          // pressure coefficient
+            U[e->no - 1] = pbl->U(*e, muNodes, tau(ig)); // velocity
+            rho[e->no - 1] = pbl->Rho(U[e->no - 1]);     // density
+            M[e->no - 1] = pbl->Mach(U[e->no - 1]);      // Mach number
+            Cp[e->no - 1] = pbl->Cp(U[e->no - 1]);       // pressure coefficient
             ig += 1;
         }
     }
@@ -219,7 +212,7 @@ void Solver::computeLoad()
             body->cLoadZ[i] = 0;
             for (auto e : body->neMap.at(body->nodes[i]))
             {
-                Eigen::Vector3d cLoad = Cp[e->no - 1] * e->getVMem().getVol() * e->normal() / e->nodes.size();
+                Eigen::Vector3d cLoad = Cp[e->no - 1] * e->getVol() * e->getNormal() / e->nodes.size();
                 body->cLoadX[i] += cLoad(0);
                 body->cLoadY[i] += cLoad(1);
                 body->cLoadZ[i] += cLoad(2);
@@ -229,11 +222,11 @@ void Solver::computeLoad()
         // Compute aerodynamic load coefficients (i.e. Load / pressure / surface)
         // compute coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
         body->Cm = 0;
-        Eigen::Vector3d Cf(0, 0, 0);        
+        Eigen::Vector3d Cf(0, 0, 0);
         for (auto e : body->tag->elems)
         {
-            Eigen::Vector3d l = e->getVMem().getCg() - pbl->xR;                       // lever arm
-            Eigen::Vector3d cf = Cp[e->no - 1] * e->getVMem().getVol() * e->normal(); // force coefficient
+            Eigen::Vector3d l = e->getCg() - pbl->xR;                          // lever arm
+            Eigen::Vector3d cf = Cp[e->no - 1] * e->getVol() * e->getNormal(); // force coefficient
             Cf -= cf;
             body->Cm -= (l(2) * cf(0) - l(0) * cf(2)) / (pbl->sR * pbl->cR); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
         }
diff --git a/fpm/src/fWake.cpp b/fpm/src/fWake.cpp
index 0b407616b333aaf83bdb651d62bcbc28cccfa0e8..f8ba0428e60f7afb6da475b3855f285a23cd0718 100644
--- a/fpm/src/fWake.cpp
+++ b/fpm/src/fWake.cpp
@@ -25,6 +25,10 @@ using namespace fpm;
 
 Wake::Wake(std::shared_ptr<MshData> _msh, std::string const &name, Tag *const &wTag) : Group(_msh, name)
 {
+    // Initialize element memory
+    for (auto e : tag->elems)
+        e->initValues(false);
+
     // Create set of TE edges
     std::unordered_set<Edge *, EquEdge, EquEdge> teEdges;
     for (auto e : tag->elems)
@@ -50,7 +54,7 @@ Wake::Wake(std::shared_ptr<MshData> _msh, std::string const &name, Tag *const &w
             if (it != teEdges.end())
             {
                 // check normal direction
-                if (e->normal().dot((*it)->el0->normal()) > 0)
+                if (e->getNormal().dot((*it)->el0->getNormal()) > 0)
                     (*it)->el1 = e;
                 else
                     (*it)->el2 = e;
diff --git a/fpm/tests/n12t.py b/fpm/tests/n12t.py
index fa1d8432039902951a9bc220bd499c472223a4e6..64a73a8b193c7c8f0079499c51e6c0035f089524 100644
--- a/fpm/tests/n12t.py
+++ b/fpm/tests/n12t.py
@@ -27,7 +27,7 @@ def main():
     p['sym'] = 0
     p['spn'] = 5 # true span if sym=0, half span if sym=1
     p['tspn'] = 3 # true span if sym=0, half span if sym=1
-    p['pars'] = {'symY' : p['sym'], 'wSpn' : p['spn'], 'wNc' : 100, 'wNs' : 10, 'tSpn' : p['tspn'], 'tNc' : 100, 'tNs' : 10}
+    p['pars'] = {'symY' : p['sym'], 'wSpn' : p['spn'], 'wNc' : 100, 'wNs' : 10, 'tSpn' : p['tspn'], 'tNc' : 100, 'tNs' : 5}
     # flow parameters
     p['aoa'] = 3 * math.pi / 180 
     p['aos'] = 0
@@ -39,10 +39,10 @@ def main():
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
     if(p['aoa'] == 3 * math.pi / 180 and p['mach'] == 0):
-        tests.add(CTest('CL', solver.Cl, 0.273, 5e-2))
-        tests.add(CTest('CD', solver.Cd, 0.004, 5e-2))
+        tests.add(CTest('CL', solver.Cl, 0.241, 5e-2))
+        tests.add(CTest('CD', solver.Cd, 0.0045, 5e-2))
         tests.add(CTest('CS', solver.Cs, -0.000, 5e-2))
-        tests.add(CTest('CM', solver.Cm, -0.277, 5e-2))
+        tests.add(CTest('CM', solver.Cm, -0.113, 5e-2))
     tests.run()
     
     # eof
diff --git a/run.py b/run.py
new file mode 100755
index 0000000000000000000000000000000000000000..86fb30f36862285d3009c12d5061d010654da023
--- /dev/null
+++ b/run.py
@@ -0,0 +1,36 @@
+#!/usr/bin/env python3
+# -*- coding: utf8 -*-
+# test encoding: à-é-è-ô-ï-€
+
+# Copyright 2020 University of Liège
+# 
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+# 
+#     http://www.apache.org/licenses/LICENSE-2.0
+# 
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# Calls fwk.wutils.run to execute a script as if dartflo was installed
+
+def main():
+    import os.path, sys
+    # adds fwk/tbox to the python path
+    thisdir = os.path.split(os.path.abspath(__file__))[0]
+    fwkdir = os.path.abspath(os.path.join(thisdir, 'ext', 'amfe'))
+    if not os.path.isdir(fwkdir):
+        raise Exception('fpm/ext/amfe not found!\n')
+    sys.path.append(fwkdir)
+    # adds "." to the pythonpath
+    sys.path.append(thisdir)
+
+    import fwk.wutils as wu
+    wu.run()
+
+if __name__ == "__main__":
+    main()
diff --git a/run_fpm.py b/run_fpm.py
deleted file mode 100755
index a63a0ecceddc9604c6a5e03e53a03a1ec11a1082..0000000000000000000000000000000000000000
--- a/run_fpm.py
+++ /dev/null
@@ -1,26 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf8 -*-
-# test encoding: à-é-è-ô-ï-€
-#
-# runs a test as if it was installed
-#   - fixes the python path in a dev environment
-#   - creates a workspace folder
-
-if __name__ == "__main__":
-    import sys
-    import os
-
-    thisdir = os.path.split(os.path.abspath(__file__))[0]
-    
-    # look for fwk
-    fwkdir = os.path.abspath(os.path.join(thisdir, '..', 'waves'))
-    if not os.path.isdir(fwkdir):
-        raise Exception(
-            "'waves' not found - clone it next to fpm (from https://gitlab.uliege.be/am-dept/waves)")
-    sys.path.append(fwkdir)
-
-    # adds "." to the pythonpath (after waves so that waves is found first)
-    sys.path.append(thisdir)
-
-    import run  # maybe we should put the code of run in fwk....
-    run.main(thisdir)