Merge branch '12-add-prestress-for-fibers' into 'master'

Resolve "Add prestress for fibers"

Closes #12

See merge request !12

apps/CMakeLists.txt

@@ -13,7 +13,7 @@ foreach(exe ${APPS})
 endforeach()
 
 # Register tests
-foreach(num 23_1 23_2 23_3 24_1 24_2)
+foreach(num 23_1 23_2 23_3 23_3_prestress 24_1 24_2)
     add_test(NAME run_hansbo_${num}
         COMMAND structured2d_p2 run ${PROJECT_SOURCE_DIR}/experiments/hansbo_reproduction/fig${num}.yml
 )

dune/structures/eulerbernoulli.hh

@@ -14,6 +14,7 @@
 #include <dune/geometry/affinegeometry.hh>
 #include <dune/pdelab.hh>
 #include <dune/structures/elasticity.hh>
+#include <dune/structures/fiberparameters.hh>
 #include <dune/structures/material.hh>
 #include <dune/structures/parametrizedcurves.hh>
 
@@ -118,13 +119,8 @@ public:
       fibre_parametrizations.push_back(
         std::make_shared<StraightFibre<2>>(fibre));
 
-      // Parse the Youngs modulus of the fibre.
-      auto modulus = fibre["youngs_modulus"].as<double>();
-      fibre_modulus.push_back(modulus);
-
-      // Parse the fibre radii
-      auto radius = fibre["radius"].as<double>();
-      fibre_radii.push_back(radius);
+      // Construct fiber parameters from config node
+      fiber_parameters.emplace_back(fibre);
     }
     std::cout << "Parsed a total of " << fibre_parametrizations.size()
               << " fibres from the configuration." << std::endl;
@@ -275,11 +271,13 @@ public:
             force[k] += local_coefficient_force_vector(lfsv.child(k), i)
                         * basis.function(i);
 
-        // Extract physical parameter of the fibre
-        auto E = fibre_modulus[fibindex];
-        auto d = fibre_radii[fibindex];
-        auto A = d;
-        auto I = (d * d * d) / 12.0;
+        // Extract physical parameters of the fibre
+        const auto& fiber_param = fiber_parameters[fibindex];
+        const auto E = fiber_param.youngs_modulus;
+        const auto d = fiber_param.radius;
+        const auto prestress = fiber_param.prestress;
+        const auto A = d;
+        const auto I = (d * d * d) / 12.0;
 
         // The needed tangential derivative quantities. These expressions are
         // generated using the generate_tangential_derivatives Python script.
@@ -421,18 +419,18 @@ public:
                              + jit[0][0] * basis.hessian(i, 0, 0))))
                 * t[0];
 
-          r.accumulate(lfsu.child(0),
-                       i,
-                       (E * A * dtut * dtvt_0 + E * I * dt2un * dt2vn_0
-                        - A * force[0] * basis.function(i))
-                         * ip.weight()
-                         * global_linegeo.integrationElement(ip.position()));
-          r.accumulate(lfsu.child(1),
-                       i,
-                       (E * A * dtut * dtvt_1 + E * I * dt2un * dt2vn_1
-                        - A * force[1] * basis.function(i))
-                         * ip.weight()
-                         * global_linegeo.integrationElement(ip.position()));
+          r.accumulate(
+            lfsu.child(0),
+            i,
+            (E * A * dtut * dtvt_0 + E * I * dt2un * dt2vn_0
+             + A * prestress * dtvt_0 - A * force[0] * basis.function(i))
+              * ip.weight() * global_linegeo.integrationElement(ip.position()));
+          r.accumulate(
+            lfsu.child(1),
+            i,
+            (E * A * dtut * dtvt_1 + E * I * dt2un * dt2vn_1
+             + A * prestress * dtvt_1 - A * force[1] * basis.function(i))
+              * ip.weight() * global_linegeo.integrationElement(ip.position()));
         }
       }
     }
@@ -553,10 +551,10 @@ public:
       auto jit_n = cellgeo_n.jacobianInverseTransposed(pos_n);
 
       // Extract physical parameter of the fibre
-      auto E = fibre_modulus[fibindex];
-      auto d = fibre_radii[fibindex];
-      auto A = d;
-      auto I = (d * d * d) / 12.0;
+      const auto& fiber_param = fiber_parameters[fibindex];
+      const auto E = fiber_param.youngs_modulus;
+      const auto d = fiber_param.radius;
+      const auto I = (d * d * d) / 12.0;
 
       // Compute the penalty factor
       const auto h_F = std::min(cellgeo_s.volume(), cellgeo_n.volume())
@@ -907,8 +905,7 @@ private:
   std::vector<std::shared_ptr<FibreParametrizationBase<2>>>
     fibre_parametrizations;
   std::vector<FibreGridIntersection> fibre_intersections;
-  std::vector<double> fibre_modulus;
-  std::vector<double> fibre_radii;
+  std::vector<FiberParameters<double>> fiber_parameters;
 
   // Store a force vector - code adapted of how the generated code does it.
   // That is not the nicest way of doing it, but we do not need additional

