nonLinCF.h

#ifndef NONLINCF_H
#define NONLINCF_H


#include <fem.hpp>
#include <solve.hpp>
#include <ngstd.hpp>
#include <vector>
#include <cmath>
#include <memory>
#include <iostream>
#include <cmath>
#include <cstring>




#include "../basicClasses/ngData_class.h"
#include "../basicClasses/magCurve.h"


using namespace ngstd;
using namespace ngfem;
using namespace ngsolve;
using namespace ngcomp;

using namespace std::complex_literals;

// #define _Pi 3.141592653589793

template <unsigned int DIM>
std::pair<unsigned int, unsigned int> findNearestNeighbour(const vector<vector<MappedIntegrationPoint<DIM, DIM>> > & mips, const BaseMappedIntegrationPoint& mip);

template <typename T, unsigned int DIM>
class nonLinCF:public ngScalar_class<T, DIM>
{ 
            
    public:
        nonLinCF(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, const vector<vector<MappedIntegrationPoint<DIM, DIM>> > & mips,
                        const shared_ptr<CoefficientFunction>& ptrInput, 
                        vector<double> H_KL, 
                        vector<double> B_KL, 
                        const std::string cFieldIn="H", const std::string cFieldOut="B", 
                        const shared_ptr<CoefficientFunction>& ptrMask=nullptr, unsigned int order=1, bool usePreisBiro=false);
        nonLinCF(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
                        const shared_ptr<CoefficientFunction>& ptrInput, 
                        vector<double> H_KL, 
                        vector<double> B_KL, 
                        const std::string cFieldIn="H", const std::string cFieldOut="B", 
                        const shared_ptr<CoefficientFunction>& ptrMask=nullptr, unsigned int order=1, bool usePreisBiro=false);




        ~nonLinCF() = default;

        using ngScalar_class<T, DIM>::mips;
        using ngScalar_class<T, DIM>::data;



        bool HasConverged(double aMeanError, double aMaxError, bool print_it=true);

        virtual void Update();
        virtual void Update(const shared_ptr<CoefficientFunction>& ptrInput);
        
        void UpdateEnergyDensity();
        const shared_ptr<ngScalar_class<T, DIM>> GetEnergyDensity();
        const shared_ptr<ngScalar_class<T, DIM>> GetIterationError();

        

        nonLinCF<T, DIM> Integrate();
        nonLinCF<T, DIM> Differntiate();

        unsigned int getOrder() const;
        KL GetKL();

        shared_ptr<CoefficientFunction> GetInputCF() const;
        void SetInputCF(shared_ptr<CoefficientFunction>) ;


        NGS_DLL_HEADER AutoDiff<1,double> operator() (AutoDiff<1,double> x) const;

        virtual void Evaluate (const BaseMappedIntegrationRule & ir, 
                        BareSliceMatrix<AutoDiff<1,double>> values) const{


            // cout << ir.Size() << " " << 2*values.Dist() << " " << &values(0).Value() << endl;
            // SliceMatrix<double> hvalues(ir.Size(), 1, 2*values.Dist(), &values(0).Value());
            // Evaluate (ir, hvalues);
            for (size_t i = 0; i < ir.Size(); i++){
                for (size_t j = this->Dimension(); j-- > 0; ){
                    values(i,j) = this->operator()(values(i,j));
                }
            }
        
        }


        void setConcavePoints(double Hconcave, double Bconcave){
			concaveH = Hconcave;
			concaveB = Bconcave;
		}

        CoefficientFunction* operator()(std::shared_ptr<CoefficientFunction> ptrIn);
        T operator()(T );

      

    protected:


        shared_ptr<CoefficientFunction> ptrInput;
        
        const std::string sFieldIn;        
        const std::string sFieldOut;
        bool usePreisBiro;


        double concaveH;
        double concaveB;

        // data
        shared_ptr<ngScalar_class<double, DIM>> ptrError;
        shared_ptr<ngScalar_class<T, DIM>> ptrEnergyDens;

        
        std::shared_ptr<KL>  oKL;


        double meanError;
        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, 
                    const vector<double> omega_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);
    
        nonLinCF_harmonicBalance(const nonLinCF_harmonicBalance&);

        
        nonLinCF_harmonicBalance copy(){
            return nonLinCF_harmonicBalance(*this);
        }


        virtual void Update();
        virtual void Update(const vector<shared_ptr<CoefficientFunction>> ptrInput);

        
        vector<double> getOmegas()const{
            return omega_i;
        }
        
        void setOmegas(const vector<double> val){
            omega_i = val;
        }

        vector<double> getTi() const{
            return  ti;
        }

        vector<std::shared_ptr<nonLinCF<T, DIM>>> getParts() const{
            return nonLin_parts;
        }

        std::shared_ptr<nonLinCF<T, DIM>> operator [] (const unsigned int o) const{
            return nonLin_parts[o];
        }


    private:
        vector<shared_ptr<CoefficientFunction>> ptrInput_i;
        vector<double> omega_i;
        vector<double> ti;
        vector<std::shared_ptr<nonLinCF<T, DIM>>> nonLin_parts;
        
};

#endif /* NONLINCF_H */