use mc33 in tdcs driver

duneuro/common/udg_solver.hh

@@ -20,7 +20,7 @@ namespace duneuro
   template <class ST, int compartments, int degree, class P, class DF, class RF, class JF>
   struct UDGSolverTraits {
     using SubTriangulation = ST;
-    using FundamentalGridView = typename ST::BaseT::GridView;
+    using FundamentalGridView = typename ST::GridView;
     static const int dimension = FundamentalGridView::dimension;
     using Problem = P;
     using FunctionSpace = UDGQkMultiPhaseSpace<FundamentalGridView, RF, degree, compartments>;

duneuro/meeg/udg_meeg_driver.hh

@@ -52,7 +52,9 @@ namespace duneuro
         , domain_(levelSetGridView_, data.levelSetData, config.sub("domain"))
         , subTriangulation_(std::make_shared<typename Traits::SubTriangulation>(
               fundamentalGridView_, levelSetGridView_, domain_.getDomainConfiguration(),
-              config.get<bool>("udg.force_refinement", false)))
+              config.get<bool>("udg.force_refinement", false),
+              config.get<std::string>("udg.tpmc_type", "full") == "full" ? tpmc::fullTPMC :
+                                                                           tpmc::simpleMC))
         , elementSearch_(std::make_shared<typename Traits::ElementSearch>(fundamentalGridView_))
         , solver_(
               std::make_shared<typename Traits::Solver>(subTriangulation_, config.sub("solver")))

duneuro/tes/udg_tdcs_driver.hh

@@ -3,7 +3,8 @@
 
 #include <dune/common/parametertree.hh>
 
-#include <dune/udg/simpletpmctriangulation.hh>
+#include <dune/udg/mc33triangulation.hh>
+#include <dune/udg/pdelab/gridfunction.hh>
 
 #include <duneuro/common/make_dof_vector.hh>
 #include <duneuro/common/structured_grid_utilities.hh>
@@ -13,7 +14,7 @@
 #include <duneuro/tes/tdcs_driver_interface.hh>
 #include <duneuro/tes/tdcs_patch_udg_parameter.hh>
 #include <duneuro/tes/udg_tdcs_solver.hh>
-#include <duneuro/udg/simpletpmc_domain.hh>
+#include <duneuro/udg/mc33_domain.hh>
 
 namespace duneuro
 {
@@ -21,7 +22,8 @@ namespace duneuro
   struct UDGTDCSDriverTraits {
     using Grid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim>>;
     using GridView = typename Grid::LevelGridView;
-    using SubTriangulation = Dune::UDG::SimpleTpmcTriangulation<GridView, GridView>;
+    using DD = Dune::UDG::DomainDefinition<GridView, GridView>;
+    using SubTriangulation = Dune::UDG::MarchingCube33SubTriangulation<DD>;
     using Problem = TDCSPatchUDGParameter<GridView>;
     using Solver = UDGSolver<SubTriangulation, compartments, degree, Problem>;
 
@@ -47,18 +49,17 @@ namespace duneuro
         , grid_(make_structured_grid<dim>(config.sub("volume_conductor.grid")))
         , fundamentalGridView_(grid_->levelGridView(0))
         , levelSetGridView_(grid_->levelGridView(grid_->maxLevel()))
-        , domain_(levelSetGridView_, data.levelSetData, config.sub("domain"))
+        , levelSets_(make_level_sets(levelSetGridView_, config.sub("domain")))
+        , domain_(fundamentalGridView_, levelSetGridView_, levelSets_, config.sub("domain"))
         , subTriangulation_(std::make_shared<typename Traits::SubTriangulation>(
-              fundamentalGridView_, levelSetGridView_, domain_.getDomainConfiguration(),
-              config.get<bool>("udg.force_refinement", false),
-              config.get<std::string>("udg.tpmc_type", "full") == "full" ? tpmc::fullTPMC :
-                                                                           tpmc::simpleMC))
+              domain_.getDomainConfiguration(), config.sub("udg")))
         , problem_(std::make_shared<typename Traits::Problem>(
               config_.get<std::vector<double>>("solver.conductivities"), patchSet))
         , solver_(std::make_shared<typename Traits::Solver>(subTriangulation_, problem_,
                                                             config.sub("solver")))
         , conductivities_(config.get<std::vector<double>>("solver.conductivities"))
     {
+      subTriangulation_->finalize();
     }
 
     virtual std::unique_ptr<Function> makeDomainFunction() const override
