From 602d5cebae55a004ff274f41536005c04455e627 Mon Sep 17 00:00:00 2001
From: Seweryn Kokot <sewkokot@gmail.com>
Date: Mon, 29 Nov 2021 12:41:15 +0100
Subject: [PATCH] HarmonicSS: use a load period instead of a circular frequency

---
 .../integrator/HarmonicSteadyState.cpp        | 26 ++++++++++---------
 SRC/analysis/integrator/HarmonicSteadyState.h |  4 +--
 2 files changed, 16 insertions(+), 14 deletions(-)

diff --git a/SRC/analysis/integrator/HarmonicSteadyState.cpp b/SRC/analysis/integrator/HarmonicSteadyState.cpp
index df7e641164..fbee15e8cb 100644
--- a/SRC/analysis/integrator/HarmonicSteadyState.cpp
+++ b/SRC/analysis/integrator/HarmonicSteadyState.cpp
@@ -61,10 +61,10 @@ void* OPS_HarmonicSteadyState()
 	return 0;
     }
 
-	double prads = 0;
+	double period = 0;
 	numData = 1;
-	if(OPS_GetDoubleInput(&numData,&prads) < 0) {
-	    opserr<<"WARNING failed to read double prads\n";
+	if(OPS_GetDoubleInput(&numData,&period) < 0) {
+	    opserr<<"WARNING failed to read double period\n";
 	    return 0;
 	}
 
@@ -82,12 +82,12 @@ void* OPS_HarmonicSteadyState()
 	}
     }
 
-    return new HarmonicSteadyState(lambda,prads,numIter,mLambda[0],mLambda[1]);
+    return new HarmonicSteadyState(lambda,period,numIter,mLambda[0],mLambda[1]);
 }
 
-HarmonicSteadyState::HarmonicSteadyState(double dLambda, double prads, int numIncr, double min, double max,  int classtag)
+HarmonicSteadyState::HarmonicSteadyState(double dLambda, double period, int numIncr, double min, double max,  int classtag)
     : StaticIntegrator(classtag),
- deltaLambda(dLambda), loadCircFreq(prads),
+ deltaLambda(dLambda), loadPeriod(period),
  specNumIncrStep(numIncr), numIncrLastStep(numIncr),
 	  dLambdaMin(min), dLambdaMax(max), gradNumber(0), sensitivityFlag(0)
 {
@@ -175,7 +175,7 @@ HarmonicSteadyState::sendSelf(int cTag,
 {
   Vector data(6);
   data(0) = deltaLambda;
-  data(1) = loadCircFreq;
+  data(1) = loadPeriod;
   data(2) = specNumIncrStep;
   data(3) = numIncrLastStep;
   data(4) = dLambdaMin;
@@ -199,7 +199,7 @@ HarmonicSteadyState::recvSelf(int cTag,
       return -1;
   }
   deltaLambda = data(0);
-  loadCircFreq = data(1);
+  loadPeriod = data(1);
   specNumIncrStep = data(2);
   numIncrLastStep = data(3);
   dLambdaMin = data(4);
@@ -217,7 +217,7 @@ HarmonicSteadyState::Print(OPS_Stream &s, int flag)
 	double currentLambda = theModel->getCurrentDomainTime();
 	s << "\t HarmonicSteadyState - currentLambda: " << currentLambda;
 	s << "  deltaLambda: " << deltaLambda << endln;
-	s << "  loadCircFreq: " << loadCircFreq << " rad/s" << endln;
+	s << "  Load Period: " << loadPeriod << endln;
     } else
 	s << "\t HarmonicSteadyState - no associated AnalysisModel\n";
 
@@ -226,19 +226,21 @@ HarmonicSteadyState::Print(OPS_Stream &s, int flag)
 int
 HarmonicSteadyState::formEleTangent(FE_Element *theEle)
 {
+  static const double twoPi = 2*3.1415926535897932;
+  static double twoPiSquareOverPeriodSquare = twoPi*twoPi/(loadPeriod*loadPeriod);
   if (statusFlag == CURRENT_TANGENT) {
     theEle->zeroTangent();
     theEle->addKtToTang();
-	theEle->addMtoTang(-loadCircFreq*loadCircFreq);
+	theEle->addMtoTang(-twoPiSquareOverPeriodSquare);
   } else if (statusFlag == INITIAL_TANGENT) {
     theEle->zeroTangent();
     theEle->addKiToTang();
-	theEle->addMtoTang(-loadCircFreq*loadCircFreq);
+	theEle->addMtoTang(-twoPiSquareOverPeriodSquare);
   } else if (statusFlag == HALL_TANGENT)  {
     theEle->zeroTangent();
     theEle->addKtToTang(cFactor);
     theEle->addKiToTang(iFactor);
-	theEle->addMtoTang(-loadCircFreq*loadCircFreq);
+	theEle->addMtoTang(-twoPiSquareOverPeriodSquare);
   }
 
   return 0;
diff --git a/SRC/analysis/integrator/HarmonicSteadyState.h b/SRC/analysis/integrator/HarmonicSteadyState.h
index 924417d516..1e6042e831 100644
--- a/SRC/analysis/integrator/HarmonicSteadyState.h
+++ b/SRC/analysis/integrator/HarmonicSteadyState.h
@@ -45,7 +45,7 @@ class ReliabilityDomain;
 class HarmonicSteadyState : public StaticIntegrator
 {
   public:
-    HarmonicSteadyState(double deltaLambda, double loadCircFreq, int numIncr,
+    HarmonicSteadyState(double deltaLambda, double loadPeriod, int numIncr,
 			   double minLambda, double maxlambda,
 			   int classtag=INTEGRATOR_TAGS_HarmonicSteadyState);
 
@@ -80,7 +80,7 @@ class HarmonicSteadyState : public StaticIntegrator
 
   private:
     double deltaLambda;  // dlambda at step (i-1)
-    double loadCircFreq; // excitation circular frequency (p) in rad/s
+    double loadPeriod; // load period in seconds (p = 2pi/T)
 
     double specNumIncrStep, numIncrLastStep; // Jd & J(i-1)
     double dLambdaMin, dLambdaMax; // min & max values for dlambda at step (i)