ngPreisachVector.cc

#include "ngPreisachVector.h"



/* --------------------------------------------------------------------------- *
*  -----------   Constructor                                                   *
*  --------------------------------------------------------------------------- */

template<unsigned int DIM>
ngPreisachVector<DIM>::ngPreisachVector(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
                    const shared_ptr<CoefficientFunction>& ptrInput, 
                    const shared_ptr<EverettMatrix>& pE, const shared_ptr<distributions>& ptrDist,
                    char cField, const shared_ptr<CoefficientFunction>& ptrMask):
                    CoefficientFunction(DIM), ptrInput{ptrInput}, ptrEverettMatrix{pE}, ptrMask{ptrMask}{
    
    if(cField == 'H')
        inverse=false;
    else
        inverse=true;

    
    unsigned int num=0;
    // create a container for every element
    mips.resize(mesh->GetNE());
    preisachVectorPoints.resize(mesh->GetNE());

    std::shared_ptr<preisachVector<DIM>> ptrPreisachP = std::make_shared<preisachVector<DIM>>(pE, ptrDist);

    ptrPreisachP->calcEnergyDensity();  // necessary to calc normEnergy and copy it for all
    cout << "create VectorPreisachPoint num:\t" << endl;

    //parFor(mips.size(), [&](int el){
    // for(el = 0; el != mips.size(); ++el){
    parFor(mips.size(), [&](int el) {
        unsigned ip;
        // create a container for every integration point
	    MappedIntegrationRule<DIM, DIM> mappedrule(intrule, mesh->GetTrafo(el, global_alloc), global_alloc);
	    mips[el].resize(intrule.GetNIP());
        preisachVectorPoints[el].resize(intrule.GetNIP());
	    for (ip = 0;ip < intrule.GetNIP();++ip) {
		    MappedIntegrationPoint<DIM, DIM> mip(intrule[ip], mesh->GetTrafo(el, global_alloc));
		    mips[el][ip] = { intrule[ip], mesh->GetTrafo(el, global_alloc) };
            
            //cout <<"(" <<  mips[el][ip].GetPoint()[0] << ", " << mips[el][ip].GetPoint()[1] << ", "  << mips[el][ip].GetPoint()[2] << ")" << endl;
            // cout << ptrMask->Evaluate(mips[el][ip]) << endl;
            if(ptrMask == nullptr || ptrMask->Evaluate(mips[el][ip])){
                //preisachVectorPoints[el][ip] = std::make_shared<preisachVector>(pE, ptrDist);
                preisachVectorPoints[el][ip] = std::make_shared<preisachVector<DIM>>(ptrPreisachP); // is faster to copy results of a single
                ++num;

                cout<< "\r" << num;

            }
            else{
                preisachVectorPoints[el].resize(0);
            }
	    }
    });
    cout << "\n";
    cout << "created " << CountVectorPreisachPoints() << " VectorPreisachPoints" << endl;

    Demagnetise();
} 


/* --------------------------------------------------------------------------- *
*  -----------   Member Funcitons                                              *
*  --------------------------------------------------------------------------- */
template<unsigned int DIM>
void ngPreisachVector<DIM>::RedefineCF(shared_ptr<CoefficientFunction> CF){
    ptrInput = CF;
}


template<unsigned int DIM>
double ngPreisachVector<DIM>::Evaluate(const BaseMappedIntegrationPoint& mip) const{
    return 0;
}


template<unsigned int DIM>
void  ngPreisachVector<DIM>::Evaluate(const BaseMappedIntegrationPoint& mip,  FlatVector<> result) const{
    if(!inverse)
        ptrFluxDens->Evaluate(mip, result);
    else
        ptrFieldStrength->Evaluate(mip, result);
    
}

