72 double const t,
double const dt, std::vector<double>
const& local_x,
73 std::vector<double>
const& ,
74 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
75 std::vector<double>& local_b_data)
77 auto const local_matrix_size = local_x.size();
79 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
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);
89 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
90 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
91 auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
92 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
94 auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
95 cap_pressure_matrix_index, cap_pressure_matrix_index);
98 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
101 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
102 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
104 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
105 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
107 auto Klpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
108 cap_pressure_matrix_index, cap_pressure_matrix_index);
110 auto Bg = local_b.template segment<nonwet_pressure_size>(
111 nonwet_pressure_matrix_index);
114 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
116 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
117 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
118 auto const& gas_phase = medium.phase(
"Gas");
120 unsigned const n_integration_points =
121 _integration_method.getNumberOfPoints();
126 for (
unsigned ip = 0; ip < n_integration_points; ip++)
128 auto const& w = _ip_data[ip].integration_weight;
129 auto const& N = _ip_data[ip].N;
130 auto const& dNdx = _ip_data[ip].dNdx;
134 double pc_int_pt = 0.;
135 double pn_int_pt = 0.;
139 _pressure_wet[ip] = pn_int_pt - pc_int_pt;
144 variables.
temperature = _process_data.temperature(t, pos)[0];
146 gas_phase.property(MPL::PropertyType::molar_mass)
147 .template value<double>(variables, pos, t, dt);
149 auto const rho_nonwet =
150 gas_phase.property(MPL::PropertyType::density)
151 .template value<double>(variables, pos, t, dt);
153 auto const rho_wet = liquid_phase.property(MPL::PropertyType::density)
154 .template value<double>(variables, pos, t, dt);
155 auto& Sw = _saturation[ip];
156 Sw = medium.property(MPL::PropertyType::saturation)
157 .template value<double>(variables, pos, t, dt);
161 medium.property(MPL::PropertyType::saturation)
162 .template dValue<double>(
163 variables, MPL::Variable::capillary_pressure, pos, t, dt);
165 auto const porosity =
166 medium.property(MPL::PropertyType::porosity)
167 .template value<double>(variables, pos, t, dt);
169 auto const drhononwet_dpn =
170 gas_phase.property(MPL::PropertyType::density)
171 .template dValue<double>(
172 variables, MPL::Variable::gas_phase_pressure, pos, t, dt);
174 auto const k_rel_wet =
175 medium.property(MPL::PropertyType::relative_permeability)
176 .template value<double>(variables, pos, t, dt);
177 auto const k_rel_nonwet =
180 MPL::PropertyType::relative_permeability_nonwetting_phase)
181 .template value<double>(variables, pos, t, dt);
183 auto const mu_nonwet =
184 gas_phase.property(MPL::PropertyType::viscosity)
185 .template value<double>(variables, pos, t, dt);
187 auto const lambda_nonwet = k_rel_nonwet / mu_nonwet;
189 auto const mu_wet = liquid_phase.property(MPL::PropertyType::viscosity)
190 .template value<double>(variables, pos, t, dt);
192 auto const lambda_wet = k_rel_wet / mu_wet;
195 medium.property(MPL::PropertyType::permeability)
196 .template value<double>(variables, pos, t, dt));
198 Mgp.noalias() += porosity * (1 - Sw) * drhononwet_dpn * M;
199 Mgpc.noalias() += -porosity * rho_nonwet * dSw_dpc * M;
200 Mlpc.noalias() += porosity * dSw_dpc * rho_wet * M;
202 laplace_operator.noalias() = dNdx.transpose() * K * dNdx * w;
204 Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator;
205 Klp.noalias() += rho_wet * lambda_wet * laplace_operator;
206 Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator;
208 if (_process_data.has_gravity)
210 auto const& b = _process_data.specific_body_force;
214 rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator;
215 Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator;
218 if (_process_data.has_mass_lumping)
220 for (
unsigned row = 0; row < Mgpc.cols(); row++)
222 for (
unsigned column = 0; column < Mgpc.cols(); column++)
226 Mgpc(row, row) += Mgpc(row, column);
227 Mgpc(row, column) = 0.0;
228 Mgp(row, row) += Mgp(row, column);
229 Mgp(row, column) = 0.0;
230 Mlpc(row, row) += Mlpc(row, column);
231 Mlpc(row, column) = 0.0;