diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h index 750f733d9bbb83c0713d5c3e6767299e5df31646..35f1400e6dabe148991750bc70345998bd98c364 100644 --- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h +++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h @@ -85,17 +85,6 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure, ip_data.eps_prev.resize(kelvin_vector_size); ip_data.sigma_eff_prev.resize(kelvin_vector_size); - ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType< - DisplacementDim, displacement_size>::Zero(DisplacementDim, - displacement_size); - for (int i = 0; i < DisplacementDim; ++i) - { - ip_data.N_u_op - .template block<1, displacement_size / DisplacementDim>( - i, i * displacement_size / DisplacementDim) - .noalias() = sm_u.N; - } - ip_data.N_u = sm_u.N; ip_data.dNdx_u = sm_u.dNdx; @@ -208,8 +197,6 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; - auto const& N_u_op = _ip_data[ip].N_u_op; - auto const& N_u = _ip_data[ip].N_u; auto const& dNdx_u = _ip_data[ip].dNdx_u; @@ -314,7 +301,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, double const rho = rho_sr * (1 - porosity) + porosity * rho_fr; local_rhs.template segment<displacement_size>(displacement_index) .noalias() -= - (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w; + (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w; // // displacement equation, pressure part @@ -711,8 +698,6 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; - auto const& N_u_op = _ip_data[ip].N_u_op; - auto const& N_u = _ip_data[ip].N_u; auto const& dNdx_u = _ip_data[ip].dNdx_u; @@ -778,7 +763,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, double const rho = rho_sr * (1 - porosity) + porosity * rho_fr; local_rhs.noalias() -= (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) - - N_u_op.transpose() * rho * b) * + N_u_op(N_u).transpose() * rho * b) * w; } } diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h index f3139b7e0db2a5ba6735a6147e1d5e19ce8a0575..af485aa1033420552076d8488e8bc053677c42dd 100644 --- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h +++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h @@ -17,6 +17,7 @@ #include "LocalAssemblerInterface.h" #include "MaterialLib/PhysicalConstant.h" #include "MaterialLib/SolidModels/LinearElasticIsotropic.h" +#include "MathLib/EigenBlockMatrixView.h" #include "MathLib/KelvinVector.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "NumLib/DOF/DOFTableUtil.h" @@ -48,10 +49,6 @@ struct IntegrationPointData final { } - typename ShapeMatricesTypeDisplacement::template MatrixType< - DisplacementDim, NPoints * DisplacementDim> - N_u_op; /**< for interpolation of the displacement vector, whereas N_u - interpolates one component (scalar) */ typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev; typename BMatricesType::KelvinVectorType eps, eps_prev; @@ -180,6 +177,10 @@ public: using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>; + static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView< + DisplacementDim, + typename ShapeMatricesTypeDisplacement::NodalRowVectorType>; + HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const&) = delete; HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler&&) = delete; diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h index 64971c5bfc50759cc9a1207aa8fa83bdbf4a9eff..86c64530ab6ce0a1857896b849312300cb066b98 100644 --- a/ProcessLib/RichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h @@ -43,9 +43,6 @@ struct IntegrationPointData final sigma_eff_prev.resize(kelvin_vector_size); } - typename ShapeMatrixTypeDisplacement::template MatrixType< - DisplacementDim, NPoints * DisplacementDim> - N_u_op; typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev; typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev; typename BMatricesType::KelvinVectorType eps, eps_prev; diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index e69b8227e66a0a02127da66d3dd5ea500bd3c941..42500cc7132ed44abbec5851bb6b8c3fb5d9218b 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -19,6 +19,7 @@ #include "MaterialLib/MPL/Medium.h" #include "MaterialLib/MPL/Utils/FormEigenTensor.h" #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" +#include "MathLib/EigenBlockMatrixView.h" #include "MathLib/KelvinVector.h" #include "NumLib/Function/Interpolation.h" #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h" @@ -162,17 +163,6 @@ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, _integration_method.getWeightedPoint(ip).getWeight() * sm_u.integralMeasure * sm_u.detJ; - ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType< - DisplacementDim, displacement_size>::Zero(DisplacementDim, - displacement_size); - for (int i = 0; i < DisplacementDim; ++i) - { - ip_data.N_u_op - .template block<1, displacement_size / DisplacementDim>( - i, i * displacement_size / DisplacementDim) - .noalias() = sm_u.N; - } - ip_data.N_u = sm_u.N; ip_data.dNdx_u = sm_u.dNdx; @@ -432,8 +422,6 @@ void RichardsMechanicsLocalAssembler< x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; - auto const& N_u_op = _ip_data[ip].N_u_op; - auto const& N_u = _ip_data[ip].N_u; auto const& dNdx_u = _ip_data[ip].dNdx_u; @@ -645,7 +633,7 @@ void RichardsMechanicsLocalAssembler< // double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR; rhs.template segment<displacement_size>(displacement_index).noalias() -= - (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w; + (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w; // // pressure equation, pressure part. @@ -789,8 +777,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; - auto const& N_u_op = _ip_data[ip].N_u_op; - auto const& N_u = _ip_data[ip].N_u; auto const& dNdx_u = _ip_data[ip].dNdx_u; @@ -984,7 +970,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR; local_rhs.template segment<displacement_size>(displacement_index) .noalias() -= - (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w; + (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w; // // displacement equation, pressure part @@ -1007,7 +993,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, .template block<displacement_size, pressure_size>( displacement_index, pressure_index) .noalias() += - N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w; + N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w; // For the swelling stress with double structure model the corresponding // Jacobian u-p entry would be required, but it does not improve diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h index c458948496de1eb58f4c4d1951a5328e9d432aec..2a28693ed49762590830c52af9b2a429ef4a1dda 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h @@ -17,6 +17,7 @@ #include "LocalAssemblerInterface.h" #include "MaterialLib/MPL/VariableType.h" #include "MaterialLib/SolidModels/LinearElasticIsotropic.h" +#include "MathLib/EigenBlockMatrixView.h" #include "MathLib/KelvinVector.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "NumLib/DOF/DOFTableUtil.h" @@ -72,6 +73,10 @@ public: using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>; + static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView< + DisplacementDim, + typename ShapeMatricesTypeDisplacement::NodalRowVectorType>; + RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler const&) = delete; RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler&&) = delete; diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h index b7ec52b90f60c57ee9c61cda81b6c04cd2a8374d..779ed47fcf582fb5ca41e5edd18b13a6063759ef 100644 --- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -18,6 +18,7 @@ #include "LocalAssemblerInterface.h" #include "MaterialLib/PhysicalConstant.h" #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" +#include "MathLib/EigenBlockMatrixView.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "NumLib/Extrapolation/ExtrapolatableElement.h" #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h" @@ -107,6 +108,9 @@ public: using IpData = IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>; + static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView< + DisplacementDim, typename ShapeMatricesType::NodalRowVectorType>; + SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const&) = delete; SmallDeformationLocalAssembler(SmallDeformationLocalAssembler&&) = delete; @@ -272,20 +276,6 @@ public: auto const& N = _ip_data[ip].N; auto const& dNdx = _ip_data[ip].dNdx; - typename ShapeMatricesType::template MatrixType<DisplacementDim, - displacement_size> - N_u_op = ShapeMatricesType::template MatrixType< - DisplacementDim, - displacement_size>::Zero(DisplacementDim, - displacement_size); - for (int i = 0; i < DisplacementDim; ++i) - { - N_u_op - .template block<1, displacement_size / DisplacementDim>( - i, i * displacement_size / DisplacementDim) - .noalias() = N; - } - auto const x_coord = NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element, N); @@ -336,7 +326,7 @@ public: auto const rho = _process_data.solid_density(t, x_position)[0]; local_b.noalias() -= - (B.transpose() * sigma - N_u_op.transpose() * rho * b) * w; + (B.transpose() * sigma - N_u_op(N).transpose() * rho * b) * w; local_Jac.noalias() += B.transpose() * C * B * w; } } diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index fed694cbf1a138e2097def2ff82509f225c68182..1ceeda5f4ae5a989d95f5c372c383614805251c4 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -50,9 +50,6 @@ struct IntegrationPointData final sigma_eff_prev.resize(kelvin_vector_size); } - typename ShapeMatrixTypeDisplacement::template MatrixType< - DisplacementDim, NPoints * DisplacementDim> - N_u_op; typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev; typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev; typename BMatricesType::KelvinVectorType eps, eps_prev; diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index 58849432443a78bd11d13b3efef6abf483328b86..ec0601f539b9b5ae19cd2cef65dcf26005f45067 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -17,6 +17,7 @@ #include "MaterialLib/MPL/Utils/FormEigenTensor.h" #include "MaterialLib/PhysicalConstant.h" #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h" +#include "MathLib/EigenBlockMatrixView.h" #include "MathLib/KelvinVector.h" #include "NumLib/Function/Interpolation.h" #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h" @@ -74,17 +75,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure, _integration_method.getWeightedPoint(ip).getWeight() * sm_u.integralMeasure * sm_u.detJ; - ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType< - DisplacementDim, displacement_size>::Zero(DisplacementDim, - displacement_size); - for (int i = 0; i < DisplacementDim; ++i) - { - ip_data.N_u_op - .template block<1, displacement_size / DisplacementDim>( - i, i * displacement_size / DisplacementDim) - .noalias() = sm_u.N; - } - ip_data.N_u = sm_u.N; ip_data.dNdx_u = sm_u.dNdx; @@ -1089,7 +1079,6 @@ void TH2MLocalAssembler< auto const& gradNpT = gradNp.transpose().eval(); auto const& gradNTT = gradNT.transpose().eval(); - auto const& Nu_op = ip.N_u_op; auto const& w = ip.integration_weight; auto const& m = Invariants::identity2; @@ -1349,7 +1338,8 @@ void TH2MLocalAssembler< KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_s_L * m * Np) * w; - fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w; + fU.noalias() -= + (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w; if (_process_data.apply_mass_lumping) { @@ -1515,7 +1505,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure, auto const& gradNpT = gradNp.transpose().eval(); auto const& gradNTT = gradNT.transpose().eval(); - auto const& Nu_op = ip.N_u_op; auto const& w = ip.integration_weight; auto const& m = Invariants::identity2; @@ -2047,7 +2036,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure, .noalias() += BuT * ip_cv.C * Bu * w; // fU_1 - fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w; + fU.noalias() -= + (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w; // KuT local_Jac @@ -2059,8 +2049,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure, local_Jac .template block<displacement_size, temperature_size>( displacement_index, temperature_index) - .noalias() += Nu_op.transpose() * ip_cv.drho_dT * b * Nu_op * w; - */ + .noalias() += N_u_op(Nu).transpose() * ip_cv.drho_dT * b * + N_u_op(Nu).transpose() * w; + */ if (_process_data.apply_mass_lumping) { diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h index 09fc3854c29383c87557b2dc7c8435641aa6d1b1..3eb21370d8d44317ee58e95f18c4bf1513ad3b75 100644 --- a/ProcessLib/TH2M/TH2MFEM.h +++ b/ProcessLib/TH2M/TH2MFEM.h @@ -18,6 +18,7 @@ #include "LocalAssemblerInterface.h" #include "MaterialLib/PhysicalConstant.h" #include "MaterialLib/SolidModels/LinearElasticIsotropic.h" +#include "MathLib/EigenBlockMatrixView.h" #include "MathLib/KelvinVector.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "NumLib/DOF/DOFTableUtil.h" @@ -73,6 +74,10 @@ public: MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>; + static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView< + DisplacementDim, + typename ShapeMatricesTypeDisplacement::NodalRowVectorType>; + using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>; TH2MLocalAssembler(TH2MLocalAssembler const&) = delete; diff --git a/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h index e7a31d16bea9540f6bb165498f9e22ff4da4470e..a4a910a70fae1cb8288ed4c8e17978cbc037cd30 100644 --- a/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h +++ b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h @@ -26,10 +26,6 @@ template <typename ShapeMatrixTypeDisplacement, typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints> struct IntegrationPointData final { - typename ShapeMatrixTypeDisplacement::template MatrixType< - DisplacementDim, NPoints * DisplacementDim> - N_u_op; - typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u; typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u; diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h index d925a4bd769a467e49e9b2fa8c5ba611eb9ea4ab..43c1b899258b119817b77761a20c5bcabd93dc63 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h @@ -69,17 +69,6 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction, integration_method.getWeightedPoint(ip).getWeight() * sm_u.integralMeasure * sm_u.detJ; - ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType< - DisplacementDim, displacement_size>::Zero(DisplacementDim, - displacement_size); - for (int i = 0; i < DisplacementDim; ++i) - { - ip_data.N_u_op - .template block<1, displacement_size / DisplacementDim>( - i, i * displacement_size / DisplacementDim) - .noalias() = sm_u.N; - } - ip_data.N_u = sm_u.N; ip_data.dNdx_u = sm_u.dNdx; @@ -316,7 +305,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, typename ConstitutiveTraits::OutputData& output_data) const { auto const& N_u = ip_data.N_u; - auto const& N_u_op = ip_data.N_u_op; auto const& dNdx_u = ip_data.dNdx_u; // N and dNdx are used for both p and T variables @@ -403,7 +391,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, block_u(out.res).noalias() = B.transpose() * sigma_total - static_cast<int>(this->process_data_.apply_body_force_for_deformation) * - N_u_op.transpose() * CD.grav_data.volumetric_body_force; + N_u_op(N_u).transpose() * CD.grav_data.volumetric_body_force; // Storage matrices out.storage_p_a_p.noalias() = CD.eq_p_data.storage_p_a_p_X_NTN * NTN; @@ -447,7 +435,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, B.transpose() * CD.s_mech_data.J_uT_BT_K_N * N; block_up(out.Jac).noalias() = B.transpose() * CD.s_mech_data.J_up_BT_K_N * N + - N_u_op.transpose() * CD.grav_data.J_up_HT_V_N * N; + N_u_op(N_u).transpose() * CD.grav_data.J_up_HT_V_N * N; block_uu(out.Jac).noalias() = B.transpose() * CD.s_mech_data.stiffness_tensor * B; diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h index 45b30aef7df0dc4dc0de90d852a54b21d66422ee..49f4768d7d2d25aefb72da58e9cd928ecd6a7c1e 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h @@ -15,6 +15,7 @@ #include "IntegrationPointData.h" #include "LocalAssemblerInterface.h" +#include "MathLib/EigenBlockMatrixView.h" #include "MathLib/KelvinVector.h" #include "NumLib/Fem/InitShapeMatrices.h" #include "NumLib/Fem/Integration/GenericIntegrationMethod.h" @@ -66,6 +67,10 @@ public: using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>; + static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView< + DisplacementDim, + typename ShapeMatricesTypeDisplacement::NodalRowVectorType>; + ThermoRichardsMechanicsLocalAssembler( ThermoRichardsMechanicsLocalAssembler const&) = delete; ThermoRichardsMechanicsLocalAssembler(