diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 6c4950e92020aed7acce5db681f66b2aac712317..1f9f86beafbc864ef460bafbd686894a6f685453 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -53,6 +53,9 @@
#ifdef OGS_BUILD_PROCESS_HEATCONDUCTION
#include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
#endif
+#ifdef OGS_BUILD_PROCESS_HEATTRANSPORTBHE
+#include "ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.h"
+#endif
#ifdef OGS_BUILD_PROCESS_HYDROMECHANICS
#include "ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h"
#endif
@@ -371,6 +374,17 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
}
else
#endif
+#ifdef OGS_BUILD_PROCESS_HEATTRANSPORTBHE
+ if (type == "HEAT_TRANSPORT_BHE")
+ {
+ process =
+ ProcessLib::HeatTransportBHE::createHeatTransportBHEProcess(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config, _curves);
+ }
+ else
+#endif
#ifdef OGS_BUILD_PROCESS_HYDROMECHANICS
if (type == "HYDRO_MECHANICS")
{
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b7287fa978964c837fbdae946c768143d58c3da5..456c89de82e6053f7d3b94cba1b9d8cda8faa644 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -90,6 +90,7 @@ set(ProcessesList
GroundwaterFlow
HT
HeatConduction
+ HeatTransportBHE
HydroMechanics
LiquidFlow
LIE
diff --git a/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b87ec86a71c6c38e298bf4beb3ba1b2af4d726fe
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.cpp
@@ -0,0 +1,113 @@
+/**
+ * \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
+ *
+ */
+
+#include "BHEBottomDirichletBoundaryCondition.h"
+
+#include <algorithm>
+#include <logog/include/logog.hpp>
+#include <vector>
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+BHEBottomDirichletBoundaryCondition::BHEBottomDirichletBoundaryCondition(
+ GlobalIndexType const global_idx_T_in_bottom,
+ GlobalIndexType const global_idx_T_out_bottom,
+ MeshLib::Mesh const& bulk_mesh,
+ std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
+ int const variable_id,
+ int const component_id)
+ : _bulk_mesh(bulk_mesh)
+{
+ DBUG(
+ "Found %d nodes for BHE bottom Dirichlet BCs for the variable %d "
+ "and "
+ "component %d",
+ vec_outflow_bc_nodes.size(), variable_id, component_id);
+
+ MeshLib::MeshSubset bc_mesh_subset{_bulk_mesh, vec_outflow_bc_nodes};
+
+ // create memory to store Tout values
+ _T_in_values.ids.clear();
+ _T_in_values.values.clear();
+
+ _bc_values.ids.clear();
+ _bc_values.values.clear();
+
+ // convert mesh node ids to global index for the given component
+ _bc_values.ids.reserve(bc_mesh_subset.getNumberOfNodes());
+ _bc_values.values.reserve(bc_mesh_subset.getNumberOfNodes());
+
+ const auto g_T_out_idx = global_idx_T_out_bottom;
+ if (g_T_out_idx >= 0)
+ {
+ _bc_values.ids.emplace_back(g_T_out_idx);
+ _bc_values.values.emplace_back(298.15);
+ }
+
+ const auto g_T_in_idx = global_idx_T_in_bottom;
+
+ if (g_T_in_idx >= 0)
+ {
+ _T_in_values.ids.emplace_back(g_T_in_idx);
+ _T_in_values.values.emplace_back(298.15);
+ }
+}
+
+void BHEBottomDirichletBoundaryCondition::getEssentialBCValues(
+ const double /*t*/, GlobalVector const& x,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+{
+ bc_values.ids.clear();
+ bc_values.values.clear();
+
+ bc_values.ids.resize(_bc_values.ids.size());
+ bc_values.values.resize(_bc_values.values.size());
+
+ const std::size_t n_nodes = _T_in_values.ids.size();
+ for (std::size_t i = 0; i < n_nodes; i++)
+ {
+ bc_values.ids[i] = _bc_values.ids[i];
+ // here, the outflow temperature is always
+ // the same as the inflow temperature
+ // get the inflow temperature from here.
+ bc_values.values[i] = x[_T_in_values.ids[i]];
+ }
+}
+
+// update new values and corresponding indices.
+void BHEBottomDirichletBoundaryCondition::preTimestep(const double /*t*/,
+ const GlobalVector& x)
+{
+ // At the bottom of each BHE, the outflow temperature
+ // is the same as the inflow temperature.
+ // Here the task is to get the inflow temperature and
+ // save it locally
+ auto const n_nodes = _bc_values.ids.size();
+ for (std::size_t i = 0; i < n_nodes; i++)
+ {
+ // read the T_out
+ _T_in_values.values[i] = x[_T_in_values.ids[i]];
+ }
+}
+
+std::unique_ptr<BHEBottomDirichletBoundaryCondition>
+createBHEBottomDirichletBoundaryCondition(
+ GlobalIndexType global_idx_T_in_bottom,
+ GlobalIndexType global_idx_T_out_bottom, MeshLib::Mesh const& bulk_mesh,
+ std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
+ int const variable_id, int const component_id)
+{
+ DBUG("Constructing BHEBottomDirichletBoundaryCondition from config.");
+
+ return std::make_unique<BHEBottomDirichletBoundaryCondition>(
+ global_idx_T_in_bottom, global_idx_T_out_bottom, bulk_mesh,
+ vec_outflow_bc_nodes, variable_id, component_id);
+}
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..063eac02a09b4d864baea5aea30f7e9796e82529
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h
@@ -0,0 +1,51 @@
+/**
+ * \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
+
+#include "BoundaryCondition.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/IndexValueVector.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace ProcessLib
+{
+class BHEBottomDirichletBoundaryCondition final : public BoundaryCondition
+{
+public:
+ BHEBottomDirichletBoundaryCondition(
+ GlobalIndexType global_idx_T_in_bottom,
+ GlobalIndexType global_idx_T_out_bottom,
+ MeshLib::Mesh const& bulk_mesh,
+ std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
+ int const variable_id,
+ int const component_id);
+
+ void getEssentialBCValues(
+ const double t, GlobalVector const& x,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
+
+ void preTimestep(const double t, const GlobalVector& x) override;
+
+private:
+ MeshLib::Mesh const& _bulk_mesh;
+
+ NumLib::IndexValueVector<GlobalIndexType> _bc_values;
+
+ NumLib::IndexValueVector<GlobalIndexType> _T_in_values;
+};
+
+std::unique_ptr<BHEBottomDirichletBoundaryCondition>
+createBHEBottomDirichletBoundaryCondition(
+ GlobalIndexType global_idx_T_in_bottom,
+ GlobalIndexType global_idx_T_out_bottom, MeshLib::Mesh const& bulk_mesh,
+ std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
+ int const variable_id, int const component_id);
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b152909b53e648ed4ba08001d252996e21153086
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
@@ -0,0 +1,117 @@
+/**
+ * \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
+ *
+ */
+
+#include "BHEInflowDirichletBoundaryCondition.h"
+
+#include <algorithm>
+#include <logog/include/logog.hpp>
+#include <vector>
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
+ GlobalIndexType global_idx_T_in_top,
+ GlobalIndexType global_idx_T_out_top,
+ MeshLib::Mesh const& bc_mesh,
+ std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
+ int const variable_id,
+ int const component_id,
+ std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
+ pt_bhe)
+ : _bc_mesh(bc_mesh), _pt_bhe(pt_bhe)
+{
+ DBUG(
+ "Found %d nodes for BHE Inflow Dirichlet BCs for the variable %d "
+ "and "
+ "component %d",
+ vec_inflow_bc_nodes.size(), variable_id, component_id);
+
+ MeshLib::MeshSubset bc_mesh_subset{_bc_mesh, vec_inflow_bc_nodes};
+
+ // create memory to store Tout values
+ _T_out_values.clear();
+ _T_out_indices.clear();
+
+ _bc_values.ids.clear();
+ _bc_values.values.clear();
+
+ // convert mesh node ids to global index for the given component
+ assert(bc_mesh_subset.getNumberOfNodes() == 1);
+ _bc_values.ids.reserve(bc_mesh_subset.getNumberOfNodes());
+ _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;
+
+ if (g_idx_T_in >= 0 && g_idx_T_out >= 0)
+ {
+ _T_out_indices.emplace_back(g_idx_T_out);
+ _T_out_values.emplace_back(320.0 /*using initial value*/);
+ _bc_values.ids.emplace_back(g_idx_T_in);
+ _bc_values.values.emplace_back(320.0 /*using initial value*/);
+ }
+}
+
+void BHEInflowDirichletBoundaryCondition::getEssentialBCValues(
+ const double t, GlobalVector const& x,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+{
+ bc_values.ids.clear();
+ bc_values.values.clear();
+
+ bc_values.ids.resize(_bc_values.ids.size());
+ bc_values.values.resize(_bc_values.values.size());
+
+ const std::size_t n_nodes = _T_out_values.size();
+
+ for (std::size_t i = 0; i < n_nodes; i++)
+ {
+ bc_values.ids[i] = _bc_values.ids[i];
+ // here call the corresponding BHE functions
+ auto const tmp_T_out = x[_T_out_indices[i]];
+ bc_values.values[i] = _pt_bhe->getTinByTout(tmp_T_out, t);
+ }
+}
+
+// update new values and corresponding indices.
+void BHEInflowDirichletBoundaryCondition::preTimestep(const double /*t*/,
+ const GlobalVector& x)
+{
+ // for each BHE, the inflow temperature is dependent on
+ // the ouflow temperature of the BHE.
+ // Here the task is to get the outflow temperature and
+ // save it locally
+ auto const n_nodes = _bc_values.ids.size();
+ for (std::size_t i = 0; i < n_nodes; i++)
+ {
+ // read the T_out
+ _T_out_values[i] = x[_T_out_indices[i]];
+ }
+}
+
+std::unique_ptr<BHEInflowDirichletBoundaryCondition>
+createBHEInflowDirichletBoundaryCondition(
+ GlobalIndexType global_idx_T_in_top, GlobalIndexType global_idx_T_out_top,
+ MeshLib::Mesh const& bc_mesh,
+ std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
+ int const variable_id, int const component_id,
+ std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
+ pt_bhe)
+{
+ 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,
+ variable_id, component_id, pt_bhe);
+}
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..53c2afc0f4f784cfcee5c674e98f5128bf83ac29
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
@@ -0,0 +1,60 @@
+/**
+ * \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
+
+#include "BoundaryCondition.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/IndexValueVector.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace ProcessLib
+{
+class BHEInflowDirichletBoundaryCondition final : public BoundaryCondition
+{
+public:
+ BHEInflowDirichletBoundaryCondition(
+ GlobalIndexType global_idx_T_in_top,
+ GlobalIndexType global_idx_T_out_top,
+ MeshLib::Mesh const& bc_mesh,
+ std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
+ int const variable_id,
+ int const component_id,
+ std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
+ pt_bhe);
+
+ void getEssentialBCValues(
+ const double t, GlobalVector const& x,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
+
+ void preTimestep(const double t, const GlobalVector& x) override;
+
+private:
+ MeshLib::Mesh const& _bc_mesh;
+
+ /// Stores the results of the outflow temperatures per boundary node.
+ std::vector<double> _T_out_values;
+ std::vector<GlobalIndexType> _T_out_indices;
+
+ NumLib::IndexValueVector<GlobalIndexType> _bc_values;
+
+ std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
+ _pt_bhe;
+};
+
+std::unique_ptr<BHEInflowDirichletBoundaryCondition>
+createBHEInflowDirichletBoundaryCondition(
+ GlobalIndexType global_idx_T_in_top, GlobalIndexType global_idx_T_out_top,
+ MeshLib::Mesh const& bc_mesh,
+ std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
+ int const variable_id, int const component_id,
+ std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
+ pt_bhe);
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1a55cc0c94d3f7f0588cff3c1c901db71b8e3b93
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
@@ -0,0 +1,101 @@
+/**
+ * \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
+ */
+
+#include "MeshUtils.h"
+
+#include "BaseLib/Algorithm.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshSearch/NodeSearch.h"
+#include "MeshLib/Node.h"
+
+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)
+{
+ // get vectors of matrix elements and BHE elements
+ vec_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);
+ }
+
+ // get BHE material IDs
+ auto opt_material_ids(
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
+ 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());
+
+ // 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++)
+ {
+ const auto bhe_mat_id = vec_BHE_mat_IDs[bhe_id];
+ std::vector<MeshLib::Element*>& vec_elements = vec_BHE_elements[bhe_id];
+ std::copy_if(all_BHE_elements.begin(), all_BHE_elements.end(),
+ std::back_inserter(vec_elements),
+ [&](MeshLib::Element* e) {
+ return (*opt_material_ids)[e->getID()] == bhe_mat_id;
+ });
+ DBUG("-> found %d elements on the BHE_%d", vec_elements.size(), bhe_id);
+ }
+
+ // 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++)
+ {
+ std::vector<MeshLib::Node*>& vec_nodes = vec_BHE_nodes[bhe_id];
+ for (MeshLib::Element* e : vec_BHE_elements[bhe_id])
+ {
+ for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
+ {
+ vec_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
+ }
+ }
+ BaseLib::makeVectorUnique(
+ vec_nodes, [](MeshLib::Node* node1, MeshLib::Node* node2) {
+ return node1->getID() < node2->getID();
+ });
+ DBUG("-> found %d nodes on the BHE_%d", vec_nodes.size(), bhe_id);
+ }
+
+ // get all the nodes into the pure soil nodes vector
+ vec_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);
+ }
+
+ // final count of 3 types of elements
+ // They are
+ // (i) soil,
+ // (ii) BHE
+ DBUG("-> found total %d soil elements and %d BHE elements",
+ vec_soil_elements.size(),
+ all_BHE_elements.size());
+}
+} // end of namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac75431d2e0057279b536fc2a13ba6b2f4d32154
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
@@ -0,0 +1,42 @@
+/**
+ * \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
+
+#include <vector>
+
+namespace MeshLib
+{
+class Element;
+class Mesh;
+class Node;
+} // namespace MeshLib
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+/**
+ * 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);
+} // end of namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/CMakeLists.txt b/ProcessLib/HeatTransportBHE/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2d9def64402e62d9154925303a0ade8532756dab
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/CMakeLists.txt
@@ -0,0 +1,8 @@
+APPEND_SOURCE_FILES(SOURCES)
+APPEND_SOURCE_FILES(SOURCES BHE)
+APPEND_SOURCE_FILES(SOURCES LocalAssemblers)
+
+add_library(HeatTransportBHE ${SOURCES})
+target_link_libraries(HeatTransportBHE PUBLIC ProcessLib)
+
+# include(Tests.cmake)
diff --git a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ece955939dc2f849f42eeb98efa69634d8d2c30d
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
@@ -0,0 +1,306 @@
+/**
+ * \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
+ *
+ */
+
+#include "CreateHeatTransportBHEProcess.h"
+
+#include "BHE/BHEAbstract.h"
+#include "BHE/BHE_1U.h"
+#include "BHE/BHE_2U.h"
+#include "BHE/BHE_CXA.h"
+#include "BHE/BHE_CXC.h"
+#include "BHE/CreateBHE1U.h"
+#include "BHE/CreateBHE2U.h"
+#include "BHE/CreateBHECXA.h"
+#include "BHE/CreateBHECXC.h"
+#include "BaseLib/Algorithm.h"
+#include "HeatTransportBHEProcess.h"
+#include "HeatTransportBHEProcessData.h"
+#include "MaterialLib/Fluid/Density/CreateFluidDensityModel.h"
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/Fluid/SpecificHeatCapacity/CreateSpecificFluidHeatCapacityModel.h"
+#include "MaterialLib/Fluid/ThermalConductivity/CreateFluidThermalConductivityModel.h"
+#include "MaterialLib/Fluid/Viscosity/CreateViscosityModel.h"
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+std::unique_ptr<Process> createHeatTransportBHEProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves)
+{
+ //! \ogs_file_param{prj__processes__process__type}
+ config.checkConfigParameter("type", "HEAT_TRANSPORT_BHE");
+
+ DBUG("Create HeatTransportBHE Process.");
+
+ // Process variable.
+
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+ process_variables;
+
+ // reading primary variables for each
+ // BHE----------------------------------------------------------
+ auto range =
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__process_variables__process_variable}
+ pv_config.getConfigParameterList<std::string>("process_variable");
+ std::vector<std::reference_wrapper<ProcessVariable>> per_process_variables;
+
+ for (std::string const& pv_name : range)
+ {
+ if (pv_name != "temperature_soil" &&
+ pv_name.find("temperature_BHE") != 0)
+ {
+ OGS_FATAL(
+ "Found a process variable name '%s'. It should be "
+ "'temperature_soil' or 'temperature_BHE_X'");
+ }
+ auto variable = std::find_if(variables.cbegin(), variables.cend(),
+ [&pv_name](ProcessVariable const& v) {
+ return v.getName() == pv_name;
+ });
+
+ if (variable == variables.end())
+ {
+ OGS_FATAL(
+ "Could not find process variable '%s' in the provided "
+ "variables "
+ "list for config tag <%s>.",
+ pv_name.c_str(), "process_variable");
+ }
+ DBUG("Found process variable \'%s\' for config tag <%s>.",
+ variable->getName().c_str(), "process_variable");
+
+ per_process_variables.emplace_back(
+ const_cast<ProcessVariable&>(*variable));
+ }
+ process_variables.push_back(std::move(per_process_variables));
+ // end of reading primary variables for each
+ // BHE----------------------------------------------------------
+
+ // solid phase thermal conductivity parameter.
+ auto& thermal_conductivity_solid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special
+ //! prj__processes__process__HEAT_TRANSPORT_BHE__thermal_conductivity_solid}
+ "thermal_conductivity_solid", parameters, 1);
+
+ DBUG("Use \'%s\' as solid phase thermal conductivity parameter.",
+ thermal_conductivity_solid.name.c_str());
+
+ // solid phase thermal conductivity parameter.
+ auto& thermal_conductivity_fluid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__thermal_conductivity_fluid}
+ "thermal_conductivity_fluid", parameters, 1);
+
+ DBUG("Use \'%s\' as fluid phase thermal conductivity parameter.",
+ thermal_conductivity_fluid.name.c_str());
+
+ // gas phase thermal conductivity parameter.
+ auto& thermal_conductivity_gas = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__thermal_conductivity_gas}
+ "thermal_conductivity_gas", parameters, 1);
+
+ DBUG("Use \'%s\' as gas phase thermal conductivity parameter.",
+ thermal_conductivity_gas.name.c_str());
+
+ // solid phase heat capacity parameter.
+ auto& heat_capacity_solid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__heat_capacity_solid}
+ "heat_capacity_solid", parameters, 1);
+
+ DBUG("Use \'%s\' as solid phase heat capacity parameter.",
+ heat_capacity_solid.name.c_str());
+
+ // fluid phase heat capacity parameter.
+ auto& heat_capacity_fluid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__heat_capacity_fluid}
+ "heat_capacity_fluid", parameters, 1);
+
+ DBUG("Use \'%s\' as fluid phase heat capacity parameter.",
+ heat_capacity_fluid.name.c_str());
+
+ // gas phase heat capacity parameter.
+ auto& heat_capacity_gas = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__heat_capacity_gas}
+ "heat_capacity_gas", parameters, 1);
+
+ DBUG("Use \'%s\' as gas phase heat capacity parameter.",
+ heat_capacity_gas.name.c_str());
+
+ // solid phase density parameter.
+ auto& density_solid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__density_solid}
+ "density_solid", parameters, 1);
+
+ DBUG("Use \'%s\' as solid phase density parameter.",
+ density_solid.name.c_str());
+
+ // fluid phase density parameter.
+ auto& density_fluid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__density_fluid}
+ "density_fluid", parameters, 1);
+
+ DBUG("Use \'%s\' as fluid phase density parameter.",
+ density_fluid.name.c_str());
+
+ // gas phase density parameter.
+ auto& density_gas = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__HEAT_TRANSPORT_BHE__density_gas}
+ "density_gas", parameters, 1);
+
+ DBUG("Use \'%s\' as gas phase density parameter.",
+ density_gas.name.c_str());
+
+ // get the refrigerant properties from fluid property class
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__fluid}
+ auto const& fluid_config = config.getConfigSubtree("fluid");
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_density}
+ auto const& rho_conf = fluid_config.getConfigSubtree("refrigerant_density");
+ auto bhe_refrigerant_density =
+ MaterialLib::Fluid::createFluidDensityModel(rho_conf);
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_viscosity}
+ auto const& mu_conf =
+ fluid_config.getConfigSubtree("refrigerant_viscosity");
+ auto bhe_refrigerant_viscosity =
+ MaterialLib::Fluid::createViscosityModel(mu_conf);
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_specific_heat_capacity}
+ auto const& cp_conf =
+ fluid_config.getConfigSubtree("refrigerant_specific_heat_capacity");
+ auto bhe_refrigerant_heat_capacity =
+ MaterialLib::Fluid::createSpecificFluidHeatCapacityModel(cp_conf);
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_thermal_conductivity}
+ auto const& lambda_conf =
+ fluid_config.getConfigSubtree("refrigerant_thermal_conductivity");
+ auto bhe_regrigerant_heat_conductivity =
+ MaterialLib::Fluid::createFluidThermalConductivityModel(lambda_conf);
+
+ // reading BHE
+ // parameters--------------------------------------------------------------
+ std::vector<std::unique_ptr<BHE::BHEAbstract>> vec_BHEs;
+ // BHE::BHE_Net BHE_network;
+
+ // now read the BHE configurations
+ auto const& bhe_configs =
+ //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers}
+ config.getConfigSubtree("borehole_heat_exchangers");
+
+ for (
+ auto const& bhe_conf :
+ //! \ogs_file_param
+ //! prj__processes__process___HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger}
+ bhe_configs.getConfigSubtreeList("borehole_heat_exchanger"))
+ {
+ // read in the parameters
+ using namespace BHE;
+ const std::string bhe_type_str =
+ bhe_conf.getConfigParameter<std::string>("bhe_type");
+
+ // convert BHE type
+ if (bhe_type_str == "BHE_TYPE_1U")
+ {
+ // initialize the 1U type BHE
+ BHE::BHE_1U* m_bhe_1u =
+ BHE::CreateBHE1U(bhe_conf,
+ curves,
+ bhe_refrigerant_density,
+ bhe_refrigerant_viscosity,
+ bhe_refrigerant_heat_capacity,
+ bhe_regrigerant_heat_conductivity);
+
+ vec_BHEs.emplace_back(std::make_unique<BHE_1U>(*m_bhe_1u));
+ }
+ else if (bhe_type_str == "BHE_TYPE_2U")
+ {
+ // initialize the 2U type BHE
+ BHE::BHE_2U* m_bhe_2u =
+ BHE::CreateBHE2U(bhe_conf,
+ curves,
+ bhe_refrigerant_density,
+ bhe_refrigerant_viscosity,
+ bhe_refrigerant_heat_capacity,
+ bhe_regrigerant_heat_conductivity);
+
+ vec_BHEs.emplace_back(std::make_unique<BHE_2U>(*m_bhe_2u));
+ }
+ else if (bhe_type_str == "BHE_TYPE_CXC")
+ {
+ // initialize the CXC type BHE
+ BHE::BHE_CXC* m_bhe_CXC =
+ BHE::CreateBHECXC(bhe_conf,
+ curves,
+ bhe_refrigerant_density,
+ bhe_refrigerant_viscosity,
+ bhe_refrigerant_heat_capacity,
+ bhe_regrigerant_heat_conductivity);
+
+ vec_BHEs.emplace_back(std::make_unique<BHE_CXC>(*m_bhe_CXC));
+ }
+ else if (bhe_type_str == "BHE_TYPE_CXA")
+ {
+ // initialize the CXA type BHE
+ BHE::BHE_CXA* m_bhe_CXA =
+ BHE::CreateBHECXA(bhe_conf,
+ curves,
+ bhe_refrigerant_density,
+ bhe_refrigerant_viscosity,
+ bhe_refrigerant_heat_capacity,
+ bhe_regrigerant_heat_conductivity);
+
+ vec_BHEs.emplace_back(std::make_unique<BHE_CXA>(*m_bhe_CXA));
+ }
+ }
+ // end of reading BHE
+ // parameters-------------------------------------------------------
+
+ HeatTransportBHEProcessData process_data{thermal_conductivity_solid,
+ thermal_conductivity_fluid,
+ thermal_conductivity_gas,
+ heat_capacity_solid,
+ heat_capacity_fluid,
+ heat_capacity_gas,
+ density_solid,
+ density_fluid,
+ density_gas,
+ std::move(vec_BHEs)};
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"HeatTransportBHE_Temperature"});
+
+ ProcessLib::createSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+
+ return std::make_unique<HeatTransportBHEProcess>(
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(process_data),
+ std::move(secondary_variables), std::move(named_function_caller));
+}
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.h b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..e6ee4f3d9db41ec449cc3313f2cd238f4b65d1cd
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.h
@@ -0,0 +1,30 @@
+/**
+ * \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
+
+#include <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+std::unique_ptr<Process> createHeatTransportBHEProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves);
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a8f1a0c8a3477dad8ee81aa169949e6534a57b08
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
@@ -0,0 +1,419 @@
+/**
+ * \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
+ *
+ */
+
+#include "HeatTransportBHEProcess.h"
+
+#include <cassert>
+#include <iostream>
+#include "ProcessLib/HeatTransportBHE/BHE/MeshUtils.h"
+#include "ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h"
+
+#include "ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h"
+#include "ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h"
+
+#include "ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h"
+#include "ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+HeatTransportBHEProcess::HeatTransportBHEProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ HeatTransportBHEProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller)
+ : 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))
+{
+ 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())
+ {
+ 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());
+ }
+
+ // 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++)
+ {
+ // 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;
+ }
+
+ MeshLib::PropertyVector<int> const* material_ids(
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
+ _process_data._mesh_prop_materialIDs = material_ids;
+}
+
+void HeatTransportBHEProcess::constructDofTable()
+{
+ // Create single component dof in every of the mesh's nodes.
+ _mesh_subset_all_nodes =
+ std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
+
+ //------------------------------------------------------------
+ // prepare mesh subsets to define DoFs
+ //------------------------------------------------------------
+ // first all the soil nodes
+ _mesh_subset_pure_soil_nodes =
+ std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_pure_soil_nodes);
+
+ std::vector<int> vec_n_BHE_unknowns;
+ vec_n_BHE_unknowns.reserve(_vec_BHE_nodes.size());
+ // the BHE nodes need to be cherry-picked from the vector
+ for (unsigned i = 0; i < _vec_BHE_nodes.size(); i++)
+ {
+ _mesh_subset_BHE_nodes.push_back(
+ std::make_unique<MeshLib::MeshSubset const>(_mesh,
+ _vec_BHE_nodes[i]));
+ int n_BHE_unknowns =
+ _process_data._vec_BHE_property[i]->getNumUnknowns();
+ vec_n_BHE_unknowns.emplace_back(n_BHE_unknowns);
+ }
+
+ // Collect the mesh subsets in a vector.
+ // All the soil nodes has 1 temperature variable
+ std::vector<MeshLib::MeshSubset> all_mesh_subsets{
+ *_mesh_subset_pure_soil_nodes};
+
+ // All the BHE nodes have additinal variables
+ std::size_t count = 0;
+ for (auto& ms : _mesh_subset_BHE_nodes)
+ {
+ std::generate_n(std::back_inserter(all_mesh_subsets),
+ // Here the number of components equals to
+ // the number of unknowns on the BHE
+ vec_n_BHE_unknowns[count],
+ [&]() { return *ms; });
+ count++;
+ }
+
+ std::vector<int> vec_n_components;
+ // This is the soil temperature for first mesh subset.
+ // 1, because for the soil part there is just one variable which is the soil
+ // temperature.
+ vec_n_components.push_back(1);
+ // now the BHE subsets
+ for (std::size_t i = 0; i < _vec_BHE_mat_IDs.size(); i++)
+ {
+ // Here the number of components equals to
+ // the number of unknowns on the BHE
+ vec_n_components.push_back(vec_n_BHE_unknowns[i]);
+ }
+
+ 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++)
+ {
+ vec_var_elements.push_back(&_vec_BHE_elements[i]);
+ }
+
+ _local_to_global_index_map =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets),
+ vec_n_components,
+ vec_var_elements,
+ NumLib::ComponentOrder::BY_COMPONENT);
+
+ // in case of debugging the dof table, activate the following line
+ // std::cout << *_local_to_global_index_map << "\n";
+}
+
+void HeatTransportBHEProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ const int process_id = 0;
+
+ // this process can only run with 3-dimensional mesh
+ ProcessLib::HeatTransportBHE::createLocalAssemblers<
+ 3, /*mesh.getDimension(),*/
+ HeatTransportBHELocalAssemblerSoil, HeatTransportBHELocalAssemblerBHE>(
+ mesh.getElements(), dof_table, _local_assemblers,
+ _process_data._vec_ele_connected_BHE_IDs,
+ _process_data._vec_BHE_property, mesh.isAxiallySymmetric(),
+ integration_order, _process_data);
+
+ // create BHE boundary conditions
+ // for each BHE, one BC on the top
+ // and one BC at the bottom.
+ std::vector<std::unique_ptr<BoundaryCondition>> bc_collections =
+ createBHEBoundaryConditionTopBottom();
+ 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]));
+}
+
+void HeatTransportBHEProcess::assembleConcreteProcess(const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble 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::assemble, _local_assemblers,
+ dof_table, t, x, M, K, b, _coupled_solutions);
+}
+
+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)
+{
+ 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);
+}
+
+void HeatTransportBHEProcess::computeSecondaryVariableConcrete(
+ const double t, GlobalVector const& x)
+{
+ DBUG("Compute heat flux for HeatTransportBHE process.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &HeatTransportBHELocalAssemblerInterface::computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, t, x,
+ _coupled_solutions);
+}
+
+std::vector<std::unique_ptr<BoundaryCondition>>
+HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom()
+{
+ /**
+ * boundary conditions
+ */
+ std::vector<std::unique_ptr<BoundaryCondition>> bcs;
+
+ 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;
+ // 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)
+ {
+ 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");
+ }
+ }
+
+ return bcs;
+}
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..217e4da211538bd1e1b2a59f75c329ed2a4d215e
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
@@ -0,0 +1,98 @@
+/**
+ * \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
+
+#include "HeatTransportBHEProcessData.h"
+#include "ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h"
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+class HeatTransportBHEProcess final : public Process
+{
+public:
+ HeatTransportBHEProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+ jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ HeatTransportBHEProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller);
+
+ //! \name ODESystem interface
+ //! @{
+ bool isLinear() const override { return false; }
+
+ void computeSecondaryVariableConcrete(double const t,
+ GlobalVector const& x) override;
+
+private:
+ void constructDofTable() override;
+
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override;
+
+ void assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) override;
+
+ void 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) override;
+
+ std::vector<std::unique_ptr<BoundaryCondition>>
+ createBHEBoundaryConditionTopBottom();
+
+ HeatTransportBHEProcessData _process_data;
+
+ 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;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..28224822f0ecdfa0bc48b733deaa84ad6822011a
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h
@@ -0,0 +1,90 @@
+/**
+ * \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
+
+#include "MeshLib/PropertyVector.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+template <typename T>
+struct Parameter;
+
+namespace HeatTransportBHE
+{
+using namespace BHE;
+
+struct HeatTransportBHEProcessData
+{
+ HeatTransportBHEProcessData(
+ Parameter<double> const& thermal_conductivity_solid_,
+ Parameter<double> const& thermal_conductivity_fluid_,
+ Parameter<double> const& thermal_conductivity_gas_,
+ Parameter<double> const& heat_capacity_solid_,
+ Parameter<double> const& heat_capacity_fluid_,
+ Parameter<double> const& heat_capacity_gas_,
+ Parameter<double> const& density_solid_,
+ Parameter<double> const& density_fluid_,
+ Parameter<double> const& density_gas_,
+ std::vector<std::unique_ptr<BHE::BHEAbstract>>&& vec_BHEs_)
+ : thermal_conductivity_solid(thermal_conductivity_solid_),
+ thermal_conductivity_fluid(thermal_conductivity_fluid_),
+ thermal_conductivity_gas(thermal_conductivity_gas_),
+ heat_capacity_solid(heat_capacity_solid_),
+ heat_capacity_fluid(heat_capacity_fluid_),
+ heat_capacity_gas(heat_capacity_gas_),
+ density_solid(density_solid_),
+ density_fluid(density_fluid_),
+ density_gas(density_gas_),
+ _vec_BHE_property(std::move(vec_BHEs_))
+ {
+ }
+
+ HeatTransportBHEProcessData(HeatTransportBHEProcessData&& other) = default;
+
+ //! Copies are forbidden.
+ HeatTransportBHEProcessData(HeatTransportBHEProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(HeatTransportBHEProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(HeatTransportBHEProcessData&&) = delete;
+
+ // ! thermal conductivity values for the three phases
+ Parameter<double> const& thermal_conductivity_solid;
+ Parameter<double> const& thermal_conductivity_fluid;
+ Parameter<double> const& thermal_conductivity_gas;
+
+ // ! heat capacity values for the three phases
+ Parameter<double> const& heat_capacity_solid;
+ Parameter<double> const& heat_capacity_fluid;
+ Parameter<double> const& heat_capacity_gas;
+
+ // ! density values for the three phases
+ Parameter<double> const& density_solid;
+ Parameter<double> const& density_fluid;
+ Parameter<double> const& density_gas;
+
+ MeshLib::PropertyVector<int> const* _mesh_prop_materialIDs = nullptr;
+ std::vector<std::size_t> _map_materialID_to_BHE_ID;
+
+ std::vector<std::unique_ptr<BHE::BHEAbstract>> _vec_BHE_property;
+
+ // a table of connected BHE IDs for each element
+ std::vector<std::vector<std::size_t>> _vec_ele_connected_BHE_IDs;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h
new file mode 100644
index 0000000000000000000000000000000000000000..03f6444c70de73ac825591fc99951107afb98d8c
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h
@@ -0,0 +1,86 @@
+/**
+ * \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
+
+#include <vector>
+
+#include <logog/include/logog.hpp>
+
+#include "LocalDataInitializer.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace detail
+{
+template <
+ int GlobalDim,
+ template <typename, typename, int> class LocalAssemblerSoilImplementation,
+ template <typename, typename, int> class LocalAssemblerBHEImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ const std::vector<std::vector<std::size_t>>& vec_ele_connected_BHE_IDs,
+ const std::vector<std::unique_ptr<BHE::BHEAbstract>>& vec_BHE_property,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ // Shape matrices initializer
+ using LocalDataInitializer = LocalDataInitializer<
+ LocalAssemblerInterface, LocalAssemblerSoilImplementation,
+ LocalAssemblerBHEImplementation, 3, ExtraCtorArgs...>;
+
+ DBUG("Create local assemblers for the HeatTransportBHE process.");
+ // Populate the vector of local assemblers.
+ local_assemblers.resize(mesh_elements.size());
+
+ LocalDataInitializer initializer(dof_table);
+
+ DBUG("Calling local assembler builder for all mesh elements.");
+ GlobalExecutor::transformDereferenced(
+ initializer, mesh_elements, local_assemblers, vec_ele_connected_BHE_IDs,
+ vec_BHE_property, std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+} // namespace detail
+
+/*! Creates local assemblers for each element of the given \c mesh.
+ *
+ * \tparam LocalAssemblerImplementation the individual local assembler type
+ * \tparam LocalAssemblerInterface the general local assembler interface
+ * \tparam ExtraCtorArgs types of additional constructor arguments.
+ * Those arguments will be passed to the constructor of
+ * \c LocalAssemblerImplementation.
+ *
+ * The first two template parameters cannot be deduced from the arguments.
+ * Therefore they always have to be provided manually.
+ */
+template <
+ int GlobalDim,
+ template <typename, typename, int> class LocalAssemblerSoilImplementation,
+ template <typename, typename, int> class LocalAssemblerBHEImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ DBUG("Create local assemblers for the HeatTransportBHE process.");
+
+ detail::createLocalAssemblers<GlobalDim, LocalAssemblerSoilImplementation,
+ LocalAssemblerBHEImplementation>(
+ dof_table, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
new file mode 100644
index 0000000000000000000000000000000000000000..f7ea07d91062fd54fbce295eec09fb9b6a299b76
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
@@ -0,0 +1,106 @@
+/**
+ * \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
+
+#include <vector>
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+#include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h"
+
+#include "HeatTransportBHEProcessAssemblerInterface.h"
+#include "IntegrationPointDataBHE.h"
+#include "SecondaryData.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
+class HeatTransportBHELocalAssemblerBHE
+ : public HeatTransportBHELocalAssemblerInterface
+{
+public:
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, BHE_Dim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+ // Using dynamic size, because the number of unknowns in the BHE is runtime
+ // value.
+ using BheLocalMatrixType =
+ typename ShapeMatricesType::template MatrixType<Eigen::Dynamic,
+ Eigen::Dynamic>;
+ HeatTransportBHELocalAssemblerBHE(
+ HeatTransportBHELocalAssemblerBHE const&) = delete;
+ HeatTransportBHELocalAssemblerBHE(HeatTransportBHELocalAssemblerBHE&&) =
+ delete;
+
+ 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,
+ HeatTransportBHEProcessData& process_data);
+
+ void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
+ std::vector<double>& /*local_M_data*/,
+ 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
+ {
+ }
+
+ void postTimestepConcrete(std::vector<double> const& /*local_x*/) override;
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _secondary_data.N[integration_point];
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+private:
+ HeatTransportBHEProcessData& _process_data;
+
+ std::vector<
+ IntegrationPointDataBHE<ShapeMatricesType>,
+ Eigen::aligned_allocator<IntegrationPointDataBHE<ShapeMatricesType>>>
+ _ip_data;
+
+ IntegrationMethod _integration_method;
+
+ std::size_t const element_id;
+
+ SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+
+ BheLocalMatrixType _R_matrix;
+
+ BheLocalMatrixType _R_s_matrix;
+
+ BheLocalMatrixType _R_pi_s_matrix;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
+
+#include "HeatTransportBHELocalAssemblerBHE_impl.h"
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..5285b4b2785a62419eda825da18d9d556ea7fd19
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
@@ -0,0 +1,226 @@
+/**
+ * \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
+
+#include <Eigen/Eigen>
+#include "HeatTransportBHELocalAssemblerBHE.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+using namespace BHE;
+
+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,
+ HeatTransportBHEProcessData& process_data)
+ : HeatTransportBHELocalAssemblerInterface(
+ ShapeFunction::NPOINTS * BHE_Dim, // no intersection
+ dofIndex_to_localIndex),
+ _process_data(process_data),
+ _integration_method(integration_order),
+ element_id(e.getID())
+{
+ // need to make sure that the BHE elements are one-dimensional
+ assert(e.getDimension() == 1);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N.resize(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
+ BHE_Dim>(e, is_axially_symmetric,
+ _integration_method);
+
+ auto mat_id = (*_process_data._mesh_prop_materialIDs)[e.getID()];
+ auto BHE_id = _process_data._map_materialID_to_BHE_ID[mat_id];
+
+ SpatialPosition x_position;
+ x_position.setElementID(element_id);
+
+ // ip data initialization
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ 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);
+ auto const& sm = shape_matrices[ip];
+ auto& ip_data = _ip_data[ip];
+ ip_data.integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ sm.integralMeasure * sm.detJ;
+ ip_data.N = sm.N;
+ ip_data.dNdx = sm.dNdx;
+
+ _secondary_data.N[ip] = sm.N;
+ }
+
+ const int BHE_n_unknowns = _ip_data[0]._bhe_instance.getNumUnknowns();
+ const int NumRMatrixRows = ShapeFunction::NPOINTS * BHE_n_unknowns;
+ _R_matrix.setZero(NumRMatrixRows, NumRMatrixRows);
+ _R_pi_s_matrix.setZero(NumRMatrixRows, ShapeFunction::NPOINTS);
+ _R_s_matrix.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+ // formulate the local BHE R matrix
+ Eigen::MatrixXd matBHE_loc_R =
+ Eigen::MatrixXd::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+ for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < BHE_n_unknowns;
+ idx_bhe_unknowns++)
+ {
+ matBHE_loc_R.setZero();
+ // Loop over Gauss points
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+
+ auto const& N = _ip_data[ip].N;
+ auto const& w = _ip_data[ip].integration_weight;
+
+ // get coefficient of R matrix for corresponding BHE.
+ auto R_coeff = _process_data._vec_BHE_property[BHE_id]
+ ->getBoundaryHeatExchangeCoeff(idx_bhe_unknowns);
+
+ // calculate mass matrix for current unknown
+ matBHE_loc_R += N.transpose() * R_coeff * N * w;
+ } // end of loop over integration point
+
+ // The following assembly action is according to Diersch (2013) FEFLOW
+ // book please refer to M.127 and M.128 on page 955 and 956
+ _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
+ // std::string sep = "\n----------------------------------------\n";
+ // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
+ // std::cout << "_R_matrix: \n" << sep;
+ // std::cout << _R_matrix.format(CleanFmt) << sep;
+ // std::cout << "_R_s_matrix: \n" << sep;
+ // std::cout << _R_s_matrix.format(CleanFmt) << sep;
+ // std::cout << "_R_pi_s_matrix: \n" << sep;
+ // std::cout << _R_pi_s_matrix.format(CleanFmt) << sep;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
+void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
+ BHE_Dim>::
+ assemble(
+ double const /*t*/, std::vector<double> const& local_x,
+ std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+ std::vector<double>& /*local_b_data*/) // local b vector is not touched
+{
+ auto const local_matrix_size = local_x.size();
+ const int BHE_n_unknowns = _ip_data[0]._bhe_instance.getNumUnknowns();
+ // plus one because the soil temperature is included in local_x
+ assert(local_matrix_size == ShapeFunction::NPOINTS * (BHE_n_unknowns + 1));
+
+ auto local_M = MathLib::createZeroedMatrix<BheLocalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<BheLocalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition x_position;
+ x_position.setElementID(element_id);
+
+ int shift = 0;
+ static const int local_BHE_matrix_size =
+ ShapeFunction::NPOINTS * BHE_n_unknowns;
+ static const int shift_start = ShapeFunction::NPOINTS;
+
+ // the mass and conductance matrix terms
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto& ip_data = _ip_data[ip];
+
+ auto const& w = ip_data.integration_weight;
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+
+ // looping over all unknowns.
+ for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < BHE_n_unknowns;
+ idx_bhe_unknowns++)
+ {
+ // get coefficient of mass from corresponding BHE.
+ auto& mass_coeff = ip_data._vec_mass_coefficients[idx_bhe_unknowns];
+ auto& laplace_mat = ip_data._vec_mat_Laplace[idx_bhe_unknowns];
+ auto& advection_vec =
+ ip_data._vec_Advection_vectors[idx_bhe_unknowns];
+
+ // calculate shift.
+ shift = shift_start + ShapeFunction::NPOINTS * idx_bhe_unknowns;
+ // local M
+ local_M
+ .template block<ShapeFunction::NPOINTS, ShapeFunction::NPOINTS>(
+ shift, shift)
+ .noalias() += N.transpose() * mass_coeff * N * w;
+
+ // local K
+ // laplace part
+ local_K
+ .block(shift, shift, ShapeFunction::NPOINTS,
+ ShapeFunction::NPOINTS)
+ .noalias() += dNdx.transpose() * laplace_mat * dNdx * w;
+ // advection part
+ local_K
+ .block(shift, shift, ShapeFunction::NPOINTS,
+ ShapeFunction::NPOINTS)
+ .noalias() +=
+ N.transpose() * advection_vec.transpose() * dNdx * w;
+ }
+ }
+
+ // add the R matrix to local_K
+ local_K.block(shift_start, shift_start, local_BHE_matrix_size,
+ local_BHE_matrix_size) += _R_matrix;
+
+ // add the R_pi_s matrix to local_K
+ local_K.block(shift_start, 0, local_BHE_matrix_size, shift_start) +=
+ _R_pi_s_matrix;
+ local_K.block(0, shift_start, shift_start, local_BHE_matrix_size) +=
+ _R_pi_s_matrix.transpose();
+
+ // add the R_s matrix to local_K
+ local_K.block(0, 0, shift_start, shift_start) += 2.0 * _R_s_matrix;
+
+ // debugging
+ // std::string sep =
+ // "\n----------------------------------------\n";
+ // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
+ // std::cout << local_K.format(CleanFmt) << sep;
+ // std::cout << local_M.format(CleanFmt) << sep;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
+void HeatTransportBHELocalAssemblerBHE<
+ ShapeFunction, IntegrationMethod,
+ BHE_Dim>::postTimestepConcrete(std::vector<double> const& /*local_x*/)
+{
+}
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
new file mode 100644
index 0000000000000000000000000000000000000000..3ee83392e2c0674c856465e1fdbf1e399875bd22
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
@@ -0,0 +1,113 @@
+/**
+ * \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
+
+#include <vector>
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h"
+
+#include "HeatTransportBHEProcessAssemblerInterface.h"
+#include "IntegrationPointDataSoil.h"
+#include "SecondaryData.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+const unsigned NUM_NODAL_DOF_SOIL = 1;
+
+template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
+class HeatTransportBHELocalAssemblerSoil
+ : public HeatTransportBHELocalAssemblerInterface
+{
+public:
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ HeatTransportBHELocalAssemblerSoil(
+ HeatTransportBHELocalAssemblerSoil const&) = delete;
+ HeatTransportBHELocalAssemblerSoil(HeatTransportBHELocalAssemblerSoil&&) =
+ delete;
+
+ 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,
+ HeatTransportBHEProcessData& process_data);
+
+ void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
+ std::vector<double>& /*local_M_data*/,
+ 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
+ {
+ auto const& N = _secondary_data.N[integration_point];
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+private:
+ HeatTransportBHEProcessData& _process_data;
+
+ std::vector<
+ IntegrationPointDataSoil<ShapeMatricesType>,
+ Eigen::aligned_allocator<IntegrationPointDataSoil<ShapeMatricesType>>>
+ _ip_data;
+
+ IntegrationMethod const _integration_method;
+
+ std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
+ _shape_matrices;
+
+ std::size_t const element_id;
+
+ bool const _is_axially_symmetric;
+
+ SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
+
+#include "HeatTransportBHELocalAssemblerSoil_impl.h"
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..3cf4209eefe010d583da12f566ae84e8f0c5f8a5
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
@@ -0,0 +1,137 @@
+/**
+ * \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
+
+#include "HeatTransportBHELocalAssemblerSoil.h"
+
+#include <valarray>
+#include <vector>
+
+#include <Eigen/Eigen>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h"
+
+// #include "IntegrationPointDataMatrix.h"
+#include "HeatTransportBHEProcessAssemblerInterface.h"
+#include "SecondaryData.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
+HeatTransportBHELocalAssemblerSoil<ShapeFunction,
+ IntegrationMethod,
+ 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,
+ HeatTransportBHEProcessData& process_data)
+ : HeatTransportBHELocalAssemblerInterface(
+ ShapeFunction::NPOINTS * GlobalDim, dofIndex_to_localIndex),
+ _process_data(process_data),
+ _integration_method(integration_order),
+ element_id(e.getID()),
+ _is_axially_symmetric(is_axially_symmetric)
+{
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N.resize(n_integration_points);
+
+ _shape_matrices = initShapeMatrices<ShapeFunction,
+ ShapeMatricesType,
+ IntegrationMethod,
+ GlobalDim>(
+ e, is_axially_symmetric, _integration_method);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
+void HeatTransportBHELocalAssemblerSoil<
+ ShapeFunction,
+ IntegrationMethod,
+ GlobalDim>::assemble(double const t,
+ std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ 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);
+
+ auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(element_id);
+
+ 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 const k_f = _process_data.thermal_conductivity_fluid(t, pos)[0];
+ // auto const k_g = _process_data.thermal_conductivity_gas(t, pos)[0];
+ auto const k_s = _process_data.thermal_conductivity_solid(t, pos)[0];
+
+ // auto const heat_capacity_f = _process_data.heat_capacity_fluid(t,
+ // pos)[0]; auto const heat_capacity_g =
+ // _process_data.heat_capacity_gas(t, pos)[0];
+ auto const heat_capacity_s =
+ _process_data.heat_capacity_solid(t, pos)[0];
+
+ // auto const density_f = _process_data.density_fluid(t, pos)[0];
+ // auto const density_g = _process_data.density_gas(t, pos)[0];
+ auto const density_s = _process_data.density_solid(t, pos)[0];
+
+ // 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;
+
+ // assemble Mass matrix
+ local_M.noalias() += sm.N.transpose() * density_s * heat_capacity_s *
+ sm.N * sm.detJ * wp.getWeight() *
+ sm.integralMeasure;
+ }
+
+ // debugging
+ // std::string sep = "\n----------------------------------------\n";
+ // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
+ // 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
new file mode 100644
index 0000000000000000000000000000000000000000..2af5e87a70b990a0cb760e2a32c51aa3055a188b
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h
@@ -0,0 +1,64 @@
+/**
+ * \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
+
+#include <utility>
+#include <vector>
+
+#include "BaseLib/Error.h"
+
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+class HeatTransportBHELocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+public:
+ HeatTransportBHELocalAssemblerInterface() : _dofIndex_to_localIndex{} {}
+ 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:
+ std::vector<unsigned> const _dofIndex_to_localIndex;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h
new file mode 100644
index 0000000000000000000000000000000000000000..3f09811c48612ad01e62dd6ce2095f579f723931
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h
@@ -0,0 +1,74 @@
+/**
+ * \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
+
+#include <Eigen/Eigen>
+
+#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+using namespace BHE;
+
+template <typename ShapeMatrixType>
+struct IntegrationPointDataBHE final
+{
+ explicit IntegrationPointDataBHE(BHEAbstract& bhe_instance)
+ : _bhe_instance(bhe_instance)
+ {
+ // depending on the type of BHE
+ const int unknown_size = _bhe_instance.getNumUnknowns();
+ // initialization
+ _vec_mass_coefficients.resize(unknown_size);
+ for (int i = 0; i < unknown_size; i++)
+ {
+ Eigen::MatrixXd mat_laplace(3, 3);
+ mat_laplace.setZero();
+ _vec_mat_Laplace.push_back(mat_laplace);
+ Eigen::VectorXd vec_advection(3);
+ vec_advection.setZero();
+ _vec_Advection_vectors.push_back(vec_advection);
+ }
+
+ // set parameter values
+ for (int j = 0; j < unknown_size; j++)
+ {
+ // mass matrix coefficients
+ _vec_mass_coefficients[j] = _bhe_instance.getMassCoeff(j);
+ // laplace matrix
+ _bhe_instance.getLaplaceMatrix(j, _vec_mat_Laplace[j]);
+ // advection vector
+ _bhe_instance.getAdvectionVector(j, _vec_Advection_vectors[j]);
+ }
+ }
+
+ BHEAbstract& _bhe_instance;
+
+ typename ShapeMatrixType::NodalRowVectorType N;
+ typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
+ double integration_weight;
+
+ // mass coefficients, length depending on the type of BHE
+ std::vector<double> _vec_mass_coefficients;
+
+ // Laplace matrices
+ std::vector<Eigen::MatrixXd> _vec_mat_Laplace;
+
+ // Advection vectors
+ std::vector<Eigen::VectorXd> _vec_Advection_vectors;
+
+ void pushBackState() {}
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
new file mode 100644
index 0000000000000000000000000000000000000000..97abc53d0dc8d7a9091c545c44147f8aaa9db837
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
@@ -0,0 +1,35 @@
+/**
+ * \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
+
+#include <memory>
+#include <vector>
+
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+template <typename ShapeMatrixType>
+struct IntegrationPointDataSoil final
+{
+ explicit IntegrationPointDataSoil() {}
+
+ typename ShapeMatrixType::NodalRowVectorType N;
+ typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
+ double integration_weight;
+
+ void pushBackState() {}
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
new file mode 100644
index 0000000000000000000000000000000000000000..f395d3551890211b73d693412b5a2fb014f50709
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
@@ -0,0 +1,378 @@
+/**
+ * \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
+
+#include <functional>
+#include <memory>
+#include <type_traits>
+#include <typeindex>
+#include <typeinfo>
+#include <unordered_map>
+
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+
+#ifndef OGS_MAX_ELEMENT_DIM
+static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
+#endif
+
+#ifndef OGS_MAX_ELEMENT_ORDER
+static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
+#endif
+
+// The following macros decide which element types will be compiled, i.e.
+// which element types will be available for use in simulations.
+
+#ifdef OGS_ENABLE_ELEMENT_SIMPLEX
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 1u
+#else
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_CUBOID
+#define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1
+#else
+#define ENABLED_ELEMENT_TYPE_CUBOID 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PRISM
+#define ENABLED_ELEMENT_TYPE_PRISM 1u << 2
+#else
+#define ENABLED_ELEMENT_TYPE_PRISM 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PYRAMID
+#define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3
+#else
+#define ENABLED_ELEMENT_TYPE_PYRAMID 0u
+#endif
+
+// Dependent element types.
+// Faces of tets, pyramids and prisms are triangles
+#define ENABLED_ELEMENT_TYPE_TRI \
+ ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+ (ENABLED_ELEMENT_TYPE_PRISM))
+// Faces of hexes, pyramids and prisms are quads
+#define ENABLED_ELEMENT_TYPE_QUAD \
+ ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+ (ENABLED_ELEMENT_TYPE_PRISM))
+
+// All enabled element types
+#define OGS_ENABLED_ELEMENTS \
+ ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
+ (ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
+
+// Include only what is needed (Well, the conditions are not sharp).
+#if OGS_ENABLED_ELEMENTS != 0
+#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
+#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
+#endif
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+/// The LocalDataInitializer is a functor creating a local assembler data with
+/// corresponding to the mesh element type shape functions and calling
+/// initialization of the new local assembler data.
+/// For example for MeshLib::Quad a local assembler data with template argument
+/// NumLib::ShapeQuad4 is created.
+template <typename LocalAssemblerInterface,
+ template <typename, typename, int> class LocalAssemblerDataSoil,
+ template <typename, typename, int> class LocalAssemblerDataBHE,
+ int GlobalDim, typename... ConstructorArgs>
+class LocalDataInitializer final
+{
+public:
+ using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+
+ explicit LocalDataInitializer(
+ NumLib::LocalToGlobalIndexMap const& dof_table)
+ : _dof_table(dof_table)
+ {
+ // REMARKS: At the moment, only a 3D mesh (soil) with 1D elements (BHE)
+ // are supported.
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Quad))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad4>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Hex))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeHex8>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Quad8))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
+ _builder[std::type_index(typeid(MeshLib::Quad9))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Hex20))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
+#endif
+
+ // /// Simplices ////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Tri))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTri3>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Tet))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTet4>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tri6))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Tet10))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
+#endif
+
+ // /// Prisms ////////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Prism))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePrism6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Prism15))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
+#endif
+
+ // /// Pyramids //////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Pyramid))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePyra5>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+ OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
+ makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
+#endif
+ // /// Lines ///////////////////////////////////
+
+#if OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
+ _builder[std::type_index(typeid(MeshLib::Line))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeLine2>();
+#endif
+
+#if OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+ _builder[std::type_index(typeid(MeshLib::Line3))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
+#endif
+ }
+
+ /// Sets the provided \c data_ptr to the newly created local assembler data.
+ ///
+ /// \attention
+ /// The index \c id is not necessarily the mesh item's id. Especially when
+ /// having multiple meshes it will differ from the latter.
+ void operator()(
+ std::size_t const id,
+ MeshLib::Element const& mesh_item,
+ LADataIntfPtr& data_ptr,
+ const std::vector<
+ std::vector<std::size_t>>& /*vec_ele_connected_BHE_IDs*/,
+ const std::vector<
+ std::unique_ptr<BHE::BHEAbstract>>& /*vec_BHE_property*/,
+ ConstructorArgs&&... args) const
+ {
+ auto const type_idx = std::type_index(typeid(mesh_item));
+ auto const it = _builder.find(type_idx);
+
+ if (it == _builder.end())
+ {
+ OGS_FATAL(
+ "You are trying to build a local assembler for an unknown mesh "
+ "element type (%s)."
+ " Maybe you have disabled this mesh element type in your build "
+ "configuration or this process requires higher order elements.",
+ 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)
+ {
+ // this is a BHE element
+ dofIndex_to_localIndex.resize(n_local_dof);
+ unsigned dof_id = 0;
+ unsigned local_id = 0;
+ for (auto i : varIDs)
+ {
+ auto const n_global_components =
+ _dof_table.getNumberOfElementComponents(i);
+ for (unsigned j = 0; j < n_global_components; j++)
+ {
+ auto const& ms = _dof_table.getMeshSubset(i, j);
+ auto const mesh_id = ms.getMeshID();
+ for (unsigned k = 0; k < mesh_item.getNumberOfNodes(); k++)
+ {
+ MeshLib::Location l(mesh_id,
+ MeshLib::MeshItemType::Node,
+ mesh_item.getNodeIndex(k));
+ auto global_index = _dof_table.getGlobalIndex(l, i, j);
+ if (global_index != NumLib::MeshComponentMap::nop)
+ {
+ dofIndex_to_localIndex[dof_id++] = local_id;
+ }
+ local_id++;
+ }
+ }
+ }
+ }
+
+ data_ptr = it->second(mesh_item, varIDs.size(), n_local_dof,
+ 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&&...)>;
+
+ template <typename ShapeFunction>
+ using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
+ typename ShapeFunction::MeshElement>::IntegrationMethod;
+
+ // local assembler builder implementations.
+ template <typename ShapeFunction>
+ using LADataSoil =
+ LocalAssemblerDataSoil<ShapeFunction, IntegrationMethod<ShapeFunction>,
+ GlobalDim>;
+
+ template <typename ShapeFunction>
+ using LADataBHE =
+ LocalAssemblerDataBHE<ShapeFunction, IntegrationMethod<ShapeFunction>,
+ GlobalDim>;
+ /// A helper forwarding to the correct version of makeLocalAssemblerBuilder
+ /// depending whether the global dimension is less than the shape function's
+ /// dimension or not.
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder()
+ {
+ return makeLocalAssemblerBuilder<ShapeFunction>(
+ static_cast<std::integral_constant<
+ bool, (GlobalDim >= ShapeFunction::DIM)>*>(nullptr));
+ }
+
+ /// Mapping of element types to local assembler constructors.
+ std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+ NumLib::LocalToGlobalIndexMap const& _dof_table;
+
+private:
+ /// Generates a function that creates a new LocalAssembler of type
+ /// LAData<ShapeFunction>. Only functions with shape function's dimension
+ /// less or equal to the global dimension are instantiated, e.g. following
+ /// combinations of shape functions and global dimensions: (Line2, 1),
+ /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
+ template <typename ShapeFunction>
+ 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) {
+ if (e.getDimension() == GlobalDim) // soil elements
+ {
+ return LADataIntfPtr{new LADataSoil<ShapeFunction>{
+ e, local_matrix_size, dofIndex_to_localIndex,
+ std::forward<ConstructorArgs>(args)...}};
+ }
+
+ return LADataIntfPtr{new LADataBHE<ShapeFunction>{
+ // BHE elements
+ e, local_matrix_size, dofIndex_to_localIndex,
+ std::forward<ConstructorArgs>(args)...}};
+ };
+ }
+
+ /// Returns nullptr for shape functions whose dimensions are less than the
+ /// global dimension.
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder(std::false_type*)
+ {
+ return nullptr;
+ }
+}; // namespace HeatTransportBHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
+
+#undef ENABLED_ELEMENT_TYPE_SIMPLEX
+#undef ENABLED_ELEMENT_TYPE_CUBOID
+#undef ENABLED_ELEMENT_TYPE_PYRAMID
+#undef ENABLED_ELEMENT_TYPE_PRISM
+#undef ENABLED_ELEMENT_TYPE_TRI
+#undef ENABLED_ELEMENT_TYPE_QUAD
+#undef OGS_ENABLED_ELEMENTS
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/SecondaryData.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/SecondaryData.h
new file mode 100644
index 0000000000000000000000000000000000000000..4a7484928e3e4ead1059862392fe1be081d3d89f
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/SecondaryData.h
@@ -0,0 +1,30 @@
+/**
+ * \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
+
+#include <vector>
+
+#include <Eigen/Eigen>
+
+#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+/// Used by for extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+ std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
+};
+} // namespace HeatTransportBHE
+} // namespace ProcessLib