51{
52 auto const local_matrix_size = local_x.size();
53
54 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
55
57 local_M_data, local_matrix_size, local_matrix_size);
59 local_K_data, local_matrix_size, local_matrix_size);
61 local_b_data, local_matrix_size);
62
63 auto Mgp =
64 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
66 auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
68
69 auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
71
73 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
74
75 auto Kgp =
76 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
78
79 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
81
82 auto Klpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
84
85 auto Bg = local_b.template segment<nonwet_pressure_size>(
87
88 auto Bl =
90
92 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
93 auto const& gas_phase = medium.phase("Gas");
94
95 unsigned const n_integration_points =
97
100
101 for (unsigned ip = 0; ip < n_integration_points; ip++)
102 {
104 auto const& w =
_ip_data[ip].integration_weight;
105 auto const& dNdx =
_ip_data[ip].dNdx;
106 auto const& M =
_ip_data[ip].massOperator;
107
108 double pc_int_pt = 0.;
109 double pn_int_pt = 0.;
111 pc_int_pt);
112
115
120 gas_phase.property(MPL::PropertyType::molar_mass)
121 .template value<double>(variables, pos, t, dt);
122
123 auto const rho_nonwet =
124 gas_phase.property(MPL::PropertyType::density)
125 .template value<double>(variables, pos, t, dt);
126
127 auto const rho_wet = liquid_phase.property(MPL::PropertyType::density)
128 .template value<double>(variables, pos, t, dt);
130 Sw = medium.property(MPL::PropertyType::saturation)
131 .template value<double>(variables, pos, t, dt);
132
133 auto const dSw_dpc =
134 medium.property(MPL::PropertyType::saturation)
135 .template dValue<double>(
136 variables, MPL::Variable::capillary_pressure, pos, t, dt);
137
139 medium.property(MPL::PropertyType::porosity)
140 .template value<double>(variables, pos, t, dt);
141
142 auto const drhononwet_dpn =
143 gas_phase.property(MPL::PropertyType::density)
144 .template dValue<double>(
145 variables, MPL::Variable::gas_phase_pressure, pos, t, dt);
146
147 auto const k_rel_wet =
148 medium.property(MPL::PropertyType::relative_permeability)
149 .template value<double>(variables, pos, t, dt);
150 auto const k_rel_nonwet =
151 medium
152 .property(
153 MPL::PropertyType::relative_permeability_nonwetting_phase)
154 .template value<double>(variables, pos, t, dt);
155
156 auto const mu_nonwet =
157 gas_phase.property(MPL::PropertyType::viscosity)
158 .template value<double>(variables, pos, t, dt);
159
160 auto const lambda_nonwet = k_rel_nonwet / mu_nonwet;
161
162 auto const mu_wet = liquid_phase.property(MPL::PropertyType::viscosity)
163 .template value<double>(variables, pos, t, dt);
164
165 auto const lambda_wet = k_rel_wet / mu_wet;
166
168 medium.property(MPL::PropertyType::permeability)
169 .template value<double>(variables, pos, t, dt));
170
171 Mgp.noalias() +=
porosity * (1 - Sw) * drhononwet_dpn * M;
172 Mgpc.noalias() += -
porosity * rho_nonwet * dSw_dpc * M;
173 Mlpc.noalias() +=
porosity * dSw_dpc * rho_wet * M;
174
175 laplace_operator.noalias() = dNdx.transpose() * K * dNdx * w;
176
177 Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator;
178 Klp.noalias() += rho_wet * lambda_wet * laplace_operator;
179 Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator;
180
182 {
184
186 Bg.noalias() +=
187 rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator;
188 Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator;
189 }
190 }
192 {
193 for (unsigned row = 0; row < Mgpc.cols(); row++)
194 {
195 for (unsigned column = 0; column < Mgpc.cols(); column++)
196 {
197 if (row != column)
198 {
199 Mgpc(row, row) += Mgpc(row, column);
200 Mgpc(row, column) = 0.0;
201 Mgp(row, row) += Mgp(row, column);
202 Mgp(row, column) = 0.0;
203 Mlpc(row, row) += Mlpc(row, column);
204 Mlpc(row, column) = 0.0;
205 }
206 }
207 }
208 }
209}
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
double gas_phase_pressure
double capillary_pressure
std::size_t getID() const
Returns the ID of the element.
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
static const int nonwet_pressure_matrix_index
static const int cap_pressure_matrix_index
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
typename ShapeMatricesType::NodalVectorType NodalVectorType
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 &)
Eigen::VectorXd const specific_body_force
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
bool const has_mass_lumping
Enables lumping of the mass matrix.
ParameterLib::Parameter< double > const & temperature