3d generalization

dune/amdis/OperatorTerm.hpp

@@ -164,6 +164,24 @@ namespace AMDiS
   
 // -----------------------------------------------------------------------------
   
+  namespace Impl 
+  {
+    template <class F, int dow>
+    using CoordsFunctorResult
+      = std::result_of<F(Dune::FieldVector<double, dow>)>;
+      
+    template <class F, int dow>
+    static constexpr bool isCallableDow()
+    { return Dune::Functions::Concept::isCallable<F, Dune::FieldVector<double, dow>>(); }
+      
+    template <class F>
+    using CoordsFunctorResult_t = typename std::conditional_t<
+      isCallableDow<F, 1>(), CoordsFunctorResult<F, 1>, std::conditional_t<
+      isCallableDow<F, 2>(), CoordsFunctorResult<F, 2>, std::conditional_t<
+      isCallableDow<F, 3>(), CoordsFunctorResult<F, 3>, void> > >::type;
+  }
+  
+  
   
   /// An expression that evaluates to the current coordinate of a dof or 
   /// quadrature point with given index.
@@ -171,7 +189,7 @@ namespace AMDiS
   class CoordsTerm
   {
   public:
-    using value_type = typename std::result_of<Functor(Dune::FieldVector<double, 2>)>::type;
+    using value_type = Impl::CoordsFunctorResult_t<Functor>;
     
     template <class F,
       class = std::enable_if_t<std::is_same<Functor, std::decay_t<F>>::value> >

dune/amdis/TermGenerator.hpp

@@ -27,9 +27,13 @@ namespace AMDiS
     using TermCategory =
         std::conditional_t< std::is_arithmetic<T>::value,            
                           _constant,
+        std::conditional_t< Dune::Functions::Concept::isFunction<T, double(Dune::FieldVector<double, 1>)>(),   
+                      _functor,
         std::conditional_t< Dune::Functions::Concept::isFunction<T, double(Dune::FieldVector<double, 2>)>(),   
                       _functor,
-                      _default > >;
+        std::conditional_t< Dune::Functions::Concept::isFunction<T, double(Dune::FieldVector<double, 3>)>(),   
+                      _functor,
+                      _default > > > >;
 
 
     /// Generator class that takes a type as argument

init/ellipt.json.2d

@@ -3,27 +3,24 @@
 
   "elliptMesh": {
     "macro file name":    "./macro/macro.stand.2d",
-    "global refinements": 10
+    "global refinements": 0
   },
 
   "ellipt": {
     "mesh":  "elliptMesh",
-    "dim":                    2,
-    "components":             1,
-    "polynomial degree[0]":   1,
+    "names": "u",
  
-    "solver":                 "cg",
     "solver" : {
+      "name":               "cg",
       "max iteration":      1000,
-      "tolerance":      1e-8,
-      "info":           10,
+      "absolute tolerance": 1e-6,
+      "info":               1,
       "left precon":        "diag"
     },
 
     "output": {
       "filename":         "output/ellipt.2d",
-      "ParaView format": 1,
-      "ParaView mode": 1
+      "output directory": "output"
     }
   }
 }

src/CMakeLists.txt

 
-set(projects "heat" "pfc" "stokes" "navier_stokes")
+set(projects2d "ellipt" "heat" "pfc" "stokes" "navier_stokes")
 
-foreach(project ${projects})
+foreach(project ${projects2d})
     add_executable(${project} ${project}.cc)
     add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project})
     target_link_dune_default_libraries(${project})
     target_link_libraries(${project} "duneamdis")
     target_compile_definitions(${project} PRIVATE DIM=2 DOW=2)
 endforeach()
+
+set(projects3d "ellipt")
+
+foreach(project ${projects3d})
+    add_executable(${project} ${project}.cc)
+    add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 ${project})
+    target_link_dune_default_libraries(${project})
+    target_link_libraries(${project} "duneamdis")
+    target_compile_definitions(${project} PRIVATE DIM=3 DOW=3)
+endforeach()
\ No newline at end of file

src/ellipt.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <iostream>
+
+#include <dune/amdis/AMDiS.hpp>
+#include <dune/amdis/AdaptInstationary.hpp>
+#include <dune/amdis/ProblemInstat.hpp>
+#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/Literals.hpp>
+
+// 1 component with polynomial degree 1
+using ElliptParam   = ProblemStatTraits<DIM, DOW, 1>;
+using ElliptProblem = ProblemStat<ElliptParam>;
+using ElliptProblemInstat = ProblemInstat<ElliptParam>;
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+  
+  ElliptProblem prob("ellipt");
+  prob.initialize(INIT_ALL);
+  
+  ElliptProblemInstat probInstat("ellipt", prob);
+  probInstat.initialize(INIT_UH_OLD);
+  
+  AdaptInfo adaptInfo("adapt");
+  
+  using Op = ElliptProblem::OperatorType;
+  Op opL, opForce; 
+  
+  opL.addSOT(1.0);
+  prob.addMatrixOperator(opL, 0, 0);
+  
+  opForce.addZOT( eval([](auto const& x) { return -1.0; }) );
+  prob.addVectorOperator(opForce, 0);
+  
+  
+  // set boundary condition
+  auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
+  auto dbcValues = [](auto const& x){ return 0.0; }; // set value
+  prob.addDirichletBC(predicate, 0, 0, dbcValues);
+      
+  *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);
+  
+  AMDiS::finalize();
+  return 0;
+}
\ No newline at end of file

src/heat.cc

@@ -19,27 +19,27 @@
 using namespace AMDiS;
 
 
-class HeatParam
-{
-  template <class... Bs>
-  using FeSpaceTuple = std::tuple<Bs...>;
-  
-  template <class Mesh, int deg>
-  using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
-  
-public:
-  static constexpr int dim = DIM;
-  static constexpr int dimworld = DOW;
-  static constexpr int nComponents = 1;
-  
-  // default values
-  using Mesh = Dune::YaspGrid<dim>;
-  using FeSpaces = FeSpaceTuple<Lagrange<Mesh, 2>>;
-};
+// class HeatParam
+// {
+//   template <class... Bs>
+//   using FeSpaceTuple = std::tuple<Bs...>;
+//   
+//   template <class Mesh, int deg>
+//   using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
+//   
+// public:
+//   static constexpr int dim = DIM;
+//   static constexpr int dimworld = DOW;
+//   static constexpr int nComponents = 1;
+//   
+//   // default values
+//   using Mesh = Dune::YaspGrid<dim>;
+//   using FeSpaces = FeSpaceTuple<Lagrange<Mesh, 2>>;
+// };
 
 
 // 1 component with polynomial degree 1
-// using HeatParam   = ProblemStatTraits<DIM, DOW, 1>;
+using HeatParam   = ProblemStatTraits<DIM, DOW, 1>;
 using HeatProblem = ProblemStat<HeatParam>;
 using HeatProblemInstat = ProblemInstat<HeatParam>;