From ffd0307d3ee9108d65fe78eec821c8d4b3291c02 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 7 Jun 2010 14:49:46 +0000
Subject: [PATCH] Added more support for Rosenbrock method.

---
 AMDiS/ChangeLog                               |  5 ++
 AMDiS/src/time/RosenbrockAdaptInstationary.cc |  2 +
 AMDiS/src/time/RosenbrockStationary.cc        | 44 ++++++++++++++
 AMDiS/src/time/RosenbrockStationary.h         | 57 ++++++++++++++-----
 4 files changed, 94 insertions(+), 14 deletions(-)

diff --git a/AMDiS/ChangeLog b/AMDiS/ChangeLog
index 5b37c9d5..f122da0d 100644
--- a/AMDiS/ChangeLog
+++ b/AMDiS/ChangeLog
@@ -3,3 +3,8 @@ Rev 1210:
 	  "FeSpace" geändert. Dies betrifft insbesondere "getFESpace()",
 	  das nun "getFeSpace()" heißt.
 	* ProblemVec::getRHS() -> ProblemVec::getRhs()
+Rev 1231:
+	* Im Operator-Konstruktor brauchen die Flags, daher MATRIX_OPERATOR und
+	  VECTOR_OPERATOR nicht mehr angegeben werden. Sollten diese gesetzt sein,
+	  so compiliert das Programm noch, aber es wird eine entsprechende 
+	  Fehlermeldung ausgegeben und das Programm bricht ab.
diff --git a/AMDiS/src/time/RosenbrockAdaptInstationary.cc b/AMDiS/src/time/RosenbrockAdaptInstationary.cc
index 523235ff..f8a3aecb 100644
--- a/AMDiS/src/time/RosenbrockAdaptInstationary.cc
+++ b/AMDiS/src/time/RosenbrockAdaptInstationary.cc
@@ -32,6 +32,8 @@ namespace AMDiS {
     problemStat->setRosenbrockMethod(rosenbrockMethod);
 
     adaptInfo->setRosenbrockMode(true);
+
+    problemStat->setTauGamma(&tauGamma);
   }
 
 
diff --git a/AMDiS/src/time/RosenbrockStationary.cc b/AMDiS/src/time/RosenbrockStationary.cc
index 40b42e1e..c0f66018 100644
--- a/AMDiS/src/time/RosenbrockStationary.cc
+++ b/AMDiS/src/time/RosenbrockStationary.cc
@@ -31,6 +31,10 @@ namespace AMDiS {
   void RosenbrockStationary::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
 					       bool asmMatrix, bool asmVector)
   {        
+    FUNCNAME("RosenbrockStationary::buildAfterCoarsen()");
+
+    TEST_EXIT(tauPtr)("No tau pointer defined in stationary problem!\n");
+
     if (first) {
       first = false;
       *unVec = *solution;
@@ -79,4 +83,44 @@ namespace AMDiS {
   void RosenbrockStationary::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
   {}
 
+
+  void RosenbrockStationary::addOperator(Operator &op, int row, int col, 
+					 double *factor, double *estFactor)
+  {
+    FUNCNAME("RosenbrockStationary::addOperator()");
+
+    TEST_EXIT(op.getUhOld() == NULL)("UhOld is not allowed to be set!\n");
+
+    op.setUhOld(stageSolution->getDOFVector(col));
+    ProblemVec::addVectorOperator(op, row);
+  }
+  
+
+  void RosenbrockStationary::addJacobianOperator(Operator &op, int row, int col, 
+						 double *factor, double *estFactor)
+  {
+    FUNCNAME("RosenbrockStationary::addJacobianOperator()");
+    
+    TEST_EXIT(factor == NULL)("Not yet implemented!\n");
+    TEST_EXIT(estFactor == NULL)("Not yet implemented!\n");
+
+    ProblemVec::addMatrixOperator(op, row, col, &minusOne, &minusOne);
+  }
+
+
+  void RosenbrockStationary::addTimeOperator(int row, int col)
+  {
+    FUNCNAME("RosenbrockStationary::addTimeOperator()");
+
+    TEST_EXIT(tauGamma)("This should not happen!\n");
+
+    Operator *op = new Operator(componentSpaces[row], componentSpaces[col]);
+    op->addZeroOrderTerm(new Simple_ZOT);
+    ProblemVec::addMatrixOperator(op, row, col, tauGamma, tauGamma);
+
+    Operator *opRhs = new Operator(componentSpaces[row]);
+    opRhs->addZeroOrderTerm(new VecAtQP_ZOT(phiSum, new IdFunc()));
+    ProblemVec::addVectorOperator(opRhs, row, &minusOne, &minusOne);
+  }
+
 }
diff --git a/AMDiS/src/time/RosenbrockStationary.h b/AMDiS/src/time/RosenbrockStationary.h
index 4a69d968..ccc4bc05 100644
--- a/AMDiS/src/time/RosenbrockStationary.h
+++ b/AMDiS/src/time/RosenbrockStationary.h
@@ -32,29 +32,36 @@ namespace AMDiS {
   class RosenbrockStationary : public ProblemVec
   {
   public:
+    class IdFunc : public AbstractFunction<double, double>
+    {
+    public:
+      IdFunc()
+	: AbstractFunction<double, double>(1)
+      {}
+      
+      double operator()(const double& u) const 
+      {
+	return u;
+      }
+    };
+
+
     RosenbrockStationary(std::string name)
       : ProblemVec(name),
-	first(true)
+	first(true),
+	minusOne(-1.0),
+	tauPtr(NULL),
+	tauGamma(NULL)
     {}
 
-    DOFVector<double>* getStageSolution(int i)
-    {
-      return stageSolution->getDOFVector(i);
-    }
-
     DOFVector<double>* getUnVec(int i)
     {
       return unVec->getDOFVector(i);
     }
 
-    DOFVector<double>* getPhiSum()
-    {
-      return phiSum;
-    }
-
-    void setTauPtr(double *ptr)
+    DOFVector<double>* getStageSolution(int i)
     {
-      tauPtr = ptr;
+      return stageSolution->getDOFVector(i);
     }
 
     void acceptTimestep();
@@ -70,6 +77,24 @@ namespace AMDiS {
       init();
     }
 
+    void addOperator(Operator &op, int row, int col, 
+		     double *factor = NULL, double *estFactor = NULL);
+
+    void addJacobianOperator(Operator &op, int row, int col, 
+			     double *factor = NULL, double *estFactor = NULL);
+
+    void addTimeOperator(int i, int j);
+
+    void setTau(double *ptr)
+    {
+      tauPtr = ptr;
+    }
+
+    void setTauGamma(double *ptr)
+    {
+      tauGamma = ptr;
+    }
+
   protected:
     void init();
 
@@ -82,9 +107,13 @@ namespace AMDiS {
 
     DOFVector<double> *phiSum, *tmpDof;
 
+    bool first;
+
+    double minusOne;
+
     double *tauPtr;
 
-    bool first;
+    double *tauGamma;
   };
 
 }
-- 
GitLab