Merge pull request #27371 from dewenyushu/cp_euler_angle

Supply Euler angles from aux-variables for CP

modules/solid_mechanics/include/materials/crystal_plasticity/ComputeElasticityTensorCP.h

@@ -52,4 +52,8 @@ class ComputeElasticityTensorCP : public ComputeElasticityTensor
 
   /// flag for user-defined rotation matrix, supplied in input file
   bool _user_provided_rotation_matrix;
+
+  // The coupled Euler angles component variables
+  unsigned int _n_euler_angle_vars;
+  const std::vector<const VariableValue *> _euler_angle_vars;
 };

modules/solid_mechanics/src/materials/crystal_plasticity/ComputeElasticityTensorCP.C

@@ -21,6 +21,8 @@ ComputeElasticityTensorCP::validParams()
                                   "The ElementReadPropertyFile "
                                   "GeneralUserObject to read element "
                                   "specific property values from file");
+  params.addCoupledVar("euler_angle_variables",
+                       "Vector of coupled variables representing the Euler angles' components.");
   return params;
 }
 
@@ -31,7 +33,9 @@ ComputeElasticityTensorCP::ComputeElasticityTensorCP(const InputParameters & par
                                : nullptr),
     _Euler_angles_mat_prop(declareProperty<RealVectorValue>("Euler_angles")),
     _crysrot(declareProperty<RankTwoTensor>(_base_name + "crysrot")),
-    _R(_Euler_angles)
+    _R(_Euler_angles),
+    _n_euler_angle_vars(coupledComponents("euler_angle_variables")),
+    _euler_angle_vars(coupledValues("euler_angle_variables"))
 {
   // the base class guarantees constant in time, but in this derived class the
   // tensor will rotate over time once plastic deformation sets in
@@ -56,6 +60,15 @@ ComputeElasticityTensorCP::ComputeElasticityTensorCP(const InputParameters & par
        (parameters.isParamSetByUser("euler_angle_3"))))
     mooseError("Bunge Euler angle information and the rotation_matrix cannot both be specified. "
                "Provide only one type of orientation information in the input file.");
+
+  // Check if source of Euler angle values has a conflict
+  if (_read_prop_user_object && _n_euler_angle_vars)
+    paramError("euler_angle_variables",
+               "Euler angles cannot be supplied from both coupled variables and auxiliary input "
+               "file in the option `read_prop_user_object`.");
+
+  if (isCoupled("euler_angle_variables") && _n_euler_angle_vars != 3)
+    paramError("euler_angle_variables", "The Euler angles should have three components.");
 }
 
 void
@@ -67,6 +80,12 @@ ComputeElasticityTensorCP::assignEulerAngles()
     _Euler_angles_mat_prop[_qp](1) = _read_prop_user_object->getData(_current_elem, 1);
     _Euler_angles_mat_prop[_qp](2) = _read_prop_user_object->getData(_current_elem, 2);
   }
+  else if (_n_euler_angle_vars)
+  {
+    _Euler_angles_mat_prop[_qp](0) = (*_euler_angle_vars[0])[_qp];
+    _Euler_angles_mat_prop[_qp](1) = (*_euler_angle_vars[1])[_qp];
+    _Euler_angles_mat_prop[_qp](2) = (*_euler_angle_vars[2])[_qp];
+  }
   else
     _Euler_angles_mat_prop[_qp] = _Euler_angles;
 

modules/solid_mechanics/test/tests/crystal_plasticity/euler_angles/euler_ang_file.txt

@@ -0,0 +1 @@
+../stress_update_material_based/euler_ang_file.txt

modules/solid_mechanics/test/tests/crystal_plasticity/euler_angles/euler_angle_auxvar.i

@@ -0,0 +1,179 @@
+[GlobalParams]
+  displacements = 'disp_x disp_y disp_z'
+[]
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  elem_type = HEX8
+[]
+
+[AuxVariables]
+  [euler_angle_1]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [euler_angle_2]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [euler_angle_3]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  # Euler angles aux variable to check the correctness of value assignments
+  [check_euler_angle_1]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [check_euler_angle_2]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [check_euler_angle_3]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+[]
+
+[Physics]
+  [SolidMechanics]
+    [QuasiStatic]
+      [all]
+        strain = FINITE
+        add_variables = true
+        generate_output = stress_zz
+      []
+    []
+  []
+[]
+
+[AuxKernels]
+  [euler_angle_1]
+    type = FunctionAux
+    variable = euler_angle_1
+    function = '10*t'
+  []
+  [euler_angle_2]
+    type = FunctionAux
+    variable = euler_angle_2
+    function = '20*t'
+  []
+  [euler_angle_3]
+    type = FunctionAux
+    variable = euler_angle_3
+    function = '30*t'
+  []
+  # output Euler angles material property to check correctness of value assignment
+  [mat_euler_angle_1]
+    type = MaterialRealVectorValueAux
+    variable = check_euler_angle_1
+    property = 'Euler_angles'
+    component = 0
+   []
+   [mat_euler_angle_2]
+    type = MaterialRealVectorValueAux
+    variable = check_euler_angle_2
+    property = 'Euler_angles'
+    component = 1
+   []
+   [mat_euler_angle_3]
+    type = MaterialRealVectorValueAux
+    variable = check_euler_angle_3
+    property = 'Euler_angles'
+    component = 2
+   []
+[]
+
+[BCs]
+  [Periodic]
+    [all]
+      variable = 'disp_x'
+      auto_direction = 'z'
+    []
+  []
+  [fix_x]
+    type = DirichletBC
+    variable = disp_x
+    boundary = 'left'
+    value = 0
+  []
+  [fix_y]
+    type = DirichletBC
+    variable = disp_y
+    boundary = 'bottom'
+    value = 0
+  []
+  [fix_z]
+    type = DirichletBC
+    variable = disp_z
+    boundary = 'back'
+    value = 0
+  []
+  [tdisp]
+    type = FunctionDirichletBC
+    variable = disp_z
+    boundary = 'front'
+    function = '0.01*t'
+  []
+[]
+
+[Materials]
+  [elasticity_tensor]
+    type = ComputeElasticityTensorCP
+    C_ijkl = '1.684e5 1.214e5 1.214e5 1.684e5 1.214e5 1.684e5 0.754e5 0.754e5 0.754e5'
+    fill_method = symmetric9
+    euler_angle_variables = 'euler_angle_1 euler_angle_2 euler_angle_3'
+  []
+  [stress]
+    type = ComputeMultipleCrystalPlasticityStress
+    crystal_plasticity_models = 'trial_xtalpl'
+    tan_mod_type = exact
+  []
+  [trial_xtalpl]
+    type = CrystalPlasticityKalidindiUpdate
+    number_slip_systems = 12
+    slip_sys_file_name = input_slip_sys.txt
+  []
+[]
+
+[Postprocessors]
+  [check_euler_angle_1]
+    type = ElementAverageValue
+    variable = check_euler_angle_1
+  []
+  [check_euler_angle_2]
+    type = ElementAverageValue
+    variable = check_euler_angle_2
+  []
+  [check_euler_angle_3]
+    type = ElementAverageValue
+    variable = check_euler_angle_3
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'PJFNK'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = ' lu '
+  nl_abs_tol = 1e-10
+  nl_rel_tol = 1e-10
+  nl_abs_step_tol = 1e-10
+
+  dt = 0.1
+  dtmin = 0.01
+  end_time = 0.5
+[]
+
+[Outputs]
+  csv = true
+[]