Loading [MathJax]/extensions/tex2jax.js
OGS
ComponentTransportFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <numeric>
14#include <vector>
15
36
37namespace ProcessLib
38{
39namespace ComponentTransport
40{
41template <typename GlobalDimNodalMatrixType>
43{
44 IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_,
45 double const& integration_weight_)
46 : dNdx(dNdx_), integration_weight(integration_weight_)
47 {
48 }
49
51 GlobalDimNodalMatrixType const dNdx;
52 double const integration_weight;
53
54 // -1 indicates that no chemical reaction takes place in the element to
55 // which the integration point belongs.
57
58 double porosity = std::numeric_limits<double>::quiet_NaN();
59 double porosity_prev = std::numeric_limits<double>::quiet_NaN();
61};
62
66{
67public:
69
70 virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
71
73 std::size_t const mesh_item_id,
74 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
75 std::vector<GlobalVector*> const& x, double const t)
76 {
77 std::vector<double> local_x_vec;
78
79 auto const n_processes = x.size();
80 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
81 {
82 auto const indices =
83 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
84 assert(!indices.empty());
85 auto const local_solution = x[process_id]->get(indices);
86 local_x_vec.insert(std::end(local_x_vec),
87 std::begin(local_solution),
88 std::end(local_solution));
89 }
90 auto const local_x = MathLib::toVector(local_x_vec);
91
93 }
94
96 std::size_t const mesh_item_id,
97 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
98 std::vector<GlobalVector*> const& x, double const t, double const dt)
99 {
100 std::vector<double> local_x_vec;
101
102 auto const n_processes = x.size();
103 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
104 {
105 auto const indices =
106 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
107 assert(!indices.empty());
108 auto const local_solution = x[process_id]->get(indices);
109 local_x_vec.insert(std::end(local_x_vec),
110 std::begin(local_solution),
111 std::end(local_solution));
112 }
113 auto const local_x = MathLib::toVector(local_x_vec);
114
115 setChemicalSystemConcrete(local_x, t, dt);
116 }
117
119 std::size_t const mesh_item_id,
120 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
121 std::vector<GlobalVector*> const& x, double const t, double const dt,
122 GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, int const process_id)
123 {
124 std::vector<double> local_x_vec;
125
126 auto const n_processes = x.size();
127 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++pcs_id)
128 {
129 auto const indices =
130 NumLib::getIndices(mesh_item_id, *dof_tables[pcs_id]);
131 assert(!indices.empty());
132 auto const local_solution = x[pcs_id]->get(indices);
133 local_x_vec.insert(std::end(local_x_vec),
134 std::begin(local_solution),
135 std::end(local_solution));
136 }
137 auto const local_x = MathLib::toVector(local_x_vec);
138
139 auto const indices =
140 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
141 auto const num_r_c = indices.size();
142
143 std::vector<double> local_M_data;
144 local_M_data.reserve(num_r_c * num_r_c);
145 std::vector<double> local_K_data;
146 local_K_data.reserve(num_r_c * num_r_c);
147 std::vector<double> local_b_data;
148 local_b_data.reserve(num_r_c);
149
150 assembleReactionEquationConcrete(t, dt, local_x, local_M_data,
151 local_K_data, local_b_data,
152 process_id);
153
154 auto const r_c_indices =
156 if (!local_M_data.empty())
157 {
158 auto const local_M =
159 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
160 M.add(r_c_indices, local_M);
161 }
162 if (!local_K_data.empty())
163 {
164 auto const local_K =
165 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
166 K.add(r_c_indices, local_K);
167 }
168 if (!local_b_data.empty())
169 {
170 b.add(indices, local_b_data);
171 }
172 }
173
174 virtual void postSpeciationCalculation(std::size_t const ele_id,
175 double const t, double const dt) = 0;
176
178 std::size_t const ele_id) = 0;
179
180 virtual std::vector<double> const& getIntPtDarcyVelocity(
181 const double t,
182 std::vector<GlobalVector*> const& x,
183 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
184 std::vector<double>& cache) const = 0;
185
186 virtual std::vector<double> const& getIntPtMolarFlux(
187 const double t, std::vector<GlobalVector*> const& x,
188 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
189 std::vector<double>& cache, int const component_id) const = 0;
190
191private:
193 Eigen::VectorXd const& /*local_x*/, double const /*t*/) = 0;
194
195 virtual void setChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/,
196 double const /*t*/,
197 double const /*dt*/) = 0;
198
200 double const t, double const dt, Eigen::VectorXd const& local_x,
201 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
202 std::vector<double>& local_b_data, int const transport_process_id) = 0;
203};
204
205template <typename ShapeFunction, int GlobalDim>
207{
208 // When monolithic scheme is adopted, nodal pressure and nodal concentration
209 // are accessed by vector index.
210 static const int pressure_index = 0;
211 const int temperature_index = -1;
213
214 static const int pressure_size = ShapeFunction::NPOINTS;
215 static const int temperature_size = ShapeFunction::NPOINTS;
216 static const int concentration_size =
217 ShapeFunction::NPOINTS; // per component
218
221
223 typename ShapeMatricesType::template MatrixType<pressure_size,
226 typename ShapeMatricesType::template VectorType<pressure_size>;
227
229 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
230 using LocalVectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
231
234
239
240public:
242 MeshLib::Element const& element,
243 std::size_t const local_matrix_size,
244 NumLib::GenericIntegrationMethod const& integration_method,
245 bool is_axially_symmetric,
246 ComponentTransportProcessData const& process_data,
247 std::vector<std::reference_wrapper<ProcessVariable>> const&
248 transport_process_variables)
249 : temperature_index(process_data.isothermal ? -1
250 : ShapeFunction::NPOINTS),
251 first_concentration_index(process_data.isothermal
252 ? ShapeFunction::NPOINTS
253 : 2 * ShapeFunction::NPOINTS),
254 _element(element),
255 _process_data(process_data),
256 _integration_method(integration_method),
257 _transport_process_variables(transport_process_variables)
258 {
259 (void)local_matrix_size;
260
261 unsigned const n_integration_points =
263 _ip_data.reserve(n_integration_points);
264
267
268 double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
269
270 auto const shape_matrices =
272 GlobalDim>(element, is_axially_symmetric,
274 auto const& medium =
276 for (unsigned ip = 0; ip < n_integration_points; ip++)
277 {
278 _ip_data.emplace_back(
279 shape_matrices[ip].dNdx,
281 shape_matrices[ip].integralMeasure *
282 shape_matrices[ip].detJ * aperture_size);
283
284 _ip_data[ip].porosity =
286 .template initialValue<double>(
287 pos, std::numeric_limits<double>::quiet_NaN() /*t*/);
288
289 _ip_data[ip].pushBackState();
290 }
291 }
292
293 void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
294 {
296 // chemical system index map
297 auto& chemical_system_index_map =
299
300 unsigned const n_integration_points =
302 for (unsigned ip = 0; ip < n_integration_points; ip++)
303 {
304 _ip_data[ip].chemical_system_id =
305 chemical_system_index_map.empty()
306 ? 0
307 : chemical_system_index_map.back() + 1;
308 chemical_system_index_map.push_back(
309 _ip_data[ip].chemical_system_id);
310 }
311 }
312
313 void initializeChemicalSystemConcrete(Eigen::VectorXd const& local_x,
314 double const t) override
315 {
317
318 auto const& medium =
320
323
324 auto const& Ns =
326 .NsHigherOrder<typename ShapeFunction::MeshElement>();
327
328 unsigned const n_integration_points =
330
331 for (unsigned ip = 0; ip < n_integration_points; ip++)
332 {
333 auto& ip_data = _ip_data[ip];
334 auto const& N = Ns[ip];
335 auto const& chemical_system_id = ip_data.chemical_system_id;
336
337 // set position with N as the shape matrix at the current
338 // integration point
340 NumLib::interpolateCoordinates<ShapeFunction,
342 N)));
343
344 auto const n_component = _transport_process_variables.size();
345 std::vector<double> C_int_pt(n_component);
346 for (unsigned component_id = 0; component_id < n_component;
347 ++component_id)
348 {
349 auto const concentration_index =
351 component_id * concentration_size;
352 auto const local_C =
353 local_x.template segment<concentration_size>(
354 concentration_index);
355
357 C_int_pt[component_id]);
358 }
359
361 ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
362 medium, pos, t);
363 }
364 }
365
366 void setChemicalSystemConcrete(Eigen::VectorXd const& local_x,
367 double const t, double dt) override
368 {
370
371 auto const& medium =
373
376
379
380 auto const& Ns =
382 .NsHigherOrder<typename ShapeFunction::MeshElement>();
383
384 unsigned const n_integration_points =
386
387 for (unsigned ip = 0; ip < n_integration_points; ip++)
388 {
389 auto& ip_data = _ip_data[ip];
390 auto const& N = Ns[ip];
391 auto& porosity = ip_data.porosity;
392 auto const& porosity_prev = ip_data.porosity_prev;
393 auto const& chemical_system_id = ip_data.chemical_system_id;
394
395 auto const n_component = _transport_process_variables.size();
396 std::vector<double> C_int_pt(n_component);
397 for (unsigned component_id = 0; component_id < n_component;
398 ++component_id)
399 {
400 auto const concentration_index =
402 component_id * concentration_size;
403 auto const local_C =
404 local_x.template segment<concentration_size>(
405 concentration_index);
406
408 C_int_pt[component_id]);
409 }
410
411 {
412 vars_prev.porosity = porosity_prev;
413
414 porosity =
416 ? porosity_prev
417 : medium
418 ->property(
420 .template value<double>(vars, vars_prev, pos, t,
421 dt);
422
423 vars.porosity = porosity;
424 }
425
427 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
428 }
429 }
430
431 void postSpeciationCalculation(std::size_t const ele_id, double const t,
432 double const dt) override
433 {
435 {
436 return;
437 }
438
439 auto const& medium = *_process_data.media_map.getMedium(ele_id);
440
442 pos.setElementID(ele_id);
443
444 for (auto& ip_data : _ip_data)
445 {
446 ip_data.porosity = ip_data.porosity_prev;
447
449 ->updateVolumeFractionPostReaction(ip_data.chemical_system_id,
450 medium, pos,
451 ip_data.porosity, t, dt);
452
454 ip_data.chemical_system_id, medium, ip_data.porosity);
455 }
456 }
457
458 void assemble(double const t, double const dt,
459 std::vector<double> const& local_x,
460 std::vector<double> const& /*local_x_prev*/,
461 std::vector<double>& local_M_data,
462 std::vector<double>& local_K_data,
463 std::vector<double>& local_b_data) override
464 {
465 auto const local_matrix_size = local_x.size();
466 // Nodal DOFs include pressure
467 int const num_nodal_dof = 1 + _transport_process_variables.size();
468 // This assertion is valid only if all nodal d.o.f. use the same shape
469 // matrices.
470 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
471
473 local_M_data, local_matrix_size, local_matrix_size);
475 local_K_data, local_matrix_size, local_matrix_size);
477 local_b_data, local_matrix_size);
478
479 // Get block matrices
480 auto Kpp = local_K.template block<pressure_size, pressure_size>(
482 auto Mpp = local_M.template block<pressure_size, pressure_size>(
484 auto Bp = local_b.template segment<pressure_size>(pressure_index);
485
486 auto local_p = Eigen::Map<const NodalVectorType>(
487 &local_x[pressure_index], pressure_size);
488
489 auto const& b =
492
493 auto const number_of_components = num_nodal_dof - 1;
494 for (int component_id = 0; component_id < number_of_components;
495 ++component_id)
496 {
497 /* Partitioned assembler matrix
498 * | pp | pc1 | pc2 | pc3 |
499 * |-----|-----|-----|-----|
500 * | c1p | c1c1| 0 | 0 |
501 * |-----|-----|-----|-----|
502 * | c2p | 0 | c2c2| 0 |
503 * |-----|-----|-----|-----|
504 * | c3p | 0 | 0 | c3c3|
505 */
506 auto concentration_index =
507 pressure_size + component_id * concentration_size;
508
509 auto KCC =
510 local_K.template block<concentration_size, concentration_size>(
511 concentration_index, concentration_index);
512 auto MCC =
513 local_M.template block<concentration_size, concentration_size>(
514 concentration_index, concentration_index);
515 auto MCp =
516 local_M.template block<concentration_size, pressure_size>(
517 concentration_index, pressure_index);
518 auto MpC =
519 local_M.template block<pressure_size, concentration_size>(
520 pressure_index, concentration_index);
521
522 auto local_C = Eigen::Map<const NodalVectorType>(
523 &local_x[concentration_index], concentration_size);
524
525 assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
526 MCC, MCp, MpC, Kpp, Mpp, Bp);
527
529 {
530 auto const stoichiometric_matrix =
533
534 assert(stoichiometric_matrix);
535
536 for (Eigen::SparseMatrix<double>::InnerIterator it(
537 *stoichiometric_matrix, component_id);
538 it;
539 ++it)
540 {
541 auto const stoichiometric_coefficient = it.value();
542 auto const coupled_component_id = it.row();
543 auto const kinetic_prefactor =
545 ->getKineticPrefactor(coupled_component_id);
546
547 auto const concentration_index =
548 pressure_size + component_id * concentration_size;
549 auto const coupled_concentration_index =
551 coupled_component_id * concentration_size;
552 auto KCmCn = local_K.template block<concentration_size,
554 concentration_index, coupled_concentration_index);
555
556 // account for the coupling between components
557 assembleKCmCn(component_id, t, dt, KCmCn,
558 stoichiometric_coefficient,
559 kinetic_prefactor);
560 }
561 }
562 }
563 }
564
566 GlobalDimVectorType const& b, int const component_id, double const t,
567 double const dt,
568 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
569 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
570 Eigen::Ref<LocalBlockMatrixType> KCC,
571 Eigen::Ref<LocalBlockMatrixType> MCC,
572 Eigen::Ref<LocalBlockMatrixType> MCp,
573 Eigen::Ref<LocalBlockMatrixType> MpC,
574 Eigen::Ref<LocalBlockMatrixType> Kpp,
575 Eigen::Ref<LocalBlockMatrixType> Mpp,
576 Eigen::Ref<LocalSegmentVectorType> Bp)
577 {
578 unsigned const n_integration_points =
580
583
585
586 // Get material properties
587 auto const& medium =
589 // Select the only valid for component transport liquid phase.
590 auto const& phase = medium.phase("AqueousLiquid");
591
592 // Assume that the component name is the same as the process variable
593 // name. Components are shifted by one because the first one is always
594 // pressure.
595 auto const& component = phase.component(
596 _transport_process_variables[component_id].get().getName());
597
598 LocalBlockMatrixType KCC_Laplacian =
599 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
600
601 std::vector<GlobalDimVectorType> ip_flux_vector;
602 double average_velocity_norm = 0.0;
604 {
605 ip_flux_vector.reserve(n_integration_points);
606 }
607
608 auto const& Ns =
610 .NsHigherOrder<typename ShapeFunction::MeshElement>();
611
612 for (unsigned ip(0); ip < n_integration_points; ++ip)
613 {
614 auto& ip_data = _ip_data[ip];
615 auto const& dNdx = ip_data.dNdx;
616 auto const& N = Ns[ip];
617 auto const& w = ip_data.integration_weight;
618 auto& porosity = ip_data.porosity;
619
620 double C_int_pt = 0.0;
621 double p_int_pt = 0.0;
622
623 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
624 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
625
626 // set position with N as the shape matrix at the current
627 // integration point
629 NumLib::interpolateCoordinates<ShapeFunction,
631 N)));
632
633 vars.concentration = C_int_pt;
634 vars.liquid_phase_pressure = p_int_pt;
635
636 // update according to a particular porosity model
638 .template value<double>(vars, pos, t, dt);
639 vars.porosity = porosity;
640
641 auto const& retardation_factor =
643 .template value<double>(vars, pos, t, dt);
644
645 auto const& solute_dispersivity_transverse = medium.template value<
646 double>(
648
649 auto const& solute_dispersivity_longitudinal =
650 medium.template value<double>(
652 longitudinal_dispersivity);
653
654 // Use the fluid density model to compute the density
655 // TODO (renchao): concentration of which component as the argument
656 // for calculation of fluid density
657 auto const density =
659 .template value<double>(vars, pos, t, dt);
660
661 auto const decay_rate =
663 .template value<double>(vars, pos, t, dt);
664
665 auto const& pore_diffusion_coefficient =
668 .value(vars, pos, t, dt));
669
672 vars, pos, t, dt));
673
674 // Use the viscosity model to compute the viscosity
676 .template value<double>(vars, pos, t, dt);
677
678 GlobalDimMatrixType const K_over_mu = K / mu;
679 GlobalDimVectorType const velocity =
681 ? GlobalDimVectorType(-K_over_mu *
682 (dNdx * p_nodal_values - density * b))
683 : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
684
685 const double drho_dp =
687 .template dValue<double>(
688 vars,
690 pos, t, dt);
691
692 const double drho_dC =
694 .template dValue<double>(
696 t, dt);
697
698 GlobalDimMatrixType const hydrodynamic_dispersion =
701 pore_diffusion_coefficient, velocity, porosity,
702 solute_dispersivity_transverse,
703 solute_dispersivity_longitudinal);
704
705 const double R_times_phi(retardation_factor * porosity);
706 GlobalDimVectorType const mass_density_flow = velocity * density;
707 auto const N_t_N = (N.transpose() * N).eval();
708
710 {
711 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
712 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
713 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
714 }
715 else
716 {
717 ip_flux_vector.emplace_back(mass_density_flow);
718 average_velocity_norm += velocity.norm();
719 }
720 MCC.noalias() += N_t_N * (R_times_phi * density * w);
721 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
722 KCC_Laplacian.noalias() +=
723 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
724
725 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
726
727 // Calculate Mpp, Kpp, and bp in the first loop over components
728 if (component_id == 0)
729 {
730 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
731 Kpp.noalias() +=
732 dNdx.transpose() * K_over_mu * dNdx * (density * w);
733
735 {
736 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
737 (density * density * w);
738 }
739 }
740 }
741
743 {
745 typename ShapeFunction::MeshElement>(
747 _ip_data,
749 ip_flux_vector,
750 average_velocity_norm /
751 static_cast<double>(n_integration_points),
752 KCC_Laplacian);
753 }
754
755 KCC.noalias() += KCC_Laplacian;
756 }
757
758 void assembleKCmCn(int const component_id, double const t, double const dt,
759 Eigen::Ref<LocalBlockMatrixType> KCmCn,
760 double const stoichiometric_coefficient,
761 double const kinetic_prefactor)
762 {
763 unsigned const n_integration_points =
765
768
770
771 auto const& medium =
773 auto const& phase = medium.phase("AqueousLiquid");
774 auto const& component = phase.component(
775 _transport_process_variables[component_id].get().getName());
776
777 auto const& Ns =
779 .NsHigherOrder<typename ShapeFunction::MeshElement>();
780
781 for (unsigned ip(0); ip < n_integration_points; ++ip)
782 {
783 auto& ip_data = _ip_data[ip];
784 auto const& w = ip_data.integration_weight;
785 auto const& N = Ns[ip];
786 auto& porosity = ip_data.porosity;
787
788 // set position with N as the shape matrix at the current
789 // integration point
791 NumLib::interpolateCoordinates<ShapeFunction,
793 N)));
794
795 auto const retardation_factor =
797 .template value<double>(vars, pos, t, dt);
798
800 .template value<double>(vars, pos, t, dt);
801
802 auto const density =
804 .template value<double>(vars, pos, t, dt);
805
806 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
807 kinetic_prefactor * retardation_factor *
808 porosity * density * N;
809 }
810 }
811
812 void assembleForStaggeredScheme(double const t, double const dt,
813 Eigen::VectorXd const& local_x,
814 Eigen::VectorXd const& local_x_prev,
815 int const process_id,
816 std::vector<double>& local_M_data,
817 std::vector<double>& local_K_data,
818 std::vector<double>& local_b_data) override
819 {
820 if (process_id == _process_data.hydraulic_process_id)
821 {
822 assembleHydraulicEquation(t, dt, local_x, local_x_prev,
823 local_M_data, local_K_data, local_b_data);
824 }
825 else if (process_id == _process_data.thermal_process_id)
826 {
827 assembleHeatTransportEquation(t, dt, local_x, local_x_prev,
828 local_M_data, local_K_data,
829 local_b_data);
830 }
831 else
832 {
833 // Go for assembling in an order of transport process id.
834 assembleComponentTransportEquation(t, dt, local_x, local_x_prev,
835 local_M_data, local_K_data,
836 local_b_data, process_id);
837 }
838 }
839
840 void assembleHydraulicEquation(double const t,
841 double const dt,
842 Eigen::VectorXd const& local_x,
843 Eigen::VectorXd const& local_x_prev,
844 std::vector<double>& local_M_data,
845 std::vector<double>& local_K_data,
846 std::vector<double>& local_b_data)
847 {
848 auto const local_p =
849 local_x.template segment<pressure_size>(pressure_index);
850 auto const local_C = local_x.template segment<concentration_size>(
852 auto const local_C_prev =
853 local_x_prev.segment<concentration_size>(first_concentration_index);
854
855 NodalVectorType local_T = getLocalTemperature(t, local_x);
856
858 local_M_data, pressure_size, pressure_size);
860 local_K_data, pressure_size, pressure_size);
862 local_b_data, pressure_size);
863
864 unsigned const n_integration_points =
866
869
870 auto const& b =
873
874 auto const& medium =
876 auto const& phase = medium.phase("AqueousLiquid");
877
880
881 auto const& Ns =
883 .NsHigherOrder<typename ShapeFunction::MeshElement>();
884
885 for (unsigned ip(0); ip < n_integration_points; ++ip)
886 {
887 auto& ip_data = _ip_data[ip];
888 auto const& dNdx = ip_data.dNdx;
889 auto const& w = ip_data.integration_weight;
890 auto const& N = Ns[ip];
891 auto& porosity = ip_data.porosity;
892 auto const& porosity_prev = ip_data.porosity_prev;
893
894 double const C_int_pt = N.dot(local_C);
895 double const p_int_pt = N.dot(local_p);
896 double const T_int_pt = N.dot(local_T);
897
898 vars.concentration = C_int_pt;
899 vars.liquid_phase_pressure = p_int_pt;
900 vars.temperature = T_int_pt;
901
902 // porosity
903 {
904 vars_prev.porosity = porosity_prev;
905
906 porosity =
908 ? porosity_prev
910 .template value<double>(vars, vars_prev, pos, t,
911 dt);
912
913 vars.porosity = porosity;
914 }
915
916 // Use the fluid density model to compute the density
917 // TODO: Concentration of which component as one of arguments for
918 // calculation of fluid density
919 auto const density =
921 .template value<double>(vars, pos, t, dt);
922
925 vars, pos, t, dt));
926
927 // Use the viscosity model to compute the viscosity
929 .template value<double>(vars, pos, t, dt);
930
931 GlobalDimMatrixType const K_over_mu = K / mu;
932
933 const double drho_dp =
935 .template dValue<double>(
936 vars,
938 pos, t, dt);
939 const double drho_dC =
941 .template dValue<double>(
943 t, dt);
944
945 // matrix assembly
946 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
947 local_K.noalias() +=
948 w * dNdx.transpose() * density * K_over_mu * dNdx;
949
951 {
952 local_b.noalias() +=
953 w * density * density * dNdx.transpose() * K_over_mu * b;
954 }
955
956 // coupling term
957 {
958 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
959
960 local_b.noalias() -=
961 w * N.transpose() * porosity * drho_dC * C_dot;
962 }
963 }
964 }
965
966 void assembleHeatTransportEquation(double const t, double const dt,
967 Eigen::VectorXd const& local_x,
968 Eigen::VectorXd const& /*local_x_prev*/,
969 std::vector<double>& local_M_data,
970 std::vector<double>& local_K_data,
971 std::vector<double>& /*local_b_data*/)
972 {
973 assert(local_x.size() ==
975
976 auto const local_p =
977 local_x.template segment<pressure_size>(pressure_index);
978 auto const local_T =
979 local_x.template segment<temperature_size>(temperature_index);
980
982 local_M_data, temperature_size, temperature_size);
984 local_K_data, temperature_size, temperature_size);
985
987 pos.setElementID(this->_element.getID());
988
989 auto const& process_data = this->_process_data;
990 auto const& medium =
991 *process_data.media_map.getMedium(this->_element.getID());
992 auto const& liquid_phase = medium.phase("AqueousLiquid");
993
994 auto const& b =
997
999
1000 unsigned const n_integration_points =
1002
1003 std::vector<GlobalDimVectorType> ip_flux_vector;
1004 double average_velocity_norm = 0.0;
1005 ip_flux_vector.reserve(n_integration_points);
1006
1007 auto const& Ns =
1009 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1010
1011 for (unsigned ip(0); ip < n_integration_points; ip++)
1012 {
1013 auto const& ip_data = this->_ip_data[ip];
1014 auto const& dNdx = ip_data.dNdx;
1015 auto const& w = ip_data.integration_weight;
1016 auto const& N = Ns[ip];
1017
1018 double p_at_xi = 0.;
1019 NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
1020 double T_at_xi = 0.;
1021 NumLib::shapeFunctionInterpolate(local_T, N, T_at_xi);
1022
1023 vars.temperature = T_at_xi;
1024 vars.liquid_phase_pressure = p_at_xi;
1025
1026 vars.liquid_saturation = 1.0;
1027
1028 auto const porosity =
1030 .template value<double>(vars, pos, t, dt);
1031 vars.porosity = porosity;
1032
1033 // Use the fluid density model to compute the density
1034 auto const fluid_density =
1035 liquid_phase
1037 .template value<double>(vars, pos, t, dt);
1038 vars.density = fluid_density;
1039 auto const specific_heat_capacity_fluid =
1040 liquid_phase
1042 .template value<double>(vars, pos, t, dt);
1043
1044 // Assemble mass matrix
1045 local_M.noalias() += w *
1047 vars, porosity, fluid_density,
1048 specific_heat_capacity_fluid, pos, t, dt) *
1049 N.transpose() * N;
1050
1051 // Assemble Laplace matrix
1052 auto const viscosity =
1053 liquid_phase
1055 .template value<double>(vars, pos, t, dt);
1056
1057 auto const intrinsic_permeability =
1059 medium
1060 .property(
1062 .value(vars, pos, t, dt));
1063
1064 GlobalDimMatrixType const K_over_mu =
1065 intrinsic_permeability / viscosity;
1066 GlobalDimVectorType const velocity =
1067 process_data.has_gravity
1068 ? GlobalDimVectorType(-K_over_mu *
1069 (dNdx * local_p - fluid_density * b))
1070 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1071
1072 GlobalDimMatrixType const thermal_conductivity_dispersivity =
1074 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1075 pos, t, dt);
1076
1077 local_K.noalias() +=
1078 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1079
1080 ip_flux_vector.emplace_back(velocity * fluid_density *
1081 specific_heat_capacity_fluid);
1082 average_velocity_norm += velocity.norm();
1083 }
1084
1086 process_data.stabilizer, this->_ip_data,
1087 _process_data.shape_matrix_cache, ip_flux_vector,
1088 average_velocity_norm / static_cast<double>(n_integration_points),
1089 local_K);
1090 }
1091
1093 double const t, double const dt, Eigen::VectorXd const& local_x,
1094 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_M_data,
1095 std::vector<double>& local_K_data,
1096 std::vector<double>& /*local_b_data*/, int const transport_process_id)
1097 {
1098 assert(static_cast<int>(local_x.size()) ==
1101 static_cast<int>(_transport_process_variables.size()) +
1103
1104 auto const local_p =
1105 local_x.template segment<pressure_size>(pressure_index);
1106
1107 NodalVectorType local_T = getLocalTemperature(t, local_x);
1108
1109 auto const local_C = local_x.template segment<concentration_size>(
1111 (transport_process_id - (_process_data.isothermal ? 1 : 2)) *
1113 auto const local_p_prev =
1114 local_x_prev.segment<pressure_size>(pressure_index);
1115
1117 local_M_data, concentration_size, concentration_size);
1119 local_K_data, concentration_size, concentration_size);
1120
1121 LocalBlockMatrixType KCC_Laplacian =
1122 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1123
1124 unsigned const n_integration_points =
1126
1127 std::vector<GlobalDimVectorType> ip_flux_vector;
1128 double average_velocity_norm = 0.0;
1130 {
1131 ip_flux_vector.reserve(n_integration_points);
1132 }
1133
1136
1137 auto const& b =
1140
1143
1144 auto const& medium =
1146 auto const& phase = medium.phase("AqueousLiquid");
1147 auto const component_id =
1148 transport_process_id - (_process_data.isothermal ? 1 : 2);
1149 auto const& component = phase.component(
1150 _transport_process_variables[component_id].get().getName());
1151
1152 auto const& Ns =
1154 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1155
1156 for (unsigned ip(0); ip < n_integration_points; ++ip)
1157 {
1158 auto& ip_data = _ip_data[ip];
1159 auto const& dNdx = ip_data.dNdx;
1160 auto const& w = ip_data.integration_weight;
1161 auto const& N = Ns[ip];
1162 auto& porosity = ip_data.porosity;
1163 auto const& porosity_prev = ip_data.porosity_prev;
1164
1165 double const C_int_pt = N.dot(local_C);
1166 double const p_int_pt = N.dot(local_p);
1167 double const T_int_pt = N.dot(local_T);
1168
1169 vars.concentration = C_int_pt;
1170 vars.liquid_phase_pressure = p_int_pt;
1171 vars.temperature = T_int_pt;
1172
1174 {
1175 vars.temperature = N.dot(local_T);
1176 }
1177
1178 // porosity
1179 {
1180 vars_prev.porosity = porosity_prev;
1181
1182 porosity =
1184 ? porosity_prev
1186 .template value<double>(vars, vars_prev, pos, t,
1187 dt);
1188
1189 vars.porosity = porosity;
1190 }
1191
1192 auto const& retardation_factor =
1194 .template value<double>(vars, pos, t, dt);
1195
1196 auto const& solute_dispersivity_transverse = medium.template value<
1197 double>(
1199 auto const& solute_dispersivity_longitudinal =
1200 medium.template value<double>(
1202 longitudinal_dispersivity);
1203
1204 // Use the fluid density model to compute the density
1205 auto const density =
1207 .template value<double>(vars, pos, t, dt);
1208 auto const decay_rate =
1210 .template value<double>(vars, pos, t, dt);
1211
1212 auto const& pore_diffusion_coefficient =
1215 .value(vars, pos, t, dt));
1216
1219 vars, pos, t, dt));
1220 // Use the viscosity model to compute the viscosity
1222 .template value<double>(vars, pos, t, dt);
1223
1224 GlobalDimMatrixType const K_over_mu = K / mu;
1225 GlobalDimVectorType const velocity =
1227 ? GlobalDimVectorType(-K_over_mu *
1228 (dNdx * local_p - density * b))
1229 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1230
1231 GlobalDimMatrixType const hydrodynamic_dispersion =
1234 pore_diffusion_coefficient, velocity, porosity,
1235 solute_dispersivity_transverse,
1236 solute_dispersivity_longitudinal);
1237
1238 double const R_times_phi = retardation_factor * porosity;
1239 auto const N_t_N = (N.transpose() * N).eval();
1240
1242 {
1243 const double drho_dC =
1245 .template dValue<double>(
1247 pos, t, dt);
1248 local_M.noalias() +=
1249 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1250 }
1251
1252 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1253
1254 // coupling term
1256 {
1257 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1258
1259 const double drho_dp =
1261 .template dValue<double>(vars,
1263 liquid_phase_pressure,
1264 pos, t, dt);
1265
1266 local_K.noalias() +=
1267 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1268 dNdx.transpose() * velocity * N * (density * w);
1269 }
1270 else
1271 {
1272 ip_flux_vector.emplace_back(velocity * density);
1273 average_velocity_norm += velocity.norm();
1274 }
1275 local_K.noalias() +=
1276 N_t_N * (decay_rate * R_times_phi * density * w);
1277
1278 KCC_Laplacian.noalias() += dNdx.transpose() *
1279 hydrodynamic_dispersion * dNdx *
1280 (density * w);
1281 }
1282
1284 {
1286 typename ShapeFunction::MeshElement>(
1288 _process_data.shape_matrix_cache, ip_flux_vector,
1289 average_velocity_norm /
1290 static_cast<double>(n_integration_points),
1291 KCC_Laplacian);
1292 }
1293 local_K.noalias() += KCC_Laplacian;
1294 }
1295
1297 double const t, double const dt, Eigen::VectorXd const& local_x,
1298 Eigen::VectorXd const& local_x_prev, int const process_id,
1299 std::vector<double>& local_b_data,
1300 std::vector<double>& local_Jac_data) override
1301 {
1302 if (process_id == _process_data.hydraulic_process_id)
1303 {
1304 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1305 local_b_data, local_Jac_data);
1306 }
1307 else
1308 {
1309 int const component_id = process_id - 1;
1311 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1312 component_id);
1313 }
1314 }
1315
1317 double const t, double const dt, Eigen::VectorXd const& local_x,
1318 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1319 std::vector<double>& local_Jac_data)
1320 {
1321 auto const p = local_x.template segment<pressure_size>(pressure_index);
1322 auto const c = local_x.template segment<concentration_size>(
1324
1325 auto const p_prev = local_x_prev.segment<pressure_size>(pressure_index);
1326 auto const c_prev =
1327 local_x_prev.segment<concentration_size>(first_concentration_index);
1328
1330 local_Jac_data, pressure_size, pressure_size);
1332 local_b_data, pressure_size);
1333
1334 unsigned const n_integration_points =
1336
1339 auto const& b =
1342
1343 auto const& medium =
1345 auto const& phase = medium.phase("AqueousLiquid");
1346
1349
1350 auto const& Ns =
1352 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1353
1354 for (unsigned ip(0); ip < n_integration_points; ++ip)
1355 {
1356 auto& ip_data = _ip_data[ip];
1357 auto const& dNdx = ip_data.dNdx;
1358 auto const& w = ip_data.integration_weight;
1359 auto const& N = Ns[ip];
1360 auto& phi = ip_data.porosity;
1361 auto const& phi_prev = ip_data.porosity_prev;
1362
1363 double const p_ip = N.dot(p);
1364 double const c_ip = N.dot(c);
1365
1366 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1367
1368 vars.liquid_phase_pressure = p_ip;
1369 vars.concentration = c_ip;
1370
1371 // porosity
1372 {
1373 vars_prev.porosity = phi_prev;
1374
1376 ? phi_prev
1378 .template value<double>(vars, vars_prev, pos, t,
1379 dt);
1380
1381 vars.porosity = phi;
1382 }
1383
1384 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1385 .template value<double>(vars, pos, t, dt);
1386
1389 vars, pos, t, dt));
1390
1392 .template value<double>(vars, pos, t, dt);
1393
1394 auto const drho_dp =
1396 .template dValue<double>(
1397 vars,
1399 pos, t, dt);
1400 auto const drho_dc =
1402 .template dValue<double>(
1404 t, dt);
1405
1406 // matrix assembly
1407 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1408 w * dNdx.transpose() * rho * k / mu * dNdx;
1409
1410 local_rhs.noalias() -=
1411 w * N.transpose() * phi *
1412 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1413 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1414
1416 {
1417 local_rhs.noalias() +=
1418 w * rho * dNdx.transpose() * k / mu * rho * b;
1419 }
1420 }
1421 }
1422
1424 double const t, double const dt, Eigen::VectorXd const& local_x,
1425 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1426 std::vector<double>& local_Jac_data, int const component_id)
1427 {
1428 auto const concentration_index =
1430
1431 auto const p = local_x.template segment<pressure_size>(pressure_index);
1432 auto const c =
1433 local_x.template segment<concentration_size>(concentration_index);
1434 auto const c_prev =
1435 local_x_prev.segment<concentration_size>(concentration_index);
1436
1439 {
1441 }
1442
1444 local_Jac_data, concentration_size, concentration_size);
1446 local_b_data, concentration_size);
1447
1448 LocalBlockMatrixType KCC_Laplacian =
1449 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1450
1451 unsigned const n_integration_points =
1453
1454 std::vector<GlobalDimVectorType> ip_flux_vector;
1455 double average_velocity_norm = 0.0;
1456 ip_flux_vector.reserve(n_integration_points);
1457
1460
1461 auto const& b =
1464
1467
1468 auto const& medium =
1470 auto const& phase = medium.phase("AqueousLiquid");
1471 auto const& component = phase.component(
1472 _transport_process_variables[component_id].get().getName());
1473
1474 auto const& Ns =
1476 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1477
1478 for (unsigned ip(0); ip < n_integration_points; ++ip)
1479 {
1480 auto& ip_data = _ip_data[ip];
1481 auto const& dNdx = ip_data.dNdx;
1482 auto const& w = ip_data.integration_weight;
1483 auto const& N = Ns[ip];
1484 auto& phi = ip_data.porosity;
1485 auto const& phi_prev = ip_data.porosity_prev;
1486
1487 double const p_ip = N.dot(p);
1488 double const c_ip = N.dot(c);
1489
1490 vars.liquid_phase_pressure = p_ip;
1491 vars.concentration = c_ip;
1492
1494 {
1495 vars.temperature = N.dot(T);
1496 }
1497
1498 // porosity
1499 {
1500 vars_prev.porosity = phi_prev;
1501
1503 ? phi_prev
1505 .template value<double>(vars, vars_prev, pos, t,
1506 dt);
1507
1508 vars.porosity = phi;
1509 }
1510
1511 auto const R =
1513 .template value<double>(vars, pos, t, dt);
1514
1515 auto const alpha_T = medium.template value<double>(
1517 auto const alpha_L = medium.template value<double>(
1519
1520 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1521 .template value<double>(vars, pos, t, dt);
1522 // first-order decay constant
1523 auto const alpha =
1525 .template value<double>(vars, pos, t, dt);
1526
1529 .value(vars, pos, t, dt));
1530
1533 vars, pos, t, dt));
1535 .template value<double>(vars, pos, t, dt);
1536 // Darcy flux
1537 GlobalDimVectorType const q =
1539 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1540 : GlobalDimVectorType(-k / mu * dNdx * p);
1541
1543 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1544 alpha_L);
1545
1546 // matrix assembly
1547 local_Jac.noalias() +=
1548 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1549
1550 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1551
1552 auto const cdot = (c - c_prev) / dt;
1553 local_rhs.noalias() -=
1554 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1555
1556 ip_flux_vector.emplace_back(q * rho);
1557 average_velocity_norm += q.norm();
1558 }
1559
1562 _process_data.shape_matrix_cache, ip_flux_vector,
1563 average_velocity_norm / static_cast<double>(n_integration_points),
1564 KCC_Laplacian);
1565
1566 local_rhs.noalias() -= KCC_Laplacian * c;
1567
1568 local_Jac.noalias() += KCC_Laplacian;
1569 }
1570
1572 double const t, double const dt, Eigen::VectorXd const& local_x,
1573 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1574 std::vector<double>& local_b_data,
1575 int const transport_process_id) override
1576 {
1577 auto const local_C = local_x.template segment<concentration_size>(
1579 (transport_process_id - 1) * concentration_size);
1580
1582 local_M_data, concentration_size, concentration_size);
1584 local_K_data, concentration_size, concentration_size);
1586 local_b_data, concentration_size);
1587
1588 unsigned const n_integration_points =
1590
1593
1596
1597 auto const& medium =
1599 auto const component_id = transport_process_id - 1;
1600
1601 auto const& Ns =
1603 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1604
1605 for (unsigned ip(0); ip < n_integration_points; ++ip)
1606 {
1607 auto& ip_data = _ip_data[ip];
1608 auto const w = ip_data.integration_weight;
1609 auto const& N = Ns[ip];
1610 auto& porosity = ip_data.porosity;
1611 auto const& porosity_prev = ip_data.porosity_prev;
1612 auto const chemical_system_id = ip_data.chemical_system_id;
1613
1614 double C_int_pt = 0.0;
1615 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1616
1617 vars.concentration = C_int_pt;
1618
1619 auto const porosity_dot = (porosity - porosity_prev) / dt;
1620
1621 // porosity
1622 {
1623 vars_prev.porosity = porosity_prev;
1624
1625 porosity =
1627 ? porosity_prev
1629 .template value<double>(vars, vars_prev, pos, t,
1630 dt);
1631 }
1632
1633 local_M.noalias() += w * N.transpose() * porosity * N;
1634
1635 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1636
1637 if (chemical_system_id == -1)
1638 {
1639 continue;
1640 }
1641
1642 auto const C_post_int_pt =
1644 component_id, chemical_system_id);
1645
1646 local_b.noalias() +=
1647 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1648 }
1649 }
1650
1651 std::vector<double> const& getIntPtDarcyVelocity(
1652 const double t,
1653 std::vector<GlobalVector*> const& x,
1654 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
1655 std::vector<double>& cache) const override
1656 {
1657 assert(x.size() == dof_table.size());
1658
1659 auto const n_processes = x.size();
1660 std::vector<std::vector<double>> local_x;
1661 local_x.reserve(n_processes);
1662
1663 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1664 {
1665 auto const indices =
1666 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1667 assert(!indices.empty());
1668 local_x.push_back(x[process_id]->get(indices));
1669 }
1670
1671 // only one process, must be monolithic.
1672 if (n_processes == 1)
1673 {
1674 auto const local_p = Eigen::Map<const NodalVectorType>(
1675 &local_x[0][pressure_index], pressure_size);
1676 auto const local_C = Eigen::Map<const NodalVectorType>(
1678 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1679 }
1680
1681 // multiple processes, must be staggered.
1682 {
1683 constexpr int pressure_process_id = 0;
1684 constexpr int concentration_process_id = 1;
1685 auto const local_p = Eigen::Map<const NodalVectorType>(
1686 &local_x[pressure_process_id][0], pressure_size);
1687 auto const local_C = Eigen::Map<const NodalVectorType>(
1688 &local_x[concentration_process_id][0], concentration_size);
1689 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1690 }
1691 }
1692
1693 std::vector<double> const& calculateIntPtDarcyVelocity(
1694 const double t,
1695 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
1696 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
1697 std::vector<double>& cache) const
1698 {
1699 auto const n_integration_points =
1701
1702 cache.clear();
1703 auto cache_mat = MathLib::createZeroedMatrix<
1704 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1705 cache, GlobalDim, n_integration_points);
1706
1709
1710 auto const& b =
1713
1715
1716 auto const& medium =
1718 auto const& phase = medium.phase("AqueousLiquid");
1719
1720 auto const& Ns =
1722 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1723
1724 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1725 {
1726 auto const& ip_data = _ip_data[ip];
1727 auto const& dNdx = ip_data.dNdx;
1728 auto const& N = Ns[ip];
1729 auto const& porosity = ip_data.porosity;
1730
1731 double C_int_pt = 0.0;
1732 double p_int_pt = 0.0;
1733
1734 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1735 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1736
1737 vars.concentration = C_int_pt;
1738 vars.liquid_phase_pressure = p_int_pt;
1739 vars.porosity = porosity;
1740
1741 // TODO (naumov) Temporary value not used by current material
1742 // models. Need extension of secondary variables interface.
1743 double const dt = std::numeric_limits<double>::quiet_NaN();
1746 vars, pos, t, dt));
1748 .template value<double>(vars, pos, t, dt);
1749 GlobalDimMatrixType const K_over_mu = K / mu;
1750
1751 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1753 {
1754 auto const rho_w =
1756 .template value<double>(vars, pos, t, dt);
1757 // here it is assumed that the vector b is directed 'downwards'
1758 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1759 }
1760 }
1761
1762 return cache;
1763 }
1764
1765 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
1766 const unsigned integration_point) const override
1767 {
1769 typename ShapeFunction::MeshElement>()[integration_point];
1770
1771 // assumes N is stored contiguously in memory
1772 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1773 }
1774
1775 Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
1776 double const t,
1777 std::vector<double> const& local_x) const override
1778 {
1779 auto const local_p = Eigen::Map<const NodalVectorType>(
1780 &local_x[pressure_index], pressure_size);
1781 auto const local_C = Eigen::Map<const NodalVectorType>(
1783
1784 // Eval shape matrices at given point
1785 // Note: Axial symmetry is set to false here, because we only need dNdx
1786 // here, which is not affected by axial symmetry.
1787 auto const shape_matrices =
1789 GlobalDim>(
1790 _element, false /*is_axially_symmetric*/,
1791 std::array{pnt_local_coords})[0];
1792
1795 auto const& b =
1798
1800
1801 auto const& medium =
1803 auto const& phase = medium.phase("AqueousLiquid");
1804
1805 // local_x contains the local concentration and pressure values
1806 double c_int_pt;
1807 NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
1808 vars.concentration = c_int_pt;
1809
1810 double p_int_pt;
1811 NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
1812 vars.liquid_phase_pressure = p_int_pt;
1813
1814 // TODO (naumov) Temporary value not used by current material models.
1815 // Need extension of secondary variables interface.
1816 double const dt = std::numeric_limits<double>::quiet_NaN();
1819 vars, pos, t, dt));
1820
1822 .template value<double>(vars, pos, t, dt);
1823 GlobalDimMatrixType const K_over_mu = K / mu;
1824
1825 GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
1826 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1827 .template value<double>(vars, pos, t, dt);
1829 {
1830 q += K_over_mu * rho_w * b;
1831 }
1832 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1833 flux.head<GlobalDim>() = rho_w * q;
1834 return flux;
1835 }
1836
1838 double const t,
1839 double const /*dt*/,
1840 Eigen::VectorXd const& local_x,
1841 Eigen::VectorXd const& /*local_x_prev*/) override
1842 {
1843 auto const local_p =
1844 local_x.template segment<pressure_size>(pressure_index);
1845 auto const local_C = local_x.template segment<concentration_size>(
1847
1848 std::vector<double> ele_velocity;
1849 calculateIntPtDarcyVelocity(t, local_p, local_C, ele_velocity);
1850
1851 auto const n_integration_points =
1853 auto const ele_velocity_mat =
1854 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
1855
1856 auto const ele_id = _element.getID();
1857 Eigen::Map<LocalVectorType>(
1858 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
1859 GlobalDim) =
1860 ele_velocity_mat.rowwise().sum() / n_integration_points;
1861 }
1862
1864 std::size_t const ele_id) override
1865 {
1866 auto const n_integration_points =
1868
1870 {
1871 auto const& medium = *_process_data.media_map.getMedium(ele_id);
1872
1873 for (auto& ip_data : _ip_data)
1874 {
1875 ip_data.porosity = ip_data.porosity_prev;
1876
1878 ->updatePorosityPostReaction(ip_data.chemical_system_id,
1879 medium, ip_data.porosity);
1880 }
1881
1883 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
1884 [](double const s, auto const& ip)
1885 { return s + ip.porosity; }) /
1886 n_integration_points;
1887 }
1888
1889 std::vector<GlobalIndexType> chemical_system_indices;
1890 chemical_system_indices.reserve(n_integration_points);
1891 std::transform(_ip_data.begin(), _ip_data.end(),
1892 std::back_inserter(chemical_system_indices),
1893 [](auto const& ip_data)
1894 { return ip_data.chemical_system_id; });
1895
1897 ele_id, chemical_system_indices);
1898 }
1899
1900 std::vector<double> const& getIntPtMolarFlux(
1901 const double t, std::vector<GlobalVector*> const& x,
1902 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
1903 std::vector<double>& cache, int const component_id) const override
1904 {
1905 std::vector<double> local_x_vec;
1906
1907 auto const n_processes = x.size();
1908 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1909 {
1910 auto const indices =
1911 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
1912 assert(!indices.empty());
1913 auto const local_solution = x[process_id]->get(indices);
1914 local_x_vec.insert(std::end(local_x_vec),
1915 std::begin(local_solution),
1916 std::end(local_solution));
1917 }
1918 auto const local_x = MathLib::toVector(local_x_vec);
1919
1920 auto const p = local_x.template segment<pressure_size>(pressure_index);
1921 auto const c = local_x.template segment<concentration_size>(
1923
1924 auto const n_integration_points =
1926
1927 cache.clear();
1928 auto cache_mat = MathLib::createZeroedMatrix<
1929 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1930 cache, GlobalDim, n_integration_points);
1931
1934
1935 auto const& b =
1938
1940
1941 auto const& medium =
1943 auto const& phase = medium.phase("AqueousLiquid");
1944
1945 auto const& component = phase.component(
1946 _transport_process_variables[component_id].get().getName());
1947
1948 auto const& Ns =
1950 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1951
1952 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1953 {
1954 auto const& ip_data = _ip_data[ip];
1955 auto const& dNdx = ip_data.dNdx;
1956 auto const& N = Ns[ip];
1957 auto const& phi = ip_data.porosity;
1958
1959 double const p_ip = N.dot(p);
1960 double const c_ip = N.dot(c);
1961
1962 vars.concentration = c_ip;
1963 vars.liquid_phase_pressure = p_ip;
1964 vars.porosity = phi;
1965
1966 double const dt = std::numeric_limits<double>::quiet_NaN();
1967
1970 vars, pos, t, dt));
1972 .template value<double>(vars, pos, t, dt);
1973 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1974 .template value<double>(vars, pos, t, dt);
1975
1976 // Darcy flux
1977 GlobalDimVectorType const q =
1979 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1980 : GlobalDimVectorType(-k / mu * dNdx * p);
1981
1982 auto const alpha_T = medium.template value<double>(
1984 auto const alpha_L = medium.template value<double>(
1988 .value(vars, pos, t, dt));
1989
1990 // Hydrodynamic dispersion
1992 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1993 alpha_L);
1994
1995 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
1996 }
1997
1998 return cache;
1999 }
2000
2001 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
2002 Eigen::VectorXd const& /*local_x_prev*/,
2003 double const /*t*/, double const /*dt*/,
2004 int const /*process_id*/) override
2005 {
2006 unsigned const n_integration_points =
2008
2009 for (unsigned ip = 0; ip < n_integration_points; ip++)
2010 {
2011 _ip_data[ip].pushBackState();
2012 }
2013 }
2014
2015private:
2018
2020 std::vector<std::reference_wrapper<ProcessVariable>> const
2022
2023 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>> _ip_data;
2024
2026 MaterialPropertyLib::VariableArray const& vars, const double porosity,
2027 const double fluid_density, const double specific_heat_capacity_fluid,
2028 ParameterLib::SpatialPosition const& pos, double const t,
2029 double const dt)
2030 {
2031 auto const& medium =
2032 *_process_data.media_map.getMedium(this->_element.getID());
2033 auto const& solid_phase = medium.phase("Solid");
2034
2035 auto const specific_heat_capacity_solid =
2036 solid_phase
2037 .property(
2039 .template value<double>(vars, pos, t, dt);
2040
2041 auto const solid_density =
2042 solid_phase.property(MaterialPropertyLib::PropertyType::density)
2043 .template value<double>(vars, pos, t, dt);
2044
2045 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2046 fluid_density * specific_heat_capacity_fluid * porosity;
2047 }
2048
2051 const double fluid_density, const double specific_heat_capacity_fluid,
2052 const GlobalDimVectorType& velocity,
2053 ParameterLib::SpatialPosition const& pos, double const t,
2054 double const dt)
2055 {
2056 auto const& medium =
2058
2059 auto thermal_conductivity =
2061 medium
2062 .property(
2064 .value(vars, pos, t, dt));
2065
2066 auto const thermal_dispersivity_transversal =
2067 medium
2069 thermal_transversal_dispersivity)
2070 .template value<double>();
2071
2072 auto const thermal_dispersivity_longitudinal =
2073 medium
2075 thermal_longitudinal_dispersivity)
2076 .template value<double>();
2077
2078 // Thermal conductivity is moved outside and zero matrix is passed
2079 // instead due to multiplication with fluid's density times specific
2080 // heat capacity.
2081 return thermal_conductivity +
2082 fluid_density * specific_heat_capacity_fluid *
2085 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2086 velocity, 0 /* phi */, thermal_dispersivity_transversal,
2087 thermal_dispersivity_longitudinal);
2088 }
2089
2091 Eigen::VectorXd const& local_x)
2092 {
2093 NodalVectorType local_T;
2095 {
2097 {
2099 _element, t);
2100 }
2101 else
2102 {
2103 local_T = NodalVectorType::Zero(temperature_size);
2104 }
2105 }
2106 else
2107 {
2108 local_T =
2109 local_x.template segment<temperature_size>(temperature_index);
2110 }
2111 return local_T;
2112 }
2113};
2114
2115} // namespace ComponentTransport
2116} // namespace ProcessLib
GlobalMatrix::IndexType GlobalIndexType
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
virtual void updatePorosityPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, double &)
virtual double getKineticPrefactor(std::size_t reaction_id) const
virtual void computeSecondaryVariable(std::size_t const, std::vector< GlobalIndexType > const &)
virtual void setChemicalSystemConcrete(std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const *, MaterialPropertyLib::VariableArray const &, ParameterLib::SpatialPosition const &, double const, double const)
virtual void updateVolumeFractionPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const, double const, double const)
std::vector< GlobalIndexType > chemical_system_index_map
virtual Eigen::SparseMatrix< double > const * getStoichiometricMatrix() const
virtual double getConcentration(int const, GlobalIndexType const) const
virtual void initializeChemicalSystemConcrete(std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const)
Phase const & phase(std::size_t index) const
Definition Medium.cpp:33
Property const & property(PropertyType const &p) const
Definition Phase.cpp:53
Component const & component(std::size_t const &index) const
Definition Phase.cpp:33
int add(IndexType row, IndexType col, double val)
Definition EigenMatrix.h:87
Global vector based on Eigen vector.
Definition EigenVector.h:26
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:77
constexpr double getWeight() const
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
auto const & NsHigherOrder() const
void setCoordinates(MathLib::Point3d const &coordinates)
void setElementID(std::size_t element_id)
virtual std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual void setChemicalSystemConcrete(Eigen::VectorXd const &, double const, double const)=0
virtual void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt)=0
virtual void computeReactionRelatedSecondaryVariable(std::size_t const ele_id)=0
void initializeChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t)
void setChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
virtual void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id)=0
void assembleReactionEquation(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, int const process_id)
virtual void initializeChemicalSystemConcrete(Eigen::VectorXd const &, double const)=0
virtual std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache, int const component_id) const =0
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
void assembleWithJacobianHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleHeatTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > LocalMatrixType
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< double > &cache, int const component_id) const override
void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt) override
NumLib::GenericIntegrationMethod const & _integration_method
void initializeChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t) override
void computeReactionRelatedSecondaryVariable(std::size_t const ele_id) override
std::vector< double > const & calculateIntPtDarcyVelocity(const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, std::vector< double > &cache) const
void assembleBlockMatrices(GlobalDimVectorType const &b, int const component_id, double const t, double const dt, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< LocalBlockMatrixType > KCC, Eigen::Ref< LocalBlockMatrixType > MCC, Eigen::Ref< LocalBlockMatrixType > MCp, Eigen::Ref< LocalBlockMatrixType > MpC, Eigen::Ref< LocalBlockMatrixType > Kpp, Eigen::Ref< LocalBlockMatrixType > Mpp, Eigen::Ref< LocalSegmentVectorType > Bp)
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
ComponentTransportProcessData const & _process_data
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
void computeSecondaryVariableConcrete(double const t, double const, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
std::vector< std::reference_wrapper< ProcessVariable > > const _transport_process_variables
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void assembleComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
NodalVectorType getLocalTemperature(double const t, Eigen::VectorXd const &local_x)
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
GlobalDimMatrixType getThermalConductivityDispersivity(MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void assembleWithJacobianComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, int const component_id)
Eigen::Matrix< double, Eigen::Dynamic, 1 > LocalVectorType
double getHeatEnergyCoefficient(MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id) override
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, ComponentTransportProcessData const &process_data, std::vector< std::reference_wrapper< ProcessVariable > > const &transport_process_variables)
typename ShapeMatricesType::template VectorType< pressure_size > LocalSegmentVectorType
void assembleKCmCn(int const component_id, double const t, double const dt, Eigen::Ref< LocalBlockMatrixType > KCmCn, double const stoichiometric_coefficient, double const kinetic_prefactor)
void setChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t, double dt) override
typename ShapeMatricesType::NodalVectorType NodalVectorType
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ retardation_factor
specify retardation factor used in component transport process.
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
void assembleAdvectionMatrix(IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
virtual Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > getNodalValuesOnElement(MeshLib::Element const &element, double const t) const
Returns a matrix of values for all nodes of the given element.
Definition Parameter.h:164
std::vector< Eigen::VectorXd > const projected_specific_body_force_vectors
Projected specific body force vector: R * R^T * b.
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
NumLib::ShapeMatrixCache shape_matrix_cache
caches for each mesh element type the shape matrix
IntegrationPointData(GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)