25 namespace RichardsMechanics
27 template <
int DisplacementDim>
29 double const t,
SpatialPosition const& x_position,
int const material_id,
34 if (permeability.rows() == DisplacementDim)
38 if (permeability.rows() == 1)
40 return Eigen::Matrix<double, DisplacementDim,
41 DisplacementDim>::Identity() *
46 "Intrinsic permeability dimension is neither %d nor one, but %dx%d.",
47 DisplacementDim, permeability.rows(), permeability.cols());
50 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
51 typename IntegrationMethod,
int DisplacementDim>
53 ShapeFunctionPressure, IntegrationMethod,
55 RichardsMechanicsLocalAssembler(
58 bool const is_axially_symmetric,
59 unsigned const integration_order,
61 : _process_data(process_data),
62 _integration_method(integration_order),
64 _is_axially_symmetric(is_axially_symmetric)
66 unsigned const n_integration_points =
69 _ip_data.reserve(n_integration_points);
72 auto const shape_matrices_u =
75 DisplacementDim>(e, is_axially_symmetric,
78 auto const shape_matrices_p =
80 IntegrationMethod, DisplacementDim>(
83 auto const& solid_material =
88 for (
unsigned ip = 0; ip < n_integration_points; ip++)
90 _ip_data.emplace_back(solid_material);
92 auto const& sm_u = shape_matrices_u[ip];
95 sm_u.integralMeasure * sm_u.detJ;
97 ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
100 for (
int i = 0; i < DisplacementDim; ++i)
102 .template block<1, displacement_size / DisplacementDim>(
106 ip_data.N_u = sm_u.N;
107 ip_data.dNdx_u = sm_u.dNdx;
109 ip_data.N_p = shape_matrices_p[ip].N;
110 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
116 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
117 typename IntegrationMethod,
int DisplacementDim>
119 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
121 std::vector<double>
const& local_x,
122 std::vector<double>& local_M_data,
123 std::vector<double>& local_K_data,
124 std::vector<double>& local_rhs_data)
129 Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
134 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
139 typename ShapeMatricesTypeDisplacement::template MatrixType<
146 typename ShapeMatricesTypeDisplacement::template MatrixType<
153 typename ShapeMatricesTypeDisplacement::template VectorType<
158 auto const material_id =
164 unsigned const n_integration_points =
166 for (
unsigned ip = 0; ip < n_integration_points; ip++)
169 auto const& w =
_ip_data[ip].integration_weight;
171 auto const& N_u_op =
_ip_data[ip].N_u_op;
174 auto const& dNdx_u =
_ip_data[ip].dNdx_u;
177 auto const& dNdx_p =
_ip_data[ip].dNdx_p;
185 ShapeFunctionDisplacement::NPOINTS,
193 auto& S_L =
_ip_data[ip].saturation;
195 auto const alpha =
_process_data.biot_coefficient(t, x_position)[0];
196 auto const rho_SR =
_process_data.solid_density(t, x_position)[0];
197 auto const K_SR =
_process_data.solid_bulk_modulus(t, x_position)[0];
198 auto const K_LR =
_process_data.fluid_bulk_modulus(t, x_position)[0];
201 material_id, t, x_position, -p_cap_ip,
temperature, p_cap_ip);
202 auto const rho_LR =
_process_data.flow_material->getFluidDensity(
207 DisplacementDim>::value>::identity2;
210 material_id, t, x_position, -p_cap_ip,
temperature, p_cap_ip);
212 double const dS_L_dp_cap =
214 material_id, t, x_position, -p_cap_ip,
temperature, S_L);
220 double const mu =
_process_data.flow_material->getFluidViscosity(
224 intrinsicPermeability<DisplacementDim>(
226 (rho_LR * k_rel / mu);
231 eps.noalias() = B *
u;
233 auto C =
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
239 K.template block<displacement_size, displacement_size>(
241 .noalias() += B.transpose() * C * B * w;
245 N_u_op.transpose() * rho * b * w;
250 double const a0 = S_L * (alpha -
porosity) / K_SR;
252 double const specific_storage =
253 dS_L_dp_cap * (p_cap_ip * a0 -
porosity) +
257 .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
261 .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
264 dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
271 .noalias() -= B.transpose() * alpha * S_L * identity2 * N_p * w;
276 M.template block<pressure_size, displacement_size>(
pressure_index,
278 .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
279 identity2.transpose() * B * w;
283 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
284 typename IntegrationMethod,
int DisplacementDim>
286 ShapeFunctionPressure, IntegrationMethod,
289 std::vector<double>
const& local_xdot,
290 const double ,
const double ,
291 std::vector<double>& ,
292 std::vector<double>& ,
293 std::vector<double>& local_rhs_data,
294 std::vector<double>& local_Jac_data)
299 Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
304 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
309 Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
313 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
318 typename ShapeMatricesTypeDisplacement::template MatrixType<
325 typename ShapeMatricesTypeDisplacement::template VectorType<
330 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
334 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
337 typename ShapeMatricesTypeDisplacement::template MatrixType<
339 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
343 typename ShapeMatricesTypeDisplacement::template MatrixType<
345 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
350 auto const material_id =
356 unsigned const n_integration_points =
358 for (
unsigned ip = 0; ip < n_integration_points; ip++)
361 auto const& w =
_ip_data[ip].integration_weight;
363 auto const& N_u_op =
_ip_data[ip].N_u_op;
366 auto const& dNdx_u =
_ip_data[ip].dNdx_u;
369 auto const& dNdx_p =
_ip_data[ip].dNdx_p;
377 ShapeFunctionDisplacement::NPOINTS,
389 for (
int i = 0; i < u_ip.size(); ++i)
392 u.segment(i * ShapeFunctionDisplacement::NPOINTS,
393 ShapeFunctionDisplacement::NPOINTS),
394 N_u, u_ip.coeffRef(i));
398 auto& S_L =
_ip_data[ip].saturation;
399 auto const& sigma_eff =
_ip_data[ip].sigma_eff;
401 auto const alpha =
_process_data.biot_coefficient(t, x_position)[0];
402 auto const rho_SR =
_process_data.solid_density(t, x_position)[0];
403 auto const K_SR =
_process_data.solid_bulk_modulus(t, x_position)[0];
404 auto const K_LR =
_process_data.fluid_bulk_modulus(t, x_position)[0];
408 material_id, t, x_position, -p_cap_ip,
temperature, p_cap_ip);
409 auto const rho_LR =
_process_data.flow_material->getFluidDensity(
414 DisplacementDim>::value>::identity2;
417 material_id, t, x_position, -p_cap_ip,
temperature, p_cap_ip);
419 double const dS_L_dp_cap =
421 material_id, t, x_position, -p_cap_ip,
temperature, S_L);
423 double const d2S_L_dp_cap_2 =
425 material_id, t, x_position, -p_cap_ip,
temperature, S_L);
431 double const mu =
_process_data.flow_material->getFluidViscosity(
435 intrinsicPermeability<DisplacementDim>(
442 eps.noalias() = B *
u;
444 auto C =
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
448 .template block<displacement_size, displacement_size>(
450 .noalias() += B.transpose() * C * B * w;
455 (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
460 Kup.noalias() += B.transpose() * alpha * S_L * identity2 * N_p * w;
473 .template block<displacement_size, pressure_size>(
475 .noalias() -= B.transpose() * alpha *
476 (S_L + p_cap_ip * dS_L_dp_cap) * identity2 * N_p * w;
479 .template block<displacement_size, pressure_size>(
482 N_u_op.transpose() *
porosity * rho_LR * dS_L_dp_cap * b * N_p * w;
486 Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
487 identity2.transpose() * B * w;
492 laplace_p.noalias() +=
493 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
495 double const a0 = (alpha -
porosity) / K_SR;
496 double const specific_storage =
497 dS_L_dp_cap * (p_cap_ip * S_L * a0 -
porosity) +
500 double const dspecific_storage_dp_cap =
501 d2S_L_dp_cap_2 * (p_cap_ip * S_L * a0 -
porosity) +
503 (
porosity / K_LR + a0 * 3 * S_L + dS_L_dp_cap * p_cap_ip * a0);
505 storage_p.noalias() +=
506 N_p.transpose() * rho_LR * specific_storage * N_p * w;
511 .noalias() += N_p.transpose() * rho_LR * p_cap_dot_ip *
512 dspecific_storage_dp_cap * N_p * w;
532 double const dk_rel_dS_l =
533 _process_data.flow_material->getRelativePermeabilityDerivative(
536 grad_p_cap = -dNdx_p * p_L;
540 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
541 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
546 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
547 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
549 local_rhs.template segment<pressure_size>(
pressure_index).noalias() +=
550 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
557 .noalias() += laplace_p + storage_p / dt;
563 .noalias() = Kpu / dt;
566 local_rhs.template segment<pressure_size>(
pressure_index).noalias() -=
567 laplace_p * p_L + storage_p * p_L_dot + Kpu * u_dot;
571 .noalias() += Kup * p_L;
574 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
575 typename IntegrationMethod,
int DisplacementDim>
577 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
580 GlobalVector
const& ,
582 std::vector<double>& cache)
const 584 static const int kelvin_vector_size =
586 auto const num_intpts =
_ip_data.size();
590 double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
591 cache, kelvin_vector_size, num_intpts);
593 for (
unsigned ip = 0; ip < num_intpts; ++ip)
595 auto const& sigma =
_ip_data[ip].sigma_eff;
603 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
604 typename IntegrationMethod,
int DisplacementDim>
606 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
609 GlobalVector
const& ,
611 std::vector<double>& cache)
const 613 auto const kelvin_vector_size =
615 auto const num_intpts =
_ip_data.size();
619 double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
620 cache, kelvin_vector_size, num_intpts);
622 for (
unsigned ip = 0; ip < num_intpts; ++ip)
631 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
632 typename IntegrationMethod,
int DisplacementDim>
634 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
640 std::vector<double>& cache)
const 642 auto const num_intpts =
_ip_data.size();
645 assert(!indices.empty());
646 auto const local_x = current_solution.get(indices);
650 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
651 cache, DisplacementDim, num_intpts);
654 Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
658 auto const material_id =
661 unsigned const n_integration_points =
666 for (
unsigned ip = 0; ip < n_integration_points; ip++)
676 intrinsicPermeability<DisplacementDim>(
680 auto const rho_LR =
_process_data.flow_material->getFluidDensity(
685 auto const& dNdx_p =
_ip_data[ip].dNdx_p;
686 cache_matrix.col(ip).noalias() =
687 -K_over_mu * dNdx_p * p_L - rho_LR * K_over_mu * b;
693 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
694 typename IntegrationMethod,
int DisplacementDim>
696 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
699 GlobalVector
const& ,
701 std::vector<double>& cache)
const 703 auto const num_intpts =
_ip_data.size();
707 Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>>(cache, 1,
710 for (
unsigned ip = 0; ip < num_intpts; ++ip)
712 cache_mat[ip] =
_ip_data[ip].saturation;
718 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
719 typename IntegrationMethod,
int DisplacementDim>
721 ShapeFunctionPressure, IntegrationMethod,
724 const double ,
const std::vector<double>& ,
725 const double ,
const double ,
726 std::vector<double>& ,
727 std::vector<double>& ,
728 std::vector<double>& ,
729 std::vector<double>& ,
732 OGS_FATAL(
"RichardsMechanics; The staggered scheme is not implemented.");
735 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
736 typename IntegrationMethod,
int DisplacementDim>
738 ShapeFunctionPressure, IntegrationMethod,
741 const double ,
const std::vector<double>& ,
742 const double ,
const double ,
743 std::vector<double>& ,
744 std::vector<double>& ,
745 std::vector<double>& ,
746 std::vector<double>& ,
749 OGS_FATAL(
"RichardsMechanics; The staggered scheme is not implemented.");
752 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
753 typename IntegrationMethod,
int DisplacementDim>
755 ShapeFunctionPressure, IntegrationMethod,
759 const std::vector<double>& local_xdot,
760 const double dxdot_dx,
762 std::vector<double>& local_M_data,
763 std::vector<double>& local_K_data,
764 std::vector<double>& local_b_data,
765 std::vector<double>& local_Jac_data,
772 t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
773 local_b_data, local_Jac_data, local_coupled_solutions);
779 t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
780 local_b_data, local_Jac_data, local_coupled_solutions);
783 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
784 typename IntegrationMethod,
int DisplacementDim>
786 ShapeFunctionPressure, IntegrationMethod,
790 bool const use_monolithic_scheme)
792 const int displacement_offset =
796 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
804 for (
int ip = 0; ip < n_integration_points; ip++)
808 auto const& dNdx_u =
_ip_data[ip].dNdx_u;
817 ShapeFunctionDisplacement::NPOINTS,
822 eps.noalias() = B *
u;
824 _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
829 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
830 typename IntegrationMethod,
int DisplacementDim>
832 ShapeFunctionPressure, IntegrationMethod,
835 std::vector<double>
const& local_x)
838 Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
841 auto const material_id =
847 unsigned const n_integration_points =
850 double saturation_avg = 0;
852 for (
unsigned ip = 0; ip < n_integration_points; ip++)
860 auto& S_L =
_ip_data[ip].saturation;
862 material_id, t, x_position, -p_cap_ip,
temperature, p_cap_ip);
863 saturation_avg += S_L;
865 saturation_avg /= n_integration_points;
870 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
void assembleWithJacobianForDeformationEquations(double const t, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, LocalCoupledSolutions const &local_coupled_solutions)
Kelvin vector dimensions for given displacement dimension.
std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
void assembleWithJacobian(double const t, std::vector< double > const &local_x, std::vector< double > const &local_xdot, const double, const double, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
Eigen::Matrix< double, DisplacementDim, DisplacementDim > intrinsicPermeability(double const t, SpatialPosition const &x_position, int const material_id, RichardsFlow::RichardsFlowMaterialProperties const &material)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
void computeSecondaryVariableConcrete(double const t, std::vector< double > const &local_x) override
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
static const int displacement_size
bool const _is_axially_symmetric
void setElementID(std::size_t element_id)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
void postNonLinearSolverConcrete(std::vector< double > const &local_x, double const t, bool const use_monolithic_scheme) override
IntegrationMethod _integration_method
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Eigen::MatrixXd const & getPermeability(const int material_id, const double t, const ProcessLib::SpatialPosition &pos, const int dim) const
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
std::vector< double > const & getIntPtDarcyVelocity(const double t, GlobalVector const ¤t_solution, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< double > &cache) const override
void assemble(double const t, std::vector< double > const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_rhs_data) override
std::vector< double > const & getIntPtSaturation(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &) const override
void assembleWithJacobianForPressureEquations(double const t, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, LocalCoupledSolutions const &local_coupled_solutions)
void assembleWithJacobianForStaggeredScheme(double const t, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, LocalCoupledSolutions const &local_coupled_solutions) override
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
static const int displacement_index
VectorType< GlobalDim > GlobalDimVectorType
std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
static const int pressure_size
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
MeshLib::Element const & _element
static const int pressure_index
virtual std::size_t getID() const final
Returns the ID of the element.
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
#define OGS_FATAL(fmt,...)
RichardsMechanicsProcessData< DisplacementDim > & _process_data
void setIntegrationPoint(unsigned integration_point)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType