41 namespace TwoPhaseFlowWithPP
45 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
46 void TwoPhaseFlowWithPPLocalAssembler<
47 ShapeFunction, IntegrationMethod,
48 GlobalDim>::assemble(
double const t,
double const dt,
49 std::vector<double>
const& local_x,
50 std::vector<double>
const& ,
51 std::vector<double>& local_M_data,
52 std::vector<double>& local_K_data,
53 std::vector<double>& local_b_data)
55 auto const local_matrix_size = local_x.size();
57 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
59 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
60 local_M_data, local_matrix_size, local_matrix_size);
61 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
62 local_K_data, local_matrix_size, local_matrix_size);
63 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
64 local_b_data, local_matrix_size);
67 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
68 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
69 auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
70 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
72 auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
73 cap_pressure_matrix_index, cap_pressure_matrix_index);
76 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
79 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
80 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
82 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
83 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
85 auto Klpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
86 cap_pressure_matrix_index, cap_pressure_matrix_index);
88 auto Bg = local_b.template segment<nonwet_pressure_size>(
89 nonwet_pressure_matrix_index);
92 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
94 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
95 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
96 auto const& gas_phase = medium.phase(
"Gas");
98 unsigned const n_integration_points =
99 _integration_method.getNumberOfPoints();
104 for (
unsigned ip = 0; ip < n_integration_points; ip++)
107 auto const& w = _ip_data[ip].integration_weight;
108 auto const& dNdx = _ip_data[ip].dNdx;
109 auto const& M = _ip_data[ip].massOperator;
111 double pc_int_pt = 0.;
112 double pn_int_pt = 0.;
116 _pressure_wet[ip] = pn_int_pt - pc_int_pt;
119 variables[
static_cast<int>(MPL::Variable::phase_pressure)] = pn_int_pt;
123 _process_data.temperature(t, pos)[0];
126 .template value<double>(variables, pos, t, dt);
128 auto const rho_nonwet =
130 .template value<double>(variables, pos, t, dt);
133 .template value<double>(variables, pos, t, dt);
134 auto& Sw = _saturation[ip];
136 .template value<double>(variables, pos, t, dt);
140 .template dValue<double>(
145 .template value<double>(variables, pos, t, dt);
147 auto const drhononwet_dpn =
149 .template dValue<double>(
150 variables, MPL::Variable::phase_pressure, pos, t, dt);
152 auto const k_rel_wet =
154 .template value<double>(variables, pos, t, dt);
155 auto const k_rel_nonwet =
159 .template value<double>(variables, pos, t, dt);
161 auto const mu_nonwet =
163 .template value<double>(variables, pos, t, dt);
165 auto const lambda_nonwet = k_rel_nonwet / mu_nonwet;
168 .template value<double>(variables, pos, t, dt);
170 auto const lambda_wet = k_rel_wet / mu_wet;
172 auto const K = MPL::formEigenTensor<GlobalDim>(
174 .template value<double>(variables, pos, t, dt));
176 Mgp.noalias() +=
porosity * (1 - Sw) * drhononwet_dpn * M;
177 Mgpc.noalias() += -
porosity * rho_nonwet * dSw_dpc * M;
178 Mlpc.noalias() +=
porosity * dSw_dpc * rho_wet * M;
180 laplace_operator.noalias() = dNdx.transpose() * K * dNdx * w;
182 Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator;
183 Klp.noalias() += rho_wet * lambda_wet * laplace_operator;
184 Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator;
186 if (_process_data.has_gravity)
188 auto const& b = _process_data.specific_body_force;
192 rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator;
193 Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator;
196 if (_process_data.has_mass_lumping)
198 for (
unsigned row = 0; row < Mgpc.cols(); row++)
200 for (
unsigned column = 0; column < Mgpc.cols(); column++)
204 Mgpc(row, row) += Mgpc(row, column);
205 Mgpc(row, column) = 0.0;
206 Mgp(row, row) += Mgp(row, column);
207 Mgp(row, column) = 0.0;
208 Mlpc(row, row) += Mlpc(row, column);
209 Mlpc(row, column) = 0.0;
Definition of the PiecewiseLinearInterpolation class.
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
typename ShapeMatricesType::NodalVectorType NodalVectorType
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
@ relative_permeability_nonwetting_phase
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
const unsigned NUM_NODAL_DOF