Merge pull request #360 from yulee/ExpressNewton

ExpressNewton algorithm

SRC/analysis/algorithm/equiSolnAlgo/ExpressNewton.cpp

@@ -0,0 +1,163 @@
+/* ****************************************************************** **
+**    OpenSees - Open System for Earthquake Engineering Simulation    **
+**          Pacific Earthquake Engineering Research Center            **
+**                                                                    **
+**                                                                    **
+** (C) Copyright 1999, The Regents of the University of California    **
+** All Rights Reserved.                                               **
+**                                                                    **
+** Commercial use of this program without express permission of the   **
+** University of California, Berkeley, is strictly prohibited.  See   **
+** file 'COPYRIGHT'  in main directory for information on usage and   **
+** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
+**                                                                    **
+** Developed by:                                                      **
+**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
+**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
+**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
+**                                                                    **
+** ****************************************************************** */
+
+// $Revision: 1.0 $
+// $Date: 2019-01-28 17:53:01 $
+// $Source: /usr/local/cvs/OpenSees/SRC/analysis/algorithm/equiSolnAlgo/ExpressNewton.cpp,v $
+
+// Written: Yuli Huang (yulee@berkeley.edu)
+// Created: 02/2020
+// Revision: A
+//
+// Description: This file contains the implementation for the ExpressNewton class.
+// ExpressNewton is a class which performs a ExpressNewton solution algorithm
+// to solve the equations.
+//
+// Reference:
+// Junjie Xu, Yuli Huang, Zhe Qu,
+// An efficient and unconditionally stable numerical algorithm for nonlinear structural dynamics
+// International Journal for Numerical Methods in Engineering,
+// https://doi.org/10.1002/nme.6456.
+// (https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6456)
+
+// What: "@(#)ExpressNewton.h, revA"
+
+#include <ExpressNewton.h>
+#include <AnalysisModel.h>
+#include <StaticAnalysis.h>
+#include <StaticIntegrator.h>
+#include <LinearSOE.h>
+#include <Vector.h>
+#include <Channel.h>
+#include <FEM_ObjectBroker.h>
+#include <ConvergenceTest.h>
+#include <ID.h>
+
+#include <Timer.h>
+#include <elementAPI.h>
+#include <string>
+
+// Constructor
+ExpressNewton::ExpressNewton(int ni, double km, int tg, int fo)
+  :EquiSolnAlgo(EquiALGORITHM_TAGS_ExpressNewton), nIter(ni), factorOnce(fo) 
+{
+  if (tg == INITIAL_TANGENT) {
+    kMultiplier1 = km;
+    kMultiplier2 = 0.0;
+  }
+  else {
+    kMultiplier1 = 0.0;
+    kMultiplier2 = km;
+  }
+}
+
+// Destructor
+ExpressNewton::~ExpressNewton()
+{
+
+}
+
+int 
+ExpressNewton::solveCurrentStep(void)
+{
+    // set up some pointers and check they are valid
+    // NOTE this could be taken away if we set Ptrs as protecetd in superclass
+
+    AnalysisModel *theAnalysisModel = this->getAnalysisModelPtr(); 
+    LinearSOE  *theSOE = this->getLinearSOEptr();
+    IncrementalIntegrator *theIntegrator = this->getIncrementalIntegratorPtr();
+
+    if ((theAnalysisModel == 0) || (theIntegrator ==0 ) || (theSOE == 0)){
+	opserr << "WARNING ExpressNewton::solveCurrentStep() -";
+	opserr << "setLinks() has not been called.\n";
+	return -5;
+    }
+
+	if (factorOnce != 2) {
+		if (theIntegrator->formTangent(HALL_TANGENT, kMultiplier1, kMultiplier2) < 0) {
+		  opserr << "WARNING ExpressNewton::solveCurrentStep() -";
+		  opserr << "the Integrator failed in formTangent()\n";
+		  return -1;
+		}
+		if (factorOnce == 1)
+			factorOnce = 2;
+    }
+
+    for (int iter = 0; iter <nIter; ++iter)
+    {
+    if (theIntegrator->formUnbalance() < 0) {
+	opserr << "WARNING ExpressNewton::solveCurrentStep() -";
+	opserr << "the Integrator failed in formUnbalance()\n";	
+	return -2;
+    }
+
+    if (theSOE->solve() < 0) {
+	opserr << "WARNING ExpressNewton::solveCurrentStep() -";
+	opserr << "the LinearSOE failed in solve()\n";	
+	return -3;
+    }
+
+    if (theIntegrator->update(theSOE->getX()) < 0) {
+	opserr << "WARNING ExpressNewton::solveCurrentStep() -";
+	opserr << "the Integrator failed in update()\n";	
+	return -4;
+    }
+    }
+
+    return 0;
+}
+
+int
+ExpressNewton::setConvergenceTest(ConvergenceTest *theNewTest)
+{
+  return 0;
+}
+
+int
+ExpressNewton::sendSelf(int cTag, Channel &theChannel)
+{
+  static Vector data(4);
+  data(0) = nIter;
+  data(1) = kMultiplier1;
+  data(1) = kMultiplier2;
+  data(2) = factorOnce;
+  return theChannel.sendVector(this->getDbTag(), cTag, data);
+
+
+}
+
+int
+ExpressNewton::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
+{
+  static Vector data(4);
+  theChannel.recvVector(this->getDbTag(), cTag, data);
+  nIter = int(data(0));
+  kMultiplier1 = data(1);
+  kMultiplier2 = data(2);
+  factorOnce = int(data(3));
+  return 0;
+}
+
+
+void
+ExpressNewton::Print(OPS_Stream &s, int flag)
+{
+    s << "\t ExpressNewton algorithm";
+}

