pass operators by value

dune/amdis/Operator.hpp

@@ -42,16 +42,16 @@ namespace AMDiS
      * evaluation of the scalar product `c = C(phi(lambda))*v(lambda)` with lambda a local coordinate.
      **/
     template <class Term>
-    Self& addZOT(Term const& c)
+    Operator& addZOT(Term const& c)
     {
       return addZOTImpl(toTerm(c));
     }
 
     template <class... Args>
-    static shared_ptr<Self> zot(Args&&... args)
+    static Operator zot(Args&&... args)
     {
-      auto op = make_shared<Self>();
-      op->addZOT(std::forward<Args>(args)...);
+      Self op;
+      op.addZOT(std::forward<Args>(args)...);
       return op;
     }
 
@@ -61,7 +61,7 @@ namespace AMDiS
     /// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
     /// The coefficient \p b must be a vector expression.
     template <class Term>
-    Self& addFOT(Term const& b, FirstOrderType type)
+    Operator& addFOT(Term const& b, FirstOrderType type)
     {
       return addFOTImpl(toTerm(b), type);
     }
@@ -71,16 +71,16 @@ namespace AMDiS
     /// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
     /// The coefficient \p b must be a scalar expression.
     template <class Term>
-    Self& addFOT(Term const& b, std::size_t i, FirstOrderType type)
+    Operator& addFOT(Term const& b, std::size_t i, FirstOrderType type)
     {
       return addFOTImpl(toTerm(b), i, type);
     }
 
     template <class... Args>
-    static shared_ptr<Self> fot(Args&&... args)
+    static Operator fot(Args&&... args)
     {
-      auto op = make_shared<Self>();
-      op->addFOT(std::forward<Args>(args)...);
+      Self op;
+      op.addFOT(std::forward<Args>(args)...);
       return op;
     }
 
@@ -88,7 +88,7 @@ namespace AMDiS
     /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
     /// where \p A can be a matrix or a scalar expression.
     template <class Term>
-    Self& addSOT(Term const& A)
+    Operator& addSOT(Term const& A)
     {
       return addSOTImpl(toTerm(A));
     }
@@ -98,16 +98,16 @@ namespace AMDiS
     /// corresponds to the derivative component of the test-function and the
     /// second index \p _j to the derivative component of the trial function.
     template <class Term>
-    Self& addSOT(Term const& A, std::size_t i, std::size_t j)
+    Operator& addSOT(Term const& A, std::size_t i, std::size_t j)
     {
       return addSOTImpl(toTerm(A), i, j);
     }
 
     template <class... Args>
-    static shared_ptr<Self> sot(Args&&... args)
+    static Operator sot(Args&&... args)
     {
-      auto op = make_shared<Self>();
-      op->addSOT(std::forward<Args>(args)...);
+      Self op;
+      op.addSOT(std::forward<Args>(args)...);
       return op;
     }
 

dune/amdis/ProblemStat.hpp

@@ -142,10 +142,8 @@ namespace AMDiS
     void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
                            double* factor = nullptr, double* estFactor = nullptr);
 
-    void addMatrixOperator(std::map< std::pair<int,int>,
-                                     std::shared_ptr<ElementOperator> > ops);
-    void addMatrixOperator(std::map< std::pair<int,int>,
-                                     std::pair<std::shared_ptr<ElementOperator>, bool> > ops);
+    void addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops);
+    void addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator, bool> > ops);
     /** @} */
 
 
@@ -159,8 +157,8 @@ namespace AMDiS
     void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i,
                            double* factor = nullptr, double* estFactor = nullptr);
 
-    void addVectorOperator(std::map< int, std::shared_ptr<ElementOperator> > ops);
-    void addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops);
+    void addVectorOperator(std::map< int, ElementOperator> ops);
+    void addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops);
     /** @} */
 
     /// Adds a Dirichlet boundary condition

dune/amdis/ProblemStat.inc.hpp

@@ -126,24 +126,24 @@ namespace AMDiS
 
   template <class Traits>
   void ProblemStatSeq<Traits>::