dune/structures/fiberparameters.hh

+#ifndef DUNE_STRUCTURES_FIBERPARAMETERS_HH
+#define DUNE_STRUCTURES_FIBERPARAMETERS_HH
+
+#include <yaml-cpp/yaml.h>
+
+/// A storage for material parameters of a fiber
+template<typename Range = double>
+struct FiberParameters
+{
+  Range radius;
+  Range youngs_modulus;
+  Range prestress;
+
+  /// Default constructor
+  FiberParameters() = default;
+  /// Default destructor
+  ~FiberParameters() = default;
+
+  /// Constructor from explicit input values
+  FiberParameters(const Range radius_,
+                  const Range youngs_modulus_,
+                  const Range prestress_)
+    : radius(radius_)
+    , youngs_modulus(youngs_modulus_)
+    , prestress(prestress_)
+  {
+  }
+
+  /// Constructor from a YAML config node
+  FiberParameters(const YAML::Node& config)
+    : radius(config["radius"].as<Range>())
+    , youngs_modulus(config["youngs_modulus"].as<Range>())
+    , prestress(config["prestress"].as<Range>())
+  {
+  }
+};
+
+#endif // DUNE_STRUCTURES_FIBERPARAMETERS_HH

dune/structures/reinforcedoperator.hh

@@ -142,6 +142,11 @@ public:
                    "          default: 1e5                              \n"
                    "          meta:                                     \n"
                    "            title: Young's modulus                  \n"
+                   "        prestress:                                  \n"
+                   "          type: float                               \n"
+                   "          default: 0.0                              \n"
+                   "          meta:                                     \n"
+                   "            title: Isotropic fiber prestress        \n"
                    "    meta:                                           \n"
                    "      title: Fibres                                 \n");
     return data;

experiments/hansbo_reproduction/fig23_1.yml

@@ -53,7 +53,7 @@ solver:
           _blockname: vis_vonmises_0
           continuous: true
           vector: 0
-      filename: output
+      filename: fig23_1
       instationary: false
       intervals: 1
       path: vtk

experiments/hansbo_reproduction/fig23_2.yml

@@ -72,7 +72,7 @@ solver:
             _blockname: vis_vonmises_0
             continuous: true
             vector: 0
-        filename: output
+        filename: fig23_2
         instationary: false
         intervals: 1
         path: vtk

experiments/hansbo_reproduction/fig23_3.yml

@@ -71,7 +71,7 @@ solver:
             _blockname: vis_stress_1
             vector: 0
             index: 0
-        filename: output
+        filename: fig23_3
         instationary: false
         intervals: 1
         path: vtk

experiments/hansbo_reproduction/fig23_3_prestress.yml

+solver:
+    vectors:
+      - Displacement Field_0
+      - Body Force_0
+      - Traction Force_0
+    grid:
+      N:
+        - 40
+        - 10
+      lowerleft:
+        - 0
+        - 0
+      refinement: 0
+      upperright:
+        - 4
+        - 1
+      _blockname: grid_0
+    blocks:
+      - _type: material
+        _blockname: material_0
+        materials:
+          -
+            group: 0
+            model: linear
+            poisson_ratio: 0.3333
+            youngs_modulus: 300
+      - _type: interpolation
+        _blockname: interpolation_0
+        functions: 0
+        vector: 0
+      - _type: interpolation
+        _blockname: interpolation_1
+        functions: 0, -1
+        vector: 1
+      - _type: constraints
+        _blockname: constraints_0
+        functions: x < 1e-8
+        vector: 0
+      - _type: linearsolver
+        _blockname: linearsolver_0
+        operator:
+          reinforced_operator:
+            _blockname: reinforced_operator_0
+            fibres:
+              -
+                end:
+                  - 4.5
+                  - 0.501
+                radius: 0.1
+                start:
+                  - -0.5
+                  - 0.501
+                youngs_modulus: 1000000
+                prestress: 100000
+            vector: 0
+        vector: 0
+      - _type: visualization
+        _blockname: visualization_0
+        blocks:
+          - _type: vis_vector
+            _blockname: vis_vector_0
+            vector: 0
+          - _type: vis_vonmises
+            _blockname: vis_vonmises_0
+            continuous: true
+            vector: 0
+          # Default index: 1
+          - _type: vis_stress
+            _blockname: vis_stress_0
+            vector: 0
+          - _type: vis_stress
+            _blockname: vis_stress_1
+            vector: 0
+            index: 0
+        filename: fig23_3_prestress
+        instationary: false
+        intervals: 1
+        path: vtk

experiments/hansbo_reproduction/fig24_1.yml

@@ -53,7 +53,7 @@ solver:
             _blockname: vis_vonmises_0
             continuous: true
             vector: 0
-        filename: output
+        filename: fig24_1
         instationary: false
         intervals: 1
         path: vtk

experiments/hansbo_reproduction/fig24_2.yml

@@ -72,7 +72,7 @@ solver:
             _blockname: vis_vonmises_0
             continuous: true
             vector: 0
-        filename: output
+        filename: fig24_2
         instationary: false
         intervals: 1
         path: vtk