SRC/analysis/algorithm/equiSolnAlgo/ExpressNewton.h

@@ -0,0 +1,72 @@
+/* ****************************************************************** **
+**    OpenSees - Open System for Earthquake Engineering Simulation    **
+**          Pacific Earthquake Engineering Research Center            **
+**                                                                    **
+**                                                                    **
+** (C) Copyright 1999, The Regents of the University of California    **
+** All Rights Reserved.                                               **
+**                                                                    **
+** Commercial use of this program without express permission of the   **
+** University of California, Berkeley, is strictly prohibited.  See   **
+** file 'COPYRIGHT'  in main directory for information on usage and   **
+** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
+**                                                                    **
+** Developed by:                                                      **
+**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
+**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
+**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
+**                                                                    **
+** ****************************************************************** */
+
+// $Revision: 1.0 $
+// $Date: 2019-01-28 17:53:01 $
+// $Source: /usr/local/cvs/OpenSees/SRC/analysis/algorithm/equiSolnAlgo/ExpressNewton.h,v $
+
+// Written: Yuli Huang (yulee@berkeley.edu)
+// Created: 02/2020
+// Revision: A
+//
+// Description: This file contains the definition for the ExpressNewton class.
+// ExpressNewton is a class which performs a ExpressNewton solution algorithm
+// to solve the equations.
+//
+// Reference:
+// Junjie Xu, Yuli Huang, Zhe Qu,
+// An efficient and unconditionally stable numerical algorithm for nonlinear structural dynamics
+// International Journal for Numerical Methods in Engineering,
+// https://doi.org/10.1002/nme.6456.
+// (https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6456)
+
+// What: "@(#)ExpressNewton.h, revA"
+
+#ifndef ExpressNewton_h
+#define ExpressNewton_h
+
+#include <EquiSolnAlgo.h>
+
+class ExpressNewton: public EquiSolnAlgo
+{
+  public:
+    ExpressNewton(int nIter, double kMultiplier, int tagent, int factorOnce);
+    ~ExpressNewton();
+
+    int solveCurrentStep(void);
+    int setConvergenceTest(ConvergenceTest *theNewTest);
+
+    virtual int sendSelf(int commitTag, Channel &theChannel);
+    virtual int recvSelf(int commitTag, Channel &theChannel, 
+			 FEM_ObjectBroker &theBroker);
+
+    void Print(OPS_Stream &s, int flag =0);    
+
+  protected:
+
+  private:
+    int factorOnce;
+    int nIter;
+    double kMultiplier1, kMultiplier2;
+};
+
+#endif
+
+

SRC/classTags.h

@@ -888,6 +888,7 @@
 #define EquiALGORITHM_TAGS_InitialNewton          13
 #define EquiALGORITHM_TAGS_ElasticAlgorithm 14
 #define EquiALGORITHM_TAGS_NewtonHallM 15
+#define EquiALGORITHM_TAGS_ExpressNewton 16
 
 #define ACCELERATOR_TAGS_Krylov		1
 #define ACCELERATOR_TAGS_Secant		2

SRC/tcl/commands.cpp

@@ -143,6 +143,7 @@ OPS_Stream *opserrPtr = &sserr;
 #include <KrylovNewton.h>
 #include <PeriodicNewton.h>
 #include <AcceleratedNewton.h>
+#include <ExpressNewton.h>
 
 // accelerators
 #include <RaphsonAccelerator.h>
@@ -4022,6 +4023,27 @@ specifyAlgorithm(ClientData clientData, Tcl_Interp *interp, int argc,
       theNewAlgo = new NewtonLineSearch(*theTest, theLineSearch); 
   }
 
+  else if (strcmp(argv[1],"ExpressNewton") == 0) {
+    int nIter = 2, factorOnce = 0, formTangent = CURRENT_TANGENT;
+    double kMultiplier = 1.0;
+	  if (argc >= 3 && Tcl_GetInt(interp, argv[2], &nIter) != TCL_OK)
+      return TCL_ERROR;
+	  if (argc >= 4 && Tcl_GetDouble(interp, argv[3], &kMultiplier) != TCL_OK)
+      return TCL_ERROR;
+    int count = 4;
+    while (argc > count) {
+      if ((strcmp(argv[count],"-initialTangent") == 0) || (strcmp(argv[count],"-InitialTangent") == 0)) {
+        formTangent = INITIAL_TANGENT;
+      } else if ((strcmp(argv[count],"-currentTangent") == 0) || (strcmp(argv[count],"-CurrentTangent") ==0 )) {
+        formTangent = CURRENT_TANGENT;
+      } else if ((strcmp(argv[count],"-factorOnce") == 0) || (strcmp(argv[count],"-FactorOnce") ==0 )) {
+        factorOnce = 1;
+      }
+      count++;
+    }
+    theNewAlgo = new ExpressNewton(nIter,kMultiplier,formTangent,factorOnce);
+  }
+
   else {
     opserr << "WARNING No EquiSolnAlgo type " << argv[1] << " exists\n";
       return TCL_ERROR;