Changed: Generalize ReactionsOnly class for internal force calculation.

Added: SIMbase::assembleForces() and IntegrandBase::setSecondaryInt().
Added: virtual method for calculation/print of interface force resultants.

src/ASM/IntegrandBase.h

@@ -84,6 +84,8 @@ public:
   virtual void initResultPoints(double, bool = false) {}
   //! \brief Initializes the global node number mapping for current patch.
   virtual void initNodeMap(const std::vector<int>&) {}
+  //! \brief Assigns a secondary integral to be computed (for reaction forces).
+  virtual void setSecondaryInt(GlobalIntegral* = nullptr) {}
   //! \brief Returns the system quantity to be integrated by \a *this.
   virtual GlobalIntegral& getGlobalInt(GlobalIntegral* gq) const;
 

src/ASM/ReactionsOnly.C

@@ -18,30 +18,34 @@
 #include "ProcessAdm.h"
 
 
-ReactionsOnly::ReactionsOnly (Vector& rf, const SAM* sam,
-                              const ProcessAdm& adm, Vector* sf)
+ReactionsOnly::ReactionsOnly (const SAM* sam, const ProcessAdm& adm,
+                              Vector* rf, Vector* sf)
   : mySam(sam), myAdm(adm), R(rf), S(sf)
 {
-  mySam->initForAssembly(b,&R);
+  mySam->initForAssembly(b,R);
 }
 
 
 void ReactionsOnly::initialize (bool)
 {
   b.init();
-  R.fill(0.0);
+  if (R)
+    R->fill(0.0);
 }
 
 
 bool ReactionsOnly::finalize (bool)
 {
- if (myAdm.dd.isPartitioned())
-   myAdm.allReduceAsSum(R);
+  if (R)
+  {
+    if (myAdm.dd.isPartitioned())
+      myAdm.allReduceAsSum(*R);
 
-  R *= -1.0;
+    (*R) *= -1.0;
 #if SP_DEBUG > 2
-  std::cout <<"\nReaction forces:"<< R;
+    std::cout <<"\nReaction forces:"<< *R;
 #endif
+  }
 
   if (!S)
     return true;
@@ -60,7 +64,7 @@ bool ReactionsOnly::finalize (bool)
 bool ReactionsOnly::assemble (const LocalIntegral* elmObj, int elmId)
 {
   const ElmMats* elMat = dynamic_cast<const ElmMats*>(elmObj);
-  if (elMat && mySam->assembleSystem(b,elMat->getRHSVector(),elmId,&R))
+  if (elMat && mySam->assembleSystem(b,elMat->getRHSVector(),elmId,R))
     return true;
 
   std::cerr <<" *** ReactionsOnly::assemble: Failure for element "<< elmId
@@ -69,23 +73,22 @@ bool ReactionsOnly::assemble (const LocalIntegral* elmObj, int elmId)
 }
 
 
-/*!
-  Returns \e true if the patch \a pidx have any dirichlet boundary conditions.
-*/
-
 bool ReactionsOnly::haveContributions (size_t pidx,
                                        const PropertyVec& pvec) const
 {
-  auto&& isConstrained = [this,pidx](const Property& p)
+  // Lambda function checking if a property have force contribution on a patch.
+  auto&& contributed = [this,pidx](const Property& p)
   {
     if (p.patch != pidx)
       return false;
     else if (S)
       return p.pcode == Property::OTHER;
+    else if (!R)
+      return false;
 
     return (p.pcode >= Property::DIRICHLET &&
             p.pcode <= Property::DIRICHLET_ANASOL);
   };
 
-  return std::find_if(pvec.begin(),pvec.end(),isConstrained) != pvec.end();
+  return std::find_if(pvec.begin(),pvec.end(),contributed) != pvec.end();
 }

src/ASM/ReactionsOnly.h

@@ -7,7 +7,7 @@
 //!
 //! \author Knut Morten Okstad / SINTEF
 //!
-//! \brief Global integral sub-class for reaction force integration.
+//! \brief Global integral class for reaction- and interface force calculation.
 //!
 //==============================================================================
 
@@ -19,7 +19,7 @@
 
 
 /*!
-  \brief Class for assembly of reaction forces only.
+  \brief Class for assembly of reaction- and interface forces.
 
   \details This class can be used for linear problems which normally consists
   of a single assembly loop. For non-linear problems, the reaction forces for
@@ -27,18 +27,21 @@
   the solution vector of the previous iteration is used. In linear problems we
   therefore need a separate assembly loop where only the reaction forces are
   calculated. This class is provided to facilitate such calculations.
+
+  The class can also be used to calculate interface force resultants,
+  for boundaries not associated with dirichlet conditions.
 */
 
 class ReactionsOnly : public GlobalIntegral
 {
 public:
   //! \brief The constructor initializes the data members.
-  //! \param[in] rf The reaction force vector to be assembled
   //! \param[in] sam Data for FE assembly management
   //! \param[in] adm Parallell processing administrator
+  //! \param[in] rf Reaction force vector to be assembled
   //! \param[in] sf Internal force vector to be assembled
-  ReactionsOnly(Vector& rf, const SAM* sam, const ProcessAdm& adm,
-                Vector* sf = nullptr);
+  ReactionsOnly(const SAM* sam, const ProcessAdm& adm,
+                Vector* rf = nullptr, Vector* sf = nullptr);
   //! \brief Empty destructor.
   virtual ~ReactionsOnly() {}
 
@@ -49,10 +52,10 @@ public:
 
   //! \brief Adds a LocalIntegral object into a corresponding global object.
   //! \param[in] elmObj The local integral object to add into \a *this.
-  //! \param[in] elmId Global number of the element associated with elmObj
+  //! \param[in] elmId Global number of the element associated with \a elmObj
   virtual bool assemble(const LocalIntegral* elmObj, int elmId);
 
-  //! \brief Returns \e false if no contributions from the specified patch.
+  //! \brief Returns \e true if the patch \a pidx have any force contributions.
   //! \param[in] pidx 1-based patch index
   //! \param[in] pvec Physical property mapping
   virtual bool haveContributions(size_t pidx,
@@ -62,9 +65,9 @@ private:
   const SAM*        mySam; //!< Data for FE assembly management
   const ProcessAdm& myAdm; //!< Parallel processing administrator
 
-  Vector&   R; //!< Nodal reaction forces
-  StdVector b; //!< Dummy right-hand-side vector
-  Vector*   S; //!< Internal force vector
+  StdVector b; //!< Internal right-hand-side vector used in the assembly process
+  Vector*   R; //!< Nodal reaction forces
+  Vector*   S; //!< Nodal internal forces
 };
 
 #endif

src/ASM/TimeDomain.h

@@ -31,6 +31,9 @@ struct TimeDomain
   //! \brief Default constructor.
   explicit TimeDomain(int i = 0, bool f = true) : it(i), first(f)
   { t = dt = dtn = CFL = 0.0; }
+  //! \brief Constructor for linear problems with fixed (non-zero) time.
+  explicit TimeDomain(double t0) : t(t0), it(0), first(true)
+  { dt = dtn = CFL = 0.0; }
 };
 
 #endif

src/SIM/SIMbase.C

@@ -22,6 +22,7 @@
 #endif
 #include "IntegrandBase.h"
 #include "AlgEqSystem.h"
+#include "ReactionsOnly.h"
 #include "SystemMatrix.h"
 #include "LinSolParams.h"
 #include "EigSolver.h"
@@ -2479,3 +2480,26 @@ void SIMbase::dumpEqSys ()
       utl::zero_print_tol = old_tol;
     }
 }
