diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
index 7e62ce2747195fc04aeffac93c38182dc90e2c2d..553b14174b1b9280b43252c8fa755a0600d49ba6 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
@@ -10,8 +10,6 @@
  */
 
 #include "BHECommonCoaxial.h"
-#include "FlowAndTemperatureControl.h"
-#include "ThermoMechanicalFlowProperties.h"
 
 namespace ProcessLib
 {
@@ -82,6 +80,42 @@ BHECommonCoaxial::pipeAdvectionVectors() const
              // grout g, Eq. 28 and Eq. 25
              {0, 0, 0}}};
 }
+
+std::array<double, BHECommonCoaxial::number_of_unknowns>
+BHECommonCoaxial::CrossSectionAreas() const
+{
+    auto S = cross_section_areas_coaxial();
+    return {S[0], S[1], S[2]};
+}
+
+std::array<double, BHECommonCoaxial::number_of_unknowns>
+BHECommonCoaxial::calcThermalResistances(double const Nu_inner_pipe,
+                                         double const Nu_annulus_pipe)
+{
+    // thermal resistances due to advective flow of refrigerant in the pipes
+    auto const R_advective = calculateAdvectiveThermalResistance(
+        _pipes.inner_pipe, _pipes.outer_pipe, refrigerant, Nu_inner_pipe,
+        Nu_annulus_pipe);
+
+    // thermal resistance due to thermal conductivity of the pipe wall material
+    auto const R_conductive = calculatePipeWallThermalResistance(
+        _pipes.inner_pipe, _pipes.outer_pipe);
+
+    // thermal resistance due to the grout transition and grout-soil exchange.
+    auto const R = calculateGroutAndGroutSoilExchangeThermalResistance(
+        _pipes.outer_pipe, grout, borehole_geometry.diameter);
+
+    // thermal resistance due to grout-soil exchange
+    double const R_gs = R.grout_soil;
+
+    // Eq. 56 and 57
+    double const R_ff = R_advective.inner_pipe_coaxial + R_advective.a_annulus +
+                        R_conductive.inner_pipe_coaxial;
+    double const R_fg =
+        R_advective.b_annulus + R_conductive.annulus + R.conductive_b;
+
+    return getThermalResistances(R_gs, R_ff, R_fg);
+}
 }  // namespace BHE
 }  // namespace HeatTransportBHE
 }  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
index 8d26285b7e9a4b0b432a3ac6ceadd9dedd5370c1..25be41e8046fd3996c4af576a7fe6c51a76372ef 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
@@ -12,6 +12,7 @@
 #include <Eigen/Eigen>
 #include "BHECommon.h"
 #include "FlowAndTemperatureControl.h"
+#include "Physics.h"
 #include "PipeConfigurationCoaxial.h"
 #include "ThermoMechanicalFlowProperties.h"
 
@@ -34,15 +35,18 @@ public:
     {
     }
 
+    static constexpr int number_of_unknowns = 3;
+    static constexpr int number_of_grout_zones = 1;
+
     double thermalResistance(int const unknown_index) const
     {
         return _thermal_resistances[unknown_index];
-    };
+    }
 
     double updateFlowRateAndTemperature(double T_out, double current_time);
 
-    static constexpr int number_of_unknowns = 3;
-    static constexpr int number_of_grout_zones = 1;
+    std::array<double, number_of_unknowns> calcThermalResistances(
+        double const Nu_inner_pipe, double const Nu_annulus_pipe);
 
     std::array<double, number_of_unknowns> pipeHeatCapacities() const;
 
@@ -53,11 +57,15 @@ public:
 
     std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
         const;
+    double cross_section_area_inner_pipe, cross_section_area_annulus,
+        cross_section_area_grout;
 
-protected:
-    virtual std::array<double, number_of_unknowns> calcThermalResistances(
-        double const Nu_inner_pipe, double const Nu_annulus_pipe) = 0;
+    std::array<double, number_of_unknowns> CrossSectionAreas() const;
 
+    virtual std::array<double, number_of_unknowns> cross_section_areas_coaxial()
+        const = 0;
+
+protected:
     void updateHeatTransferCoefficients(double const flow_rate)
     {
         auto const tm_flow_properties_annulus =
@@ -78,15 +86,26 @@ protected:
                 flow_rate);
 
         _flow_velocity_inner = tm_flow_properties.velocity;
+
         _thermal_resistances =
             calcThermalResistances(tm_flow_properties.nusselt_number,
                                    tm_flow_properties_annulus.nusselt_number);
