diff --git a/src/Elliptic/BoundaryConditions/AnalyticSolution.hpp b/src/Elliptic/BoundaryConditions/AnalyticSolution.hpp
index 81e58d3945ee..8d9958fa63bd 100644
--- a/src/Elliptic/BoundaryConditions/AnalyticSolution.hpp
+++ b/src/Elliptic/BoundaryConditions/AnalyticSolution.hpp
@@ -4,16 +4,13 @@
 #pragma once
 
 #include <cstddef>
+#include <memory>
 #include <ostream>
 #include <pup.h>
 #include <string>
 #include <vector>
 
-#include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
-#include "DataStructures/Tensor/Slice.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
-#include "Domain/Structure/Direction.hpp"
-#include "Domain/Structure/IndexToSliceAt.hpp"
 #include "Domain/Tags.hpp"
 #include "Domain/Tags/FaceNormal.hpp"
 #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
@@ -21,13 +18,14 @@
 #include "Elliptic/BoundaryConditions/Tags.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
 #include "Options/String.hpp"
-#include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
+#include "Parallel/Tags/Metavariables.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
+#include "Utilities/CallWithDynamicType.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
-#include "Utilities/TypeTraits/IsA.hpp"
 
 namespace elliptic::BoundaryConditions {
 
@@ -61,14 +59,29 @@ class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
   using Base = BoundaryCondition<Dim>;
 
  public:
+  struct Solution {
+    using type = std::unique_ptr<elliptic::analytic_data::AnalyticSolution>;
+    static constexpr Options::String help = {
+        "The analytic solution to impose on the boundary"};
+  };
+
   using options =
-      tmpl::list<elliptic::OptionTags::BoundaryConditionType<FieldTags>...>;
+      tmpl::list<Solution,
+                 elliptic::OptionTags::BoundaryConditionType<FieldTags>...>;
   static constexpr Options::String help =
       "Boundary conditions from the analytic solution";
 
   AnalyticSolution() = default;
-  AnalyticSolution(const AnalyticSolution&) = default;
-  AnalyticSolution& operator=(const AnalyticSolution&) = default;
+  AnalyticSolution(const AnalyticSolution& rhs) : Base(rhs) { *this = rhs; }
+  AnalyticSolution& operator=(const AnalyticSolution& rhs) {
+    if (rhs.solution_ != nullptr) {
+      solution_ = rhs.solution_->get_clone();
+    } else {
+      solution_ = nullptr;
+    }
+    boundary_condition_types_ = rhs.boundary_condition_types_;
+    return *this;
+  }
   AnalyticSolution(AnalyticSolution&&) = default;
   AnalyticSolution& operator=(AnalyticSolution&&) = default;
   ~AnalyticSolution() = default;
@@ -81,11 +94,13 @@ class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
 
   /// Select which `elliptic::BoundaryConditionType` to apply for each field
   explicit AnalyticSolution(
+      std::unique_ptr<elliptic::analytic_data::AnalyticSolution> solution,
       // This pack expansion repeats the type `elliptic::BoundaryConditionType`
       // for each system field
       const typename elliptic::OptionTags::BoundaryConditionType<
           FieldTags>::type... boundary_condition_types)
-      : boundary_condition_types_{boundary_condition_types...} {}
+      : solution_(std::move(solution)),
+        boundary_condition_types_{boundary_condition_types...} {}
 
   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
       const override {
@@ -110,41 +125,35 @@ class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
   }
 
   using argument_tags =
-      tmpl::list<::Tags::AnalyticSolutionsBase, domain::Tags::Mesh<Dim>,
-                 domain::Tags::Direction<Dim>,
+      tmpl::list<Parallel::Tags::Metavariables,
+                 domain::Tags::Coordinates<Dim, Frame::Inertial>,
                  ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<
                      Dim, Frame::Inertial>>>;
-  using volume_tags =
-      tmpl::list<::Tags::AnalyticSolutionsBase, domain::Tags::Mesh<Dim>>;
+  using volume_tags = tmpl::list<Parallel::Tags::Metavariables>;
 
-  template <typename OptionalAnalyticSolutions>
+  template <typename Metavariables>
   void apply(const gsl::not_null<typename FieldTags::type*>... fields,
              const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
-             const OptionalAnalyticSolutions& optional_analytic_solutions,
-             const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
+             const Metavariables& /*meta*/,
+             const tnsr::I<DataVector, Dim>& face_inertial_coords,
              const tnsr::i<DataVector, Dim>& face_normal) const {
-    const auto& analytic_solutions = [&optional_analytic_solutions]()
-        -> const auto& {
-      if constexpr (tt::is_a_v<std::optional, OptionalAnalyticSolutions>) {
-        if (not optional_analytic_solutions.has_value()) {
-          ERROR_NO_TRACE(
-              "Trying to impose boundary conditions from an analytic solution, "
-              "but no analytic solution is available. You probably selected "
-              "the 'AnalyticSolution' boundary condition but chose to solve a "
-              "problem that has no analytic solution. If this is the case, you "
-              "should probably select a different boundary condition.");
-        }
-        return *optional_analytic_solutions;
-      } else {
-        return optional_analytic_solutions;
-      }
-    }
-    ();
-    const size_t slice_index =
-        index_to_slice_at(volume_mesh.extents(), direction);
-    const auto impose_boundary_condition = [this, &analytic_solutions,
-                                            &volume_mesh, &direction,
-                                            &slice_index, &face_normal](
+    // Retrieve variables for both Dirichlet and Neumann conditions, then decide
+    // which to impose. We could also retrieve either the field for the flux for
+    // each field individually based on the selection, but that would incur the
+    // overhead of calling into the analytic solution multiple times and
+    // possibly computing temporary quantities multiple times. This performance
+    // consideration is probably irrelevant because the boundary conditions are
+    // only evaluated once at the beginning of the solve.
+    using analytic_tags = tmpl::list<FieldTags..., FluxTags...>;
+    using factory_classes =
+        typename Metavariables::factory_creation::factory_classes;
+    const auto solution_vars = call_with_dynamic_type<
+        tuples::tagged_tuple_from_typelist<analytic_tags>,
+        tmpl::at<factory_classes, elliptic::analytic_data::AnalyticSolution>>(
+        solution_.get(), [&face_inertial_coords](const auto* const derived) {
+          return derived->variables(face_inertial_coords, analytic_tags{});
+        });
+    const auto impose_boundary_condition = [this, &solution_vars, &face_normal](
                                                auto field_tag_v,
                                                auto flux_tag_v,
                                                const auto field,
@@ -154,18 +163,11 @@ class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
       switch (get<elliptic::Tags::BoundaryConditionType<field_tag>>(
           boundary_condition_types_)) {
         case elliptic::BoundaryConditionType::Dirichlet:
-          data_on_slice(
-              field,
-              get<::Tags::detail::AnalyticImpl<field_tag>>(analytic_solutions),
-              volume_mesh.extents(), direction.dimension(), slice_index);
+          *field = get<field_tag>(solution_vars);
           break;
         case elliptic::BoundaryConditionType::Neumann:
-          normal_dot_flux(
-              n_dot_flux, face_normal,
-              data_on_slice(get<::Tags::detail::AnalyticImpl<flux_tag>>(
-                                analytic_solutions),
-                            volume_mesh.extents(), direction.dimension(),
-                            slice_index));
+          normal_dot_flux(n_dot_flux, face_normal,
+                          get<flux_tag>(solution_vars));
           break;
         default:
           ERROR("Unsupported boundary condition type: "
@@ -213,16 +215,7 @@ class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
   void pup(PUP::er& p) override;
 
  private:
-  friend bool operator==(const AnalyticSolution& lhs,
-                         const AnalyticSolution& rhs) {
-    return lhs.boundary_condition_types_ == rhs.boundary_condition_types_;
-  }
-
-  friend bool operator!=(const AnalyticSolution& lhs,
-                         const AnalyticSolution& rhs) {
-    return not(lhs == rhs);
-  }
-
+  std::unique_ptr<elliptic::analytic_data::AnalyticSolution> solution_{nullptr};
   tuples::TaggedTuple<elliptic::Tags::BoundaryConditionType<FieldTags>...>
       boundary_condition_types_{};
 };
@@ -232,6 +225,7 @@ template <typename System, size_t Dim, typename... FieldTags,
 void AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
                       tmpl::list<FluxTags...>>::pup(PUP::er& p) {
   Base::pup(p);
+  p | solution_;
   p | boundary_condition_types_;
 }
 
diff --git a/tests/InputFiles/Elasticity/BentBeam.yaml b/tests/InputFiles/Elasticity/BentBeam.yaml
index 15eb18367d6c..0a4462a66b85 100644
--- a/tests/InputFiles/Elasticity/BentBeam.yaml
+++ b/tests/InputFiles/Elasticity/BentBeam.yaml
@@ -20,7 +20,7 @@ ResourceInfo:
   AvoidGlobalProc0: false
   Singletons: Auto
 
-Background:
+Background: &solution
   BentBeam:
     Length: 2.
     Height: 1.
@@ -44,6 +44,7 @@ DomainCreator:
     TimeDependence: None
     BoundaryCondition:
       AnalyticSolution:
+        Solution: *solution
         Displacement: Dirichlet
 
 Discretization:
diff --git a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml
index 075662f9e791..7560e13c9ff5 100644
--- a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml
+++ b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml
@@ -22,7 +22,7 @@ ResourceInfo:
   AvoidGlobalProc0: false
   Singletons: Auto
 
-Background:
+Background: &solution
   HalfSpaceMirror:
     BeamWidth: 0.177
     Material: &fused_silica
@@ -54,12 +54,15 @@ DomainCreator:
     BoundaryConditions:
       LowerZ:
         AnalyticSolution:
+          Solution: *solution
           Displacement: Dirichlet
       UpperZ:
         AnalyticSolution:
+          Solution: *solution
           Displacement: Dirichlet
       Mantle:
         AnalyticSolution:
+          Solution: *solution
           Displacement: Dirichlet
 
 Discretization:
diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml
index 9cd4acfbaae8..fc074aef9d80 100644
--- a/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml
+++ b/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml
@@ -23,7 +23,7 @@ ResourceInfo:
   AvoidGlobalProc0: false
   Singletons: Auto
 
-Background:
+Background: &solution
   ProductOfSinusoids:
     WaveNumbers: [1]
 
@@ -44,9 +44,11 @@ DomainCreator:
     BoundaryConditions:
       LowerBoundary:
         AnalyticSolution:
+          Solution: *solution
           Field: Dirichlet
       UpperBoundary:
         AnalyticSolution:
+          Solution: *solution
           Field: Neumann
 
 Discretization:
diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml
index 36c06bb76b6b..8f6b3c97a663 100644
--- a/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml
+++ b/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml
@@ -21,7 +21,7 @@ ResourceInfo:
   AvoidGlobalProc0: false
   Singletons: Auto
 
-Background:
+Background: &solution
   ProductOfSinusoids:
     WaveNumbers: [1, 1]
 
@@ -39,6 +39,7 @@ DomainCreator:
     TimeDependence: None
     BoundaryCondition:
       AnalyticSolution:
+        Solution: *solution
         Field: Dirichlet
 
 Discretization:
diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml
index f9dc2ca7c269..4ac8c31956bd 100644
--- a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml
+++ b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml
@@ -22,7 +22,7 @@ ResourceInfo:
   AvoidGlobalProc0: false
   Singletons: Auto
 
-Background:
+Background: &solution
   ProductOfSinusoids:
     WaveNumbers: [1, 1, 1]
 
@@ -40,12 +40,15 @@ DomainCreator:
     TimeDependence: None
     BoundaryConditionInX:
       AnalyticSolution:
+        Solution: *solution
         Field: Dirichlet
     BoundaryConditionInY:
       AnalyticSolution:
+        Solution: *solution
         Field: Dirichlet
     BoundaryConditionInZ:
       AnalyticSolution:
+        Solution: *solution
         Field: Dirichlet
 
 Discretization:
diff --git a/tests/InputFiles/Xcts/KerrSchild.yaml b/tests/InputFiles/Xcts/KerrSchild.yaml
index 229230f9fa4e..90ec788c0b37 100644
--- a/tests/InputFiles/Xcts/KerrSchild.yaml
+++ b/tests/InputFiles/Xcts/KerrSchild.yaml
@@ -22,7 +22,7 @@ ResourceInfo:
   AvoidGlobalProc0: false
   Singletons: Auto
 
-Background:
+Background: &solution
   KerrSchild:
     Mass: 1.
     Spin: [0., 0., 0.]
@@ -37,6 +37,7 @@ DomainCreator:
     Interior:
       ExciseWithBoundaryCondition:
         AnalyticSolution:
+          Solution: *solution
           ConformalFactor: Dirichlet
           LapseTimesConformalFactor: Dirichlet
           ShiftExcess: Dirichlet
@@ -50,6 +51,7 @@ DomainCreator:
     TimeDependentMaps: None
     OuterBoundaryCondition:
       AnalyticSolution:
+        Solution: *solution
         ConformalFactor: Dirichlet
         LapseTimesConformalFactor: Dirichlet
         ShiftExcess: Dirichlet
diff --git a/tests/InputFiles/Xcts/TovStar.yaml b/tests/InputFiles/Xcts/TovStar.yaml
index 53142b22e528..f34034f6cd94 100644
--- a/tests/InputFiles/Xcts/TovStar.yaml
+++ b/tests/InputFiles/Xcts/TovStar.yaml
@@ -22,7 +22,7 @@ ResourceInfo:
   AvoidGlobalProc0: false
   Singletons: Auto
 
-Background:
+Background: &solution
   TovStar:
     CentralDensity: 1.e-3
     EquationOfState:
@@ -49,6 +49,7 @@ DomainCreator:
     TimeDependentMaps: None
     OuterBoundaryCondition:
       AnalyticSolution:
+        Solution: *solution
         ConformalFactor: Dirichlet
         LapseTimesConformalFactor: Dirichlet
         ShiftExcess: Dirichlet
diff --git a/tests/Unit/Elliptic/BoundaryConditions/Test_AnalyticSolution.cpp b/tests/Unit/Elliptic/BoundaryConditions/Test_AnalyticSolution.cpp
index 2e9c4e182843..410b221688fe 100644
--- a/tests/Unit/Elliptic/BoundaryConditions/Test_AnalyticSolution.cpp
+++ b/tests/Unit/Elliptic/BoundaryConditions/Test_AnalyticSolution.cpp
@@ -14,7 +14,6 @@
 #include "DataStructures/DataBox/DataBoxTag.hpp"
 #include "DataStructures/DataBox/Prefixes.hpp"
 #include "DataStructures/DataVector.hpp"
-#include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
 #include "Domain/FaceNormal.hpp"
@@ -32,8 +31,11 @@
 #include "NumericalAlgorithms/Spectral/Basis.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
+#include "Options/Protocols/FactoryCreation.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
 #include "Utilities/Literals.hpp"
+#include "Utilities/ProtocolHelpers.hpp"
+#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
 #include "Utilities/TMPL.hpp"
 
 namespace elliptic::BoundaryConditions {
@@ -58,27 +60,72 @@ struct System {
   using primal_fluxes = tmpl::list<FluxTag<Dim, 1>, FluxTag<Dim, 2>>;
 };
 
+template <size_t Dim>
+using test_tags = tmpl::list<ScalarFieldTag<1>, ScalarFieldTag<2>,
+                             FluxTag<Dim, 1>, FluxTag<Dim, 2>>;
+
+template <size_t Dim>
+struct TestSolution : elliptic::analytic_data::AnalyticSolution {
+  using options = tmpl::list<>;
+  static constexpr Options::String help{"A solution."};
+  TestSolution() = default;
+  TestSolution(const TestSolution&) = default;
+  TestSolution& operator=(const TestSolution&) = default;
+  TestSolution(TestSolution&&) = default;
+  TestSolution& operator=(TestSolution&&) = default;
+  ~TestSolution() override = default;
+  explicit TestSolution(CkMigrateMessage* m)
+      : elliptic::analytic_data::AnalyticSolution(m) {}
+  using PUP::able::register_constructor;
+  WRAPPED_PUPable_decl_template(TestSolution);  // NOLINT
+
+  std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
+      const override {
+    return std::make_unique<TestSolution>(*this);
+  }
+
+  tuples::tagged_tuple_from_typelist<test_tags<Dim>> variables(
+      const tnsr::I<DataVector, Dim>& x, test_tags<Dim> /*meta*/) const {
+    // Create arbitrary analytic solution data
+    Variables<test_tags<Dim>> result{x.begin()->size()};
+    std::iota(result.data(), result.data() + result.size(), 1.);
+    return {get<ScalarFieldTag<1>>(result), get<ScalarFieldTag<2>>(result),
+            get<FluxTag<Dim, 1>>(result), get<FluxTag<Dim, 2>>(result)};
+  }
+};
+
+template <size_t Dim>
+PUP::able::PUP_ID TestSolution<Dim>::my_PUP_ID = 0;  // NOLINT
+
+template <size_t Dim>
+struct Metavariables {
+  struct factory_creation
+      : tt::ConformsTo<Options::protocols::FactoryCreation> {
+    using factory_classes = tmpl::map<
+        tmpl::pair<elliptic::BoundaryConditions::BoundaryCondition<Dim>,
+                   tmpl::list<elliptic::BoundaryConditions::AnalyticSolution<
+                       System<Dim>>>>,
+        tmpl::pair<elliptic::analytic_data::AnalyticSolution,
+                   tmpl::list<TestSolution<Dim>>>>;
+  };
+};
+
 template <size_t Dim>
 void test_analytic_solution() {
   CAPTURE(Dim);
   // Test factory-creation
-  const auto created =
-      TestHelpers::test_factory_creation<BoundaryCondition<Dim>,
-                                         AnalyticSolution<System<Dim>>>(
+  register_factory_classes_with_charm<Metavariables<Dim>>();
+  const auto created = serialize_and_deserialize(
+      TestHelpers::test_creation<std::unique_ptr<BoundaryCondition<Dim>>,
+                                 Metavariables<Dim>>(
           "AnalyticSolution:\n"
+          "  Solution: TestSolution\n"
           "  Field1: Dirichlet\n"
-          "  Field2: Neumann");
+          "  Field2: Neumann"));
   REQUIRE(dynamic_cast<const AnalyticSolution<System<Dim>>*>(created.get()) !=
           nullptr);
   const auto& boundary_condition =
       dynamic_cast<const AnalyticSolution<System<Dim>>&>(*created);
-  {
-    INFO("Semantics");
-    test_serialization(boundary_condition);
-    test_copy_semantics(boundary_condition);
-    auto move_boundary_condition = boundary_condition;
-    test_move_semantics(std::move(move_boundary_condition), boundary_condition);
-  }
   {
     INFO("Properties");
     CHECK(boundary_condition.boundary_condition_types() ==
@@ -90,29 +137,20 @@ void test_analytic_solution() {
     INFO("Test applying the boundary conditions");
     const Mesh<Dim> volume_mesh{3, Spectral::Basis::Legendre,
                                 Spectral::Quadrature::GaussLobatto};
-    const size_t volume_num_points = volume_mesh.number_of_grid_points();
     const auto direction = Direction<Dim>::lower_xi();
     const auto face_mesh = volume_mesh.slice_away(direction.dimension());
     const size_t face_num_points = face_mesh.number_of_grid_points();
+    tnsr::I<DataVector, Dim> face_inertial_coords{face_num_points, 0.};
     tnsr::i<DataVector, Dim> face_normal{face_num_points, 0.};
     get<0>(face_normal) = -2.;
-    // Create arbitrary analytic solution data
-    Variables<tmpl::list<::Tags::detail::AnalyticImpl<ScalarFieldTag<1>>,
-                         ::Tags::detail::AnalyticImpl<ScalarFieldTag<2>>,
-                         ::Tags::detail::AnalyticImpl<FluxTag<Dim, 1>>,
-                         ::Tags::detail::AnalyticImpl<FluxTag<Dim, 2>>>>
-        analytic_solutions{volume_num_points};
-    std::iota(analytic_solutions.data(),
-              analytic_solutions.data() + analytic_solutions.size(), 1.);
     const auto box = db::create<db::AddSimpleTags<
-        domain::Tags::Mesh<Dim>,
-        ::Tags::AnalyticSolutions<
-            tmpl::list<ScalarFieldTag<1>, ScalarFieldTag<2>, FluxTag<Dim, 1>,
-                       FluxTag<Dim, 2>>>,
-        domain::Tags::Faces<Dim, domain::Tags::Direction<Dim>>,
+        Parallel::Tags::MetavariablesImpl<Metavariables<Dim>>,
+        domain::Tags::Faces<Dim,
+                            domain::Tags::Coordinates<Dim, Frame::Inertial>>,
         domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>>(
-        volume_mesh, std::make_optional(std::move(analytic_solutions)),
-        DirectionMap<Dim, Direction<Dim>>{{direction, direction}},
+        Metavariables<Dim>{},
+        DirectionMap<Dim, tnsr::I<DataVector, Dim>>{
+            {direction, face_inertial_coords}},
         DirectionMap<Dim, tnsr::i<DataVector, Dim>>{{direction, face_normal}});
     Variables<tmpl::list<ScalarFieldTag<1>, ScalarFieldTag<2>,
                          ::Tags::NormalDotFlux<ScalarFieldTag<1>>,
@@ -131,9 +169,9 @@ void test_analytic_solution() {
       if constexpr (Dim == 1) {
         return {1.};
       } else if constexpr (Dim == 2) {
-        return {1., 4., 7.};
+        return {1., 2., 3.};
       } else if constexpr (Dim == 3) {
-        return {1., 4., 7., 10., 13., 16., 19., 22., 25.};
+        return {1., 2., 3., 4., 5., 6., 7., 8., 9.};
       }
     }();
     CHECK_ITERABLE_APPROX(get(get<ScalarFieldTag<1>>(vars)),
@@ -145,11 +183,11 @@ void test_analytic_solution() {
     // Imposed Neumann conditions on field 2
     const auto expected_neumann_field = []() -> DataVector {
       if constexpr (Dim == 1) {
-        return {-20.};
+        return {-8.};
       } else if constexpr (Dim == 2) {
-        return {-74., -80., -86.};
+        return {-26., -28., -30.};
       } else if constexpr (Dim == 3) {
-        return {-272., -278., -284., -290., -296., -302., -308., -314., -320.};
+        return {-92., -94., -96., -98., -100., -102., -104., -106., -108.};
       }
     }();
     CHECK_ITERABLE_APPROX(
diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp
index 2623e4811df4..9a5f5315d1e6 100644
--- a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp
+++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp
@@ -131,7 +131,7 @@ make_boundary_condition(
     const elliptic::BoundaryConditionType boundary_condition_type) {
   return std::make_unique<
       elliptic::BoundaryConditions::AnalyticSolution<System>>(
-      boundary_condition_type);
+      nullptr, boundary_condition_type);
 }
 
 // Generate some random element-centered subdomain data on each element
diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp
index 31d35dbe10a0..8a8c28127964 100644
--- a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp
+++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp
@@ -161,6 +161,8 @@ struct Metavariables {
         tmpl::pair<
             elliptic::BoundaryConditions::BoundaryCondition<System::volume_dim>,
             tmpl::list<elliptic::BoundaryConditions::AnalyticSolution<System>>>,
+        tmpl::pair<elliptic::analytic_data::AnalyticSolution,
+                   tmpl::list<AnalyticSolution>>,
         tmpl::pair<AnalyticSolution, tmpl::list<AnalyticSolution>>>;
   };
 };
@@ -443,9 +445,11 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {{4}},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Neumann)};
       const ElementId<1> left_id{0, {{{2, 0}}}};
       const ElementId<1> midleft_id{0, {{{2, 1}}}};
@@ -524,9 +528,11 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {{12}},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Neumann)};
       Approx analytic_solution_aux_approx =
           Approx::custom().epsilon(1.e-8).scale(M_PI);
@@ -562,6 +568,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {{3, 2}},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           nullptr};
       const ElementId<2> northwest_id{0, {{{1, 0}, {1, 1}}}};
@@ -627,6 +634,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {{12, 12}},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           nullptr};
       Approx analytic_solution_aux_approx =
@@ -657,12 +665,15 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {{2, 3, 4}},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           nullptr};
       const ElementId<3> self_id{0, {{{1, 0}, {1, 0}, {1, 0}}}};
@@ -751,12 +762,15 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {{12, 12, 12}},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet),
           nullptr};
       Approx analytic_solution_aux_approx =
@@ -797,6 +811,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet)};
       const ElementId<2> lowerleft_id{0, {{{0, 0}, {1, 0}}}};
       const ElementId<2> upperleft_id{0, {{{0, 0}, {1, 1}}}};
@@ -881,6 +896,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           {},
           std::make_unique<
               elliptic::BoundaryConditions::AnalyticSolution<system>>(
+              analytic_solution.get_clone(),
               elliptic::BoundaryConditionType::Dirichlet)};
       Approx analytic_solution_aux_approx =
           Approx::custom().epsilon(1.e-12).scale(M_PI);