+
+
+bool SIMbase::assembleForces (const Vector& solution, double t0,
+                              Vector* R, Vector* S)
+{
+  if (!myProblem) return false;
+
+  // Assign secondary integral to the integrand
+  myProblem->setSecondaryInt(new ReactionsOnly(mySam,adm,R,S));
+
+  // Temporarily nullify myEqSys such that the secondary integral will be used
+  AlgEqSystem* tmpEqSys = myEqSys;
+  myEqSys = nullptr;
+
+  // Integrate the reaction and/or interface forces
+  bool ok = this->setMode(SIM::RHS_ONLY) && this->assembleSystem(t0,{solution});
+
+  // Release the secondary integral and restore the system matrix pointer
+  myProblem->setSecondaryInt();
+  myEqSys = tmpEqSys;
+
+  return ok;
+}

src/SIM/SIMbase.h

@@ -278,11 +278,12 @@ public:
                               bool newLHSmatrix = true, bool poorConvg = false);
 
   //! \brief Administers assembly of the linear equation system.
+  //! \param[in] t0 Time for evaluation of time-dependent property functions
   //! \param[in] pSol Primary solution vectors in DOF-order
   //!
   //! \details Use this version for linear/stationary problems only.
-  bool assembleSystem(const Vectors& pSol = Vectors())
-  { return this->assembleSystem(TimeDomain(),pSol); }
+  bool assembleSystem(double t0 = 0.0, const Vectors& pSol = Vectors())
+  { return this->assembleSystem(TimeDomain(t0),pSol); }
 
   //! \brief Extracts the assembled load vector for inspection/visualization.
   //! \param[out] loadVec Global load vector in DOF-order
