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 8658 : PressureBaseTempl<is_ad>::validParams()
20 : {
21 8658 : InputParameters params = PressureBaseParent<is_ad>::validParams();
22 17316 : params.addDeprecatedParam<unsigned int>(
23 : "component", "The component for the pressure", "This parameter is no longer necessary");
24 17316 : params.addParam<bool>("use_displaced_mesh", true, "Whether to use the displaced mesh.");
25 8658 : return params;
26 0 : }
27 :
28 : template <bool is_ad>
29 : InputParameters
30 9530 : PressureBaseTempl<is_ad>::actionParams()
31 : {
32 9530 : auto params = emptyInputParameters();
33 19060 : params.addRequiredCoupledVar("displacements",
34 : "The string of displacements suitable for the problem statement");
35 9530 : return params;
36 0 : }
37 :
38 : template <bool is_ad>
39 3140 : PressureBaseTempl<is_ad>::PressureBaseTempl(const InputParameters & parameters)
40 : : PressureBaseParent<is_ad>(parameters),
41 3140 : _ndisp(this->coupledComponents("displacements")),
42 3140 : _component(
43 18056 : [this, ¶meters]()
44 : {
45 11718 : for (unsigned int i = 0; i < _ndisp; ++i)
46 17156 : _disp_var.push_back(this->coupled("displacements", i));
47 :
48 5874 : for (unsigned int i = 0; i < _ndisp; ++i)
49 : {
50 5874 : if (_var.number() == _disp_var[i])
51 : {
52 6280 : if (parameters.isParamSetByUser("component") &&
53 3314 : 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 3140 : return i;
62 : }
63 : }
64 :
65 0 : this->paramError("variable", "The BC variable should be a displacement variable.");
66 6280 : }())
67 : {
68 3140 : }
69 :
70 : template <bool is_ad>
71 : void
72 3016 : PressureBaseTempl<is_ad>::initialSetup()
73 : {
74 3016 : auto boundary_ids = this->boundaryIDs();
75 : std::set<SubdomainID> block_ids;
76 6254 : for (auto bndry_id : boundary_ids)
77 : {
78 3238 : auto bids = _mesh.getBoundaryConnectedBlocks(bndry_id);
79 3238 : block_ids.insert(bids.begin(), bids.end());
80 : }
81 :
82 3016 : _coord_type = _fe_problem.getCoordSystem(*block_ids.begin());
83 6050 : for (auto blk_id : block_ids)
84 : {
85 3034 : if (_coord_type != _fe_problem.getCoordSystem(blk_id))
86 0 : mooseError("The Pressure BC requires subdomains to have the same coordinate system.");
87 : }
88 3016 : }
89 :
90 : template <bool is_ad>
91 : GenericReal<is_ad>
92 40862256 : PressureBaseTempl<is_ad>::computeQpResidual()
93 : {
94 45295392 : return computePressure() * (_normals[_qp](_component) * _test[_i][_qp]);
95 : }
96 :
97 2492 : PressureBase::PressureBase(const InputParameters & parameters)
98 : : PressureBaseTempl<false>(parameters),
99 2492 : _q_dxi(nullptr),
100 2492 : _q_deta(nullptr),
101 2492 : _phi_dxi(nullptr),
102 2492 : _phi_deta(nullptr),
103 2492 : _use_displaced_mesh(this->template getParam<bool>("use_displaced_mesh")),
104 4984 : _fe(libMesh::n_threads())
105 : {
106 2492 : }
107 :
108 : Real
109 13393656 : PressureBase::computeFaceStiffness(const unsigned int local_j, const unsigned int coupled_component)
110 : {
111 : //
112 : // Note that this approach will break down for shell elements, i.e.,
113 : // topologically 2D elements in 3D space with pressure loads on
114 : // the faces.
115 : //
116 13393656 : const Real phi_dxi = (*_phi_dxi)[local_j][_qp];
117 13393656 : const Real phi_deta = _phi_deta ? (*_phi_deta)[local_j][_qp] : 0;
118 :
119 13393656 : const RealGradient & dqdxi = (*_q_dxi)[_qp];
120 : const RealGradient out_of_plane(0, 0, 1);
121 13393656 : const RealGradient & dqdeta = _q_deta ? (*_q_deta)[_qp] : out_of_plane;
122 : // Here, b is dqdxi (cross) dqdeta
123 : // Then, normal is b/length(b)
124 13393656 : RealGradient b(dqdxi(1) * dqdeta(2) - dqdxi(2) * dqdeta(1),
125 13393656 : dqdxi(2) * dqdeta(0) - dqdxi(0) * dqdeta(2),
126 13393656 : dqdxi(0) * dqdeta(1) - dqdxi(1) * dqdeta(0));
127 13393656 : const Real inv_length = 1 / (b * _normals[_qp]);
128 :
129 13393656 : const unsigned int i = _component;
130 : const unsigned int j = coupled_component;
131 :
132 : // const int posneg[3][3] = {{0, -1, 1}, {1, 0, -1}, {-1, 1, 0}};
133 13393656 : const int posneg = 1 - (j + (2 - (i + 1) % 3)) % 3;
134 :
135 : // const unsigned int index[3][3] = {{0, 2, 1}, {2, 1, 0}, {1, 0, 2}};
136 13393656 : const unsigned int index = 2 - (j + (i + 2) % 3) % 3;
137 :
138 13393656 : const Real variation_b = posneg * (phi_deta * dqdxi(index) - phi_dxi * dqdeta(index));
139 :
140 : Real rz_term = 0;
141 13393656 : if (_coord_type == Moose::COORD_RZ && j == _subproblem.getAxisymmetricRadialCoord())
142 : {
143 31036 : rz_term = _normals[_qp](i) * _phi[_j][_qp] / _q_point[_qp](0);
144 : }
145 :
146 13393656 : return computePressure() * _test[_i][_qp] * (inv_length * variation_b + rz_term);
147 : }
148 :
149 : Real
150 56698284 : PressureBase::computeStiffness(const unsigned int coupled_component)
151 : {
152 56698284 : if (_ndisp > 1)
153 : {
154 56698188 : const std::map<unsigned int, unsigned int>::iterator j_it = _node_map.find(_j);
155 56698188 : if (_test[_i][_qp] == 0 || j_it == _node_map.end())
156 : return 0;
157 :
158 13393656 : return computeFaceStiffness(j_it->second, coupled_component);
159 : }
160 :
161 96 : else if (_coord_type == Moose::COORD_RSPHERICAL)
162 : {
163 96 : return computePressure() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] *
164 96 : (2 / _q_point[_qp](0));
165 : }
166 :
167 0 : if (_coord_type == Moose::COORD_RZ)
168 : {
169 0 : return computePressure() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] /
170 0 : _q_point[_qp](0);
171 : }
172 :
173 : return 0;
174 : }
175 :
176 : Real
177 28483012 : PressureBase::computeQpJacobian()
178 : {
179 28483012 : if (_use_displaced_mesh)
180 27251140 : return computeStiffness(_component);
181 :
182 : return 0;
183 : }
184 :
185 : Real
186 30897896 : PressureBase::computeQpOffDiagJacobian(const unsigned int jvar_num)
187 : {
188 30897896 : if (_use_displaced_mesh)
189 59022816 : for (unsigned int j = 0; j < _ndisp; ++j)
190 59022048 : if (jvar_num == _disp_var[j])
191 29447144 : return computeStiffness(j);
192 :
193 : return 0;
194 : }
195 :
196 : void
197 842604 : PressureBase::precalculateQpJacobian()
198 : {
199 842604 : if (_ndisp == 1)
200 : return;
201 :
202 842580 : if (_fe[_tid] == nullptr)
203 : {
204 1722 : const unsigned int dim = _sys.mesh().dimension() - 1;
205 1722 : QBase * const & qrule = _assembly.writeableQRuleFace();
206 1722 : _fe[_tid] = FEBase::build(dim, _var.feType());
207 1722 : _fe[_tid]->attach_quadrature_rule(qrule);
208 : }
209 :
210 842580 : _q_dxi = &_fe[_tid]->get_dxyzdxi();
211 842580 : _phi_dxi = &_fe[_tid]->get_dphidxi();
212 842580 : if (_coord_type == Moose::COORD_XYZ)
213 : {
214 832048 : _q_deta = &_fe[_tid]->get_dxyzdeta();
215 832048 : _phi_deta = &_fe[_tid]->get_dphideta();
216 : }
217 :
218 842580 : _fe[_tid]->reinit(_current_side_elem);
219 :
220 842580 : if (_coord_type == Moose::COORD_XYZ)
221 : {
222 832048 : if (_q_deta->empty())
223 40392 : _q_deta = nullptr;
224 832048 : if (_phi_deta->empty())
225 40392 : _phi_deta = nullptr;
226 : }
227 :
228 : // Compute node map (given elem node, supply face node)
229 : _node_map.clear();
230 842580 : const unsigned int num_node_elem = _current_elem->n_nodes();
231 842580 : const Node * const * elem_nodes = _current_elem->get_nodes();
232 842580 : const unsigned int num_node_face = _current_side_elem->n_nodes();
233 842580 : const Node * const * face_nodes = _current_side_elem->get_nodes();
234 : unsigned int num_found = 0;
235 6145616 : for (unsigned i = 0; i < num_node_elem; ++i)
236 : {
237 23251944 : for (unsigned j = 0; j < num_node_face; ++j)
238 20446464 : if (**(elem_nodes + i) == **(face_nodes + j))
239 : {
240 3340136 : _node_map[i] = j;
241 3340136 : ++num_found;
242 3340136 : break;
243 : }
244 6145616 : if (num_found == num_node_face)
245 : break;
246 : }
247 : }
248 :
249 : void
250 411056 : PressureBase::precalculateQpOffDiagJacobian(const MooseVariableFEBase & /*jvar*/)
251 : {
252 411056 : precalculateQpJacobian();
253 411056 : }
254 :
255 : template class PressureBaseTempl<false>;
256 : template class PressureBaseTempl<true>;
|
