75 {
76 auto const local_matrix_size = local_x.size();
77
78
79 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
80
81 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
82 local_M_data, local_matrix_size, local_matrix_size);
83 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
84 local_K_data, local_matrix_size, local_matrix_size);
85 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
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
102
103 auto p_nodal_values = Eigen::Map<const NodalVectorType>(
105
106 auto const& medium =
108 auto const& liquid_phase = medium.phase("AqueousLiquid");
109 auto const& solid_phase = medium.phase("Solid");
110
111 auto const& b =
112 process_data
113 .projected_specific_body_force_vectors[this->
_element.
getID()];
114
116
117 unsigned const n_integration_points =
119
120 std::vector<GlobalDimVectorType> ip_flux_vector;
121 double average_velocity_norm = 0.0;
122 ip_flux_vector.reserve(n_integration_points);
123
124 auto const& Ns =
125 process_data.shape_matrix_cache
126 .template NsHigherOrder<typename ShapeFunction::MeshElement>();
127
128 for (unsigned ip(0); ip < n_integration_points; ip++)
129 {
131
132 auto const& ip_data = this->
_ip_data[ip];
133 auto const& dNdx = ip_data.dNdx;
134 auto const&
N = Ns[ip];
135 auto const& w = ip_data.integration_weight;
136
137 double T_int_pt = 0.0;
138 double p_int_pt = 0.0;
139
141
144
146
147
148 auto const specific_storage =
150 .template value<double>(vars, pos, t, dt);
151
154 .template value<double>(vars, pos, t, dt);
156
157 auto const intrinsic_permeability =
158 MaterialPropertyLib::formEigenTensor<GlobalDim>(
159 medium
160 .property(
162 .value(vars, pos, t, dt));
163
164 auto const specific_heat_capacity_fluid =
165 liquid_phase
167 .template value<double>(vars, pos, t, dt);
168
169
171 liquid_phase
173 .template value<double>(vars, pos, t, dt);
174
176
178 liquid_phase
180 .template value<double>(vars, pos, t, dt);
182
184 process_data.has_gravity
186 fluid_density * b))
188
189
192 vars, fluid_density, specific_heat_capacity_fluid, velocity,
193 pos, t, dt);
194
195 KTT.noalias() +=
196 dNdx.transpose() * thermal_conductivity_dispersivity * dNdx * w;
197
198 ip_flux_vector.emplace_back(velocity * fluid_density *
199 specific_heat_capacity_fluid);
200 average_velocity_norm += velocity.norm();
201
202 Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
203 MTT.noalias() += w *
205 vars, porosity, fluid_density,
206 specific_heat_capacity_fluid, pos, t, dt) *
208 Mpp.noalias() += w *
N.transpose() * specific_storage *
N;
209 if (process_data.has_gravity)
210 {
212 }
213
214
215 }
216
217 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
218 process_data.stabilizer, this->_ip_data,
219 process_data.shape_matrix_cache, ip_flux_vector,
220 average_velocity_norm / static_cast<double>(n_integration_points),
221 KTT);
222 }
double liquid_phase_pressure
std::size_t getID() const
Returns the ID of the element.
unsigned getNumberOfPoints() const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
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
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
double fluid_density(const double p, const double T, const double x)