added harmonic balance to nonLinCF

C++/nonLin/libNonLin.cc

@@ -22,9 +22,7 @@
 
 
 template<typename T, unsigned int DIM>
-void declare_nonLinCF(py::module &main_m){
-    py::module_ m = main_m.def_submodule("nonLin", "nonLinear CoefficientFunction");
-
+void declare_nonLinCF(py::module &m){
 
     std::string pyclass_name = std::string("nonLinCF") + std::to_string(DIM);
 
@@ -69,16 +67,52 @@ void declare_nonLinCF(py::module &main_m){
 }
 
 
-void ExportReclibNonLin(py::module &m)
+template<typename T, unsigned int DIM>
+void declare_nonLinCF_harmonicBalance(py::module &m){
+
+
+    std::string pyclass_name = std::string("nonLinCF_harmonicBalance") + std::to_string(DIM);
+
+    py::class_<nonLinCF_harmonicBalance<T, DIM>, shared_ptr<nonLinCF_harmonicBalance<T, DIM>>, nonLinCF<T, DIM>>
+    (m, pyclass_name.c_str(), "nonLinear CoefficientFunction\n\
+    mesh, intrule, ptrInput ... scalar input CF, points of x axis, points of y axis, cFieldIn='H'/'B', cFieldOut = 'H/B/mu/mudiff', mask=1, order=1/3")
+	.def(py::init(
+			[](shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
+                        const vector<shared_ptr<CoefficientFunction>> ptrInput_i, 
+                        vector<double> H_KL, 
+                        vector<double> B_KL,
+                        vector<double> ti,  
+                        const std::string cFieldIn, const std::string cFieldOut, const shared_ptr<CoefficientFunction>& ptrMask, unsigned int order, 
+                        bool usePreisBiro)
+	{
+		return new nonLinCF_harmonicBalance<T, DIM>(mesh, intrule, ptrInput_i, H_KL, B_KL, ti, cFieldIn, cFieldOut, ptrMask, order, usePreisBiro);
+	}), py::arg("mesh"), py::arg("IntegrationRule"), py::arg("InputFunctions"), py::arg("H_KL"), py::arg("B_KL"), py::arg("ti"), py::arg("field_in") = "H", py::arg("field_out") = "B",py::arg("mask")=nullptr, py::arg("order")=3, py::arg("usePreisBiro")=false)
+    .def("Update",  overload_cast_<>()(&nonLinCF_harmonicBalance<T, DIM>::Update), "Updates the output CF according to the current input CF")
+    .def("Update",  overload_cast_<const vector<shared_ptr<CoefficientFunction>> >()(&nonLinCF_harmonicBalance<T, DIM>::Update), "Updates the output CF according to the given input CF")
+    .def("getParts", &nonLinCF_harmonicBalance<T, DIM>::getParts)
+    .def("getTi", &nonLinCF_harmonicBalance<T, DIM>::getTi)
+    .def_property("omega0",  &nonLinCF_harmonicBalance<T, DIM>::getOmega, &nonLinCF_harmonicBalance<T, DIM>::setOmega)
+
+    ;
+}
+
+
+void ExportReclibNonLin(py::module &main_m)
 {
 
     // --------------------------------------------------------------------------- *
     // -----------   NonLin CF KL                                                 
     // --------------------------------------------------------------------------- *
+    py::module_ m = main_m.def_submodule("nonLin", "nonLinear CoefficientFunction");
 
     declare_nonLinCF<double, 1>(m);
     declare_nonLinCF<double, 2>(m);
     declare_nonLinCF<double, 3>(m);
+
+    declare_nonLinCF_harmonicBalance<double, 1>(m);
+    declare_nonLinCF_harmonicBalance<double, 2>(m);
+    declare_nonLinCF_harmonicBalance<double, 3>(m);
+
 }
 
 

C++/nonLin/nonLinCF.cc

@@ -43,7 +43,6 @@ nonLinCF<T, DIM>::nonLinCF(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrul
     // Update();
 }; 
 
-
 /* --------------------------------------------------------------------------- *
 *  -----------   Member Funcitons                                              *
 *  --------------------------------------------------------------------------- */
@@ -56,7 +55,6 @@ void nonLinCF<T, DIM>::findConcavePoints(){
 			while(oKL->GetData()[i].first == 0){
 				++i;
 			}
-			cout << oKL->GetData().size() << endl;
 			double max = oKL->GetData()[i].second/oKL->GetData()[i].first;
 			for(i = i+1; i!= oKL->GetData().size(); ++i){
 				if(oKL->GetData()[i].second/oKL->GetData()[i].first > max){
@@ -68,7 +66,7 @@ void nonLinCF<T, DIM>::findConcavePoints(){
 			concaveB = std::fabs(oKL->GetData()[ind_max].second);
 			concaveH = std::fabs(oKL->GetData()[ind_max].first);
 
-			cout << "concaveB" << concaveB << " " << concaveH << endl;
+			
 		}
 
 
@@ -408,6 +406,183 @@ unsigned int nonLinCF<T, DIM>::getOrder() const{
 }
 
 
+
+// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+template <typename T, unsigned int DIM>
+nonLinCF_harmonicBalance<T, DIM>::nonLinCF_harmonicBalance(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
+            const vector<shared_ptr<CoefficientFunction>> ptrInput_i, 
+            vector<double> H_KL, 
+            vector<double> B_KL, 
+            vector<double> ti, 
+            const std::string cFieldIn, const std::string cFieldOut, 
+            const shared_ptr<CoefficientFunction>& ptrMask, unsigned int order, bool usePreisBiro):
+            nonLinCF<T, DIM>{mesh, intrule, nullptr, H_KL, B_KL, cFieldIn, cFieldOut, ptrMask, order, usePreisBiro}, 
+            ptrInput_i{ptrInput_i}, ti{ti}, omega_0{2.0*M_PI/(ti.back() - ti.front())}{
+                unsigned int i;
+
+                for(i = 0; i!= ti.size(); ++i){
+                    nonLin_parts.push_back(std::make_shared<nonLinCF<T, DIM>>(mesh, intrule, this->mips, nullptr, H_KL, B_KL, cFieldIn, cFieldOut, ptrMask, order, usePreisBiro ));
+                }
+                Update();
+            }
+
+template <typename T, unsigned int DIM>
+void nonLinCF_harmonicBalance<T, DIM>::Update(){
+    Update(ptrInput_i);
+}
+
+
+
+
+
+template <typename T, unsigned int DIM>
+void nonLinCF_harmonicBalance<T, DIM>::Update(const vector<shared_ptr<CoefficientFunction>> ptrInput_i){
+
+
+    if(ptrInput_i.size() == 0)
+        cout << __PRETTY_FUNCTION__ << ": invalid input pointer" << endl;
+
+
+ 
+    parFor(this->mips.size(), [&](int el) {
+        unsigned int ip;
+        std::complex<double> val_in_c;
+        Vector<Complex> fV_tmp(1);
+        
+        double val_in;
+        double out_old;
+        unsigned int i;
+        unsigned int j;
+
+        for(ip = 0; ip != this->mips[el].size(); ++ip){
+
+
+            // --------------------- update parts -----------------
+            for(i = 0; i!= ti.size(); ++i){
+
+                val_in_c = 0;
+                for(j = 0; j!= ptrInput_i.size(); ++j){
+                    ptrInput_i[j]->Evaluate(this->mips[el][ip], fV_tmp);
+
+                    val_in_c += fV_tmp(0) * std::exp( j * 1i * omega_0 * ti[i]);
+                }
+
+                val_in = std::abs(val_in_c.real());
+
+                if(!std::strcmp(this->sFieldIn.c_str(), "H")){       
+                    
+                    if(!std::strcmp(this->sFieldOut.c_str(), "B"))
+                        nonLin_parts[i]->Set(el, ip, this->oKL->EvaluateH(val_in));
+                    else if(!std::strcmp(this->sFieldOut.c_str(), "mu")){
+                        if (this->usePreisBiro && std::abs(val_in) <= this->concaveH){
+                            // double	b = h*pastMu[el][p];
+                            // double hnew = bh->GetHfromB(b);
+                            // return b / hnew;
+                            val_in = this->oKL->EvaluateB(val_in * nonLin_parts[i]->Get(el, ip).getxx());
+                        }   
+                        if(abs(val_in) < 1e-6)
+                            nonLin_parts[i]->Set(el, ip, this->oKL->getZeroMu());
+                        else{
+
+                            nonLin_parts[i]->Set(el, ip, this->oKL->EvaluateH(val_in)/val_in);
+                        }
+                    }
+                    else if(!std::strcmp(this->sFieldOut.c_str(), "mu_diff")){
+                        if (this->usePreisBiro && std::abs(val_in) <= this->concaveH){
+                            val_in = this->oKL->EvaluateB(val_in * nonLin_parts[i]->Get(el, ip).getxx());
+                        }   
+                        nonLin_parts[i]->Set(el, ip, this->oKL->DifferentiateH(val_in));
+                    }
+                    else{
+                        throw std::invalid_argument("invalid input parameter: valid parameters are (H, B, mu, mu_diff)");
+                    }
+                }
+
+                else {
+                    if(!std::strcmp(this->sFieldOut.c_str(), "H"))
+                        nonLin_parts[i]->Set(el, ip, this->oKL->EvaluateB(val_in));
+                    else if(!std::strcmp(this->sFieldOut.c_str(), "nu")){
+                        if (this->usePreisBiro && std::abs(val_in) > this->concaveB){
+                            val_in = this->oKL->EvaluateH(val_in * nonLin_parts[i]->Get(el, ip).getxx());
+                        }   
+                        if(abs(this->oKL->EvaluateB(val_in)) < 1e-6)
+                            nonLin_parts[i]->Set(el, ip, 1/(this->oKL->getZeroMu()));
+                        else
+                            nonLin_parts[i]->Set(el, ip, this->oKL->EvaluateB(val_in)/val_in);
+                    }
+
+                    else if(!std::strcmp(this->sFieldOut.c_str(), "nu_diff")){
+                        if (this->usePreisBiro && std::abs(val_in) > this->concaveB){
+                            val_in = this->oKL->EvaluateH(val_in * nonLin_parts[i]->Get(el, ip).getxx());
+                        }   
+                        nonLin_parts[i]->Set(el, ip, this->oKL->DifferentiateB(val_in));
+                    }
+                    else{
+                        throw std::invalid_argument("invalid input parameter: valid parameters are (H, B, nu, nu_diff)");
+                    }
+
+                }
+            
+            }
+
+            
+
+            // --------------------- update global -----------------
+            out_old = this->Get(el, ip).xx();
+            val_in = 0;
+            for(i = 0; i!= ti.size()-1; ++i){
+                val_in += nonLin_parts[i]->Get(el, ip).xx() * (ti[i+1] - ti[i]);
+            }    
+            this->Set(el, ip, val_in * 1.0/(ti.back() - ti.front()));
+
+
+            // --------------------- error calc -----------------
+            if (this->ptrError != nullptr){
+                // cout << "error calculation" << endl;
+                this->ptrError->Set(el, ip,  this->Get(el, ip).xx() - out_old);
+            }
+        }
+
+    });
+
+
+
+
+
+
+
+    this->maxError = 0;
+    double sumOut = 0;
+    this->meanError = 0;
+
+    
+    if (this->ptrError != nullptr){
+        int el, ip;
+        for(el = 0; el != this->mips.size(); ++el){
+            for(ip = 0; ip != this->mips[el].size(); ++ip){                
+                this->maxError = std::max(this->ptrError->Get(el, ip).norm()/this->Get(el, ip).norm(), this->maxError);
+                this->meanError += this->ptrError->Get(el, ip).norm();
+                sumOut += this->Get(el, ip).norm();
+            }
+        }
+        
+        if (sumOut > 1e-14) {
+            this->meanError = this->meanError / sumOut;
+        }
+        else {
+            this->meanError = 0;
+        }
+    }
+
+    this->meanError *= 100;
+    this->maxError *= 100;
+}
+
 template class nonLinCF<double, 1>;
 template class nonLinCF<double, 2>;
 template class nonLinCF<double, 3>;
+
+template class nonLinCF_harmonicBalance<double, 1>;
+template class nonLinCF_harmonicBalance<double, 2>;
+template class nonLinCF_harmonicBalance<double, 3>;
+

C++/nonLin/nonLinCF.h

@@ -13,6 +13,8 @@
 #include <cstring>
 
 
+
+
 #include "../basicClasses/ngData_class.h"
 #include "../basicClasses/magCurve.h"
 
@@ -22,6 +24,8 @@ using namespace ngfem;
 using namespace ngsolve;
 using namespace ngcomp;
 
+using namespace std::complex_literals;
+
 // #define _Pi 3.141592653589793
 
 template <unsigned int DIM>
@@ -46,6 +50,8 @@ class nonLinCF:public ngScalar_class<T, DIM>
                         const shared_ptr<CoefficientFunction>& ptrMask=nullptr, unsigned int order=1, bool usePreisBiro=false);
 
 
+
+
         ~nonLinCF() = default;
 
         using ngScalar_class<T, DIM>::mips;
@@ -55,8 +61,8 @@ class nonLinCF:public ngScalar_class<T, DIM>
 
         bool HasConverged(double aMeanError, double aMaxError, bool print_it=true);
 
-        void Update();
-        void Update(const shared_ptr<CoefficientFunction>& ptrInput);
+        virtual void Update();
+        virtual void Update(const shared_ptr<CoefficientFunction>& ptrInput);
         
         void UpdateEnergyDensity();
         const shared_ptr<ngScalar_class<T, DIM>> GetEnergyDensity();
@@ -102,7 +108,7 @@ class nonLinCF:public ngScalar_class<T, DIM>
 
       
 
-    private:
+    protected:
 
 
         shared_ptr<CoefficientFunction> ptrInput;
@@ -127,8 +133,52 @@ class nonLinCF:public ngScalar_class<T, DIM>
         double maxError;
 
 
+
+    private:
         void findConcavePoints();
 };
 
+template<typename T, unsigned int DIM>
+class nonLinCF_harmonicBalance : public nonLinCF<T, DIM>{
+
+    public:
+        nonLinCF_harmonicBalance(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
+                    const vector<shared_ptr<CoefficientFunction>> ptrInput_i, 
+                    vector<double> H_KL, 
+                    vector<double> B_KL, 
+                    vector<double> ti, 
+                    const std::string cFieldIn="H", const std::string cFieldOut="B", 
+                    const shared_ptr<CoefficientFunction>& ptrMask=nullptr, unsigned int order=1, bool usePreisBiro=false);
+    
+
+        virtual void Update();
+        virtual void Update(const vector<shared_ptr<CoefficientFunction>> ptrInput);
+
+        
+        double getOmega()const{
+            return omega_0;
+        }
+        
+        void setOmega(const double val){
+            omega_0 = val;
+        }
+
+        vector<double> getTi() const{
+            return  ti;
+        }
+
+        vector<std::shared_ptr<nonLinCF<T, DIM>>> getParts() const{
+            return nonLin_parts;
+        }
+
+
+    private:
+        vector<shared_ptr<CoefficientFunction>> ptrInput_i;
+        vector<double> ti;
+        double omega_0 = 0;
+        vector<std::shared_ptr<nonLinCF<T, DIM>>> nonLin_parts;
+        
+};
+
 #endif /* NONLINCF_H */