70 double const t,
double const dt, std::vector<double>
const& local_x,
71 std::vector<double>
const& ,
72 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
73 std::vector<double>& )
75 assert(local_x.size() == ShapeFunction::NPOINTS);
79 local_M_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
81 local_K_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
86 auto const& medium = *_process_data.media_map.getMedium(_element_id);
87 auto const& solid_phase = medium.phase(
"Solid");
88 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
92 unsigned const n_integration_points =
93 _integration_method.getNumberOfPoints();
95 for (
unsigned ip = 0; ip < n_integration_points; ip++)
97 auto& ip_data = _ip_data[ip];
98 auto const& N = ip_data.N;
99 auto const& dNdx = ip_data.dNdx;
100 auto const& w = ip_data.integration_weight;
102 double T_int_pt = 0.0;
108 auto const density_s =
110 .template value<double>(vars, pos, t, dt);
112 auto const heat_capacity_s =
116 .template value<double>(vars, pos, t, dt);
118 auto const density_f =
120 .template value<double>(vars, pos, t, dt);
122 auto const heat_capacity_f =
126 .template value<double>(vars, pos, t, dt);
128 auto const porosity =
130 .template value<double>(vars, pos, t, dt);
132 auto const velocity =
135 .template value<Eigen::Vector3d>(vars, pos, t, dt);
138 auto const thermal_conductivity =
143 .value(vars, pos, t, dt));
145 auto thermal_conductivity_dispersivity = thermal_conductivity;
147 double const velocity_magnitude = velocity.norm();
149 if (velocity_magnitude >= std::numeric_limits<double>::epsilon())
151 auto const thermal_dispersivity_longitudinal =
154 thermal_longitudinal_dispersivity)
155 .template value<double>();
156 auto const thermal_dispersivity_transversal =
159 thermal_transversal_dispersivity)
160 .template value<double>();
162 auto const thermal_dispersivity =
163 density_f * heat_capacity_f *
164 (thermal_dispersivity_transversal * velocity_magnitude *
165 Eigen::Matrix3d::Identity() +
166 (thermal_dispersivity_longitudinal -
167 thermal_dispersivity_transversal) /
168 velocity_magnitude * velocity * velocity.transpose());
169 thermal_conductivity_dispersivity += thermal_dispersivity;
174 (dNdx.transpose() * thermal_conductivity_dispersivity * dNdx +
175 N.transpose() * velocity.transpose() * dNdx * density_f *
180 local_M.noalias() += N.transpose() * N * w *
181 (density_s * heat_capacity_s * (1 - porosity) +
182 density_f * heat_capacity_f * porosity);
185 if (_process_data._mass_lumping)
188 if (_process_data.mass_lumping_soil_elements[_element_id])
190 local_M = local_M.colwise().sum().eval().asDiagonal();
203 double const t,
double const dt, std::vector<double>
const& local_x,
204 std::vector<double>
const& local_x_prev,
205 std::vector<double>& local_rhs_data, std::vector<double>& local_Jac_data)
207 assert(local_x.size() == ShapeFunction::NPOINTS);
208 auto const local_matrix_size = local_x.size();
211 Eigen::Map<NodalVectorType const>(local_x.data(), local_matrix_size);
212 auto x_prev = Eigen::Map<NodalVectorType const>(local_x_prev.data(),
216 local_Jac_data, local_matrix_size, local_matrix_size);
218 local_rhs_data, local_matrix_size);
220 std::vector<double> local_M_data;
221 std::vector<double> local_K_data;
222 assemble(t, dt, local_x, local_x_prev, local_M_data, local_K_data,
227 local_M_data, local_matrix_size, local_matrix_size);
229 local_K_data, local_matrix_size, local_matrix_size);
232 local_Jac.noalias() += local_K + local_M / dt;
233 local_rhs.noalias() -= local_K * x + local_M * (x - x_prev) / dt;