-  addMatrixOperator(std::map< std::pair<int,int>, std::shared_ptr<ElementOperator> > ops)
+  addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
   {
     for (auto op : ops){
       int r = op.first.first;
       int c = op.first.second;
-      matrixOperators[r][c].emplace_back(op.second);
+      matrixOperators[r][c].emplace_back(std::make_shared<ElementOperator>(op.second));
       matrixChanging[r][c] = true;
     }
   }
 
   template <class Traits>
   void ProblemStatSeq<Traits>::
-  addMatrixOperator(std::map< std::pair<int,int>, std::pair<std::shared_ptr<ElementOperator>,bool> > ops)
+  addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
   {
     for (auto op : ops) {
       int r = op.first.first;
       int c = op.first.second;
-      matrixOperators[r][c].emplace_back(op.second.first);
+      matrixOperators[r][c].emplace_back(std::make_shared<ElementOperator>(op.second.first));
       matrixChanging[r][c] = matrixChanging[r][c] || op.second.second;
     }
   }
@@ -168,20 +168,20 @@ namespace AMDiS
 
   template <class Traits>
   void ProblemStatSeq<Traits>::
-  addVectorOperator(std::map< int, std::shared_ptr<ElementOperator> > ops)
+  addVectorOperator(std::map< int, ElementOperator> ops)
   {
     for (auto op : ops) {
-      vectorOperators[op.first].emplace_back(op.second);
+      vectorOperators[op.first].emplace_back(std::make_shared<ElementOperator>(op.second));
       vectorChanging[op.first] = true;
     }
   }
 
   template <class Traits>
   void ProblemStatSeq<Traits>::
-  addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops)
+  addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
   {
     for (auto op : ops) {
-      vectorOperators[op.first].emplace_back(op.second.first);
+      vectorOperators[op.first].emplace_back(std::make_shared<ElementOperator>(op.second.first));
       vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
     }
   }

src/MTLPfcPrecon.hpp

@@ -3,8 +3,8 @@
 #include "BlockPreconditioner.hpp"
 #include <dune/amdis/common/Literals.hpp>
 
-namespace AMDiS 
-{ 
+namespace AMDiS
+{
   template <class Matrix, class Vector>
   class MTLPfcPrecon : public BlockPreconditioner<Matrix, Vector>
   {
@@ -12,104 +12,104 @@ namespace AMDiS
     using Self = MTLPfcPrecon;
     using Super = BlockPreconditioner<Matrix, Vector>;
     using PreconBase = typename Super::PreconBase;
-    
+
     using MTLMatrix = BaseMatrix_t<Matrix>;
     using MTLVector = Vector;
     using LinearSolverType = LinearSolverInterface<MTLMatrix, MTLVector>;
-    
+
   public:
     struct Creator : CreatorInterface<PreconBase>
     {
       virtual std::shared_ptr<PreconBase> create() override
       {
         precon = std::make_shared<Self>();
-        return precon; 
+        return precon;
       }
-      
+
       std::shared_ptr<Self> precon = NULL;
-    };  
-    
+    };
+
   public:
     MTLPfcPrecon()
       : Super()
     {
-      std::string solverNameM = "cg", 
-                  solverNameMpL = "cg", 
+      std::string solverNameM = "cg",
+                  solverNameMpL = "cg",
                   solverNameMpL2 = "cg";
       Parameters::get("precon_pfc->M->solver", solverNameM);
       Parameters::get("precon_pfc->MpL->solver", solverNameMpL);
       Parameters::get("precon_pfc->MpL2->solver", solverNameMpL2);
-      
+
       CreatorInterfaceName<LinearSolverType>* creator;
-      
+
       creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameM, "precon_pfc->M->solver") );
       solverM = creator->create("precon_pfc->M");
-      
+
       creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL, "precon_pfc->MpL->solver") );
       solverMpL = creator->create("precon_pfc->MpL");
-      
+
       creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL2, "precon_pfc->MpL2->solver") );
       solverMpL2 = creator->create("precon_pfc->MpL2");
     }
-    
+
     virtual ~MTLPfcPrecon()
     {
       exit();
     }
-      
+
     void setData(double* tau_, double M0_ = 1.0)
     {
       tau = tau_;
       M0 = M0_;
     }
-    
+
     // implementation in MTLPfcPrecon.inc.hpp
     virtual void init(Matrix const& A) override;
-    
+
     virtual void exit() override {}
-    
+
     // implementation in MTLPfcPrecon.inc.hpp
-    virtual void solve(MTLVector const& b, MTLVector& x) const override;  
-      
-                       
+    virtual void solve(MTLVector const& b, MTLVector& x) const override;
+
+
     template <typename VectorX, typename VectorB>
     void solveM(VectorX& x, const VectorB& b) const
     {
       SolverInfo solverInfo("precon_pfc->M");
       solverM->solve(getM(), x, b, solverInfo);
     }
-    
+
     template <typename VectorX, typename VectorB>
     void solveMpL(VectorX& x, const VectorB& b) const
     {
       SolverInfo solverInfo("precon_pfc->MpL");
       solverMpL->solve(MpL, x, b, solverInfo);
     }
-    
+
     template <typename VectorX, typename VectorB>
     void solveMpL2(VectorX& x, const VectorB& b) const
     {
       SolverInfo solverInfo("precon_pfc->MpL2");
       solverMpL2->solve(MpL2, x, b, solverInfo);
     }
-    
-    MTLMatrix const& getM()  const { return matM  ? *matM  : (*Super::A)(2_c, 2_c); } // < u, v >
-    MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(1_c, 0_c); } // < M0*tau*grad(u), grad(v) >
-    MTLMatrix const& getL()  const { return matL  ? *matL  : (*Super::A)(2_c, 1_c); } // < grad(u), grad(v) >
-    
+
+    MTLMatrix const& getM()  const { return matM  ? *matM  : (*Super::A)(_2, _2); } // < u, v >
+    MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(_1, _0); } // < M0*tau*grad(u), grad(v) >
+    MTLMatrix const& getL()  const { return matL  ? *matL  : (*Super::A)(_2, _1); } // < grad(u), grad(v) >
+
     double getTau() const { return *tau; }
-    
-  protected:     
+
+  protected:
     double* tau = NULL;
     double M0 = 1.0;
-    
+
     MTLMatrix* matM = NULL;
     MTLMatrix* matL = NULL;
     MTLMatrix* matL0 = NULL;
-    
+
     MTLMatrix MpL;
     MTLMatrix MpL2;
-    
+
     mutable MTLVector y0;
     mutable MTLVector y1;
     mutable MTLVector tmp;
@@ -118,7 +118,7 @@ namespace AMDiS
     std::shared_ptr<LinearSolverType> solverMpL;
     std::shared_ptr<LinearSolverType> solverMpL2;
   };
-  
+
 } // end namespace AMDiS
 
 #include "MTLPfcPrecon.inc.hpp"

src/ellipt.cc

@@ -13,9 +13,9 @@
 using namespace AMDiS;
 
 // 1 component with polynomial degree 1
-using ElliptParam   = DefaultTraitsMesh<Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>, 1>;
+using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
+using ElliptParam   = DefaultTraitsMesh<Grid, 1>;
 using ElliptProblem = ProblemStat<ElliptParam>;
-using ElliptProblemInstat = ProblemInstat<ElliptParam>;
 
 int main(int argc, char** argv)
 {
@@ -24,9 +24,6 @@ int main(int argc, char** argv)
   ElliptProblem prob("ellipt");
   prob.initialize(INIT_ALL);
 
-  ElliptProblemInstat probInstat("ellipt", prob);
-  probInstat.initialize(INIT_UH_OLD);
-
   AdaptInfo adaptInfo("adapt");
 
   using Op = ElliptProblem::ElementOperator;
@@ -35,7 +32,7 @@ int main(int argc, char** argv)
   opL.addSOT(1.0);
   prob.addMatrixOperator(opL, 0, 0);
 
-  opForce.addZOT( eval([](auto const& x) { return -1.0; }) );
+  opForce.addZOT([](auto const& x) { return -1.0; });
   prob.addVectorOperator(opForce, 0);
 
 

src/gridview.cc

@@ -25,13 +25,12 @@ int main(int argc, char** argv)
 {
   AMDiS::init(argc, argv);
 
-
   // Create grid from structured grid
   std::array<unsigned int, 2> n = {{4, 4}};
   Dune::FieldVector<double, 2> lower = {{0.0, 0.0}};
   Dune::FieldVector<double, 2> upper = {{1.0, 1.0}};
   auto grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, n);
-
+  // NOTE: can not be used with AlbertaGrid
 
   ElliptProblem prob("ellipt", *grid);
   prob.initialize(INIT_ALL);
@@ -63,13 +62,6 @@ int main(int argc, char** argv)
   *prob.getSolution() = 0.0; // maybe not necessary
 
   prob.buildAfterCoarsen(adaptInfo, Flag(0));
-
-  // write matrix to file
-  {
-    mtl::io::matrix_market_ostream out("matrix.mtx");
-    out << prob.getSystemMatrix()->getMatrix<0,0>();
-  }
-
   prob.solve(adaptInfo);
   prob.writeFiles(adaptInfo, true);
 

src/heat.cc

@@ -58,10 +58,10 @@ int main(int argc, char** argv)
   opL.addSOT( constant(1.0) );
   prob.addMatrixOperator(opL, 0, 0);
 
-  opTimeRhs.addZOT( valueOf(prob.getSolution(0_c)) );
+  opTimeRhs.addZOT( valueOf(prob.getSolution(_0)) );
   prob.addVectorOperator(opTimeRhs, 0, probInstat.getInvTau());
 
-  opForce.addZOT( eval([](auto const& x) { return -1.0; }) );
+  opForce.addZOT([](auto const& x) { return -1.0; });
   prob.addVectorOperator(opForce, 0);
 
 

src/navier_stokes.cc

@@ -50,38 +50,39 @@ int main(int argc, char** argv)
     using Op = StokesProblem::ElementOperator;
     forEach(range_<0,DOW>, [&](auto const _i)
     {
-	// <1/tau * u_i, v_i>
+	      // <1/tau * u_i, v_i>
         auto opTime = Op::zot( density );
         auto opTimeOld = Op::zot( valueOf(prob.getSolution(_i), density) );
-        prob.addMatrixOperator(*opTime, _i, _i, probInstat.getInvTau());
-        prob.addVectorOperator(*opTimeOld,  _i, probInstat.getInvTau());
+        prob.addMatrixOperator(opTime, _i, _i, probInstat.getInvTau());
+        prob.addVectorOperator(opTimeOld,  _i, probInstat.getInvTau());
 
-#if STRATEGY == 1
         // <(u * nabla)u_i^old, v_i>
         forEach(range_<0, DOW>, [&](auto const _j)
         {
           auto opNonlin = Op::zot( derivativeOf(prob.getSolution(_i), _j, density) );
-          prob.addMatrixOperator(*opNonlin, _i, _j);
+          prob.addMatrixOperator(opNonlin, _i, _j);
         });
-#elif STRATEGY == 2
+
         // <(u^old * nabla)u_i, v_i>
         forEach(range_<0, DOW>, [&](auto const _j)
         {
           auto opNonlin = Op::fot( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
-          prob.addMatrixOperator(*opNonlin, _i, _i);
+          prob.addMatrixOperator(opNonlin, _i, _i);
         });
-#else
+
         // <(u^old * nabla)u_i^old, v_i>
         forEach(range_<0, DOW>, [&](auto const _j)
         {
-          auto opNonlin = Op::zot( density * valueOf(prob.getSolution(_j)) * derivativeOf(prob.getSolution(_i), _j) );
-          prob.addVectorOperator(*opNonlin, _i);
+          auto opNonlin = Op::zot( func([density](double uj, double dj_ui)
+            {
+              return density * uj * dj_ui;
+            }, valueOf(prob.getSolution(_j)), derivativeOf(prob.getSolution(_i), _j)) );
+          prob.addVectorOperator(opNonlin, _i);
         });
-#endif
 
-	auto opL = std::make_pair( Op::sot( viscosity ), false ); // <viscosity*grad(u_i), grad(v_i)>
-	auto opP = std::make_pair( Op::fot( 1.0, _i, GRD_PSI ), false); // <p, d_i(v_i)>
-	auto opDiv = std::make_pair( Op::fot( 1.0, _i, GRD_PHI ), false); // <d_i(u_i), q>
+        auto opL = std::make_pair( Op::sot( viscosity ), false ); // <viscosity*grad(u_i), grad(v_i)>
+        auto opP = std::make_pair( Op::fot( 1.0, _i, GRD_PSI ), false); // <p, d_i(v_i)>
+        auto opDiv = std::make_pair( Op::fot( 1.0, _i, GRD_PHI ), false); // <d_i(u_i), q>
 
         const int i = _i;
         prob.addMatrixOperator({{
@@ -105,7 +106,7 @@ int main(int argc, char** argv)
 
     // set boundary conditions
     for (size_t i = 0; i < DOW; ++i)
-	prob.addDirichletBC(box, i, i, zero);
+	    prob.addDirichletBC(box, i, i, zero);
 
     prob.addDirichletBC(left, 0, 0, zero);
     prob.addDirichletBC(left, 1, 1, parabolic);

src/pfc.cc

@@ -57,8 +57,8 @@ int main(int argc, char** argv)
 
     // < [-(1+r) - 3*psi^2]*u, v > + < 2*grad(u), grad(v) >
     auto opLhs01_ = Op::zot( -(1.0 + r) );
-    opLhs01_->addZOT( valueOfFunc(Psi, [](auto psi) { return -3.0 * Math::pow<2>(psi); }, 2) );
-    opLhs01_->addSOT( 2.0 );
+    opLhs01_.addZOT( valueOfFunc(Psi, [](auto psi) { return -3.0 * Math::pow<2>(psi); }, 2) );
+    opLhs01_.addSOT( 2.0 );
     auto opLhs01 = std::make_pair(opLhs01_, true);
 
     // < -2*psi^3, v >

src/stokes.cc

@@ -31,24 +31,24 @@ int main(int argc, char** argv)
     using Op = StokesProblem::ElementOperator;
     for (int i = 0; i < DOW; ++i) {
         // <viscosity*grad(u_i), grad(v_i)>
-        Op* opL = new Op;
-        opL->addSOT( constant(viscosity) );
-        prob.addMatrixOperator(*opL, i, i);
+        Op opL;
+        opL.addSOT( constant(viscosity) );
+        prob.addMatrixOperator(opL, i, i);
 
         // <p, d_i(v_i)>
-        Op* opB = new Op;
-        opB->addFOT( constant(1.0), i, GRD_PSI );
-        prob.addMatrixOperator(*opB, i, DOW);
+        Op opB;
+        opB.addFOT( constant(1.0), i, GRD_PSI );
+        prob.addMatrixOperator(opB, i, DOW);
 
         // <d_i(u_i), q>
-        Op* opDiv = new Op;
-        opDiv->addFOT( constant(1.0), i, GRD_PHI );
-        prob.addMatrixOperator(*opDiv, DOW, i);
+        Op opDiv;
+        opDiv.addFOT( constant(1.0), i, GRD_PHI );
+        prob.addMatrixOperator(opDiv, DOW, i);
     }
 
-    Op* opZero = new Op;
-    opZero->addZOT( constant(0.0) );
-    prob.addMatrixOperator(*opZero, DOW, DOW);
+    Op* opZero;
+    opZero.addZOT( constant(0.0) );
+    prob.addMatrixOperator(opZero, DOW, DOW);
 
 
     // define boundary regions