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 : #include "FEProblemBase.h"
16 :
17 : template <bool is_ad>
18 : InputParameters
19 10082 : PressureBaseTempl<is_ad>::validParams()
20 : {
21 10082 : InputParameters params = PressureBaseParent<is_ad>::validParams();
22 20164 : params.addDeprecatedParam<unsigned int>(
23 : "component", "The component for the pressure", "This parameter is no longer necessary");
24 20164 : params.addParam<bool>("use_displaced_mesh", true, "Whether to use the displaced mesh.");
25 10082 : return params;
26 0 : }
27 :
28 : template <bool is_ad>
29 : InputParameters
30 11104 : PressureBaseTempl<is_ad>::actionParams()
31 : {
32 11104 : auto params = emptyInputParameters();
33 22208 : params.addRequiredCoupledVar("displacements",
34 : "The string of displacements suitable for the problem statement");
35 11104 : return params;
36 0 : }
37 :
38 : template <bool is_ad>
39 3652 : PressureBaseTempl<is_ad>::PressureBaseTempl(const InputParameters & parameters)
40 : : PressureBaseParent<is_ad>(parameters),
41 3652 : _ndisp(this->coupledComponents("displacements")),
42 3652 : _component(
43 10956 : [this, ¶meters]()
44 : {
45 13612 : for (unsigned int i = 0; i < _ndisp; ++i)
46 19920 : _disp_var.push_back(this->coupled("displacements", i));
47 :
48 6825 : for (unsigned int i = 0; i < _ndisp; ++i)
49 : {
50 6825 : if (_var.number() == _disp_var[i])
51 : {
52 7304 : if (parameters.isParamSetByUser("component") &&
53 3859 : i != this->template getParam<unsigned int>("component"))
54 0 : this->paramError(
55 : "component",
56 : "The component this BC is acting on is now inferred from the position "
57 : "of the `variable` in the `displacements` variable list. The explicitly "
58 : "specified component value is at odds with teh automatically inferred "
59 : "value. The `component` parameter has been deprecated. Please double "
60 : "check your input for potential mestakes.");
61 3652 : return i;
62 : }
63 : }
64 :
65 0 : this->paramError("variable", "The BC variable should be a displacement variable.");
66 7304 : }())
67 : {
68 3652 : }
69 :
70 : template <bool is_ad>
71 : void
72 3509 : PressureBaseTempl<is_ad>::initialSetup()
73 : {
74 3509 : auto boundary_ids = this->boundaryIDs();
75 : std::set<SubdomainID> block_ids;
76 7252 : for (auto bndry_id : boundary_ids)
77 : {
78 3743 : auto bids = _mesh.getBoundaryConnectedBlocks(bndry_id);
79 3743 : block_ids.insert(bids.begin(), bids.end());
80 : }
81 :
82 3509 : if (block_ids.size())
83 3483 : _coord_type = _fe_problem.getCoordSystem(*block_ids.begin());
84 : else
85 : {
86 : mooseInfo(
87 : "No connected blocks were found, the coordinate system type is obtained from the mesh.");
88 26 : _coord_type = _fe_problem.mesh().getUniqueCoordSystem();
89 : }
90 7013 : for (auto blk_id : block_ids)
91 : {
92 3504 : if (_coord_type != _fe_problem.getCoordSystem(blk_id))
93 0 : mooseError("The Pressure BC requires subdomains to have the same coordinate system.");
94 : }
95 3509 : }
96 :
97 : template <bool is_ad>
98 : GenericReal<is_ad>
99 50738900 : PressureBaseTempl<is_ad>::computeQpResidual()
100 : {
101 55933004 : return computePressure() * (_normals[_qp](_component) * _test[_i][_qp]);
102 : }
103 :
104 2930 : PressureBase::PressureBase(const InputParameters & parameters)
105 : : PressureBaseTempl<false>(parameters),
106 2930 : _q_dxi(nullptr),
107 2930 : _q_deta(nullptr),
108 2930 : _phi_dxi(nullptr),
109 2930 : _phi_deta(nullptr),
110 2930 : _use_displaced_mesh(this->template getParam<bool>("use_displaced_mesh")),
111 5860 : _fe(libMesh::n_threads())
112 : {
113 2930 : }
114 :
115 : Real
116 17917813 : PressureBase::computeFaceStiffness(const unsigned int local_j, const unsigned int coupled_component)
117 : {
118 : //
119 : // Note that this approach will break down for shell elements, i.e.,
120 : // topologically 2D elements in 3D space with pressure loads on
121 : // the faces.
122 : //
123 17917813 : const Real phi_dxi = (*_phi_dxi)[local_j][_qp];
124 17917813 : const Real phi_deta = _phi_deta ? (*_phi_deta)[local_j][_qp] : 0;
125 :
126 17917813 : const RealGradient & dqdxi = (*_q_dxi)[_qp];
127 : const RealGradient out_of_plane(0, 0, 1);
128 17917813 : const RealGradient & dqdeta = _q_deta ? (*_q_deta)[_qp] : out_of_plane;
129 : // Here, b is dqdxi (cross) dqdeta
130 : // Then, normal is b/length(b)
131 17917813 : RealGradient b(dqdxi(1) * dqdeta(2) - dqdxi(2) * dqdeta(1),
132 17917813 : dqdxi(2) * dqdeta(0) - dqdxi(0) * dqdeta(2),
133 17917813 : dqdxi(0) * dqdeta(1) - dqdxi(1) * dqdeta(0));
134 17917813 : const Real inv_length = 1 / (b * _normals[_qp]);
135 :
136 17917813 : const unsigned int i = _component;
137 : const unsigned int j = coupled_component;
138 :
139 : // const int posneg[3][3] = {{0, -1, 1}, {1, 0, -1}, {-1, 1, 0}};
140 17917813 : const int posneg = 1 - (j + (2 - (i + 1) % 3)) % 3;
141 :
142 : // const unsigned int index[3][3] = {{0, 2, 1}, {2, 1, 0}, {1, 0, 2}};
143 17917813 : const unsigned int index = 2 - (j + (i + 2) % 3) % 3;
144 :
145 17917813 : const Real variation_b = posneg * (phi_deta * dqdxi(index) - phi_dxi * dqdeta(index));
146 :
147 : Real rz_term = 0;
148 17917813 : if (_coord_type == Moose::COORD_RZ && j == _subproblem.getAxisymmetricRadialCoord())
149 : {
150 43814 : rz_term = _normals[_qp](i) * _phi[_j][_qp] / _q_point[_qp](0);
151 : }
152 :
153 17917813 : return computePressure() * _test[_i][_qp] * (inv_length * variation_b + rz_term);
154 : }
155 :
156 : Real
157 75494007 : PressureBase::computeStiffness(const unsigned int coupled_component)
158 : {
159 75494007 : if (_ndisp > 1)
160 : {
161 75493879 : const std::map<unsigned int, unsigned int>::iterator j_it = _node_map.find(_j);
162 75493879 : if (_test[_i][_qp] == 0 || j_it == _node_map.end())
163 : return 0;
164 :
165 17917813 : return computeFaceStiffness(j_it->second, coupled_component);
166 : }
167 :
168 128 : else if (_coord_type == Moose::COORD_RSPHERICAL)
169 : {
170 128 : return computePressure() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] *
171 128 : (2 / _q_point[_qp](0));
172 : }
173 :
174 0 : if (_coord_type == Moose::COORD_RZ)
175 : {
176 0 : return computePressure() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] /
177 0 : _q_point[_qp](0);
178 : }
179 :
180 : return 0;
181 : }
182 :
183 : Real
184 38068945 : PressureBase::computeQpJacobian()
185 : {
186 38068945 : if (_use_displaced_mesh)
187 36546577 : return computeStiffness(_component);
188 :
189 : return 0;
190 : }
191 :
192 : Real
193 40713446 : PressureBase::computeQpOffDiagJacobian(const unsigned int jvar_num)
194 : {
195 40713446 : if (_use_displaced_mesh)
196 78041870 : for (unsigned int j = 0; j < _ndisp; ++j)
197 78040718 : if (jvar_num == _disp_var[j])
198 38947430 : return computeStiffness(j);
199 :
200 : return 0;
201 : }
202 :
203 : void
204 1136799 : PressureBase::precalculateQpJacobian()
205 : {
206 1136799 : if (_ndisp == 1)
207 : return;
208 :
209 1136767 : if (_fe[_tid] == nullptr)
210 : {
211 2117 : const unsigned int dim = _sys.mesh().dimension() - 1;
212 2117 : QBase * const & qrule = _assembly.writeableQRuleFace();
213 2117 : _fe[_tid] = FEBase::build(dim, _var.feType());
214 2117 : _fe[_tid]->attach_quadrature_rule(qrule);
215 : }
216 :
217 1136767 : _q_dxi = &_fe[_tid]->get_dxyzdxi();
218 1136767 : _phi_dxi = &_fe[_tid]->get_dphidxi();
219 1136767 : if (_coord_type == Moose::COORD_XYZ)
220 : {
221 1121889 : _q_deta = &_fe[_tid]->get_dxyzdeta();
222 1121889 : _phi_deta = &_fe[_tid]->get_dphideta();
223 : }
224 :
225 1136767 : _fe[_tid]->reinit(_current_side_elem);
226 :
227 1136767 : if (_coord_type == Moose::COORD_XYZ)
228 : {
229 1121889 : if (_q_deta->empty())
230 57878 : _q_deta = nullptr;
231 1121889 : if (_phi_deta->empty())
232 57878 : _phi_deta = nullptr;
233 : }
234 :
235 : // Compute node map (given elem node, supply face node)
236 : _node_map.clear();
237 1136767 : const unsigned int num_node_elem = _current_elem->n_nodes();
238 1136767 : const Node * const * elem_nodes = _current_elem->get_nodes();
239 1136767 : const unsigned int num_node_face = _current_side_elem->n_nodes();
240 1136767 : const Node * const * face_nodes = _current_side_elem->get_nodes();
241 : unsigned int num_found = 0;
242 8200499 : for (unsigned i = 0; i < num_node_elem; ++i)
243 : {
244 30757137 : for (unsigned j = 0; j < num_node_face; ++j)
245 27044027 : if (**(elem_nodes + i) == **(face_nodes + j))
246 : {
247 4487389 : _node_map[i] = j;
248 4487389 : ++num_found;
249 4487389 : break;
250 : }
251 8200499 : if (num_found == num_node_face)
252 : break;
253 : }
254 : }
255 :
256 : void
257 552110 : PressureBase::precalculateQpOffDiagJacobian(const MooseVariableFEBase & /*jvar*/)
258 : {
259 552110 : precalculateQpJacobian();
260 552110 : }
261 :
262 : template class PressureBaseTempl<false>;
263 : template class PressureBaseTempl<true>;
|
