diff --git a/include/IO/property/reader.hpp b/include/IO/property/reader.hpp
index 6832e8c5..ac589f60 100644
--- a/include/IO/property/reader.hpp
+++ b/include/IO/property/reader.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
+#include "IO/reader.hpp"
 #include "compute/interface.hpp"
 #include "enumerations/interface.hpp"
-#include "IO/reader.hpp"
 
 namespace specfem {
 namespace IO {
diff --git a/include/IO/property/reader.tpp b/include/IO/property/reader.tpp
index 60253ebb..4fcd104a 100644
--- a/include/IO/property/reader.tpp
+++ b/include/IO/property/reader.tpp
@@ -1,18 +1,21 @@
 #pragma once
 
+#include "IO/property/reader.hpp"
 #include "compute/assembly/assembly.hpp"
 #include "enumerations/dimension.hpp"
 #include "enumerations/medium.hpp"
 #include "kokkos_abstractions.h"
 #include "point/properties.hpp"
-#include "IO/property/reader.hpp"
 #include <Kokkos_Core.hpp>
 
 template <typename InputLibrary>
-specfem::IO::property_reader<InputLibrary>::property_reader(const std::string input_folder): input_folder(input_folder) {}
+specfem::IO::property_reader<InputLibrary>::property_reader(
+    const std::string input_folder)
+    : input_folder(input_folder) {}
 
 template <typename InputLibrary>
-void specfem::IO::property_reader<InputLibrary>::read(specfem::compute::assembly &assembly) {
+void specfem::IO::property_reader<InputLibrary>::read(
+    specfem::compute::assembly &assembly) {
   auto &properties = assembly.properties;
 
   using DomainView =
@@ -21,32 +24,67 @@ void specfem::IO::property_reader<InputLibrary>::read(specfem::compute::assembly
   typename InputLibrary::File file(input_folder + "/Properties");
 
   {
-    typename InputLibrary::Group elastic = file.openGroup("/ElasticIsotropic");
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSVIsotropic");
+
+    elastic.openDataset("rho", properties.elastic_sv_isotropic.h_rho).read();
+    elastic.openDataset("mu", properties.elastic_sv_isotropic.h_mu).read();
+    elastic
+        .openDataset("lambdaplus2mu",
+                     properties.elastic_sv_isotropic.h_lambdaplus2mu)
+        .read();
+  }
+
+  {
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSHIsotropic");
 
-    elastic.openDataset("rho", properties.elastic_isotropic.h_rho).read();
-    elastic.openDataset("mu", properties.elastic_isotropic.h_mu).read();
-    elastic.openDataset("lambdaplus2mu", properties.elastic_isotropic.h_lambdaplus2mu).read();
+    elastic.openDataset("rho", properties.elastic_sh_isotropic.h_rho).read();
+    elastic.openDataset("mu", properties.elastic_sh_isotropic.h_mu).read();
+    elastic
+        .openDataset("lambdaplus2mu",
+                     properties.elastic_sh_isotropic.h_lambdaplus2mu)
+        .read();
   }
 
   {
-    typename InputLibrary::Group elastic = file.openGroup("/ElasticAnisotropic");
-
-    elastic.openDataset("rho", properties.elastic_anisotropic.h_rho).read();
-    elastic.openDataset("c11", properties.elastic_anisotropic.h_c11).read();
-    elastic.openDataset("c13", properties.elastic_anisotropic.h_c13).read();
-    elastic.openDataset("c15", properties.elastic_anisotropic.h_c15).read();
-    elastic.openDataset("c33", properties.elastic_anisotropic.h_c33).read();
-    elastic.openDataset("c35", properties.elastic_anisotropic.h_c35).read();
-    elastic.openDataset("c55", properties.elastic_anisotropic.h_c55).read();
-    elastic.openDataset("c12", properties.elastic_anisotropic.h_c12).read();
-    elastic.openDataset("c23", properties.elastic_anisotropic.h_c23).read();
-    elastic.openDataset("c25", properties.elastic_anisotropic.h_c25).read();
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSVAnisotropic");
+
+    elastic.openDataset("rho", properties.elastic_sv_anisotropic.h_rho).read();
+    elastic.openDataset("c11", properties.elastic_sv_anisotropic.h_c11).read();
+    elastic.openDataset("c13", properties.elastic_sv_anisotropic.h_c13).read();
+    elastic.openDataset("c15", properties.elastic_sv_anisotropic.h_c15).read();
+    elastic.openDataset("c33", properties.elastic_sv_anisotropic.h_c33).read();
+    elastic.openDataset("c35", properties.elastic_sv_anisotropic.h_c35).read();
+    elastic.openDataset("c55", properties.elastic_sv_anisotropic.h_c55).read();
+    elastic.openDataset("c12", properties.elastic_sv_anisotropic.h_c12).read();
+    elastic.openDataset("c23", properties.elastic_sv_anisotropic.h_c23).read();
+    elastic.openDataset("c25", properties.elastic_sv_anisotropic.h_c25).read();
+  }
+
+  {
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSHAnisotropic");
+
+    elastic.openDataset("rho", properties.elastic_sh_anisotropic.h_rho).read();
+    elastic.openDataset("c11", properties.elastic_sh_anisotropic.h_c11).read();
+    elastic.openDataset("c13", properties.elastic_sh_anisotropic.h_c13).read();
+    elastic.openDataset("c15", properties.elastic_sh_anisotropic.h_c15).read();
+    elastic.openDataset("c33", properties.elastic_sh_anisotropic.h_c33).read();
+    elastic.openDataset("c35", properties.elastic_sh_anisotropic.h_c35).read();
+    elastic.openDataset("c55", properties.elastic_sh_anisotropic.h_c55).read();
+    elastic.openDataset("c12", properties.elastic_sh_anisotropic.h_c12).read();
+    elastic.openDataset("c23", properties.elastic_sh_anisotropic.h_c23).read();
+    elastic.openDataset("c25", properties.elastic_sh_anisotropic.h_c25).read();
   }
 
   {
     typename InputLibrary::Group acoustic = file.openGroup("/Acoustic");
 
-    acoustic.openDataset("rho_inverse", properties.acoustic_isotropic.h_rho_inverse).read();
+    acoustic
+        .openDataset("rho_inverse", properties.acoustic_isotropic.h_rho_inverse)
+        .read();
     acoustic.openDataset("kappa", properties.acoustic_isotropic.h_kappa).read();
   }
 
diff --git a/include/IO/property/writer.hpp b/include/IO/property/writer.hpp
index 69fdecd1..8f0ff219 100644
--- a/include/IO/property/writer.hpp
+++ b/include/IO/property/writer.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
+#include "IO/writer.hpp"
 #include "compute/interface.hpp"
 #include "enumerations/interface.hpp"
-#include "IO/writer.hpp"
 
 namespace specfem {
 namespace IO {
@@ -37,5 +37,5 @@ template <typename OutputLibrary> class property_writer : public writer {
 private:
   std::string output_folder; ///< Path to output folder
 };
-} // namespace writer
+} // namespace IO
 } // namespace specfem
diff --git a/include/IO/property/writer.tpp b/include/IO/property/writer.tpp
index fdac6130..dede5fb6 100644
--- a/include/IO/property/writer.tpp
+++ b/include/IO/property/writer.tpp
@@ -29,22 +29,24 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
   const int ngllz = mesh.points.ngllz;
   const int ngllx = mesh.points.ngllx;
 
-  int n_elastic_isotropic;
-  int n_elastic_anisotropic;
+  int n_elastic_sv_isotropic;
+  int n_elastic_sh_isotropic;
+  int n_elastic_sv_anisotropic;
+  int n_elastic_sh_anisotropic;
   int n_acoustic;
 
   {
-    typename OutputLibrary::Group elastic = file.createGroup("/ElasticIsotropic");
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSVIsotropic");
 
     const auto element_indices = element_types.get_elements_on_host(
         specfem::element::medium_tag::elastic_sv,
         specfem::element::property_tag::isotropic);
-    n_elastic_isotropic = element_indices.size();
+        n_elastic_sv_isotropic = element_indices.size();
 
-    DomainView x("xcoordinates_elastic_isotropic", n_elastic_isotropic, ngllz, ngllx);
-    DomainView z("zcoordinates_elastic_isotropic", n_elastic_isotropic, ngllz, ngllx);
+    DomainView x("xcoordinates_elastic_isotropic", n_elastic_sv_isotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_isotropic", n_elastic_sv_isotropic, ngllz, ngllx);
 
-    for (int i = 0; i < n_elastic_isotropic; i++) {
+    for (int i = 0; i < n_elastic_sv_isotropic; i++) {
       const int ispec = element_indices(i);
       for (int iz = 0; iz < ngllz; iz++) {
         for (int ix = 0; ix < ngllx; ix++) {
@@ -57,23 +59,88 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
     elastic.createDataset("X", x).write();
     elastic.createDataset("Z", z).write();
 
-    elastic.createDataset("rho", properties.elastic_isotropic.h_rho).write();
-    elastic.createDataset("mu", properties.elastic_isotropic.h_mu).write();
-    elastic.createDataset("lambdaplus2mu", properties.elastic_isotropic.h_lambdaplus2mu).write();
+    elastic.createDataset("rho", properties.elastic_sv_isotropic.h_rho).write();
+    elastic.createDataset("mu", properties.elastic_sv_isotropic.h_mu).write();
+    elastic.createDataset("lambdaplus2mu", properties.elastic_sv_isotropic.h_lambdaplus2mu).write();
   }
 
   {
-    typename OutputLibrary::Group elastic = file.createGroup("/ElasticAnisotropic");
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSHIsotropic");
+
+    const auto element_indices = element_types.get_elements_on_host(
+        specfem::element::medium_tag::elastic_sh,
+        specfem::element::property_tag::isotropic);
+        n_elastic_sh_isotropic = element_indices.size();
+
+    DomainView x("xcoordinates_elastic_isotropic", n_elastic_sh_isotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_isotropic", n_elastic_sh_isotropic, ngllz, ngllx);
+
+    for (int i = 0; i < n_elastic_sh_isotropic; i++) {
+      const int ispec = element_indices(i);
+      for (int iz = 0; iz < ngllz; iz++) {
+        for (int ix = 0; ix < ngllx; ix++) {
+          x(i, iz, ix) = mesh.points.h_coord(0, ispec, iz, ix);
+          z(i, iz, ix) = mesh.points.h_coord(1, ispec, iz, ix);
+        }
+      }
+    }
+
+    elastic.createDataset("X", x).write();
+    elastic.createDataset("Z", z).write();
+
+    elastic.createDataset("rho", properties.elastic_sh_isotropic.h_rho).write();
+    elastic.createDataset("mu", properties.elastic_sh_isotropic.h_mu).write();
+    elastic.createDataset("lambdaplus2mu", properties.elastic_sh_isotropic.h_lambdaplus2mu).write();
+  }
+
+  {
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSVAnisotropic");
 
     const auto element_indices = element_types.get_elements_on_host(
         specfem::element::medium_tag::elastic_sv,
         specfem::element::property_tag::anisotropic);
-    n_elastic_anisotropic = element_indices.size();
+    n_elastic_sv_anisotropic = element_indices.size();
+
+    DomainView x("xcoordinates_elastic_anisotropic", n_elastic_sv_anisotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_anisotropic", n_elastic_sv_anisotropic, ngllz, ngllx);
+
+    for (int i = 0; i < n_elastic_sv_anisotropic; i++) {
+      const int ispec = element_indices(i);
+      for (int iz = 0; iz < ngllz; iz++) {
+        for (int ix = 0; ix < ngllx; ix++) {
+          x(i, iz, ix) = mesh.points.h_coord(0, ispec, iz, ix);
+          z(i, iz, ix) = mesh.points.h_coord(1, ispec, iz, ix);
+        }
+      }
+    }
+
+    elastic.createDataset("X", x).write();
+    elastic.createDataset("Z", z).write();
+
+    elastic.createDataset("rho", properties.elastic_sv_anisotropic.h_rho).write();
+    elastic.createDataset("c11", properties.elastic_sv_anisotropic.h_c11).write();
+    elastic.createDataset("c13", properties.elastic_sv_anisotropic.h_c13).write();
+    elastic.createDataset("c15", properties.elastic_sv_anisotropic.h_c15).write();
+    elastic.createDataset("c33", properties.elastic_sv_anisotropic.h_c33).write();
+    elastic.createDataset("c35", properties.elastic_sv_anisotropic.h_c35).write();
+    elastic.createDataset("c55", properties.elastic_sv_anisotropic.h_c55).write();
+    elastic.createDataset("c12", properties.elastic_sv_anisotropic.h_c12).write();
+    elastic.createDataset("c23", properties.elastic_sv_anisotropic.h_c23).write();
+    elastic.createDataset("c25", properties.elastic_sv_anisotropic.h_c25).write();
+  }
+
+  {
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSHAnisotropic");
+
+    const auto element_indices = element_types.get_elements_on_host(
+        specfem::element::medium_tag::elastic_sh,
+        specfem::element::property_tag::anisotropic);
+    n_elastic_sh_anisotropic = element_indices.size();
 
-    DomainView x("xcoordinates_elastic_anisotropic", n_elastic_anisotropic, ngllz, ngllx);
-    DomainView z("zcoordinates_elastic_anisotropic", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView x("xcoordinates_elastic_anisotropic", n_elastic_sh_anisotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_anisotropic", n_elastic_sh_anisotropic, ngllz, ngllx);
 
-    for (int i = 0; i < n_elastic_anisotropic; i++) {
+    for (int i = 0; i < n_elastic_sh_anisotropic; i++) {
       const int ispec = element_indices(i);
       for (int iz = 0; iz < ngllz; iz++) {
         for (int ix = 0; ix < ngllx; ix++) {
@@ -86,16 +153,16 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
     elastic.createDataset("X", x).write();
     elastic.createDataset("Z", z).write();
 
-    elastic.createDataset("rho", properties.elastic_anisotropic.h_rho).write();
-    elastic.createDataset("c11", properties.elastic_anisotropic.h_c11).write();
-    elastic.createDataset("c13", properties.elastic_anisotropic.h_c13).write();
-    elastic.createDataset("c15", properties.elastic_anisotropic.h_c15).write();
-    elastic.createDataset("c33", properties.elastic_anisotropic.h_c33).write();
-    elastic.createDataset("c35", properties.elastic_anisotropic.h_c35).write();
-    elastic.createDataset("c55", properties.elastic_anisotropic.h_c55).write();
-    elastic.createDataset("c12", properties.elastic_anisotropic.h_c12).write();
-    elastic.createDataset("c23", properties.elastic_anisotropic.h_c23).write();
-    elastic.createDataset("c25", properties.elastic_anisotropic.h_c25).write();
+    elastic.createDataset("rho", properties.elastic_sh_anisotropic.h_rho).write();
+    elastic.createDataset("c11", properties.elastic_sh_anisotropic.h_c11).write();
+    elastic.createDataset("c13", properties.elastic_sh_anisotropic.h_c13).write();
+    elastic.createDataset("c15", properties.elastic_sh_anisotropic.h_c15).write();
+    elastic.createDataset("c33", properties.elastic_sh_anisotropic.h_c33).write();
+    elastic.createDataset("c35", properties.elastic_sh_anisotropic.h_c35).write();
+    elastic.createDataset("c55", properties.elastic_sh_anisotropic.h_c55).write();
+    elastic.createDataset("c12", properties.elastic_sh_anisotropic.h_c12).write();
+    elastic.createDataset("c23", properties.elastic_sh_anisotropic.h_c23).write();
+    elastic.createDataset("c25", properties.elastic_sh_anisotropic.h_c25).write();
   }
 
   {
@@ -124,7 +191,7 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
     acoustic.createDataset("kappa", properties.acoustic_isotropic.h_kappa).write();
   }
 
-  assert(n_elastic_isotropic + n_elastic_anisotropic + n_acoustic == nspec);
+  assert(n_elastic_sv_isotropic + n_elastic_sv_anisotropic + n_elastic_sh_isotropic + n_elastic_sh_anisotropic + n_acoustic == nspec);
 
   std::cout << "Properties written to " << output_folder << "/Properties"
             << std::endl;
diff --git a/include/IO/seismogram/reader.hpp b/include/IO/seismogram/reader.hpp
index 5b9b3f04..8586c027 100644
--- a/include/IO/seismogram/reader.hpp
+++ b/include/IO/seismogram/reader.hpp
@@ -1,8 +1,8 @@
 #ifndef SPECFEM_READER_SEISMOGRAM_HPP
 #define SPECFEM_READER_SEISMOGRAM_HPP
 
-#include "enumerations/specfem_enums.hpp"
 #include "IO/reader.hpp"
+#include "enumerations/specfem_enums.hpp"
 
 namespace specfem {
 namespace forcing_function {
@@ -17,13 +17,13 @@ class seismogram_reader {
 public:
   seismogram_reader(){};
   seismogram_reader(const char *filename,
-             const specfem::enums::seismogram::format type,
-             specfem::kokkos::HostView2d<type_real> source_time_function)
+                    const specfem::enums::seismogram::format type,
+                    specfem::kokkos::HostView2d<type_real> source_time_function)
       : filename(filename), type(type),
         source_time_function(source_time_function) {}
   seismogram_reader(const std::string &filename,
-             const specfem::enums::seismogram::format type,
-             specfem::kokkos::HostView2d<type_real> source_time_function)
+                    const specfem::enums::seismogram::format type,
+                    specfem::kokkos::HostView2d<type_real> source_time_function)
       : filename(filename), type(type),
         source_time_function(source_time_function) {}
   void read();
diff --git a/include/IO/wavefield/writer.hpp b/include/IO/wavefield/writer.hpp
index a22fd719..46029f90 100644
--- a/include/IO/wavefield/writer.hpp
+++ b/include/IO/wavefield/writer.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
+#include "IO/writer.hpp"
 #include "compute/interface.hpp"
 #include "enumerations/interface.hpp"
-#include "IO/writer.hpp"
 
 namespace specfem {
 namespace IO {
diff --git a/include/IO/writer.hpp b/include/IO/writer.hpp
index 735d2ef6..f98a04de 100644
--- a/include/IO/writer.hpp
+++ b/include/IO/writer.hpp
@@ -23,5 +23,5 @@ class writer {
   virtual void write(specfem::compute::assembly &assembly){};
 };
 
-} // namespace writer
+} // namespace IO
 } // namespace specfem
diff --git a/include/compute/element_types/element_types.hpp b/include/compute/element_types/element_types.hpp
index bbd97c98..aced2d3e 100644
--- a/include/compute/element_types/element_types.hpp
+++ b/include/compute/element_types/element_types.hpp
@@ -112,6 +112,7 @@ struct element_types {
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(MEDIUM_TAG_VARIABLES,
                                  WHERE(DIMENSION_TAG_DIM2)
                                      WHERE(MEDIUM_TAG_ELASTIC_SV,
+                                           MEDIUM_TAG_ELASTIC_SH,
                                            MEDIUM_TAG_ACOUSTIC))
 
 #undef MEDIUM_TAG_VARIABLES
@@ -128,6 +129,7 @@ struct element_types {
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(MATERIAL_SYSTEMS_VARIABLE_NAMES,
                                       WHERE(DIMENSION_TAG_DIM2)
                                           WHERE(MEDIUM_TAG_ELASTIC_SV,
+                                                MEDIUM_TAG_ELASTIC_SH,
                                                 MEDIUM_TAG_ACOUSTIC)
                                               WHERE(PROPERTY_TAG_ISOTROPIC,
                                                     PROPERTY_TAG_ANISOTROPIC))
@@ -146,7 +148,8 @@ struct element_types {
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       ELEMENT_TYPES_VARIABLE_NAMES,
       WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
+          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH,
+                MEDIUM_TAG_ACOUSTIC)
               WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC)
                   WHERE(BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
                         BOUNDARY_TAG_STACEY,
diff --git a/include/compute/impl/value_containers.hpp b/include/compute/impl/value_containers.hpp
index 95d9c25c..806f52b6 100644
--- a/include/compute/impl/value_containers.hpp
+++ b/include/compute/impl/value_containers.hpp
@@ -30,11 +30,19 @@ struct value_containers {
 
   containers_type<specfem::element::medium_tag::elastic_sv,
                   specfem::element::property_tag::isotropic>
-      elastic_isotropic; ///< Elastic isotropic material values
+      elastic_sv_isotropic; ///< Elastic isotropic material values
 
   containers_type<specfem::element::medium_tag::elastic_sv,
                   specfem::element::property_tag::anisotropic>
-      elastic_anisotropic; ///< Elastic isotropic material values
+      elastic_sv_anisotropic; ///< Elastic isotropic material values
+
+  containers_type<specfem::element::medium_tag::elastic_sh,
+                  specfem::element::property_tag::isotropic>
+      elastic_sh_isotropic; ///< Elastic isotropic material values
+
+  containers_type<specfem::element::medium_tag::elastic_sh,
+                  specfem::element::property_tag::anisotropic>
+      elastic_sh_anisotropic; ///< Elastic isotropic material values
 
   containers_type<specfem::element::medium_tag::acoustic,
                   specfem::element::property_tag::isotropic>
@@ -62,17 +70,27 @@ struct value_containers {
       get_container() const {
     if constexpr ((MediumTag == specfem::element::medium_tag::elastic_sv) &&
                   (PropertyTag == specfem::element::property_tag::isotropic)) {
-      return elastic_isotropic;
+      return elastic_sv_isotropic;
     } else if constexpr ((MediumTag ==
                           specfem::element::medium_tag::elastic_sv) &&
                          (PropertyTag ==
                           specfem::element::property_tag::anisotropic)) {
-      return elastic_anisotropic;
+      return elastic_sv_anisotropic;
     } else if constexpr ((MediumTag ==
                           specfem::element::medium_tag::acoustic) &&
                          (PropertyTag ==
                           specfem::element::property_tag::isotropic)) {
       return acoustic_isotropic;
+    } else if constexpr ((MediumTag ==
+                          specfem::element::medium_tag::elastic_sh) &&
+                         (PropertyTag ==
+                          specfem::element::property_tag::isotropic)) {
+      return elastic_sh_isotropic;
+    } else if constexpr ((MediumTag ==
+                          specfem::element::medium_tag::elastic_sh) &&
+                         (PropertyTag ==
+                          specfem::element::property_tag::anisotropic)) {
+      return elastic_sh_anisotropic;
     } else {
       static_assert("Material type not implemented");
     }
@@ -83,14 +101,18 @@ struct value_containers {
    *
    */
   void copy_to_host() {
-    elastic_isotropic.copy_to_host();
-    elastic_anisotropic.copy_to_host();
+    elastic_sv_isotropic.copy_to_host();
+    elastic_sv_anisotropic.copy_to_host();
+    elastic_sh_isotropic.copy_to_host();
+    elastic_sh_anisotropic.copy_to_host();
     acoustic_isotropic.copy_to_host();
   }
 
   void copy_to_device() {
-    elastic_isotropic.copy_to_device();
-    elastic_anisotropic.copy_to_device();
+    elastic_sv_isotropic.copy_to_device();
+    elastic_sv_anisotropic.copy_to_device();
+    elastic_sh_isotropic.copy_to_device();
+    elastic_sh_anisotropic.copy_to_device();
     acoustic_isotropic.copy_to_device();
   }
 };
diff --git a/include/compute/properties/properties.hpp b/include/compute/properties/properties.hpp
index 2addf56a..3c81501d 100644
--- a/include/compute/properties/properties.hpp
+++ b/include/compute/properties/properties.hpp
@@ -23,7 +23,7 @@ namespace compute {
  *
  */
 struct properties
-    : public impl::value_containers<specfem::medium::material_properties> {
+    : public impl::value_containers<specfem::medium::properties_container> {
   /**
    * @name Constructors
    */
@@ -59,12 +59,12 @@ struct properties
    */
   void copy_to_host() {
     impl::value_containers<
-        specfem::medium::material_properties>::copy_to_host();
+        specfem::medium::properties_container>::copy_to_host();
   }
 
   void copy_to_device() {
     impl::value_containers<
-        specfem::medium::material_properties>::copy_to_device();
+        specfem::medium::properties_container>::copy_to_device();
   }
 };
 
diff --git a/include/compute/receivers/receivers.hpp b/include/compute/receivers/receivers.hpp
index a06636f9..2163be40 100644
--- a/include/compute/receivers/receivers.hpp
+++ b/include/compute/receivers/receivers.hpp
@@ -391,6 +391,7 @@ struct receivers : public impl::StationIterator,
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(RECEIVER_INDICES_VARIABLE_NAME,
                                       WHERE(DIMENSION_TAG_DIM2)
                                           WHERE(MEDIUM_TAG_ELASTIC_SV,
+                                                MEDIUM_TAG_ELASTIC_SH,
                                                 MEDIUM_TAG_ACOUSTIC)
                                               WHERE(PROPERTY_TAG_ISOTROPIC,
                                                     PROPERTY_TAG_ANISOTROPIC))
diff --git a/include/compute/sources/sources.hpp b/include/compute/sources/sources.hpp
index 72b6966d..ef474629 100644
--- a/include/compute/sources/sources.hpp
+++ b/include/compute/sources/sources.hpp
@@ -160,6 +160,7 @@ struct sources {
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(SOURCE_MEDIUM_DECLARATION,
                                  WHERE(DIMENSION_TAG_DIM2)
                                      WHERE(MEDIUM_TAG_ELASTIC_SV,
+                                           MEDIUM_TAG_ELASTIC_SH,
                                            MEDIUM_TAG_ACOUSTIC))
 
 #undef SOURCE_MEDIUM_DECLARATION
@@ -208,7 +209,8 @@ struct sources {
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       SOURCE_INDICES_VARIABLES_NAME,
       WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
+          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH,
+                MEDIUM_TAG_ACOUSTIC)
               WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC)
                   WHERE(BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
                         BOUNDARY_TAG_STACEY,
@@ -310,8 +312,8 @@ load_on_device(const IteratorIndexType iterator_index,
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
       SOURCE_MEDIUM_LOAD_ON_DEVICE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC))
 
 #undef SOURCE_MEDIUM_LOAD_ON_DEVICE
 
@@ -388,8 +390,8 @@ void load_on_host(const IteratorIndexType iterator_index,
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
       SOURCE_MEDIUM_LOAD_ON_HOST,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC))
 
 #undef SOURCE_MEDIUM_LOAD_ON_HOST
 
diff --git a/include/medium/dim2/acoustic/isotropic/properties_container.hpp b/include/medium/dim2/acoustic/isotropic/properties_container.hpp
index 9c76c2b6..a0b5a333 100644
--- a/include/medium/dim2/acoustic/isotropic/properties_container.hpp
+++ b/include/medium/dim2/acoustic/isotropic/properties_container.hpp
@@ -12,7 +12,7 @@ struct properties_container<specfem::element::medium_tag::acoustic,
                             specfem::element::property_tag::isotropic> {
 
   constexpr static auto dimension = specfem::dimension::type::dim2;
-  constexpr static auto value_type = specfem::element::medium_tag::acoustic;
+  constexpr static auto medium_tag = specfem::element::medium_tag::acoustic;
   constexpr static auto property_type =
       specfem::element::property_tag::isotropic;
 
@@ -37,6 +37,17 @@ struct properties_container<specfem::element::medium_tag::acoustic,
         kappa("specfem::compute::properties::kappa", nspec, ngllz, ngllx),
         h_kappa(Kokkos::create_mirror_view(kappa)) {}
 
+  properties_container(
+      const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
+      const int ngllz, const int ngllx,
+      const specfem::mesh::materials &materials, const bool has_gll_model,
+      const specfem::kokkos::HostView1d<int> property_index_mapping)
+      : properties_container(elements.extent(0), ngllz, ngllx) {
+
+    impl::constructor(elements, ngllz, ngllx, materials, has_gll_model,
+                      property_index_mapping, *this);
+  }
+
   template <
       typename PointProperties,
       typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
@@ -46,7 +57,7 @@ struct properties_container<specfem::element::medium_tag::acoustic,
 
     static_assert(PointProperties::dimension == dimension,
                   "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
+    static_assert(PointProperties::medium_tag == medium_tag,
                   "Medium tag mismatch");
     static_assert(PointProperties::property_tag == property_type,
                   "Property tag mismatch");
@@ -71,7 +82,7 @@ struct properties_container<specfem::element::medium_tag::acoustic,
 
     static_assert(PointProperties::dimension == dimension,
                   "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
+    static_assert(PointProperties::medium_tag == medium_tag,
                   "Medium tag mismatch");
     static_assert(PointProperties::property_tag == property_type,
                   "Property tag mismatch");
@@ -105,7 +116,7 @@ struct properties_container<specfem::element::medium_tag::acoustic,
 
     static_assert(PointProperties::dimension == dimension,
                   "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
+    static_assert(PointProperties::medium_tag == medium_tag,
                   "Medium tag mismatch");
     static_assert(PointProperties::property_tag == property_type,
                   "Property tag mismatch");
@@ -130,7 +141,7 @@ struct properties_container<specfem::element::medium_tag::acoustic,
 
     static_assert(PointProperties::dimension == dimension,
                   "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
+    static_assert(PointProperties::medium_tag == medium_tag,
                   "Medium tag mismatch");
     static_assert(PointProperties::property_tag == property_type,
                   "Property tag mismatch");
@@ -173,7 +184,7 @@ struct properties_container<specfem::element::medium_tag::acoustic,
 
     static_assert(PointProperties::dimension == dimension,
                   "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
+    static_assert(PointProperties::medium_tag == medium_tag,
                   "Medium tag mismatch");
     static_assert(PointProperties::property_tag == property_type,
                   "Property tag mismatch");
@@ -194,7 +205,7 @@ struct properties_container<specfem::element::medium_tag::acoustic,
 
     static_assert(PointProperties::dimension == dimension,
                   "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
+    static_assert(PointProperties::medium_tag == medium_tag,
                   "Medium tag mismatch");
     static_assert(PointProperties::property_tag == property_type,
                   "Property tag mismatch");
diff --git a/include/medium/dim2/elastic/anisotropic/kernels_container.hpp b/include/medium/dim2/elastic/anisotropic/kernels_container.hpp
index 40349e89..24a7ec16 100644
--- a/include/medium/dim2/elastic/anisotropic/kernels_container.hpp
+++ b/include/medium/dim2/elastic/anisotropic/kernels_container.hpp
@@ -9,10 +9,11 @@ namespace specfem {
 namespace medium {
 
 template <>
-class kernels_container<specfem::element::medium_tag::elastic_sv,
+class kernels_container<specfem::element::medium_tag::elastic,
                         specfem::element::property_tag::anisotropic> {
 public:
-  constexpr static auto value_type = specfem::element::medium_tag::elastic_sv;
+  constexpr static auto dimension = specfem::dimension::type::dim2;
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic;
   constexpr static auto property_type =
       specfem::element::property_tag::anisotropic;
   int nspec;
@@ -59,15 +60,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     initialize();
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void load_device_kernels(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -82,19 +78,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     kernels.c55 = c55(ispec, iz, ix);
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void load_device_kernels(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      PointKernelType &kernels) const {
+      const specfem::point::simd_index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
-
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = specfem::datatype::simd<type_real, true>;
+    using mask_type = simd::mask_type;
+    using tag_type = simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -118,15 +109,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_from(&c55(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
-  void
-  load_host_kernels(specfem::point::index<PointKernelType::dimension> &index,
-                    PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+  void load_host_kernels(
+      specfem::point::index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -141,19 +127,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     kernels.c55 = h_c55(ispec, iz, ix);
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   void load_host_kernels(
-      specfem::point::simd_index<PointKernelType::dimension> &index,
-      PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      specfem::point::simd_index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -177,15 +158,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_from(&h_c55(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void update_kernels_on_device(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -200,19 +176,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     c55(ispec, iz, ix) = kernels.c55;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void update_kernels_on_device(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -236,15 +207,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_to(&c55(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   void update_kernels_on_host(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -259,19 +225,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     h_c55(ispec, iz, ix) = kernels.c55;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   void update_kernels_on_host(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -294,15 +255,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_to(&h_c55(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void add_kernels_on_device(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -317,19 +273,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     c55(ispec, iz, ix) += kernels.c55;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void add_kernels_on_device(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = typename specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -337,7 +288,7 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     const int iz = index.iz;
     const int ix = index.ix;
 
-    simd_type lhs;
+    simd::datatype lhs;
 
     Kokkos::Experimental::where(mask, lhs).copy_from(&rho(ispec, iz, ix),
                                                      tag_type());
@@ -382,15 +333,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
                                                    tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   void add_kernels_on_host(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -405,19 +351,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     h_c55(ispec, iz, ix) += kernels.c55;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   void add_kernels_on_host(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -425,7 +366,7 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     const int iz = index.iz;
     const int ix = index.ix;
 
-    simd_type lhs;
+    simd::datatype lhs;
 
     Kokkos::Experimental::where(mask, lhs).copy_from(&h_rho(ispec, iz, ix),
                                                      tag_type());
@@ -507,5 +448,38 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
   }
 };
 
+template <>
+class kernels_container<specfem::element::medium_tag::elastic_sv,
+                        specfem::element::property_tag::anisotropic>
+    : public kernels_container<specfem::element::medium_tag::elastic,
+                               specfem::element::property_tag::anisotropic> {
+
+public:
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+
+  kernels_container() = default;
+
+  kernels_container(const int nspec, const int ngllz, const int ngllx)
+      : kernels_container<specfem::element::medium_tag::elastic,
+                          specfem::element::property_tag::anisotropic>(
+            nspec, ngllz, ngllx) {}
+};
+
+template <>
+class kernels_container<specfem::element::medium_tag::elastic_sh,
+                        specfem::element::property_tag::anisotropic>
+    : public kernels_container<specfem::element::medium_tag::elastic,
+                               specfem::element::property_tag::anisotropic> {
+
+public:
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh;
+
+  kernels_container() = default;
+
+  kernels_container(const int nspec, const int ngllz, const int ngllx)
+      : kernels_container<specfem::element::medium_tag::elastic,
+                          specfem::element::property_tag::anisotropic>(
+            nspec, ngllz, ngllx) {}
+};
 } // namespace medium
 } // namespace specfem
diff --git a/include/medium/dim2/elastic/anisotropic/properties_container.hpp b/include/medium/dim2/elastic/anisotropic/properties_container.hpp
index 3f3531ff..cac5e444 100644
--- a/include/medium/dim2/elastic/anisotropic/properties_container.hpp
+++ b/include/medium/dim2/elastic/anisotropic/properties_container.hpp
@@ -8,11 +8,11 @@ namespace specfem {
 namespace medium {
 
 template <>
-struct properties_container<specfem::element::medium_tag::elastic_sv,
+struct properties_container<specfem::element::medium_tag::elastic,
                             specfem::element::property_tag::anisotropic> {
 
   constexpr static auto dimension = specfem::dimension::type::dim2;
-  constexpr static auto value_type = specfem::element::medium_tag::elastic_sv;
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic;
   constexpr static auto property_type =
       specfem::element::property_tag::anisotropic;
 
@@ -69,19 +69,10 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
         c25("specfem::compute::properties::c25", nspec, ngllz, ngllx),
         h_c25(Kokkos::create_mirror_view(c25)) {}
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
-  KOKKOS_FORCEINLINE_FUNCTION void
-  load_device_properties(const specfem::point::index<dimension> &index,
-                         PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  KOKKOS_FORCEINLINE_FUNCTION void load_device_properties(
+      const specfem::point::index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, false>
+          &property) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -102,21 +93,12 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     property.rho_vs = sqrt(property.rho * property.c55);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
-  KOKKOS_FORCEINLINE_FUNCTION void
-  load_device_properties(const specfem::point::simd_index<dimension> &index,
-                         PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
-
-    using simd = typename PointProperties::simd;
+  KOKKOS_FORCEINLINE_FUNCTION void load_device_properties(
+      const specfem::point::simd_index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, true>
+          &property) const {
+
+    using simd = typename specfem::datatype::simd<type_real, true>;
     using mask_type = typename simd::mask_type;
     using tag_type = typename simd::tag_type;
 
@@ -151,19 +133,10 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     property.rho_vs = Kokkos::sqrt(property.rho * property.c55);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
-  inline void
-  load_host_properties(const specfem::point::index<dimension> &index,
-                       PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void load_host_properties(
+      const specfem::point::index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, false>
+          &property) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -184,21 +157,12 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     property.rho_vs = sqrt(property.rho * property.c55);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
-  inline void
-  load_host_properties(const specfem::point::simd_index<dimension> &index,
-                       PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void load_host_properties(
+      const specfem::point::simd_index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, true>
+          &property) const {
 
-    using simd = typename PointProperties::simd;
+    using simd = typename specfem::datatype::simd<type_real, true>;
     using mask_type = typename simd::mask_type;
     using tag_type = typename simd::tag_type;
 
@@ -257,18 +221,10 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     Kokkos::deep_copy(h_c25, c25);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
-  inline void assign(const specfem::point::index<dimension> &index,
-                     const PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    // static_assert(PointProperties::medium_tag == value_type,
-    //               "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void
+  assign(const specfem::point::index<dimension> &index,
+         const specfem::point::properties<dimension, medium_tag, property_type,
+                                          false> &property) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -286,20 +242,12 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     h_c25(ispec, iz, ix) = property.c25;
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
-  inline void assign(const specfem::point::simd_index<dimension> &index,
-                     const PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    // static_assert(PointProperties::medium_tag == value_type,
-    //               "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void
+  assign(const specfem::point::simd_index<dimension> &index,
+         const specfem::point::properties<dimension, medium_tag, property_type,
+                                          true> &property) const {
 
-    using simd = typename PointProperties::simd;
+    using simd = typename specfem::datatype::simd<type_real, true>;
     using mask_type = typename simd::mask_type;
     using tag_type = typename simd::tag_type;
 
@@ -332,5 +280,58 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
   }
 };
 
+template <>
+struct properties_container<specfem::element::medium_tag::elastic_sv,
+                            specfem::element::property_tag::anisotropic>
+    : public properties_container<specfem::element::medium_tag::elastic,
+                                  specfem::element::property_tag::anisotropic> {
+
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+
+  properties_container() = default;
+
+  properties_container(const int nspec, const int ngllz, const int ngllx)
+      : properties_container<specfem::element::medium_tag::elastic,
+                             specfem::element::property_tag::anisotropic>(
+            nspec, ngllz, ngllx){};
+
+  properties_container(
+      const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
+      const int ngllz, const int ngllx,
+      const specfem::mesh::materials &materials, const bool has_gll_model,
+      const specfem::kokkos::HostView1d<int> property_index_mapping)
+      : properties_container(elements.extent(0), ngllz, ngllx) {
+
+    impl::constructor(elements, ngllz, ngllx, materials, has_gll_model,
+                      property_index_mapping, *this);
+  }
+};
+
+template <>
+struct properties_container<specfem::element::medium_tag::elastic_sh,
+                            specfem::element::property_tag::anisotropic>
+    : public properties_container<specfem::element::medium_tag::elastic,
+                                  specfem::element::property_tag::anisotropic> {
+
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh;
+  properties_container() = default;
+
+  properties_container(const int nspec, const int ngllz, const int ngllx)
+      : properties_container<specfem::element::medium_tag::elastic,
+                             specfem::element::property_tag::anisotropic>(
+            nspec, ngllz, ngllx){};
+
+  properties_container(
+      const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
+      const int ngllz, const int ngllx,
+      const specfem::mesh::materials &materials, const bool has_gll_model,
+      const specfem::kokkos::HostView1d<int> property_index_mapping)
+      : properties_container(elements.extent(0), ngllz, ngllx) {
+
+    impl::constructor(elements, ngllz, ngllx, materials, has_gll_model,
+                      property_index_mapping, *this);
+  }
+};
+
 } // namespace medium
 } // namespace specfem
diff --git a/include/medium/dim2/elastic/isotropic/kernels_container.hpp b/include/medium/dim2/elastic/isotropic/kernels_container.hpp
index 09a96f42..78a6c8d7 100644
--- a/include/medium/dim2/elastic/isotropic/kernels_container.hpp
+++ b/include/medium/dim2/elastic/isotropic/kernels_container.hpp
@@ -10,10 +10,11 @@ namespace specfem {
 namespace medium {
 
 template <>
-class kernels_container<specfem::element::medium_tag::elastic_sv,
+class kernels_container<specfem::element::medium_tag::elastic,
                         specfem::element::property_tag::isotropic> {
 public:
-  constexpr static auto value_type = specfem::element::medium_tag::elastic_sv;
+  constexpr static auto dimension = specfem::dimension::type::dim2;
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic;
   constexpr static auto property_type =
       specfem::element::property_tag::isotropic;
   int nspec;
@@ -56,15 +57,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     initialize();
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void load_device_kernels(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -78,19 +74,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     kernels.beta = beta(ispec, iz, ix);
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void load_device_kernels(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      PointKernelType &kernels) const {
+      const specfem::point::simd_index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
-
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = typename specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -112,15 +103,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_from(&beta(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
-  void
-  load_host_kernels(specfem::point::index<PointKernelType::dimension> &index,
-                    PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+  void load_host_kernels(
+      specfem::point::index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -134,19 +120,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     kernels.beta = h_beta(ispec, iz, ix);
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   void load_host_kernels(
-      specfem::point::simd_index<PointKernelType::dimension> &index,
-      PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      specfem::point::simd_index<dimension> &index,
+      specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = typename specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -168,15 +149,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_from(&h_beta(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void update_kernels_on_device(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -190,19 +166,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     beta(ispec, iz, ix) = kernels.beta;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void update_kernels_on_device(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = typename specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -224,15 +195,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_to(&beta(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   void update_kernels_on_host(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -246,19 +212,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     h_beta(ispec, iz, ix) = kernels.beta;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   void update_kernels_on_host(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = typename specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -280,15 +241,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
         .copy_to(&h_beta(ispec, iz, ix), tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void add_kernels_on_device(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -302,19 +258,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     beta(ispec, iz, ix) += kernels.beta;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   KOKKOS_INLINE_FUNCTION void add_kernels_on_device(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = typename specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -322,7 +273,7 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     const int iz = index.iz;
     const int ix = index.ix;
 
-    simd_type lhs;
+    simd::datatype lhs;
 
     Kokkos::Experimental::where(mask, lhs).copy_from(&rho(ispec, iz, ix),
                                                      tag_type());
@@ -361,15 +312,10 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
                                                    tag_type());
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<!PointKernelType::simd::using_simd, int> = 0>
   void add_kernels_on_host(
-      const specfem::point::index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, false>
+          &kernels) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -383,19 +329,14 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     h_beta(ispec, iz, ix) += kernels.beta;
   }
 
-  template <
-      typename PointKernelType,
-      typename std::enable_if_t<PointKernelType::simd::using_simd, int> = 0>
   void add_kernels_on_host(
-      const specfem::point::simd_index<PointKernelType::dimension> &index,
-      const PointKernelType &kernels) const {
-
-    static_assert(PointKernelType::medium_tag == value_type);
-    static_assert(PointKernelType::property_tag == property_type);
+      const specfem::point::simd_index<dimension> &index,
+      const specfem::point::kernels<dimension, medium_tag, property_type, true>
+          &kernels) const {
 
-    using simd_type = typename PointKernelType::simd::datatype;
-    using mask_type = typename PointKernelType::simd::mask_type;
-    using tag_type = typename PointKernelType::simd::tag_type;
+    using simd = typename specfem::datatype::simd<type_real, true>;
+    using mask_type = typename simd::mask_type;
+    using tag_type = typename simd::tag_type;
 
     mask_type mask([&](std::size_t lane) { return index.mask(lane); });
 
@@ -403,7 +344,7 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
     const int iz = index.iz;
     const int ix = index.ix;
 
-    simd_type lhs;
+    simd::datatype lhs;
 
     Kokkos::Experimental::where(mask, lhs).copy_from(&h_rho(ispec, iz, ix),
                                                      tag_type());
@@ -476,5 +417,38 @@ class kernels_container<specfem::element::medium_tag::elastic_sv,
   }
 };
 
+template <>
+class kernels_container<specfem::element::medium_tag::elastic_sv,
+                        specfem::element::property_tag::isotropic>
+    : public kernels_container<specfem::element::medium_tag::elastic,
+                               specfem::element::property_tag::isotropic> {
+public:
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+
+  kernels_container() = default;
+
+  kernels_container(const int nspec, const int ngllz, const int ngllx)
+      : kernels_container<specfem::element::medium_tag::elastic,
+                          specfem::element::property_tag::isotropic>(
+            nspec, ngllz, ngllx) {}
+};
+
+template <>
+class kernels_container<specfem::element::medium_tag::elastic_sh,
+                        specfem::element::property_tag::isotropic>
+    : public kernels_container<specfem::element::medium_tag::elastic,
+                               specfem::element::property_tag::isotropic> {
+
+public:
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh;
+
+  kernels_container() = default;
+
+  kernels_container(const int nspec, const int ngllz, const int ngllx)
+      : kernels_container<specfem::element::medium_tag::elastic,
+                          specfem::element::property_tag::isotropic>(
+            nspec, ngllz, ngllx) {}
+};
+
 } // namespace medium
 } // namespace specfem
diff --git a/include/medium/dim2/elastic/isotropic/properties_container.hpp b/include/medium/dim2/elastic/isotropic/properties_container.hpp
index 4537f5ab..d7f80772 100644
--- a/include/medium/dim2/elastic/isotropic/properties_container.hpp
+++ b/include/medium/dim2/elastic/isotropic/properties_container.hpp
@@ -8,11 +8,11 @@ namespace specfem {
 namespace medium {
 
 template <>
-struct properties_container<specfem::element::medium_tag::elastic_sv,
+struct properties_container<specfem::element::medium_tag::elastic,
                             specfem::element::property_tag::isotropic> {
 
   constexpr static auto dimension = specfem::dimension::type::dim2;
-  constexpr static auto value_type = specfem::element::medium_tag::elastic_sv;
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic;
   constexpr static auto property_type =
       specfem::element::property_tag::isotropic;
 
@@ -41,19 +41,10 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
                       ngllz, ngllx),
         h_lambdaplus2mu(Kokkos::create_mirror_view(lambdaplus2mu)) {}
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
-  KOKKOS_FORCEINLINE_FUNCTION void
-  load_device_properties(const specfem::point::index<dimension> &index,
-                         PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  KOKKOS_FORCEINLINE_FUNCTION void load_device_properties(
+      const specfem::point::index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, false>
+          &property) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -67,21 +58,12 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     property.rho_vs = sqrt(property.rho * property.mu);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
-  KOKKOS_FORCEINLINE_FUNCTION void
-  load_device_properties(const specfem::point::simd_index<dimension> &index,
-                         PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
-
-    using simd = typename PointProperties::simd;
+  KOKKOS_FORCEINLINE_FUNCTION void load_device_properties(
+      const specfem::point::simd_index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, true>
+          &property) const {
+
+    using simd = specfem::datatype::simd<type_real, true>;
     using mask_type = typename simd::mask_type;
     using tag_type = typename simd::tag_type;
 
@@ -103,19 +85,10 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     property.rho_vs = Kokkos::sqrt(property.rho * property.mu);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
-  inline void
-  load_host_properties(const specfem::point::index<dimension> &index,
-                       PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void load_host_properties(
+      const specfem::point::index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, false>
+          &property) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -129,21 +102,12 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     property.rho_vs = sqrt(property.rho * property.mu);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
-  inline void
-  load_host_properties(const specfem::point::simd_index<dimension> &index,
-                       PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    static_assert(PointProperties::medium_tag == value_type,
-                  "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void load_host_properties(
+      const specfem::point::simd_index<dimension> &index,
+      specfem::point::properties<dimension, medium_tag, property_type, true>
+          &property) const {
 
-    using simd = typename PointProperties::simd;
+    using simd = specfem::datatype::simd<type_real, true>;
     using mask_type = typename simd::mask_type;
     using tag_type = typename simd::tag_type;
 
@@ -177,18 +141,10 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     Kokkos::deep_copy(h_lambdaplus2mu, lambdaplus2mu);
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
-  inline void assign(const specfem::point::index<dimension> &index,
-                     const PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    // static_assert(PointProperties::medium_tag == value_type,
-    //               "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void
+  assign(const specfem::point::index<dimension> &index,
+         const specfem::point::properties<dimension, medium_tag, property_type,
+                                          false> &property) const {
 
     const int ispec = index.ispec;
     const int iz = index.iz;
@@ -199,20 +155,12 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
     h_lambdaplus2mu(ispec, iz, ix) = property.lambdaplus2mu;
   }
 
-  template <
-      typename PointProperties,
-      typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
-  inline void assign(const specfem::point::simd_index<dimension> &index,
-                     const PointProperties &property) const {
-
-    static_assert(PointProperties::dimension == dimension,
-                  "Dimension mismatch");
-    // static_assert(PointProperties::medium_tag == value_type,
-    //               "Medium tag mismatch");
-    static_assert(PointProperties::property_tag == property_type,
-                  "Property tag mismatch");
+  inline void
+  assign(const specfem::point::simd_index<dimension> &index,
+         const specfem::point::properties<dimension, medium_tag, property_type,
+                                          true> &property) const {
 
-    using simd = typename PointProperties::simd;
+    using simd = specfem::datatype::simd<type_real, true>;
     using mask_type = typename simd::mask_type;
     using tag_type = typename simd::tag_type;
 
@@ -231,5 +179,58 @@ struct properties_container<specfem::element::medium_tag::elastic_sv,
   }
 };
 
+template <>
+struct properties_container<specfem::element::medium_tag::elastic_sv,
+                            specfem::element::property_tag::isotropic>
+    : public properties_container<specfem::element::medium_tag::elastic,
+                                  specfem::element::property_tag::isotropic> {
+
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+
+  properties_container() = default;
+
+  properties_container(const int nspec, const int ngllz, const int ngllx)
+      : properties_container<specfem::element::medium_tag::elastic,
+                             specfem::element::property_tag::isotropic>(
+            nspec, ngllz, ngllx){};
+
+  properties_container(
+      const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
+      const int ngllz, const int ngllx,
+      const specfem::mesh::materials &materials, const bool has_gll_model,
+      const specfem::kokkos::HostView1d<int> property_index_mapping)
+      : properties_container(elements.extent(0), ngllz, ngllx) {
+
+    impl::constructor(elements, ngllz, ngllx, materials, has_gll_model,
+                      property_index_mapping, *this);
+  }
+};
+
+template <>
+struct properties_container<specfem::element::medium_tag::elastic_sh,
+                            specfem::element::property_tag::isotropic>
+    : public properties_container<specfem::element::medium_tag::elastic,
+                                  specfem::element::property_tag::isotropic> {
+
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh;
+  properties_container() = default;
+
+  properties_container(const int nspec, const int ngllz, const int ngllx)
+      : properties_container<specfem::element::medium_tag::elastic,
+                             specfem::element::property_tag::isotropic>(
+            nspec, ngllz, ngllx){};
+
+  properties_container(
+      const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
+      const int ngllz, const int ngllx,
+      const specfem::mesh::materials &materials, const bool has_gll_model,
+      const specfem::kokkos::HostView1d<int> property_index_mapping)
+      : properties_container(elements.extent(0), ngllz, ngllx) {
+
+    impl::constructor(elements, ngllz, ngllx, materials, has_gll_model,
+                      property_index_mapping, *this);
+  }
+};
+
 } // namespace medium
 } // namespace specfem
diff --git a/include/medium/material_properties.hpp b/include/medium/material_properties.hpp
index 5c146cb1..4ed38985 100644
--- a/include/medium/material_properties.hpp
+++ b/include/medium/material_properties.hpp
@@ -9,54 +9,54 @@
 namespace specfem {
 namespace medium {
 
-template <specfem::element::medium_tag type,
-          specfem::element::property_tag property>
-struct material_properties
-    : public specfem::medium::properties_container<type, property> {
-  constexpr static auto value_type = type;
-  constexpr static auto property_type = property;
-  constexpr static auto dimension = specfem::dimension::type::dim2;
-
-  material_properties() = default;
-
-  material_properties(
-      const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
-      const int ngllz, const int ngllx,
-      const specfem::mesh::materials &materials, const bool has_gll_model,
-      const specfem::kokkos::HostView1d<int> property_index_mapping)
-      : specfem::medium::properties_container<type, property>(
-            elements.extent(0), ngllz, ngllx) {
-
-    const int nelement = elements.extent(0);
-    int count = 0;
-    for (int i = 0; i < nelement; ++i) {
-      const int ispec = elements(i);
-      property_index_mapping(ispec) = count;
-      if (!has_gll_model) {
-        for (int iz = 0; iz < ngllz; ++iz) {
-          for (int ix = 0; ix < ngllx; ++ix) {
-            // Get the material at index from mesh::materials
-            auto material =
-                std::get<specfem::medium::material<type, property> >(
-                    materials[ispec]);
-
-            // Assign the material property to the property container
-            auto point_property = material.get_properties();
-            this->assign(specfem::point::index<dimension>(count, iz, ix),
-                         point_property);
-          }
-        }
-      }
-      count++;
-    }
-
-    if (!has_gll_model) {
-      this->copy_to_device();
-    }
-
-    return;
-  }
-};
+// template <specfem::element::medium_tag type,
+//           specfem::element::property_tag property>
+// struct material_properties
+//     : public specfem::medium::properties_container<type, property> {
+//   constexpr static auto value_type = type;
+//   constexpr static auto property_type = property;
+//   constexpr static auto dimension = specfem::dimension::type::dim2;
+
+//   material_properties() = default;
+
+//   material_properties(
+//       const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
+//       const int ngllz, const int ngllx,
+//       const specfem::mesh::materials &materials, const bool has_gll_model,
+//       const specfem::kokkos::HostView1d<int> property_index_mapping)
+//       : specfem::medium::properties_container<type, property>(
+//             elements.extent(0), ngllz, ngllx) {
+
+//     const int nelement = elements.extent(0);
+//     int count = 0;
+//     for (int i = 0; i < nelement; ++i) {
+//       const int ispec = elements(i);
+//       property_index_mapping(ispec) = count;
+//       if (!has_gll_model) {
+//         for (int iz = 0; iz < ngllz; ++iz) {
+//           for (int ix = 0; ix < ngllx; ++ix) {
+//             // Get the material at index from mesh::materials
+//             auto material =
+//                 std::get<specfem::medium::material<type, property> >(
+//                     materials[ispec]);
+
+//             // Assign the material property to the property container
+//             auto point_property = material.get_properties();
+//             this->assign(specfem::point::index<dimension>(count, iz, ix),
+//                          point_property);
+//           }
+//         }
+//       }
+//       count++;
+//     }
+
+//     if (!has_gll_model) {
+//       this->copy_to_device();
+//     }
+
+//     return;
+//   }
+// };
 
 } // namespace medium
 } // namespace specfem
diff --git a/include/medium/properties_container.hpp b/include/medium/properties_container.hpp
index 2b3c9d89..82d390c7 100644
--- a/include/medium/properties_container.hpp
+++ b/include/medium/properties_container.hpp
@@ -5,6 +5,50 @@
 namespace specfem {
 namespace medium {
 
+namespace impl {
+template <typename PropertiesContainer>
+void constructor(
+    const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
+    const int ngllz, const int ngllx, const specfem::mesh::materials &materials,
+    const bool has_gll_model,
+    const specfem::kokkos::HostView1d<int> property_index_mapping,
+    PropertiesContainer &properties) {
+
+  constexpr auto dimension = PropertiesContainer::dimension;
+  constexpr auto type = PropertiesContainer::medium_tag;
+  constexpr auto property = PropertiesContainer::property_type;
+
+  const int nelement = elements.extent(0);
+  int count = 0;
+  for (int i = 0; i < nelement; ++i) {
+    const int ispec = elements(i);
+    property_index_mapping(ispec) = count;
+    if (!has_gll_model) {
+      for (int iz = 0; iz < ngllz; ++iz) {
+        for (int ix = 0; ix < ngllx; ++ix) {
+          // Get the material at index from mesh::materials
+          auto material = std::get<specfem::medium::material<type, property> >(
+              materials[ispec]);
+
+          // Assign the material property to the property container
+          auto point_property = material.get_properties();
+          properties.assign(specfem::point::index<dimension>(count, iz, ix),
+                            point_property);
+        }
+      }
+    }
+    count++;
+  }
+
+  if (!has_gll_model) {
+    properties.copy_to_device();
+  }
+
+  return;
+}
+
+} // namespace impl
+
 template <specfem::element::medium_tag type,
           specfem::element::property_tag property>
 struct properties_container {
diff --git a/include/parameter_parser/interface.hpp b/include/parameter_parser/interface.hpp
index 5fe0b32e..269e91e4 100644
--- a/include/parameter_parser/interface.hpp
+++ b/include/parameter_parser/interface.hpp
@@ -7,8 +7,8 @@
 #include "run_setup.hpp"
 #include "setup.hpp"
 // #include "solver/interface.hpp"
-#include "time_scheme/interface.hpp"
 #include "IO/seismogram/writer.hpp"
 #include "IO/wavefield/writer.hpp"
+#include "time_scheme/interface.hpp"
 
 #endif
diff --git a/include/parameter_parser/writer/kernel.hpp b/include/parameter_parser/writer/kernel.hpp
index 9c84bca9..7aab6685 100644
--- a/include/parameter_parser/writer/kernel.hpp
+++ b/include/parameter_parser/writer/kernel.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
-#include "enumerations/simulation.hpp"
 #include "IO/reader.hpp"
 #include "IO/writer.hpp"
+#include "enumerations/simulation.hpp"
 #include "yaml-cpp/yaml.h"
 
 namespace specfem {
diff --git a/include/parameter_parser/writer/wavefield.hpp b/include/parameter_parser/writer/wavefield.hpp
index 707bda4d..0f98c9d5 100644
--- a/include/parameter_parser/writer/wavefield.hpp
+++ b/include/parameter_parser/writer/wavefield.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
-#include "enumerations/simulation.hpp"
 #include "IO/reader.hpp"
 #include "IO/writer.hpp"
+#include "enumerations/simulation.hpp"
 #include "yaml-cpp/yaml.h"
 
 namespace specfem {
diff --git a/include/point/kernels.hpp b/include/point/kernels.hpp
index 2c14f54f..7683114e 100644
--- a/include/point/kernels.hpp
+++ b/include/point/kernels.hpp
@@ -30,7 +30,7 @@ struct kernels;
  */
 template <bool UseSIMD>
 struct kernels<specfem::dimension::type::dim2,
-               specfem::element::medium_tag::elastic_sv,
+               specfem::element::medium_tag::elastic,
                specfem::element::property_tag::isotropic, UseSIMD> {
 public:
   /**
@@ -50,7 +50,7 @@ struct kernels<specfem::dimension::type::dim2,
    *
    */
   ///@{
-  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic;
   constexpr static auto property_tag =
       specfem::element::property_tag::isotropic;
   constexpr static auto dimension = specfem::dimension::type::dim2;
@@ -135,6 +135,120 @@ struct kernels<specfem::dimension::type::dim2,
 };
 // end elastic isotropic
 
+template <bool UseSIMD>
+struct kernels<specfem::dimension::type::dim2,
+               specfem::element::medium_tag::elastic_sv,
+               specfem::element::property_tag::isotropic, UseSIMD>
+    : public kernels<specfem::dimension::type::dim2,
+                     specfem::element::medium_tag::elastic,
+                     specfem::element::property_tag::isotropic, UseSIMD> {
+
+public:
+  /**
+   * @name Typedefs
+   *
+   */
+  ///@{
+  using simd =
+      typename specfem::datatype::simd<type_real, UseSIMD>; ///< SIMD type
+  using value_type =
+      typename simd::datatype; ///< Underlying data type to store the kernels
+  ///@}
+
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+  /**
+   * @brief Default constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels() = default;
+
+  /**
+   * @brief single value constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type value)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::isotropic, UseSIMD>(value) {}
+
+  /**
+   * @brief Constructor
+   *
+   * @param rho \f$ K_{\rho} \f$
+   * @param mu \f$ K_{\mu} \f$
+   * @param kappa \f$ K_{\kappa} \f$
+   * @param rhop \f$ K_{\rho'} \f$
+   * @param alpha \f$ K_{\alpha} \f$
+   * @param beta \f$ K_{\beta} \f$
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type rho, const value_type mu, const value_type kappa,
+          const value_type rhop, const value_type alpha, const value_type beta)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::isotropic, UseSIMD>(
+            rho, mu, kappa, rhop, alpha, beta) {}
+};
+
+template <bool UseSIMD>
+struct kernels<specfem::dimension::type::dim2,
+               specfem::element::medium_tag::elastic_sh,
+               specfem::element::property_tag::isotropic, UseSIMD>
+    : public kernels<specfem::dimension::type::dim2,
+                     specfem::element::medium_tag::elastic,
+                     specfem::element::property_tag::isotropic, UseSIMD> {
+
+public:
+  /**
+   * @name Typedefs
+   *
+   */
+  ///@{
+  using simd =
+      typename specfem::datatype::simd<type_real, UseSIMD>; ///< SIMD type
+  using value_type =
+      typename simd::datatype; ///< Underlying data type to store the kernels
+  ///@}
+
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh;
+  /**
+   * @brief Default constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels() = default;
+
+  /**
+   * @brief single value constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type value)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::isotropic, UseSIMD>(value) {}
+
+  /**
+   * @brief Constructor
+   *
+   * @param rho \f$ K_{\rho} \f$
+   * @param mu \f$ K_{\mu} \f$
+   * @param kappa \f$ K_{\kappa} \f$
+   * @param rhop \f$ K_{\rho'} \f$
+   * @param alpha \f$ K_{\alpha} \f$
+   * @param beta \f$ K_{\beta} \f$
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type rho, const value_type mu, const value_type kappa,
+          const value_type rhop, const value_type alpha, const value_type beta)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::isotropic, UseSIMD>(
+            rho, mu, kappa, rhop, alpha, beta) {}
+};
+
 /**
  * @brief Template specialization for the kernels struct for 2D elastic
  * anisotropic elements
@@ -143,7 +257,7 @@ struct kernels<specfem::dimension::type::dim2,
  */
 template <bool UseSIMD>
 struct kernels<specfem::dimension::type::dim2,
-               specfem::element::medium_tag::elastic_sv,
+               specfem::element::medium_tag::elastic,
                specfem::element::property_tag::anisotropic, UseSIMD> {
 public:
   /**
@@ -163,7 +277,7 @@ struct kernels<specfem::dimension::type::dim2,
    *
    */
   ///@{
-  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic;
   constexpr static auto property_tag =
       specfem::element::property_tag::anisotropic;
   constexpr static auto dimension = specfem::dimension::type::dim2;
@@ -248,6 +362,123 @@ struct kernels<specfem::dimension::type::dim2,
   KOKKOS_FUNCTION
   bool operator!=(const value_type value) { return !(*this == value); }
 };
+
+template <bool UseSIMD>
+struct kernels<specfem::dimension::type::dim2,
+               specfem::element::medium_tag::elastic_sv,
+               specfem::element::property_tag::anisotropic, UseSIMD>
+    : public kernels<specfem::dimension::type::dim2,
+                     specfem::element::medium_tag::elastic,
+                     specfem::element::property_tag::anisotropic, UseSIMD> {
+
+public:
+  /**
+   * @name Typedefs
+   *
+   */
+  ///@{
+  using simd =
+      typename specfem::datatype::simd<type_real, UseSIMD>; ///< SIMD type
+  using value_type =
+      typename simd::datatype; ///< Underlying data type to store the kernels
+  ///@}
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv;
+  /**
+   * @brief Default constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels() = default;
+
+  /**
+   * @brief single value constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type value)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::anisotropic, UseSIMD>(value) {}
+
+  /**
+   * @brief Constructor
+   *
+   * @param rho \f$ K_{\rho} \f$
+   * @param c11 \f$ K_{c_{11}} \f$
+   * @param c13 \f$ K_{c_{13}} \f$
+   * @param c15 \f$ K_{c_{15}} \f$
+   * @param c33 \f$ K_{c_{33}} \f$
+   * @param c35 \f$ K_{c_{35}} \f$
+   * @param c55 \f$ K_{c_{55}} \f$
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type rho, const value_type c11, const value_type c13,
+          const value_type c15, const value_type c33, const value_type c35,
+          const value_type c55)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::anisotropic, UseSIMD>(
+            rho, c11, c13, c15, c33, c35, c55) {}
+};
+
+template <bool UseSIMD>
+struct kernels<specfem::dimension::type::dim2,
+               specfem::element::medium_tag::elastic_sh,
+               specfem::element::property_tag::anisotropic, UseSIMD>
+    : public kernels<specfem::dimension::type::dim2,
+                     specfem::element::medium_tag::elastic,
+                     specfem::element::property_tag::anisotropic, UseSIMD> {
+
+public:
+  /**
+   * @name Typedefs
+   *
+   */
+  ///@{
+  using simd =
+      typename specfem::datatype::simd<type_real, UseSIMD>; ///< SIMD type
+  using value_type =
+      typename simd::datatype; ///< Underlying data type to store the kernels
+  ///@}
+
+  constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh;
+  /**
+   * @brief Default constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels() = default;
+
+  /**
+   * @brief single value constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type value)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::anisotropic, UseSIMD>(value) {}
+
+  /**
+   * @brief Constructor
+   *
+   * @param rho \f$ K_{\rho} \f$
+   * @param c11 \f$ K_{c_{11}} \f$
+   * @param c13 \f$ K_{c_{13}} \f$
+   * @param c15 \f$ K_{c_{15}} \f$
+   * @param c33 \f$ K_{c_{33}} \f$
+   * @param c35 \f$ K_{c_{35}} \f$
+   * @param c55 \f$ K_{c_{55}} \f$
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type rho, const value_type c11, const value_type c13,
+          const value_type c15, const value_type c33, const value_type c35,
+          const value_type c55)
+      : kernels<specfem::dimension::type::dim2,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::anisotropic, UseSIMD>(
+            rho, c11, c13, c15, c33, c35, c55) {}
+};
 // end elastic anisotropic
 
 /**
diff --git a/src/IO/kernel/writer.cpp b/src/IO/kernel/writer.cpp
index dd3181b1..4110b80a 100644
--- a/src/IO/kernel/writer.cpp
+++ b/src/IO/kernel/writer.cpp
@@ -5,6 +5,8 @@
 
 // Explicit instantiation
 
-template class specfem::IO::kernel_writer<specfem::IO::HDF5<specfem::IO::write> >;
+template class specfem::IO::kernel_writer<
+    specfem::IO::HDF5<specfem::IO::write> >;
 
-template class specfem::IO::kernel_writer<specfem::IO::ASCII<specfem::IO::write> >;
+template class specfem::IO::kernel_writer<
+    specfem::IO::ASCII<specfem::IO::write> >;
diff --git a/src/IO/mesh.cpp b/src/IO/mesh.cpp
index 3ad22eb5..8a24385a 100644
--- a/src/IO/mesh.cpp
+++ b/src/IO/mesh.cpp
@@ -196,16 +196,24 @@ specfem::IO::read_mesh(const std::string filename,
 
   const auto l_elastic_sv_isotropic =
       mesh.materials.elastic_sv_isotropic.material_properties;
+  const auto l_elastic_sh_isotropic =
+      mesh.materials.elastic_sh_isotropic.material_properties;
   const auto l_acoustic_isotropic =
       mesh.materials.acoustic_isotropic.material_properties;
 
   const auto l_elastic_sv_anisotropic =
       mesh.materials.elastic_sv_anisotropic.material_properties;
+  const auto l_elastic_sh_anisotropic =
+      mesh.materials.elastic_sh_anisotropic.material_properties;
 
   for (const auto material : l_elastic_sv_isotropic) {
     mpi->cout(material.print());
   }
 
+  for (const auto material : l_elastic_sh_isotropic) {
+    mpi->cout(material.print());
+  }
+
   for (const auto material : l_acoustic_isotropic) {
     mpi->cout(material.print());
   }
@@ -214,9 +222,13 @@ specfem::IO::read_mesh(const std::string filename,
     mpi->cout(material.print());
   }
 
-  assert(l_elastic_sv_isotropic.size() + l_acoustic_isotropic.size() +
-             l_elastic_sv_anisotropic.size() ==
-         mesh.materials.n_materials);
+  for (const auto material : l_elastic_sh_anisotropic) {
+    mpi->cout(material.print());
+  }
+
+  assert((l_elastic_sv_isotropic.size() + l_elastic_sh_isotropic.size() +
+          l_acoustic_isotropic.size() + l_elastic_sv_anisotropic.size() +
+          l_elastic_sh_anisotropic.size()) == mesh.materials.n_materials);
 
   mesh.tags = specfem::mesh::tags<specfem::dimension::type::dim2>(
       mesh.materials, mesh.boundaries);
diff --git a/src/IO/property/reader.cpp b/src/IO/property/reader.cpp
index 8e1cb19a..ca7bdcb4 100644
--- a/src/IO/property/reader.cpp
+++ b/src/IO/property/reader.cpp
@@ -5,7 +5,8 @@
 #include "IO/reader.hpp"
 
 // Explicit instantiation
-template class specfem::IO::property_reader<specfem::IO::HDF5<specfem::IO::read> >;
+template class specfem::IO::property_reader<
+    specfem::IO::HDF5<specfem::IO::read> >;
 
 template class specfem::IO::property_reader<
     specfem::IO::ASCII<specfem::IO::read> >;
diff --git a/src/compute/compute_kernels.cpp b/src/compute/compute_kernels.cpp
index 2a70c000..73e2c3f6 100644
--- a/src/compute/compute_kernels.cpp
+++ b/src/compute/compute_kernels.cpp
@@ -36,16 +36,26 @@ specfem::compute::kernels::kernels(
       specfem::element::property_tag::isotropic>(
       acoustic_elements, ngllz, ngllx, h_property_index_mapping);
 
-  elastic_isotropic = specfem::medium::material_kernels<
+  elastic_sv_isotropic = specfem::medium::material_kernels<
       specfem::element::medium_tag::elastic_sv,
       specfem::element::property_tag::isotropic>(
       elastic_isotropic_elements, ngllz, ngllx, h_property_index_mapping);
 
-  elastic_anisotropic = specfem::medium::material_kernels<
+  elastic_sh_isotropic = specfem::medium::material_kernels<
+      specfem::element::medium_tag::elastic_sh,
+      specfem::element::property_tag::isotropic>(
+      elastic_isotropic_elements, ngllz, ngllx, h_property_index_mapping);
+
+  elastic_sv_anisotropic = specfem::medium::material_kernels<
       specfem::element::medium_tag::elastic_sv,
       specfem::element::property_tag::anisotropic>(
       elastic_anisotropic_elements, ngllz, ngllx, h_property_index_mapping);
 
+  elastic_sh_anisotropic = specfem::medium::material_kernels<
+      specfem::element::medium_tag::elastic_sh,
+      specfem::element::property_tag::anisotropic>(
+      elastic_anisotropic_elements, ngllz, ngllx, h_property_index_mapping);
+
   Kokkos::deep_copy(property_index_mapping, h_property_index_mapping);
 
   return;
diff --git a/src/compute/compute_properties.cpp b/src/compute/compute_properties.cpp
index 4319e5e4..b69312e9 100644
--- a/src/compute/compute_properties.cpp
+++ b/src/compute/compute_properties.cpp
@@ -15,13 +15,23 @@ specfem::compute::properties::properties(
   this->h_property_index_mapping =
       Kokkos::create_mirror_view(property_index_mapping);
 
-  const auto elastic_isotropic_elements = element_types.get_elements_on_host(
+  const auto elastic_sv_isotropic_elements = element_types.get_elements_on_host(
       specfem::element::medium_tag::elastic_sv,
       specfem::element::property_tag::isotropic);
 
-  const auto elastic_anisotropic_elements = element_types.get_elements_on_host(
-      specfem::element::medium_tag::elastic_sv,
-      specfem::element::property_tag::anisotropic);
+  const auto elastic_sv_anisotropic_elements =
+      element_types.get_elements_on_host(
+          specfem::element::medium_tag::elastic_sv,
+          specfem::element::property_tag::anisotropic);
+
+  const auto elastic_sh_isotropic_elements = element_types.get_elements_on_host(
+      specfem::element::medium_tag::elastic_sh,
+      specfem::element::property_tag::isotropic);
+
+  const auto elastic_sh_anisotropic_elements =
+      element_types.get_elements_on_host(
+          specfem::element::medium_tag::elastic_sh,
+          specfem::element::property_tag::anisotropic);
 
   const auto acoustic_elements = element_types.get_elements_on_host(
       specfem::element::medium_tag::acoustic,
@@ -31,22 +41,34 @@ specfem::compute::properties::properties(
     h_property_index_mapping(ispec) = -1;
   }
 
-  acoustic_isotropic = specfem::medium::material_properties<
+  acoustic_isotropic = specfem::medium::properties_container<
       specfem::element::medium_tag::acoustic,
       specfem::element::property_tag::isotropic>(
       acoustic_elements, ngllz, ngllx, materials, has_gll_model,
       h_property_index_mapping);
 
-  elastic_isotropic = specfem::medium::material_properties<
+  elastic_sv_isotropic = specfem::medium::properties_container<
       specfem::element::medium_tag::elastic_sv,
       specfem::element::property_tag::isotropic>(
-      elastic_isotropic_elements, ngllz, ngllx, materials, has_gll_model,
+      elastic_sv_isotropic_elements, ngllz, ngllx, materials, has_gll_model,
       h_property_index_mapping);
 
-  elastic_anisotropic = specfem::medium::material_properties<
+  elastic_sv_anisotropic = specfem::medium::properties_container<
       specfem::element::medium_tag::elastic_sv,
       specfem::element::property_tag::anisotropic>(
-      elastic_anisotropic_elements, ngllz, ngllx, materials, has_gll_model,
+      elastic_sv_anisotropic_elements, ngllz, ngllx, materials, has_gll_model,
+      h_property_index_mapping);
+
+  elastic_sh_isotropic = specfem::medium::properties_container<
+      specfem::element::medium_tag::elastic_sh,
+      specfem::element::property_tag::isotropic>(
+      elastic_sh_isotropic_elements, ngllz, ngllx, materials, has_gll_model,
+      h_property_index_mapping);
+
+  elastic_sh_anisotropic = specfem::medium::properties_container<
+      specfem::element::medium_tag::elastic_sh,
+      specfem::element::property_tag::anisotropic>(
+      elastic_sh_anisotropic_elements, ngllz, ngllx, materials, has_gll_model,
       h_property_index_mapping);
 
   Kokkos::deep_copy(property_index_mapping, h_property_index_mapping);
diff --git a/src/compute/compute_receivers.cpp b/src/compute/compute_receivers.cpp
index 383776c4..a373d691 100644
--- a/src/compute/compute_receivers.cpp
+++ b/src/compute/compute_receivers.cpp
@@ -112,9 +112,9 @@ specfem::compute::receivers::receivers(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       COUNT_RECEIVERS_PER_MATERIAL_SYSTEM,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef COUNT_RECEIVERS_PER_MATERIAL_SYSTEM
 
@@ -145,9 +145,9 @@ specfem::compute::receivers::receivers(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       ALLOCATE_RECEIVERS_PER_MATERIAL_SYSTEM,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef ALLOCATE_RECEIVERS_PER_MATERIAL_SYSTEM
 
@@ -187,9 +187,9 @@ specfem::compute::receivers::receivers(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       ASSIGN_RECEIVERS_PER_MATERIAL_SYSTEM,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef ASSIGN_RECEIVERS_PER_MATERIAL_SYSTEM
 
@@ -216,9 +216,10 @@ specfem::compute::receivers::get_indices_on_host(
   }
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
-      RETURN_VALUE, WHERE(DIMENSION_TAG_DIM2) WHERE(MEDIUM_TAG_ELASTIC_SV,
-                                                    MEDIUM_TAG_ACOUSTIC)
-                        WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      RETURN_VALUE,
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef RETURN_VALUE
 }
@@ -240,9 +241,10 @@ specfem::compute::receivers::get_indices_on_device(
   }
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
-      RETURN_VALUE, WHERE(DIMENSION_TAG_DIM2) WHERE(MEDIUM_TAG_ELASTIC_SV,
-                                                    MEDIUM_TAG_ACOUSTIC)
-                        WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      RETURN_VALUE,
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef RETURN_VALUE
 }
diff --git a/src/compute/compute_sources.cpp b/src/compute/compute_sources.cpp
index 7ef17520..9b63e58e 100644
--- a/src/compute/compute_sources.cpp
+++ b/src/compute/compute_sources.cpp
@@ -101,8 +101,9 @@ specfem::compute::sources::sources(
           sources, element_types, mesh);
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
-      SORT_SOURCES_PER_MEDIUM, WHERE(DIMENSION_TAG_DIM2) WHERE(
-                                   MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+      SORT_SOURCES_PER_MEDIUM,
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC))
 
 #undef SORT_SOURCES_PER_MEDIUM
 
@@ -119,9 +120,10 @@ specfem::compute::sources::sources(
                            GET_NAME(MEDIUM_TAG))                               \
           .size();
 
-  CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
-      COUNT_SOURCES, WHERE(DIMENSION_TAG_DIM2)
-                         WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+  CALL_MACRO_FOR_ALL_MEDIUM_TAGS(COUNT_SOURCES, WHERE(DIMENSION_TAG_DIM2)
+                                                    WHERE(MEDIUM_TAG_ELASTIC_SV,
+                                                          MEDIUM_TAG_ELASTIC_SH,
+                                                          MEDIUM_TAG_ACOUSTIC))
 
 #undef COUNT_SOURCES
 
@@ -184,7 +186,8 @@ specfem::compute::sources::sources(
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
       ASSIGN_MEMBERS, WHERE(DIMENSION_TAG_DIM2)
-                          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+                          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH,
+                                MEDIUM_TAG_ACOUSTIC))
 
 #undef ASSIGN_MEMBERS
 
@@ -229,11 +232,11 @@ specfem::compute::sources::sources(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       COUNT_SOURCES_PER_ELEMENT_TYPE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef COUNT_SOURCES_PER_ELEMENT_TYPE
 
@@ -339,11 +342,11 @@ specfem::compute::sources::sources(
   // Creating the views for all the sources
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       ALLOCATE_SOURCES_PER_ELEMENT_TYPE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef ALLOCATE_SOURCES_PER_ELEMENT_TYPE
 
@@ -494,11 +497,11 @@ specfem::compute::sources::sources(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       ASSIGN_SOURCES_PER_ELEMENT_TYPE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef ASSIGN_SOURCES_PER_ELEMENT_TYPE
 
@@ -556,11 +559,11 @@ specfem::compute::sources::get_sources_on_host(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       RETURN_VALUE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef RETURN_VALUE
 }
@@ -616,11 +619,11 @@ specfem::compute::sources::get_sources_on_device(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       RETURN_VALUE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef RETURN_VALUE
 }
diff --git a/src/compute/element_types/element_types.cpp b/src/compute/element_types/element_types.cpp
index 97d23efa..d80f0891 100644
--- a/src/compute/element_types/element_types.cpp
+++ b/src/compute/element_types/element_types.cpp
@@ -28,8 +28,9 @@ specfem::compute::element_types::element_types(
   }
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
-      COUNT_MEDIUM_TAG_INDICES, WHERE(DIMENSION_TAG_DIM2) WHERE(
-                                    MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+      COUNT_MEDIUM_TAG_INDICES,
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC))
 
 #undef COUNT_MEDIUM_TAG_INDICES
 
@@ -46,8 +47,8 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
       ALLOCATE_MEDIUM_TAG_VIEWS,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC))
 
 #undef ALLOCATE_MEDIUM_TAG_VIEWS
 
@@ -72,8 +73,8 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
       ASSIGN_MEDIUM_TAG_INDICES,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC))
 
 #undef ASSIGN_MEDIUM_TAG_INDICES
 
@@ -91,9 +92,9 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       COUNT_MATERIAL_SYSTEM_INDICES,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef COUNT_MATERIAL_SYSTEM_INDICES
 
@@ -113,9 +114,9 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       ALLOCATE_MATERIAL_SYSTEM_VIEWS,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef ALLOCATE_MATERIAL_SYSTEM_VIEWS
 
@@ -145,9 +146,9 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       ASSIGN_MATERIAL_SYSTEM_INDICES,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef ASSIGN_MATERIAL_SYSTEM_INDICES
 
@@ -169,11 +170,11 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       COUNT_ELEMENT_TYPES_INDICES,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef COUNT_ELEMENT_TYPES_INDICES
 
@@ -195,11 +196,11 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       ALLOCATE_ELEMENT_TYPES_VIEWS,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef ALLOCATE_ELEMENT_TYPES_VIEWS
 
@@ -233,11 +234,11 @@ specfem::compute::element_types::element_types(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       ASSIGN_ELEMENT_TYPES_INDICES,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef ASSIGN_ELEMENT_TYPES_INDICES
 }
@@ -254,9 +255,12 @@ specfem::compute::element_types::get_elements_on_host(
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
       RETURN_VARIABLE, WHERE(DIMENSION_TAG_DIM2)
-                           WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+                           WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH,
+                                 MEDIUM_TAG_ACOUSTIC))
 
 #undef RETURN_VARIABLE
+
+  throw std::runtime_error("Medium tag not found");
 }
 
 Kokkos::View<int *, Kokkos::DefaultExecutionSpace>
@@ -271,9 +275,12 @@ specfem::compute::element_types::get_elements_on_device(
 
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(
       RETURN_VARIABLE, WHERE(DIMENSION_TAG_DIM2)
-                           WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC))
+                           WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH,
+                                 MEDIUM_TAG_ACOUSTIC))
 
 #undef RETURN_VARIABLE
+
+  throw std::runtime_error("Medium tag not found");
 }
 
 Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace>
@@ -291,11 +298,13 @@ specfem::compute::element_types::get_elements_on_host(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       RETURN_VARIABLE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef RETURN_VARIABLE
+
+  throw std::runtime_error("Medium tag or property tag not found");
 }
 
 Kokkos::View<int *, Kokkos::DefaultExecutionSpace>
@@ -313,11 +322,13 @@ specfem::compute::element_types::get_elements_on_device(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       RETURN_VARIABLE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef RETURN_VARIABLE
+
+  throw std::runtime_error("Medium tag or property tag not found");
 }
 
 Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace>
@@ -337,13 +348,16 @@ specfem::compute::element_types::get_elements_on_host(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       RETURN_VARIABLE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef RETURN_VARIABLE
+
+  throw std::runtime_error(
+      "Medium tag, property tag or boundary tag not found");
 }
 
 Kokkos::View<int *, Kokkos::DefaultExecutionSpace>
@@ -363,11 +377,14 @@ specfem::compute::element_types::get_elements_on_device(
 
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       RETURN_VARIABLE,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
-                  BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
-                  BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC) WHERE(
+              BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
+              BOUNDARY_TAG_STACEY, BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET))
 
 #undef RETURN_VARIABLE
+
+  throw std::runtime_error(
+      "Medium tag, property tag or boundary tag not found");
 }
diff --git a/src/parameter_parser/writer/kernel.cpp b/src/parameter_parser/writer/kernel.cpp
index 589f68a7..b0dfa893 100644
--- a/src/parameter_parser/writer/kernel.cpp
+++ b/src/parameter_parser/writer/kernel.cpp
@@ -43,13 +43,11 @@ specfem::runtime_configuration::kernel::instantiate_kernel_writer() const {
       [&]() -> std::shared_ptr<specfem::IO::writer> {
     if (this->simulation_type == specfem::simulation::type::combined) {
       if (this->output_format == "HDF5") {
-        return std::make_shared<
-            specfem::IO::kernel_writer<specfem::IO::HDF5<specfem::IO::write> > >(
-            this->output_folder);
+        return std::make_shared<specfem::IO::kernel_writer<
+            specfem::IO::HDF5<specfem::IO::write> > >(this->output_folder);
       } else if (this->output_format == "ASCII") {
-        return std::make_shared<
-            specfem::IO::kernel_writer<specfem::IO::ASCII<specfem::IO::write> > >(
-            this->output_folder);
+        return std::make_shared<specfem::IO::kernel_writer<
+            specfem::IO::ASCII<specfem::IO::write> > >(this->output_folder);
       } else {
         throw std::runtime_error("Unknown wavefield format");
       }
diff --git a/src/parameter_parser/writer/property.cpp b/src/parameter_parser/writer/property.cpp
index 63aa3305..11affbdd 100644
--- a/src/parameter_parser/writer/property.cpp
+++ b/src/parameter_parser/writer/property.cpp
@@ -46,13 +46,11 @@ specfem::runtime_configuration::property::instantiate_property_writer() const {
       return nullptr;
     }
     if (this->output_format == "HDF5") {
-      return std::make_shared<
-          specfem::IO::property_writer<specfem::IO::HDF5<specfem::IO::write> > >(
-          this->output_folder);
+      return std::make_shared<specfem::IO::property_writer<
+          specfem::IO::HDF5<specfem::IO::write> > >(this->output_folder);
     } else if (this->output_format == "ASCII") {
-      return std::make_shared<
-          specfem::IO::property_writer<specfem::IO::ASCII<specfem::IO::write> > >(
-          this->output_folder);
+      return std::make_shared<specfem::IO::property_writer<
+          specfem::IO::ASCII<specfem::IO::write> > >(this->output_folder);
     } else {
       throw std::runtime_error("Unknown model format");
     }
@@ -74,9 +72,8 @@ specfem::runtime_configuration::property::instantiate_property_reader() const {
           specfem::IO::property_reader<specfem::IO::HDF5<specfem::IO::read> > >(
           this->output_folder);
     } else if (this->output_format == "ASCII") {
-      return std::make_shared<
-          specfem::IO::property_reader<specfem::IO::ASCII<specfem::IO::read> > >(
-          this->output_folder);
+      return std::make_shared<specfem::IO::property_reader<
+          specfem::IO::ASCII<specfem::IO::read> > >(this->output_folder);
     } else {
       throw std::runtime_error("Unknown model format");
     }
diff --git a/src/parameter_parser/writer/wavefield.cpp b/src/parameter_parser/writer/wavefield.cpp
index a09f7d85..edb6a14c 100644
--- a/src/parameter_parser/writer/wavefield.cpp
+++ b/src/parameter_parser/writer/wavefield.cpp
@@ -70,9 +70,8 @@ specfem::runtime_configuration::wavefield::instantiate_wavefield_reader()
       [&]() -> std::shared_ptr<specfem::IO::reader> {
     if (this->simulation_type == specfem::simulation::type::combined) {
       if (this->output_format == "HDF5") {
-        return std::make_shared<
-            specfem::IO::wavefield_reader<specfem::IO::HDF5<specfem::IO::read> > >(
-            this->output_folder);
+        return std::make_shared<specfem::IO::wavefield_reader<
+            specfem::IO::HDF5<specfem::IO::read> > >(this->output_folder);
       } else if (this->output_format == "ASCII") {
         return std::make_shared<specfem::IO::wavefield_reader<
             specfem::IO::ASCII<specfem::IO::read> > >(this->output_folder);
diff --git a/src/quadrature/gll/gll_utils.cpp b/src/quadrature/gll/gll_utils.cpp
index 135f7649..be0e7d1c 100644
--- a/src/quadrature/gll/gll_utils.cpp
+++ b/src/quadrature/gll/gll_utils.cpp
@@ -59,7 +59,8 @@ void specfem::quadrature::gll::gll_utils::jacg(HostMirror1d xjac, const int np,
 
     if (np - j < 0 || np - j > np - 1) {
       std::ostringstream oss;
-      oss << "ERROR :Index xjac is out of range: xjac.extent(0) == " << np << " & np-j =" << np - j;
+      oss << "ERROR :Index xjac is out of range: xjac.extent(0) == " << np
+          << " & np-j =" << np - j;
       throw std::runtime_error(oss.str());
     }
     xjac(np - j) = x;
diff --git a/src/source_time_function/external.cpp b/src/source_time_function/external.cpp
index ffac2fac..dbcd6338 100644
--- a/src/source_time_function/external.cpp
+++ b/src/source_time_function/external.cpp
@@ -1,7 +1,7 @@
 #include "source_time_function/external.hpp"
+#include "IO/seismogram/reader.hpp"
 #include "enumerations/specfem_enums.hpp"
 #include "kokkos_abstractions.h"
-#include "IO/seismogram/reader.hpp"
 #include <fstream>
 #include <tuple>
 #include <vector>
diff --git a/tests/unit-tests/assembly/kernels/kernels.cpp b/tests/unit-tests/assembly/kernels/kernels.cpp
index 1973db64..394827e6 100644
--- a/tests/unit-tests/assembly/kernels/kernels.cpp
+++ b/tests/unit-tests/assembly/kernels/kernels.cpp
@@ -16,10 +16,9 @@ std::string get_error_message(
 
 template <>
 std::string get_error_message(
-    const specfem::point::kernels<specfem::dimension::type::dim2,
-                                  specfem::element::medium_tag::elastic_sv,
-                                  specfem::element::property_tag::isotropic,
-                                  false> &point_kernel,
+    const specfem::point::kernels<
+        specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+        specfem::element::property_tag::isotropic, false> &point_kernel,
     const type_real value) {
   std::ostringstream message;
 
@@ -39,9 +38,38 @@ template <>
 std::string get_error_message(
     const specfem::point::kernels<specfem::dimension::type::dim2,
                                   specfem::element::medium_tag::elastic_sv,
-                                  specfem::element::property_tag::anisotropic,
+                                  specfem::element::property_tag::isotropic,
+                                  false> &point_kernel,
+    const type_real value) {
+
+  return get_error_message(
+      static_cast<specfem::point::kernels<
+          specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+          specfem::element::property_tag::isotropic, false> >(point_kernel),
+      value);
+}
+
+template <>
+std::string get_error_message(
+    const specfem::point::kernels<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::elastic_sh,
+                                  specfem::element::property_tag::isotropic,
                                   false> &point_kernel,
     const type_real value) {
+
+  return get_error_message(
+      static_cast<specfem::point::kernels<
+          specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+          specfem::element::property_tag::isotropic, false> >(point_kernel),
+      value);
+}
+
+template <>
+std::string get_error_message(
+    const specfem::point::kernels<
+        specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+        specfem::element::property_tag::anisotropic, false> &point_kernel,
+    const type_real value) {
   std::ostringstream message;
 
   message << "\n\t Expected: " << value;
@@ -57,6 +85,34 @@ std::string get_error_message(
   return message.str();
 }
 
+template <>
+std::string get_error_message(
+    const specfem::point::kernels<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::elastic_sv,
+                                  specfem::element::property_tag::anisotropic,
+                                  false> &point_kernel,
+    const type_real value) {
+  return get_error_message(
+      static_cast<specfem::point::kernels<
+          specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+          specfem::element::property_tag::anisotropic, false> >(point_kernel),
+      value);
+}
+
+template <>
+std::string get_error_message(
+    const specfem::point::kernels<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::elastic_sh,
+                                  specfem::element::property_tag::anisotropic,
+                                  false> &point_kernel,
+    const type_real value) {
+  return get_error_message(
+      static_cast<specfem::point::kernels<
+          specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+          specfem::element::property_tag::anisotropic, false> >(point_kernel),
+      value);
+}
+
 template <>
 std::string get_error_message(
     const specfem::point::kernels<
@@ -98,7 +154,7 @@ specfem::point::kernels<specfem::dimension::type::dim2,
 get_point_kernel(const int ispec, const int iz, const int ix,
                  const specfem::compute::kernels &kernels) {
 
-  const auto elastic_isotropic = kernels.elastic_isotropic;
+  const auto elastic_isotropic = kernels.elastic_sv_isotropic;
 
   const int ispec_l = kernels.h_property_index_mapping(ispec);
 
@@ -117,6 +173,32 @@ get_point_kernel(const int ispec, const int iz, const int ix,
   return point_kernel;
 }
 
+template <>
+specfem::point::kernels<specfem::dimension::type::dim2,
+                        specfem::element::medium_tag::elastic_sh,
+                        specfem::element::property_tag::isotropic, false>
+get_point_kernel(const int ispec, const int iz, const int ix,
+                 const specfem::compute::kernels &kernels) {
+
+  const auto elastic_isotropic = kernels.elastic_sh_isotropic;
+
+  const int ispec_l = kernels.h_property_index_mapping(ispec);
+
+  specfem::point::kernels<specfem::dimension::type::dim2,
+                          specfem::element::medium_tag::elastic_sh,
+                          specfem::element::property_tag::isotropic, false>
+      point_kernel;
+
+  point_kernel.rho = elastic_isotropic.h_rho(ispec_l, iz, ix);
+  point_kernel.mu = elastic_isotropic.h_mu(ispec_l, iz, ix);
+  point_kernel.kappa = elastic_isotropic.h_kappa(ispec_l, iz, ix);
+  point_kernel.rhop = elastic_isotropic.h_rhop(ispec_l, iz, ix);
+  point_kernel.alpha = elastic_isotropic.h_alpha(ispec_l, iz, ix);
+  point_kernel.beta = elastic_isotropic.h_beta(ispec_l, iz, ix);
+
+  return point_kernel;
+}
+
 template <>
 specfem::point::kernels<specfem::dimension::type::dim2,
                         specfem::element::medium_tag::elastic_sv,
@@ -142,6 +224,31 @@ get_point_kernel(
   return point_kernel_l;
 }
 
+template <>
+specfem::point::kernels<specfem::dimension::type::dim2,
+                        specfem::element::medium_tag::elastic_sh,
+                        specfem::element::property_tag::isotropic, false>
+get_point_kernel(
+    const int lane,
+    const specfem::point::kernels<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::elastic_sh,
+                                  specfem::element::property_tag::isotropic,
+                                  true> &point_kernel) {
+  specfem::point::kernels<specfem::dimension::type::dim2,
+                          specfem::element::medium_tag::elastic_sh,
+                          specfem::element::property_tag::isotropic, false>
+      point_kernel_l;
+
+  point_kernel_l.rho = point_kernel.rho[lane];
+  point_kernel_l.mu = point_kernel.mu[lane];
+  point_kernel_l.kappa = point_kernel.kappa[lane];
+  point_kernel_l.rhop = point_kernel.rhop[lane];
+  point_kernel_l.alpha = point_kernel.alpha[lane];
+  point_kernel_l.beta = point_kernel.beta[lane];
+
+  return point_kernel_l;
+}
+
 template <>
 specfem::point::kernels<specfem::dimension::type::dim2,
                         specfem::element::medium_tag::elastic_sv,
@@ -149,7 +256,7 @@ specfem::point::kernels<specfem::dimension::type::dim2,
 get_point_kernel(const int ispec, const int iz, const int ix,
                  const specfem::compute::kernels &kernels) {
 
-  const auto elastic_anisotropic = kernels.elastic_anisotropic;
+  const auto elastic_anisotropic = kernels.elastic_sv_anisotropic;
 
   const int ispec_l = kernels.h_property_index_mapping(ispec);
 
@@ -169,6 +276,33 @@ get_point_kernel(const int ispec, const int iz, const int ix,
   return point_kernel;
 }
 
+template <>
+specfem::point::kernels<specfem::dimension::type::dim2,
+                        specfem::element::medium_tag::elastic_sh,
+                        specfem::element::property_tag::anisotropic, false>
+get_point_kernel(const int ispec, const int iz, const int ix,
+                 const specfem::compute::kernels &kernels) {
+
+  const auto elastic_anisotropic = kernels.elastic_sh_anisotropic;
+
+  const int ispec_l = kernels.h_property_index_mapping(ispec);
+
+  specfem::point::kernels<specfem::dimension::type::dim2,
+                          specfem::element::medium_tag::elastic_sh,
+                          specfem::element::property_tag::anisotropic, false>
+      point_kernel;
+
+  point_kernel.rho = elastic_anisotropic.h_rho(ispec_l, iz, ix);
+  point_kernel.c11 = elastic_anisotropic.h_c11(ispec_l, iz, ix);
+  point_kernel.c13 = elastic_anisotropic.h_c13(ispec_l, iz, ix);
+  point_kernel.c15 = elastic_anisotropic.h_c15(ispec_l, iz, ix);
+  point_kernel.c33 = elastic_anisotropic.h_c33(ispec_l, iz, ix);
+  point_kernel.c35 = elastic_anisotropic.h_c35(ispec_l, iz, ix);
+  point_kernel.c55 = elastic_anisotropic.h_c55(ispec_l, iz, ix);
+
+  return point_kernel;
+}
+
 template <>
 specfem::point::kernels<specfem::dimension::type::dim2,
                         specfem::element::medium_tag::elastic_sv,
@@ -195,6 +329,32 @@ get_point_kernel(
   return point_kernel_l;
 }
 
+template <>
+specfem::point::kernels<specfem::dimension::type::dim2,
+                        specfem::element::medium_tag::elastic_sh,
+                        specfem::element::property_tag::anisotropic, false>
+get_point_kernel(
+    const int lane,
+    const specfem::point::kernels<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::elastic_sh,
+                                  specfem::element::property_tag::anisotropic,
+                                  true> &point_kernel) {
+  specfem::point::kernels<specfem::dimension::type::dim2,
+                          specfem::element::medium_tag::elastic_sh,
+                          specfem::element::property_tag::anisotropic, false>
+      point_kernel_l;
+
+  point_kernel_l.rho = point_kernel.rho[lane];
+  point_kernel_l.c11 = point_kernel.c11[lane];
+  point_kernel_l.c13 = point_kernel.c13[lane];
+  point_kernel_l.c15 = point_kernel.c15[lane];
+  point_kernel_l.c33 = point_kernel.c33[lane];
+  point_kernel_l.c35 = point_kernel.c35[lane];
+  point_kernel_l.c55 = point_kernel.c55[lane];
+
+  return point_kernel_l;
+}
+
 template <>
 specfem::point::kernels<specfem::dimension::type::dim2,
                         specfem::element::medium_tag::acoustic,
@@ -511,9 +671,9 @@ void test_kernels(specfem::compute::assembly &assembly) {
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       TEST_STORE_AND_ADD,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef TEST_STORE_AND_ADD
 }
diff --git a/tests/unit-tests/assembly/properties/properties.cpp b/tests/unit-tests/assembly/properties/properties.cpp
index 242c309e..6ccae91f 100644
--- a/tests/unit-tests/assembly/properties/properties.cpp
+++ b/tests/unit-tests/assembly/properties/properties.cpp
@@ -163,7 +163,7 @@ specfem::point::properties<specfem::dimension::type::dim2,
 get_point_property(const int ispec, const int iz, const int ix,
                    const specfem::compute::properties &properties) {
 
-  const auto elastic_isotropic = properties.elastic_isotropic;
+  const auto elastic_isotropic = properties.elastic_sv_isotropic;
 
   const int ispec_l = properties.h_property_index_mapping(ispec);
 
@@ -180,6 +180,30 @@ get_point_property(const int ispec, const int iz, const int ix,
   return point_property;
 }
 
+template <>
+specfem::point::properties<specfem::dimension::type::dim2,
+                           specfem::element::medium_tag::elastic_sh,
+                           specfem::element::property_tag::isotropic, false>
+get_point_property(const int ispec, const int iz, const int ix,
+                   const specfem::compute::properties &properties) {
+
+  const auto elastic_isotropic = properties.elastic_sh_isotropic;
+
+  const int ispec_l = properties.h_property_index_mapping(ispec);
+
+  specfem::point::properties<specfem::dimension::type::dim2,
+                             specfem::element::medium_tag::elastic_sh,
+                             specfem::element::property_tag::isotropic, false>
+      point_property;
+
+  point_property.rho = elastic_isotropic.h_rho(ispec_l, iz, ix);
+  point_property.mu = elastic_isotropic.h_mu(ispec_l, iz, ix);
+  point_property.lambdaplus2mu =
+      elastic_isotropic.h_lambdaplus2mu(ispec_l, iz, ix);
+
+  return point_property;
+}
+
 template <>
 specfem::point::properties<specfem::dimension::type::dim2,
                            specfem::element::medium_tag::elastic_sv,
@@ -202,6 +226,28 @@ get_point_property(
   return point_property_l;
 }
 
+template <>
+specfem::point::properties<specfem::dimension::type::dim2,
+                           specfem::element::medium_tag::elastic_sh,
+                           specfem::element::property_tag::isotropic, false>
+get_point_property(
+    const int lane,
+    const specfem::point::properties<specfem::dimension::type::dim2,
+                                     specfem::element::medium_tag::elastic_sh,
+                                     specfem::element::property_tag::isotropic,
+                                     true> &point_property) {
+  specfem::point::properties<specfem::dimension::type::dim2,
+                             specfem::element::medium_tag::elastic_sh,
+                             specfem::element::property_tag::isotropic, false>
+      point_property_l;
+
+  point_property_l.rho = point_property.rho[lane];
+  point_property_l.mu = point_property.mu[lane];
+  point_property_l.lambdaplus2mu = point_property.lambdaplus2mu[lane];
+
+  return point_property_l;
+}
+
 template <>
 specfem::point::properties<specfem::dimension::type::dim2,
                            specfem::element::medium_tag::elastic_sv,
@@ -209,7 +255,7 @@ specfem::point::properties<specfem::dimension::type::dim2,
 get_point_property(const int ispec, const int iz, const int ix,
                    const specfem::compute::properties &properties) {
 
-  const auto elastic_anisotropic = properties.elastic_anisotropic;
+  const auto elastic_anisotropic = properties.elastic_sv_anisotropic;
 
   const int ispec_l = properties.h_property_index_mapping(ispec);
 
@@ -229,6 +275,33 @@ get_point_property(const int ispec, const int iz, const int ix,
   return point_property;
 }
 
+template <>
+specfem::point::properties<specfem::dimension::type::dim2,
+                           specfem::element::medium_tag::elastic_sh,
+                           specfem::element::property_tag::anisotropic, false>
+get_point_property(const int ispec, const int iz, const int ix,
+                   const specfem::compute::properties &properties) {
+
+  const auto elastic_anisotropic = properties.elastic_sh_anisotropic;
+
+  const int ispec_l = properties.h_property_index_mapping(ispec);
+
+  specfem::point::properties<specfem::dimension::type::dim2,
+                             specfem::element::medium_tag::elastic_sh,
+                             specfem::element::property_tag::anisotropic, false>
+      point_property;
+
+  point_property.rho = elastic_anisotropic.h_rho(ispec_l, iz, ix);
+  point_property.c11 = elastic_anisotropic.h_c11(ispec_l, iz, ix);
+  point_property.c13 = elastic_anisotropic.h_c13(ispec_l, iz, ix);
+  point_property.c15 = elastic_anisotropic.h_c15(ispec_l, iz, ix);
+  point_property.c33 = elastic_anisotropic.h_c33(ispec_l, iz, ix);
+  point_property.c35 = elastic_anisotropic.h_c35(ispec_l, iz, ix);
+  point_property.c55 = elastic_anisotropic.h_c55(ispec_l, iz, ix);
+
+  return point_property;
+}
+
 template <>
 specfem::point::properties<specfem::dimension::type::dim2,
                            specfem::element::medium_tag::elastic_sv,
@@ -255,6 +328,32 @@ get_point_property(
   return point_property_l;
 }
 
+template <>
+specfem::point::properties<specfem::dimension::type::dim2,
+                           specfem::element::medium_tag::elastic_sh,
+                           specfem::element::property_tag::anisotropic, false>
+get_point_property(
+    const int lane,
+    const specfem::point::properties<
+        specfem::dimension::type::dim2,
+        specfem::element::medium_tag::elastic_sh,
+        specfem::element::property_tag::anisotropic, true> &point_property) {
+  specfem::point::properties<specfem::dimension::type::dim2,
+                             specfem::element::medium_tag::elastic_sh,
+                             specfem::element::property_tag::anisotropic, false>
+      point_property_l;
+
+  point_property_l.rho = point_property.rho[lane];
+  point_property_l.c11 = point_property.c11[lane];
+  point_property_l.c13 = point_property.c13[lane];
+  point_property_l.c15 = point_property.c15[lane];
+  point_property_l.c33 = point_property.c33[lane];
+  point_property_l.c35 = point_property.c35[lane];
+  point_property_l.c55 = point_property.c55[lane];
+
+  return point_property_l;
+}
+
 template <>
 specfem::point::properties<specfem::dimension::type::dim2,
                            specfem::element::medium_tag::acoustic,
@@ -339,14 +438,12 @@ void check_point_properties(
 
 template <bool using_simd>
 void check_point_properties(
-    const specfem::point::properties<specfem::dimension::type::dim2,
-                                     specfem::element::medium_tag::elastic_sv,
-                                     specfem::element::property_tag::isotropic,
-                                     using_simd> &p1,
-    const specfem::point::properties<specfem::dimension::type::dim2,
-                                     specfem::element::medium_tag::elastic_sv,
-                                     specfem::element::property_tag::isotropic,
-                                     using_simd> &p2,
+    const specfem::point::properties<
+        specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+        specfem::element::property_tag::isotropic, using_simd> &p1,
+    const specfem::point::properties<
+        specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
+        specfem::element::property_tag::isotropic, using_simd> &p2,
     const int &n_simd_elements) {
   check_eq<using_simd>(p1.rho, p2.rho, n_simd_elements);
   check_eq<using_simd>(p1.mu, p2.mu, n_simd_elements);
@@ -366,12 +463,10 @@ void check_point_properties(
 template <bool using_simd>
 void check_point_properties(
     const specfem::point::properties<
-        specfem::dimension::type::dim2,
-        specfem::element::medium_tag::elastic_sv,
+        specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
         specfem::element::property_tag::anisotropic, using_simd> &p1,
     const specfem::point::properties<
-        specfem::dimension::type::dim2,
-        specfem::element::medium_tag::elastic_sv,
+        specfem::dimension::type::dim2, specfem::element::medium_tag::elastic,
         specfem::element::property_tag::anisotropic, using_simd> &p2,
     const int &n_simd_elements) {
   check_eq<using_simd>(p1.rho, p2.rho, n_simd_elements);
@@ -747,9 +842,9 @@ void test_properties(
 
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       TEST_STORE_AND_LOAD,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef TEST_STORE_AND_LOAD
 
@@ -761,9 +856,9 @@ void test_properties(
   // stage 5: check if properties are correctly written and read
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(
       TEST_COMPUTE_TO_MESH,
-      WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
-              WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
+      WHERE(DIMENSION_TAG_DIM2) WHERE(
+          MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH, MEDIUM_TAG_ACOUSTIC)
+          WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC))
 
 #undef TEST_COMPUTE_TO_MESH
   // check_compute_to_mesh<specfem::element::medium_tag::elastic_sv,
diff --git a/tests/unit-tests/assembly/test_config.yaml b/tests/unit-tests/assembly/test_config.yaml
index 60d620bd..3c4853e4 100644
--- a/tests/unit-tests/assembly/test_config.yaml
+++ b/tests/unit-tests/assembly/test_config.yaml
@@ -1,51 +1,89 @@
 Tests:
-  # Test 1:
-  - name : "Test 1: Simple mesh with flat topography"
+
+  - name : "Test 1: Simple mesh with flat topography (P_SV wave)"
     description: >
       Testing mesh reader on a simple mesh with flat topography.
     config:
       nproc : 1
+      elastic_wave: "P_SV"
     databases:
       mesh : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/database.bin"
       sources : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/sources.yaml"
       stations : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/STATIONS"
-  # Test 2:
-  - name : "Test 2: Simple mesh with curved topography"
+
+  - name : "Test 2: Simple mesh with flat topography (SH wave)"
+    description: >
+      Testing mesh reader on a simple mesh with flat topography.
+    config:
+      nproc : 1
+      elastic_wave: "SH"
+    databases:
+      mesh : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/database.bin"
+      sources : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/sources.yaml"
+      stations : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/STATIONS"
+
+  - name : "Test 3: Simple mesh with curved topography"
     description: >
       Testing mesh reader on a simple mesh with curved topography.
     config:
       nproc : 1
+      elastic_wave: "P_SV"
     databases:
       mesh : "../../../tests/unit-tests/data/mesh/simple_mesh_curved_topography/database.bin"
       sources : "../../../tests/unit-tests/data/mesh/simple_mesh_curved_topography/sources.yaml"
       stations : "../../../tests/unit-tests/data/mesh/simple_mesh_curved_topography/STATIONS"
-  # Test 3:
-  - name : "Test 3: Simple mesh with flat ocean bottom"
+
+  - name : "Test 4: Simple mesh with flat ocean bottom"
     description: >
       Testing mesh reader on a simple mesh with flat ocean bottom. Fluid-solid interface is at the ocean bottom.
     config:
       nproc : 1
+      elastic_wave: "P_SV"
     databases:
       mesh : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_flat_ocean_bottom/database.bin"
       sources : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_flat_ocean_bottom/sources.yaml"
       stations : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_flat_ocean_bottom/STATIONS"
-  # Test 4:
-  - name : "Test 4: Simple mesh with curved ocean bottom"
+
+  - name : "Test 5: Simple mesh with curved ocean bottom"
     description: >
       Testing mesh reader on a simple mesh with curved ocean bottom. Fluid-solid interface is at the ocean bottom.
     config:
       nproc : 1
+      elastic_wave: "P_SV"
     databases:
       mesh : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_curved_ocean_bottom/database.bin"
       sources : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_curved_ocean_bottom/sources.yaml"
       stations : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_curved_ocean_bottom/STATIONS"
-  # Test 5:
-  - name : "Test 5: Gmesh Example"
+
+  - name : "Test 6: Gmesh Example"
     description: >
       Testing mesh reader on a Gmesh example mesh.
     config:
       nproc : 1
+      elastic_wave: "P_SV"
     databases:
       mesh : "../../../tests/unit-tests/data/mesh/Gmesh_Example_Stacey/database.bin"
       sources : "../../../tests/unit-tests/data/mesh/Gmesh_Example_Stacey/sources.yaml"
       stations : "../../../tests/unit-tests/data/mesh/Gmesh_Example_Stacey/STATIONS"
+
+  - name : "Test 7: Homogeneous Elastic Anisotropic Material (P_SV wave)"
+    description: >
+      Testing the reading of a homogeneous elastic anisotropic material.
+    config:
+      nproc : 1
+      elastic_wave: "P_SV"
+    databases:
+      mesh : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/database.bin"
+      sources : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/sources.yaml"
+      stations : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/STATIONS"
+
+  - name : "Test 8: Homogeneous Elastic Anisotropic Material (SH wave)"
+    description: >
+      Testing the reading of a homogeneous elastic anisotropic material.
+    config:
+      nproc : 1
+      elastic_wave: "SH"
+    databases:
+      mesh : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/database.bin"
+      sources : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/sources.yaml"
+      stations : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/STATIONS"
diff --git a/tests/unit-tests/assembly/test_fixture/test_fixture.cpp b/tests/unit-tests/assembly/test_fixture/test_fixture.cpp
index 4ff4889c..401955f2 100644
--- a/tests/unit-tests/assembly/test_fixture/test_fixture.cpp
+++ b/tests/unit-tests/assembly/test_fixture/test_fixture.cpp
@@ -32,8 +32,9 @@ ASSEMBLY::ASSEMBLY() {
   for (auto &Test : Tests) {
     const auto [database_file, sources_file, stations_file] =
         Test.get_databases();
-    const auto mesh = specfem::IO::read_mesh(
-        database_file, specfem::enums::elastic_wave::p_sv, mpi);
+
+    const auto wave = Test.get_elastic_wave();
+    const auto mesh = specfem::IO::read_mesh(database_file, wave, mpi);
 
     this->Meshes.push_back(mesh);
 
diff --git a/tests/unit-tests/assembly/test_fixture/test_fixture.hpp b/tests/unit-tests/assembly/test_fixture/test_fixture.hpp
index f8a8caef..6ad1c5be 100644
--- a/tests/unit-tests/assembly/test_fixture/test_fixture.hpp
+++ b/tests/unit-tests/assembly/test_fixture/test_fixture.hpp
@@ -40,18 +40,48 @@ struct database {
   std::string stations;
 };
 
+struct config {
+public:
+  config() : nproc(1), elastic_wave("P_SV"){};
+  config(const YAML::Node &Node) {
+    nproc = Node["nproc"].as<int>();
+    elastic_wave = Node["elastic_wave"].as<std::string>();
+  }
+
+  int get_nproc() { return nproc; }
+  specfem::enums::elastic_wave get_elastic_wave() {
+    if (elastic_wave == "P_SV")
+      return specfem::enums::elastic_wave::p_sv;
+    else if (elastic_wave == "SH")
+      return specfem::enums::elastic_wave::sh;
+    else
+      throw std::runtime_error("Elastic wave type not supported");
+  }
+
+private:
+  int nproc;
+  std::string elastic_wave;
+};
+
 struct Test {
 public:
   Test(const YAML::Node &Node) {
     name = Node["name"].as<std::string>();
     description = Node["description"].as<std::string>();
     YAML::Node databases = Node["databases"];
+    YAML::Node configuration = Node["config"];
     try {
       database = test_configuration::database(databases);
     } catch (std::runtime_error &e) {
       throw std::runtime_error("Error in test configuration: " + name + "\n" +
                                e.what());
     }
+    try {
+      config = test_configuration::config(configuration);
+    } catch (std::runtime_error &e) {
+      throw std::runtime_error("Error in test configuration: " + name + "\n" +
+                               e.what());
+    }
     return;
   }
 
@@ -59,9 +89,16 @@ struct Test {
     return database.get_databases();
   }
 
+  int get_nproc() { return config.get_nproc(); }
+
+  specfem::enums::elastic_wave get_elastic_wave() {
+    return config.get_elastic_wave();
+  }
+
   std::string name;
   std::string description;
   test_configuration::database database;
+  test_configuration::config config;
 };
 } // namespace test_configuration
 
diff --git a/tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/sources.yaml b/tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/sources.yaml
new file mode 100644
index 00000000..2a674e3d
--- /dev/null
+++ b/tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/sources.yaml
@@ -0,0 +1,13 @@
+number-of-sources: 1
+sources:
+  - force:
+      x : 2500.0
+      z : 2500.0
+      source_surf: false
+      angle : 0.0
+      vx : 0.0
+      vz : 0.0
+      Ricker:
+        factor: 1e10
+        tshift: 0.0
+        f0: 10.0