[Problem] 3D: Use an asymmetric weak patch

src/enumparser.cc

@@ -56,6 +56,17 @@ Config::stateModel StringToEnum<Config::stateModel>::convert(
   DUNE_THROW(Dune::Exception, "failed to parse enum");
 }
 
+Config::PatchType StringToEnum<Config::PatchType>::convert(
+    std::string const &s) {
+  if (s == "Rectangular")
+    return Config::Rectangular;
+
+  if (s == "Trapezoidal")
+    return Config::Trapezoidal;
+
+  DUNE_THROW(Dune::Exception, "failed to parse enum");
+}
+
 Solver::VerbosityMode StringToEnum<Solver::VerbosityMode>::convert(
     std::string const &s) {
   if (s == "full")
@@ -73,4 +84,5 @@ Solver::VerbosityMode StringToEnum<Solver::VerbosityMode>::convert(
 template std::istream &operator>>(std::istream &lhs, Config::FrictionModel &);
 template std::istream &operator>>(std::istream &lhs, Config::stateModel &);
 template std::istream &operator>>(std::istream &lhs, Config::scheme &);
+template std::istream &operator>>(std::istream &lhs, Config::PatchType &);
 template std::istream &operator>>(std::istream &lhs, Solver::VerbosityMode &);

src/enumparser.hh

@@ -23,6 +23,10 @@ template <> struct StringToEnum<Config::scheme> {
   static Config::scheme convert(std::string const &s);
 };
 
+template <> struct StringToEnum<Config::PatchType> {
+  static Config::PatchType convert(std::string const &s);
+};
+
 template <> struct StringToEnum<Solver::VerbosityMode> {
   static Solver::VerbosityMode convert(std::string const &s);
 };

src/enums.hh

@@ -5,6 +5,7 @@ struct Config {
   enum FrictionModel { Truncated, Regularised };
   enum stateModel { AgeingLaw, SlipLaw };
   enum scheme { Newmark, BackwardEuler };
+  enum PatchType { Rectangular, Trapezoidal };
 };
 
 #endif

src/sand-wedge-data/parset-3D.cfg

@@ -2,6 +2,9 @@
 [boundary.friction]
 smallestDiameter= 2e-2  # [m]
 
+[boundary.friction.weakening]
+patchType       = Trapezoidal
+
 [timeSteps]
 refinementTolerance = 1e-5
 

src/sand-wedge-data/weakpatch.hh

 #ifndef SRC_SAND_WEDGE_DATA_WEAKPATCH_HH
 #define SRC_SAND_WEDGE_DATA_WEAKPATCH_HH
 
-template <class LocalVector> ConvexPolyhedron<LocalVector> getWeakPatch() {
+template <class LocalVector>
+ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset) {
   ConvexPolyhedron<LocalVector> weakPatch;
 #if MY_DIM == 3
   weakPatch.vertices.resize(4);
@@ -11,6 +12,16 @@ template <class LocalVector> ConvexPolyhedron<LocalVector> getWeakPatch() {
     weakPatch.vertices[k][2] = -MyGeometry::depth / 2.0;
     weakPatch.vertices[k + 2][2] = MyGeometry::depth / 2.0;
   }
+  switch (parset.get<Config::PatchType>("patchType")) {
+    case Config::Rectangular:
+      break;
+    case Config::Trapezoidal:
+      weakPatch.vertices[1][0] += 0.05 * MyGeometry::lengthScale;
+      weakPatch.vertices[3][0] -= 0.05 * MyGeometry::lengthScale;
+      break;
+    default:
+      assert(false);
+  }
 #else
   weakPatch.vertices.resize(2);
   weakPatch.vertices[0] = MyGeometry::X;

src/sand-wedge.cc

@@ -119,7 +119,8 @@ int main(int argc, char *argv[]) {
     using ScalarMatrix = MyAssembler::ScalarMatrix;
     using ScalarVector = MyAssembler::ScalarVector;
 
-    auto weakPatch = getWeakPatch<LocalVector>();
+    auto const weakPatch =
+        getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"));
 
     // {{{ Set up grid
     GridConstructor<Grid> gridConstructor;
@@ -140,8 +141,8 @@ int main(int argc, char *argv[]) {
     std::cout << "min diameter: " << minDiameter << std::endl;
     std::cout << "max diameter: " << maxDiameter << std::endl;
 
-    GridView const leafView = grid->leafGridView();
-    size_t const leafVertexCount = leafView.size(dims);
+    auto const leafView = grid->leafGridView();
+    auto const leafVertexCount = leafView.size(dims);
 
     std::cout << "Number of DOFs: " << leafVertexCount << std::endl;