template<unsigned int DIM>
void ngPreisachVector<DIM>::UpdatePast(){
    UpdatePast(ptrInput);
}
template<unsigned int DIM>
void ngPreisachVector<DIM>::UpdatePast(const shared_ptr<CoefficientFunction>& ptrIn){

    parFor(mips.size(), [&](unsigned int el) {
        double tmp[DIM] = {0};
        FlatVector fV_in(DIM, tmp);
        vec_tD tD_in;
    //for(int el = 0; el != mips.size(); ++el){

        unsigned int ip;
        
        for(ip = 0; ip < mips[el].size(); ++ip){
        //parFor(mips[el].size(), [&](int ip){
            

            //cout <<"value for UpdatePast " << value << "pointer to input:" << ptrInput << endl;
            // if a preisach model has been created - area == iron

            if(preisachVectorPoints[el].size() != 0){
                if(!inverse){
                    // evaluate with FlatVector
                    ptrIn->Evaluate(mips[el][ip], fV_in);

                    // cast it to vec_tD
                    tD_in.first = fV_in(0);
                    tD_in.second = fV_in(1);
                    tD_in.third = fV_in(2);

                    // add it as vec_tD
                    preisachVectorPoints[el][ip]->addH(tD_in);
                }
                else{

                    cout << "inverse model not implemented yet" << endl;
                    exit(-1);
    
                }

                // update Container

                if(updateB)ptrFluxDens->Set(el, ip, preisachVectorPoints[el][ip]->getB());
                if(updateH)ptrFieldStrength->Set(el, ip, preisachVectorPoints[el][ip]->getH()); 
                if(updateMu)ptrMat->Set(el, ip, preisachVectorPoints[el][ip]->getMu()); 
                if(updateMuDiff)ptrMatDiff->Set(el, ip, preisachVectorPoints[el][ip]->getMuDiff());     

                
            }
        }
    });


}
template<unsigned int DIM>
void ngPreisachVector<DIM>::Update(){  
    Update(ptrInput);
}
template<unsigned int DIM>
void ngPreisachVector<DIM>::Update(const shared_ptr<CoefficientFunction>& ptrIn){  

    // cout << "UPDATE of all Values" << endl;

    // unsigned int el;
    
    
    // update all preisach points
    parFor(mips.size(), [&](int el) {
        unsigned int ip;
        triple<vec_tD, mat_t<double, DIM, DIM>, mat_t<double, DIM, DIM>> val_pilot;
        double tmp[DIM] = {0};
        FlatVector fV_in(DIM, tmp);
        vec_tD tD_in;
        
    //for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip < mips[el].size(); ++ip){
            
            // if a preisach model has been created - area == iron
            if(preisachVectorPoints[el].size() != 0){

                // MuDiffold = MuDiff;
                if(!inverse){
                    // evaluate with FlatVector
                    ptrIn->Evaluate(mips[el][ip], fV_in);

                    // cast it to vec_tD
                    tD_in.first = fV_in(0);
                    tD_in.second = fV_in(1);
                    tD_in.third = fV_in(2);

                    
                    // pilot
                    val_pilot = preisachVectorPoints[el][ip]->pilotH_full(tD_in);


                    // cout << val_pilot.second.xx << endl;
                    
                }
                else{
                    cout << __FILE__ << ": inverse model not implemented yet" << endl;
                    exit(-1);
                }      
                // cout << "update containers" << endl;
                // save values for Error calculations
                if(updateB)ptrFluxDens->Set(el, ip, val_pilot.first);
                if(updateH)ptrFieldStrength->Set(el, ip, tD_in);
                if(updateMu)ptrMat->Set(el, ip, val_pilot.second);
                if(updateMuDiff)ptrMatDiff->Set(el, ip, val_pilot.third);
                // cout <<"Update\t" << el << ", "<< ip << ":\t" <<  ptrMat->Get(el, ip).xx << " \t\t" << val_pilot.second.xx << endl;


            }
        }
    });


}
template<unsigned int DIM>
void ngPreisachVector<DIM>::UpdateEnergyDensity(){  

    //cout << "UPDATE of all Values" << endl;

    unsigned int el;
    unsigned int ip;
    
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip < mips[el].size(); ++ip){
                
            // if a preisach model has been created - area == iron
            if(preisachVectorPoints[el][ip] != nullptr){
                // save values
                ptrEnergyDens->Set(el, ip, preisachVectorPoints[el][ip]->calcEnergyDensity());
            }
        }
    }
}
template<unsigned int DIM>
void ngPreisachVector<DIM>::Demagnetise(){
    int el_copy = -1 , ip_copy = -1;
    unsigned int el, ip;

    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                if(el_copy == -1 && ip_copy == -1){
                    preisachVectorPoints[el][ip]->demagnetise();

                    el_copy = el;
                    ip_copy = ip;                    
                }
                else{
                    *preisachVectorPoints[el][ip] = *preisachVectorPoints[el_copy][ip_copy];
                }
                if(ptrFieldStrength != nullptr)
                    ptrFieldStrength->Set(el, ip, preisachVectorPoints[el][ip]->getH());
                if(ptrFluxDens != nullptr)
                    ptrFluxDens->Set(el, ip, preisachVectorPoints[el][ip]->getB());
                if(ptrMat != nullptr)
                    ptrMat->Set(el, ip, preisachVectorPoints[el][ip]->getMu());
                if(ptrMatDiff != nullptr)
                    ptrMatDiff->Set(el, ip, preisachVectorPoints[el][ip]->getMuDiff());
                if(ptrEnergyDens != nullptr)
                    ptrEnergyDens->Set(el, ip, preisachVectorPoints[el][ip]->calcEnergyDensity());
            }
        }
    }

}

