Add shear and compression wave velocities to Linear Elastic Material Properties for #692

include/materials/linear_elastic.h

@@ -77,6 +77,16 @@ class LinearElastic : public Material<Tdim> {
   double poisson_ratio_{std::numeric_limits<double>::max()};
   //! Bulk modulus
   double bulk_modulus_{std::numeric_limits<double>::max()};
+  //! S-Wave Velocity
+  double s_wave_velocity_{std::numeric_limits<double>::max()};
+  //! P-Wave Velocity
+  double p_wave_velocity_{std::numeric_limits<double>::max()};
+  //! Layer Thickness
+  double layer_thickness_{std::numeric_limits<double>::max()};
+  //! P Spring Coefficient
+  double p_spring_coeff_{std::numeric_limits<double>::max()};
+  //! S Spring Coefficient
+  double s_spring_coeff_{std::numeric_limits<double>::max()};
 };  // LinearElastic class
 }  // namespace mpm
 

include/materials/linear_elastic.tcc

@@ -9,9 +9,34 @@ mpm::LinearElastic<Tdim>::LinearElastic(unsigned id,
         material_properties.at("youngs_modulus").template get<double>();
     poisson_ratio_ =
         material_properties.at("poisson_ratio").template get<double>();
+
     // Calculate bulk modulus
     bulk_modulus_ = youngs_modulus_ / (3.0 * (1. - 2. * poisson_ratio_));
 
+    // Special material properties
+    if (material_properties.contains("earthquake")) {
+      bool earthquake =
+          material_properties.at("earthquake").template get<bool>();
+
+      if (earthquake) {
+        // Calculate constrained and shear modulus
+        double constrained_modulus =
+            youngs_modulus_ * (1. - poisson_ratio_) /
+            ((1. + poisson_ratio_) * (1. - 2. * poisson_ratio_));
+        double shear_modulus = youngs_modulus_ / (2.0 * (1. + poisson_ratio_));
+
+        // Calculate wave velocities
+        p_wave_velocity_ = sqrt(constrained_modulus / density_);
+        s_wave_velocity_ = sqrt(shear_modulus / density_);
+
+        // Calculate spring coefficients
+        layer_thickness_ =
+            material_properties.at("layer_thickness").template get<double>();
+        p_spring_coeff_ = constrained_modulus / layer_thickness_;
+        s_spring_coeff_ = shear_modulus / layer_thickness_;
+      }
+    }
+
     properties_ = material_properties;
     // Set elastic tensor
     this->compute_elastic_tensor();

tests/materials/linear_elastic_test.cc

@@ -147,6 +147,19 @@ TEST_CASE("LinearElastic is checked in 2D", "[material][linear_elastic][2D]") {
     REQUIRE(stress(4) == Approx(0.00000000000000e+00).epsilon(Tolerance));
     REQUIRE(stress(5) == Approx(0.00000000000000e+00).epsilon(Tolerance));
   }
+
+  SECTION("LinearElastic check properties earthquake") {
+    unsigned id = 0;
+    jmaterial["earthquake"] = true;
+    jmaterial["layer_thickness"] = 10000;
+
+    auto material =
+        Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
+            "LinearElastic2D", std::move(id), jmaterial);
+
+    REQUIRE(material->template property<double>("layer_thickness") ==
+            Approx(jmaterial["layer_thickness"]).epsilon(Tolerance));
+  }
 }
 
 //! Check linearelastic class in 3D