From 04961a91b1aa14a08fbe4c36334288915c91d5ef Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Tue, 6 Mar 2018 17:09:10 +0100
Subject: [PATCH] .

---
 src/assemblers_tmpl.cc                      |  2 +-
 src/explicitgrid.hh                         |  8 ++
 src/gridselector.hh                         | 30 +++++++
 src/hdf5/patchinfo-writer_tmpl.cc           |  8 +-
 src/hdf5/surface-writer_tmpl.cc             |  2 +-
 src/multi-body-problem-data/mygrids_tmpl.cc |  6 +-
 src/multi-body-problem.cc                   | 93 ++++++++++++---------
 7 files changed, 99 insertions(+), 50 deletions(-)

diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc
index c2384ddc..e068ad5e 100644
--- a/src/assemblers_tmpl.cc
+++ b/src/assemblers_tmpl.cc
@@ -4,4 +4,4 @@
 
 #include "explicitgrid.hh"
 
-template class MyAssembler<GridView, MY_DIM>;
+template class MyAssembler<LeafGridView, MY_DIM>;
diff --git a/src/explicitgrid.hh b/src/explicitgrid.hh
index e69de29b..b538e24d 100644
--- a/src/explicitgrid.hh
+++ b/src/explicitgrid.hh
@@ -0,0 +1,8 @@
+#ifndef SRC_EXPLICITGRID_HH
+#define SRC_EXPLICITGRID_HH
+
+#include "gridselector.hh"
+
+using LeafGridView = Grid::LeafGridView;
+using LevelGridView = Grid::LevelGridView;
+#endif
diff --git a/src/gridselector.hh b/src/gridselector.hh
index e69de29b..94fb9d85 100644
--- a/src/gridselector.hh
+++ b/src/gridselector.hh
@@ -0,0 +1,30 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#define WANT_ALUGRID 0
+#define WANT_UG 1
+
+#define WANT_GRID WANT_UG
+
+#if WANT_GRID == WANT_ALUGRID
+
+#if !HAVE_ALUGRID
+#error ALUGRID was requested but not found
+#endif
+#include <dune/grid/alugrid.hh>
+using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>;
+
+#elif WANT_GRID == WANT_UG
+
+#if !HAVE_UG
+#error UG was requested but not found
+#endif
+#include <dune/grid/uggrid.hh>
+using Grid = Dune::UGGrid<MY_DIM>;
+
+#else
+
+#error requested a grid that does not exist
+
+#endif
diff --git a/src/hdf5/patchinfo-writer_tmpl.cc b/src/hdf5/patchinfo-writer_tmpl.cc
index ef346e10..d2f3c014 100644
--- a/src/hdf5/patchinfo-writer_tmpl.cc
+++ b/src/hdf5/patchinfo-writer_tmpl.cc
@@ -4,11 +4,11 @@
 #include "../program_state.hh"
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<GridView, double>;
+using P1Basis = P1NodalBasis<LeafGridView, double>;
 using MyFunction = BasisGridFunction<P1Basis, Vector>;
 
-template class GridEvaluator<LocalVector, GridView>;
+template class GridEvaluator<LocalVector, LeafGridView>;
 template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
+GridEvaluator<LocalVector, LeafGridView>::evaluate(MyFunction const &f) const;
 
-template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;
+template class PatchInfoWriter<MyProgramState, P1Basis, LeafGridView>;
diff --git a/src/hdf5/surface-writer_tmpl.cc b/src/hdf5/surface-writer_tmpl.cc
index e4570aa7..516b59a7 100644
--- a/src/hdf5/surface-writer_tmpl.cc
+++ b/src/hdf5/surface-writer_tmpl.cc
@@ -5,4 +5,4 @@
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 
-template class SurfaceWriter<MyProgramState, GridView>;
+template class SurfaceWriter<MyProgramState, LeafGridView>;
diff --git a/src/multi-body-problem-data/mygrids_tmpl.cc b/src/multi-body-problem-data/mygrids_tmpl.cc
index 99ae7cd6..2fee2d0e 100644
--- a/src/multi-body-problem-data/mygrids_tmpl.cc
+++ b/src/multi-body-problem-data/mygrids_tmpl.cc
@@ -8,11 +8,11 @@
 
 template class GridsConstructor<Grid>;
 
-template struct MyFaces<GridView>;
+template struct MyFaces<LeafGridView>;
 template struct MyFaces<LevelGridView>;
 