/* --------------------------------------------------------------------------- *
*  -----------   Properties                                                    *
*  --------------------------------------------------------------------------- */
template<unsigned int DIM>
shared_ptr<CoefficientFunction> ngPreisachVector<DIM>::GetCF(){
    return ptrInput;
}

/* --------------------------------------------------------------------------- *
*  -----------   Member Functions                                              *
*  --------------------------------------------------------------------------- */
// vector<array<double, 2>> ngPreisachVector<DIM>::GetKLAt(const BaseMappedIntegrationPoint& mip) const{
//     unsigned int el = mip.GetTransformation().GetElementNr();
//     if(preisachVectorPoints[el][0] != nullptr){
//         preisachVectorPoints[el][0]->setUpdateKL(true);
//         return preisachVectorPoints[el][0]->getKL(); 
//     }
//     else{
//         cout << "no preisach Model implemented at that point" << endl;
//         vector<array<double,2>> a;
//         return a;
//     }

// }
template<unsigned int DIM>
void ngPreisachVector<DIM>::SetUsePreviousForMuDiffOnly(const bool b){
    unsigned int el, ip;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                preisachVectorPoints[el][ip]->setUsePreviousForMuDiffOnly(b);
            }
        }
    }
            
}template<unsigned int DIM>
bool ngPreisachVector<DIM>::GetUsePreviousForMuDiffOnly()const{
    unsigned int el, ip;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                return preisachVectorPoints[el][ip]->getUsePreviousForMuDiffOnly();
            }
        }
    }
    return false;

}
template<unsigned int DIM>
void ngPreisachVector<DIM>::SetUsePerfectDemag(const bool b){

    unsigned int el_copy = -1 , ip_copy = -1;
    unsigned int el, ip;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                if(el_copy == -1 && ip_copy == -1){
                    preisachVectorPoints[el][ip]->setUsePerfectDemag(b);
                    el_copy = el;
                    ip_copy = ip;
                }
                else{
                    *preisachVectorPoints[el][ip] = *preisachVectorPoints[el_copy][ip_copy];
                }
            }
        }
    }
            
}template<unsigned int DIM>
bool ngPreisachVector<DIM>::GetUsePerfectDemag()const{
    unsigned int el, ip;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                return preisachVectorPoints[el][ip]->getUsePerfectDemag();
            }
        }
    }
    return false;

}template<unsigned int DIM>
void ngPreisachVector<DIM>::SetDifferential(const bool b){

    unsigned int el_copy = -1 , ip_copy = -1;
    unsigned int el, ip;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                if(el_copy == -1 && ip_copy == -1){
                    preisachVectorPoints[el][ip]->setDifferential(b);
                    el_copy = el;
                    ip_copy = ip;
                }
                else{
                    *preisachVectorPoints[el][ip] = *preisachVectorPoints[el_copy][ip_copy];
                }
            }
        }
    }
            
}template<unsigned int DIM>
bool ngPreisachVector<DIM>::GetDifferential()const{
    unsigned int el, ip;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                return preisachVectorPoints[el][ip]->getDifferential();
            }
        }
    }

    return false;
}

template<unsigned int DIM>
unsigned int ngPreisachVector<DIM>::GetDataSize()const{
    unsigned int el, ip;

    unsigned int ret = 0;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != mips.size(); ++el){
        for(ip = 0; ip != mips[el].size(); ++ip){
            if(preisachVectorPoints[el][ip] != nullptr){
                ret += preisachVectorPoints[el][ip]->getDataSize();
            }
        }
    }

    ret += ptrEverettMatrix->size() * sizeof(double);
    if(ptrFieldStrength != nullptr) ret += ptrFieldStrength->GetDataSize();
    if(ptrFluxDens != nullptr) ret += ptrFluxDens->GetDataSize();
    if(ptrMat != nullptr) ret += ptrMat->GetDataSize();
    if(ptrMatDiff != nullptr) ret += ptrMatDiff->GetDataSize();
    if(ptrEnergyDens != nullptr) ret += ptrMatDiff->GetDataSize();


    ret += sizeof(bool) * 5;        // auto updates

    return ret;
}
template<unsigned int DIM>
unsigned int ngPreisachVector<DIM>::CountVectorPreisachPoints() const{
    unsigned int el, ip;
    unsigned int num = 0;
    //parFor(mips.size(), [&](int el){
    for(el = 0; el != preisachVectorPoints.size(); ++el){
        for(ip = 0; ip != preisachVectorPoints[el].size(); ++ip){
            num += 1;
        }
    }

    return num;

}

