From a4b38fff9de21c295ccc0e7579fa031ac482953a Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Tue, 21 Dec 2021 19:10:50 +0100
Subject: [PATCH] added: support for Integrand::AVERAGE flag in LR integrands

---
 src/ASM/LR/ASMu2D.C | 26 ++++++++++++++++++++++++++
 src/ASM/LR/ASMu3D.C | 32 ++++++++++++++++++++++++++++++++
 2 files changed, 58 insertions(+)

diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C
index 4d33a4b9..569bc772 100644
--- a/src/ASM/LR/ASMu2D.C
+++ b/src/ASM/LR/ASMu2D.C
@@ -1187,6 +1187,32 @@ bool ASMu2D::integrate (Integrand& integrand,
         dXidu[1] = el->vmax() - el->vmin();
       }
 
+      if (integrand.getIntegrandType() & Integrand::AVERAGE)
+      {
+        // --- Compute average value of basis functions over the element -----
+
+        fe.Navg.resize(el->support().size(), true);
+        double area = 0.0;
+        int ip = (iel-1)*nGP*nGP;
+        for (int j = 0; j < nGP; j++)
+          for (int i = 0; i < nGP; i++, ip++)
+          {
+            SplineUtils::extractBasis(spline1[ip],fe.N,dNdu);
+
+            // Compute Jacobian determinant of coordinate mapping
+            // and multiply by weight of current integration point
+            double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
+            double weight = dA*wg[i]*wg[j];
+
+            // Numerical quadrature
+            fe.Navg.add(fe.N,detJac*weight);
+            area += detJac*weight;
+          }
+
+        // Divide by element area
+        fe.Navg /= area;
+      }
+
       // Initialize element quantities
       LocalIntegral* A = integrand.getLocalIntegral(el->support().size(),fe.iel);
       if (!integrand.initElement(MNPC[iel-1],fe,X,nRed*nRed,*A))
diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C
index 14f18673..379d7c39 100644
--- a/src/ASM/LR/ASMu3D.C
+++ b/src/ASM/LR/ASMu3D.C
@@ -1089,6 +1089,38 @@ bool ASMu3D::integrate (Integrand& integrand,
         dXidu[2] = el->getParmin(2);
       }
 
+      if (integrand.getIntegrandType() & Integrand::AVERAGE)
+      {
+        // --- Compute average value of basis functions over the element -----
+
+        int ip = (iel-1)*nGP*nGP*nGP + firstIp;
+        fe.Navg.resize(el->support().size(),true);
+        double vol = 0.0;
+        int ig = 1;
+        for (int k = 0; k < nGP; k++)
+          for (int j = 0; j < nGP; j++)
+            for (int i = 0; i < nGP; i++, ++ip, ++ig)
+            {
+              B.fillColumn(1, BN.getColumn(ig));
+              B.fillColumn(2, BdNdu.getColumn(ig)*2.0/du);
+              B.fillColumn(3, BdNdv.getColumn(ig)*2.0/dv);
+              B.fillColumn(4, BdNdw.getColumn(ig)*2.0/dw);
+              this->evaluateBasis(fe, dNdu, C, B);
+
+              // Compute Jacobian determinant of coordinate mapping
+              // and multiply by weight of current integration point
+              double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
+              double weight = dV*wg[i]*wg[j]*wg[k];
+
+              // Numerical quadrature
+              fe.Navg.add(fe.N,detJac*weight);
+              vol += detJac*weight;
+            }
+
+        // Divide by element volume
+        fe.Navg /= vol;
+      }
+
       // Initialize element quantities
       LocalIntegral* A = integrand.getLocalIntegral(el->support().size(),fe.iel);
       if (!integrand.initElement(MNPC[iel-1],fe,X,nRed*nRed*nRed,*A))