Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Pressure.h"
11 :
12 : #include "Assembly.h"
13 : #include "Function.h"
14 : #include "MooseError.h"
15 :
16 : registerMooseObject("TensorMechanicsApp", Pressure);
17 : registerMooseObject("TensorMechanicsApp", ADPressure);
18 :
19 : template <bool is_ad>
20 : InputParameters
21 2874 : PressureTempl<is_ad>::validParams()
22 : {
23 2874 : InputParameters params = PressureParent<is_ad>::validParams();
24 2874 : params.addClassDescription("Applies a pressure on a given boundary in a given direction");
25 5748 : params.addDeprecatedParam<unsigned int>(
26 : "component", "The component for the pressure", "This parameter is no longer necessary");
27 5748 : params.addRequiredCoupledVar("displacements",
28 : "The string of displacements suitable for the problem statement");
29 5748 : params.addDeprecatedParam<Real>("constant",
30 : "The magnitude to use in computing the pressure",
31 : "Use 'factor' in place of 'constant'");
32 5748 : params.addParam<Real>("factor", 1.0, "The magnitude to use in computing the pressure");
33 5748 : params.addParam<FunctionName>("function", "The function that describes the pressure");
34 5748 : params.addParam<PostprocessorName>("postprocessor",
35 : "Postprocessor that will supply the pressure value");
36 5748 : params.addParam<Real>("alpha", 0.0, "alpha parameter required for HHT time integration scheme");
37 2874 : params.set<bool>("use_displaced_mesh") = true;
38 2874 : return params;
39 0 : }
40 :
41 : template <bool is_ad>
42 1437 : PressureTempl<is_ad>::PressureTempl(const InputParameters & parameters)
43 : : PressureParent<is_ad>(parameters),
44 1437 : _component(libMesh::invalid_uint),
45 1437 : _ndisp(this->coupledComponents("displacements")),
46 4264 : _factor(parameters.isParamSetByUser("factor") ? this->template getParam<Real>("factor")
47 1484 : : parameters.isParamSetByUser("constant") ? this->template getParam<Real>("constant")
48 : : 1.0),
49 4094 : _function(this->isParamValid("function") ? &this->getFunction("function") : NULL),
50 1437 : _postprocessor(
51 2874 : this->isParamValid("postprocessor") ? &this->getPostprocessorValue("postprocessor") : NULL),
52 2874 : _alpha(this->template getParam<Real>("alpha")),
53 1437 : _fe_side(_assembly.getFEFace(_var.feType(), _sys.mesh().dimension())),
54 1437 : _q_dxi(nullptr),
55 1437 : _q_deta(nullptr),
56 1437 : _phi_dxi(nullptr),
57 1437 : _phi_deta(nullptr),
58 2874 : _use_displaced_mesh(this->template getParam<bool>("use_displaced_mesh")),
59 4311 : _fe(libMesh::n_threads())
60 : {
61 2827 : if (parameters.isParamSetByUser("factor") && parameters.isParamSetByUser("constant"))
62 0 : mooseError("Error in " + _name + ". Cannot set 'factor' and 'constant'.");
63 :
64 5355 : for (unsigned int i = 0; i < _ndisp; ++i)
65 : {
66 3918 : _disp_var.push_back(this->coupled("displacements", i));
67 3918 : if (_var.number() == _disp_var[i])
68 : {
69 1437 : _component = i;
70 2874 : if (parameters.isParamSetByUser("component") &&
71 1497 : _component != this->template getParam<unsigned int>("component"))
72 0 : mooseError("Incompatibility between component and displacements in " + _name);
73 : }
74 : }
75 1437 : if (_component == libMesh::invalid_uint)
76 0 : mooseError("Problem with displacements in " + _name);
77 1437 : }
78 :
79 : template <bool is_ad>
80 : void
81 1373 : PressureTempl<is_ad>::initialSetup()
82 : {
83 1373 : auto boundary_ids = this->boundaryIDs();
84 : std::set<SubdomainID> block_ids;
85 2857 : for (auto bndry_id : boundary_ids)
86 : {
87 1484 : auto bids = _mesh.getBoundaryConnectedBlocks(bndry_id);
88 1484 : block_ids.insert(bids.begin(), bids.end());
89 : }
90 :
91 1373 : _coord_type = _fe_problem.getCoordSystem(*block_ids.begin());
92 2755 : for (auto blk_id : block_ids)
93 : {
94 1382 : if (_coord_type != _fe_problem.getCoordSystem(blk_id))
95 0 : mooseError("The Pressure BC requires subdomains to have the same coordinate system.");
96 : }
97 1373 : }
98 :
99 : template <bool is_ad>
100 : GenericReal<is_ad>
101 17844854 : PressureTempl<is_ad>::computeQpResidual()
102 : {
103 20179038 : return computeFactor() * (_normals[_qp](_component) * _test[_i][_qp]);
104 : }
105 :
106 : template <bool is_ad>
107 : GenericReal<is_ad>
108 23587970 : PressureTempl<is_ad>::computeFactor() const
109 : {
110 23587970 : GenericReal<is_ad> factor = _factor;
111 :
112 23587970 : if (_function)
113 22480264 : factor *= _function->value(_t + _alpha * _dt, _q_point[_qp]);
114 :
115 23587970 : if (_postprocessor)
116 0 : factor *= *_postprocessor;
117 :
118 23587970 : return factor;
119 : }
120 :
121 : Real
122 5743068 : Pressure::computeFaceStiffness(const unsigned int local_j, const unsigned int coupled_component)
123 : {
124 : //
125 : // Note that this approach will break down for shell elements, i.e.,
126 : // topologically 2D elements in 3D space with pressure loads on
127 : // the faces.
128 : //
129 5743068 : const Real phi_dxi = (*_phi_dxi)[local_j][_qp];
130 5743068 : const Real phi_deta = _phi_deta ? (*_phi_deta)[local_j][_qp] : 0;
131 :
132 5743068 : const RealGradient & dqdxi = (*_q_dxi)[_qp];
133 : const RealGradient out_of_plane(0, 0, 1);
134 5743068 : const RealGradient & dqdeta = _q_deta ? (*_q_deta)[_qp] : out_of_plane;
135 : // Here, b is dqdxi (cross) dqdeta
136 : // Then, normal is b/length(b)
137 5743068 : RealGradient b(dqdxi(1) * dqdeta(2) - dqdxi(2) * dqdeta(1),
138 5743068 : dqdxi(2) * dqdeta(0) - dqdxi(0) * dqdeta(2),
139 5743068 : dqdxi(0) * dqdeta(1) - dqdxi(1) * dqdeta(0));
140 5743068 : const Real inv_length = 1 / (b * _normals[_qp]);
141 :
142 5743068 : const unsigned int i = _component;
143 : const unsigned int j = coupled_component;
144 :
145 : // const int posneg[3][3] = {{0, -1, 1}, {1, 0, -1}, {-1, 1, 0}};
146 5743068 : const int posneg = 1 - (j + (2 - (i + 1) % 3)) % 3;
147 :
148 : // const unsigned int index[3][3] = {{0, 2, 1}, {2, 1, 0}, {1, 0, 2}};
149 5743068 : const unsigned int index = 2 - (j + (i + 2) % 3) % 3;
150 :
151 5743068 : const Real variation_b = posneg * (phi_deta * dqdxi(index) - phi_dxi * dqdeta(index));
152 :
153 : Real rz_term = 0;
154 5743068 : if (_coord_type == Moose::COORD_RZ && j == _subproblem.getAxisymmetricRadialCoord())
155 : {
156 14510 : rz_term = _normals[_qp](i) * _phi[_j][_qp] / _q_point[_qp](0);
157 : }
158 :
159 5743068 : return computeFactor() * _test[_i][_qp] * (inv_length * variation_b + rz_term);
160 : }
161 :
162 : Real
163 24518838 : Pressure::computeStiffness(const unsigned int coupled_component)
164 : {
165 24518838 : if (_ndisp > 1)
166 : {
167 24518790 : const std::map<unsigned int, unsigned int>::iterator j_it = _node_map.find(_j);
168 24518790 : if (_test[_i][_qp] == 0 || j_it == _node_map.end())
169 : return 0;
170 :
171 5743068 : return computeFaceStiffness(j_it->second, coupled_component);
172 : }
173 :
174 48 : else if (_coord_type == Moose::COORD_RSPHERICAL)
175 : {
176 48 : return computeFactor() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] *
177 48 : (2 / _q_point[_qp](0));
178 : }
179 :
180 0 : if (_coord_type == Moose::COORD_RZ)
181 : {
182 0 : return computeFactor() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] /
183 0 : _q_point[_qp](0);
184 : }
185 :
186 : return 0;
187 : }
188 :
189 : Real
190 11662418 : Pressure::computeQpJacobian()
191 : {
192 11662418 : if (_use_displaced_mesh)
193 11385938 : return computeStiffness(_component);
194 :
195 : return 0;
196 : }
197 :
198 : Real
199 13194724 : Pressure::computeQpOffDiagJacobian(const unsigned int jvar_num)
200 : {
201 13194724 : if (_use_displaced_mesh)
202 26335896 : for (unsigned int j = 0; j < _ndisp; ++j)
203 26335512 : if (jvar_num == _disp_var[j])
204 13132900 : return computeStiffness(j);
205 :
206 : return 0;
207 : }
208 :
209 : void
210 345294 : Pressure::precalculateQpJacobian()
211 : {
212 345294 : if (_ndisp == 1)
213 : return;
214 :
215 345282 : if (_fe[_tid] == nullptr)
216 : {
217 763 : const unsigned int dim = _sys.mesh().dimension() - 1;
218 763 : QBase * const & qrule = _assembly.writeableQRuleFace();
219 763 : _fe[_tid] = FEBase::build(dim, _var.feType());
220 763 : _fe[_tid]->attach_quadrature_rule(qrule);
221 : }
222 :
223 345282 : _q_dxi = &_fe[_tid]->get_dxyzdxi();
224 345282 : _phi_dxi = &_fe[_tid]->get_dphidxi();
225 345282 : if (_coord_type == Moose::COORD_XYZ)
226 : {
227 340304 : _q_deta = &_fe[_tid]->get_dxyzdeta();
228 340304 : _phi_deta = &_fe[_tid]->get_dphideta();
229 : }
230 :
231 345282 : _fe[_tid]->reinit(_current_side_elem);
232 :
233 345282 : if (_coord_type == Moose::COORD_XYZ)
234 : {
235 340304 : if (_q_deta->empty())
236 19444 : _q_deta = nullptr;
237 340304 : if (_phi_deta->empty())
238 19444 : _phi_deta = nullptr;
239 : }
240 :
241 : // Compute node map (given elem node, supply face node)
242 : _node_map.clear();
243 345282 : const unsigned int num_node_elem = _current_elem->n_nodes();
244 345282 : const Node * const * elem_nodes = _current_elem->get_nodes();
245 345282 : const unsigned int num_node_face = _current_side_elem->n_nodes();
246 345282 : const Node * const * face_nodes = _current_side_elem->get_nodes();
247 : unsigned int num_found = 0;
248 2478800 : for (unsigned i = 0; i < num_node_elem; ++i)
249 : {
250 9413052 : for (unsigned j = 0; j < num_node_face; ++j)
251 8302080 : if (**(elem_nodes + i) == **(face_nodes + j))
252 : {
253 1367828 : _node_map[i] = j;
254 1367828 : ++num_found;
255 1367828 : break;
256 : }
257 2478800 : if (num_found == num_node_face)
258 : break;
259 : }
260 : }
261 :
262 : void
263 170344 : Pressure::precalculateQpOffDiagJacobian(const MooseVariableFEBase & /*jvar*/)
264 : {
265 170344 : precalculateQpJacobian();
266 170344 : }
267 :
268 : template class PressureTempl<false>;
269 : template class PressureTempl<true>;
|
