18 #include "libmesh/quadrature.h" 19 #include "libmesh/utility.h" 20 #include "libmesh/enum_quadrature_type.h" 21 #include "libmesh/fe_type.h" 22 #include "libmesh/string_to_enum.h" 23 #include "libmesh/quadrature_gauss.h" 33 "rotations",
"The rotations appropriate for the simulation geometry and coordinate system");
36 "The displacements appropriate for the simulation geometry and coordinate system");
39 "Thickness of the shell. Can be supplied as either a number or a variable name.");
41 "Quadrature order in out of plane direction");
43 "large_strain",
false,
"Set to true to turn on finite strain calculations.");
46 "Vector that is projected onto an element to define the " 47 "direction for the first local coordinate direction e1.");
53 _nrot(coupledComponents(
"rotations")),
54 _ndisp(coupledComponents(
"displacements")),
57 _thickness(coupledValue(
"thickness")),
58 _large_strain(getParam<bool>(
"large_strain")),
62 _nonlinear_sys(_fe_problem.getNonlinearSystemBase(0)),
69 _node_normal_old(getMaterialPropertyOldByName<
RealVectorValue>(
"node_normal")),
84 _local_transformation_matrix(),
85 _local_transformation_matrix_old(),
86 _covariant_transformation_matrix(),
87 _covariant_transformation_matrix_old(),
88 _contravariant_transformation_matrix(),
89 _contravariant_transformation_matrix_old(),
90 _total_global_strain(),
91 _sol(_nonlinear_sys.currentSolution()),
92 _sol_old(_nonlinear_sys.solutionOld()),
93 _ref_first_local_dir(getParam<
RealVectorValue>(
"reference_first_local_direction"))
99 "ADComputeIncrementalShellStrain: The number of variables supplied in 'displacements' " 100 "must be 3 and that in 'rotations' must be 2.");
104 "ADComputeIncrementalShellStrain: Shell element is implemented only for linear elements.");
109 "The norm of the defined first local axis is less than the acceptable telerance (1e-3)");
112 for (
unsigned int i = 0; i <
_ndisp; ++i)
124 _t_qrule = std::make_unique<libMesh::QGauss>(
125 1, Utility::string_to_enum<Order>(getParam<std::string>(
"through_thickness_order")));
152 for (
unsigned int i = 0; i <
_t_points.size(); ++i)
155 &declareADProperty<RankTwoTensor>(
"strain_increment_t_points_" + std::to_string(i));
157 &declareADProperty<RankTwoTensor>(
"total_strain_t_points_" + std::to_string(i));
159 &getMaterialPropertyOldByName<RankTwoTensor>(
"total_strain_t_points_" + std::to_string(i));
160 _B[i] = &declareADProperty<DenseMatrix<Real>>(
"B_t_points_" + std::to_string(i));
161 _B_old[i] = &getMaterialPropertyOldByName<DenseMatrix<Real>>(
"B_t_points_" + std::to_string(i));
162 _ge[i] = &declareADProperty<RankTwoTensor>(
"ge_t_points_" + std::to_string(i));
163 _ge_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
"ge_t_points_" + std::to_string(i));
164 _J_map[i] = &declareADProperty<Real>(
"J_mapping_t_points_" + std::to_string(i));
165 _J_map_old[i] = &getMaterialPropertyOldByName<Real>(
"J_mapping_t_points_" + std::to_string(i));
166 _dxyz_dxi[i] = &declareADProperty<RealVectorValue>(
"dxyz_dxi_t_points_" + std::to_string(i));
168 &getMaterialPropertyOldByName<RealVectorValue>(
"dxyz_dxi_t_points_" + std::to_string(i));
169 _dxyz_deta[i] = &declareADProperty<RealVectorValue>(
"dxyz_deta_t_points_" + std::to_string(i));
171 &getMaterialPropertyOldByName<RealVectorValue>(
"dxyz_deta_t_points_" + std::to_string(i));
173 &declareADProperty<RealVectorValue>(
"dxyz_dzeta_t_points_" + std::to_string(i));
175 &getMaterialPropertyOldByName<RealVectorValue>(
"dxyz_dzeta_t_points_" + std::to_string(i));
178 &declareProperty<RankTwoTensor>(
"local_transformation_t_points_" + std::to_string(i));
180 "local_transformation_t_points_" + std::to_string(i));
182 &declareProperty<RankTwoTensor>(
"covariant_transformation_t_points_" + std::to_string(i));
184 "covariant_transformation_t_points_" + std::to_string(i));
186 "contravariant_transformation_t_points_" + std::to_string(i));
188 "contravariant_transformation_t_points_" + std::to_string(i));
190 &declareProperty<RankTwoTensor>(
"total_global_strain_t_points_" + std::to_string(i));
207 "ADComputeIncrementalShellStrain: Shell element is implemented only for 2D elements");
209 mooseError(
"ADComputeIncrementalShellStrain: Shell element needs to have exactly four nodes.");
210 if (
_qrule->get_points().size() != 4)
211 mooseError(
"ADComputeIncrementalShellStrain: Shell element needs to have exactly four " 212 "quadrature points.");
227 FEType fe_type(Utility::string_to_enum<Order>(
"First"),
228 Utility::string_to_enum<FEFamily>(
"LAGRANGE"));
232 _phi_map = fe->get_fe_map().get_phi_map();
234 for (
unsigned int i = 0; i < 4; ++i)
237 for (
unsigned int i = 0; i <
_2d_points.size(); ++i)
258 for (
unsigned int i = 0; i <
_2d_points.size(); ++i)
263 for (
unsigned int temp1 = 0; temp1 < 5; ++temp1)
266 for (
unsigned int temp2 = 0; temp2 < 20; ++temp2)
280 for (
unsigned int ii = 0; ii < 3; ++ii)
281 for (
unsigned int jj = 0; jj < 3; ++jj)
300 FEType fe_type(Utility::string_to_enum<Order>(
"First"),
301 Utility::string_to_enum<FEFamily>(
"LAGRANGE"));
305 _phi_map = fe->get_fe_map().get_phi_map();
307 for (
unsigned int i = 0; i <
_nodes.size(); ++i)
313 normal /= normal.
norm();
315 for (
unsigned int k = 0;
k < 4; ++
k)
323 for (
unsigned int t = 0; t <
_t_points.size(); ++t)
340 for (
unsigned int i = 0; i <
_2d_points.size(); ++i)
351 for (
unsigned int k = 0;
k <
_nodes.size(); ++
k)
366 for (
unsigned int i = 0; i <
_2d_points.size(); ++i)
387 gmn(1, 0) = gmn(0, 1);
388 gmn(2, 0) = gmn(0, 2);
389 gmn(2, 1) = gmn(1, 2);
398 Real normx = std::sqrt(J(0, 0) * J(0, 0) + J(0, 1) * J(0, 1) + J(0, 2) * J(0, 2));
399 Real normy = std::sqrt(J(1, 0) * J(1, 0) + J(1, 1) * J(1, 1) + J(1, 2) * J(1, 2));
400 Real normz = std::sqrt(J(2, 0) * J(2, 0) + J(2, 1) * J(2, 1) + J(2, 2) * J(2, 2));
415 (*_transformation_matrix)[i] = J;
444 (*
_ge[
j])[i](0, 0) = gi0 *
_e1;
445 (*
_ge[
j])[i](0, 1) = gi0 *
_e2;
446 (*
_ge[
j])[i](0, 2) = gi0 *
_e3;
447 (*
_ge[
j])[i](1, 0) = gi1 *
_e1;
448 (*
_ge[
j])[i](1, 1) = gi1 *
_e2;
449 (*
_ge[
j])[i](1, 2) = gi1 *
_e3;
450 (*
_ge[
j])[i](2, 0) = gi2 *
_e1;
451 (*
_ge[
j])[i](2, 1) = gi2 *
_e2;
452 (*
_ge[
j])[i](2, 2) = gi2 *
_e3;
476 mooseError(
"The reference first local axis must not be perpendicular to any of the shell " 492 for (
unsigned int k = 0;
k <
_nodes.size(); ++
k)
501 for (
unsigned int k = 0;
k <
_nodes.size(); ++
k)
509 for (
unsigned int i = 0; i <
_2d_points.size(); ++i)
513 (*
_B[
j])[i].resize(5, 20);
515 for (
unsigned int k = 0;
k <
_nodes.size(); ++
k)
544 (*
_B[
j])[i](2, 12 +
k) =
547 (*
_B[
j])[i](2, 16 +
k) =
612 for (
unsigned int j = 0;
j < 4; ++
j)
617 for (
unsigned int i = 0; i <
_ndisp; ++i)
622 if (ADReal::do_derivatives)
627 for (
unsigned int i = 0; i <
_nrot; ++i)
632 if (ADReal::do_derivatives)
std::vector< ADMaterialProperty< RankTwoTensor > * > _strain_increment
Strain increment in the covariant coordinate system.
NonlinearSystemBase & _nonlinear_sys
Reference to the nonlinear system object.
RankTwoTensorTempl< T > inverse() const
std::vector< const MaterialProperty< DenseMatrix< Real > > * > _B_old
Old B_matrix for small strain.
const MaterialProperty< RealVectorValue > & _node_normal_old
Material property storing the old normal to the element at the 4 nodes.
FEProblemBase & _fe_problem
const QBase *const & _qrule
std::vector< unsigned int > _rot_num
Variable numbers corresponding to the rotational variables.
std::vector< ADRealVectorValue > _v2
First tangential vectors at nodes.
std::vector< std::vector< Real > > _dphidxi_map
Derivatives of shape functions w.r.t isoparametric coordinates xi.
const NumericVector< Number > & _sol_old
ADRealVectorValue _x2
simulation variables
auto norm() const -> decltype(std::norm(Real()))
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
virtual void zero() override final
std::vector< MaterialProperty< RankTwoTensor > * > _contravariant_transformation_matrix
Contravariant base vector matrix material property to transform strain.
std::vector< const MaterialProperty< Real > * > _J_map_old
Old material property containing jacobian of transformation.
std::vector< unsigned int > _disp_num
Variable numbers corresponding to the displacement variables.
std::vector< const MaterialProperty< RealVectorValue > * > _dxyz_deta_old
Old derivative of global x, y and z w.r.t isoparametric coordinate eta.
unsigned int number() const
std::vector< std::vector< Real > > _dphideta_map
Derivatives of shape functions w.r.t isoparametric coordinates eta.
std::vector< ADMaterialProperty< DenseMatrix< Real > > * > _B
B_matrix for small strain.
static const std::string component
registerMooseObject("SolidMechanicsApp", ADComputeIncrementalShellStrain)
std::vector< ADMaterialProperty< RankTwoTensor > * > _total_strain
Total strain increment in the covariant coordinate system.
virtual void initQpStatefulProperties() override
std::vector< ADMaterialProperty< RealVectorValue > * > _dxyz_dzeta
Derivative of global x, y and z w.r.t isoparametric coordinate zeta.
static constexpr std::size_t dim
std::vector< const MaterialProperty< RealVectorValue > * > _dxyz_dxi_old
Old derivative of global x, y and z w.r.t isoparametric coordinate xi.
MooseVariable * getVar(const std::string &var_name, unsigned int comp)
const std::vector< double > y
const VariableValue & _thickness
Coupled variable for the shell thickness.
unsigned int _nrot
Number of coupled rotational variables.
virtual void computeLocalCoordinates()
Computes the element's local coordinates and store in _e1, _e2 and _e3.
std::vector< const MaterialProperty< RankTwoTensor > * > _covariant_transformation_matrix_old
virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) override
virtual void computeBMatrix()
Computes the B matrix that connects strains to nodal displacements and rotations. ...
static RankTwoTensorTempl initializeFromRows(const libMesh::TypeVector< ADReal > &row0, const libMesh::TypeVector< ADReal > &row1, const libMesh::TypeVector< ADReal > &row2)
std::vector< ADMaterialProperty< RealVectorValue > * > _dxyz_deta
Derivative of global x, y and z w.r.t isoparametric coordinate eta.
std::vector< std::vector< unsigned int > > _soln_rot_index
Indices of solution vector corresponding to rotation DOFs in 2 directions at the 4 nodes...
virtual void computeGMatrix()
Computes the transformation matrix from natural coordinates to local cartesian coordinates for elasti...
std::vector< const MaterialProperty< RankTwoTensor > * > _ge_old
Old ge matrix for elasticity tensor conversion.
std::vector< const MaterialProperty< RankTwoTensor > * > _contravariant_transformation_matrix_old
static InputParameters validParams()
const std::vector< double > x
Real f(Real x)
Test function for Brents method.
std::vector< const Node * > _nodes
Vector storing pointers to the nodes of the shell element.
std::vector< Point > _t_points
Quadrature points in the out of plane direction in isoparametric coordinate system.
ADDenseVector _strain_vector
Vector that stores the strain in the the 2 axial and 3 shear directions.
unsigned int _ndisp
Number of coupled displacement variables.
std::vector< ADMaterialProperty< RealVectorValue > * > _dxyz_dxi
Derivative of global x, y and z w.r.t isoparametric coordinate xi.
ADDenseVector _soln_vector
Vector that stores the incremental solution at all the 20 DOFs in the 4 noded element.
ADMaterialProperty< RealVectorValue > & _node_normal
Material property storing the normal to the element at the 4 nodes. Stored as a material property for...
unsigned int number() const
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
virtual void computeProperties() override
std::vector< const MaterialProperty< RankTwoTensor > * > _local_transformation_matrix_old
bool hasSecondOrderElements()
std::unique_ptr< libMesh::QGauss > _t_qrule
Quadrature rule in the out of plane direction.
std::vector< ADRealVectorValue > _v1
First tangential vectors at nodes.
const RealVectorValue _ref_first_local_dir
std::vector< std::vector< unsigned int > > _soln_disp_index
Indices of solution vector corresponding to displacement DOFs in 3 directions at the 4 nodes...
static InputParameters validParams()
std::vector< MaterialProperty< RankTwoTensor > * > _local_transformation_matrix
Transformation matrix to map stresses from global coordinate to the element's local coordinate...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void updateGVectors()
Updates the vectors required for shear locking computation for finite rotations.
std::vector< ADMaterialProperty< RankTwoTensor > * > _ge
ge matrix for elasticity tensor conversion
std::vector< std::vector< Real > > _phi_map
Shape function value.
IntRange< T > make_range(T beg, T end)
std::vector< MaterialProperty< RankTwoTensor > * > _total_global_strain
Total strain in global coordinate system.
void mooseError(Args &&... args) const
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
const NumericVector< Number > *const & _sol
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< MaterialProperty< RankTwoTensor > * > _covariant_transformation_matrix
Covariant base vector matrix material property to transform stress.
std::vector< const MaterialProperty< RealVectorValue > * > _dxyz_dzeta_old
Old derivative of global x, y and z w.r.t isoparametric coordinate zeta.
virtual void computeSolnVector()
Computes the 20x1 soln vector and its derivatives for each shell element.
std::vector< Point > _2d_points
Quadrature points in the in-plane direction in isoparametric coordinate system.
std::vector< ADMaterialProperty< Real > * > _J_map
Material property containing jacobian of transformation.
RankTwoTensor _unrotated_total_strain
static const std::string k
const Elem *const & _current_elem
ADComputeIncrementalShellStrain(const InputParameters ¶meters)
ADMaterialProperty< RankTwoTensor > * _transformation_matrix
Rotation matrix material property.
virtual void computeNodeNormal()
Computes the node normal at each node.
std::vector< const MaterialProperty< RankTwoTensor > * > _total_strain_old
Old total strain increment in the covariant coordinate system.