From 7736de71ca1703d10fc378efac2c4f908e2c7e9a Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Thu, 10 Feb 2022 13:22:47 +0100
Subject: [PATCH] Added: SIMgeneric::getInterfaceForces()

---
 src/SIM/SIMgeneric.C | 26 ++++++++++++++++++++++++++
 src/SIM/SIMgeneric.h |  8 ++++++++
 2 files changed, 34 insertions(+)

diff --git a/src/SIM/SIMgeneric.C b/src/SIM/SIMgeneric.C
index 4381963f..cb66a8af 100644
--- a/src/SIM/SIMgeneric.C
+++ b/src/SIM/SIMgeneric.C
@@ -17,6 +17,7 @@
 #include "IntegrandBase.h"
 #include "Utilities.h"
 #include "Function.h"
+#include "SAM.h"
 #include "IFEM.h"
 
 
@@ -59,6 +60,31 @@ Vector SIMgeneric::getSolution (const Vector& psol, const double* par,
 }
 
 
+Vector SIMgeneric::getInterfaceForces (const Vector& sf,
+                                       const RealArray& weights, int code) const
+{
+  Vector force(nsd);
+  if (!mySam)
+    return force;
+
+  IntVec glbNodes;
+  this->getBoundaryNodes(code,glbNodes);
+
+  for (int inod : glbNodes)
+  {
+    double w = inod <= (int)weights.size() ? weights[inod-1] : 1.0;
+    std::pair<int,int> dof = mySam->getNodeDOFs(inod);
+    for (unsigned char i = 0; i < nsd; i++, dof.first++)
+      if (dof.first <= dof.second && dof.first < (int)sf.size())
+        force[i] += w*sf(dof.first);
+      else
+        break;
+  }
+
+  return force;
+}
+
+
 int SIMgeneric::evalPoint (const double* xi, Vec3& X, double* param,
                            int patch, bool global) const
 {
diff --git a/src/SIM/SIMgeneric.h b/src/SIM/SIMgeneric.h
index bd996f59..648cbfe4 100644
--- a/src/SIM/SIMgeneric.h
+++ b/src/SIM/SIMgeneric.h
@@ -57,6 +57,14 @@ public:
   int evalPoint(const double* xi, Vec3& X, double* param = nullptr,
                 int patch = 1, bool global = false) const;
 
+  //! \brief Extracts the force resultant on a specified interface.
+  //! \param[in] sf Internal nodal forces
+  //! \param[in] weights Nodal weights
+  //! \param[in] code Code indentifying the interface to extract forces for
+  //! \return The interface force resultant
+  Vector getInterfaceForces(const Vector& sf,
+                            const RealArray& weights, int code) const;
+
   //! \brief Returns the element that contains a specified spatial point.
   //! \param[in] param The parameters of the point in the knot-span domain
   //! \param[in] patch 1-based patch index containing the point