diff --git a/C++/basicClasses/libBasicClasses.cc b/C++/basicClasses/libBasicClasses.cc
index c7f578126f38d439169b6ec6cb3de13b6000c1a0..694eb32fbafae987ed820194bd307451a8d5b809 100644
--- a/C++/basicClasses/libBasicClasses.cc
+++ b/C++/basicClasses/libBasicClasses.cc
@@ -8,6 +8,7 @@
 #include "../basicTypes/datas.h"
 #include "ngData_class.h"
 #include "ngData_class_extended.h"
+#include "../preisach/ngPreisachScalar.h"
 
 
 using basicDataTypes::triple;
@@ -219,37 +220,37 @@ void declare_vec_t(py::module &m) {
 }
 
 
-
-template<typename T, unsigned int DIM>
-void declare_ngData_class(py::module &m, std::string name, std::string suffix=""){
+// declaration in libMain.h
+template<typename T, unsigned int DIM, typename CL_PARENT>
+void declare_ngData_class(py::module &m, std::string name, std::string suffix){
 
     //Data Class
     std::string pyclass_name = std::string("ng") + name +std::string("_class") + std::to_string(DIM) + suffix;
-    py::class_<ngData_class<T, DIM>, shared_ptr<ngData_class<T, DIM>> , CoefficientFunction>
+    py::class_<ngData_class<T, DIM, CL_PARENT>, shared_ptr<ngData_class<T, DIM, CL_PARENT>> , CoefficientFunction>
     (m, pyclass_name.c_str())
     .def(py::init(
         [](shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, const shared_ptr<CoefficientFunction>& ptrMask)
     {
-        return new ngData_class<T, DIM>(mesh, intrule, ptrMask);
+        return new ngData_class<T, DIM, CL_PARENT>(mesh, intrule, ptrMask);
     }), py::arg("mesh"), py::arg("intrule"), py::arg("mask")=nullptr) 
-    .def("Save", &ngData_class<T, DIM>::Save, py::arg("filename"), py::arg("binary")=true, "save all values in the container to a file.")
-    .def("Load", &ngData_class<T, DIM>::Load, py::arg("filename"), py::arg("skipTest")=1, py::arg("binary")=true, "load the values for the container from a file. The arguments of the constructor have to fit. skipTest will make loading faster.")
-    .def("Set", overload_cast_<const shared_ptr<CoefficientFunction>&>()(&ngData_class<T, DIM>::Set), py::arg("CF"), "set the values according to a Coefficient function")
-    .def_property_readonly("numBytes", &ngData_class<T, DIM>::GetNumBytes, "approximation of the current memory demand")
-    .def("__str__", &ngData_class<T, DIM>::toString)
-    .def("__repr__", &ngData_class<T, DIM>::toString)
-    .def("__sum__", &ngData_class<T, DIM>::Sum)
-    .def("Sum", &ngData_class<T, DIM>::Sum)
-    .def("SerialElement", &ngData_class<T, DIM>::serialElement)
-    .def("GetData", &ngData_class<T, DIM>::getData)
-    .def("GetMips", &ngData_class<T, DIM>::getMips)
-    .def("GetMipsAsMatrix", &ngData_class<T, DIM>::getMipsAsMatrix)
-    .def("GetRotatedValues", &ngData_class<T, DIM>::GetRotatedValues, py::arg("angle"), py::arg("axis")=2, "rotate all values around the axis in the center")
-    .def_property("autoUpdate", &ngData_class<T, DIM>::GetAutoUpdate, &ngData_class<T, DIM>::SetAutoUpdate, "Set the additional field to mark the object to be updated")
-    .def_property("derivative", &ngData_class<T, DIM>::GetDerivative, &ngData_class<T, DIM>::SetDerivative, "Set the derivative used for AutoDiff")
-    .def("copy", &ngData_class<T, DIM>::copy)
-    .def_property_readonly("numEntries", &ngData_class<T, DIM>::GetNumEntries, "number of allocated objects of type T")
-    .def_property("vec", &ngData_class<T, DIM>::GetVec, &ngData_class<T, DIM>::SetVec , "vec")
+    .def("Save", &ngData_class<T, DIM, CL_PARENT>::Save, py::arg("filename"), py::arg("binary")=true, "save all values in the container to a file.")
+    .def("Load", &ngData_class<T, DIM, CL_PARENT>::Load, py::arg("filename"), py::arg("skipTest")=1, py::arg("binary")=true, "load the values for the container from a file. The arguments of the constructor have to fit. skipTest will make loading faster.")
+    .def("Set", overload_cast_<const shared_ptr<CoefficientFunction>&>()(&ngData_class<T, DIM, CL_PARENT>::Set), py::arg("CF"), "set the values according to a Coefficient function")
+    .def_property_readonly("numBytes", &ngData_class<T, DIM, CL_PARENT>::GetNumBytes, "approximation of the current memory demand")
+    .def("__str__", &ngData_class<T, DIM, CL_PARENT>::toString)
+    .def("__repr__", &ngData_class<T, DIM, CL_PARENT>::toString)
+    .def("__sum__", &ngData_class<T, DIM, CL_PARENT>::Sum)
+    .def("Sum", &ngData_class<T, DIM, CL_PARENT>::Sum)
+    .def("SerialElement", &ngData_class<T, DIM, CL_PARENT>::serialElement)
+    .def("GetData", &ngData_class<T, DIM, CL_PARENT>::getData)
+    .def("GetMips", &ngData_class<T, DIM, CL_PARENT>::getMips)
+    .def("GetMipsAsMatrix", &ngData_class<T, DIM, CL_PARENT>::getMipsAsMatrix)
+    .def("GetRotatedValues", &ngData_class<T, DIM, CL_PARENT>::GetRotatedValues, py::arg("angle"), py::arg("axis")=2, "rotate all values around the axis in the center")
+    .def_property("autoUpdate", &ngData_class<T, DIM, CL_PARENT>::GetAutoUpdate, &ngData_class<T, DIM, CL_PARENT>::SetAutoUpdate, "Set the additional field to mark the object to be updated")
+    .def_property("derivative", &ngData_class<T, DIM, CL_PARENT>::GetDerivative, &ngData_class<T, DIM, CL_PARENT>::SetDerivative, "Set the derivative used for AutoDiff")
+    .def("copy", &ngData_class<T, DIM, CL_PARENT>::copy)
+    .def_property_readonly("numEntries", &ngData_class<T, DIM, CL_PARENT>::GetNumEntries, "number of allocated objects of type T")
+    .def_property("vec", &ngData_class<T, DIM, CL_PARENT>::GetVec, &ngData_class<T, DIM, CL_PARENT>::SetVec , "vec")
     
     // destroys convergence
     // .def("__add__", [](ngData_class<T, DIM> *instance, const ngData_class<T, DIM>& o){return instance->operator+(o);})
@@ -375,6 +376,16 @@ void declare_PreisachVectorPilots_t(py::module &m){
 }
 
 
+template<typename T>
+void declare_DummyPilot(py::module &m, std::string name){
+    //return type of pV pilots
+    std::string pyclass_name = std::string("__dummyPilot_") + name + std::string("__");
+    py::class_<dummyPilot<T>, shared_ptr<dummyPilot<T>>>
+    (m, pyclass_name.c_str())
+    ;
+}
+
+
 void ExportReclibBasicClasses(py::module &main_m){
     py::module_ m = main_m.def_submodule("basicClasses", "This is a submodule with all basicClasses and types.");
 
@@ -384,12 +395,18 @@ void ExportReclibBasicClasses(py::module &main_m){
 
     // declare_vvector_t<double>(m, "");
     // declare_vvector_t<std::complex<double>>(m, "complex");
+    declare_DummyPilot<double>(m, "double");
+
 
     declare_ngData_class<scal_t<double>, 1>(m, "Scalar");
     declare_ngData_class<scal_t<double>, 2>(m, "Scalar");
     declare_ngData_class<scal_t<double>, 3>(m, "Scalar");
     declare_ngData_class<scal_t<std::complex<double>>, 3>(m, "ScalarComplex");
 
+    declare_ngData_class<vec_t<double, 1>, 1, ngPreisachScalar<1> >(m, "Vector", "PreisachScalar");
+    declare_ngData_class<vec_t<double, 1>, 2, ngPreisachScalar<2> >(m, "Vector", "PreisachScalar");
+    declare_ngData_class<vec_t<double, 1>, 3, ngPreisachScalar<3> >(m, "Vector", "PreisachScalar");
+
     
 
     // declare_ngData_class<vec_t<double, 1>, 1>(m, "Vector"); //same as <scal_t<double>, 1>
diff --git a/C++/basicClasses/ngData_class.cc b/C++/basicClasses/ngData_class.cc
index d0104755dcd9ae481f14251a2d3b9a8c5ccb088c..810f1a64a9d18852cefc17b53982b96f6d5562c2 100644
--- a/C++/basicClasses/ngData_class.cc
+++ b/C++/basicClasses/ngData_class.cc
@@ -1,14 +1,15 @@
 #include "ngData_class.h"
+#include "../preisach/ngPreisachScalar.h"
 
-
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH>::ngData_class():CoefficientFunction( T::Nr * T::Nc, basicDataTypes::is_complex<typename T::value_type>{}){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT>::ngData_class():CoefficientFunction( T::Nr * T::Nc, basicDataTypes::is_complex<typename T::value_type>{}){
 
 }
 
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH>::ngData_class(shared_ptr<MeshAccess> mesh_in, IntegrationRule intrule_in, const vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>> >& mips_in): 
-                CoefficientFunction( T::Nr * T::Nc, basicDataTypes::is_complex<typename T::value_type>{}), mips{mips_in}, mesh{mesh_in}, intrule{intrule_in}{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT>::ngData_class(shared_ptr<MeshAccess> mesh_in, IntegrationRule intrule_in, const vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>> >& mips_in, 
+                CL_PARENT * parent, typename T::value_type (CL_PARENT::*pilot) (const BaseMappedIntegrationPoint& mip, const typename T::value_type) const): 
+                CoefficientFunction( T::Nr * T::Nc, basicDataTypes::is_complex<typename T::value_type>{}), mips{mips_in}, mesh{mesh_in}, intrule{intrule_in}, parent{parent}, pilot{pilot}{
 
     // cout << __PRETTY_FUNCTION__ << endl;
     if constexpr(T::Nc > 1){
@@ -42,9 +43,10 @@ ngData_class<T, DIM_MESH>::ngData_class(shared_ptr<MeshAccess> mesh_in, Integrat
 };
 
 
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH>::ngData_class(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, const shared_ptr<CoefficientFunction>& ptrMask):
-            CoefficientFunction( T::Nr * T::Nc, basicDataTypes::is_complex<typename T::value_type>{}), mesh{mesh}, intrule{intrule}{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT>::ngData_class(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, const shared_ptr<CoefficientFunction>& ptrMask, 
+            CL_PARENT * parent, typename T::value_type (CL_PARENT::*pilot) (const BaseMappedIntegrationPoint& mip, const typename T::value_type) const):
+            CoefficientFunction( T::Nr * T::Nc, basicDataTypes::is_complex<typename T::value_type>{}), mesh{mesh}, intrule{intrule}, parent{parent}, pilot{pilot}{
 
     // cout << __PRETTY_FUNCTION__ << endl;
 
@@ -63,10 +65,10 @@ ngData_class<T, DIM_MESH>::ngData_class(shared_ptr<MeshAccess>& mesh, Integratio
 }
 
 // copy const
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH>::ngData_class(const ngData_class& o):
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT>::ngData_class(const ngData_class& o):
     CoefficientFunction( T::Nr * T::Nc, basicDataTypes::is_complex<typename T::value_type>{}), 
-    mips{o.mips}, mesh{o.mesh}, intrule{o.intrule}, data{o.data}, num{o.num}, autoUpdate{o.autoUpdate}{            // cout << __PRETTY_FUNCTION__ << endl;
+    mips{o.mips}, mesh{o.mesh}, intrule{o.intrule}, data{o.data}, num{o.num}, autoUpdate{o.autoUpdate}, parent{parent}, pilot{pilot}{            // cout << __PRETTY_FUNCTION__ << endl;
     if constexpr(T::Nc > 1){
         // Set CF dimensions 
         Array<int> dims(2, global_alloc);
@@ -77,8 +79,8 @@ ngData_class<T, DIM_MESH>::ngData_class(const ngData_class& o):
 }
 
 // copy assignment 
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH> ngData_class<T, DIM_MESH>::operator =(const ngData_class& o){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT> ngData_class<T, DIM_MESH, CL_PARENT>::operator =(const ngData_class& o){
             
     // cout << __PRETTY_FUNCTION__ << endl;
     
@@ -100,13 +102,13 @@ ngData_class<T, DIM_MESH> ngData_class<T, DIM_MESH>::operator =(const ngData_cla
     return *this;
 }
 
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH> ngData_class<T, DIM_MESH>::copy() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT> ngData_class<T, DIM_MESH, CL_PARENT>::copy() const{
     return *this;
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::generateMipsFromMesh(const shared_ptr<CoefficientFunction> ptrMask){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::generateMipsFromMesh(const shared_ptr<CoefficientFunction> ptrMask){
     if(mesh == nullptr){
         throw std::runtime_error("mesh has to be set first");
     }
@@ -136,7 +138,7 @@ void ngData_class<T, DIM_MESH>::generateMipsFromMesh(const shared_ptr<Coefficien
             mips[el][ip] = { (intrule)[ip], mesh->GetTrafo(ElementId(VorB, el), global_alloc) };
 
             // mask
-            if(ptrMask != nullptr && ptrMask->Evaluate(mips[el][0]) == 0){
+            if(ptrMask && ptrMask->Evaluate(mips[el][0]) == 0){
                 mips[el].clear();
                 break;
             }
@@ -148,8 +150,8 @@ void ngData_class<T, DIM_MESH>::generateMipsFromMesh(const shared_ptr<Coefficien
 }
         
         // operators
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH>& ngData_class<T, DIM_MESH>::operator +=(const ngData_class& o){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT>& ngData_class<T, DIM_MESH, CL_PARENT>::operator +=(const ngData_class& o){
     unsigned int el, ip;
     for(el = 0; el != o.mips.size(); ++el){
         for (ip = 0;ip < o.mips[el].size(); ip++) {
@@ -160,9 +162,8 @@ ngData_class<T, DIM_MESH>& ngData_class<T, DIM_MESH>::operator +=(const ngData_c
     return *this;   
 }
 
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH>& ngData_class<T, DIM_MESH>::operator -=(const ngData_class& o){    
-    cout << "sub" << endl;
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT>& ngData_class<T, DIM_MESH, CL_PARENT>::operator -=(const ngData_class& o){    
     unsigned int el, ip;
     for(el = 0; el != o.mips.size(); ++el){
         for (ip = 0;ip < o.mips[el].size(); ip++) {
@@ -172,15 +173,15 @@ ngData_class<T, DIM_MESH>& ngData_class<T, DIM_MESH>::operator -=(const ngData_c
     return *this;   
 }
 
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH> ngData_class<T, DIM_MESH>::operator +(const ngData_class& o) const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT> ngData_class<T, DIM_MESH, CL_PARENT>::operator +(const ngData_class& o) const{
     ngData_class ret(*this);
     ret += o;
     return ret;
 }
 
-template<typename T, unsigned int DIM_MESH>
-ngData_class<T, DIM_MESH> ngData_class<T, DIM_MESH>::operator -(const ngData_class& o) const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+ngData_class<T, DIM_MESH, CL_PARENT> ngData_class<T, DIM_MESH, CL_PARENT>::operator -(const ngData_class& o) const{
     ngData_class ret(*this);
     ret -= o;
     return ret;
@@ -203,35 +204,35 @@ ngData_class<T, DIM_MESH> ngData_class<T, DIM_MESH>::operator -(const ngData_cla
         // }
 
 
-template<typename T, unsigned int DIM_MESH>
-std::map<std::pair<const unsigned int, const unsigned int>, T>& ngData_class<T, DIM_MESH>::getData(){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+std::map<std::pair<const unsigned int, const unsigned int>, T>& ngData_class<T, DIM_MESH, CL_PARENT>::getData(){
     return data;
     
 }
 
-template<typename T, unsigned int DIM_MESH>
-std::map<std::pair<const unsigned int, const unsigned int>, T> ngData_class<T, DIM_MESH>::getDataConst() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+std::map<std::pair<const unsigned int, const unsigned int>, T> ngData_class<T, DIM_MESH, CL_PARENT>::getDataConst() const{
     return data;
     
 }
 
-template<typename T, unsigned int DIM_MESH>
-vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>>>& ngData_class<T, DIM_MESH>::getMips(){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>> >& ngData_class<T, DIM_MESH, CL_PARENT>::getMips(){
     return mips;
 }
 
-template<typename T, unsigned int DIM_MESH>
-shared_ptr<MeshAccess>& ngData_class<T, DIM_MESH>::getMesh(){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+shared_ptr<MeshAccess>& ngData_class<T, DIM_MESH, CL_PARENT>::getMesh(){
     return mesh;
 }
 
-template<typename T, unsigned int DIM_MESH>
-IntegrationRule& ngData_class<T, DIM_MESH>::getIntrule(){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+IntegrationRule& ngData_class<T, DIM_MESH, CL_PARENT>::getIntrule(){
     return intrule;
 }
 
-template<typename T, unsigned int DIM_MESH>
-vector<vector<double>> ngData_class<T, DIM_MESH>::getMipsAsMatrix(){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+vector<vector<double>> ngData_class<T, DIM_MESH, CL_PARENT>::getMipsAsMatrix(){
     vector<vector<double>> tmp;
     unsigned int el, ip;
 
@@ -257,8 +258,8 @@ vector<vector<double>> ngData_class<T, DIM_MESH>::getMipsAsMatrix(){
 }
 
 
-template<typename T, unsigned int DIM_MESH>
-auto ngData_class<T, DIM_MESH>::getState() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+auto ngData_class<T, DIM_MESH, CL_PARENT>::getState() const{
     throw std::runtime_error("not implemented");
     // return std::make_tuple(, instance->getPtrMask(), instance->getNum());
     return std::make_tuple(data, mesh);
@@ -266,16 +267,16 @@ auto ngData_class<T, DIM_MESH>::getState() const{
 
 }
 
-template<typename T, unsigned int DIM_MESH>
-T& ngData_class<T, DIM_MESH>::serialElement(const unsigned int ind){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+T& ngData_class<T, DIM_MESH, CL_PARENT>::serialElement(const unsigned int ind){
     // data per element = T::Nc * T::Nr
     auto it = data.begin();
     std::advance(it, ind);
     return it->second;
 }
 
-template<typename T, unsigned int DIM_MESH>
-const T ngData_class<T, DIM_MESH>::Get(const unsigned int el, const unsigned int ip) const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+const T ngData_class<T, DIM_MESH, CL_PARENT>::Get(const unsigned int el, const unsigned int ip) const{
     if( el < mips.size() && ip < mips[el].size() ){
         return data.at(std::make_pair(el, ip));
     }
@@ -285,8 +286,8 @@ const T ngData_class<T, DIM_MESH>::Get(const unsigned int el, const unsigned int
     }
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Set(const unsigned int el, const unsigned int ip, const T& val){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Set(const unsigned int el, const unsigned int ip, const T& val){
     if( el < mips.size() && ip < mips[el].size() ){
         data[std::make_pair(el, ip)] = val;
     }
@@ -296,8 +297,8 @@ void ngData_class<T, DIM_MESH>::Set(const unsigned int el, const unsigned int ip
 }
 
 // special case for double
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Set(const unsigned int el, const unsigned int ip, const double& val){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Set(const unsigned int el, const unsigned int ip, const double& val){
     if( el < mips.size() && ip < mips[el].size() ){
         data[std::make_pair(el, ip)][0][0] = val;
     }
@@ -306,8 +307,8 @@ void ngData_class<T, DIM_MESH>::Set(const unsigned int el, const unsigned int ip
     }
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Set(const shared_ptr<CoefficientFunction>& ptrIn){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Set(const shared_ptr<CoefficientFunction>& ptrIn){
     parFor(mips.size(), [&](int el) {
         unsigned int ip, i, l;
         typename T::value_type tmp[T::Nr * T::Nc];
@@ -349,13 +350,13 @@ void ngData_class<T, DIM_MESH>::Set(const shared_ptr<CoefficientFunction>& ptrIn
     
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Set(const ngData_class<T, DIM_MESH>& ptrIn){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Set(const ngData_class<T, DIM_MESH, CL_PARENT>& ptrIn){
     this->operator=(ptrIn);
 }
 
-template<typename T, unsigned int DIM_MESH>
-std::string ngData_class<T, DIM_MESH>::toString() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+std::string ngData_class<T, DIM_MESH, CL_PARENT>::toString() const{
     std::stringstream s;
 
     s << *this;
@@ -363,8 +364,8 @@ std::string ngData_class<T, DIM_MESH>::toString() const{
     return s.str();
 }
 
-template<typename T, unsigned int DIM_MESH>
-unsigned int ngData_class<T, DIM_MESH>::GetNumBytes() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+unsigned int ngData_class<T, DIM_MESH, CL_PARENT>::GetNumBytes() const{
     
     if constexpr(!basicDataTypes::is_vvector_t<T>{}){
         return num * sizeof(typename T::value_type) * T::Nr * T::Nc;
@@ -377,18 +378,18 @@ unsigned int ngData_class<T, DIM_MESH>::GetNumBytes() const{
         
 } 
 
-template<typename T, unsigned int DIM_MESH>
-unsigned int ngData_class<T, DIM_MESH>::GetNumEntries() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+unsigned int ngData_class<T, DIM_MESH, CL_PARENT>::GetNumEntries() const{
     return num;
 }
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::SetNumEntries(unsigned int o){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::SetNumEntries(unsigned int o){
     num = o;
 }
 
 // for now, just around xyz-axis
-template<typename T, unsigned int DIM_MESH>
-std::shared_ptr<ngData_class<T, DIM_MESH>> ngData_class<T, DIM_MESH>::GetRotatedValues(const double phi, unsigned int axis){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+std::shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> ngData_class<T, DIM_MESH, CL_PARENT>::GetRotatedValues(const double phi, unsigned int axis){
 
     // if constexpr(basicDataTypes::is_complex<typename T::value_type>{}){
     //     throw std::runtime_error("not implemented for complex values yet");
@@ -414,7 +415,7 @@ std::shared_ptr<ngData_class<T, DIM_MESH>> ngData_class<T, DIM_MESH>::GetRotated
         R = mat_t<double, DIM_MESH, DIM_MESH>{1};
     }
 
-    std::shared_ptr<ngData_class<T, DIM_MESH>> data_ret = std::make_shared<ngData_class<T, DIM_MESH>>(mesh, intrule, mips);
+    std::shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> data_ret = std::make_shared<ngData_class<T, DIM_MESH, CL_PARENT>>(mesh, intrule, mips);
 
 
 
@@ -482,30 +483,30 @@ std::shared_ptr<ngData_class<T, DIM_MESH>> ngData_class<T, DIM_MESH>::GetRotated
 
 }
 
-template<typename T, unsigned int DIM_MESH>
-bool ngData_class<T, DIM_MESH>::GetAutoUpdate() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+bool ngData_class<T, DIM_MESH, CL_PARENT>::GetAutoUpdate() const{
     return autoUpdate;
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::SetAutoUpdate(const bool val){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::SetAutoUpdate(const bool val){
     autoUpdate = val;
 }
 
-template<typename T, unsigned int DIM_MESH>
-shared_ptr<ngData_class<T, DIM_MESH>> ngData_class<T, DIM_MESH>::GetDerivative() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> ngData_class<T, DIM_MESH, CL_PARENT>::GetDerivative() const{
     return derivative;
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::SetDerivative(const shared_ptr<ngData_class<T, DIM_MESH>> val){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::SetDerivative(const shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> val){
     derivative = val;
 }
 
 
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Evaluate(const BaseMappedIntegrationPoint& ip, FlatVector<Complex> result) const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Evaluate(const BaseMappedIntegrationPoint& ip, FlatVector<Complex> result) const{
 
     if constexpr(basicDataTypes::is_complex<typename T::value_type>{}){
         double tmp[2 * (T::Nc * T::Nr)];
@@ -534,8 +535,8 @@ void ngData_class<T, DIM_MESH>::Evaluate(const BaseMappedIntegrationPoint& ip, F
 
 
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Evaluate(const BaseMappedIntegrationPoint& mip,  FlatVector<> result) const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Evaluate(const BaseMappedIntegrationPoint& mip,  FlatVector<> result) const{
 
     // cout << __PRETTY_FUNCTION__ << " enter" << endl;
     if constexpr(!basicDataTypes::is_vvector_t<T>{}){
@@ -581,8 +582,8 @@ void ngData_class<T, DIM_MESH>::Evaluate(const BaseMappedIntegrationPoint& mip,
     }
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Evaluate (const BaseMappedIntegrationRule & ir, 
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Evaluate (const BaseMappedIntegrationRule & ir, 
                         BareSliceMatrix<AutoDiff<1,double>> values) const{
 
             // typename T::value_type tmp[this->Dimension()] = {};
@@ -613,16 +614,22 @@ void ngData_class<T, DIM_MESH>::Evaluate (const BaseMappedIntegrationRule & ir,
                         // values(i, j).Value(0) = this->Evaluate(ir[i]);
 
                         // double eps = 1e-5;
-                        double val = this->Evaluate(ir[i]);
+                        double val=0;
+                        if (parent && pilot){
+
+                            // cout << values(i, j).Value() << " " << values(i, j).DValue(0) << endl;
+                            val = (parent->*pilot)(ir[i], values(i, j).Value() + values(i, j).DValue(0));
+                            // this->Evaluate(ir[i]);
+                        }
 
                         
-                        double dval2;
-                        if (derivative == nullptr){
-                            dval2 = this->Evaluate(ir[i]);
-                            cout << "warning: no derivative set for dataClass" << endl;
+                        double dval2=0;
+                        if (derivative){
+                            dval2 = derivative->Evaluate(ir[i]);
+
                         }
                         else{
-                            dval2 = derivative->Evaluate(ir[i]);
+                            cout << "warning: no derivative set for dataClass" << endl;
                         }
                         // cout << "dval = " << dval << " =?= " << dval2 << endl;
 
@@ -639,8 +646,8 @@ void ngData_class<T, DIM_MESH>::Evaluate (const BaseMappedIntegrationRule & ir,
                 }
             }        
 }
-template<typename T, unsigned int DIM_MESH>
-double ngData_class<T, DIM_MESH>::Evaluate(const BaseMappedIntegrationPoint& mip) const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+double ngData_class<T, DIM_MESH, CL_PARENT>::Evaluate(const BaseMappedIntegrationPoint& mip) const{
     // cout << __PRETTY_FUNCTION__ << " enter" << endl;
 
     if constexpr(!basicDataTypes::is_vvector_t<T>{}){
@@ -654,8 +661,8 @@ double ngData_class<T, DIM_MESH>::Evaluate(const BaseMappedIntegrationPoint& mip
     }
 }
 
-template<typename T, unsigned int DIM_MESH>
-typename T::value_type ngData_class<T, DIM_MESH>::Sum() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+typename T::value_type ngData_class<T, DIM_MESH, CL_PARENT>::Sum() const{
     typename T::value_type sum{};
     for (int el = 0; el < mips.size(); el++) {
         // create a container for every integration point
@@ -675,8 +682,8 @@ typename T::value_type ngData_class<T, DIM_MESH>::Sum() const{
     return sum;
 }
 
-template<typename T, unsigned int DIM_MESH>
-std::vector<typename T::value_type>  ngData_class<T, DIM_MESH>::GetVec() const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+std::vector<typename T::value_type>  ngData_class<T, DIM_MESH, CL_PARENT>::GetVec() const{
     std::vector<typename T::value_type> ret{};
     for (int el = 0; el < mips.size(); el++) {
         // create a container for every integration point
@@ -701,8 +708,8 @@ std::vector<typename T::value_type>  ngData_class<T, DIM_MESH>::GetVec() const{
     return ret;
 }
 
-template<typename T, unsigned int DIM_MESH>
-void  ngData_class<T, DIM_MESH>::SetVec(const std::vector<typename T::value_type>& val){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void  ngData_class<T, DIM_MESH, CL_PARENT>::SetVec(const std::vector<typename T::value_type>& val){
     
     int n = 0;
     for (int el = 0; el < mips.size(); el++) {
@@ -729,15 +736,15 @@ void  ngData_class<T, DIM_MESH>::SetVec(const std::vector<typename T::value_type
     }
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::setData(std::map<std::pair<const unsigned int, const unsigned int>, T> o){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::setData(std::map<std::pair<const unsigned int, const unsigned int>, T> o){
     data = o;
 }
 
 
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Save(std::string fileName, bool binary) const{
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Save(std::string fileName, bool binary) const{
 
     std::ofstream file(fileName, std::ios::out | std::ios::binary) ;
     if (file.is_open()){
@@ -774,8 +781,8 @@ void ngData_class<T, DIM_MESH>::Save(std::string fileName, bool binary) const{
     }
 }
 
-template<typename T, unsigned int DIM_MESH>
-void ngData_class<T, DIM_MESH>::Load(std::string filename, const unsigned int iSkip, bool binary){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+void ngData_class<T, DIM_MESH, CL_PARENT>::Load(std::string filename, const unsigned int iSkip, bool binary){
     
     if(!binary){
         size_t pos = 0;
@@ -960,8 +967,8 @@ void ngData_class<T, DIM_MESH>::Load(std::string filename, const unsigned int iS
 
 
 
-template<typename T, unsigned int DIM_MESH>
-std::ostream& operator << (std::ostream& o, const ngData_class<T, DIM_MESH>& m){
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT>
+std::ostream& operator << (std::ostream& o, const ngData_class<T, DIM_MESH, CL_PARENT>& m){
 
     unsigned int el, ip, i, l;
 
@@ -1011,6 +1018,13 @@ template class ngData_class<basicDataTypes::vec_t<double, 2>, 2>;
 template class ngData_class<basicDataTypes::vec_t<double, 2>, 1>;  // 2D VPM for 1D Mesh
 template class ngData_class<basicDataTypes::vec_t<double, 3>, 3>;
 template class ngData_class<basicDataTypes::vec_t<double, 2>, 3>;  // 2D VPM for 3D Mesh
+// new
+template class ngData_class<vec_t<double, 1>, 1, ngPreisachScalar<1> >;
+template class ngData_class<vec_t<double, 1>, 2, ngPreisachScalar<2> >;
+template class ngData_class<vec_t<double, 1>, 3, ngPreisachScalar<3> >;
+// template class ngData_class<vec_t<double, 1>, 3, ngPreisachScalar<3> >;
+
+
 
 template class ngData_class<basicDataTypes::mat_t<double, 1, 1>, 1>;
 template class ngData_class<basicDataTypes::mat_t<double, 1, 1>, 2>;
diff --git a/C++/basicClasses/ngData_class.h b/C++/basicClasses/ngData_class.h
index 6ef54a3b786c2ddb310e006df1eaf1edddcc0eae..17a358a06f26c4d9759993b0acf905cd64b53c6a 100644
--- a/C++/basicClasses/ngData_class.h
+++ b/C++/basicClasses/ngData_class.h
@@ -30,12 +30,22 @@ using basicDataTypes::vec_t;
 using basicDataTypes::vvector_t;
 using basicDataTypes::scal_t;
 
+template<typename T>
+class dummyPilot{
+    public:
+    dummyPilot(){};
+    T pilot(const BaseMappedIntegrationPoint& mip, const T){
+        return T{};
+    }
+};
+
 
-template<typename T, unsigned int DIM_MESH>
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT=dummyPilot<typename T::value_type>>
 class ngData_class : public CoefficientFunction{
     protected:
         // data
         vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>> > mips;
+        // std::map<size_t, MappedIntegrationPoint<DIM_MESH, DIM_MESH> > mips;
         shared_ptr<MeshAccess> mesh;
         IntegrationRule intrule;
         // shared_ptr<CoefficientFunction> ptrMask;
@@ -48,15 +58,20 @@ class ngData_class : public CoefficientFunction{
         std::map<std::pair<const unsigned int, const unsigned int>, T> data;
         unsigned int num = 0;
 
-        shared_ptr<ngData_class<T, DIM_MESH>> derivative = nullptr;
+        shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> derivative = nullptr;
 
 
+        CL_PARENT* parent = nullptr;
+        typename T::value_type (CL_PARENT::*pilot) (const BaseMappedIntegrationPoint& mip, const typename T::value_type) const= nullptr;
+
         
     public:
         ngData_class();
-        ngData_class(shared_ptr<MeshAccess> mesh_in, IntegrationRule intrule_in, const vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>> >& mips_in);
+        ngData_class(shared_ptr<MeshAccess> mesh_in, IntegrationRule intrule_in, const vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>> >& mips_in, 
+         CL_PARENT * parent = nullptr, typename T::value_type (CL_PARENT::*pilot) (const BaseMappedIntegrationPoint& mip, const typename T::value_type) const = nullptr);
 
-        ngData_class(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, const shared_ptr<CoefficientFunction>& ptrMask=nullptr);
+        ngData_class(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, const shared_ptr<CoefficientFunction>& ptrMask=nullptr,
+         CL_PARENT * parent = nullptr, typename T::value_type (CL_PARENT::*pilot) (const BaseMappedIntegrationPoint& mip, const typename T::value_type) const = nullptr);
 
         // copy const
         ngData_class(const ngData_class& o);
@@ -98,7 +113,7 @@ class ngData_class : public CoefficientFunction{
 
         std::map<std::pair<const unsigned int, const unsigned int>, T>& getData();
         std::map<std::pair<const unsigned int, const unsigned int>, T> getDataConst() const;
-        vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>>>& getMips();
+        vector<vector<MappedIntegrationPoint<DIM_MESH, DIM_MESH>> >& getMips();
 
         shared_ptr<MeshAccess>& getMesh();
 
@@ -132,17 +147,17 @@ class ngData_class : public CoefficientFunction{
         void SetNumEntries(unsigned int o);
 
         // for now, just around xyz-axis
-        std::shared_ptr<ngData_class<T, DIM_MESH>> GetRotatedValues(const double phi, unsigned int axis=2);
+        std::shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> GetRotatedValues(const double phi, unsigned int axis=2);
 
 
         [[nodiscard]] bool GetAutoUpdate() const;
         void SetAutoUpdate(const bool val);
 
-        [[nodiscard]] shared_ptr<ngData_class<T, DIM_MESH>> GetDerivative() const;
-        void SetDerivative(const shared_ptr<ngData_class<T, DIM_MESH>> val);
+        [[nodiscard]] shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> GetDerivative() const;
+        void SetDerivative(const shared_ptr<ngData_class<T, DIM_MESH, CL_PARENT>> val);
 
-        template<typename iT, unsigned int iDIM_MESH>
-        friend std::ostream& operator << (std::ostream& o, const ngData_class<iT, iDIM_MESH>& m);
+        template<typename iT, unsigned int iDIM_MESH, typename iCL_PARENT>
+        friend std::ostream& operator << (std::ostream& o, const ngData_class<iT, iDIM_MESH, iCL_PARENT>& m);
         
 
         virtual void Evaluate(const BaseMappedIntegrationPoint& ip, FlatVector<Complex> result) const;
@@ -169,12 +184,13 @@ class ngData_class : public CoefficientFunction{
 };
 
 
-template<typename T, unsigned int DIM>
-using ngScalar_class = ngData_class<scal_t<T>, DIM>;
-template<typename T, unsigned int DIM>
-using ngVector_class = ngData_class<vec_t<T, DIM>, DIM>;
-template<typename T, unsigned int DIM>
-using ngMat_class = ngData_class<mat_t<T, DIM, DIM>, DIM>;
+
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT=dummyPilot<T>>
+using ngScalar_class = ngData_class<scal_t<T>, DIM_MESH, CL_PARENT>;
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT=dummyPilot<T>>
+using ngVector_class = ngData_class<vec_t<T, DIM_MESH>, DIM_MESH,  CL_PARENT>;
+template<typename T, unsigned int DIM_MESH, typename CL_PARENT=dummyPilot<T>>
+using ngMat_class = ngData_class<mat_t<T, DIM_MESH, DIM_MESH>, DIM_MESH,  CL_PARENT>;
 
 
 
diff --git a/C++/libMain.h b/C++/libMain.h
index 6f7fbf27fc57e23dbd5c5fbcabfe0277189aeb16..a93ced5cd274d1c842a222a19179b95936938e74 100644
--- a/C++/libMain.h
+++ b/C++/libMain.h
@@ -36,6 +36,9 @@ void ExportReclibPhiFunctions(py::module &m);
 
 #ifdef CMAKE_BASICCLASSES
 void ExportReclibBasicClasses(py::module &m);
+
+template<typename T, unsigned int DIM, typename CL_PARENT=dummyPilot<typename T::value_type>>
+void declare_ngData_class(py::module &m, std::string name, std::string suffix="");
 #endif
 
 #ifdef CMAKE_NONLINEAR
diff --git a/C++/nonLin/nonLinCF.cc b/C++/nonLin/nonLinCF.cc
index ea838e18b57420f9dacb47c58406a77958a825f8..b4969cc324a85a723ffa876f063c8ea44f6dad54 100644
--- a/C++/nonLin/nonLinCF.cc
+++ b/C++/nonLin/nonLinCF.cc
@@ -22,7 +22,7 @@ nonLinCF<T, DIM>::nonLinCF(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrul
     oKL = std::make_shared<KL>(H_KL, B_KL, order, true);
     findConcavePoints();
     // cout << "order!" << endl;
-    Update();
+    // Update();
 }; 
 
 template <typename T, unsigned int DIM>
@@ -40,7 +40,7 @@ nonLinCF<T, DIM>::nonLinCF(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrul
     oKL = std::make_shared<KL>(H_KL, B_KL, order, true);
     findConcavePoints();
     // cout << "order!" << endl;
-    Update();
+    // Update();
 }; 
 
 
@@ -363,33 +363,23 @@ nonLinCF<T, DIM> nonLinCF<T, DIM>::Differntiate(){
 template <typename T, unsigned int DIM>
 AutoDiff<1,double> nonLinCF<T, DIM>:: operator() (AutoDiff<1,double> x) const
   {
-    // double eps = 1e-5;
-    double val = this->oKL->EvaluateH(x.Value());
-    // double valr = this->oKL->EvaluateH(x.Value()+eps);
-    // double vall = this->oKL->EvaluateH(x.Value()-eps);
-    // double dval = (valr-vall) / (2*eps);
 
-    if(strcmp(sFieldOut.c_str(), "B")){
-        throw std::runtime_error("to use NewtonMetho use H_KL, mu_KL and sFieldOut=B ");
-    }
+    // cout << x.Value() << " " << x.DValue(0) << endl;
 
+    double x_in = x.Value() + x.DValue(0);
+    double val = this->oKL->EvaluateH(x_in);
     
-    // cout << "x " << x.Value() << " " << x.DValue(0) << endl;
+    if(strcmp(sFieldOut.c_str(), "B")){
+        throw std::runtime_error("to use NewtonMethod use H_KL, mu_KL and sFieldOut=B ");
+    }
 
-    double dval2 = this->oKL->DifferentiateH(x.Value());
-    // cout << "dval = " << dval << " =?= " << dval2 << endl;
 
-    AutoDiff<1,double> res(val);
-    // if (x.DValue(0) != 0){
-    //     res.DValue(0) = dval2;// * x.DValue(0);
-    // }
-    // else{
-    res.DValue(0) = dval2;
-    // }
+    double dval = this->oKL->DifferentiateH(x_in);
 
     
-    
-    // cout << "r " << res.Value() << " " << res.DValue(0) << endl;
+
+    AutoDiff<1,double> res(val);
+    res.DValue(0) = dval * x.DValue(0);
 
 
     return res;
diff --git a/C++/preisach/libPreisachScalar.cc b/C++/preisach/libPreisachScalar.cc
index 89089c9035468527799dac33d28e09e2b2fd30fc..e6a31827e3fc9edc15a29692006a061b893f41ee 100644
--- a/C++/preisach/libPreisachScalar.cc
+++ b/C++/preisach/libPreisachScalar.cc
@@ -97,6 +97,9 @@ void declare_ngPreisachScalar(py::module &m) {
 
     ;
 
+    // new
+    // declare_ngData_class<vec_t<double, 1>, DIM, ngPreisachScalar<DIM> >(m, "Vector", "PreisachScalar");
+
 }
 
 
@@ -115,6 +118,10 @@ void ExportReclibPreisachScalar(py::module &main_m){
     declare_ngPreisachScalar<2>(m);
     declare_ngPreisachScalar<3>(m);
 
+
+
+
+
     // --------------------------------------------------------------------------- *
     //  -----------   Fundamental Preisach Operators (no CF)                        *
     //  --------------------------------------------------------------------------- 
diff --git a/C++/preisach/ngPreisachScalar.cc b/C++/preisach/ngPreisachScalar.cc
index d8f2b148685973be54a8b8f2cfdac3afdd123df6..7dacb7b720d1017a800e417ee1ae28d729dfef3c 100644
--- a/C++/preisach/ngPreisachScalar.cc
+++ b/C++/preisach/ngPreisachScalar.cc
@@ -204,6 +204,8 @@ void ngPreisachScalar<DIM>::UpdatePast(const shared_ptr<CoefficientFunction>& pt
 
 }
 
+
+
 template <unsigned int DIM>
 void ngPreisachScalar<DIM>::Update(){  
     Update(ptrInput);
@@ -738,11 +740,11 @@ const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetReactive
 }
 
 template <unsigned int DIM>
-const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetB(const bool autoupdate){
+const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ngPreisachScalar<DIM>::GetB(const bool autoupdate){
     //----------- MEMORY DEMANDING ---------------
     // flux density
     if(ptrFluxDens == nullptr){
-        ptrFluxDens = std::make_shared<ngScalar_class<double, DIM>>(mesh, intrule, mips);
+        ptrFluxDens = std::make_shared<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>>(mesh, intrule, mips, this, &ngPreisachScalar<DIM>::pilotForward<TYPES_PILOT::OUTPUT>);
         ptrFluxDens->SetAutoUpdate(autoupdate);
         Update();
     }
@@ -750,30 +752,30 @@ const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetB(const
 }
 
 template <unsigned int DIM>
-const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetH(const bool autoupdate){
+const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ngPreisachScalar<DIM>::GetH(const bool autoupdate){
     // field strength
     if(ptrFieldStrength == nullptr){
-        ptrFieldStrength = std::make_shared<ngScalar_class<double, DIM>>(mesh, intrule, mips);
+        ptrFieldStrength = std::make_shared<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>>(mesh, intrule, mips, this, &ngPreisachScalar<DIM>::pilotForward<TYPES_PILOT::INPUT>);
         ptrFieldStrength->SetAutoUpdate(autoupdate);
         Update();
     }
     return ptrFieldStrength;
 }
 template <unsigned int DIM>
-const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetMu(const bool autoupdate){
+const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ngPreisachScalar<DIM>::GetMu(const bool autoupdate){
     // diff_mu
     if(ptrMu == nullptr){
-        ptrMu = std::make_shared<ngScalar_class<double, DIM>>(mesh, intrule, mips);
+        ptrMu = std::make_shared<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>>(mesh, intrule, mips, this, &ngPreisachScalar<DIM>::pilotForward<TYPES_PILOT::MAT>);
         ptrMu->SetAutoUpdate(autoupdate);
         Update();
     }
     return ptrMu;
 }
 template <unsigned int DIM>
-const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetMuDiff(const bool autoupdate){
+const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ngPreisachScalar<DIM>::GetMuDiff(const bool autoupdate){
     // diff_mu
     if(ptrMuDiff == nullptr){
-        ptrMuDiff = std::make_shared<ngScalar_class<double, DIM>>(mesh, intrule, mips);
+        ptrMuDiff = std::make_shared<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>>(mesh, intrule, mips, this, &ngPreisachScalar<DIM>::pilotForward<TYPES_PILOT::MAT_DIFF>);
         ptrMuDiff->SetAutoUpdate(autoupdate);
         Update();
     }
@@ -781,10 +783,10 @@ const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetMuDiff(c
 }
 
 template <unsigned int DIM>
-const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetMatDiffDiff(const bool autoupdate){
+const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ngPreisachScalar<DIM>::GetMatDiffDiff(const bool autoupdate){
     // diffdiff_mu
     if(ptrMatDiffDiff == nullptr){
-        ptrMatDiffDiff = std::make_shared<ngScalar_class<double, DIM>>(mesh, intrule, mips);
+        ptrMatDiffDiff = std::make_shared<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>>(mesh, intrule, mips, this, &ngPreisachScalar<DIM>::pilotForward<TYPES_PILOT::MAT_DIFFDIFF>);
         ptrMatDiffDiff->SetAutoUpdate(true);        
         Update();
 
@@ -858,6 +860,9 @@ unsigned int ngPreisachScalar<DIM>::GetNumBytes()const{
     return ret;
 }
 
+
+
+
 template class ngPreisachScalar<1>;
 template class ngPreisachScalar<2>;
 template class ngPreisachScalar<3>;
diff --git a/C++/preisach/ngPreisachScalar.h b/C++/preisach/ngPreisachScalar.h
index 3ed1cb69a9da3dd3f0cf2dde29cca1208023c3af..ce0e8e070ba2524519a5a25aafff4cbd7c2dbc49 100644
--- a/C++/preisach/ngPreisachScalar.h
+++ b/C++/preisach/ngPreisachScalar.h
@@ -28,6 +28,13 @@ using namespace ngcomp;
 
 using preisachScalar_ns::preisach;
 
+enum class TYPES_PILOT { 
+                    INPUT,
+                    OUTPUT,
+                    MAT, 
+                    MAT_DIFF, 
+                    MAT_DIFFDIFF };
+
 
 template<unsigned int DIM>
 class ngPreisachScalar:public CoefficientFunction
@@ -51,11 +58,11 @@ class ngPreisachScalar:public CoefficientFunction
         ~ngPreisachScalar() = default;
 
 
-        const shared_ptr<ngScalar_class<double, DIM>> GetB(const bool autoupdate=true);
-        const shared_ptr<ngScalar_class<double, DIM>> GetH(const bool autoupdate=true);
-        const shared_ptr<ngScalar_class<double, DIM>> GetMu(const bool autoupdate=true);
-        const shared_ptr<ngScalar_class<double, DIM>> GetMuDiff(const bool autoupdate=true);
-        const shared_ptr<ngScalar_class<double, DIM>> GetMatDiffDiff(const bool autoupdate=true);
+        const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> GetB(const bool autoupdate=true);
+        const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> GetH(const bool autoupdate=true);
+        const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> GetMu(const bool autoupdate=true);
+        const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> GetMuDiff(const bool autoupdate=true);
+        const shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> GetMatDiffDiff(const bool autoupdate=true);
         const shared_ptr<ngScalar_class<double, DIM>> GetNu(const bool autoupdate=true);
         const shared_ptr<ngScalar_class<double, DIM>> GetNuDiff(const bool autoupdate=true);
         const shared_ptr<ngScalar_class<double, DIM>> GetEnergyDensity(const bool autoupdate=false);
@@ -70,6 +77,31 @@ class ngPreisachScalar:public CoefficientFunction
         void Demagnetise();
         bool HasConverged(double aMeanError, double aMaxError) const;
 
+
+        template<TYPES_PILOT TYPE_P>
+        double pilotForward (const BaseMappedIntegrationPoint& mip, const double val) const{
+            // cout << __PRETTY_FUNCTION__ << endl;
+
+            auto pilot = GetPreisach(mip)->pilotForward(val);
+            
+            if constexpr(TYPE_P == TYPES_PILOT::OUTPUT){
+                return pilot->getOutputVal_interpolated();
+            }
+            else if constexpr(TYPE_P == TYPES_PILOT::INPUT){
+                return pilot->getInputVal_interpolated();
+            }
+            else if constexpr(TYPE_P == TYPES_PILOT::MAT){
+                return pilot->getMat();
+            }
+            else if constexpr(TYPE_P == TYPES_PILOT::MAT_DIFF){
+                return pilot->getMatDiff();
+            }
+            else if constexpr(TYPE_P == TYPES_PILOT::MAT_DIFFDIFF){
+                return pilot->calcMatDiffDiff(true);
+            }
+            return 0;
+        }
+
         void Update();
         void Update(const shared_ptr<CoefficientFunction>& ptr);
         void UpdateEnergyDensity();
@@ -133,15 +165,16 @@ class ngPreisachScalar:public CoefficientFunction
         
         vector<vector<MappedIntegrationPoint<DIM, DIM>> > mips;
         vector<vector<shared_ptr<preisach> > > preisachScalarPoints;
+        // std::unordered_map<*BaseMappedIntegrationPoint, shared_ptr<preisach>>  preisachScalarPointsTes;
         std::unordered_map<unsigned int, shared_ptr<preisach>> preisachScalarContainer;
         std::unordered_map<shared_ptr<preisach>, std::pair<unsigned int, unsigned int> > elIpResponsibleForUpdatingPreisach;
         shared_ptr<EverettMatrix> ptrEverettMatrix;
 
-        shared_ptr<ngScalar_class<double, DIM>> ptrFluxDens;
-        shared_ptr<ngScalar_class<double, DIM>> ptrFieldStrength;
-        shared_ptr<ngScalar_class<double, DIM>> ptrMu;
-        shared_ptr<ngScalar_class<double, DIM>> ptrMuDiff;
-        shared_ptr<ngScalar_class<double, DIM>> ptrMatDiffDiff;
+        shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ptrFluxDens;
+        shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ptrFieldStrength;
+        shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ptrMu;
+        shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ptrMuDiff;
+        shared_ptr<ngScalar_class<double, DIM, ngPreisachScalar<DIM>>> ptrMatDiffDiff;
         shared_ptr<ngScalar_class<double, DIM>> ptrNu;
         shared_ptr<ngScalar_class<double, DIM>> ptrNuDiff;
         shared_ptr<ngScalar_class<double, DIM>> ptrEnergyDens;
diff --git a/C++/preisach/ngPreisachVector.cc b/C++/preisach/ngPreisachVector.cc
index 2fca1f0616fb4e04b501c0c4dcfeb2ff6836a499..17b93bfad2181e4376ac90a4c6d47dca101b48c0 100644
--- a/C++/preisach/ngPreisachVector.cc
+++ b/C++/preisach/ngPreisachVector.cc
@@ -266,7 +266,7 @@ void ngPreisachVector<DIM_MESH, DIM_PREISACH>::Update(const shared_ptr<Coefficie
                         // evaluate with FlatVector
                         ptrIn->Evaluate(mips[el][ip], fV_in);
 
-                        cout << el<< " " << ip << " " <<tmp[0] << ", " << tmp[1] << endl;
+                        // cout << el<< " " << ip << " " <<tmp[0] << ", " << tmp[1] << endl;
                         // cast it to vec_t
                         UNUSED(preisachVectorPoints[el][ip]->pilotForward(vec_t<double, DIM_PREISACH>(tmp, DIM_PREISACH)));