-template MyFaces<GridView> GridsConstructor<Grid>::constructFaces(
-    GridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+template MyFaces<LeafGridView> GridsConstructor<Grid>::constructFaces(
+    LeafGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
 
 template MyFaces<LevelGridView> GridsConstructor<Grid>::constructFaces(
     LevelGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 7533b579..83d79227 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -120,9 +120,15 @@ int main(int argc, char *argv[]) {
 
     auto const parset = getParameters(argc, argv);
 
-    using MyAssembler = MyAssembler<LeafGridView, dims>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, dims>>;
+    using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
+    using DeformedGrid = DeformedGridComplex::DeformedGridType;
+
+    using DefLeafGridView = DeformedGrid::LeafGridView;
+    using DefLevelGridView = DeformedGrid::LevelGridView;
+
+    using MyAssembler = MyAssembler<DefLeafGridView, dims>;
     using Matrix = MyAssembler::Matrix;
-    using Vector = MyAssembler::Vector;
     using LocalVector = Vector::block_type;
     using ScalarMatrix = MyAssembler::ScalarMatrix;
     using ScalarVector = MyAssembler::ScalarVector;
@@ -154,32 +160,18 @@ int main(int argc, char *argv[]) {
     }
 #endif
 
-    // set up (deformed) grids
+    // set up grids
     GridsConstructor<Grid> gridConstructor(cuboidGeometries);
     auto& grids = gridConstructor.getGrids();
 
-    using DeformedGridComplex = DeformedContinuaComplex<Grid, Vector>;
-    using DeformedGrid = DeformedGridComplex::DeformedGridType;
-    Dune::Contact::DeformedContinuaComplex<Grid, Vector> deformedGrids(grids);
+    std::vector<ConvexPolyhedron<LocalVector>> weakPatches(grids.size());
 
-    using LeafGridView = DeformedGrid::LeafGridView;
-    using LevelGridView = DeformedGrid::LevelGridView;
-    std::vector<std::shared_ptr<LeafGridView>> leafViews(deformedGrids.size());
-    std::vector<std::shared_ptr<LevelGridView>> levelViews(deformedGrids.size());
-    std::vector<int> leafVertexCounts(deformedGrids.size());
-
-    using LeafFaces = MyFaces<LeafGridView>;
-    using LevelFaces = MyFaces<LevelGridView>;
-    std::vector<std::shared_ptr<LeafFaces>> leafFaces(deformedGrids.size());
-    std::vector<std::shared_ptr<LevelFaces>> levelFaces(deformedGrids.size());
-
-    std::vector<ConvexPolyhedron<LocalVector>> weakPatches(deformedGrids.size());
-
-    for (size_t i=0; i<deformedGrids.size(); i++) {
-        auto const &cuboidGeometry = *cuboidGeometries[i];
+    for (size_t i=0; i<grids.size(); i++) {
+        const auto& cuboidGeometry = *cuboidGeometries[i];
 
         // define weak patch and refine grid
-        weakPatches[i] = getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"), cuboidGeometry);
+        auto& weakPatch = weakPatches[i];
+        weakPatch = getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"), cuboidGeometry);
         refine(*grids[i], weakPatch, parset.get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
 
         writeToVTK(grids[i]->leafGridView(), "", "grid" + std::to_string(i));
@@ -195,19 +187,38 @@ int main(int argc, char *argv[]) {
         }
         std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
         std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
+    }
 
-        leafViews[i] = std::make_shared<LeafGridView>(grids[i]->leafGridView());
-        levelViews[i] = std::make_shared<LevelGridView>(grids[i]->levelGridView(0));
+    // set up deformed grids
+    DeformedGridComplex deformedGridComplex;
+    for (size_t i=0; i<grids.size(); i++) {
+        deformedGridComplex.addGrid(*grids[i], "body" + std::to_string(i));
+    }
+    const auto deformedGrids = deformedGridComplex.gridVector();
+
+    std::vector<std::shared_ptr<DefLeafGridView>> leafViews(deformedGrids.size());
+    std::vector<std::shared_ptr<DefLevelGridView>> levelViews(deformedGrids.size());
+    std::vector<int> leafVertexCounts(deformedGrids.size());
+
+    using LeafFaces = MyFaces<DefLeafGridView>;
+    using LevelFaces = MyFaces<DefLevelGridView>;
+    std::vector<std::shared_ptr<LeafFaces>> leafFaces(deformedGrids.size());
+    std::vector<std::shared_ptr<LevelFaces>> levelFaces(deformedGrids.size());
+
+    for (size_t i=0; i<deformedGrids.size(); i++) {
+        leafViews[i] = std::make_shared<DefLeafGridView>(deformedGrids[i]->leafGridView());
+        levelViews[i] = std::make_shared<DefLevelGridView>(deformedGrids[i]->levelGridView(0));
 
         leafVertexCounts[i] = leafViews[i]->size(dims);
 
+        const auto& cuboidGeometry = *cuboidGeometries[i];
         leafFaces[i] = std::make_shared<LeafFaces>(*leafViews[i], cuboidGeometry);
         levelFaces[i] = std::make_shared<LevelFaces>(*levelViews[i], cuboidGeometry);
     }
 
     // set up contact couplings
-    using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
-    using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+    using LeafBoundaryPatch = BoundaryPatch<DefLeafGridView>;
+    using LevelBoundaryPatch = BoundaryPatch<DefLevelGridView>;
 
     using CouplingPair = Dune::Contact::CouplingPair<DeformedGrid, DeformedGrid, field_type>;
     std::vector<std::shared_ptr<CouplingPair>> couplings(bodyCount-1);
@@ -222,13 +233,13 @@ int main(int argc, char *argv[]) {
     }
 
     // set up dune-contact nBodyAssembler
-    std::vector<const DeformedGrid*> const_grids(bodyCount);
+    /*std::vector<const DeformedGrid*> const_grids(bodyCount);
     for (size_t i=0; i<const_grids.size(); i++) {
-        const_grids[i] = grids[i].get();
-    }
+        const_grids[i] = &deformedGrids.grid("body" + std::to_string(i));
+    }*/
 
     Dune::Contact::NBodyAssembler<DeformedGrid, Vector> nBodyAssembler(bodyCount, bodyCount-1);
-    nBodyAssembler.setGrids(const_grids);
+    nBodyAssembler.setGrids(deformedGrids);
 
     for (size_t i=0; i<couplings.size(); i++) {
         nBodyAssembler.setCoupling(*couplings[i], i);
@@ -238,9 +249,9 @@ int main(int argc, char *argv[]) {
     nBodyAssembler.assembleObstacle();
 
     // set up boundary conditions
-    std::vector<BoundaryPatch<LeafGridView>> neumannBoundaries(bodyCount);
-    std::vector<BoundaryPatch<LeafGridView>> frictionalBoundaries(bodyCount);
-    std::vector<BoundaryPatch<LeafGridView>> surfaces(bodyCount);
+    std::vector<BoundaryPatch<DefLeafGridView>> neumannBoundaries(bodyCount);
+    std::vector<BoundaryPatch<DefLeafGridView>> frictionalBoundaries(bodyCount);
+    std::vector<BoundaryPatch<DefLeafGridView>> surfaces(bodyCount);
 
     // Dirichlet boundary
     std::vector<Dune::BitSetVector<dims>> noNodes(bodyCount);
@@ -258,17 +269,17 @@ int main(int argc, char *argv[]) {
 
       // Neumann boundary
       auto& neumannBoundary = neumannBoundaries[i];
-      neumannBoundary = BoundaryPatch<LeafGridView>(*leafViews[i]);
+      neumannBoundary = BoundaryPatch<DefLeafGridView>(*leafViews[i]);
 
       // friction boundary
       auto& frictionalBoundary = frictionalBoundaries[i];
       if (i==0) {
-        frictionalBoundary = BoundaryPatch<LeafGridView>(leafFaces[i]->upper);
+        frictionalBoundary = BoundaryPatch<DefLeafGridView>(leafFaces[i]->upper);
       } else if (i==bodyCount-1) {
-        frictionalBoundary = BoundaryPatch<LeafGridView>(leafFaces[i]->lower);
+        frictionalBoundary = BoundaryPatch<DefLeafGridView>(leafFaces[i]->lower);
       } else {
-          frictionalBoundary = BoundaryPatch<LeafGridView>(leafFaces[i]->lower);
-          frictionalBoundary.addPatch(BoundaryPatch<LeafGridView>(leafFaces[i]->upper));
+          frictionalBoundary = BoundaryPatch<DefLeafGridView>(leafFaces[i]->lower);
+          frictionalBoundary.addPatch(BoundaryPatch<DefLeafGridView>(leafFaces[i]->upper));
       }
       //surfaces[i] = BoundaryPatch<GridView>(myFaces.upper);
 
@@ -395,7 +406,7 @@ int main(int argc, char *argv[]) {
     std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction(weakPatches.size());
     for (size_t i=0; i<frictionInfo.size(); i++) {
         frictionInfo[i] = std::make_shared<MyGlobalFrictionData<LocalVector>>(parset.sub("boundary.friction"), weakPatches[i]);
-        globalFriction[i] = assemblers[i]->assembleFrictionNonlinearity(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), frictionalBoundary, frictionInfo, programState.weightedNormalStress);
+        globalFriction[i] = assemblers[i]->assembleFrictionNonlinearity(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), frictionalBoundaries[i], *frictionInfo[i], programState.weightedNormalStress[i]);
         globalFriction[i]->updateAlpha(programState.alpha[i]);
     }
 
@@ -453,7 +464,7 @@ int main(int argc, char *argv[]) {
     report(true);
     */
 
-
+/*
     // Set up TNNMG solver
     using NonlinearFactory = SolverFactory<
         dims,
@@ -514,7 +525,7 @@ int main(int argc, char *argv[]) {
         break;
       }
     }
-
+*/
   } catch (Dune::Exception &e) {
     Dune::derr << "Dune reported error: " << e << std::endl;
   } catch (std::exception &e) {
-- 
GitLab