diff --git a/unit/src/TestSubdomains.C b/unit/src/TestSubdomains.C
new file mode 100644
index 00000000..95a7a2c3
--- /dev/null
+++ b/unit/src/TestSubdomains.C
@@ -0,0 +1,122 @@
+#include "coefficients.h"
+#include "mfem.hpp"
+#include "gtest/gtest.h"
+
+double
+scalar_f(const mfem::Vector & x)
+{
+  return x.Sum();
+}
+
+void
+vector_f(const mfem::Vector & x, mfem::Vector & f)
+{
+  for (int i = 0; i < x.Size(); ++i)
+  {
+    f[i] = sin(x[i]);
+  }
+}
+
+TEST(CheckSetup, PWCoefficient)
+{
+  mfem::Vector one_vec(3);
+  one_vec = 1.0;
+
+  platypus::Subdomain air("air", 1);
+  air._scalar_coefficients.Register("electrical_conductivity",
+                                    std::make_shared<mfem::ConstantCoefficient>(1e-4));
+  air._scalar_coefficients.Register("magnetic_permeability",
+                                    std::make_shared<mfem::ConstantCoefficient>(M_PI * 4.0e-7));
+  air._vector_coefficients.Register("vector_function",
+                                    std::make_shared<mfem::VectorConstantCoefficient>(one_vec));
+
+  platypus::Subdomain plate("plate", 2);
+  plate._scalar_coefficients.Register("electrical_conductivity",
+                                      std::make_shared<mfem::FunctionCoefficient>(scalar_f));
+  plate._scalar_coefficients.Register("magnetic_permeability",
+                                      std::make_shared<mfem::ConstantCoefficient>(M_PI * 4.0e-7));
+  plate._vector_coefficients.Register(
+      "vector_function", std::make_shared<mfem::VectorFunctionCoefficient>(3, vector_f));
+
+  platypus::Subdomain coil1("coil1", 3);
+  coil1._scalar_coefficients.Register("electrical_conductivity",
+                                      std::make_shared<mfem::ConstantCoefficient>(1e4));
+  coil1._scalar_coefficients.Register("magnetic_permeability",
+                                      std::make_shared<mfem::ConstantCoefficient>(M_PI * 4.0e-7));
+  coil1._vector_coefficients.Register(
+      "vector_function", std::make_shared<mfem::VectorFunctionCoefficient>(3, vector_f));
+
+  platypus::Subdomain coil2("coil2", 4);
+  coil2._scalar_coefficients.Register("electrical_conductivity",
+                                      std::make_shared<mfem::ConstantCoefficient>(1e4));
+  coil2._scalar_coefficients.Register("magnetic_permeability",
+                                      std::make_shared<mfem::ConstantCoefficient>(M_PI * 4.0e-7));
+  coil2._vector_coefficients.Register(
+      "vector_function", std::make_shared<mfem::VectorFunctionCoefficient>(3, vector_f));
+
+  platypus::Subdomain coil3("coil3", 5);
+  coil3._scalar_coefficients.Register("electrical_conductivity",
+                                      std::make_shared<mfem::ConstantCoefficient>(1e4));
+  coil3._scalar_coefficients.Register("magnetic_permeability",
+                                      std::make_shared<mfem::ConstantCoefficient>(M_PI * 4.0e-7));
+  coil3._vector_coefficients.Register(
+      "vector_function", std::make_shared<mfem::VectorFunctionCoefficient>(3, vector_f));
+
+  platypus::Subdomain coil4("coil4", 6);
+  coil4._scalar_coefficients.Register("electrical_conductivity",
+                                      std::make_shared<mfem::ConstantCoefficient>(1e4));
+  coil4._scalar_coefficients.Register("magnetic_permeability",
+                                      std::make_shared<mfem::ConstantCoefficient>(M_PI * 4.0e-7));
+  coil4._vector_coefficients.Register(
+      "vector_function", std::make_shared<mfem::VectorFunctionCoefficient>(3, vector_f));
+
+  platypus::Coefficients coefficients(
+      std::vector<platypus::Subdomain>({air, plate, coil1, coil2, coil3, coil4}));
+
+  coefficients._scalars.Register("I", std::make_shared<mfem::ConstantCoefficient>(1.0));
+  coefficients._vectors.Register("const_vector",
+                                 std::make_shared<mfem::VectorConstantCoefficient>(one_vec));
+
+  EXPECT_EQ(coefficients._subdomains.size(), 6);
+
+  EXPECT_TRUE(coefficients._scalars.Has("electrical_conductivity"));
+  EXPECT_TRUE(coefficients._scalars.Has("magnetic_permeability"));
+  EXPECT_TRUE(coefficients._scalars.Has("I"));
+
+  EXPECT_TRUE(coefficients._vectors.Has("vector_function"));
+  EXPECT_TRUE(coefficients._vectors.Has("const_vector"));
+}
+
+TEST(CheckData, Coefficients)
+{
+
+  platypus::Subdomain wire("wire", 1);
+  wire._scalar_coefficients.Register("property_one",
+                                     std::make_shared<mfem::ConstantCoefficient>(1.0));
+  wire._scalar_coefficients.Register("property_two",
+                                     std::make_shared<mfem::ConstantCoefficient>(150.0));
+
+  platypus::Subdomain air("air", 2);
+  air._scalar_coefficients.Register("property_one",
+                                    std::make_shared<mfem::ConstantCoefficient>(26.0));
+  air._scalar_coefficients.Register("property_two",
+                                    std::make_shared<mfem::ConstantCoefficient>(152.0));
+
+  platypus::Coefficients coefficients(std::vector<platypus::Subdomain>({wire, air}));
+
+  // Verify predefined values
+  mfem::IsoparametricTransformation t;
+  mfem::IntegrationPoint ip;
+
+  mfem::Coefficient * pw = coefficients._scalars.Get("property_one");
+  t.Attribute = 1;
+  EXPECT_DOUBLE_EQ(pw->Eval(t, ip), 1.0);
+  t.Attribute = 2;
+  EXPECT_DOUBLE_EQ(pw->Eval(t, ip), 26.0);
+
+  pw = coefficients._scalars.Get("property_two");
+  t.Attribute = 1;
+  EXPECT_DOUBLE_EQ(pw->Eval(t, ip), 150.0);
+  t.Attribute = 2;
+  EXPECT_DOUBLE_EQ(pw->Eval(t, ip), 152.0);
+}
\ No newline at end of file
diff --git a/unit/src/TestVariables.C b/unit/src/TestVariables.C
deleted file mode 100644
index ceb5115e..00000000
--- a/unit/src/TestVariables.C
+++ /dev/null
@@ -1,24 +0,0 @@
-#include "problem_builder.h"
-#include "mfem.hpp"
-#include "gtest/gtest.h"
-
-TEST(CheckSetup, variables)
-{
-  mfem::Mesh mesh("../../data/beam-tet.mesh", 1, 1);
-
-  auto problem_builder = std::make_unique<platypus::TimeDomainProblemBuilder>();
-
-  auto pmesh = std::make_shared<mfem::ParMesh>(MPI_COMM_WORLD, mesh);
-
-  problem_builder->SetMesh(pmesh);
-  problem_builder->AddFESpace(std::string("HCurl"), std::string("ND_3D_P2"));
-  problem_builder->AddGridFunction(std::string("vector_potential"), std::string("HCurl"));
-
-  auto problem = problem_builder->ReturnProblem();
-
-  auto stored_gf = problem->_gridfunctions.Get("vector_potential");
-  auto stored_fespace = problem->_fespaces.Get("HCurl");
-
-  EXPECT_EQ(stored_fespace->GetVSize(), stored_gf->ParFESpace()->GetVSize());
-  EXPECT_GT(stored_fespace->GetVSize(), 0);
-}