OGS
HeatTransportBHELocalAssemblerBHE-impl.h
Go to the documentation of this file.
1 
11 #pragma once
12 
16 
17 namespace ProcessLib
18 {
19 namespace HeatTransportBHE
20 {
21 template <typename ShapeFunction, typename IntegrationMethod, typename BHEType>
24  BHEType const& bhe,
25  bool const is_axially_symmetric,
26  unsigned const integration_order,
27  HeatTransportBHEProcessData& process_data)
28  : _process_data(process_data),
29  _integration_method(integration_order),
30  _bhe(bhe),
31  _element_id(e.getID())
32 {
33  // need to make sure that the BHE elements are one-dimensional
34  assert(e.getDimension() == 1);
35 
36  unsigned const n_integration_points =
37  _integration_method.getNumberOfPoints();
38 
39  _ip_data.reserve(n_integration_points);
40  _secondary_data.N.resize(n_integration_points);
41 
42  auto const shape_matrices =
44  3 /* GlobalDim */>(e, is_axially_symmetric,
46 
47  // ip data initialization
48  for (unsigned ip = 0; ip < n_integration_points; ip++)
49  {
50  auto const& sm = shape_matrices[ip];
51  // create the class IntegrationPointDataBHE in place
52  _ip_data.push_back(
53  {sm.N, sm.dNdx,
54  _integration_method.getWeightedPoint(ip).getWeight() *
55  sm.integralMeasure * sm.detJ});
56 
57  _secondary_data.N[ip] = sm.N;
58  }
59 
60  // calculate the element direction vector
61  auto const p0 =
62  Eigen::Map<Eigen::Vector3d const>(e.getNode(0)->getCoords(), 3);
63  auto const p1 =
64  Eigen::Map<Eigen::Vector3d const>(e.getNode(1)->getCoords(), 3);
65 
66  _element_direction = (p1 - p0).normalized();
67 
71  static constexpr int max_num_thermal_exchange_terms = 5;
72  // formulate the local BHE R matrix
73  for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < bhe_unknowns;
74  idx_bhe_unknowns++)
75  {
76  typename ShapeMatricesType::template MatrixType<
78  matBHE_loc_R = ShapeMatricesType::template MatrixType<
82  // Loop over Gauss points
83  for (unsigned ip = 0; ip < n_integration_points; ip++)
84  {
85  auto const& N = _ip_data[ip].N;
86  auto const& w = _ip_data[ip].integration_weight;
87 
88  auto const& R = _bhe.thermalResistance(idx_bhe_unknowns);
89  // calculate mass matrix for current unknown
90  matBHE_loc_R += N.transpose() * N * (1 / R) * w;
91  } // end of loop over integration point
92 
93  // The following assembly action is according to Diersch (2013) FEFLOW
94  // book please refer to M.127 and M.128 on page 955 and 956
95  // The if check is absolutely necessary because
96  // (i) In the CXA and CXC case, there are 3 exchange terms,
97  // and it is the same as the number of unknowns;
98  // (ii) In the 1U case, there are 4 exchange terms,
99  // and it is again same as the number of unknowns;
100  // (iii) In the 2U case, there are 5 exchange terms,
101  // and it is less than the number of unknowns (8).
102  if (idx_bhe_unknowns < max_num_thermal_exchange_terms)
103  {
104  _bhe.template assembleRMatrices<ShapeFunction::NPOINTS>(
105  idx_bhe_unknowns, matBHE_loc_R, _R_matrix, _R_pi_s_matrix,
106  _R_s_matrix);
107  }
108  } // end of loop over BHE unknowns
109 
110  // debugging
111  // std::string sep =
112  // "\n----------------------------------------\n";
113  // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
114  // std::cout << "_R_matrix: \n" << sep;
115  // std::cout << _R_matrix.format(CleanFmt) << sep;
116  // std::cout << "_R_s_matrix: \n" << sep;
117  // std::cout << _R_s_matrix.format(CleanFmt) << sep;
118  // std::cout << "_R_pi_s_matrix: \n" << sep;
119  // std::cout << _R_pi_s_matrix.format(CleanFmt) << sep;
120 }
121 
122 template <typename ShapeFunction, typename IntegrationMethod, typename BHEType>
123 void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
124  BHEType>::
125  assemble(
126  double const /*t*/, double const /*dt*/,
127  std::vector<double> const& /*local_x*/,
128  std::vector<double> const& /*local_xdot*/,
129  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
130  std::vector<double>& /*local_b_data*/) // local b vector is not touched
131 {
132  auto local_M = MathLib::createZeroedMatrix<BheLocalMatrixType>(
133  local_M_data, local_matrix_size, local_matrix_size);
134  auto local_K = MathLib::createZeroedMatrix<BheLocalMatrixType>(
135  local_K_data, local_matrix_size, local_matrix_size);
136 
137  unsigned const n_integration_points =
138  _integration_method.getNumberOfPoints();
139 
140  auto const& pipe_heat_capacities = _bhe.pipeHeatCapacities();
141  auto const& pipe_heat_conductions = _bhe.pipeHeatConductions();
142  auto const& pipe_advection_vectors =
143  _bhe.pipeAdvectionVectors(_element_direction);
144  auto const& cross_section_areas = _bhe.crossSectionAreas();
145 
146  // the mass and conductance matrix terms
147  for (unsigned ip = 0; ip < n_integration_points; ip++)
148  {
149  auto& ip_data = _ip_data[ip];
150 
151  auto const& w = ip_data.integration_weight;
152  auto const& N = ip_data.N;
153  auto const& dNdx = ip_data.dNdx;
154 
155  // looping over all unknowns.
156  for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < bhe_unknowns;
157  idx_bhe_unknowns++)
158  {
159  // get coefficient of mass from corresponding BHE.
160  auto const& mass_coeff = pipe_heat_capacities[idx_bhe_unknowns];
161  auto const& lambda = pipe_heat_conductions[idx_bhe_unknowns];
162  auto const& advection_vector =
163  pipe_advection_vectors[idx_bhe_unknowns];
164  auto const& A = cross_section_areas[idx_bhe_unknowns];
165 
166  int const single_bhe_unknowns_index =
167  bhe_unknowns_index +
168  single_bhe_unknowns_size * idx_bhe_unknowns;
169  // local M
170  local_M
171  .template block<single_bhe_unknowns_size,
172  single_bhe_unknowns_size>(
173  single_bhe_unknowns_index, single_bhe_unknowns_index)
174  .noalias() += N.transpose() * N * mass_coeff * A * w;
175 
176  // local K
177  // laplace part
178  local_K
179  .template block<single_bhe_unknowns_size,
180  single_bhe_unknowns_size>(
181  single_bhe_unknowns_index, single_bhe_unknowns_index)
182  .noalias() += dNdx.transpose() * dNdx * lambda * A * w;
183  // advection part
184  local_K
185  .template block<single_bhe_unknowns_size,
186  single_bhe_unknowns_size>(
187  single_bhe_unknowns_index, single_bhe_unknowns_index)
188  .noalias() +=
189  N.transpose() * advection_vector.transpose() * dNdx * A * w;
190  }
191  }
192 
193  // add the R matrix to local_K
194  local_K.template block<bhe_unknowns_size, bhe_unknowns_size>(
195  bhe_unknowns_index, bhe_unknowns_index) += _R_matrix;
196 
197  // add the R_pi_s matrix to local_K
198  local_K
199  .template block<bhe_unknowns_size, soil_temperature_size>(
200  bhe_unknowns_index, soil_temperature_index)
201  .noalias() += _R_pi_s_matrix;
202  local_K
203  .template block<soil_temperature_size, bhe_unknowns_size>(
204  soil_temperature_index, bhe_unknowns_index)
205  .noalias() += _R_pi_s_matrix.transpose();
206 
207  // add the R_s matrix to local_K
208  local_K
209  .template block<soil_temperature_size, soil_temperature_size>(
210  soil_temperature_index, soil_temperature_index)
211  .noalias() += _bhe.number_of_grout_zones * _R_s_matrix;
212 
213  // debugging
214  // std::string sep = "\n----------------------------------------\n";
215  // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
216  // std::cout << local_K.format(CleanFmt) << sep;
217  // std::cout << local_M.format(CleanFmt) << sep;
218 }
219 } // namespace HeatTransportBHE
220 } // namespace ProcessLib
const T * getCoords() const
Definition: TemplatePoint.h:75
virtual const Node * getNode(unsigned idx) const =0
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
ShapeMatricesType::template MatrixType< soil_temperature_size, soil_temperature_size > _R_s_matrix
ShapeMatricesType::template MatrixType< bhe_unknowns_size, soil_temperature_size > _R_pi_s_matrix
std::vector< IntegrationPointDataBHE< ShapeMatricesType >, Eigen::aligned_allocator< IntegrationPointDataBHE< ShapeMatricesType > > > _ip_data
ShapeMatricesType::template MatrixType< bhe_unknowns_size, bhe_unknowns_size > _R_matrix
HeatTransportBHELocalAssemblerBHE(HeatTransportBHELocalAssemblerBHE const &)=delete
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:28