diff --git a/C++/nonLin/libNonLin.cc b/C++/nonLin/libNonLin.cc
index 3c57814c598c54e253b8ba837c07d72b724b3dbe..5c5d63075c0d5df6926a12d2802cf49b23fc325f 100644
--- a/C++/nonLin/libNonLin.cc
+++ b/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);
+
}
diff --git a/C++/nonLin/nonLinCF.cc b/C++/nonLin/nonLinCF.cc
index b4969cc324a85a723ffa876f063c8ea44f6dad54..bf1cdee84f59821f31994a5b446a048534576bf2 100644
--- a/C++/nonLin/nonLinCF.cc
+++ b/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>;
+
diff --git a/C++/nonLin/nonLinCF.h b/C++/nonLin/nonLinCF.h
index a9c1e7d80e57c984d14eeede88f3ce0d8a755c19..46ffcd0af2b3614c0ca6f6453dec0400dcba2fb1 100644
--- a/C++/nonLin/nonLinCF.h
+++ b/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 */