@@ -89,12 +90,47 @@ namespace duneuro
                 fundamentalGridView_, conductivities_));
         vtkWriter.addVertexData(std::make_shared<Dune::UDG::DomainIndexUnfittedVTKGridFunction<
                                     typename Traits::GridView>>(fundamentalGridView_));
+        vtkWriter.addVertexData(std::make_shared<Dune::UDG::HostCellIndexUnfittedVTKGridFunction<
+                                    typename Traits::GridView>>(fundamentalGridView_));
         auto modeString = config.get<std::string>("mode", "volume");
         if ((modeString == "faces") || (modeString == "boundary")) {
           vtkWriter.addVertexData(std::make_shared<Dune::UDG::DomainIndexUnfittedVTKGridFunction<
                                       typename Traits::GridView>>(fundamentalGridView_, false));
         }
         vtkWriter.write(config, dataTree);
+
+        using UGF = Dune::PDELab::UnfittedDiscreteGridFunction<
+            typename Traits::Solver::Traits::FunctionSpace::GFS, typename Traits::DomainDOFVector,
+            typename Traits::SubTriangulation>;
+        UGF ugf(solver_->functionSpace().getGFS(),
+                function.cast<typename Traits::DomainDOFVector>(), *subTriangulation, 0, true);
+        std::vector<std::pair<double, double>> entries;
+        for (const auto& element : Dune::elements(fundamentalGridView_)) {
+          for (const auto& intersection : Dune::intersections(fundamentalGridView_, element)) {
+            const auto& ig = intersection.geometry();
+            if (intersection.boundary()) {
+              bool found = true;
+              for (unsigned int i = 0; i < ig.corners(); ++i) {
+                if (ig.corner(i)[0] < 1. - 1e-6) {
+                  found = false;
+                  break;
+                }
+              }
+              if (found) {
+                typename UGF::Traits::RangeType y;
+                for (unsigned int i = 0; i < ig.corners(); ++i) {
+                  ugf.evaluate(ig.inside(), intersection.inside().geometry().local(ig.corner(i)),
+                               y);
+                  entries.push_back({ig.corner(i)[1], y});
+                }
+              }
+            }
+          }
+        }
+        std::ofstream of("test_right_side.csv");
+        for (const auto& p: entries) {
+          of << p.first << " " << p.second << "\n";
+        }
       } else {
         DUNE_THROW(Dune::Exception, "Unknown format \"" << format << "\"");
       }
@@ -113,6 +149,8 @@ namespace duneuro
                 fundamentalGridView_, conductivities_));
         vtkWriter.addVertexData(std::make_shared<Dune::UDG::DomainIndexUnfittedVTKGridFunction<
                                     typename Traits::GridView>>(fundamentalGridView_));
+        vtkWriter.addVertexData(std::make_shared<Dune::UDG::HostCellIndexUnfittedVTKGridFunction<
+                                    typename Traits::GridView>>(fundamentalGridView_));
         auto modeString = config.get<std::string>("mode", "volume");
         if ((modeString == "faces") || (modeString == "boundary")) {
           vtkWriter.addVertexData(std::make_shared<Dune::UDG::DomainIndexUnfittedVTKGridFunction<
@@ -129,7 +167,8 @@ namespace duneuro
     std::unique_ptr<typename Traits::Grid> grid_;
     typename Traits::GridView fundamentalGridView_;
     typename Traits::GridView levelSetGridView_;
-    SimpleTPMCDomain<typename Traits::GridView, typename Traits::GridView> domain_;
+    std::vector<std::shared_ptr<MC33LevelSetGridFunction<typename Traits::GridView>>> levelSets_;
+    MC33Domain<typename Traits::GridView, typename Traits::GridView> domain_;
     std::shared_ptr<typename Traits::SubTriangulation> subTriangulation_;
     std::shared_ptr<typename Traits::Problem> problem_;
     std::shared_ptr<typename Traits::Solver> solver_;