changed: reorder operations in GlbNorm to make them work with multithreaded assembly

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1405 e10b68d5-8a6e-419e-a041-bce267b0401d

src/ASM/ElmMats.h

@@ -70,9 +70,6 @@ public:
     return b.front();
   }
 
-  //! \brief Virtual destruction method to clean up after numerical integration.
-  virtual void destruct() { delete this; }
-
   std::vector<Matrix> A; //!< The element coefficient matrices
   std::vector<Vector> b; //!< The element right-hand-side vectors
 

src/ASM/ElmNorm.h

@@ -14,8 +14,8 @@
 #ifndef _ELM_NORM_H
 #define _ELM_NORM_H
 
-#include <cstddef>
 #include "LocalIntegral.h"
+#include <cstddef>
 
 
 /*!
@@ -29,7 +29,17 @@ class ElmNorm : public LocalIntegral
 {
 public:
   //! \brief The constructor assigns the internal pointer.
+  //! \param[in] p Pointer to element norm values
+  //! \param[in] n Number of norm values
   ElmNorm(double* p, size_t n) : ptr(p), nnv(n) {}
+  //! \brief Alternative constructor using the internal buffer \a buf.
+  //! \param[in] n Number of norm values
+  //!
+  //! \details This constructor is used when the element norms are not requested
+  //! by the application, but are only used to assembly the global norms.
+  //! To avoid the need for a global array of element norms in that case,
+  //! an internal array is then used instead.
+  ElmNorm(size_t n) : buf(n,0.0), nnv(n) { ptr = &buf.front(); }
   //! \brief Empty destructor.
   virtual ~ElmNorm() {}
 
@@ -41,15 +51,32 @@ public:
   //! \brief Returns the number of norm values.
   size_t size() const { return nnv; }
 
+  //! \brief Returns whether the element norms are stored externally or not.
+  bool externalStorage() const { return buf.empty(); }
+
   //! \brief Virtual destruction method to clean up after numerical integration.
-  //! \details For this class we only clear the two solution vector containers
-  //! (to save memory). We do NOT delete the object itself, since it might be
-  //! needed in a second integration loop over element boundaries, for instance.
-  virtual void destruct() { vec.clear(); psol.clear(); }
+  //! \details Unless the internal buffer \a buf is used, these method only
+  //! clears the two solution vector containers (to save memory).
+  //! The object itself is then NOT deleted, since it might be used in a second
+  //! integration loop over element boundaries, for instance.
+  virtual void destruct()
+  {
+    if (buf.empty())
+    {
+      vec.clear();
+      psol.clear();
+    }
+    else
+      delete this; // The internal buffer has been used, delete the whole thing
+  }
 
 private:
-  double* ptr; //!< Pointer to the actual norm values
-  size_t  nnv; //!< Number of norm values
+  RealArray buf; //!< Internal buffer used when element norms are not requested
+  double*   ptr; //!< Pointer to the actual norm values
+  size_t    nnv; //!< Number of norm values
+
+public:
+  Vectors  psol; //!< Element-level projected solution vectors
 };
 
 #endif

src/ASM/GlbNorm.C

@@ -34,6 +34,13 @@ bool GlbNorm::assemble (const LocalIntegral* elmObj, int elmId)
   const ElmNorm* ptr = dynamic_cast<const ElmNorm*>(elmObj);
   if (!ptr) return false;
 
+  // If the element norms are requested (i.e. the internal buffer
+  // is not used) the actuall summation of element norms into the
+  // global norm is postponed to the very end (when invoked with
+  // elmId=0). This must be done like this to allow more than one
+  // loop over the elements during the norm integration.
+  if (elmId > 0 && ptr->externalStorage()) return true;
+
   ElmNorm& elVals = *const_cast<ElmNorm*>(ptr);
   for (size_t i = 0; i < elVals.size(); i++)
   {

src/ASM/GlbNorm.h

@@ -40,7 +40,7 @@ public:
   //! \brief Adds element norm quantities into the global norm object.
   //! \param[in] elmObj Pointer to the element norms to add into \a *this
   //! \param[in] elmId Global number of the element associated with \a *elmObj
-  virtual bool assemble(const LocalIntegral* elmObj, int elmId);
+  virtual bool assemble(const LocalIntegral* elmObj, int elmId = 0);
 
 private:
   //! \brief Applies the operation \a myOp on the given \a value.

src/ASM/GlobalIntegral.h

@@ -34,7 +34,7 @@ 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
-  virtual bool assemble(const LocalIntegral* elmObj, int elmId) = 0;
+  virtual bool assemble(const LocalIntegral* elmObj, int elmId) { return true; }
 };
 
 #endif

src/ASM/IntegrandBase.h

@@ -22,6 +22,8 @@ class NormBase;
 class AnaSol;
 class VTF;
 
+typedef std::vector<LocalIntegral*> LintegralVec;
+
 /*!
   \brief Base class representing a system level integrated quantity.
 */