@@ -464,20 +465,31 @@ public:
   //! \param[in] ncv Number of Arnoldi vectors (see ARPack documentation)
   //! \param[in] shift Eigenvalue shift
   //! \param[out] solution Computed eigenvalues and associated eigenvectors
-  //! \param[in] iA Index of system matrix \b A in \a myEqSys->A
-  //! \param[in] iB Index of system matrix \b B in \a myEqSys->A
+  //! \param[in] iA Index of system matrix \b A in \ref myEqSys
+  //! \param[in] iB Index of system matrix \b B in \ref myEqSys
   bool systemModes(std::vector<Mode>& solution,
                    int nev, int ncv, int iop, double shift,
                    size_t iA = 0, size_t iB = 1);
   //! \brief Performs a generalized eigenvalue analysis of the assembled system.
   //! \param[out] solution Computed eigenvalues and associated eigenvectors
-  //! \param[in] iA Index of system matrix \b A in \a myEqSys->A
-  //! \param[in] iB Index of system matrix \b B in \a myEqSys->A
+  //! \param[in] iA Index of system matrix \b A in \ref myEqSys
+  //! \param[in] iB Index of system matrix \b B in \ref myEqSys
   bool systemModes(std::vector<Mode>& solution, size_t iA = 0, size_t iB = 1)
   {
     return this->systemModes(solution,opt.nev,opt.ncv,opt.eig,opt.shift,iA,iB);
   }
 
+  //! \brief Returns whether reaction forces are to be computed or not.
+  virtual bool haveBoundaryReactions(bool = false) const { return false; }
+
+  //! \brief Assembles reaction and interface forces for specified boundaries.
+  //! \param[in] solution Current primary solution vector
+  //! \param[in] t0 Current time (for time-dependent loads)
+  //! \param[out] R Nodal reaction force container
+  //! \param[out] S Nodal interface force container
+  bool assembleForces(const Vector& solution, double t0,
+                      Vector* R, Vector* S = nullptr);
+
 
   // Post-processing methods
   // =======================

src/SIM/SIMoutput.h

@@ -297,6 +297,8 @@ public:
   virtual double getEffectivityIndex(const Vectors&, size_t, size_t) const = 0;
   //! \brief Prints integrated solution norms to the log stream.
   virtual void printNorms(const Vectors&, size_t = 36) const = 0;
+  //! \brief Prints out interface force resultants to log stream.
+  virtual void printIFforces(const Vector&, RealArray&) {}
 
   //! \brief Writes out the additional functions to VTF-file.
   virtual bool writeAddFuncs(int iStep, int& nBlock, int idBlock, double time);