Merge pull request #29180 from lynnmunday/eigenDecompMatl

Eigen decomp material

modules/solid_mechanics/doc/content/source/materials/EigenDecompositionMaterial.md

@@ -0,0 +1,17 @@
+# EigenDecompositionMaterial
+
+!syntax description /Materials/EigenDecompositionMaterial
+
+## Description
+
+This class reads in a symmetric `RankTwoTensor` material given by [!param](/Materials/EigenDecompositionMaterial/rank_two_tensor) and performs an eigendecomposition on it.  The results of the decomposition are scalar materials named `max_eigen_value`, `mid_eigen_value`, and `min_eigen_value` and the corresponding vector material properties named `max_eigen_vector`, `mid_eigen_vector`, and `min_eigen_vector`. These names can be preceeded by the [!param](/Materials/EigenDecompositionMaterial/base_name) to allow for more than one `EigenDecompositionMaterial` per block.  These material properties can be output using the Material Outputs system with the [!param](/Materials/EigenDecompositionMaterial/output_properties) as shown in the below example.  This material is useful for visualizing the elemental maximum principal stress and direction as a vector.  An error will be produced if an unsymmetric `RankTwoTensor` is decomposed which can happen for the deformation gradient.
+
+## Example Input File Syntax
+
+!listing modules/solid_mechanics/test/tests/eigen_decomp_material/prescribed_strain_3D.i block=Materials/nonADeig_decomp
+
+!syntax parameters /Materials/EigenDecompositionMaterial
+
+!syntax inputs /Materials/EigenDecompositionMaterial
+
+!syntax children /Materials/EigenDecompositionMaterial

modules/solid_mechanics/include/materials/EigenDecompositionMaterial.h

@@ -0,0 +1,44 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "Material.h"
+#include "RankTwoTensorForward.h"
+/**
+ * Perform eigendecomposition on a RankTwoTensor material property
+ */
+template <bool is_ad>
+class EigenDecompositionMaterialTempl : public Material
+{
+public:
+  static InputParameters validParams();
+
+  EigenDecompositionMaterialTempl(const InputParameters & parameters);
+
+protected:
+  virtual void computeQpProperties() override;
+  /// Base name to allow multiple tensors to be decomposed
+  const std::string _base_name;
+  /// Rank two tensor for eigen decomposition
+  const GenericMaterialProperty<RankTwoTensor, is_ad> & _tensor;
+
+  /// eigen vectors
+  GenericMaterialProperty<RealVectorValue, is_ad> & _max_eigen_vector;
+  GenericMaterialProperty<RealVectorValue, is_ad> & _mid_eigen_vector;
+  GenericMaterialProperty<RealVectorValue, is_ad> & _min_eigen_vector;
+
+  /// eigen values
+  GenericMaterialProperty<Real, is_ad> & _max_eigen_value;
+  GenericMaterialProperty<Real, is_ad> & _mid_eigen_value;
+  GenericMaterialProperty<Real, is_ad> & _min_eigen_value;
+};
+
+typedef EigenDecompositionMaterialTempl<false> EigenDecompositionMaterial;
+typedef EigenDecompositionMaterialTempl<true> ADEigenDecompositionMaterial;

modules/solid_mechanics/src/materials/EigenDecompositionMaterial.C

