From 413c652c7a422681cedb48dcd6b0448e744c2be5 Mon Sep 17 00:00:00 2001
From: Haibing Shao <haibing.shao@ufz.de>
Date: Mon, 21 Oct 2019 13:52:17 +0200
Subject: [PATCH] add a element direction vector in the BHE assembler class.

---
 .../HeatTransportBHE/BHE/BHECommonCoaxial.cpp |  3 +-
 .../HeatTransportBHE/BHE/BHECommonCoaxial.h   |  3 +-
 ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp    | 12 ++++---
 ProcessLib/HeatTransportBHE/BHE/BHE_1P.h      |  4 +--
 ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp    |  2 +-
 ProcessLib/HeatTransportBHE/BHE/BHE_1U.h      |  3 +-
 ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp    |  2 +-
 ProcessLib/HeatTransportBHE/BHE/BHE_2U.h      |  3 +-
 .../HeatTransportBHELocalAssemblerBHE-impl.h  | 31 ++++++++++---------
 .../HeatTransportBHELocalAssemblerBHE.h       |  5 +--
 10 files changed, 39 insertions(+), 29 deletions(-)

diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
index 5f3b2cc6627..f8a6d32eb6f 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
@@ -87,7 +87,8 @@ BHECommonCoaxial::pipeHeatConductions() const
 }
 
 std::array<Eigen::Vector3d, BHECommonCoaxial::number_of_unknowns>
-BHECommonCoaxial::pipeAdvectionVectors() const
+    BHECommonCoaxial::pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/) const
 {
     double const rho_r = refrigerant.density;
     double const Cp_r = refrigerant.specific_heat_capacity;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
index b033891a34e..03264584ad6 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
@@ -50,7 +50,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/)
         const;
 
     double cross_section_area_inner_pipe, cross_section_area_annulus,
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
index 9be0555074b..59a8eb1abcd 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
@@ -57,15 +57,17 @@ std::array<double, BHE_1P::number_of_unknowns> BHE_1P::pipeHeatConductions()
 }
 
 std::array<Eigen::Vector3d, BHE_1P::number_of_unknowns>
-BHE_1P::pipeAdvectionVectors() const
+BHE_1P::pipeAdvectionVectors(Eigen::Vector3d elem_direction_vec) const
 {
     double const& rho_r = refrigerant.density;
     double const& Cp_r = refrigerant.specific_heat_capacity;
+    Eigen::Vector3d adv_vector =
+        rho_r * Cp_r * _flow_velocity * elem_direction_vec;
 
-    return {{// pipe, Eq. 19
-             {0, 0, -rho_r * Cp_r * _flow_velocity},
-             // grout, Eq. 21
-             {0, 0, 0}}};
+    return {// pipe, Eq. 19
+            adv_vector,
+            // grout, Eq. 21
+            {0, 0, 0}};
 }
 
 double BHE_1P::compute_R_gs(double const chi, double const R_g)
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
index 82f9b9e2b41..d98916309b6 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
@@ -66,8 +66,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
-        const;
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d elem_direction_vec) const;
 
     template <int NPoints,
               typename SingleUnknownMatrixType,
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index c5451e452d2..915c9c20fb5 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -84,7 +84,7 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatConductions()
 }
 
 std::array<Eigen::Vector3d, BHE_1U::number_of_unknowns>
-BHE_1U::pipeAdvectionVectors() const
+    BHE_1U::pipeAdvectionVectors(Eigen::Vector3d /*elem_direction_vec*/) const
 {
     double const& rho_r = refrigerant.density;
     double const& Cp_r = refrigerant.specific_heat_capacity;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index 6702d7322a6..5b3b1d1b315 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -55,7 +55,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/)
         const;
 
     template <int NPoints,
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
index 734058abeab..ce09a6e0482 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -94,7 +94,7 @@ std::array<double, BHE_2U::number_of_unknowns> BHE_2U::pipeHeatConductions()
 }
 
 std::array<Eigen::Vector3d, BHE_2U::number_of_unknowns>
-BHE_2U::pipeAdvectionVectors() const
+    BHE_2U::pipeAdvectionVectors(Eigen::Vector3d /*elem_direction_vec*/) const
 {
     double const rho_r = refrigerant.density;
     double const Cp_r = refrigerant.specific_heat_capacity;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
index c63ff205d38..d3c69e92d8f 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
@@ -55,7 +55,8 @@ public:
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
-    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
+        Eigen::Vector3d /*elem_direction_vec*/)
         const;
 
     template <int NPoints,
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
index 33e35b1bd51..9583cd5d02b 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -57,6 +57,15 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
         _secondary_data.N[ip] = sm.N;
     }
 
+    // calculate the element direction vector
+    const double* coords_1 = e.getNode(0)->getCoords();
+    const double* coords_2 = e.getNode(1)->getCoords();
+    _element_direction_vector.setZero(3);
+    _element_direction_vector(0) = coords_2[0] - coords_1[0];
+    _element_direction_vector(1) = coords_2[1] - coords_1[1];
+    _element_direction_vector(2) = coords_2[2] - coords_1[2];
+    _element_direction_vector /= _element_direction_vector.norm();
+
     _R_matrix.setZero(bhe_unknowns_size, bhe_unknowns_size);
     _R_pi_s_matrix.setZero(bhe_unknowns_size, soil_temperature_size);
     _R_s_matrix.setZero(soil_temperature_size, soil_temperature_size);
@@ -65,12 +74,10 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
     for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < bhe_unknowns;
          idx_bhe_unknowns++)
     {
-        typename ShapeMatricesType::template MatrixType<
-            n_int_points, n_int_points>
+        typename ShapeMatricesType::template MatrixType<n_int_points,
+                                                        n_int_points>
             matBHE_loc_R = ShapeMatricesType::template MatrixType<
-                n_int_points,
-                n_int_points>::Zero(n_int_points,
-                                                n_int_points);
+                n_int_points, n_int_points>::Zero(n_int_points, n_int_points);
         // Loop over Gauss points
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
@@ -131,7 +138,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& pipe_advection_vectors = _bhe.pipeAdvectionVectors(_element_direction_vector);
     auto const& cross_section_areas = _bhe.crossSectionAreas();
 
     // the mass and conductance matrix terms
@@ -155,26 +162,22 @@ void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
             auto const& A = cross_section_areas[idx_bhe_unknowns];
 
             int const single_bhe_unknowns_index =
-                bhe_unknowns_index +
-                n_int_points * idx_bhe_unknowns;
+                bhe_unknowns_index + n_int_points * idx_bhe_unknowns;
             // local M
             local_M
-                .template block<n_int_points,
-                                n_int_points>(
+                .template block<n_int_points, n_int_points>(
                     single_bhe_unknowns_index, single_bhe_unknowns_index)
                 .noalias() += N.transpose() * N * mass_coeff * A * w;
 
             // local K
             // laplace part
             local_K
-                .template block<n_int_points,
-                                n_int_points>(
+                .template block<n_int_points, n_int_points>(
                     single_bhe_unknowns_index, single_bhe_unknowns_index)
                 .noalias() += dNdx.transpose() * dNdx * lambda * A * w;
             // advection part
             local_K
-                .template block<n_int_points,
-                                n_int_points>(
+                .template block<n_int_points, n_int_points>(
                     single_bhe_unknowns_index, single_bhe_unknowns_index)
                 .noalias() +=
                 N.transpose() * advection_vector.transpose() * dNdx * A * w;
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
index da02fb44654..d259d3884d2 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
@@ -32,8 +32,7 @@ class HeatTransportBHELocalAssemblerBHE
     static constexpr int n_int_points = ShapeFunction::NPOINTS;
     static constexpr int soil_temperature_size = ShapeFunction::NPOINTS;
     static constexpr int soil_temperature_index = 0;
-    static constexpr int bhe_unknowns_size =
-        n_int_points * bhe_unknowns;
+    static constexpr int bhe_unknowns_size = n_int_points * bhe_unknowns;
     static constexpr int bhe_unknowns_index = ShapeFunction::NPOINTS;
     static constexpr int local_matrix_size =
         soil_temperature_size + bhe_unknowns_size;
@@ -91,6 +90,8 @@ private:
 
     SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
 
+    Eigen::Vector3d _element_direction_vector;
+
     typename ShapeMatricesType::template MatrixType<bhe_unknowns_size,
                                                     bhe_unknowns_size>
         _R_matrix;
-- 
GitLab