24 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
27 Eigen::VectorXd
const& local_x,
28 Eigen::VectorXd
const& local_xdot,
30 std::vector<double>& local_M_data,
31 std::vector<double>& local_K_data,
32 std::vector<double>& local_b_data)
34 if (process_id == this->_process_data.heat_transport_process_id)
36 assembleHeatTransportEquation(t, dt, local_x, local_M_data,
41 assembleHydraulicEquation(t, dt, local_x, local_xdot, local_M_data,
42 local_K_data, local_b_data);
45 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
48 Eigen::VectorXd
const& local_x,
49 Eigen::VectorXd
const& local_xdot,
50 std::vector<double>& local_M_data,
51 std::vector<double>& local_K_data,
52 std::vector<double>& local_b_data)
55 local_x.template segment<pressure_size>(pressure_index);
57 local_x.template segment<temperature_size>(temperature_index);
59 auto const local_Tdot =
60 local_xdot.template segment<temperature_size>(temperature_index);
62 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
63 local_M_data, pressure_size, pressure_size);
64 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
65 local_K_data, pressure_size, pressure_size);
66 auto local_b = MathLib::createZeroedVector<LocalVectorType>(local_b_data,
72 auto const& process_data = this->_process_data;
73 auto const& medium = *this->_process_data.media_map->getMedium(
74 this->_element.getID());
75 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
76 auto const& solid_phase = medium.phase(
"Solid");
78 auto const& b = process_data.specific_body_force;
82 unsigned const n_integration_points =
83 this->_integration_method.getNumberOfPoints();
85 for (
unsigned ip(0); ip < n_integration_points; ip++)
89 auto const& ip_data = this->_ip_data[ip];
90 auto const& N = ip_data.N;
91 auto const& dNdx = ip_data.dNdx;
92 auto const& w = ip_data.integration_weight;
94 double p_int_pt = 0.0;
95 double T_int_pt = 0.0;
106 auto const porosity =
108 .template value<double>(vars, pos, t, dt);
111 .template value<double>(vars, pos, t, dt);
113 const double dfluid_density_dp =
115 .template dValue<double>(
122 .template value<double>(vars, pos, t, dt);
126 auto const specific_storage =
128 .template value<double>(vars, pos, t, dt);
130 auto const intrinsic_permeability =
131 MaterialPropertyLib::formEigenTensor<GlobalDim>(
133 .value(vars, pos, t, dt));
140 (porosity * dfluid_density_dp /
fluid_density + specific_storage) *
143 local_K.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
145 if (process_data.has_gravity)
151 if (!process_data.has_fluid_thermal_expansion)
158 auto const solid_thermal_expansion =
159 process_data.solid_thermal_expansion(t, pos)[0];
160 const double dfluid_density_dT =
163 .template dValue<double>(
166 double Tdot_int_pt = 0.;
168 auto const biot_constant =
169 process_data.biot_constant(t, pos)[0];
170 const double eff_thermal_expansion =
171 3.0 * (biot_constant - porosity) * solid_thermal_expansion -
173 local_b.noalias() += eff_thermal_expansion * Tdot_int_pt * w * N;
178 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
182 Eigen::VectorXd
const& local_x,
183 std::vector<double>& local_M_data,
184 std::vector<double>& local_K_data)
187 local_x.template segment<pressure_size>(pressure_index);
189 local_x.template segment<temperature_size>(temperature_index);
191 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
192 local_M_data, temperature_size, temperature_size);
193 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
194 local_K_data, temperature_size, temperature_size);
199 auto const& process_data = this->_process_data;
201 *process_data.media_map->getMedium(this->_element.getID());
202 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
204 auto const& b = process_data.specific_body_force;
207 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
211 unsigned const n_integration_points =
212 this->_integration_method.getNumberOfPoints();
214 for (
unsigned ip(0); ip < n_integration_points; ip++)
218 auto const& ip_data = this->_ip_data[ip];
219 auto const& N = ip_data.N;
220 auto const& dNdx = ip_data.dNdx;
221 auto const& w = ip_data.integration_weight;
235 auto const porosity =
237 .template value<double>(vars, pos, t, dt);
244 .template value<double>(vars, pos, t, dt);
245 auto const specific_heat_capacity_fluid =
247 .template value<double>(vars, pos, t, dt);
250 local_M.noalias() += w *
251 this->getHeatEnergyCoefficient(
253 specific_heat_capacity_fluid, pos, t, dt) *
259 .template value<double>(vars, pos, t, dt);
261 auto const intrinsic_permeability =
262 MaterialPropertyLib::formEigenTensor<GlobalDim>(
265 .value(vars, pos, t, dt));
270 process_data.has_gravity
276 this->getThermalConductivityDispersivity(
278 velocity, I, pos, t, dt);
281 w * (dNdx.transpose() * thermal_conductivity_dispersivity * dNdx +
282 N.transpose() * velocity.transpose() * dNdx *
fluid_density *
283 specific_heat_capacity_fluid);
287 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
288 std::vector<double>
const&
292 std::vector<GlobalVector*>
const& x,
293 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
294 std::vector<double>& cache)
const
296 assert(x.size() == dof_table.size());
297 auto const n_processes = dof_table.size();
299 std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes;
300 indices_of_all_coupled_processes.reserve(n_processes);
301 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
305 assert(!indices.empty());
306 indices_of_all_coupled_processes.push_back(indices);
308 auto const local_xs =
311 return this->getIntPtDarcyVelocityLocal(t, local_xs, cache);
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void assembleHeatTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data)
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
double fluid_density(const double p, const double T, const double x)
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType >> const &indices)