34{
   37 
   38    auto const local_matrix_size = local_x.size();
   39 
   40    assert(local_matrix_size == ShapeFunction::NPOINTS * 
NUM_NODAL_DOF);
 
   41 
   43        local_M_data, local_matrix_size, local_matrix_size);
   45        local_K_data, local_matrix_size, local_matrix_size);
   47        local_b_data, local_matrix_size);
   48 
   49    auto Map =
   50        local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
   52    auto Mapc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
   54    auto Max = local_M.template block<nonwet_pressure_size, contaminant_size>(
   56    auto Mat = local_M.template block<nonwet_pressure_size, temperature_size>(
   58 
   59    auto Mwp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
   61    auto Mwpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
   63    auto Mwx = local_M.template block<cap_pressure_size, contaminant_size>(
   65    auto Mwt = local_M.template block<cap_pressure_size, temperature_size>(
   67 
   68    auto Mcp = local_M.template block<contaminant_size, nonwet_pressure_size>(
   70    auto Mcpc = local_M.template block<contaminant_size, cap_pressure_size>(
   72    auto Mcx = local_M.template block<contaminant_size, contaminant_size>(
   74    auto Mct = local_M.template block<contaminant_size, temperature_size>(
   76 
   77    auto Mep = local_M.template block<temperature_size, nonwet_pressure_size>(
   79    auto Mepc = local_M.template block<temperature_size, cap_pressure_size>(
   81    auto Mex = local_M.template block<temperature_size, contaminant_size>(
   83    auto Met = local_M.template block<temperature_size, temperature_size>(
   85 
   87        NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
   88 
   89    auto Kap =
   90        local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
   92    auto Kapc = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
   94    auto Kax = local_K.template block<nonwet_pressure_size, contaminant_size>(
   96    auto Kat = local_K.template block<nonwet_pressure_size, temperature_size>(
   98 
   99    auto Kwp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
  101    auto Kwpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
  103    auto Kwx = local_K.template block<cap_pressure_size, contaminant_size>(
  105    auto Kwt = local_K.template block<cap_pressure_size, temperature_size>(
  107 
  108    auto Kcp = local_K.template block<contaminant_size, nonwet_pressure_size>(
  110    auto Kcpc = local_K.template block<contaminant_size, cap_pressure_size>(
  112    auto Kcx = local_K.template block<contaminant_size, contaminant_size>(
  114    auto Kct = local_K.template block<contaminant_size, temperature_size>(
  116 
  117    auto Kep = local_K.template block<temperature_size, nonwet_pressure_size>(
  119    auto Kepc = local_K.template block<temperature_size, cap_pressure_size>(
  121    auto Ket = local_K.template block<temperature_size, temperature_size>(
  123 
  124    auto Ba = local_b.template segment<nonwet_pressure_size>(
  126    auto Bw =
  128    auto Bc =
  130    auto Be =
  132 
  133    unsigned const n_integration_points =
  135 
  136    auto const num_nodes = ShapeFunction::NPOINTS;
  137    auto const pg_nodal_values =
  138        Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes);
  139    auto const pc_nodal_values =
  140        Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
  141 
  143 
  145        GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
  146 
  147    for (unsigned ip = 0; ip < n_integration_points; ip++)
  148    {
  150        auto const& 
N = ip_data.N;
 
  151        auto const& dNdx = ip_data.dNdx;
  152        auto const& w = ip_data.integration_weight;
  153        auto const& mass_operator = ip_data.mass_operator;
  154        auto const& diffusion_operator = ip_data.diffusion_operator;
  155        double pg_int_pt = 0.;
  156        double pc_int_pt = 0.;
  157        double Xc_int_pt = 0.;
  158        double T_int_pt = 0.;
  160                                         Xc_int_pt, T_int_pt);
  161 
  167 
  169        double const ideal_gas_constant_times_T_int_pt =
  174 
  175        auto const& medium =
  177        auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
 
  178        auto const& solid_phase = medium.phase("Solid");
  179        auto const& gas_phase = medium.phase("Gas");
  180 
  181        auto const& water_vapour_component = gas_phase.
component(
"w");
 
  182        auto const& dry_air_component = gas_phase.component("a");
  183        auto const& contaminant_vapour_component = gas_phase.component("c");
  184        auto const& dissolved_contaminant_component =
  185            liquid_phase.component("c");
  186 
  189                .template value<double>(vars, pos, t, dt);
  190 
  191        auto const water_mol_mass =
  192            water_vapour_component
  194                .template value<double>(vars, pos, t, dt);
  195        auto const air_mol_mass =
  196            dry_air_component
  198                .template value<double>(vars, pos, t, dt);
  199        auto const contaminant_mol_mass =
  200            contaminant_vapour_component
  202                .template value<double>(vars, pos, t, dt);
  203 
  204        double const Sw =
  206                .template value<double>(vars, pos, t, dt);
  207 
  210 
  211        double const dSw_dpc =
  213                .template dValue<double>(
  215                    pos, t, dt);
  216 
  217        auto const density_water =
  219                .template value<double>(vars, pos, t, dt);
  220 
  221        
  222        double const mol_density_water = density_water / water_mol_mass;
  223        double const mol_density_wet = mol_density_water;
  224 
  225        double const mol_density_nonwet =
  226            pg_int_pt / ideal_gas_constant_times_T_int_pt;
  227        double const mol_density_tot =
  228            Sw * mol_density_wet + (1 - Sw) * mol_density_nonwet;
  229 
  230        double const d_mol_density_nonwet_dpg =
  231            1 / ideal_gas_constant_times_T_int_pt;
  232        double const d_mol_density_nonwet_dT = -mol_density_nonwet / T_int_pt;
  233        double const d_mol_density_tot_dpc =
  234            (mol_density_wet - mol_density_nonwet) * dSw_dpc;
  235        double const d_mol_density_tot_dpg =
  236            (1 - Sw) * d_mol_density_nonwet_dpg;
  237        double const d_mol_density_tot_dT = (1 - Sw) * d_mol_density_nonwet_dT;
  238 
  239        
  240        double const latent_heat_evaporation =
  241            water_vapour_component
  242                .property(
  244                .template value<double>(vars, pos, t, dt);
  245 
  247 
  248        
  250        double const p_sat =
  251            water_vapour_component
  253                .template value<double>(vars, pos, t, dt);
  254        double const dp_sat_dT =
  255            water_vapour_component
  257                .template dValue<double>(
  259                    dt);
  260 
  261        
  262        double const K = std::exp(-pc_int_pt / mol_density_water /
  263                                  ideal_gas_constant_times_T_int_pt);
  264        double const dK_dT = pc_int_pt / mol_density_water /
  265                             ideal_gas_constant_times_T_int_pt / T_int_pt * K;
  266 
  267        
  268        
  269        double const p_vapour_nonwet = p_sat * K;
  270        double const d_p_vapour_nonwet_dT = dp_sat_dT * K + p_sat * dK_dT;
  271        double const d_p_vapour_nonwet_dpc =
  272            p_vapour_nonwet *
  273            (-1 / mol_density_water / ideal_gas_constant_times_T_int_pt);
  274 
  275        
  276        double const henry_contam =
  277            contaminant_vapour_component
  279                .template value<double>(vars, pos, t, dt);
  280        double d_henry_contaminant_dT =
  281            contaminant_vapour_component
  283                .template dValue<double>(
  285                    dt);
  286 
  287        
  288        double const k_c = pg_int_pt * henry_contam / mol_density_wet;
  289        double const k_w = pg_int_pt / p_vapour_nonwet;
  290 
  291        
  292        double const Ntot_c =
  293            Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
  294        
  295        double const x_contaminant_nonwet =
  296            Xc_int_pt * mol_density_tot / Ntot_c;
  297        double const x_contaminant_wet = k_c * x_contaminant_nonwet;
  298        double const x_water_wet = 1 - x_contaminant_wet;
  299        double const x_water_nonwet = x_water_wet / k_w;
  300        double const x_air_nonwet = 1 - x_water_nonwet - x_contaminant_nonwet;
  301 
  305 
  306        double const d_kc_dpg = henry_contam / mol_density_wet;
  307        double const d_kc_dT =
  308            pg_int_pt * d_henry_contaminant_dT / mol_density_wet;
  309 
  310        
  311        double const d_x_contaminant_nonwet_dpc =
  312            Xc_int_pt * (d_mol_density_tot_dpc / Ntot_c -
  313                         (mol_density_wet * k_c - mol_density_nonwet) *
  314                             dSw_dpc * mol_density_tot / Ntot_c / Ntot_c);
  315        double const d_x_contaminant_nonwet_dpg =
  316            Xc_int_pt * (d_mol_density_tot_dpg / Ntot_c -
  317                         (Sw * mol_density_wet * d_kc_dpg +
  318                          (1 - Sw) * d_mol_density_nonwet_dpg) *
  319                             mol_density_tot / Ntot_c / Ntot_c);
  320        double const d_x_contaminant_nonwet_dXc = mol_density_tot / Ntot_c;
  321        double const d_x_contaminant_nonwet_dT =
  322            Xc_int_pt * (d_mol_density_tot_dT / Ntot_c -
  323                         (Sw * mol_density_wet * d_kc_dT +
  324                          (1 - Sw) * d_mol_density_nonwet_dT) *
  325                             mol_density_tot / Ntot_c / Ntot_c);
  326 
  327        double const d_x_contaminant_wet_dpc = k_c * d_x_contaminant_nonwet_dpc;
  328        double const d_x_contaminant_wet_dpg =
  329            k_c * d_x_contaminant_nonwet_dpg + d_kc_dpg * x_contaminant_nonwet;
  330        double const d_x_contaminant_wet_dXc = k_c * d_x_contaminant_nonwet_dXc;
  331        double const d_x_contaminant_wet_dT =
  332            k_c * d_x_contaminant_nonwet_dT + d_kc_dT * x_contaminant_nonwet;
  333 
  334        double const d_x_water_wet_dpc = -d_x_contaminant_wet_dpc;
  335        double const d_x_water_wet_dpg = -d_x_contaminant_wet_dpg;
  336        double const d_x_water_wet_dXc = -d_x_contaminant_wet_dXc;
  337        double const d_x_water_wet_dT = -d_x_contaminant_wet_dT;
  338 
  339        double const d_x_water_nonwet_dpc =
  340            (d_p_vapour_nonwet_dpc * x_water_wet +
  341             p_vapour_nonwet * d_x_water_wet_dpc) /
  342            pg_int_pt;
  343        double const d_x_water_nonwet_dpg =
  344            p_vapour_nonwet * (d_x_water_wet_dpg / pg_int_pt -
  345                               x_water_wet / pg_int_pt / pg_int_pt);
  346        double const d_x_water_nonwet_dXc = d_x_water_wet_dXc / k_w;
  347        double const d_x_water_nonwet_dT =
  348            (d_p_vapour_nonwet_dT * x_water_wet +
  349             p_vapour_nonwet * d_x_water_wet_dT) /
  350            pg_int_pt;
  351 
  352        double const d_x_air_nonwet_dpc =
  353            -d_x_water_nonwet_dpc - d_x_contaminant_nonwet_dpc;
  354        double const d_x_air_nonwet_dpg =
  355            -d_x_water_nonwet_dpg - d_x_contaminant_nonwet_dpg;
  356        double const d_x_air_nonwet_dXc =
  357            -d_x_water_nonwet_dXc - d_x_contaminant_nonwet_dXc;
  358        double const d_x_air_nonwet_dT =
  359            -d_x_water_nonwet_dT - d_x_contaminant_nonwet_dT;
  360 
  361        
  362        double const density_contaminant_wet =
  363            mol_density_wet * contaminant_mol_mass * x_contaminant_wet;
  364        double const density_wet = density_water + density_contaminant_wet;
  365 
  366        double const density_air_nonwet =
  367            mol_density_nonwet * air_mol_mass * x_air_nonwet;
  368        double const density_water_nonwet =
  369            mol_density_nonwet * water_mol_mass * x_water_nonwet;
  370        double const density_contaminant_nonwet =
  371            mol_density_nonwet * contaminant_mol_mass * x_contaminant_nonwet;
  372        double const density_nonwet = density_air_nonwet +
  373                                      density_water_nonwet +
  374                                      density_contaminant_nonwet;
  375 
  376        auto const density_solid =
  378                .template value<double>(vars, pos, t, dt);
  379 
  380        
  381        double const d_density_nonwet_dpg =
  382            d_mol_density_nonwet_dpg *
  383                (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
  384                 contaminant_mol_mass * x_contaminant_nonwet) +
  385            mol_density_nonwet *
  386                (air_mol_mass * d_x_air_nonwet_dpg +
  387                 water_mol_mass * d_x_water_nonwet_dpg +
  388                 contaminant_mol_mass * d_x_contaminant_nonwet_dpg);
  389        double const d_density_nonwet_dpc =
  390            mol_density_nonwet *
  391            (air_mol_mass * d_x_air_nonwet_dpc +
  392             water_mol_mass * d_x_water_nonwet_dpc +
  393             contaminant_mol_mass * d_x_contaminant_nonwet_dpc);
  394        double const d_density_nonwet_dXc =
  395            mol_density_nonwet *
  396            (air_mol_mass * d_x_air_nonwet_dXc +
  397             water_mol_mass * d_x_water_nonwet_dXc +
  398             contaminant_mol_mass * d_x_contaminant_nonwet_dXc);
  399        double const d_density_nonwet_dT =
  400            d_mol_density_nonwet_dT *
  401                (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
  402                 contaminant_mol_mass * x_contaminant_nonwet) +
  403            mol_density_nonwet *
  404                (air_mol_mass * d_x_air_nonwet_dT +
  405                 water_mol_mass * d_x_water_nonwet_dT +
  406                 contaminant_mol_mass * d_x_contaminant_nonwet_dT);
  407 
  408        double const mol_mass_nonwet =
  409            x_water_nonwet * water_mol_mass + x_air_nonwet * air_mol_mass +
  410            x_contaminant_nonwet * contaminant_mol_mass;
  411        
  412        double const X_water_nonwet =
  413            x_water_nonwet * water_mol_mass / mol_mass_nonwet;
  414        double const X_air_nonwet =
  415            x_air_nonwet * air_mol_mass / mol_mass_nonwet;
  416        double const X_contaminant_nonwet = 1 - X_water_nonwet - X_air_nonwet;
  417 
  418        
  419        double const heat_capacity_dry_air =
  420            dry_air_component
  421                .property(
  423                .template value<double>(vars, pos, t, dt);
  424        double const heat_capacity_water_vapour =
  425            water_vapour_component
  426                .property(
  428                .template value<double>(vars, pos, t, dt);
  429        double const heat_capacity_contaminant_vapour =
  430            contaminant_vapour_component
  431                .property(
  433                .template value<double>(vars, pos, t, dt);
  434 
  435        double const heat_capacity_water =
  436            liquid_phase
  437                .property(
  439                .template value<double>(vars, pos, t, dt);
  440        double const heat_capacity_solid =
  441            solid_phase
  442                .property(
  444                .template value<double>(vars, pos, t, dt);
  445 
  446        
  447        
  448        
  449        double const enthalpy_air_nonwet =
  451            IdealGasConstant * T_int_pt / air_mol_mass;
  452        double const enthalpy_water_nonwet =
  454            latent_heat_evaporation;
  455        double const enthalpy_contaminant_nonwet =
  456            heat_capacity_contaminant_vapour *
  458            IdealGasConstant * T_int_pt / contaminant_mol_mass;
  459 
  460        
  461        double const enthalpy_nonwet =
  462            enthalpy_air_nonwet * X_air_nonwet +
  463            enthalpy_water_nonwet * X_water_nonwet +
  464            enthalpy_contaminant_nonwet * X_contaminant_nonwet;
  465        double const enthalpy_wet =
  467 
  468        
  469        double const internal_energy_nonwet =
  470            enthalpy_nonwet - pg_int_pt / density_nonwet;
  471        double const internal_energy_wet = enthalpy_wet;
  472 
  473        
  474        double const d_enthalpy_air_nonwet_dT =
  476        double const d_enthalpy_contaminant_nonwet_dT =
  477            heat_capacity_contaminant_vapour +
  479 
  480        double const d_enthalpy_nonwet_dT =
  481            heat_capacity_water * X_water_nonwet +
  482            d_enthalpy_air_nonwet_dT * X_air_nonwet +
  483            d_enthalpy_contaminant_nonwet_dT * X_contaminant_nonwet;
  484 
  485        
  486        Map.noalias() += 
porosity * (1 - Sw) *
 
  487                         (mol_density_nonwet * d_x_air_nonwet_dpg +
  488                          x_air_nonwet * d_mol_density_nonwet_dpg) *
  489                         mass_operator;
  490        Mapc.noalias() +=
  492            ((1 - Sw) * d_x_air_nonwet_dpc - x_air_nonwet * dSw_dpc) *
  493            mass_operator;
  494        Max.noalias() += 
porosity * mol_density_nonwet * (1 - Sw) *
 
  495                         d_x_air_nonwet_dXc * mass_operator;
  496        Mat.noalias() += 
porosity * (1 - Sw) *
 
  497                         (mol_density_nonwet * d_x_air_nonwet_dT +
  498                          x_air_nonwet * d_mol_density_nonwet_dT) *
  499                         mass_operator;
  500 
  501        Mwp.noalias() +=
  503            (mol_density_wet * Sw * d_x_water_wet_dpg +
  504             (1 - Sw) * x_water_nonwet * d_mol_density_nonwet_dpg +
  505             mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dpg) *
  506            mass_operator;
  507        Mwpc.noalias() +=
  509            (mol_density_wet *
  510                 (x_water_wet * dSw_dpc + Sw * d_x_water_wet_dpc) +
  511             mol_density_nonwet *
  512                 ((1 - Sw) * d_x_water_nonwet_dpc - x_water_nonwet * dSw_dpc)) *
  513            mass_operator;
  514        Mwx.noalias() +=
  516            (mol_density_wet * Sw * d_x_water_wet_dXc +
  517             mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
  518            mass_operator;
  520                         ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
  521                          (d_p_vapour_nonwet_dT - p_vapour_nonwet / T_int_pt)) *
  522                         mass_operator;
  523 
  524        Mcp.noalias() +=
  526            (mol_density_wet * Sw * d_x_contaminant_wet_dpg +
  527             (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dpg +
  528             mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dpg) *
  529            mass_operator;
  530        Mcpc.noalias() +=
  532            (mol_density_wet *
  533                 (x_contaminant_wet * dSw_dpc + Sw * d_x_contaminant_wet_dpc) +
  534             mol_density_nonwet * ((1 - Sw) * d_x_contaminant_nonwet_dpc -
  535                                   x_contaminant_nonwet * dSw_dpc)) *
  536            mass_operator;
  537        Mcx.noalias() +=
  539            (mol_density_wet * Sw * d_x_contaminant_wet_dXc +
  540             mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dXc) *
  541            mass_operator;
  542        Mct.noalias() +=
  544            (mol_density_wet * Sw * d_x_contaminant_wet_dT +
  545             (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dT +
  546             mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dT) *
  547            mass_operator;
  548 
  550                         (d_density_nonwet_dpg * enthalpy_nonwet - 1) *
  551                         (1 - Sw) * mass_operator;
  553                              (density_wet * internal_energy_wet -
  554                               density_nonwet * internal_energy_nonwet) *
  555                              dSw_dpc * mass_operator +
  556                          porosity * d_density_nonwet_dpc * enthalpy_nonwet *
  557                              (1 - Sw) * mass_operator;
  558        Mex.noalias() += 
porosity * d_density_nonwet_dXc * enthalpy_nonwet *
 
  559                         (1 - Sw) * mass_operator;
  560        Met.noalias() +=
  561            ((1 - 
porosity) * density_solid * heat_capacity_solid +
 
  562             porosity * ((1 - Sw) * (d_density_nonwet_dT * enthalpy_nonwet +
  563                                     density_nonwet * d_enthalpy_nonwet_dT) +
  564                         Sw * density_wet * heat_capacity_water)) *
  565            mass_operator;
  566 
  567        
  568        double const diffusion_coeff_water_nonwet =
  569            water_vapour_component
  571                .template value<double>(vars, pos, t, dt);
  572        double const diffusion_coeff_contaminant_nonwet =
  573            contaminant_vapour_component
  575                .template value<double>(vars, pos, t, dt);
  576        double const diffusion_coeff_contaminant_wet =
  577            dissolved_contaminant_component
  579                .template value<double>(vars, pos, t, dt);
  580 
  581        
  582        double const k_rel_nonwet =
  583            medium
  586                .template value<double>(vars, pos, t, dt);
  587        auto const mu_nonwet =
  589                .template value<double>(vars, pos, t, dt);
  590        double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
  591 
  592        
  593        double const k_rel_wet =
  594            medium
  595                .property(
  597                .template value<double>(vars, pos, t, dt);
  598        auto const mu_wet =
  600                .template value<double>(vars, pos, t, dt);
  601        double const lambda_wet = k_rel_wet / mu_wet;
  602 
  603        
  607                    .value(vars, pos, t, dt));
  608 
  609        
  611 
  612        
  616                      -lambda_wet * permeability *
  617                      (dNdx * (pg_nodal_values - pc_nodal_values) -
  618                       density_wet * b))
  620                                      (pg_nodal_values - pc_nodal_values));
  621 
  622        laplace_operator.noalias() = dNdx.transpose() * 
permeability * dNdx * w;
 
  623 
  624        
  625        auto const solute_dispersivity_transverse =
  626            medium.template value<double>(
  628        auto const solute_dispersivity_longitudinal =
  629            medium.template value<double>(
  631 
  632        double const velocity_wet_magnitude = velocity_wet.norm();
  634            velocity_wet_magnitude != 0.0
  636                      (porosity * Sw * diffusion_coeff_contaminant_wet +
  637                       solute_dispersivity_transverse *
  638                           velocity_wet_magnitude) *
  639                          I +
  640                      (solute_dispersivity_longitudinal -
  641                       solute_dispersivity_transverse) /
  642                          velocity_wet_magnitude * velocity_wet *
  643                          velocity_wet.transpose())
  645                      (
porosity * Sw * diffusion_coeff_contaminant_wet +
 
  646                       solute_dispersivity_transverse *
  647                           velocity_wet_magnitude) *
  648                      I);
  649 
  650        auto const dispersion_operator =
  651            dNdx.transpose() * hydrodynamic_dispersion * dNdx * w;
  652 
  653        
  654        
  655        Kap.noalias() +=
  656            (mol_density_nonwet * x_air_nonwet * lambda_nonwet) *
  657                laplace_operator -
  658            (1 - Sw) * 
porosity * mol_density_nonwet *
 
  659                (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpg +
  660                 diffusion_coeff_contaminant_nonwet *
  661                     d_x_contaminant_nonwet_dpg) *
  662                diffusion_operator;
  663        Kapc.noalias() +=
  664            -(1 - Sw) * porosity * mol_density_nonwet *
  665            (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpc +
  666             diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dpc) *
  667            diffusion_operator;
  668        Kax.noalias() +=
  669            -(1 - Sw) * porosity * mol_density_nonwet *
  670            (diffusion_coeff_water_nonwet * d_x_water_nonwet_dXc +
  671             diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dXc) *
  672            diffusion_operator;
  673        Kat.noalias() +=
  674            -(1 - Sw) * porosity * mol_density_nonwet *
  675            (diffusion_coeff_water_nonwet * d_x_water_nonwet_dT +
  676             diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT) *
  677            diffusion_operator;
  678 
  679        Kwp.noalias() +=
  680            (mol_density_wet * x_water_wet * lambda_wet +
  681             mol_density_nonwet * x_water_nonwet * lambda_nonwet) *
  682                laplace_operator +
  683            porosity *
  684                (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
  685                     d_x_water_wet_dpg +
  686                 (1 - Sw) * diffusion_coeff_water_nonwet * mol_density_nonwet *
  687                     d_x_water_nonwet_dpg) *
  688                diffusion_operator;
  689        Kwpc.noalias() +=
  690            -mol_density_wet * x_water_wet * lambda_wet * laplace_operator +
  692                (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
  693                     d_x_water_wet_dpc +
  694                 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
  695                     d_x_water_nonwet_dpc) *
  696                diffusion_operator;
  697        Kwx.noalias() +=
  699            (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
  700                 d_x_water_wet_dXc +
  701             (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
  702                 d_x_water_nonwet_dXc) *
  703            diffusion_operator;
  704        Kwt.noalias() +=
  706            (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
  707                 d_x_water_wet_dT +
  708             (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
  709                 d_x_water_nonwet_dT) *
  710            diffusion_operator;
  711 
  712        Kcp.noalias() +=
  713            (mol_density_wet * x_contaminant_wet * lambda_wet +
  714             mol_density_nonwet * x_contaminant_nonwet * lambda_nonwet) *
  715                laplace_operator +
  716            mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpg +
  717            porosity * (1 - Sw) * mol_density_nonwet *
  718                diffusion_coeff_contaminant_nonwet *
  719                d_x_contaminant_nonwet_dpg * diffusion_operator;
  720        Kcpc.noalias() +=
  721            -mol_density_wet * x_contaminant_wet * lambda_wet *
  722                laplace_operator +
  723            mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpc +
  724            porosity * (1 - Sw) * mol_density_nonwet *
 
  725                diffusion_coeff_contaminant_nonwet *
  726                d_x_contaminant_nonwet_dpc * diffusion_operator;
  727        Kcx.noalias() +=
  728            mol_density_wet * dispersion_operator * d_x_contaminant_wet_dXc +
  729            porosity * (1 - Sw) * mol_density_nonwet *
 
  730                diffusion_coeff_contaminant_nonwet *
  731                d_x_contaminant_nonwet_dXc * diffusion_operator;
  732        Kct.noalias() +=
  733            mol_density_wet * dispersion_operator * d_x_contaminant_wet_dT +
  734            porosity * (1 - Sw) * mol_density_nonwet *
 
  735                diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT *
  736                diffusion_operator;
  737 
  738        Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
  739                          lambda_wet * density_wet * enthalpy_wet) *
  740                         laplace_operator;
  741        Kepc.noalias() +=
  742            -lambda_wet * enthalpy_wet * density_wet * laplace_operator;
  743 
  744        if (medium.hasProperty(
  746        {
  747            auto const lambda =
  748                medium
  749                    .property(
  751                    .value(vars, pos, t, dt);
  752 
  755 
  756            Ket.noalias() +=
  757                dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
  758        }
  759        else
  760        {
  761            auto const thermal_conductivity_solid =
  762                solid_phase
  763                    .property(
  765                    .value(vars, pos, t, dt);
  766 
  767            auto const thermal_conductivity_fluid =
  768                liquid_phase
  769                    .property(
  771                    .template value<double>(vars, pos, t, dt) *
  772                Sw;
  773 
  776                    GlobalDim>(thermal_conductivity_solid,
  777                               thermal_conductivity_fluid, 
porosity);
 
  778 
  779            Ket.noalias() +=
  780                dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
  781        }
  782 
  784        {
  787            Ba.noalias() += (mol_density_nonwet * x_air_nonwet * lambda_nonwet *
  788                             density_nonwet) *
  789                            gravity_operator;
  790            Bw.noalias() +=
  791                (mol_density_wet * x_water_wet * lambda_wet * density_wet +
  792                 mol_density_nonwet * x_water_nonwet * lambda_nonwet *
  793                     density_nonwet) *
  794                gravity_operator;
  795            Bc.noalias() += (mol_density_wet * x_contaminant_wet * lambda_wet *
  796                                 density_wet +
  797                             mol_density_nonwet * x_contaminant_nonwet *
  798                                 lambda_nonwet * density_nonwet) *
  799                            gravity_operator;
  800            Be.noalias() +=
  801                (lambda_nonwet * density_nonwet * density_nonwet *
  802                     enthalpy_nonwet +
  803                 lambda_wet * density_wet * density_wet * enthalpy_wet) *
  804                gravity_operator;
  805        }  
  806    }
  808    {
  809        Map = Map.colwise().sum().eval().asDiagonal();
  810        Mapc = Mapc.colwise().sum().eval().asDiagonal();
  811        Max = Max.colwise().sum().eval().asDiagonal();
  812        Mat = Mat.colwise().sum().eval().asDiagonal();
  813        Mwp = Mwp.colwise().sum().eval().asDiagonal();
  814        Mwpc = Mwpc.colwise().sum().eval().asDiagonal();
  815        Mwx = Mwx.colwise().sum().eval().asDiagonal();
  816        Mwt = Mwt.colwise().sum().eval().asDiagonal();
  817        Mcp = Mcp.colwise().sum().eval().asDiagonal();
  818        Mcpc = Mcpc.colwise().sum().eval().asDiagonal();
  819        Mcx = Mcx.colwise().sum().eval().asDiagonal();
  820        Mct = Mct.colwise().sum().eval().asDiagonal();
  821        Mep = Mep.colwise().sum().eval().asDiagonal();
  822        Mepc = Mepc.colwise().sum().eval().asDiagonal();
  823        Mex = Mex.colwise().sum().eval().asDiagonal();
  824        Met = Met.colwise().sum().eval().asDiagonal();
  825    }  
  826}
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
Component const & component(std::size_t const &index) const
double gas_phase_pressure
double capillary_pressure
double enthalpy_of_evaporation
std::size_t getID() const
Returns the ID of the element.
static const int temperature_matrix_index
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
static const int nonwet_pressure_matrix_index
static const int contaminant_matrix_index
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
static const int cap_pressure_matrix_index
typename ShapeMatricesType::NodalVectorType NodalVectorType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
constexpr double CelsiusZeroInKelvin
Zero degrees Celsius in Kelvin.
constexpr double IdealGasConstant
Eigen::Matrix< double, GlobalDim, GlobalDim > formEffectiveThermalConductivity(MaterialPropertyLib::PropertyDataType const &solid_thermal_conductivity, const double fluid_thermal_conductivity, const double porosity)
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
@ relative_permeability_nonwetting_phase
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
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 &)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
bool const has_mass_lumping
Eigen::VectorXd const specific_body_force
MaterialPropertyLib::MaterialSpatialDistributionMap media_map