+
+        auto const cross_section_area = calculateCrossSectionAreasCoaxial(
+            _pipes.inner_pipe, _pipes.outer_pipe, borehole_geometry.area());
+
+        cross_section_area_inner_pipe = cross_section_area.inner;
+        cross_section_area_annulus = cross_section_area.annulus;
+        cross_section_area_grout = cross_section_area.grout;
     }
 
     PipeConfigurationCoaxial const _pipes;
 
     virtual std::array<double, 2> velocities() const = 0;
 
+    virtual std::array<double, number_of_unknowns> getThermalResistances(
+        double const& R_gs, double const& R_ff, double const& R_fg) const = 0;
+
     /// Here we store the thermal resistances needed for computation of the heat
     /// exchange coefficients in the governing equations of BHE.
     /// These governing equations can be found in
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index 672115c9cd8fd62036ff45395bb539fd05418b7c..24b873f51904e0dd5cf76d9ae8cad4be88572dbd 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -140,10 +140,12 @@ private:
     PipeConfiguration1U const _pipes;
 
 public:
-    std::array<double, number_of_unknowns> const cross_section_areas = {
-        {_pipes.inlet.area(), _pipes.outlet.area(),
-         borehole_geometry.area() / 2 - _pipes.inlet.area(),
-         borehole_geometry.area() / 2 - _pipes.outlet.area()}};
+    std::array<double, number_of_unknowns> CrossSectionAreas() const
+    {
+        return {{_pipes.inlet.area(), _pipes.outlet.area(),
+                 borehole_geometry.area() / 2 - _pipes.inlet.area(),
+                 borehole_geometry.area() / 2 - _pipes.outlet.area()}};
+    }
 
 private:
     void updateHeatTransferCoefficients(double const flow_rate);
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.cpp
deleted file mode 100644
index 4b3665b546bd64752d4706db2d645dce6218ac95..0000000000000000000000000000000000000000
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "BHE_CXA.h"
-#include "ThermoMechanicalFlowProperties.h"
-
-namespace ProcessLib
-{
-namespace HeatTransportBHE
-{
-namespace BHE
-{
-BHE_CXA::BHE_CXA(BoreholeGeometry const& borehole,
-                 RefrigerantProperties const& refrigerant,
-                 GroutParameters const& grout,
-                 FlowAndTemperatureControl const& flowAndTemperatureControl,
-                 PipeConfigurationCoaxial const& pipes)
-    : BHECommonCoaxial{borehole, refrigerant, grout, flowAndTemperatureControl,
-                       pipes}
-{
-    // Initialize thermal resistances.
-    auto values = apply_visitor(
-        [&](auto const& control) {
-            return control(refrigerant.reference_temperature,
-                           0. /* initial time */);
-        },
-        flowAndTemperatureControl);
-    updateHeatTransferCoefficients(values.flow_rate);
-}
-
-/// Nu_o is the Nusselt number of inner pipe, Nu_i is the Nusselt number of
-/// annulus.
-std::array<double, BHE_CXA::number_of_unknowns> BHE_CXA::calcThermalResistances(
-    double const Nu_inner_pipe, double const Nu_annulus_pipe)
-{
-    // thermal resistances due to advective flow of refrigerant in the pipes
-    auto const R_advective = calculateAdvectiveThermalResistance(
-        _pipes.inner_pipe, _pipes.outer_pipe, refrigerant,
-        Nu_inner_pipe, Nu_annulus_pipe);
-
-    // thermal resistance due to thermal conductivity of the pipe wall material
-    auto const R_conductive = calculatePipeWallThermalResistance(
-        _pipes.inner_pipe, _pipes.outer_pipe);
-
-    // thermal resistance due to the grout transition and grout-soil exchange.
-    auto const R = calculateGroutAndGroutSoilExchangeThermalResistance(
-        _pipes.outer_pipe, grout, borehole_geometry.diameter);
-
-    // thermal resistance due to grout-soil exchange
-    double const R_gs = R.grout_soil;
-
-    // Eq. 56 and 57
-    double const R_ff = R_advective.inner_pipe_coaxial + R_advective.a_annulus +
-                        R_conductive.inner_pipe_coaxial;
-    double const R_fig =
-        R_advective.b_annulus + R_conductive.annulus + R.conductive_b;
-
-    return {{R_fig, R_ff, R_gs}};
-}
-}  // namespace BHE
-}  // namespace HeatTransportBHE
-}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
index 9ed03977f6e7262000ac7339b1d95f830698da78..741711baa36d570700a589454f54a3c7c994d163 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
@@ -14,8 +14,6 @@
 #include "BaseLib/Error.h"
 
 #include "BHECommonCoaxial.h"