@@ -0,0 +1,72 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "EigenDecompositionMaterial.h"
+#include "RankTwoTensor.h"
+
+registerMooseObject("SolidMechanicsApp", EigenDecompositionMaterial);
+registerMooseObject("SolidMechanicsApp", ADEigenDecompositionMaterial);
+
+template <bool is_ad>
+InputParameters
+EigenDecompositionMaterialTempl<is_ad>::validParams()
+{
+  InputParameters params = Material::validParams();
+  params.addClassDescription("Emits material properties for the eigenvalues and eigenvectors of a "
+                             "symmetric rank two tensor.");
+  params.addRequiredParam<MaterialPropertyName>(
+      "rank_two_tensor",
+      "The name of the symmetric rank two tensor to used in eigen decomposition.");
+  params.addParam<std::string>(
+      "base_name",
+      "Optional parameter to allow multiple tensors to be decomposed on the same block.");
+  return params;
+}
+
+template <bool is_ad>
+EigenDecompositionMaterialTempl<is_ad>::EigenDecompositionMaterialTempl(
+    const InputParameters & parameters)
+  : Material(parameters),
+    _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
+    _tensor(getGenericMaterialProperty<RankTwoTensor, is_ad>("rank_two_tensor")),
+    _max_eigen_vector(
+        declareGenericProperty<RealVectorValue, is_ad>(_base_name + "max_eigen_vector")),
+    _mid_eigen_vector(
+        declareGenericProperty<RealVectorValue, is_ad>(_base_name + "mid_eigen_vector")),
+    _min_eigen_vector(
+        declareGenericProperty<RealVectorValue, is_ad>(_base_name + "min_eigen_vector")),
+    _max_eigen_value(declareGenericProperty<Real, is_ad>(_base_name + "max_eigen_value")),
+    _mid_eigen_value(declareGenericProperty<Real, is_ad>(_base_name + "mid_eigen_value")),
+    _min_eigen_value(declareGenericProperty<Real, is_ad>(_base_name + "min_eigen_value"))
+{
+  if (LIBMESH_DIM != 3)
+    mooseError("EigenDecompositionMaterial is only defined for LIBMESH_DIM=3");
+}
+
+template <bool is_ad>
+void
+EigenDecompositionMaterialTempl<is_ad>::computeQpProperties()
+{
+
+  if (!_tensor[_qp].isSymmetric())
+    mooseError("EigenDecompositionMaterial will only operate on symmetric rank two tensors.");
+
+  std::vector<GenericReal<is_ad>> eigval(3, 0.0);
+  GenericRankTwoTensor<is_ad> eigvec;
+
+  _tensor[_qp].symmetricEigenvaluesEigenvectors(eigval, eigvec);
+
+  _max_eigen_vector[_qp] = eigvec.column(2);
+  _mid_eigen_vector[_qp] = eigvec.column(1);
+  _min_eigen_vector[_qp] = eigvec.column(0);
+
+  _max_eigen_value[_qp] = eigval[2];
+  _mid_eigen_value[_qp] = eigval[1];
+  _min_eigen_value[_qp] = eigval[0];
+}

modules/solid_mechanics/test/tests/eigen_decomp_material/pipe_2d.i

@@ -0,0 +1,125 @@
+
+# pipe with outer diameter = 24", wall thickness 0.979"
+# units in MPa, meters
+[GlobalParams]
+  order = SECOND
+  family = LAGRANGE
+  displacements = 'disp_x disp_y'
+[]
+
+[Mesh]
+  coord_type = RZ
+  [pipe]
+    type = GeneratedMeshGenerator
+    dim = 2
+    xmin = 0.2799334 #12"-0.979"
+    xmax = .3048 #12"
+    ymin = 0
+    ymax = 0.00497 #0.979"/5
+    nx = 5
+    ny = 1
+    elem_type = quad9
+  []
+[]
+
+[Variables]
+  [disp_x]
+  []
+  [disp_y]
+  []
+[]
+
+[Physics]
+  [SolidMechanics]
+    [QuasiStatic]
+      [all]
+        strain = SMALL
+        add_variables = true
+        new_system = true
+        formulation = TOTAL
+        volumetric_locking_correction = false
+      []
+    []
+  []
+[]
+
+[Functions]
+  [inner_pressure]
+    type = ConstantFunction
+    value = 3
+  []
+  [big_pipe]
+    type = ParsedFunction
+    expression = -(p*11.021*11.021)/(12*12-11.021*11.021)
+    symbol_names = 'p'
+    symbol_values = 3
+  []
+[]
+
+[BCs]
+  [fixBottom]
+    type = DirichletBC
+    variable = disp_y
+    boundary = bottom
+    value = 0.0
+    preset = true
+  []
+  [Pressure]
+    [inside]
+      boundary = left
+      function = inner_pressure
+    []
+    [axial]
+      boundary = top
+      function = big_pipe
+    []
+  []
+[]
+
+[Constraints]
+  [top]
+    type = EqualValueBoundaryConstraint
+    variable = disp_y
+    secondary = top
+    penalty = 1e+14
+    formulation = penalty
+  []
+[]
+
+[Materials]
+  [elastic_tensor]
+    type = ComputeIsotropicElasticityTensor
+    youngs_modulus = 168000
+    poissons_ratio = 0.31
+  []
+  [compute_stress]
+    type = ComputeLagrangianLinearElasticStress
+  []
+  [nonADeig_decomp]
+    type = EigenDecompositionMaterial
+    rank_two_tensor = cauchy_stress
+    outputs = exodus
+    output_properties = "max_eigen_vector mid_eigen_vector min_eigen_vector "
+                        "max_eigen_valuemid_eigen_value min_eigen_value"
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = 'NEWTON'
+  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
+  petsc_options_value = 'lu       superlu_dist'
+  nl_rel_tol = 5e-8
+[]
+
+[Postprocessors]
+  [eigval_max]
+    type = ElementAverageMaterialProperty
+    mat_prop = max_eigen_value
+  []
+[]
+
+[Outputs]
+  exodus = on
+  console = true
+[]