template<unsigned int DIM>
const shared_ptr<ngScalar_class<DIM>> ngPreisachVector<DIM>::GetEnergyDensity(){
    // EnergyDensity
    if(ptrEnergyDens == nullptr){
        ptrEnergyDens = std::make_shared<ngScalar_class<DIM>>(mips, ptrMask);
    }

    return ptrEnergyDens;
}
template<unsigned int DIM>
const shared_ptr<ngVector_class<DIM>> ngPreisachVector<DIM>::GetB(){
    //----------- MEMORY DEMANDING ---------------
    // flux density
    if(ptrFluxDens == nullptr){
        ptrFluxDens = std::make_shared<ngVector_class<DIM>>(mips, ptrMask);
        updateB = true;
    }
    return ptrFluxDens;
}
template<unsigned int DIM>
const shared_ptr<ngVector_class<DIM>> ngPreisachVector<DIM>::GetH(){
    // field strength
    if(ptrFieldStrength == nullptr){
        ptrFieldStrength = std::make_shared<ngVector_class<DIM>>(mips, ptrMask);
        updateH = true;
    }
    return ptrFieldStrength;
}

template<unsigned int DIM>
const shared_ptr<ngMat_class<DIM>> ngPreisachVector<DIM>::GetMat(){


    // mu
    if(ptrMat == nullptr){
        ptrMat = std::make_shared<ngMat_class<DIM>>(mips, ptrMask);
        updateMu = true;
    }
    return ptrMat;
}

template<unsigned int DIM>
const shared_ptr<ngMat_class<DIM>> ngPreisachVector<DIM>::GetMatDiff(){


    // diff_mu
    if(ptrMatDiff == nullptr){
        ptrMatDiff = std::make_shared<ngMat_class<DIM>>(mips, ptrMask);
        updateMuDiff = true;
    }
    return ptrMatDiff;
}
// shared_ptr<CoefficientFunction> ngPreisachVector<DIM>::GetMu()const{
//     if(inverse){
//         ptrMat->SetInverse(true);
//         cout << "changed evaluation of Nu to Mu" << endl;
//     }
//     return ptrMat;
// }template<unsigned int DIM>
// shared_ptr<CoefficientFunction> ngPreisachVector<DIM>::GetNu()const{
//     if(!inverse){
//         ptrMat->SetInverse(true);
//         cout << "changed evaluation of Mu to Nu" << endl;
//     }
//     return ptrMat;
// }template<unsigned int DIM>
// shared_ptr<CoefficientFunction> ngPreisachVector<DIM>::GetMuDiff()const{
//     if(inverse){
//         ptrMatDiff->SetInverse(true);
//         cout << "changed evaluation of NuDiff to MuDiff" << endl;
//     }
//     return ptrMatDiff;
// }template<unsigned int DIM>
// shared_ptr<CoefficientFunction> ngPreisachVector<DIM>::GetNuDiff()const{
//     if(!inverse){
//         ptrMatDiff->SetInverse(true);
//         cout << "changed evaluation of MuDiff to NuDiff" << endl;
//     }
//     return ptrMatDiff;
// }

template<unsigned int DIM>
const shared_ptr<KL> ngPreisachVector<DIM>::GetInitialMagnetisationCurve()const{
    return  ptrEverettMatrix->getInitialMagnetisationCurve();
}

template<unsigned int DIM>
preisachVector<DIM> ngPreisachVector<DIM>::GetPreisachVector(const BaseMappedIntegrationPoint& mip) const{
    std::pair<unsigned int, unsigned int> id = findNearestNeighbour<DIM>(mips, mip);
    return preisachVectorPoints[id.first][id.second];
}


template class ngPreisachVector<1>;
template class ngPreisachVector<2>;
template class ngPreisachVector<3>;