-#include "FlowAndTemperatureControl.h"
-#include "PipeConfigurationCoaxial.h"
 
 namespace ProcessLib
 {
@@ -36,7 +34,6 @@ namespace BHE
  * surrounding soil is regulated through the thermal resistance values, which
  * are calculated specifically during the initialization of the class.
  */
-
 class BHE_CXA final : public BHECommonCoaxial
 {
 public:
@@ -44,7 +41,19 @@ public:
             RefrigerantProperties const& refrigerant,
             GroutParameters const& grout,
             FlowAndTemperatureControl const& flowAndTemperatureControl,
-            PipeConfigurationCoaxial const& pipes);
+            PipeConfigurationCoaxial const& pipes)
+        : BHECommonCoaxial{borehole, refrigerant, grout,
+                           flowAndTemperatureControl, pipes}
+    {
+        // Initialize thermal resistances.
+        auto values = apply_visitor(
+            [&](auto const& control) {
+                return control(refrigerant.reference_temperature,
+                               0. /* initial time */);
+            },
+            flowAndTemperatureControl);
+        updateHeatTransferCoefficients(values.flow_rate);
+    }
 
     template <int NPoints, typename SingleUnknownMatrixType,
               typename RMatrixType, typename RPiSMatrixType,
@@ -99,19 +108,25 @@ public:
     }
 
 public:
-    std::array<double, number_of_unknowns> const cross_section_areas = {
-        {_pipes.outer_pipe.area() - _pipes.inner_pipe.outsideArea(),
-         _pipes.inner_pipe.area(),
-         borehole_geometry.area() - _pipes.outer_pipe.outsideArea()}};
+    std::array<double, number_of_unknowns> cross_section_areas_coaxial()
+        const override
+    {
+        return {cross_section_area_annulus, cross_section_area_inner_pipe,
+                cross_section_area_grout};
+    }
 
 private:
-    std::array<double, number_of_unknowns> calcThermalResistances(
-        double const Nu_inner_pipe, double const Nu_annulus_pipe) override;
-
     std::array<double, 2> velocities() const override
     {
         return {_flow_velocity_annulus, _flow_velocity_inner};
     }
+
+    std::array<double, number_of_unknowns> getThermalResistances(
+        double const& R_gs, double const& R_ff,
+        double const& R_fg) const override
+    {
+        return {R_fg, R_ff, R_gs};
+    }
 };
 }  // namespace BHE
 }  // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.cpp
deleted file mode 100644
index 6f2aff4b8eecc7dd5167e48faf4015cc5e2cf727..0000000000000000000000000000000000000000
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.cpp
+++ /dev/null
@@ -1,71 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "BHE_CXC.h"
-#include "ThermoMechanicalFlowProperties.h"
-
-namespace ProcessLib
-{
-namespace HeatTransportBHE
-{
-namespace BHE
-{
-BHE_CXC::BHE_CXC(BoreholeGeometry const& borehole,
-                 RefrigerantProperties const& refrigerant,
-                 GroutParameters const& grout,
-                 FlowAndTemperatureControl const& flowAndTemperatureControl,
-                 PipeConfigurationCoaxial const& pipes)
-    : BHECommonCoaxial{borehole, refrigerant, grout, flowAndTemperatureControl,
-                       pipes}
-{
-    // Initialize thermal resistances.
-    auto values = apply_visitor(
-        [&](auto const& control) {
-            return control(refrigerant.reference_temperature,
-                           0. /* initial time */);
-        },
-        flowAndTemperatureControl);
-    updateHeatTransferCoefficients(values.flow_rate);
-}
-
-/// Nu_o is the Nusselt number of inner pipe, Nu_i is the Nusselt number of
-/// annulus.
-std::array<double, BHE_CXC::number_of_unknowns> BHE_CXC::calcThermalResistances(
-    double const Nu_inner_pipe, double const Nu_annulus_pipe)
-{
-    // thermal resistances due to advective flow of refrigerant in the pipes
-    auto const R_advective =
-        calculateAdvectiveThermalResistance(_pipes.inner_pipe,
-                                            _pipes.outer_pipe,
-                                            refrigerant,
-                                            Nu_inner_pipe,
-                                            Nu_annulus_pipe);
-
-    // thermal resistance due to thermal conductivity of the pipe wall material
-    auto const R_conductive = calculatePipeWallThermalResistance(
-        _pipes.inner_pipe, _pipes.outer_pipe);
-
-    // thermal resistance due to the grout transition grout-soil exchange.
-    auto const R = calculateGroutAndGroutSoilExchangeThermalResistance(
-        _pipes.outer_pipe, grout, borehole_geometry.diameter);
-
-    // thermal resistance due to grout-soil exchange
-    double const R_gs = R.grout_soil;
-
-    // Eq. 71 and 72
-    double const R_ff = R_advective.inner_pipe_coaxial + R_advective.a_annulus +
-                        R_conductive.inner_pipe_coaxial;
-    double const R_fog =
-        R_advective.b_annulus + R_conductive.annulus + R.conductive_b;
-
-    return {{R_ff, R_fog, R_gs}};
-}
-}  // namespace BHE
-}  // namespace HeatTransportBHE
-}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
index 7843b1223668e73139d33ac47d4d488ec8f048ff..0c7abd97e96d8f662ce542a5af2aa4deffc05234 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
@@ -14,8 +14,6 @@
 #include "BaseLib/Error.h"
 
 #include "BHECommonCoaxial.h"
