diff --git a/C++/multiscale/MSVPM.cc b/C++/multiscale/MSVPM.cc
index 573a4cc3d4a9c00d9da680352e22c35c6c07c14a..81330eead613cd3c056e774accc0adbeb90d70a4 100644
--- a/C++/multiscale/MSVPM.cc
+++ b/C++/multiscale/MSVPM.cc
@@ -8,14 +8,13 @@ namespace multiscale{
/* --------------------------------------------------------------------------- *
* ----------- MSHysteresis *
* --------------------------------------------------------------------------- */
-
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::coefHysteresisMS(
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::coefHysteresisMS(
const shared_ptr<MSWrapperVec<double, DIM_MESH>>& HorB,
const shared_ptr<EverettMatrix>& pE,
const shared_ptr<distributions>& ptrDist,
- char field, const bool differential ):
- HorB(HorB), pE{pE}, differential{differential}{
+ char field, const std::string fieldOut ):
+ HorB(HorB), fieldOut{fieldOut}{
unsigned int i, num;
if (field == 'B') {
@@ -65,7 +64,7 @@ coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::coefHysteresisMS(
- coef_local = std::make_shared<ngData_class_extended<mat_t<double, DIM_PREISACH, DIM_MS_COEF>, DIM_MESH>>(HorB->getMesh(), HorB->getIntrule(), HorB->getMips(), HorB->getLocalPoints(), HorB->getOrientation());
+ coef_local = std::make_shared<ngData_class_extended<mat_t<double, DIM_MS_COEF_R, DIM_MS_COEF_C>, DIM_MESH>>(HorB->getMesh(), HorB->getIntrule(), HorB->getMips(), HorB->getLocalPoints(), HorB->getOrientation());
ptrGlobalBCF = std::make_shared<ngData_class_extended<vec_t<double, DIM_PREISACH>, DIM_MESH>>(HorB->getMesh(), HorB->getIntrule(), HorB->getMips(), HorB->getLocalPoints(), HorB->getOrientation());
@@ -73,8 +72,28 @@ coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::coefHysteresisMS(
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-shared_ptr<CoefficientFunction> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Average(string name) {
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::coefHysteresisMS(
+ const shared_ptr<MSWrapperVec<double, DIM_MESH>>& HorB,
+ vector<vector<vector<shared_ptr<preisachVector<DIM_PREISACH>>>>> preisachVectorPoints,
+ char field, const std::string fieldOut ):
+ HorB(HorB), preisachVectorPoints{preisachVectorPoints}, fieldOut{fieldOut}{
+
+ unsigned int i, num;
+ if (field == 'B') {
+ inverse = true;
+ }
+ else {
+ inverse = false;
+ }
+
+ coef_local = std::make_shared<ngData_class_extended<mat_t<double, DIM_MS_COEF_R, DIM_MS_COEF_C>, DIM_MESH>>(HorB->getMesh(), HorB->getIntrule(), HorB->getMips(), HorB->getLocalPoints(), HorB->getOrientation());
+ ptrGlobalBCF = std::make_shared<ngData_class_extended<vec_t<double, DIM_PREISACH>, DIM_MESH>>(HorB->getMesh(), HorB->getIntrule(), HorB->getMips(), HorB->getLocalPoints(), HorB->getOrientation());
+
+}
+
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+shared_ptr<CoefficientFunction> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::Average(string name) {
if(parts.count(name) != 0)
return parts[name];
else{
@@ -87,21 +106,21 @@ shared_ptr<CoefficientFunction> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-map<string, shared_ptr<coefMSPart<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>>>& coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetParts() { return parts; }
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+map<string, shared_ptr<coefMSPart<DIM_MESH, DIM_MS_COEF_R, DIM_MS_COEF_C>>>& coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetParts() { return parts; }
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-shared_ptr<CoefficientFunction> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetParts(string name) { return Average(name); }
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+shared_ptr<CoefficientFunction> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetParts(string name) { return Average(name); }
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::RegisterPart(string name, shared_ptr<PhiFunction> phi1, shared_ptr<PhiFunction> phi2) {
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::RegisterPart(string name, shared_ptr<PhiFunction> phi1, shared_ptr<PhiFunction> phi2) {
// create part and update it
if(parts.count(name) != 0) {
cout << "part name " << name << " is registered already! Use a different name" <<endl;
return;
}
- parts[name] = make_shared<coefMSPart<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>>(coef_local, HorB->getPhi(phi1), HorB->getPhi(phi2), HorB->getWeights());
+ parts[name] = make_shared<coefMSPart<DIM_MESH, DIM_MS_COEF_R, DIM_MS_COEF_C>>(coef_local, HorB->getPhi(phi1), HorB->getPhi(phi2), HorB->getWeights());
HorB->Update();
@@ -109,8 +128,8 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::RegisterPart(string
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::SetUsePreviousForMuDiffOnly(const bool b){
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::SetUsePreviousForMuDiffOnly(const bool b){
unsigned int el, ip, l;
//parFor(mips.size(), [&](int el){
@@ -125,8 +144,8 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::SetUsePreviousForMuD
}
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetUsePreviousForMuDiffOnly()const{
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetUsePreviousForMuDiffOnly()const{
unsigned int el, ip, l;
//parFor(mips.size(), [&](int el){
for(el = 0; el != HorB->getMips().size(); ++el){
@@ -141,8 +160,8 @@ bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetUsePreviousForMuD
return false;
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::SetUsePerfectDemag(const bool b){
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::SetUsePerfectDemag(const bool b){
unsigned int el_copy = -1 , ip_copy = -1, l_copy = -1;
unsigned int el, ip, l;
@@ -166,8 +185,8 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::SetUsePerfectDemag(c
}
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetUsePerfectDemag()const{
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetUsePerfectDemag()const{
unsigned int el, ip, l;
//parFor(mips.size(), [&](int el){
for(el = 0; el != HorB->getMips().size(); ++el){
@@ -183,8 +202,8 @@ bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetUsePerfectDemag()
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-unsigned int coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::CountVectorPreisachPoints() const{
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+unsigned int coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::CountVectorPreisachPoints() const{
unsigned int el, ip;
unsigned int num = 0;
//parFor(HorB->getMips().size(), [&](int el){
@@ -198,8 +217,8 @@ unsigned int coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::CountVectorP
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Update(double scale) {
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::Update(double scale) {
//vec_tD tD_in;
@@ -230,16 +249,46 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Update(double scale)
// pilot
preisachVectorPoints[el][p][i]->pilotH(tD_in);
- if constexpr(DIM_PREISACH == DIM_MS_COEF){
- // mu
- if(differential)
- coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getPilots().fourth); // saved mu values from pilotH function
- else
- coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getPilots().third); // saved differential mu values from pilotH function
+
+
+ if constexpr(DIM_PREISACH > 1){
+ if constexpr(DIM_PREISACH == DIM_MS_COEF_C){
+ // mu
+ if(fieldOut=="mat_diff" || fieldOut=="mu_diff" || fieldOut=="nu_diff")
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getPilots().fourth); // saved differential mu values from pilotH function
+ else if(fieldOut=="mat" || fieldOut=="mu" || fieldOut=="nu")
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getPilots().third); // saved mu values from pilotH function
+ else
+ throw std::runtime_error("for MuhystMS only mu, nu, mat, mu_diff, nu_diff, mat_diff are allowed for field out");
+ }
+ else if constexpr(DIM_PREISACH == DIM_MS_COEF_R){
+ if(fieldOut=="B" || fieldOut=="Output")
+ // B
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getPilots().second); // saved B values from pilotH function
+ else if(fieldOut=="H" || fieldOut=="Input")
+ // H
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getPilots().first); // saved H values from pilotH function
+ else
+ throw std::runtime_error("for MuhystMS only B, H, Output, Input are allowed for field out");
+
+ }
+ else{
+ if(fieldOut=="loss")
+ // calcLossIncrement
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->calcLossIncrement(scale)); // saved B values from pilotH function
+ else if(fieldOut=="alternating_loss")
+ // calcLossIncrement
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->calcAlternatingLossIncrement(scale)); // saved H values from pilotH function
+ else if(fieldOut=="rotational_loss")
+ // calcLossIncrement
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->calcAlternatingLossIncrement(scale)); // saved H values from pilotH function
+ else
+ throw std::runtime_error("for MuhystMS only loss, alternating_loss, rotational_loss are allowed for field out");
+ }
}
else{
- // B
- coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getPilots().second); // saved mu values from pilotH function
+ throw std::runtime_error("1d case of coefHysteresisMS not implemented yet");
+
}
}
else{
@@ -280,8 +329,8 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Update(double scale)
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::UpdatePast(double scale) {
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::UpdatePast(double scale) {
HorB->Update();
//vec_tD tD_in;
@@ -304,21 +353,48 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::UpdatePast(double sc
}
}
- // pilot
+ // add
preisachVectorPoints[el][p][i]->addH(tD_in);
- if constexpr(DIM_PREISACH == DIM_MS_COEF){
- // mu
- if(differential)
- coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getMatDiff()); // saved mu values from pilotH function
- else
- coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getMat()); // saved differential mu values from pilotH function
+ if constexpr(DIM_PREISACH > 1){
+ if constexpr(DIM_PREISACH == DIM_MS_COEF_C){
+ // mu
+ if(fieldOut=="mat_diff" || fieldOut=="mu_diff" || fieldOut=="nu_diff")
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getMatDiff()); // saved differential mu values from pilotH function
+ else if(fieldOut=="mat" || fieldOut=="mu" || fieldOut=="nu")
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getMat()); // saved mu values from pilotH function
+ else
+ throw std::runtime_error("for MuhystMS only mu, nu, mat, mu_diff, nu_diff, mat_diff are allowed for field out");
+ }
+ else if constexpr(DIM_PREISACH == DIM_MS_COEF_R){
+ if(fieldOut=="B" || fieldOut=="Output")
+ // B
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getOutputVal()); // saved B values from pilotH function
+ else if(fieldOut=="H" || fieldOut=="Input")
+ // H
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getInputVal()); // saved H values from pilotH function
+ else
+ throw std::runtime_error("for MuhystMS only B, H, Output, Input are allowed for field out");
+
+ }
+ else{
+ if(fieldOut=="loss")
+ // calcLossIncrement
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->calcLossIncrement(scale)); // saved B values from pilotH function
+ else if(fieldOut=="alternating_loss")
+ // calcLossIncrement
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->calcAlternatingLossIncrement(scale)); // saved H values from pilotH function
+ else if(fieldOut=="rotational_loss")
+ // calcLossIncrement
+ coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->calcAlternatingLossIncrement(scale)); // saved H values from pilotH function
+ else
+ throw std::runtime_error("for MuhystMS only loss, alternating_loss, rotational_loss are allowed for field out");
+ }
}
else{
- //B
- coef_local->Set(el, p, i, preisachVectorPoints[el][p][i]->getOutputVal());
- }
-
+ throw std::runtime_error("1d case of coefHysteresisMS not implemented yet");
+
+ }
}
else{
@@ -352,8 +428,28 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::UpdatePast(double sc
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-unsigned int coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetNumBytes()const{
+
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::UpdateLossIncrements(double dt) {
+
+ HorB->Update();
+ //vec_tD tD_in;
+ parFor(preisachVectorPoints.size(), [&](int el) {
+ vec_t<double, DIM_PREISACH> tD_in;
+ //for (int el = 0; el < HorB->getMips().size(); el++) {
+ for (int p = 0; p < preisachVectorPoints[el].size(); p++) {
+ for (int i = 0;i < preisachVectorPoints[el][p].size();i++) {
+ preisachVectorPoints[el][p][i]->calcLossIncrement(dt);
+ preisachVectorPoints[el][p][i]->calcAlternatingLossIncrement(dt);
+ preisachVectorPoints[el][p][i]->calcRotationalLossIncrement(dt);
+ }
+ }
+ });
+
+}
+
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+unsigned int coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetNumBytes()const{
unsigned int el, ip, l;
unsigned int ret = 0;
@@ -368,16 +464,16 @@ unsigned int coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetNumBytes(
}
}
- ret += pE->size() * sizeof(double);
+ // ret += pE->size() * sizeof(double);
if(ptrGlobalBCF != nullptr) ret += ptrGlobalBCF->GetNumBytes();
return ret;
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Test() {
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::Test() {
- cout << "differential " << differential << endl;
+ cout << "fieldOut " << fieldOut << endl;
cout << "inverse " << inverse << endl;
cout << "reg: " << parts.size() << endl;
@@ -391,19 +487,19 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Test() {
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::isIron(const BaseMappedIntegrationPoint& mip, const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::isIron(const BaseMappedIntegrationPoint& mip, const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
return isIron(mip.GetPoint()[2], zheight_SingleSheet, iNumLayers, dFillFactor);
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::isIron(const double z, const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+bool coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::isIron(const double z, const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
return getGlobalBCF(zheight_SingleSheet, iNumLayers, dFillFactor)->isIron(z);
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::setData(const shared_ptr<const coefHysteresisMS> o){
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::setData(const shared_ptr<const coefHysteresisMS> o){
// copy all data
coef_local->setData(o->GetCoef_local()->getData());
for( auto p : o->parts){
@@ -412,8 +508,8 @@ void coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::setData(const shared
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-std::shared_ptr<ngScalar_class<double, DIM_MESH>> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::getSheetStructurCF(const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+std::shared_ptr<ngScalar_class<double, DIM_MESH>> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::getSheetStructurCF(const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
ptrStructureCF = std::make_shared<ngScalar_class<double, DIM_MESH>>(HorB->getMesh(), HorB->getIntrule(), HorB->getMips());
@@ -439,8 +535,8 @@ std::shared_ptr<ngScalar_class<double, DIM_MESH>> coefHysteresisMS<DIM_MESH, DIM
return ptrStructureCF;
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-std::shared_ptr<ngData_class_extended<vec_t<double, DIM_PREISACH>, DIM_MESH>> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::getGlobalBCF(const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+std::shared_ptr<ngData_class_extended<vec_t<double, DIM_PREISACH>, DIM_MESH>> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::getGlobalBCF(const double zheight_SingleSheet, const unsigned int iNumLayers, const double dFillFactor){
if(ptrGlobalBCF==nullptr){
@@ -453,20 +549,26 @@ std::shared_ptr<ngData_class_extended<vec_t<double, DIM_PREISACH>, DIM_MESH>> co
return ptrGlobalBCF;
}
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
-shared_ptr<MSWrapperVec<double, DIM_MESH>> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetHWrapper() const{
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
+shared_ptr<MSWrapperVec<double, DIM_MESH>> coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetHWrapper() const{
return HorB;
}
-
-template class coefHysteresisMS<1, 1, 1>;
-template class coefHysteresisMS<2, 2, 2>;
-template class coefHysteresisMS<3, 3, 3>;
-template class coefHysteresisMS<3, 2, 2>;
-
-template class coefHysteresisMS<2, 2, 1>;
-template class coefHysteresisMS<3, 3, 1>;
-template class coefHysteresisMS<3, 2, 1>;
+// mu
+template class coefHysteresisMS<1, 1, 1, 1>;
+template class coefHysteresisMS<2, 2, 2, 2>;
+template class coefHysteresisMS<3, 3, 3, 3>;
+template class coefHysteresisMS<3, 2, 2, 2>;
+
+// B
+template class coefHysteresisMS<2, 2, 2, 1>;
+template class coefHysteresisMS<3, 3, 3, 1>;
+template class coefHysteresisMS<3, 2, 2, 1>;
+
+// loss
+template class coefHysteresisMS<2, 2, 1, 1>;
+template class coefHysteresisMS<3, 3, 1, 1>;
+template class coefHysteresisMS<3, 2, 1, 1>;
}
\ No newline at end of file
diff --git a/C++/multiscale/MSVPM.h b/C++/multiscale/MSVPM.h
index cc954849c52143b919268aa75495e5790c7dfbb1..7e9d81ea65fe2a26fc69dd9854cc53668ade19aa 100644
--- a/C++/multiscale/MSVPM.h
+++ b/C++/multiscale/MSVPM.h
@@ -26,18 +26,22 @@ template<unsigned int DIM_MESH, unsigned int DIM_PREISACH>
using BHysteresisMSPart = coefMSPart<DIM_MESH, DIM_PREISACH, 1>;
-
-
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
class coefHysteresisMS {
public:
coefHysteresisMS(const shared_ptr<MSWrapperVec<double, DIM_MESH>>& HorB,
const shared_ptr<EverettMatrix>& pE,
const shared_ptr<distributions>& ptrDist,
- char field='H', const bool differential=true) ;
+ char field_in='H', const std::string fieldOut="mu") ;
+
+ coefHysteresisMS(const shared_ptr<MSWrapperVec<double, DIM_MESH>>& HorB,
+ vector<vector<vector<shared_ptr<preisachVector<DIM_PREISACH>>>>> preisachVectorPoints,
+ char field_in='H', const std::string fieldOut="mu") ;
+
+
shared_ptr<CoefficientFunction> Average(string name) ;
[[nodiscard]] shared_ptr<CoefficientFunction> GetParts(string name) ;
- [[nodiscard]] map<string, shared_ptr<coefMSPart<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>>>& GetParts();
+ [[nodiscard]] map<string, shared_ptr<coefMSPart<DIM_MESH, DIM_MS_COEF_R, DIM_MS_COEF_C>>>& GetParts();
void RegisterPart(string name, shared_ptr<PhiFunction> phi1, shared_ptr<PhiFunction> phi2) ;
@@ -47,9 +51,11 @@ public:
void UpdatePast(double scale=1) ;
+ void UpdateLossIncrements(double dt=1) ;
+
void Test() ;
- shared_ptr<ngData_class_extended<mat_t<double, DIM_PREISACH, DIM_MS_COEF>, DIM_MESH>> GetCoef_local()const{
+ shared_ptr<ngData_class_extended<mat_t<double, DIM_MS_COEF_R, DIM_MS_COEF_C>, DIM_MESH>> GetCoef_local()const{
return coef_local;
};
@@ -84,18 +90,17 @@ private:
shared_ptr<MSWrapperVec<double, DIM_MESH>> HorB;
- shared_ptr<ngData_class_extended<mat_t<double, DIM_PREISACH, DIM_MS_COEF>, DIM_MESH>> coef_local;
+ shared_ptr<ngData_class_extended<mat_t<double, DIM_MS_COEF_R, DIM_MS_COEF_C>, DIM_MESH>> coef_local;
- shared_ptr<EverettMatrix> pE;
vector<vector<vector<shared_ptr<preisachVector<DIM_PREISACH>>>>> preisachVectorPoints;
- map<string, shared_ptr<coefMSPart<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>>> parts;
+ map<string, shared_ptr<coefMSPart<DIM_MESH, DIM_MS_COEF_R, DIM_MS_COEF_C>>> parts;
std::shared_ptr<ngScalar_class<double, DIM_MESH>> ptrStructureCF;
std::shared_ptr<ngData_class_extended<vec_t<double, DIM_PREISACH>, DIM_MESH>> ptrGlobalBCF;
- bool differential;
+ const std::string fieldOut;
bool inverse;
// double mean_error = 1;
@@ -104,10 +109,10 @@ private:
template<unsigned int DIM_MESH, unsigned int DIM_PREISACH>
-using MuHysteresisMS = coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_PREISACH>;
+using MuHysteresisMS = coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_PREISACH, DIM_PREISACH>;
template<unsigned int DIM_MESH, unsigned int DIM_PREISACH>
-using BHysteresisMS = coefHysteresisMS<DIM_MESH, DIM_PREISACH, 1>;
+using BHysteresisMS = coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_PREISACH, 1>;
diff --git a/C++/multiscale/libMultiscale.cc b/C++/multiscale/libMultiscale.cc
index ef9564b2912afdcce05a2e95742c1c9eaf3729f5..e6771e2970d56423801bbd86f2eb1ec2db88ea2c 100644
--- a/C++/multiscale/libMultiscale.cc
+++ b/C++/multiscale/libMultiscale.cc
@@ -136,18 +136,23 @@ void declare_coefMSPart(py::module &m, std::string typestr) {
-template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF>
+template<unsigned int DIM_MESH, unsigned int DIM_PREISACH, unsigned int DIM_MS_COEF_R, unsigned int DIM_MS_COEF_C>
void declare_coefHysteresisMS(py::module &m, std::string typestr) {
std::string pyclass_name;
pyclass_name = std::string("");
- if constexpr(DIM_MS_COEF == DIM_PREISACH){
- pyclass_name += "Mu";
- }
- else{
- pyclass_name += "B";
+ if constexpr(DIM_PREISACH > 1){
+ if constexpr(DIM_MS_COEF_C == DIM_MS_COEF_R && DIM_MS_COEF_C == DIM_PREISACH){
+ pyclass_name += "Mu";
+ }
+ else if constexpr(DIM_MS_COEF_R == DIM_PREISACH){
+ pyclass_name += "B";
+ }
+ else{
+ pyclass_name += "loss";
+ }
}
if constexpr(DIM_MESH == DIM_PREISACH){
@@ -157,7 +162,7 @@ void declare_coefHysteresisMS(py::module &m, std::string typestr) {
pyclass_name += std::string("HysteresisMS") + std::to_string(DIM_MESH) + std::string("x") + std::to_string(DIM_PREISACH)+ typestr ;
}
- py::class_<coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>, shared_ptr<coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>>>
+ py::class_<coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>, shared_ptr<coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>>>
(m, pyclass_name.c_str(), "homogenised tensor-valued permeability for the multiscale method with the vector Preisach model.")
@@ -166,29 +171,39 @@ void declare_coefHysteresisMS(py::module &m, std::string typestr) {
[]( const shared_ptr<MSWrapperVec<double, DIM_MESH>>& HorB,
const shared_ptr<EverettMatrix>& pE,
const shared_ptr<distributions>& ptrDist,
- char field, const bool differential)
+ char field, const std::string fieldOut)
{
- return new coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>(HorB, pE, ptrDist, field, differential);
+ return new coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>(HorB, pE, ptrDist, field, fieldOut);
}), py::arg("HorB"), py::arg("EverettMatrix"),
- py::arg("PointDistribution"), py::arg("field") = "H", py::arg("differential")=false
+ py::arg("PointDistribution"), py::arg("field") = "H", py::arg("fieldOut")="mu"
+ )
+ .def(py::init(
+ []( const shared_ptr<MSWrapperVec<double, DIM_MESH>>& HorB,
+ vector<vector<vector<shared_ptr<preisachVector<DIM_PREISACH>>>>> preisachVectorPoints,
+ char field, const std::string fieldOut)
+ {
+ return new coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>(HorB, preisachVectorPoints, field, fieldOut);
+ }), py::arg("HorB"), py::arg("preisachVectorPoints"),
+ py::arg("field") = "H", py::arg("fieldOut")="mu"
)
- .def("Average", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Average, "function pointer to calculate the average a combination e.g. mu_av = CoefficientFunction([muAir, dmunonlin.Average('dmu_av')*kF + muAir*(1-kF)])")
- .def("GetParts", overload_cast_<>()(&coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetParts), "get all registered parts")
- .def("Test", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Test)
- .def("RegisterPart", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::RegisterPart, "register a combination e.g. .RegisterPart('mu_av','phi0','phi0')")
- .def("Update", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::Update, py::arg("scale")=1.0, "Update the object without changing the history of the VPMs")
- .def("UpdatePast", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::UpdatePast, py::arg("scale")=1.0, "Update the object by adding the current input CF to the history of the VPMs")
- .def("GetSheetStructurCF", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::getSheetStructurCF, "get a CF which represents the sheet structure")
- .def("GetGlobalBCF", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::getGlobalBCF, "get a visual output of the B field", py::arg("zheight_SingleSheet"), py::arg("numSheets"), py::arg("dFillFactor"))
- .def("SetUsePreviousForMuDiffOnly", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::SetUsePreviousForMuDiffOnly, "differential permeability/reluctivy should be calculated with the previous values only (backward Euler)")
- .def_property_readonly("GetUsePreviousForMuDiffOnly", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetUsePreviousForMuDiffOnly)
- .def("SetUsePerfectDemag", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::SetUsePerfectDemag, "perfect demagnetisation should be used for all SPMs")
- .def_property_readonly("GetUsePerfectDemag", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetUsePerfectDemag)
- .def_property_readonly("coefLocal", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetCoef_local)
- .def("CountVectorPreisachPoints", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::CountVectorPreisachPoints, "get the number of VPMs")
- .def("GetHWrapper", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetHWrapper, "get the Wrapper")
- .def("GetPreisachVectorPoints", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::GetPreisachVectorPoints, "get the VPMs")
- .def("SetData", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF>::setData, "copy the data of another coefHysteresisMS")
+ .def("Average", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::Average, "function pointer to calculate the average a combination e.g. mu_av = CoefficientFunction([muAir, dmunonlin.Average('dmu_av')*kF + muAir*(1-kF)])")
+ .def("GetParts", overload_cast_<>()(&coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetParts), "get all registered parts")
+ .def("Test", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::Test)
+ .def("RegisterPart", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::RegisterPart, "register a combination e.g. .RegisterPart('mu_av','phi0','phi0')")
+ .def("Update", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::Update, py::arg("scale")=1.0, "Update the object without changing the history of the VPMs")
+ .def("UpdatePast", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::UpdatePast, py::arg("scale")=1.0, "Update the object by adding the current input CF to the history of the VPMs")
+ .def("UpdateLossIncrements", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::UpdateLossIncrements, py::arg("dt")=1.0, "Update the object by adding the current input CF to the history of the VPMs")
+ .def("GetSheetStructurCF", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::getSheetStructurCF, "get a CF which represents the sheet structure")
+ .def("GetGlobalBCF", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::getGlobalBCF, "get a visual output of the B field", py::arg("zheight_SingleSheet"), py::arg("numSheets"), py::arg("dFillFactor"))
+ .def("SetUsePreviousForMuDiffOnly", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::SetUsePreviousForMuDiffOnly, "differential permeability/reluctivy should be calculated with the previous values only (backward Euler)")
+ .def_property_readonly("GetUsePreviousForMuDiffOnly", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetUsePreviousForMuDiffOnly)
+ .def("SetUsePerfectDemag", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::SetUsePerfectDemag, "perfect demagnetisation should be used for all SPMs")
+ .def_property_readonly("GetUsePerfectDemag", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetUsePerfectDemag)
+ .def_property_readonly("coefLocal", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetCoef_local)
+ .def("CountVectorPreisachPoints", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::CountVectorPreisachPoints, "get the number of VPMs")
+ .def("GetHWrapper", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetHWrapper, "get the Wrapper")
+ .def("GetPreisachVectorPoints", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::GetPreisachVectorPoints, "get the VPMs")
+ .def("SetData", &coefHysteresisMS<DIM_MESH, DIM_PREISACH, DIM_MS_COEF_R, DIM_MS_COEF_C>::setData, "copy the data of another coefHysteresisMS")
;
}
@@ -294,14 +309,19 @@ void ExportReclibMultiscale(py::module &main_m){
// * ---------------------------------------------------------------------------
// muMSHyst
- declare_coefHysteresisMS<3, 3, 3>(m, "");
- declare_coefHysteresisMS<3, 2, 2>(m, "");
- declare_coefHysteresisMS<2, 2, 2>(m, "");
+ declare_coefHysteresisMS<3, 3, 3, 3>(m, "");
+ declare_coefHysteresisMS<3, 2, 2, 2>(m, "");
+ declare_coefHysteresisMS<2, 2, 2, 2>(m, "");
// BMSHyst
- declare_coefHysteresisMS<3, 3, 1>(m, "");
- declare_coefHysteresisMS<3, 2, 1>(m, "");
- declare_coefHysteresisMS<2, 2, 1>(m, "");
+ declare_coefHysteresisMS<3, 3, 3, 1>(m, "");
+ declare_coefHysteresisMS<3, 2, 2, 1>(m, "");
+ declare_coefHysteresisMS<2, 2, 2, 1>(m, "");
+
+ // loss
+ declare_coefHysteresisMS<3, 3, 1, 1>(m, "");
+ declare_coefHysteresisMS<3, 2, 1, 1>(m, "");
+ declare_coefHysteresisMS<2, 2, 1, 1>(m, "");
declare_coefMSPart<2, 2, 2>(m, "");
diff --git a/python/MS_helper_functions.py b/python/MS_helper_functions.py
index 3efbb4e6b07895a904c0597372db54284160d403..fd646660e0d42e56d5f3dc25940f3262b1443349 100644
--- a/python/MS_helper_functions.py
+++ b/python/MS_helper_functions.py
@@ -1196,10 +1196,10 @@ class cl_MS():
if len(self.coupling_map[phi_i][phi_j].dims) == 2:
# mu_hyst
tmp = CF((tmp[0], tmp[1], 0, tmp[2], tmp[3], 0, 0, 0, np.mean([tmp[0], tmp[1], tmp[2], tmp[3]])), dims = (3,3))
- else:
+ elif self.coupling_map[phi_i][phi_j].dim > 1:
#B
tmp = CF((tmp[0], tmp[1], 0))
-
+
ret += InnerProduct(tmp * terms_u[i][0], terms_v[j][0])
else:
# normal