From e2521f6b28d289d43e673409eb1f2ce693765c9d Mon Sep 17 00:00:00 2001
From: kmo <kmo@e10b68d5-8a6e-419e-a041-bce267b0401d>
Date: Thu, 23 Feb 2012 12:02:09 +0000
Subject: [PATCH] Added projection of secondary solutions for the Poisson
 solver

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1482 e10b68d5-8a6e-419e-a041-bce267b0401d
---
 src/Integrands/Poisson.C | 60 ++++++++++++++++++++++++++++++++++------
 src/Integrands/Poisson.h |  4 +--
 2 files changed, 53 insertions(+), 11 deletions(-)

diff --git a/src/Integrands/Poisson.C b/src/Integrands/Poisson.C
index 2fd02563..f89ffa21 100644
--- a/src/Integrands/Poisson.C
+++ b/src/Integrands/Poisson.C
@@ -262,6 +262,23 @@ NormBase* Poisson::getNormIntegrand (AnaSol* asol) const
 }
 
 
+PoissonNorm::PoissonNorm (Poisson& p, VecFunc* a) : NormBase(p), anasol(a)
+{
+  nrcmp = myProblem.getNoFields(2);
+}
+
+
+size_t PoissonNorm::getNoFields () const
+{
+  size_t nf = anasol ? 4 : 2;
+  for (size_t i = 0; i < prjsol.size(); i++)
+    if (!prjsol.empty())
+      nf += anasol ? 3 : 2;
+
+  return nf;
+}
+
+
 bool PoissonNorm::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
 			   const Vec3& X) const
 {
@@ -269,31 +286,56 @@ bool PoissonNorm::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
   ElmNorm& pnorm = static_cast<ElmNorm&>(elmInt);
 
   // Evaluate the inverse constitutive matrix at this point
-  Matrix C;
-  problem.formCmatrix(C,X,true);
+  Matrix Cinv;
+  problem.formCmatrix(Cinv,X,true);
 
   // Evaluate the finite element heat flux field
-  Vector sigmah;
-  if (!problem.evalSol(sigmah,pnorm.vec.front(),fe.dNdX,X)) return false;
+  Vector sigmah, sigma, error;
+  if (!problem.evalSol(sigmah,pnorm.vec.front(),fe.dNdX,X))
+    return false;
 
   // Integrate the energy norm a(u^h,u^h)
-  pnorm[0] += sigmah.dot(C*sigmah)*fe.detJxW;
+  pnorm[0] += sigmah.dot(Cinv*sigmah)*fe.detJxW;
   // Evaluate the temperature field
   double u = pnorm.vec.front().dot(fe.N);
   // Integrate the external energy (h,u^h)
   pnorm[1] += problem.getHeat(X)*u*fe.detJxW;
 
+  size_t ip = 2;
   if (anasol)
   {
     // Evaluate the analytical heat flux
-    Vector sigma((*anasol)(X).ptr(),sigmah.size());
+    sigma.fill((*anasol)(X).ptr(),nrcmp);
     // Integrate the energy norm a(u,u)
-    pnorm[2] += sigma.dot(C*sigma)*fe.detJxW;
+    pnorm[ip++] += sigma.dot(Cinv*sigma)*fe.detJxW;
     // Integrate the error in energy norm a(u-u^h,u-u^h)
-    sigma -= sigmah;
-    pnorm[3] += sigma.dot(C*sigma)*fe.detJxW;
+    error = sigma - sigmah;
+    pnorm[ip++] += error.dot(Cinv*error)*fe.detJxW;
   }
 
+  size_t i, j;
+  for (i = 0; i < pnorm.psol.size(); i++)
+    if (!pnorm.psol[i].empty())
+    {
+      // Evaluate projected heat flux field
+      Vector sigmar(nrcmp);
+      for (j = 0; j < nrcmp; j++)
+	sigmar[j] = pnorm.psol[i].dot(fe.N,j,nrcmp);
+
+      // Integrate the energy norm a(u^r,u^r)
+      pnorm[ip++] += sigmar.dot(Cinv*sigmar)*fe.detJxW;
+      // Integrate the error in energy norm a(u^r-u^h,u^r-u^h)
+      error = sigmar - sigmah;
+      pnorm[ip++] += error.dot(Cinv*error)*fe.detJxW;
+
+      if (anasol)
+      {
+	// Integrate the error in the projected solution a(u-u^r,u-u^r)
+	error = sigma - sigmar;
+	pnorm[ip++] += error.dot(Cinv*error)*fe.detJxW;
+      }
+    }
+
   return true;
 }
 
diff --git a/src/Integrands/Poisson.h b/src/Integrands/Poisson.h
index 610b23e4..2374ed32 100644
--- a/src/Integrands/Poisson.h
+++ b/src/Integrands/Poisson.h
@@ -153,7 +153,7 @@ public:
   //! \brief The only constructor initializes its data members.
   //! \param[in] p The Poisson problem to evaluate norms for
   //! \param[in] a The analytical heat flux (optional)
-  PoissonNorm(Poisson& p, VecFunc* a = 0) : NormBase(p), anasol(a) {}
+  PoissonNorm(Poisson& p, VecFunc* a = 0);
   //! \brief Empty destructor.
   virtual ~PoissonNorm() {}
 
@@ -161,7 +161,7 @@ public:
   virtual bool hasBoundaryTerms() const { return true; }
 
   //! \brief Returns the number of norm quantities.
-  virtual size_t getNoFields() const { return anasol ? 4 : 2; }
+  virtual size_t getNoFields() const;
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions