75 {
76 auto const local_matrix_size = local_x.size();
77
78
79 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
80
82 local_M_data, local_matrix_size, local_matrix_size);
84 local_K_data, local_matrix_size, local_matrix_size);
86 local_b_data, local_matrix_size);
87
88 auto KTT = local_K.template block<temperature_size, temperature_size>(
90 auto MTT = local_M.template block<temperature_size, temperature_size>(
92 auto Kpp = local_K.template block<pressure_size, pressure_size>(
94 auto Mpp = local_M.template block<pressure_size, pressure_size>(
96 auto Bp = local_b.template block<pressure_size, 1>(
pressure_index, 0);
97
99
100 auto p_nodal_values = Eigen::Map<const NodalVectorType>(
102
103 auto const& medium =
105 auto const& liquid_phase = medium.phase("AqueousLiquid");
106 auto const& solid_phase = medium.phase("Solid");
107
108 auto const& b =
109 process_data
110 .projected_specific_body_force_vectors[this->
_element.
getID()];
111
113
114 unsigned const n_integration_points =
116
117 std::vector<GlobalDimVectorType> ip_flux_vector;
118 double average_velocity_norm = 0.0;
119 ip_flux_vector.reserve(n_integration_points);
120
121 auto const& Ns =
122 process_data.shape_matrix_cache
123 .template NsHigherOrder<typename ShapeFunction::MeshElement>();
124
125 for (unsigned ip(0); ip < n_integration_points; ip++)
126 {
127 auto const& ip_data = this->
_ip_data[ip];
128 auto const& dNdx = ip_data.dNdx;
129 auto const&
N = Ns[ip];
130 auto const& w = ip_data.integration_weight;
131
138
139 double T_int_pt = 0.0;
140 double p_int_pt = 0.0;
141
143
146
148
149
150 auto const specific_storage =
152 .template value<double>(vars, pos, t, dt);
153
156 .template value<double>(vars, pos, t, dt);
158
159 auto const intrinsic_permeability =
161 medium
162 .property(
164 .value(vars, pos, t, dt));
165
166 auto const specific_heat_capacity_fluid =
167 liquid_phase
169 .template value<double>(vars, pos, t, dt);
170
171
172 auto const fluid_density =
173 liquid_phase
175 .template value<double>(vars, pos, t, dt);
176
178
180 liquid_phase
182 .template value<double>(vars, pos, t, dt);
184
186 process_data.has_gravity
188 fluid_density * b))
190
191
194 vars, fluid_density, specific_heat_capacity_fluid, velocity,
195 pos, t, dt);
196
197 KTT.noalias() +=
198 dNdx.transpose() * thermal_conductivity_dispersivity * dNdx * w;
199
200 ip_flux_vector.emplace_back(velocity * fluid_density *
201 specific_heat_capacity_fluid);
202 average_velocity_norm += velocity.norm();
203
204 Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
205 MTT.noalias() += w *
207 vars, porosity, fluid_density,
208 specific_heat_capacity_fluid, pos, t, dt) *
210 Mpp.noalias() += w *
N.transpose() * specific_storage *
N;
211 if (process_data.has_gravity)
212 {
213 Bp += w * fluid_density * dNdx.transpose() * K_over_mu * b;
214 }
215
216
217 }
218
220 process_data.stabilizer, this->_ip_data,
221 process_data.shape_matrix_cache, ip_flux_vector,
222 average_velocity_norm / static_cast<double>(n_integration_points),
223 KTT);
224 }
double liquid_phase_pressure
std::size_t getID() const
Returns the ID of the element.
unsigned getNumberOfPoints() const
NumLib::GenericIntegrationMethod const & _integration_method
static const int temperature_index
double getHeatEnergyCoefficient(MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
static const int pressure_size
GlobalDimMatrixType getThermalConductivityDispersivity(MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
MeshLib::Element const & _element
HTProcessData const & _process_data
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
static const int pressure_index
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
void assembleAdvectionMatrix(IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)