-#include "FlowAndTemperatureControl.h"
-#include "PipeConfigurationCoaxial.h"
 
 namespace ProcessLib
 {
@@ -43,7 +41,19 @@ public:
             RefrigerantProperties const& refrigerant,
             GroutParameters const& grout,
             FlowAndTemperatureControl const& flowAndTemperatureControl,
-            PipeConfigurationCoaxial const& pipes);
+            PipeConfigurationCoaxial const& pipes)
+        : BHECommonCoaxial{borehole, refrigerant, grout,
+                           flowAndTemperatureControl, pipes}
+    {
+        // Initialize thermal resistances.
+        auto values = apply_visitor(
+            [&](auto const& control) {
+                return control(refrigerant.reference_temperature,
+                               0. /* initial time */);
+            },
+            flowAndTemperatureControl);
+        updateHeatTransferCoefficients(values.flow_rate);
+    }
 
     template <int NPoints, typename SingleUnknownMatrixType,
               typename RMatrixType, typename RPiSMatrixType,
@@ -98,19 +108,25 @@ public:
     }
 
 public:
-    std::array<double, number_of_unknowns> const cross_section_areas = {
-        {_pipes.inner_pipe.area(),
-         _pipes.outer_pipe.area() - _pipes.inner_pipe.outsideArea(),
-         borehole_geometry.area() - _pipes.outer_pipe.outsideArea()}};
+    std::array<double, number_of_unknowns> cross_section_areas_coaxial()
+        const override
+    {
+        return {cross_section_area_inner_pipe, cross_section_area_annulus,
+                cross_section_area_grout};
+    }
 
 private:
-    std::array<double, number_of_unknowns> calcThermalResistances(
-        double const Nu_inner_pipe, double const Nu_annulus_pipe) override;
-
     std::array<double, 2> velocities() const override
     {
         return {_flow_velocity_inner, _flow_velocity_annulus};
     }
+
+    std::array<double, number_of_unknowns> getThermalResistances(
+        double const& R_gs, double const& R_ff,
+        double const& R_fg) const override
+    {
+        return {R_ff, R_fg, R_gs};
+    }
 };
 }  // namespace BHE
 }  // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/Physics.h b/ProcessLib/HeatTransportBHE/BHE/Physics.h
index 8f7dc069545237d34fb235bf5e52c16f6ba5256e..89aa0dc9c8e42d2bd0a5ceba67b354c502bb5d30 100644
--- a/ProcessLib/HeatTransportBHE/BHE/Physics.h
+++ b/ProcessLib/HeatTransportBHE/BHE/Physics.h
@@ -11,6 +11,8 @@
 
 #pragma once
 
+#include "Pipe.h"
+
 namespace ProcessLib
 {
 namespace HeatTransportBHE
@@ -102,6 +104,22 @@ inline double PipeOutsideDiameter(double const pipe_diameter,
 {
     return pipe_diameter + 2. * pipe_wall_thickness;
 }
+
+struct CrossSectionAreasCoaxial
+{
+    double const inner;
+    double const annulus;
+    double const grout;
+};
+
+inline CrossSectionAreasCoaxial calculateCrossSectionAreasCoaxial(
+    Pipe const& inner_pipe, Pipe const& outer_pipe, double const borehole_area)
+{
+    double const inner = inner_pipe.area();
+    double const annulus = outer_pipe.area() - inner_pipe.outsideArea();
+    double const grout = borehole_area - outer_pipe.outsideArea();
+    return {inner, annulus, grout};
+}
 }  // end of namespace BHE
 }  // end of namespace HeatTransportBHE
 }  // end of namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
index e30497703e01df48a8226bf7f09494338aa76e25..26955a42d9e1839d282426e398f01273ccf46ed6 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -119,7 +119,7 @@ void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
     auto const& pipe_heat_capacities = _bhe.pipeHeatCapacities();
     auto const& pipe_heat_conductions = _bhe.pipeHeatConductions();
     auto const& pipe_advection_vectors = _bhe.pipeAdvectionVectors();
-    auto const& cross_section_areas = _bhe.cross_section_areas;
+    auto const& cross_section_areas = _bhe.CrossSectionAreas();
 
     // the mass and conductance matrix terms
     for (unsigned ip = 0; ip < n_integration_points; ip++)