diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 1f9f86beafbc864ef460bafbd686894a6f685453..d214974d3f142273449deaf60706d42d615a0fda 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -377,6 +377,13 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
#ifdef OGS_BUILD_PROCESS_HEATTRANSPORTBHE
if (type == "HEAT_TRANSPORT_BHE")
{
+ if (_mesh_vec[0]->getDimension() != 3)
+ {
+ OGS_FATAL(
+ "HEAT_TRANSPORT_BHE can only work with a 3-dimentional "
+ "mesh! ");
+ }
+
process =
ProcessLib::HeatTransportBHE::createHeatTransportBHEProcess(
*_mesh_vec[0], std::move(jacobian_assembler),
@@ -469,17 +476,17 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
{
case 2:
process =
- ProcessLib::PhaseField::createPhaseFieldProcess<
- 2>(*_mesh_vec[0], std::move(jacobian_assembler),
- _process_variables, _parameters,
- integration_order, process_config);
+ ProcessLib::PhaseField::createPhaseFieldProcess<2>(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
break;
case 3:
process =
- ProcessLib::PhaseField::createPhaseFieldProcess<
- 3>(*_mesh_vec[0], std::move(jacobian_assembler),
- _process_variables, _parameters,
- integration_order, process_config);
+ ProcessLib::PhaseField::createPhaseFieldProcess<3>(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
break;
}
}
diff --git a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
index b152909b53e648ed4ba08001d252996e21153086..f26e9b6623cb0bed1dc79e3e3c0a1ba109765d53 100644
--- a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
@@ -17,8 +17,7 @@
namespace ProcessLib
{
BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
- GlobalIndexType global_idx_T_in_top,
- GlobalIndexType global_idx_T_out_top,
+ std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
int const variable_id,
@@ -28,8 +27,7 @@ BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
: _bc_mesh(bc_mesh), _pt_bhe(pt_bhe)
{
DBUG(
- "Found %d nodes for BHE Inflow Dirichlet BCs for the variable %d "
- "and "
+ "Found %d nodes for BHE Inflow Dirichlet BCs for the variable %d and "
"component %d",
vec_inflow_bc_nodes.size(), variable_id, component_id);
@@ -48,8 +46,8 @@ BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
_bc_values.values.reserve(bc_mesh_subset.getNumberOfNodes());
// that might be slow, but only done once
- const auto g_idx_T_in = global_idx_T_in_top;
- const auto g_idx_T_out = global_idx_T_out_top;
+ const auto g_idx_T_in = in_out_global_indices.first;
+ const auto g_idx_T_out = in_out_global_indices.second;
if (g_idx_T_in >= 0 && g_idx_T_out >= 0)
{
@@ -99,7 +97,7 @@ void BHEInflowDirichletBoundaryCondition::preTimestep(const double /*t*/,
std::unique_ptr<BHEInflowDirichletBoundaryCondition>
createBHEInflowDirichletBoundaryCondition(
- GlobalIndexType global_idx_T_in_top, GlobalIndexType global_idx_T_out_top,
+ std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
int const variable_id, int const component_id,
@@ -108,10 +106,8 @@ createBHEInflowDirichletBoundaryCondition(
{
DBUG("Constructing BHEInflowDirichletBoundaryCondition.");
- //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Dirichlet__parameter}
-
return std::make_unique<BHEInflowDirichletBoundaryCondition>(
- global_idx_T_in_top, global_idx_T_out_top, bc_mesh, vec_inflow_bc_nodes,
+ std::move(in_out_global_indices), bc_mesh, vec_inflow_bc_nodes,
variable_id, component_id, pt_bhe);
}
-} // namespace ProcessLib
+} // namespace ProcessLib
\ No newline at end of file
diff --git a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
index 53c2afc0f4f784cfcee5c674e98f5128bf83ac29..442473b1492c8be48629c73cb78770ad4d7e7158 100644
--- a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
@@ -21,8 +21,7 @@ class BHEInflowDirichletBoundaryCondition final : public BoundaryCondition
{
public:
BHEInflowDirichletBoundaryCondition(
- GlobalIndexType global_idx_T_in_top,
- GlobalIndexType global_idx_T_out_top,
+ std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
int const variable_id,
@@ -37,6 +36,7 @@ public:
void preTimestep(const double t, const GlobalVector& x) override;
private:
+ // TODO (haibing) re-organize as the bottom BC data structure
MeshLib::Mesh const& _bc_mesh;
/// Stores the results of the outflow temperatures per boundary node.
@@ -51,7 +51,7 @@ private:
std::unique_ptr<BHEInflowDirichletBoundaryCondition>
createBHEInflowDirichletBoundaryCondition(
- GlobalIndexType global_idx_T_in_top, GlobalIndexType global_idx_T_out_top,
+ std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
int const variable_id, int const component_id,
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h b/ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h
index 9a5631db1a40c2f0cd0b3482e2b16c5f46082266..09126e877d2de5e95bad4c9db5a8e35655619aef 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h
@@ -35,6 +35,11 @@
#include "ProcessLib/Utils/ProcessUtils.h"
#include "boost/math/constants/constants.hpp"
+#include "BoreholeGeometry.h"
+#include "GroutParameters.h"
+#include "PipeParameters.h"
+#include "RefrigerantParameters.h"
+
namespace ProcessLib
{
namespace HeatTransportBHE
@@ -76,127 +81,6 @@ enum class BHE_DISCHARGE_TYPE
class BHEAbstract
{
public:
- struct BoreholeGeometry
- {
- /**
- * length/depth of the BHE
- * unit is m
- */
- double length;
-
- /**
- * diameter of the BHE
- * unit is m
- */
- double diameter;
- };
-
- struct PipeParameters
- {
- /**
- * radius of the pipline inner side
- * unit is m
- */
- double r_inner;
-
- /**
- * radius of the pipline outer side
- * unit is m
- */
- double r_outer;
-
- /**
- * pipe-in wall thickness
- * unit is m
- */
- double b_in;
-
- /**
- * pipe-out wall thickness
- * unit is m
- */
- double b_out;
-
- /**
- * thermal conductivity of the pipe wall
- * unit is kg m sec^-3 K^-1
- */
-
- double lambda_p;
-
- /**
- * thermal conductivity of the inner pipe wall
- * unit is kg m sec^-3 K^-1
- */
- double lambda_p_i;
-
- /**
- * thermal conductivity of the outer pipe wall
- * unit is kg m sec^-3 K^-1
- */
- double lambda_p_o;
- };
-
- struct RefrigerantParameters
- {
- /**
- * dynamics viscosity of the refrigerant
- * unit is kg m-1 sec-1
- */
- double mu_r;
-
- /**
- * density of the refrigerant
- * unit is kg m-3
- */
- double rho_r;
-
- /**
- * thermal conductivity of the refrigerant
- * unit is kg m sec^-3 K^-1
- */
- double lambda_r;
-
- /**
- * specific heat capacity of the refrigerant
- * unit is m^2 sec^-2 K^-1
- */
- double heat_cap_r;
-
- /**
- * longitudinal dispersivity of the
- * referigerant flow in the pipeline
- */
- double alpha_L;
- };
-
- struct GroutParameters
- {
- /**
- * density of the grout
- * unit is kg m-3
- */
- double rho_g;
-
- /**
- * porosity of the grout
- * unit is [-]
- */
- double porosity_g;
-
- /**
- * specific heat capacity of the grout
- * unit is m^2 sec^-2 K^-1
- */
- double heat_cap_g;
-
- /**
- * thermal conductivity of the grout
- * unit is kg m sec^-3 K^-1
- */
- double lambda_g;
- };
-
struct ExternallyDefinedRaRb
{
/**
@@ -300,15 +184,7 @@ public:
* initialization calcultion,
* need to be overwritten.
*/
- virtual void initialize()
- {
- calcPipeFlowVelocity();
- calcRenoldsNum();
- calcPrandtlNum();
- calcNusseltNum();
- calcThermalResistances();
- calcHeatTransferCoefficients();
- };
+ virtual void initialize() = 0;
/**
* update all parameters based on the new flow rate
@@ -317,12 +193,7 @@ public:
virtual void updateFlowRate(double new_flow_rate)
{
Q_r = new_flow_rate;
- calcPipeFlowVelocity();
- calcRenoldsNum();
- calcPrandtlNum();
- calcNusseltNum();
- calcThermalResistances();
- calcHeatTransferCoefficients();
+ initialize();
};
virtual void updateFlowRateFromCurve(double current_time) = 0;
@@ -334,28 +205,61 @@ public:
virtual void calcThermalResistances() = 0;
/**
- * Nusselt number calculation,
- * need to be overwritten.
- */
- virtual void calcNusseltNum() = 0;
-
- /**
- * Renolds number calculation,
- * need to be overwritten.
+ * flow velocity inside the pipeline
*/
- virtual void calcRenoldsNum() = 0;
+ double calcPipeFlowVelocity(double const& flow_rate,
+ double const& pipe_diameter)
+ {
+ return 4.0 * flow_rate / (PI * pipe_diameter * pipe_diameter);
+ }
- /**
- * Prandtl number calculation,
- * need to be overwritten.
- */
- virtual void calcPrandtlNum() = 0;
+ double calcPrandtlNumber(double const& viscosity,
+ double const& heat_capacity,
+ double const& heat_conductivity)
+ {
+ return viscosity * heat_capacity / heat_conductivity;
+ }
- /**
- * flow velocity inside the pipeline
- * need to be overwritten.
- */
- virtual void calcPipeFlowVelocity() = 0;
+ double calcRenoldsNumber(double const velocity_norm,
+ double const pipe_diameter,
+ double const viscosity,
+ double const density)
+ {
+ return velocity_norm * pipe_diameter / (viscosity / density);
+ };
+ double calcNusseltNumber(double const pipe_diameter,
+ double const pipe_length)
+ {
+ double tmp_Nu(0.0);
+ double gamma(0.0), xi(0.0);
+
+ if (Re < 2300.0)
+ {
+ tmp_Nu = 4.364;
+ }
+ else if (Re >= 2300.0 && Re < 10000.0)
+ {
+ gamma = (Re - 2300) / (10000 - 2300);
+
+ tmp_Nu = (1.0 - gamma) * 4.364;
+ tmp_Nu +=
+ gamma *
+ ((0.0308 / 8.0 * 1.0e4 * Pr) /
+ (1.0 + 12.7 * std::sqrt(0.0308 / 8.0) *
+ (std::pow(Pr, 2.0 / 3.0) - 1.0)) *
+ (1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0)));
+ }
+ else if (Re > 10000.0)
+ {
+ xi = pow(1.8 * std::log10(Re) - 1.5, -2.0);
+ tmp_Nu = (xi / 8.0 * Re * Re) /
+ (1.0 + 12.7 * std::sqrt(xi / 8.0) *
+ (std::pow(Pr, 2.0 / 3.0) - 1.0)) *
+ (1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0));
+ }
+
+ return tmp_Nu;
+ };
/**
* heat transfer coefficient,
@@ -409,6 +313,9 @@ public:
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
R_s_matrix) const = 0;
+ virtual std::vector<std::pair<int, int>> const&
+ inflowOutflowBcComponentIds() const = 0;
+
public:
/**
* name of the borehole heat exchanger
@@ -530,6 +437,21 @@ protected:
double Q_r;
const double PI;
+
+ /**
+ * Reynolds number
+ */
+ double Re;
+
+ /**
+ * Prandtl number
+ */
+ double Pr;
+
+ /**
+ * pipe distance
+ */
+ double omega;
};
} // end of namespace BHE
} // end of namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index debfa4ea3602ae8896a8d242375f24ba69444422..85bbd24c3d14fbb90cb1c50f1e6a2b72587716b0 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -11,6 +11,33 @@
using namespace ProcessLib::HeatTransportBHE::BHE;
+void BHE_1U::initialize()
+
+{
+ double tmp_u = calcPipeFlowVelocity(Q_r, 2.0 * pipe_param.r_inner);
+ _u(0) = tmp_u;
+ _u(1) = tmp_u;
+
+ // calculate Renolds number
+ Re = calcRenoldsNumber(std::abs(_u(0)),
+ 2.0 * pipe_param.r_inner,
+ refrigerant_param.mu_r,
+ refrigerant_param.rho_r);
+ // calculate Prandtl number
+ Pr = calcPrandtlNumber(refrigerant_param.mu_r,
+ refrigerant_param.heat_cap_r,
+ refrigerant_param.lambda_r);
+
+ // calculate Nusselt number
+ double tmp_Nu =
+ calcNusseltNumber(2.0 * pipe_param.r_inner, borehole_geometry.length);
+ _Nu(0) = tmp_Nu;
+ _Nu(1) = tmp_Nu;
+
+ calcThermalResistances();
+ calcHeatTransferCoefficients();
+}
+
/**
* calculate thermal resistance
*/
@@ -153,74 +180,6 @@ void BHE_1U::calcThermalResistances()
// -------------------------------------------------------------------------
}
-/**
- * Nusselt number calculation
- */
-void BHE_1U::calcNusseltNum()
-{
- // see Eq. 32 in Diersch_2011_CG
-
- double tmp_Nu = 0.0;
- double gamma, xi;
- double d;
- double const& L = borehole_geometry.length;
-
- d = 2.0 * pipe_param.r_inner;
-
- if (_Re < 2300.0)
- {
- tmp_Nu = 4.364;
- }
- else if (_Re >= 2300.0 && _Re < 10000.0)
- {
- gamma = (_Re - 2300) / (10000 - 2300);
-
- tmp_Nu = (1.0 - gamma) * 4.364;
- tmp_Nu += gamma * ((0.0308 / 8.0 * 1.0e4 * _Pr) /
- (1.0 + 12.7 * std::sqrt(0.0308 / 8.0) *
- (std::pow(_Pr, 2.0 / 3.0) - 1.0)) *
- (1.0 + std::pow(d / L, 2.0 / 3.0)));
- }
- else if (_Re > 10000.0)
- {
- xi = pow(1.8 * std::log10(_Re) - 1.5, -2.0);
- tmp_Nu = (xi / 8.0 * _Re * _Pr) /
- (1.0 + 12.7 * std::sqrt(xi / 8.0) *
- (std::pow(_Pr, 2.0 / 3.0) - 1.0)) *
- (1.0 + std::pow(d / L, 2.0 / 3.0));
- }
-
- _Nu(0) = tmp_Nu;
- _Nu(1) = tmp_Nu;
-}
-
-/**
- * Renolds number calculation
- */
-void BHE_1U::calcRenoldsNum()
-{
- double u_norm, d;
- double const& mu_r = refrigerant_param.mu_r;
- double const& rho_r = refrigerant_param.rho_r;
-
- u_norm = std::abs(_u(0));
- d = 2.0 * pipe_param.r_inner; // inner diameter of the pipeline
-
- _Re = u_norm * d / (mu_r / rho_r);
-}
-
-/**
- * Prandtl number calculation
- */
-void BHE_1U::calcPrandtlNum()
-{
- double const& mu_r = refrigerant_param.mu_r;
- double const& heat_cap_r = refrigerant_param.heat_cap_r;
- double const& lambda_r = refrigerant_param.lambda_r;
-
- _Pr = mu_r * heat_cap_r / lambda_r;
-}
-
/**
* calculate heat transfer coefficient
*/
@@ -232,19 +191,6 @@ void BHE_1U::calcHeatTransferCoefficients()
_PHI_gs = 1.0 / _R_gs;
}
-/**
- * flow velocity inside the pipeline
- */
-void BHE_1U::calcPipeFlowVelocity()
-{
- double tmp_u;
-
- tmp_u = Q_r / (PI * pipe_param.r_inner * pipe_param.r_inner);
-
- _u(0) = tmp_u;
- _u(1) = tmp_u;
-}
-
double BHE_1U::getMassCoeff(std::size_t idx_unknown) const
{
double const& rho_r = refrigerant_param.rho_r;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index a0d9689a2b22f5d3a9719ff84968508d1a81b079..ec9620a6b981ba17ca2f0d216a2c4ab321e898aa 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -76,6 +76,7 @@ public:
{
_u = Eigen::Vector2d::Zero();
_Nu = Eigen::Vector2d::Zero();
+ // 1U type of BHE has 2 pipes
Q_r = my_Qr;
@@ -162,6 +163,8 @@ public:
*/
std::size_t getNumUnknowns() const { return 4; }
+ void initialize();
+
void updateFlowRateFromCurve(double current_time)
{
if (use_flowrate_curve)
@@ -177,26 +180,6 @@ public:
*/
void calcThermalResistances();
- /**
- * Nusselt number calculation
- */
- void calcNusseltNum();
-
- /**
- * Renolds number calculation
- */
- void calcRenoldsNum();
-
- /**
- * Prandtl number calculation
- */
- void calcPrandtlNum();
-
- /**
- * flow velocity inside the pipeline
- */
- void calcPipeFlowVelocity();
-
/**
* calculate heat transfer coefficient
*/
@@ -247,6 +230,12 @@ public:
*/
double getTinByTout(double T_out, double current_time);
+ std::vector<std::pair<int, int>> const& inflowOutflowBcComponentIds()
+ const override
+ {
+ return _inflow_outflow_bc_component_ids;
+ }
+
/**
* required by eigen library,
* to make sure the dynamically allocated class has
@@ -255,6 +244,9 @@ public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
private:
+ std::vector<std::pair<int, int>> const _inflow_outflow_bc_component_ids = {
+ {0, 1}};
+
/**
* thermal resistances
*/
@@ -296,37 +288,21 @@ private:
double _PHI_fig, _PHI_fog, _PHI_gg, _PHI_gs;
/**
- * Reynolds number
+ * specific exchange surfaces S
*/
- double _Re;
-
+ double S_i, S_o, S_g1, S_gs;
/**
- * Prandtl number
+ * cross section area
*/
- double _Pr;
-
+ double CSA_i, CSA_o, CSA_g1, CSA_g2;
/**
* Nusselt number
*/
Eigen::Vector2d _Nu;
-
/**
* flow velocity inside the pipeline
*/
Eigen::Vector2d _u;
-
- /**
- * pipe distance
- */
- double omega;
- /**
- * specific exchange surfaces S
- */
- double S_i, S_o, S_g1, S_gs;
- /**
- * cross section area
- */
- double CSA_i, CSA_o, CSA_g1, CSA_g2;
};
} // namespace BHE
} // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
new file mode 100644
index 0000000000000000000000000000000000000000..658aa9816a3140857fddd45ddfba9b91a65272dd
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
@@ -0,0 +1,36 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct BoreholeGeometry
+{
+ /**
+ * length/depth of the BHE
+ * unit is m
+ */
+ double const length;
+
+ /**
+ * diameter of the BHE
+ * unit is m
+ */
+ double const diameter;
+};
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp b/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
index d01467d69b13d94c1363b5ea9b18a06568211750..1ec84499d054be780769287c38092a19a9889c60 100644
--- a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
@@ -246,4 +246,4 @@ BHE::BHE_1U* CreateBHE1U(
}
} // namespace BHE
} // namespace HeatTransportBHE
-} // namespace ProcessLib
+} // namespace ProcessLib
\ No newline at end of file
diff --git a/ProcessLib/HeatTransportBHE/BHE/GroutParameters.h b/ProcessLib/HeatTransportBHE/BHE/GroutParameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..c1355838f4065c9acf112bcb3fdba2dffb96641f
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/GroutParameters.h
@@ -0,0 +1,48 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct GroutParameters
+{
+ /**
+ * density of the grout
+ * unit is kg m-3
+ */
+ double const rho_g;
+
+ /**
+ * porosity of the grout
+ * unit is [-]
+ */
+ double const porosity_g;
+
+ /**
+ * specific heat capacity of the grout
+ * unit is m^2 sec^-2 K^-1
+ */
+ double const heat_cap_g;
+
+ /**
+ * thermal conductivity of the grout
+ * unit is kg m sec^-3 K^-1
+ */
+ double const lambda_g;
+};
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
index 1a55cc0c94d3f7f0588cff3c1c901db71b8e3b93..844ce563fa9e267f729e1034c3b57fe799cd8fef 100644
--- a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
@@ -18,42 +18,40 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
-void getBHEDataInMesh(
- MeshLib::Mesh const& mesh,
- std::vector<MeshLib::Element*>& vec_soil_elements,
- std::vector<int>& vec_BHE_mat_IDs,
- std::vector<std::vector<MeshLib::Element*>>& vec_BHE_elements,
- std::vector<MeshLib::Node*>& vec_soil_nodes,
- std::vector<std::vector<MeshLib::Node*>>& vec_BHE_nodes)
+BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh)
{
- // get vectors of matrix elements and BHE elements
- vec_soil_elements.reserve(mesh.getNumberOfElements());
+ BHEMeshData bheMeshData;
+ // partition all the mesh elements, and copy them into
+ // two seperate vectors, one with only matrix elements
+ // and the other only BHE elements
+ bheMeshData.soil_elements.reserve(mesh.getNumberOfElements());
std::vector<MeshLib::Element*> all_BHE_elements;
- for (MeshLib::Element* e : mesh.getElements())
- {
- // As for the first step, all elements are counted as a soil element
- // first. Those elements connected with a BHE will picked up and
- // reorganized into a seperate vector at the end of the function.
- if (e->getDimension() == mesh.getDimension())
- vec_soil_elements.push_back(e);
- else if (e->getDimension() == (unsigned int)1)
- all_BHE_elements.push_back(e);
- }
+ auto& all_mesh_elements = mesh.getElements();
+ std::partition_copy(
+ std::begin(all_mesh_elements),
+ std::end(all_mesh_elements),
+ std::back_inserter(all_BHE_elements),
+ std::back_inserter(bheMeshData.soil_elements),
+ [](MeshLib::Element* e) { return e->getDimension() == 1; });
// get BHE material IDs
- auto opt_material_ids(
- mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
+ auto opt_material_ids = MeshLib::materialIDs(mesh);
+ if (opt_material_ids == nullptr)
+ {
+ OGS_FATAL("Not able to get material IDs! ");
+ }
for (MeshLib::Element* e : all_BHE_elements)
- vec_BHE_mat_IDs.push_back((*opt_material_ids)[e->getID()]);
- BaseLib::makeVectorUnique(vec_BHE_mat_IDs);
- DBUG("-> found %d BHE material groups", vec_BHE_mat_IDs.size());
+ bheMeshData.BHE_mat_IDs.push_back((*opt_material_ids)[e->getID()]);
+ BaseLib::makeVectorUnique(bheMeshData.BHE_mat_IDs);
+ DBUG("-> found %d BHE material groups", bheMeshData.BHE_mat_IDs.size());
// create a vector of BHE elements for each group
- vec_BHE_elements.resize(vec_BHE_mat_IDs.size());
- for (unsigned bhe_id = 0; bhe_id < vec_BHE_mat_IDs.size(); bhe_id++)
+ bheMeshData.BHE_elements.resize(bheMeshData.BHE_mat_IDs.size());
+ for (unsigned bhe_id = 0; bhe_id < bheMeshData.BHE_mat_IDs.size(); bhe_id++)
{
- const auto bhe_mat_id = vec_BHE_mat_IDs[bhe_id];
- std::vector<MeshLib::Element*>& vec_elements = vec_BHE_elements[bhe_id];
+ const auto bhe_mat_id = bheMeshData.BHE_mat_IDs[bhe_id];
+ std::vector<MeshLib::Element*>& vec_elements =
+ bheMeshData.BHE_elements[bhe_id];
std::copy_if(all_BHE_elements.begin(), all_BHE_elements.end(),
std::back_inserter(vec_elements),
[&](MeshLib::Element* e) {
@@ -63,11 +61,11 @@ void getBHEDataInMesh(
}
// get a vector of BHE nodes
- vec_BHE_nodes.resize(vec_BHE_mat_IDs.size());
- for (unsigned bhe_id = 0; bhe_id < vec_BHE_mat_IDs.size(); bhe_id++)
+ bheMeshData.BHE_nodes.resize(bheMeshData.BHE_mat_IDs.size());
+ for (unsigned bhe_id = 0; bhe_id < bheMeshData.BHE_mat_IDs.size(); bhe_id++)
{
- std::vector<MeshLib::Node*>& vec_nodes = vec_BHE_nodes[bhe_id];
- for (MeshLib::Element* e : vec_BHE_elements[bhe_id])
+ std::vector<MeshLib::Node*>& vec_nodes = bheMeshData.BHE_nodes[bhe_id];
+ for (MeshLib::Element* e : bheMeshData.BHE_elements[bhe_id])
{
for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
{
@@ -82,20 +80,20 @@ void getBHEDataInMesh(
}
// get all the nodes into the pure soil nodes vector
- vec_soil_nodes.reserve(mesh.getNumberOfNodes());
+ bheMeshData.soil_nodes.reserve(mesh.getNumberOfNodes());
for (MeshLib::Node* n : mesh.getNodes())
{
// All elements are counted as a soil element
- vec_soil_nodes.push_back(n);
+ bheMeshData.soil_nodes.push_back(n);
}
- // final count of 3 types of elements
- // They are
- // (i) soil,
- // (ii) BHE
+ // finalLy counting two types of elements
+ // They are (i) soil, and (ii) BHE type of elements
DBUG("-> found total %d soil elements and %d BHE elements",
- vec_soil_elements.size(),
- all_BHE_elements.size());
+ bheMeshData.soil_elements.size(),
+ bheMeshData.BHE_elements.size());
+
+ return bheMeshData;
}
} // end of namespace HeatTransportBHE
-} // namespace ProcessLib
+} // namespace ProcessLib
\ No newline at end of file
diff --git a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
index ac75431d2e0057279b536fc2a13ba6b2f4d32154..d275b896ec2e0c41dbc02948a9544672f3f65def 100644
--- a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
+++ b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
@@ -21,22 +21,22 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
+struct BHEMeshData
+{
+ std::vector<MeshLib::Element*> soil_elements;
+ std::vector<int> BHE_mat_IDs;
+ std::vector<std::vector<MeshLib::Element*>> BHE_elements;
+ std::vector<MeshLib::Node*> soil_nodes;
+ std::vector<std::vector<MeshLib::Node*>> BHE_nodes;
+};
+
/**
* get data about fracture and matrix elements/nodes from a mesh
*
* @param mesh A mesh which includes BHE elements, i.e. 1-dimensional
* elements. It is assumed that elements forming a BHE have a distinct
* material ID.
- * @param vec_soil_elements a vector of soil elements
- * @param vec_BHE_elements a vector of BHE elements (grouped by BHE IDs)
- * @param vec_BHE_nodes a vector of BHE nodes (grouped by BHE IDs)
*/
-void getBHEDataInMesh(
- MeshLib::Mesh const& mesh,
- std::vector<MeshLib::Element*>& vec_soil_elements,
- std::vector<int>& vec_BHE_mat_IDs,
- std::vector<std::vector<MeshLib::Element*>>& vec_BHE_elements,
- std::vector<MeshLib::Node*>& vec_pure_soil_nodes,
- std::vector<std::vector<MeshLib::Node*>>& vec_BHE_nodes);
+BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh);
} // end of namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/PipeParameters.h b/ProcessLib/HeatTransportBHE/BHE/PipeParameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..1c016d7491186b836f21448f29fd25a1eaeb08cf
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/PipeParameters.h
@@ -0,0 +1,67 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct PipeParameters
+{
+ /**
+ * radius of the pipline inner side
+ * unit is m
+ */
+ double const r_inner;
+
+ /**
+ * radius of the pipline outer side
+ * unit is m
+ */
+ double const r_outer;
+
+ /**
+ * pipe-in wall thickness
+ * unit is m
+ */
+ double const b_in;
+
+ /**
+ * pipe-out wall thickness
+ * unit is m
+ */
+ double const b_out;
+
+ /**
+ * thermal conductivity of the pipe wall
+ * unit is kg m sec^-3 K^-1
+ */
+
+ double const lambda_p;
+
+ /**
+ * thermal conductivity of the inner pipe wall
+ * unit is kg m sec^-3 K^-1
+ */
+ double const lambda_p_i;
+
+ /**
+ * thermal conductivity of the outer pipe wall
+ * unit is kg m sec^-3 K^-1
+ */
+ double const lambda_p_o;
+};
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/RefrigerantParameters.h b/ProcessLib/HeatTransportBHE/BHE/RefrigerantParameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..8b46b31a3510667dc580dc1215fe298e50c552f9
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/RefrigerantParameters.h
@@ -0,0 +1,54 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct RefrigerantParameters
+{
+ /**
+ * dynamics viscosity of the refrigerant
+ * unit is kg m-1 sec-1
+ */
+ double const mu_r;
+
+ /**
+ * density of the refrigerant
+ * unit is kg m-3
+ */
+ double const rho_r;
+
+ /**
+ * thermal conductivity of the refrigerant
+ * unit is kg m sec^-3 K^-1
+ */
+ double const lambda_r;
+
+ /**
+ * specific heat capacity of the refrigerant
+ * unit is m^2 sec^-2 K^-1
+ */
+ double const heat_cap_r;
+
+ /**
+ * longitudinal dispersivity of the
+ * referigerant flow in the pipeline
+ */
+ double const alpha_L;
+};
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/CMakeLists.txt b/ProcessLib/HeatTransportBHE/CMakeLists.txt
index 2d9def64402e62d9154925303a0ade8532756dab..2838acf7603c13002f509cb6abfc37c9fe4366a9 100644
--- a/ProcessLib/HeatTransportBHE/CMakeLists.txt
+++ b/ProcessLib/HeatTransportBHE/CMakeLists.txt
@@ -5,4 +5,4 @@ APPEND_SOURCE_FILES(SOURCES LocalAssemblers)
add_library(HeatTransportBHE ${SOURCES})
target_link_libraries(HeatTransportBHE PUBLIC ProcessLib)
-# include(Tests.cmake)
+include(Tests.cmake)
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
index a8f1a0c8a3477dad8ee81aa169949e6534a57b08..5a9da59c74d1d1e0e95645bcd7817d436fd66513 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
@@ -37,39 +37,36 @@ HeatTransportBHEProcess::HeatTransportBHEProcess(
: Process(mesh, std::move(jacobian_assembler), parameters,
integration_order, std::move(process_variables),
std::move(secondary_variables), std::move(named_function_caller)),
- _process_data(std::move(process_data))
+ _process_data(std::move(process_data)),
+ _bheMeshData(getBHEDataInMesh(mesh))
{
- getBHEDataInMesh(mesh,
- _vec_pure_soil_elements,
- _vec_BHE_mat_IDs,
- _vec_BHE_elements,
- _vec_pure_soil_nodes,
- _vec_BHE_nodes);
-
- if (_vec_BHE_mat_IDs.size() != _process_data._vec_BHE_property.size())
+ if (_bheMeshData.BHE_mat_IDs.size() !=
+ _process_data._vec_BHE_property.size())
{
OGS_FATAL(
"The number of the given BHE properties (%d) are not "
"consistent"
" with the number of BHE groups in a mesh (%d).",
_process_data._vec_BHE_property.size(),
- _vec_BHE_mat_IDs.size());
+ _bheMeshData.BHE_mat_IDs.size());
}
+ auto material_ids = MeshLib::materialIDs(mesh);
+ if (material_ids == nullptr)
+ {
+ OGS_FATAL("Not able to get material IDs! ");
+ }
// create a map from a material ID to a BHE ID
auto max_BHE_mat_id =
- std::max_element(_vec_BHE_mat_IDs.begin(), _vec_BHE_mat_IDs.end());
- const std::size_t nBHEmatIDs = *max_BHE_mat_id + 1;
- _process_data._map_materialID_to_BHE_ID.resize(nBHEmatIDs);
- for (std::size_t i = 0; i < _vec_BHE_mat_IDs.size(); i++)
+ std::max_element(material_ids->begin(), material_ids->end());
+ _process_data._map_materialID_to_BHE_ID.resize(*max_BHE_mat_id + 1);
+ for (std::size_t i = 0; i < _bheMeshData.BHE_mat_IDs.size(); i++)
{
- // by default, it is assumed that the soil compartment takes material ID
- // 0 and the BHE take the successive material group.
- _process_data._map_materialID_to_BHE_ID[_vec_BHE_mat_IDs[i]] = i;
+ // fill in the map structure
+ _process_data._map_materialID_to_BHE_ID[_bheMeshData.BHE_mat_IDs[i]] =
+ i;
}
- MeshLib::PropertyVector<int> const* material_ids(
- mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
_process_data._mesh_prop_materialIDs = material_ids;
}
@@ -84,16 +81,16 @@ void HeatTransportBHEProcess::constructDofTable()
//------------------------------------------------------------
// first all the soil nodes
_mesh_subset_pure_soil_nodes =
- std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_pure_soil_nodes);
+ std::make_unique<MeshLib::MeshSubset>(_mesh, _bheMeshData.soil_nodes);
std::vector<int> vec_n_BHE_unknowns;
- vec_n_BHE_unknowns.reserve(_vec_BHE_nodes.size());
+ vec_n_BHE_unknowns.reserve(_bheMeshData.BHE_nodes.size());
// the BHE nodes need to be cherry-picked from the vector
- for (unsigned i = 0; i < _vec_BHE_nodes.size(); i++)
+ for (unsigned i = 0; i < _bheMeshData.BHE_nodes.size(); i++)
{
_mesh_subset_BHE_nodes.push_back(
- std::make_unique<MeshLib::MeshSubset const>(_mesh,
- _vec_BHE_nodes[i]));
+ std::make_unique<MeshLib::MeshSubset const>(
+ _mesh, _bheMeshData.BHE_nodes[i]));
int n_BHE_unknowns =
_process_data._vec_BHE_property[i]->getNumUnknowns();
vec_n_BHE_unknowns.emplace_back(n_BHE_unknowns);
@@ -122,7 +119,7 @@ void HeatTransportBHEProcess::constructDofTable()
// temperature.
vec_n_components.push_back(1);
// now the BHE subsets
- for (std::size_t i = 0; i < _vec_BHE_mat_IDs.size(); i++)
+ for (std::size_t i = 0; i < _bheMeshData.BHE_mat_IDs.size(); i++)
{
// Here the number of components equals to
// the number of unknowns on the BHE
@@ -132,9 +129,9 @@ void HeatTransportBHEProcess::constructDofTable()
std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
// vec_var_elements.push_back(&_vec_pure_soil_elements);
vec_var_elements.push_back(&(_mesh.getElements()));
- for (std::size_t i = 0; i < _vec_BHE_elements.size(); i++)
+ for (std::size_t i = 0; i < _bheMeshData.BHE_elements.size(); i++)
{
- vec_var_elements.push_back(&_vec_BHE_elements[i]);
+ vec_var_elements.push_back(&_bheMeshData.BHE_elements[i]);
}
_local_to_global_index_map =
@@ -153,12 +150,10 @@ void HeatTransportBHEProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
- const int process_id = 0;
-
- // this process can only run with 3-dimensional mesh
+ assert(mesh.getDimension() == 3);
ProcessLib::HeatTransportBHE::createLocalAssemblers<
- 3, /*mesh.getDimension(),*/
- HeatTransportBHELocalAssemblerSoil, HeatTransportBHELocalAssemblerBHE>(
+ 3 /* mesh dimension */, HeatTransportBHELocalAssemblerSoil,
+ HeatTransportBHELocalAssemblerBHE>(
mesh.getElements(), dof_table, _local_assemblers,
_process_data._vec_ele_connected_BHE_IDs,
_process_data._vec_BHE_property, mesh.isAxiallySymmetric(),
@@ -169,10 +164,13 @@ void HeatTransportBHEProcess::initializeConcreteProcess(
// and one BC at the bottom.
std::vector<std::unique_ptr<BoundaryCondition>> bc_collections =
createBHEBoundaryConditionTopBottom();
+
+ const int process_id = 0;
auto& current_process_BCs = _boundary_conditions[process_id];
- auto const bc_collections_size = bc_collections.size();
- for (unsigned i = 0; i < bc_collections_size; i++)
- current_process_BCs.addCreatedBC(std::move(bc_collections[i]));
+ for (auto& bc_collection : bc_collections)
+ {
+ current_process_BCs.addCreatedBC(std::move(bc_collection));
+ }
}
void HeatTransportBHEProcess::assembleConcreteProcess(const double t,
@@ -192,19 +190,12 @@ void HeatTransportBHEProcess::assembleConcreteProcess(const double t,
}
void HeatTransportBHEProcess::assembleWithJacobianConcreteProcess(
- const double t, GlobalVector const& x, GlobalVector const& xdot,
- const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b, GlobalMatrix& Jac)
+ const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/,
+ const double /*dxdot_dx*/, const double /*dx_dx*/, GlobalMatrix& /*M*/,
+ GlobalMatrix& /*K*/, GlobalVector& /*b*/, GlobalMatrix& /*Jac*/)
{
- DBUG("AssembleWithJacobian HeatTransportBHE process.");
-
- std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
- dof_table = {std::ref(*_local_to_global_index_map)};
- // Call global assembler for each local assembly item.
- GlobalExecutor::executeMemberDereferenced(
- _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
- _local_assemblers, dof_table, t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac,
- _coupled_solutions);
+ OGS_FATAL(
+ "HeatTransportBHE: analytical Jacobian assembly is not implemented");
}
void HeatTransportBHEProcess::computeSecondaryVariableConcrete(
@@ -227,189 +218,51 @@ HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom()
int const n_BHEs = static_cast<int>(_process_data._vec_BHE_property.size());
- // bulk dof table
- NumLib::LocalToGlobalIndexMap dof_table_bulk = *_local_to_global_index_map;
-
- std::vector<MeshLib::Node*> bc_top_nodes;
- std::vector<MeshLib::Node*> bc_bottom_nodes;
-
// for each BHE
for (int bhe_i = 0; bhe_i < n_BHEs; bhe_i++)
{
- // get the BHE type
- auto const& bhe_typ = _process_data._vec_BHE_property[bhe_i]->bhe_type;
+ auto const& bhe_nodes = _bheMeshData.BHE_nodes[bhe_i];
// find the variable ID
// the soil temperature is 0-th variable
// the BHE temperature is therefore bhe_i + 1
const int variable_id = bhe_i + 1;
- // find the node in mesh that are at the top
- auto const n_bhe_nodes = _vec_BHE_nodes[bhe_i].size();
- unsigned int idx_top = 0;
- unsigned int idx_bottom = _vec_BHE_nodes[bhe_i].size() - 1;
- double top_z_coord = _vec_BHE_nodes[bhe_i][idx_top]->getCoords()[2];
- double bottom_z_coord =
- _vec_BHE_nodes[bhe_i][idx_bottom]->getCoords()[2];
- for (unsigned bhe_node_i = 0; bhe_node_i < n_bhe_nodes; bhe_node_i++)
- {
- if (_vec_BHE_nodes[bhe_i][bhe_node_i]->getCoords()[2] >=
- top_z_coord)
- idx_top = bhe_node_i;
-
- if (_vec_BHE_nodes[bhe_i][bhe_node_i]->getCoords()[2] <=
- bottom_z_coord)
- idx_bottom = bhe_node_i;
- }
- bc_top_nodes.clear();
- bc_top_nodes.emplace_back(_vec_BHE_nodes[bhe_i][idx_top]);
- bc_bottom_nodes.clear();
- bc_bottom_nodes.emplace_back(_vec_BHE_nodes[bhe_i][idx_bottom]);
-
- MeshLib::MeshSubset bc_mesh_subset_top{_mesh, bc_top_nodes};
- MeshLib::Mesh const& bc_mesh_top = bc_mesh_subset_top.getMesh();
-
- MeshLib::MeshSubset bc_mesh_subset_bottom{_mesh, bc_bottom_nodes};
- MeshLib::Mesh const& bc_mesh_bottom = bc_mesh_subset_bottom.getMesh();
- auto get_global_bhe_bc_index_top = [&](std::size_t const mesh_id,
- int const component_id) {
- return dof_table_bulk.getGlobalIndex(
- {mesh_id, MeshLib::MeshItemType::Node,
- bc_top_nodes[0]->getID()},
- variable_id, component_id);
- };
- auto get_global_bhe_bc_index_bottom = [&](std::size_t const mesh_id,
- int const component_id) {
- return dof_table_bulk.getGlobalIndex(
- {mesh_id, MeshLib::MeshItemType::Node,
- bc_bottom_nodes[0]->getID()},
- variable_id, component_id);
- };
-
- // the create_BC function will be repeatedly used
- auto create_BC_top_inflow = [&](std::size_t const mesh_id,
- int const component_id_in,
- int const component_id_out) {
- auto const global_index_in =
- get_global_bhe_bc_index_top(mesh_id, component_id_in);
- auto const global_index_out =
- get_global_bhe_bc_index_top(mesh_id, component_id_out);
- return ProcessLib::createBHEInflowDirichletBoundaryCondition(
- global_index_in, global_index_out, _mesh, bc_top_nodes,
- variable_id, component_id_in,
- _process_data._vec_BHE_property[bhe_i]);
- };
- auto create_BC_bottom_outflow = [&](std::size_t const mesh_id,
- int const component_id_in,
- int const component_id_out) {
- auto const global_index_in =
- get_global_bhe_bc_index_bottom(mesh_id, component_id_in);
- auto const global_index_out =
- get_global_bhe_bc_index_bottom(mesh_id, component_id_out);
- return ProcessLib::createBHEBottomDirichletBoundaryCondition(
- global_index_in, global_index_out, _mesh, bc_top_nodes,
- variable_id, component_id_in);
- };
-
- // depending on the BHE type
- switch (bhe_typ)
+ // Bottom and top nodes w.r.t. the z coordinate.
+ auto const bottom_top_nodes = std::minmax_element(
+ begin(bhe_nodes), end(bhe_nodes),
+ [&](auto const& a, auto const& b) {
+ return a->getCoords()[2] < b->getCoords()[2];
+ });
+ auto const bc_bottom_node_id = (*bottom_top_nodes.first)->getID();
+ MeshLib::Node* const bc_top_node = *bottom_top_nodes.second;
+
+ auto get_global_bhe_bc_indices =
+ [&](std::size_t const node_id,
+ std::pair<int, int> const& in_out_component_id) {
+ return std::make_pair(
+ _local_to_global_index_map->getGlobalIndex(
+ {_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
+ variable_id, in_out_component_id.first),
+ _local_to_global_index_map->getGlobalIndex(
+ {_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
+ variable_id, in_out_component_id.second));
+ };
+
+ for (auto const& in_out_component_id :
+ _process_data._vec_BHE_property[bhe_i]
+ ->inflowOutflowBcComponentIds())
{
- case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_2U:
- {
- unsigned const component_id_T_in_1 = 0;
- unsigned const component_id_T_out_1 = 2;
-
- unsigned const component_id_T_out_2 = 3;
- unsigned const component_id_T_in_2 = 1;
-
- // there are 2 BCs on the top node
- std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top_1 =
- create_BC_top_inflow(bc_mesh_top.getID(),
- component_id_T_in_1,
- component_id_T_out_1);
- std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top_2 =
- create_BC_top_inflow(bc_mesh_top.getID(),
- component_id_T_in_2,
- component_id_T_out_2);
-
- // there are also 2 BCs on the bottom node
- std::unique_ptr<BHEBottomDirichletBoundaryCondition>
- bc_bottom_1 = create_BC_bottom_outflow(
- bc_mesh_bottom.getID(), component_id_T_in_1,
- component_id_T_out_1);
- std::unique_ptr<BHEBottomDirichletBoundaryCondition>
- bc_bottom_2 = create_BC_bottom_outflow(
- bc_mesh_bottom.getID(), component_id_T_in_2,
- component_id_T_out_2);
-
- // add bc_top and bc_bottom to the vector
- bcs.push_back(std::move(bc_top_1));
- bcs.push_back(std::move(bc_top_2));
- bcs.push_back(std::move(bc_bottom_1));
- bcs.push_back(std::move(bc_bottom_2));
- }
- break;
- case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_1U:
- {
- unsigned const component_id_T_in = 0;
- unsigned const component_id_T_out = 1;
- // there is one BC on the top node
- std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top =
- create_BC_top_inflow(bc_mesh_top.getID(), component_id_T_in,
- component_id_T_out);
-
- // there is also 1 BC on the bottom node
- std::unique_ptr<BHEBottomDirichletBoundaryCondition> bc_bottom =
- create_BC_bottom_outflow(bc_mesh_bottom.getID(),
- component_id_T_in,
- component_id_T_out);
-
- // add bc_top and bc_bottom to the vector
- bcs.push_back(std::move(bc_top));
- bcs.push_back(std::move(bc_bottom));
- }
- break;
- case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_CXC:
- {
- unsigned const component_id_T_in = 1;
- unsigned const component_id_T_out = 0;
- // there is one BC on the top node
- std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top =
- create_BC_top_inflow(bc_mesh_top.getID(), component_id_T_in,
- component_id_T_out);
-
- // there is also 1 BC on the bottom node
- std::unique_ptr<BHEBottomDirichletBoundaryCondition> bc_bottom =
- create_BC_bottom_outflow(bc_mesh_bottom.getID(),
- component_id_T_in,
- component_id_T_out);
-
- // add bc_top and bc_bottom to the vector
- bcs.push_back(std::move(bc_top));
- bcs.push_back(std::move(bc_bottom));
- }
- break;
- case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_CXA:
- {
- unsigned const component_id_T_in = 0;
- unsigned const component_id_T_out = 1;
- // there is one BC on the top node
- std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top =
- create_BC_top_inflow(bc_mesh_top.getID(), component_id_T_in,
- component_id_T_out);
-
- // there is also 1 BC on the bottom node
- std::unique_ptr<BHEBottomDirichletBoundaryCondition> bc_bottom =
- create_BC_bottom_outflow(bc_mesh_bottom.getID(),
- component_id_T_in,
- component_id_T_out);
-
- // add bc_top and bc_bottom to the vector
- bcs.push_back(std::move(bc_top));
- bcs.push_back(std::move(bc_bottom));
- }
- break;
- default:
- OGS_FATAL("WRONG BHE TYPE");
+ // Top, inflow.
+ bcs.push_back(ProcessLib::createBHEInflowDirichletBoundaryCondition(
+ get_global_bhe_bc_indices(bc_top_node->getID(),
+ in_out_component_id),
+ _mesh, {bc_top_node}, variable_id, in_out_component_id.first,
+ _process_data._vec_BHE_property[bhe_i]));
+
+ // Bottom, outflow.
+ bcs.push_back(ProcessLib::createBHEBottomDirichletBoundaryCondition(
+ get_global_bhe_bc_indices(bc_bottom_node_id,
+ in_out_component_id)));
}
}
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
index 217e4da211538bd1e1b2a59f75c329ed2a4d215e..6af02ecd02a1a2afc4e25d43fe8d87b3669d8a54 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
@@ -10,6 +10,7 @@
#pragma once
#include "HeatTransportBHEProcessData.h"
+#include "ProcessLib/HeatTransportBHE/BHE/MeshUtils.h"
#include "ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h"
#include "ProcessLib/Process.h"
@@ -17,6 +18,8 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
+struct BHEMeshData;
+
class HeatTransportBHEProcess final : public Process
{
public:
@@ -64,35 +67,18 @@ private:
std::vector<std::unique_ptr<HeatTransportBHELocalAssemblerInterface>>
_local_assemblers;
- /**
- * These are the elements that are representing BHEs
- */
- std::vector<std::vector<MeshLib::Element*>> _vec_BHE_elements;
-
- /**
- * These are the elements that are not connected with any BHE
- */
- std::vector<MeshLib::Element*> _vec_pure_soil_elements;
-
- /**
- * These are the soil nodes that are not connected with a BHE
- */
- std::vector<MeshLib::Node*> _vec_pure_soil_nodes;
-
- /**
- * Mesh nodes that are located on any BHE
- * ordered according to each BHE
- */
- std::vector<std::vector<MeshLib::Node*>> _vec_BHE_nodes;
-
std::vector<std::unique_ptr<MeshLib::MeshSubset const>>
_mesh_subset_BHE_nodes;
+
std::vector<std::unique_ptr<MeshLib::MeshSubset const>>
_mesh_subset_BHE_soil_nodes;
+
std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_pure_soil_nodes;
+
std::unique_ptr<MeshLib::MeshSubset const>
_mesh_subset_soil_nodes_connected_with_BHE;
- std::vector<int> _vec_BHE_mat_IDs;
+
+ const BHEMeshData _bheMeshData;
};
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
index f7ea07d91062fd54fbce295eec09fb9b6a299b76..0a904d254e297f626b98fa3a82232b36c72f68ba 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
@@ -42,7 +42,6 @@ public:
HeatTransportBHELocalAssemblerBHE(
MeshLib::Element const& e,
- std::size_t const local_matrix_size,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool const is_axially_symmetric,
unsigned const integration_order,
@@ -53,16 +52,6 @@ public:
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_b_data*/) override;
- void assembleWithJacobian(double const /*t*/,
- Eigen::VectorXd const& /*local_u*/,
- Eigen::VectorXd& /*local_b*/,
- Eigen::MatrixXd& /*local_J*/) override
- {
- OGS_FATAL(
- "HeatTransportBHELocalAssembler: assembly with jacobian is not "
- "implemented.");
- }
-
void preTimestepConcrete(std::vector<double> const& /*local_x*/,
double const /*t*/,
double const /*delta_t*/) override
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
index 5285b4b2785a62419eda825da18d9d556ea7fd19..5d5d8c79caa0d95a595570bee259602387f92c40 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
@@ -25,7 +25,6 @@ template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHE_Dim>::
HeatTransportBHELocalAssemblerBHE(
MeshLib::Element const& e,
- std::size_t const /*local_matrix_size*/,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool const is_axially_symmetric,
unsigned const integration_order,
@@ -62,9 +61,8 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHE_Dim>::
{
x_position.setIntegrationPoint(ip);
- IntegrationPointDataBHE<ShapeMatricesType> int_Point_Data_BHE(
- *(_process_data._vec_BHE_property[BHE_id]));
- _ip_data.emplace_back(int_Point_Data_BHE);
+ // create the class IntegrationPointDataBHE in place
+ _ip_data.emplace_back(*(_process_data._vec_BHE_property[BHE_id]));
auto const& sm = shape_matrices[ip];
auto& ip_data = _ip_data[ip];
ip_data.integration_weight =
@@ -109,7 +107,6 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHE_Dim>::
_process_data._vec_BHE_property[BHE_id]->setRMatrices(
idx_bhe_unknowns, ShapeFunction::NPOINTS, matBHE_loc_R, _R_matrix,
_R_pi_s_matrix, _R_s_matrix);
-
} // end of loop over BHE unknowns
// debugging
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
index 3ee83392e2c0674c856465e1fdbf1e399875bd22..f556f807109d5349a5f411ba72a3c22c22d7555b 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
@@ -26,8 +26,6 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
-const unsigned NUM_NODAL_DOF_SOIL = 1;
-
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
class HeatTransportBHELocalAssemblerSoil
: public HeatTransportBHELocalAssemblerInterface
@@ -45,7 +43,6 @@ public:
HeatTransportBHELocalAssemblerSoil(
MeshLib::Element const& e,
- std::size_t const local_matrix_size,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool is_axially_symmetric,
unsigned const integration_order,
@@ -56,29 +53,6 @@ public:
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_b_data*/) override;
- void assembleWithJacobian(double const /*t*/,
- std::vector<double> const& /*local_x*/,
- std::vector<double> const& /*local_xdot*/,
- const double /*dxdot_dx*/, const double /*dx_dx*/,
- std::vector<double>& /*local_M_data*/,
- std::vector<double>& /*local_K_data*/,
- std::vector<double>& /*local_b_data*/,
- std::vector<double>& /*local_Jac_data*/) override
- {
- OGS_FATAL(
- "HeatTransportBHELocalAssemblerMatrix: assembly with jacobian is "
- "not "
- "implemented.");
- }
-
- void preTimestepConcrete(std::vector<double> const& /*local_x*/,
- double const /*t*/,
- double const /*delta_t*/) override
- {
- }
-
- void postTimestepConcrete(std::vector<double> const& /*local_x*/) override;
-
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -92,8 +66,8 @@ private:
HeatTransportBHEProcessData& _process_data;
std::vector<
- IntegrationPointDataSoil<ShapeMatricesType>,
- Eigen::aligned_allocator<IntegrationPointDataSoil<ShapeMatricesType>>>
+ IntegrationPointDataSoil<NodalMatrixType>,
+ Eigen::aligned_allocator<IntegrationPointDataSoil<NodalMatrixType>>>
_ip_data;
IntegrationMethod const _integration_method;
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
index 3cf4209eefe010d583da12f566ae84e8f0c5f8a5..ffe4d351484de80c34f86c288afd8e2e7bd9fc6b 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
@@ -37,7 +37,6 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction,
GlobalDim>::
HeatTransportBHELocalAssemblerSoil(
MeshLib::Element const& e,
- std::size_t const /*local_matrix_size*/,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool const is_axially_symmetric,
unsigned const integration_order,
@@ -60,6 +59,26 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction,
IntegrationMethod,
GlobalDim>(
e, is_axially_symmetric, _integration_method);
+
+ SpatialPosition x_position;
+ x_position.setElementID(element_id);
+
+ // ip data initialization
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+
+ // create the class IntegrationPointDataBHE in place
+ _ip_data.emplace_back();
+ auto const& sm = _shape_matrices[ip];
+ auto& ip_data = _ip_data[ip];
+ double const w = _integration_method.getWeightedPoint(ip).getWeight() *
+ sm.integralMeasure * sm.detJ;
+ ip_data.NTN_product_times_w = sm.N.transpose() * sm.N * w;
+ ip_data.dNdxTdNdx_product_times_w = sm.dNdx.transpose() * sm.dNdx * w;
+
+ _secondary_data.N[ip] = sm.N;
+ }
}
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
@@ -72,14 +91,13 @@ void HeatTransportBHELocalAssemblerSoil<
std::vector<double>& local_K_data,
std::vector<double>& /*local_b_data*/)
{
- auto const local_matrix_size = local_x.size();
-
- assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF_SOIL);
+ assert(local_x.size() == ShapeFunction::NPOINTS);
+ (void)local_x; // Avoid unused arg warning.
auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
- local_M_data, local_matrix_size, local_matrix_size);
+ local_M_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
- local_K_data, local_matrix_size, local_matrix_size);
+ local_K_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -90,8 +108,9 @@ void HeatTransportBHELocalAssemblerSoil<
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
- auto const& sm = _shape_matrices[ip];
- auto const& wp = _integration_method.getWeightedPoint(ip);
+ auto& ip_data = _ip_data[ip];
+ auto const& M_ip = ip_data.NTN_product_times_w;
+ auto const& K_ip = ip_data.dNdxTdNdx_product_times_w;
// auto const k_f = _process_data.thermal_conductivity_fluid(t, pos)[0];
// auto const k_g = _process_data.thermal_conductivity_gas(t, pos)[0];
@@ -110,13 +129,10 @@ void HeatTransportBHELocalAssemblerSoil<
// for now only using the solid phase parameters
// assemble Conductance matrix
- local_K.noalias() += sm.dNdx.transpose() * k_s * sm.dNdx * sm.detJ *
- wp.getWeight() * sm.integralMeasure;
+ local_K.noalias() += K_ip * k_s;
// assemble Mass matrix
- local_M.noalias() += sm.N.transpose() * density_s * heat_capacity_s *
- sm.N * sm.detJ * wp.getWeight() *
- sm.integralMeasure;
+ local_M.noalias() += M_ip * density_s * heat_capacity_s;
}
// debugging
@@ -125,13 +141,5 @@ void HeatTransportBHELocalAssemblerSoil<
// std::cout << local_K.format(CleanFmt) << sep;
// std::cout << local_M.format(CleanFmt) << sep;
}
-
-template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
-void HeatTransportBHELocalAssemblerSoil<
- ShapeFunction,
- IntegrationMethod,
- GlobalDim>::postTimestepConcrete(std::vector<double> const& /*local_x*/)
-{
-}
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h
index 2af5e87a70b990a0cb760e2a32c51aa3055a188b..d0390c17999c0dc34459ea3f520c61e8381c5731 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h
@@ -26,38 +26,20 @@ class HeatTransportBHELocalAssemblerInterface
public NumLib::ExtrapolatableElement
{
public:
- HeatTransportBHELocalAssemblerInterface() : _dofIndex_to_localIndex{} {}
- HeatTransportBHELocalAssemblerInterface(
- std::size_t n_local_size, std::vector<unsigned> dofIndex_to_localIndex)
+ // The following function is necessary, because the BHE or Soil assemblers
+ // create their local_x memory based on the passed on dofIndex_to_localIndex
+ // Please keep it unchanged.
+ HeatTransportBHELocalAssemblerInterface(std::size_t /*n_local_size*/,
+ std::vector<unsigned>
+ dofIndex_to_localIndex)
: _dofIndex_to_localIndex(std::move(dofIndex_to_localIndex))
{
}
- void assembleWithJacobian(double const /*t*/,
- std::vector<double> const& /*local_x_*/,
- std::vector<double> const& /*local_xdot*/,
- const double /*dxdot_dx*/, const double /*dx_dx*/,
- std::vector<double>& /*local_M_data*/,
- std::vector<double>& /*local_K_data*/,
- std::vector<double>& /*local_b_data*/,
- std::vector<double>& /*local_Jac_data*/) override
- {
- OGS_FATAL(
- "HeatTransportBHELocalAssemblerInterface::assembleWithJacobian() "
- "is not implemented");
- }
-
- virtual void assembleWithJacobian(double const /*t*/,
- Eigen::VectorXd const& /*local_T*/,
- Eigen::VectorXd& /*local_b*/,
- Eigen::MatrixXd& /*local_A*/)
- {
- OGS_FATAL(
- "HeatTransportBHELocalAssemblerInterface::assembleWithJacobian() "
- "is not implemented");
- }
-
private:
+ // this dofTalbe must be kept here.
+ // When initializing the assembler, this will be needed to
+ // initialize the memory for local assemblers.
std::vector<unsigned> const _dofIndex_to_localIndex;
};
} // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
index 97abc53d0dc8d7a9091c545c44147f8aaa9db837..04e5c338af97f8b019a4369dd1ca19e54b2213d5 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
@@ -12,22 +12,15 @@
#include <memory>
#include <vector>
-#include "MaterialLib/SolidModels/MechanicsBase.h"
-
namespace ProcessLib
{
namespace HeatTransportBHE
{
-template <typename ShapeMatrixType>
+template <typename NodalMatrixType>
struct IntegrationPointDataSoil final
{
- explicit IntegrationPointDataSoil() {}
-
- typename ShapeMatrixType::NodalRowVectorType N;
- typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
- double integration_weight;
-
- void pushBackState() {}
+ NodalMatrixType NTN_product_times_w;
+ NodalMatrixType dNdxTdNdx_product_times_w;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
index f395d3551890211b73d693412b5a2fb014f50709..ac08d79591b912dbfd474cb4ece29860ecd57276 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
@@ -253,17 +253,17 @@ public:
type_idx.name());
}
- auto varIDs = _dof_table.getElementVariableIDs(id);
- auto n_local_dof = _dof_table.getNumberOfElementDOF(id);
std::vector<unsigned> dofIndex_to_localIndex;
- if (mesh_item.getDimension() < GlobalDim)
+ // Create customed dof table when it is a BHE element.
+ if (mesh_item.getDimension() == 1)
{
// this is a BHE element
+ auto const n_local_dof = _dof_table.getNumberOfElementDOF(id);
dofIndex_to_localIndex.resize(n_local_dof);
unsigned dof_id = 0;
unsigned local_id = 0;
- for (auto i : varIDs)
+ for (auto const i : _dof_table.getElementVariableIDs(id))
{
auto const n_global_components =
_dof_table.getNumberOfElementComponents(i);
@@ -276,7 +276,8 @@ public:
MeshLib::Location l(mesh_id,
MeshLib::MeshItemType::Node,
mesh_item.getNodeIndex(k));
- auto global_index = _dof_table.getGlobalIndex(l, i, j);
+ auto const global_index =
+ _dof_table.getGlobalIndex(l, i, j);
if (global_index != NumLib::MeshComponentMap::nop)
{
dofIndex_to_localIndex[dof_id++] = local_id;
@@ -287,16 +288,13 @@ public:
}
}
- data_ptr = it->second(mesh_item, varIDs.size(), n_local_dof,
- dofIndex_to_localIndex,
+ data_ptr = it->second(mesh_item, dofIndex_to_localIndex,
std::forward<ConstructorArgs>(args)...);
}
private:
using LADataBuilder = std::function<LADataIntfPtr(
MeshLib::Element const& e,
- std::size_t const n_variables,
- std::size_t const local_matrix_size,
std::vector<unsigned> const& dofIndex_to_localIndex,
ConstructorArgs&&...)>;
@@ -340,20 +338,43 @@ private:
static LADataBuilder makeLocalAssemblerBuilder(std::true_type*)
{
return [](MeshLib::Element const& e,
- std::size_t const /*n_variables*/,
- std::size_t const local_matrix_size,
std::vector<unsigned> const& dofIndex_to_localIndex,
- ConstructorArgs&&... args) {
+ ConstructorArgs&&... args) -> LADataIntfPtr {
if (e.getDimension() == GlobalDim) // soil elements
{
return LADataIntfPtr{new LADataSoil<ShapeFunction>{
- e, local_matrix_size, dofIndex_to_localIndex,
+ e, dofIndex_to_localIndex,
std::forward<ConstructorArgs>(args)...}};
}
- return LADataIntfPtr{new LADataBHE<ShapeFunction>{
+ return nullptr;
+ };
+ }
+
+ template <>
+ static LADataBuilder makeLocalAssemblerBuilder<NumLib::ShapeLine2>(
+ std::true_type*)
+ {
+ return [](MeshLib::Element const& e,
+ std::vector<unsigned> const& dofIndex_to_localIndex,
+ ConstructorArgs&&... args) -> LADataIntfPtr {
+ return LADataIntfPtr{new LADataBHE<NumLib::ShapeLine2>{
+ // BHE elements
+ e, dofIndex_to_localIndex,
+ std::forward<ConstructorArgs>(args)...}};
+ };
+ }
+
+ template <>
+ static LADataBuilder makeLocalAssemblerBuilder<NumLib::ShapeLine3>(
+ std::true_type*)
+ {
+ return [](MeshLib::Element const& e,
+ std::vector<unsigned> const& dofIndex_to_localIndex,
+ ConstructorArgs&&... args) -> LADataIntfPtr {
+ return LADataIntfPtr{new LADataBHE<NumLib::ShapeLine3>{
// BHE elements
- e, local_matrix_size, dofIndex_to_localIndex,
+ e, dofIndex_to_localIndex,
std::forward<ConstructorArgs>(args)...}};
};
}