OGS
NonLinearFbar.h
Go to the documentation of this file.
1
12#pragma once
13
14#include <boost/algorithm/string/predicate.hpp>
15#include <cmath>
16#include <string_view>
17#include <tuple>
18
19#include "GMatrix.h"
26
27namespace ProcessLib
28{
29namespace NonLinearFbar
30{
37
39 std::string_view const bar_det_f_type_name);
40
41template <int DisplacementDim, int NPOINTS, typename VectorTypeForFbar,
42 typename GradientVectorType, typename DNDX_Type>
43VectorTypeForFbar computeQVector(DNDX_Type const& dNdx,
44 GradientVectorType const& F,
45 bool const /*is_axially_symmetric*/)
46{
47 // Note: MathLib::VectorizedTensor::toTensor returns a (3 x 3) matrix by
48 // Eigen::Matrix3D. The size of the return matrix is wrong, which
49 // should be (DisplacementDim x DisplacementDim).
50 // Besides, Eigen::Matrix member inverse does not work with the mapped
51 // matrix, e.g. that by using topLeftCorner. Therefore a hard copy is
52 // needed for the following inverse, i.e. declaring F0inv as type
53 // Eigen::MatrixXd instead of type auto.
54 //
55 Eigen::MatrixXd const F_matrix =
57 .template topLeftCorner<DisplacementDim, DisplacementDim>();
58 auto const Finv = F_matrix.inverse();
59
60 VectorTypeForFbar FInvN =
61 VectorTypeForFbar::Zero(DisplacementDim * NPOINTS);
62 for (int i = 0; i < NPOINTS; ++i)
63 {
64 auto const dNidx = dNdx.col(i);
65 for (int k = 0; k < DisplacementDim; k++)
66 {
67 FInvN[k * NPOINTS + i] = Finv.col(k).dot(dNidx);
68 }
69 }
70
71 return FInvN;
72}
73
74template <int DisplacementDim, typename GradientVectorType,
75 typename VectorTypeForFbar, typename NodalVectorType,
76 typename ShapeFunction, typename ShapeMatricesType, typename IpData>
77std::tuple<double, VectorTypeForFbar> computeFBarInitialVariablesAverage(
78 std::vector<IpData, Eigen::aligned_allocator<IpData>> const& ip_data,
79 bool const compute_detF0_only, Eigen::VectorXd const& u,
80 NumLib::GenericIntegrationMethod const& integration_method,
81 MeshLib::Element const& element, bool const is_axially_symmetric)
82{
83 unsigned const n_integration_points =
84 integration_method.getNumberOfPoints();
85 // Compute the element volume
86 double volume = 0.0;
87 for (unsigned ip = 0; ip < n_integration_points; ip++)
88 {
89 auto const& w = ip_data[ip].integration_weight;
90 volume += w;
91 }
92
93 VectorTypeForFbar averaged_grad_N =
94 VectorTypeForFbar::Zero(DisplacementDim * ShapeFunction::NPOINTS);
95 NodalVectorType averaged_N_div_r;
96 if (is_axially_symmetric)
97 {
98 averaged_N_div_r = NodalVectorType::Zero(ShapeFunction::NPOINTS);
99 }
100
101 GradientVectorType F0 =
103
104 for (unsigned i = 0; i < ShapeFunction::NPOINTS; i++)
105 {
106 Eigen::Vector3d const bar_gradN =
107 NumLib::averageGradShapeFunction<DisplacementDim, ShapeFunction,
108 ShapeMatricesType, IpData>(
109 i, element, integration_method, ip_data, is_axially_symmetric) /
110 volume;
111 averaged_grad_N.template segment<DisplacementDim>(i * DisplacementDim) =
112 bar_gradN.template segment<DisplacementDim>(0);
113
114 for (int k = 0; k < DisplacementDim; k++)
115 {
116 F0.template segment<DisplacementDim>(k * DisplacementDim) +=
117 u[k * ShapeFunction::NPOINTS + i] *
118 bar_gradN.template segment<DisplacementDim>(0);
119 }
120 if (is_axially_symmetric)
121 {
122 averaged_N_div_r[i] = bar_gradN[2];
123 F0[4] += bar_gradN[2] * u[i];
124 }
125 }
126
127 if (compute_detF0_only)
128 {
130 }
131
132 VectorTypeForFbar F0InvN =
133 VectorTypeForFbar::Zero(DisplacementDim * ShapeFunction::NPOINTS);
134
135 Eigen::MatrixXd const F_matrix =
137 .template topLeftCorner<DisplacementDim, DisplacementDim>();
138 auto const Finv = F_matrix.inverse();
139
140 for (unsigned i = 0; i < ShapeFunction::NPOINTS; ++i)
141 {
142 auto const dNidx = averaged_grad_N.template segment<DisplacementDim>(
143 i * DisplacementDim);
144 for (int k = 0; k < DisplacementDim; k++)
145 {
146 F0InvN[k * ShapeFunction::NPOINTS + i] = Finv.col(k).dot(dNidx);
147 }
148 // TODO: if(is_axially_symmetric)
149 }
150
151 return {MathLib::VectorizedTensor::determinant(F0), F0InvN};
152}
153
154template <int DisplacementDim, typename GradientVectorType,
155 typename GradientMatrixType, typename VectorTypeForFbar,
156 typename ShapeFunction, typename ShapeMatricesType>
157std::tuple<double, VectorTypeForFbar> computeFBarInitialVariablesElementCenter(
158 bool const compute_detF0_only, Eigen::VectorXd const& u,
159 MeshLib::Element const& element, bool const is_axially_symmetric)
160{
161 auto const shape_matrices_at_element_center =
163 ShapeFunction, ShapeMatricesType, DisplacementDim>(
164 element, is_axially_symmetric);
165
166 auto const N_0 = shape_matrices_at_element_center.N;
167 auto const dNdx_0 = shape_matrices_at_element_center.dNdx;
168
169 // For the 2D case the 33-component is needed (and the four entries
170 // of the non-symmetric matrix); In 3d there are nine entries.
171 GradientMatrixType G0(
172 DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0),
173 ShapeFunction::NPOINTS * DisplacementDim);
174
175 auto const x_coord =
177 element, N_0);
179 dNdx_0, G0, is_axially_symmetric, N_0, x_coord);
180
181 GradientVectorType const grad_u = G0 * u;
182 GradientVectorType const F0 =
184
185 if (compute_detF0_only)
186 {
188 }
189
190 VectorTypeForFbar F0InvN =
191 computeQVector<DisplacementDim, ShapeFunction::NPOINTS,
192 VectorTypeForFbar, GradientVectorType,
193 typename ShapeMatricesType::GlobalDimNodalMatrixType>(
194 dNdx_0, F0, is_axially_symmetric);
195
196 return {MathLib::VectorizedTensor::determinant(F0), F0InvN};
197}
198
199template <int DisplacementDim>
201 bool const is_axially_symmetric)
202{
203 if (DisplacementDim == 3 || is_axially_symmetric)
204 {
207 DisplacementDim)>::identity2;
208 }
209
210 auto identity2 = MathLib::KelvinVector::Invariants<
212 DisplacementDim)>::identity2;
213
214 identity2[2] = 0.0;
215
216 return identity2;
217}
218
219template <int DisplacementDim>
221 double const alpha_p2,
223 bool const is_axially_symmetric)
224{
225 return (2.0 * eps_bar +
226 identityForF<DisplacementDim>(is_axially_symmetric)) /
227 alpha_p2;
228}
229
230} // namespace NonLinearFbar
231} // namespace ProcessLib
Definition of the Element class.
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
double determinant(Eigen::MatrixBase< Derived > const &tensor)
Computes determinant of a vectorized tensor.
Eigen::Matrix3d toTensor(Type< DisplacementDim > const &tensor)
Converts a vectorized tensor to a 3x3 matrix.
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Eigen::Vector3d averageGradShapeFunction(int const local_node_id, MeshLib::Element const &element, NumLib::GenericIntegrationMethod const &integration_method, std::vector< IpData, Eigen::aligned_allocator< IpData > > const &ip_data, const bool is_axially_symmetric)
ShapeMatricesType::ShapeMatrices initShapeMatricesAtElementCenter(MeshLib::Element const &e, bool const is_axially_symmetric)
void computeGMatrix(DNDX_Type const &dNdx, GMatrixType &g_matrix, const bool is_axially_symmetric, N_Type const &N, const double radius)
Fills a G-matrix based on given shape function dN/dx values.
Definition GMatrix.h:25
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > identityForF(bool const is_axially_symmetric)
std::tuple< double, VectorTypeForFbar > computeFBarInitialVariablesElementCenter(bool const compute_detF0_only, Eigen::VectorXd const &u, MeshLib::Element const &element, bool const is_axially_symmetric)
VectorTypeForFbar computeQVector(DNDX_Type const &dNdx, GradientVectorType const &F, bool const)
BarDetFType convertStringToDetFBarType(std::string_view const bar_det_f_type_name)
std::tuple< double, VectorTypeForFbar > computeFBarInitialVariablesAverage(std::vector< IpData, Eigen::aligned_allocator< IpData > > const &ip_data, bool const compute_detF0_only, Eigen::VectorXd const &u, NumLib::GenericIntegrationMethod const &integration_method, MeshLib::Element const &element, bool const is_axially_symmetric)
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > compute2EPlusI(double const alpha_p2, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps_bar, bool const is_axially_symmetric)