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 auto const n_component = _transport_process_variables.size();
338 std::vector<double> C_int_pt(n_component);
339 for (unsigned component_id = 0; component_id < n_component;
340 ++component_id)
341 {
342 auto const concentration_index =
344 component_id * concentration_size;
345 auto const local_C =
346 local_x.template segment<concentration_size>(
347 concentration_index);
348
350 C_int_pt[component_id]);
351 }
352
354 ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
355 medium, pos, t);
356 }
357 }
358
359 void setChemicalSystemConcrete(Eigen::VectorXd const& local_x,
360 double const t, double dt) override
361 {
363
364 auto const& medium =
366
369
372
373 auto const& Ns =
375 .NsHigherOrder<typename ShapeFunction::MeshElement>();
376
377 unsigned const n_integration_points =
379
380 for (unsigned ip = 0; ip < n_integration_points; ip++)
381 {
382 auto& ip_data = _ip_data[ip];
383 auto const& N = Ns[ip];
384 auto& porosity = ip_data.porosity;
385 auto const& porosity_prev = ip_data.porosity_prev;
386 auto const& chemical_system_id = ip_data.chemical_system_id;
387
388 auto const n_component = _transport_process_variables.size();
389 std::vector<double> C_int_pt(n_component);
390 for (unsigned component_id = 0; component_id < n_component;
391 ++component_id)
392 {
393 auto const concentration_index =
395 component_id * concentration_size;
396 auto const local_C =
397 local_x.template segment<concentration_size>(
398 concentration_index);
399
401 C_int_pt[component_id]);
402 }
403
404 {
405 vars_prev.porosity = porosity_prev;
406
407 porosity =
409 ? porosity_prev
410 : medium
411 ->property(
413 .template value<double>(vars, vars_prev, pos, t,
414 dt);
415
416 vars.porosity = porosity;
417 }
418
420 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
421 }
422 }
423
424 void postSpeciationCalculation(std::size_t const ele_id, double const t,
425 double const dt) override
426 {
428 {
429 return;
430 }
431
432 auto const& medium = *_process_data.media_map.getMedium(ele_id);
433
435 pos.setElementID(ele_id);
436
437 for (auto& ip_data : _ip_data)
438 {
439 ip_data.porosity = ip_data.porosity_prev;
440
442 ->updateVolumeFractionPostReaction(ip_data.chemical_system_id,
443 medium, pos,
444 ip_data.porosity, t, dt);
445
447 ip_data.chemical_system_id, medium, ip_data.porosity);
448 }
449 }
450
451 void assemble(double const t, double const dt,
452 std::vector<double> const& local_x,
453 std::vector<double> const& /*local_x_prev*/,
454 std::vector<double>& local_M_data,
455 std::vector<double>& local_K_data,
456 std::vector<double>& local_b_data) override
457 {
458 auto const local_matrix_size = local_x.size();
459 // Nodal DOFs include pressure
460 int const num_nodal_dof = 1 + _transport_process_variables.size();
461 // This assertion is valid only if all nodal d.o.f. use the same shape
462 // matrices.
463 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
464
466 local_M_data, local_matrix_size, local_matrix_size);
468 local_K_data, local_matrix_size, local_matrix_size);
470 local_b_data, local_matrix_size);
471
472 // Get block matrices
473 auto Kpp = local_K.template block<pressure_size, pressure_size>(
475 auto Mpp = local_M.template block<pressure_size, pressure_size>(
477 auto Bp = local_b.template segment<pressure_size>(pressure_index);
478
479 auto local_p = Eigen::Map<const NodalVectorType>(
480 &local_x[pressure_index], pressure_size);
481
482 auto const& b =
485
486 auto const number_of_components = num_nodal_dof - 1;
487 for (int component_id = 0; component_id < number_of_components;
488 ++component_id)
489 {
490 /* Partitioned assembler matrix
491 * | pp | pc1 | pc2 | pc3 |
492 * |-----|-----|-----|-----|
493 * | c1p | c1c1| 0 | 0 |
494 * |-----|-----|-----|-----|
495 * | c2p | 0 | c2c2| 0 |
496 * |-----|-----|-----|-----|
497 * | c3p | 0 | 0 | c3c3|
498 */
499 auto concentration_index =
500 pressure_size + component_id * concentration_size;
501
502 auto KCC =
503 local_K.template block<concentration_size, concentration_size>(
504 concentration_index, concentration_index);
505 auto MCC =
506 local_M.template block<concentration_size, concentration_size>(
507 concentration_index, concentration_index);
508 auto MCp =
509 local_M.template block<concentration_size, pressure_size>(
510 concentration_index, pressure_index);
511 auto MpC =
512 local_M.template block<pressure_size, concentration_size>(
513 pressure_index, concentration_index);
514
515 auto local_C = Eigen::Map<const NodalVectorType>(
516 &local_x[concentration_index], concentration_size);
517
518 assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
519 MCC, MCp, MpC, Kpp, Mpp, Bp);
520
522 {
523 auto const stoichiometric_matrix =
526
527 assert(stoichiometric_matrix);
528
529 for (Eigen::SparseMatrix<double>::InnerIterator it(
530 *stoichiometric_matrix, component_id);
531 it;
532 ++it)
533 {
534 auto const stoichiometric_coefficient = it.value();
535 auto const coupled_component_id = it.row();
536 auto const kinetic_prefactor =
538 ->getKineticPrefactor(coupled_component_id);
539
540 auto const concentration_index =
541 pressure_size + component_id * concentration_size;
542 auto const coupled_concentration_index =
544 coupled_component_id * concentration_size;
545 auto KCmCn = local_K.template block<concentration_size,
547 concentration_index, coupled_concentration_index);
548
549 // account for the coupling between components
550 assembleKCmCn(component_id, t, dt, KCmCn,
551 stoichiometric_coefficient,
552 kinetic_prefactor);
553 }
554 }
555 }
556 }
557
559 GlobalDimVectorType const& b, int const component_id, double const t,
560 double const dt,
561 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
562 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
563 Eigen::Ref<LocalBlockMatrixType> KCC,
564 Eigen::Ref<LocalBlockMatrixType> MCC,
565 Eigen::Ref<LocalBlockMatrixType> MCp,
566 Eigen::Ref<LocalBlockMatrixType> MpC,
567 Eigen::Ref<LocalBlockMatrixType> Kpp,
568 Eigen::Ref<LocalBlockMatrixType> Mpp,
569 Eigen::Ref<LocalSegmentVectorType> Bp)
570 {
571 unsigned const n_integration_points =
573
576
578
579 // Get material properties
580 auto const& medium =
582 // Select the only valid for component transport liquid phase.
583 auto const& phase = medium.phase("AqueousLiquid");
584
585 // Assume that the component name is the same as the process variable
586 // name. Components are shifted by one because the first one is always
587 // pressure.
588 auto const& component = phase.component(
589 _transport_process_variables[component_id].get().getName());
590
591 LocalBlockMatrixType KCC_Laplacian =
592 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
593
594 std::vector<GlobalDimVectorType> ip_flux_vector;
595 double average_velocity_norm = 0.0;
597 {
598 ip_flux_vector.reserve(n_integration_points);
599 }
600
601 auto const& Ns =
603 .NsHigherOrder<typename ShapeFunction::MeshElement>();
604
605 for (unsigned ip(0); ip < n_integration_points; ++ip)
606 {
607 auto& ip_data = _ip_data[ip];
608 auto const& dNdx = ip_data.dNdx;
609 auto const& N = Ns[ip];
610 auto const& w = ip_data.integration_weight;
611 auto& porosity = ip_data.porosity;
612
613 double C_int_pt = 0.0;
614 double p_int_pt = 0.0;
615
616 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
617 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
618
619 vars.concentration = C_int_pt;
620 vars.liquid_phase_pressure = p_int_pt;
621
622 // update according to a particular porosity model
624 .template value<double>(vars, pos, t, dt);
625 vars.porosity = porosity;
626
627 auto const& retardation_factor =
629 .template value<double>(vars, pos, t, dt);
630
631 auto const& solute_dispersivity_transverse = medium.template value<
632 double>(
634
635 auto const& solute_dispersivity_longitudinal =
636 medium.template value<double>(
638 longitudinal_dispersivity);
639
640 // Use the fluid density model to compute the density
641 // TODO (renchao): concentration of which component as the argument
642 // for calculation of fluid density
643 auto const density =
645 .template value<double>(vars, pos, t, dt);
646
647 auto const decay_rate =
649 .template value<double>(vars, pos, t, dt);
650
651 auto const& pore_diffusion_coefficient =
654 .value(vars, pos, t, dt));
655
658 vars, pos, t, dt));
659
660 // Use the viscosity model to compute the viscosity
662 .template value<double>(vars, pos, t, dt);
663
664 GlobalDimMatrixType const K_over_mu = K / mu;
665 GlobalDimVectorType const velocity =
667 ? GlobalDimVectorType(-K_over_mu *
668 (dNdx * p_nodal_values - density * b))
669 : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
670
671 const double drho_dp =
673 .template dValue<double>(
674 vars,
676 pos, t, dt);
677
678 const double drho_dC =
680 .template dValue<double>(
682 t, dt);
683
684 GlobalDimMatrixType const hydrodynamic_dispersion =
687 pore_diffusion_coefficient, velocity, porosity,
688 solute_dispersivity_transverse,
689 solute_dispersivity_longitudinal);
690
691 const double R_times_phi(retardation_factor * porosity);
692 GlobalDimVectorType const mass_density_flow = velocity * density;
693 auto const N_t_N = (N.transpose() * N).eval();
694
696 {
697 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
698 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
699 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
700 }
701 else
702 {
703 ip_flux_vector.emplace_back(mass_density_flow);
704 average_velocity_norm += velocity.norm();
705 }
706 MCC.noalias() += N_t_N * (R_times_phi * density * w);
707 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
708 KCC_Laplacian.noalias() +=
709 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
710
711 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
712
713 // Calculate Mpp, Kpp, and bp in the first loop over components
714 if (component_id == 0)
715 {
716 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
717 Kpp.noalias() +=
718 dNdx.transpose() * K_over_mu * dNdx * (density * w);
719
721 {
722 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
723 (density * density * w);
724 }
725 }
726 }
727
729 {
731 typename ShapeFunction::MeshElement>(
733 _ip_data,
735 ip_flux_vector,
736 average_velocity_norm /
737 static_cast<double>(n_integration_points),
738 KCC_Laplacian);
739 }
740
741 KCC.noalias() += KCC_Laplacian;
742 }
743
744 void assembleKCmCn(int const component_id, double const t, double const dt,
745 Eigen::Ref<LocalBlockMatrixType> KCmCn,
746 double const stoichiometric_coefficient,
747 double const kinetic_prefactor)
748 {
749 unsigned const n_integration_points =
751
754
756
757 auto const& medium =
759 auto const& phase = medium.phase("AqueousLiquid");
760 auto const& component = phase.component(
761 _transport_process_variables[component_id].get().getName());
762
763 auto const& Ns =
765 .NsHigherOrder<typename ShapeFunction::MeshElement>();
766
767 for (unsigned ip(0); ip < n_integration_points; ++ip)
768 {
769 auto& ip_data = _ip_data[ip];
770 auto const& w = ip_data.integration_weight;
771 auto const& N = Ns[ip];
772 auto& porosity = ip_data.porosity;
773
774 auto const retardation_factor =
776 .template value<double>(vars, pos, t, dt);
777
779 .template value<double>(vars, pos, t, dt);
780
781 auto const density =
783 .template value<double>(vars, pos, t, dt);
784
785 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
786 kinetic_prefactor * retardation_factor *
787 porosity * density * N;
788 }
789 }
790
791 void assembleForStaggeredScheme(double const t, double const dt,
792 Eigen::VectorXd const& local_x,
793 Eigen::VectorXd const& local_x_prev,
794 int const process_id,
795 std::vector<double>& local_M_data,
796 std::vector<double>& local_K_data,
797 std::vector<double>& local_b_data) override
798 {
799 if (process_id == _process_data.hydraulic_process_id)
800 {
801 assembleHydraulicEquation(t, dt, local_x, local_x_prev,
802 local_M_data, local_K_data, local_b_data);
803 }
804 else if (process_id == _process_data.thermal_process_id)
805 {
806 assembleHeatTransportEquation(t, dt, local_x, local_x_prev,
807 local_M_data, local_K_data,
808 local_b_data);
809 }
810 else
811 {
812 // Go for assembling in an order of transport process id.
813 assembleComponentTransportEquation(t, dt, local_x, local_x_prev,
814 local_M_data, local_K_data,
815 local_b_data, process_id);
816 }
817 }
818
819 void assembleHydraulicEquation(double const t,
820 double const dt,
821 Eigen::VectorXd const& local_x,
822 Eigen::VectorXd const& local_x_prev,
823 std::vector<double>& local_M_data,
824 std::vector<double>& local_K_data,
825 std::vector<double>& local_b_data)
826 {
827 auto const local_p =
828 local_x.template segment<pressure_size>(pressure_index);
829 auto const local_C = local_x.template segment<concentration_size>(
831 auto const local_C_prev =
832 local_x_prev.segment<concentration_size>(first_concentration_index);
833
834 NodalVectorType local_T = getLocalTemperature(t, local_x);
835
837 local_M_data, pressure_size, pressure_size);
839 local_K_data, pressure_size, pressure_size);
841 local_b_data, pressure_size);
842
843 unsigned const n_integration_points =
845
848
849 auto const& b =
852
853 auto const& medium =
855 auto const& phase = medium.phase("AqueousLiquid");
856
859
860 auto const& Ns =
862 .NsHigherOrder<typename ShapeFunction::MeshElement>();
863
864 for (unsigned ip(0); ip < n_integration_points; ++ip)
865 {
866 auto& ip_data = _ip_data[ip];
867 auto const& dNdx = ip_data.dNdx;
868 auto const& w = ip_data.integration_weight;
869 auto const& N = Ns[ip];
870 auto& porosity = ip_data.porosity;
871 auto const& porosity_prev = ip_data.porosity_prev;
872
873 double const C_int_pt = N.dot(local_C);
874 double const p_int_pt = N.dot(local_p);
875 double const T_int_pt = N.dot(local_T);
876
877 vars.concentration = C_int_pt;
878 vars.liquid_phase_pressure = p_int_pt;
879 vars.temperature = T_int_pt;
880
881 // porosity
882 {
883 vars_prev.porosity = porosity_prev;
884
885 porosity =
887 ? porosity_prev
889 .template value<double>(vars, vars_prev, pos, t,
890 dt);
891
892 vars.porosity = porosity;
893 }
894
895 // Use the fluid density model to compute the density
896 // TODO: Concentration of which component as one of arguments for
897 // calculation of fluid density
898 auto const density =
900 .template value<double>(vars, pos, t, dt);
901
904 vars, pos, t, dt));
905
906 // Use the viscosity model to compute the viscosity
908 .template value<double>(vars, pos, t, dt);
909
910 GlobalDimMatrixType const K_over_mu = K / mu;
911
912 const double drho_dp =
914 .template dValue<double>(
915 vars,
917 pos, t, dt);
918 const double drho_dC =
920 .template dValue<double>(
922 t, dt);
923
924 // matrix assembly
925 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
926 local_K.noalias() +=
927 w * dNdx.transpose() * density * K_over_mu * dNdx;
928
930 {
931 local_b.noalias() +=
932 w * density * density * dNdx.transpose() * K_over_mu * b;
933 }
934
935 // coupling term
936 {
937 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
938
939 local_b.noalias() -=
940 w * N.transpose() * porosity * drho_dC * C_dot;
941 }
942 }
943 }
944
945 void assembleHeatTransportEquation(double const t, double const dt,
946 Eigen::VectorXd const& local_x,
947 Eigen::VectorXd const& /*local_x_prev*/,
948 std::vector<double>& local_M_data,
949 std::vector<double>& local_K_data,
950 std::vector<double>& /*local_b_data*/)
951 {
952 assert(local_x.size() ==
954
955 auto const local_p =
956 local_x.template segment<pressure_size>(pressure_index);
957 auto const local_T =
958 local_x.template segment<temperature_size>(temperature_index);
959
961 local_M_data, temperature_size, temperature_size);
963 local_K_data, temperature_size, temperature_size);
964
966 pos.setElementID(this->_element.getID());
967
968 auto const& process_data = this->_process_data;
969 auto const& medium =
970 *process_data.media_map.getMedium(this->_element.getID());
971 auto const& liquid_phase = medium.phase("AqueousLiquid");
972
973 auto const& b =
976
978
979 unsigned const n_integration_points =
981
982 std::vector<GlobalDimVectorType> ip_flux_vector;
983 double average_velocity_norm = 0.0;
984 ip_flux_vector.reserve(n_integration_points);
985
986 auto const& Ns =
988 .NsHigherOrder<typename ShapeFunction::MeshElement>();
989
990 for (unsigned ip(0); ip < n_integration_points; ip++)
991 {
992 auto const& ip_data = this->_ip_data[ip];
993 auto const& dNdx = ip_data.dNdx;
994 auto const& w = ip_data.integration_weight;
995 auto const& N = Ns[ip];
996
997 double p_at_xi = 0.;
998 NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
999 double T_at_xi = 0.;
1000 NumLib::shapeFunctionInterpolate(local_T, N, T_at_xi);
1001
1002 vars.temperature = T_at_xi;
1003 vars.liquid_phase_pressure = p_at_xi;
1004
1005 vars.liquid_saturation = 1.0;
1006
1007 auto const porosity =
1009 .template value<double>(vars, pos, t, dt);
1010 vars.porosity = porosity;
1011
1012 // Use the fluid density model to compute the density
1013 auto const fluid_density =
1014 liquid_phase
1016 .template value<double>(vars, pos, t, dt);
1017 vars.density = fluid_density;
1018 auto const specific_heat_capacity_fluid =
1019 liquid_phase
1021 .template value<double>(vars, pos, t, dt);
1022
1023 // Assemble mass matrix
1024 local_M.noalias() += w *
1026 vars, porosity, fluid_density,
1027 specific_heat_capacity_fluid, pos, t, dt) *
1028 N.transpose() * N;
1029
1030 // Assemble Laplace matrix
1031 auto const viscosity =
1032 liquid_phase
1034 .template value<double>(vars, pos, t, dt);
1035
1036 auto const intrinsic_permeability =
1038 medium
1039 .property(
1041 .value(vars, pos, t, dt));
1042
1043 GlobalDimMatrixType const K_over_mu =
1044 intrinsic_permeability / viscosity;
1045 GlobalDimVectorType const velocity =
1046 process_data.has_gravity
1047 ? GlobalDimVectorType(-K_over_mu *
1048 (dNdx * local_p - fluid_density * b))
1049 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1050
1051 GlobalDimMatrixType const thermal_conductivity_dispersivity =
1053 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1054 pos, t, dt);
1055
1056 local_K.noalias() +=
1057 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1058
1059 ip_flux_vector.emplace_back(velocity * fluid_density *
1060 specific_heat_capacity_fluid);
1061 average_velocity_norm += velocity.norm();
1062 }
1063
1065 process_data.stabilizer, this->_ip_data,
1066 _process_data.shape_matrix_cache, ip_flux_vector,
1067 average_velocity_norm / static_cast<double>(n_integration_points),
1068 local_K);
1069 }
1070
1072 double const t, double const dt, Eigen::VectorXd const& local_x,
1073 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_M_data,
1074 std::vector<double>& local_K_data,
1075 std::vector<double>& /*local_b_data*/, int const transport_process_id)
1076 {
1077 assert(static_cast<int>(local_x.size()) ==
1080 static_cast<int>(_transport_process_variables.size()) +
1082
1083 auto const local_p =
1084 local_x.template segment<pressure_size>(pressure_index);
1085
1086 NodalVectorType local_T = getLocalTemperature(t, local_x);
1087
1088 auto const local_C = local_x.template segment<concentration_size>(
1090 (transport_process_id - (_process_data.isothermal ? 1 : 2)) *
1092 auto const local_p_prev =
1093 local_x_prev.segment<pressure_size>(pressure_index);
1094
1096 local_M_data, concentration_size, concentration_size);
1098 local_K_data, concentration_size, concentration_size);
1099
1100 LocalBlockMatrixType KCC_Laplacian =
1101 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1102
1103 unsigned const n_integration_points =
1105
1106 std::vector<GlobalDimVectorType> ip_flux_vector;
1107 double average_velocity_norm = 0.0;
1109 {
1110 ip_flux_vector.reserve(n_integration_points);
1111 }
1112
1115
1116 auto const& b =
1119
1122
1123 auto const& medium =
1125 auto const& phase = medium.phase("AqueousLiquid");
1126 auto const component_id =
1127 transport_process_id - (_process_data.isothermal ? 1 : 2);
1128 auto const& component = phase.component(
1129 _transport_process_variables[component_id].get().getName());
1130
1131 auto const& Ns =
1133 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1134
1135 for (unsigned ip(0); ip < n_integration_points; ++ip)
1136 {
1137 auto& ip_data = _ip_data[ip];
1138 auto const& dNdx = ip_data.dNdx;
1139 auto const& w = ip_data.integration_weight;
1140 auto const& N = Ns[ip];
1141 auto& porosity = ip_data.porosity;
1142 auto const& porosity_prev = ip_data.porosity_prev;
1143
1144 double const C_int_pt = N.dot(local_C);
1145 double const p_int_pt = N.dot(local_p);
1146 double const T_int_pt = N.dot(local_T);
1147
1148 vars.concentration = C_int_pt;
1149 vars.liquid_phase_pressure = p_int_pt;
1150 vars.temperature = T_int_pt;
1151
1153 {
1154 vars.temperature = N.dot(local_T);
1155 }
1156
1157 // porosity
1158 {
1159 vars_prev.porosity = porosity_prev;
1160
1161 porosity =
1163 ? porosity_prev
1165 .template value<double>(vars, vars_prev, pos, t,
1166 dt);
1167
1168 vars.porosity = porosity;
1169 }
1170
1171 auto const& retardation_factor =
1173 .template value<double>(vars, pos, t, dt);
1174
1175 auto const& solute_dispersivity_transverse = medium.template value<
1176 double>(
1178 auto const& solute_dispersivity_longitudinal =
1179 medium.template value<double>(
1181 longitudinal_dispersivity);
1182
1183 // Use the fluid density model to compute the density
1184 auto const density =
1186 .template value<double>(vars, pos, t, dt);
1187 auto const decay_rate =
1189 .template value<double>(vars, pos, t, dt);
1190
1191 auto const& pore_diffusion_coefficient =
1194 .value(vars, pos, t, dt));
1195
1198 vars, pos, t, dt));
1199 // Use the viscosity model to compute the viscosity
1201 .template value<double>(vars, pos, t, dt);
1202
1203 GlobalDimMatrixType const K_over_mu = K / mu;
1204 GlobalDimVectorType const velocity =
1206 ? GlobalDimVectorType(-K_over_mu *
1207 (dNdx * local_p - density * b))
1208 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1209
1210 GlobalDimMatrixType const hydrodynamic_dispersion =
1213 pore_diffusion_coefficient, velocity, porosity,
1214 solute_dispersivity_transverse,
1215 solute_dispersivity_longitudinal);
1216
1217 double const R_times_phi = retardation_factor * porosity;
1218 auto const N_t_N = (N.transpose() * N).eval();
1219
1221 {
1222 const double drho_dC =
1224 .template dValue<double>(
1226 pos, t, dt);
1227 local_M.noalias() +=
1228 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1229 }
1230
1231 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1232
1233 // coupling term
1235 {
1236 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1237
1238 const double drho_dp =
1240 .template dValue<double>(vars,
1242 liquid_phase_pressure,
1243 pos, t, dt);
1244
1245 local_K.noalias() +=
1246 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1247 dNdx.transpose() * velocity * N * (density * w);
1248 }
1249 else
1250 {
1251 ip_flux_vector.emplace_back(velocity * density);
1252 average_velocity_norm += velocity.norm();
1253 }
1254 local_K.noalias() +=
1255 N_t_N * (decay_rate * R_times_phi * density * w);
1256
1257 KCC_Laplacian.noalias() += dNdx.transpose() *
1258 hydrodynamic_dispersion * dNdx *
1259 (density * w);
1260 }
1261
1263 {
1265 typename ShapeFunction::MeshElement>(
1267 _process_data.shape_matrix_cache, ip_flux_vector,
1268 average_velocity_norm /
1269 static_cast<double>(n_integration_points),
1270 KCC_Laplacian);
1271 }
1272 local_K.noalias() += KCC_Laplacian;
1273 }
1274
1276 double const t, double const dt, Eigen::VectorXd const& local_x,
1277 Eigen::VectorXd const& local_x_prev, int const process_id,
1278 std::vector<double>& local_b_data,
1279 std::vector<double>& local_Jac_data) override
1280 {
1281 if (process_id == _process_data.hydraulic_process_id)
1282 {
1283 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1284 local_b_data, local_Jac_data);
1285 }
1286 else
1287 {
1288 int const component_id = process_id - 1;
1290 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1291 component_id);
1292 }
1293 }
1294
1296 double const t, double const dt, Eigen::VectorXd const& local_x,
1297 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1298 std::vector<double>& local_Jac_data)
1299 {
1300 auto const p = local_x.template segment<pressure_size>(pressure_index);
1301 auto const c = local_x.template segment<concentration_size>(
1303
1304 auto const p_prev = local_x_prev.segment<pressure_size>(pressure_index);
1305 auto const c_prev =
1306 local_x_prev.segment<concentration_size>(first_concentration_index);
1307
1309 local_Jac_data, pressure_size, pressure_size);
1311 local_b_data, pressure_size);
1312
1313 unsigned const n_integration_points =
1315
1318 auto const& b =
1321
1322 auto const& medium =
1324 auto const& phase = medium.phase("AqueousLiquid");
1325
1328
1329 auto const& Ns =
1331 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1332
1333 for (unsigned ip(0); ip < n_integration_points; ++ip)
1334 {
1335 auto& ip_data = _ip_data[ip];
1336 auto const& dNdx = ip_data.dNdx;
1337 auto const& w = ip_data.integration_weight;
1338 auto const& N = Ns[ip];
1339 auto& phi = ip_data.porosity;
1340 auto const& phi_prev = ip_data.porosity_prev;
1341
1342 double const p_ip = N.dot(p);
1343 double const c_ip = N.dot(c);
1344
1345 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1346
1347 vars.liquid_phase_pressure = p_ip;
1348 vars.concentration = c_ip;
1349
1350 // porosity
1351 {
1352 vars_prev.porosity = phi_prev;
1353
1355 ? phi_prev
1357 .template value<double>(vars, vars_prev, pos, t,
1358 dt);
1359
1360 vars.porosity = phi;
1361 }
1362
1363 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1364 .template value<double>(vars, pos, t, dt);
1365
1368 vars, pos, t, dt));
1369
1371 .template value<double>(vars, pos, t, dt);
1372
1373 auto const drho_dp =
1375 .template dValue<double>(
1376 vars,
1378 pos, t, dt);
1379 auto const drho_dc =
1381 .template dValue<double>(
1383 t, dt);
1384
1385 // matrix assembly
1386 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1387 w * dNdx.transpose() * rho * k / mu * dNdx;
1388
1389 local_rhs.noalias() -=
1390 w * N.transpose() * phi *
1391 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1392 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1393
1395 {
1396 local_rhs.noalias() +=
1397 w * rho * dNdx.transpose() * k / mu * rho * b;
1398 }
1399 }
1400 }
1401
1403 double const t, double const dt, Eigen::VectorXd const& local_x,
1404 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1405 std::vector<double>& local_Jac_data, int const component_id)
1406 {
1407 auto const concentration_index =
1409
1410 auto const p = local_x.template segment<pressure_size>(pressure_index);
1411 auto const c =
1412 local_x.template segment<concentration_size>(concentration_index);
1413 auto const c_prev =
1414 local_x_prev.segment<concentration_size>(concentration_index);
1415
1418 {
1420 }
1421
1423 local_Jac_data, concentration_size, concentration_size);
1425 local_b_data, concentration_size);
1426
1427 LocalBlockMatrixType KCC_Laplacian =
1428 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1429
1430 unsigned const n_integration_points =
1432
1433 std::vector<GlobalDimVectorType> ip_flux_vector;
1434 double average_velocity_norm = 0.0;
1435 ip_flux_vector.reserve(n_integration_points);
1436
1439
1440 auto const& b =
1443
1446
1447 auto const& medium =
1449 auto const& phase = medium.phase("AqueousLiquid");
1450 auto const& component = phase.component(
1451 _transport_process_variables[component_id].get().getName());
1452
1453 auto const& Ns =
1455 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1456
1457 for (unsigned ip(0); ip < n_integration_points; ++ip)
1458 {
1459 auto& ip_data = _ip_data[ip];
1460 auto const& dNdx = ip_data.dNdx;
1461 auto const& w = ip_data.integration_weight;
1462 auto const& N = Ns[ip];
1463 auto& phi = ip_data.porosity;
1464 auto const& phi_prev = ip_data.porosity_prev;
1465
1466 double const p_ip = N.dot(p);
1467 double const c_ip = N.dot(c);
1468
1469 vars.liquid_phase_pressure = p_ip;
1470 vars.concentration = c_ip;
1471
1473 {
1474 vars.temperature = N.dot(T);
1475 }
1476
1477 // porosity
1478 {
1479 vars_prev.porosity = phi_prev;
1480
1482 ? phi_prev
1484 .template value<double>(vars, vars_prev, pos, t,
1485 dt);
1486
1487 vars.porosity = phi;
1488 }
1489
1490 auto const R =
1492 .template value<double>(vars, pos, t, dt);
1493
1494 auto const alpha_T = medium.template value<double>(
1496 auto const alpha_L = medium.template value<double>(
1498
1499 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1500 .template value<double>(vars, pos, t, dt);
1501 // first-order decay constant
1502 auto const alpha =
1504 .template value<double>(vars, pos, t, dt);
1505
1508 .value(vars, pos, t, dt));
1509
1512 vars, pos, t, dt));
1514 .template value<double>(vars, pos, t, dt);
1515 // Darcy flux
1516 GlobalDimVectorType const q =
1518 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1519 : GlobalDimVectorType(-k / mu * dNdx * p);
1520
1522 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1523 alpha_L);
1524
1525 // matrix assembly
1526 local_Jac.noalias() +=
1527 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1528
1529 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1530
1531 auto const cdot = (c - c_prev) / dt;
1532 local_rhs.noalias() -=
1533 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1534
1535 ip_flux_vector.emplace_back(q * rho);
1536 average_velocity_norm += q.norm();
1537 }
1538
1541 _process_data.shape_matrix_cache, ip_flux_vector,
1542 average_velocity_norm / static_cast<double>(n_integration_points),
1543 KCC_Laplacian);
1544
1545 local_rhs.noalias() -= KCC_Laplacian * c;
1546
1547 local_Jac.noalias() += KCC_Laplacian;
1548 }
1549
1551 double const t, double const dt, Eigen::VectorXd const& local_x,
1552 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1553 std::vector<double>& local_b_data,
1554 int const transport_process_id) override
1555 {
1556 auto const local_C = local_x.template segment<concentration_size>(
1558 (transport_process_id - 1) * concentration_size);
1559
1561 local_M_data, concentration_size, concentration_size);
1563 local_K_data, concentration_size, concentration_size);
1565 local_b_data, concentration_size);
1566
1567 unsigned const n_integration_points =
1569
1572
1575
1576 auto const& medium =
1578 auto const component_id = transport_process_id - 1;
1579
1580 auto const& Ns =
1582 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1583
1584 for (unsigned ip(0); ip < n_integration_points; ++ip)
1585 {
1586 auto& ip_data = _ip_data[ip];
1587 auto const w = ip_data.integration_weight;
1588 auto const& N = Ns[ip];
1589 auto& porosity = ip_data.porosity;
1590 auto const& porosity_prev = ip_data.porosity_prev;
1591 auto const chemical_system_id = ip_data.chemical_system_id;
1592
1593 double C_int_pt = 0.0;
1594 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1595
1596 vars.concentration = C_int_pt;
1597
1598 auto const porosity_dot = (porosity - porosity_prev) / dt;
1599
1600 // porosity
1601 {
1602 vars_prev.porosity = porosity_prev;
1603
1604 porosity =
1606 ? porosity_prev
1608 .template value<double>(vars, vars_prev, pos, t,
1609 dt);
1610 }
1611
1612 local_M.noalias() += w * N.transpose() * porosity * N;
1613
1614 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1615
1616 if (chemical_system_id == -1)
1617 {
1618 continue;
1619 }
1620
1621 auto const C_post_int_pt =
1623 component_id, chemical_system_id);
1624
1625 local_b.noalias() +=
1626 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1627 }
1628 }
1629
1630 std::vector<double> const& getIntPtDarcyVelocity(
1631 const double t,
1632 std::vector<GlobalVector*> const& x,
1633 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
1634 std::vector<double>& cache) const override
1635 {
1636 assert(x.size() == dof_table.size());
1637
1638 auto const n_processes = x.size();
1639 std::vector<std::vector<double>> local_x;
1640 local_x.reserve(n_processes);
1641
1642 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1643 {
1644 auto const indices =
1645 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1646 assert(!indices.empty());
1647 local_x.push_back(x[process_id]->get(indices));
1648 }
1649
1650 // only one process, must be monolithic.
1651 if (n_processes == 1)
1652 {
1653 auto const local_p = Eigen::Map<const NodalVectorType>(
1654 &local_x[0][pressure_index], pressure_size);
1655 auto const local_C = Eigen::Map<const NodalVectorType>(
1657 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1658 }
1659
1660 // multiple processes, must be staggered.
1661 {
1662 constexpr int pressure_process_id = 0;
1663 constexpr int concentration_process_id = 1;
1664 auto const local_p = Eigen::Map<const NodalVectorType>(
1665 &local_x[pressure_process_id][0], pressure_size);
1666 auto const local_C = Eigen::Map<const NodalVectorType>(
1667 &local_x[concentration_process_id][0], concentration_size);
1668 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1669 }
1670 }
1671
1672 std::vector<double> const& calculateIntPtDarcyVelocity(
1673 const double t,
1674 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
1675 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
1676 std::vector<double>& cache) const
1677 {
1678 auto const n_integration_points =
1680
1681 cache.clear();
1682 auto cache_mat = MathLib::createZeroedMatrix<
1683 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1684 cache, GlobalDim, n_integration_points);
1685
1688
1689 auto const& b =
1692
1694
1695 auto const& medium =
1697 auto const& phase = medium.phase("AqueousLiquid");
1698
1699 auto const& Ns =
1701 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1702
1703 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1704 {
1705 auto const& ip_data = _ip_data[ip];
1706 auto const& dNdx = ip_data.dNdx;
1707 auto const& N = Ns[ip];
1708 auto const& porosity = ip_data.porosity;
1709
1710 double C_int_pt = 0.0;
1711 double p_int_pt = 0.0;
1712
1713 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1714 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1715
1716 vars.concentration = C_int_pt;
1717 vars.liquid_phase_pressure = p_int_pt;
1718 vars.porosity = porosity;
1719
1720 // TODO (naumov) Temporary value not used by current material
1721 // models. Need extension of secondary variables interface.
1722 double const dt = std::numeric_limits<double>::quiet_NaN();
1725 vars, pos, t, dt));
1727 .template value<double>(vars, pos, t, dt);
1728 GlobalDimMatrixType const K_over_mu = K / mu;
1729
1730 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1732 {
1733 auto const rho_w =
1735 .template value<double>(vars, pos, t, dt);
1736 // here it is assumed that the vector b is directed 'downwards'
1737 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1738 }
1739 }
1740
1741 return cache;
1742 }
1743
1744 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
1745 const unsigned integration_point) const override
1746 {
1748 typename ShapeFunction::MeshElement>()[integration_point];
1749
1750 // assumes N is stored contiguously in memory
1751 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1752 }
1753
1754 Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
1755 double const t,
1756 std::vector<double> const& local_x) const override
1757 {
1758 auto const local_p = Eigen::Map<const NodalVectorType>(
1759 &local_x[pressure_index], pressure_size);
1760 auto const local_C = Eigen::Map<const NodalVectorType>(
1762
1763 // Eval shape matrices at given point
1764 // Note: Axial symmetry is set to false here, because we only need dNdx
1765 // here, which is not affected by axial symmetry.
1766 auto const shape_matrices =
1768 GlobalDim>(
1769 _element, false /*is_axially_symmetric*/,
1770 std::array{pnt_local_coords})[0];
1771
1774 auto const& b =
1777
1779
1780 auto const& medium =
1782 auto const& phase = medium.phase("AqueousLiquid");
1783
1784 // local_x contains the local concentration and pressure values
1785 double c_int_pt;
1786 NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
1787 vars.concentration = c_int_pt;
1788
1789 double p_int_pt;
1790 NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
1791 vars.liquid_phase_pressure = p_int_pt;
1792
1793 // TODO (naumov) Temporary value not used by current material models.
1794 // Need extension of secondary variables interface.
1795 double const dt = std::numeric_limits<double>::quiet_NaN();
1798 vars, pos, t, dt));
1799
1801 .template value<double>(vars, pos, t, dt);
1802 GlobalDimMatrixType const K_over_mu = K / mu;
1803
1804 GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
1805 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1806 .template value<double>(vars, pos, t, dt);
1808 {
1809 q += K_over_mu * rho_w * b;
1810 }
1811 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1812 flux.head<GlobalDim>() = rho_w * q;
1813 return flux;
1814 }
1815
1817 double const t,
1818 double const /*dt*/,
1819 Eigen::VectorXd const& local_x,
1820 Eigen::VectorXd const& /*local_x_prev*/) override
1821 {
1822 auto const local_p =
1823 local_x.template segment<pressure_size>(pressure_index);
1824 auto const local_C = local_x.template segment<concentration_size>(
1826
1827 std::vector<double> ele_velocity;
1828 calculateIntPtDarcyVelocity(t, local_p, local_C, ele_velocity);
1829
1830 auto const n_integration_points =
1832 auto const ele_velocity_mat =
1833 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
1834
1835 auto const ele_id = _element.getID();
1836 Eigen::Map<LocalVectorType>(
1837 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
1838 GlobalDim) =
1839 ele_velocity_mat.rowwise().sum() / n_integration_points;
1840 }
1841
1843 std::size_t const ele_id) override
1844 {
1845 auto const n_integration_points =
1847
1849 {
1850 auto const& medium = *_process_data.media_map.getMedium(ele_id);
1851
1852 for (auto& ip_data : _ip_data)
1853 {
1854 ip_data.porosity = ip_data.porosity_prev;
1855
1857 ->updatePorosityPostReaction(ip_data.chemical_system_id,
1858 medium, ip_data.porosity);
1859 }
1860
1862 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
1863 [](double const s, auto const& ip)
1864 { return s + ip.porosity; }) /
1865 n_integration_points;
1866 }
1867
1868 std::vector<GlobalIndexType> chemical_system_indices;
1869 chemical_system_indices.reserve(n_integration_points);
1870 std::transform(_ip_data.begin(), _ip_data.end(),
1871 std::back_inserter(chemical_system_indices),
1872 [](auto const& ip_data)
1873 { return ip_data.chemical_system_id; });
1874
1876 ele_id, chemical_system_indices);
1877 }
1878
1879 std::vector<double> const& getIntPtMolarFlux(
1880 const double t, std::vector<GlobalVector*> const& x,
1881 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
1882 std::vector<double>& cache, int const component_id) const override
1883 {
1884 std::vector<double> local_x_vec;
1885
1886 auto const n_processes = x.size();
1887 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1888 {
1889 auto const indices =
1890 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
1891 assert(!indices.empty());
1892 auto const local_solution = x[process_id]->get(indices);
1893 local_x_vec.insert(std::end(local_x_vec),
1894 std::begin(local_solution),
1895 std::end(local_solution));
1896 }
1897 auto const local_x = MathLib::toVector(local_x_vec);
1898
1899 auto const p = local_x.template segment<pressure_size>(pressure_index);
1900 auto const c = local_x.template segment<concentration_size>(
1902
1903 auto const n_integration_points =
1905
1906 cache.clear();
1907 auto cache_mat = MathLib::createZeroedMatrix<
1908 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1909 cache, GlobalDim, n_integration_points);
1910
1913
1914 auto const& b =
1917
1919
1920 auto const& medium =
1922 auto const& phase = medium.phase("AqueousLiquid");
1923
1924 auto const& component = phase.component(
1925 _transport_process_variables[component_id].get().getName());
1926
1927 auto const& Ns =
1929 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1930
1931 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1932 {
1933 auto const& ip_data = _ip_data[ip];
1934 auto const& dNdx = ip_data.dNdx;
1935 auto const& N = Ns[ip];
1936 auto const& phi = ip_data.porosity;
1937
1938 double const p_ip = N.dot(p);
1939 double const c_ip = N.dot(c);
1940
1941 vars.concentration = c_ip;
1942 vars.liquid_phase_pressure = p_ip;
1943 vars.porosity = phi;
1944
1945 double const dt = std::numeric_limits<double>::quiet_NaN();
1946
1949 vars, pos, t, dt));
1951 .template value<double>(vars, pos, t, dt);
1952 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1953 .template value<double>(vars, pos, t, dt);
1954
1955 // Darcy flux
1956 GlobalDimVectorType const q =
1958 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1959 : GlobalDimVectorType(-k / mu * dNdx * p);
1960
1961 auto const alpha_T = medium.template value<double>(
1963 auto const alpha_L = medium.template value<double>(
1967 .value(vars, pos, t, dt));
1968
1969 // Hydrodynamic dispersion
1971 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1972 alpha_L);
1973
1974 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
1975 }
1976
1977 return cache;
1978 }
1979
1980 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
1981 Eigen::VectorXd const& /*local_x_prev*/,
1982 double const /*t*/, double const /*dt*/,
1983 int const /*process_id*/) override
1984 {
1985 unsigned const n_integration_points =
1987
1988 for (unsigned ip = 0; ip < n_integration_points; ip++)
1989 {
1990 _ip_data[ip].pushBackState();
1991 }
1992 }
1993
1994private:
1997
1999 std::vector<std::reference_wrapper<ProcessVariable>> const
2001
2002 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>> _ip_data;
2003
2005 MaterialPropertyLib::VariableArray const& vars, const double porosity,
2006 const double fluid_density, const double specific_heat_capacity_fluid,
2007 ParameterLib::SpatialPosition const& pos, double const t,
2008 double const dt)
2009 {
2010 auto const& medium =
2011 *_process_data.media_map.getMedium(this->_element.getID());
2012 auto const& solid_phase = medium.phase("Solid");
2013
2014 auto const specific_heat_capacity_solid =
2015 solid_phase
2016 .property(
2018 .template value<double>(vars, pos, t, dt);
2019
2020 auto const solid_density =
2021 solid_phase.property(MaterialPropertyLib::PropertyType::density)
2022 .template value<double>(vars, pos, t, dt);
2023
2024 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2025 fluid_density * specific_heat_capacity_fluid * porosity;
2026 }
2027
2030 const double fluid_density, const double specific_heat_capacity_fluid,
2031 const GlobalDimVectorType& velocity,
2032 ParameterLib::SpatialPosition const& pos, double const t,
2033 double const dt)
2034 {
2035 auto const& medium =
2037
2038 auto thermal_conductivity =
2040 medium
2041 .property(
2043 .value(vars, pos, t, dt));
2044
2045 auto const thermal_dispersivity_transversal =
2046 medium
2048 thermal_transversal_dispersivity)
2049 .template value<double>();
2050
2051 auto const thermal_dispersivity_longitudinal =
2052 medium
2054 thermal_longitudinal_dispersivity)
2055 .template value<double>();
2056
2057 // Thermal conductivity is moved outside and zero matrix is passed
2058 // instead due to multiplication with fluid's density times specific
2059 // heat capacity.
2060 return thermal_conductivity +
2061 fluid_density * specific_heat_capacity_fluid *
2064 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2065 velocity, 0 /* phi */, thermal_dispersivity_transversal,
2066 thermal_dispersivity_longitudinal);
2067 }
2068
2070 Eigen::VectorXd const& local_x)
2071 {
2072 NodalVectorType local_T;
2074 {
2076 {
2078 _element, t);
2079 }
2080 else
2081 {
2082 local_T = NodalVectorType::Zero(temperature_size);
2083 }
2084 }
2085 else
2086 {
2087 local_T =
2088 local_x.template segment<temperature_size>(temperature_index);
2089 }
2090 return local_T;
2091 }
2092};
2093
2094} // namespace ComponentTransport
2095} // 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:25
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:76
constexpr double getWeight() const
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
auto const & NsHigherOrder() const
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)
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_)