@@ -259,6 +261,7 @@ protected:
   Vectors prjsol; //!< Projected secondary solution vectors for current patch
 
   unsigned short int nrcmp; //!< Number of projected solution components
+  LintegralVec* lints; //!< Vector of local integrals used during norm integration
 };
 
 #endif

src/SIM/SIMbase.C

@@ -1139,6 +1139,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
 			     const Vectors& psol, const Vectors& ssol,
 			     Vector& gNorm, Matrix* eNorm)
 {
+  PROFILE1("Norm integration");
+
   NormBase* norm = myProblem->getNormIntegrand(mySol);
   if (!norm)
   {
@@ -1148,8 +1150,6 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
     return false;
   }
 
-  PROFILE1("Norm integration");
-
   myProblem->initIntegration(time);
   const Vector& primsol = psol.front();
 
@@ -1162,6 +1162,14 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
   size_t nNorms = norm->getNoFields();
   gNorm.resize(nNorms,true);
 
+#ifdef USE_OPENMP
+  // When assembling in parallel, we must always do the norm summation
+  // at the end in a serial loop, to avoid that the threads try to update the
+  // same memory address simultaneously.
+  Matrix dummy;
+  if (!eNorm) eNorm = &dummy;
+#endif
+
   // Initialize norm integral classes
   GlbNorm globalNorm(gNorm,GlbNorm::SQRT);
   LintegralVec elementNorms;
@@ -1171,7 +1179,9 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
     elementNorms.reserve(eNorm->cols());
     for (i = 0; i < eNorm->cols(); i++)
       elementNorms.push_back(new ElmNorm(eNorm->ptr(i),nNorms));
+    norm->setLocalIntegrals(&elementNorms);
   }
+  norm->setLocalIntegrals(&elementNorms);
 
   // Loop over the different material regions, integrating solution norm terms
   // for the patch domain associated with each material
@@ -1206,12 +1216,6 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
       lp = i+1;
     }
 
-  // Integrate norm contributions due to Neumann boundary conditions, if any.
-  // Note: Currently, only the global norms are considered here.
-  // The corresponding element-level norms are not stored. This is mainly
-  // because the current design only allows one loop over the elements
-  // in the element-norm calculations. Consider rivising this later.
-
   if (norm->hasBoundaryTerms())
     for (i = 0; i < myProps.size() && ok; i++)
       if (myProps[i].pcode == Property::NEUMANN)
@@ -1242,8 +1246,18 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
 	  else
 	    ok = false;
 
+  // Clean up the dynamically allocated norm objects. This will also perform
+  // the actual global norm assembly, in case the element norms are stored,
+  // and always when multi-threading is used.
+  for (i = 0; i < elementNorms.size(); i++)
+  {
+    globalNorm.assemble(elementNorms[i]);
+    delete elementNorms[i];
+  }
+  delete norm;
+
   // Add problem-dependent external norm contributions
-  if (gNorm.size() >= 2)
+  if (gNorm.size() > 1)
     gNorm(2) += this->externalEnergy(psol);
 
 #ifdef PARALLEL_PETSC
@@ -1256,10 +1270,6 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
   }
 #endif
 
-  delete norm;
-  for (i = 0; i < elementNorms.size(); i++)
-    delete elementNorms[i];
-
   return ok;
 }