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 5601 : PressureBaseTempl<is_ad>::validParams()
20 : {
21 5601 : InputParameters params = PressureBaseParent<is_ad>::validParams();
22 11202 : params.addDeprecatedParam<unsigned int>(
23 : "component", "The component for the pressure", "This parameter is no longer necessary");
24 11202 : params.addParam<bool>("use_displaced_mesh", true, "Whether to use the displaced mesh.");
25 5601 : return params;
26 0 : }
27 :
28 : template <bool is_ad>
29 : InputParameters
30 6167 : PressureBaseTempl<is_ad>::actionParams()
31 : {
32 6167 : auto params = emptyInputParameters();
33 12334 : params.addRequiredCoupledVar("displacements",
34 : "The string of displacements suitable for the problem statement");
35 6167 : return params;
36 0 : }
37 :
38 : template <bool is_ad>
39 2031 : PressureBaseTempl<is_ad>::PressureBaseTempl(const InputParameters & parameters)
40 : : PressureBaseParent<is_ad>(parameters),
41 2031 : _ndisp(this->coupledComponents("displacements")),
42 2031 : _component(
43 6093 : [this, ¶meters]()
44 : {
45 7569 : for (unsigned int i = 0; i < _ndisp; ++i)
46 11076 : _disp_var.push_back(this->coupled("displacements", i));
47 :
48 3797 : for (unsigned int i = 0; i < _ndisp; ++i)
49 : {
50 3797 : if (_var.number() == _disp_var[i])
51 : {
52 4062 : if (parameters.isParamSetByUser("component") &&
53 2151 : 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 2031 : return i;
62 : }
63 : }
64 :
65 0 : this->paramError("variable", "The BC variable should be a displacement variable.");
66 4062 : }())
67 : {
68 2031 : }
69 :
70 : template <bool is_ad>
71 : void
72 1951 : PressureBaseTempl<is_ad>::initialSetup()
73 : {
74 1951 : auto boundary_ids = this->boundaryIDs();
75 : std::set<SubdomainID> block_ids;
76 4025 : for (auto bndry_id : boundary_ids)
77 : {
78 2074 : auto bids = _mesh.getBoundaryConnectedBlocks(bndry_id);
79 2074 : block_ids.insert(bids.begin(), bids.end());
80 : }
81 :
82 1951 : if (block_ids.size())
83 1937 : _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 14 : _coord_type = _fe_problem.mesh().getUniqueCoordSystem();
89 : }
90 3900 : for (auto blk_id : block_ids)
91 : {
92 1949 : if (_coord_type != _fe_problem.getCoordSystem(blk_id))
93 0 : mooseError("The Pressure BC requires subdomains to have the same coordinate system.");
94 : }
95 1951 : }
96 :
97 : template <bool is_ad>
98 : GenericReal<is_ad>
99 29100496 : PressureBaseTempl<is_ad>::computeQpResidual()
100 : {
101 32074960 : return computePressure() * (_normals[_qp](_component) * _test[_i][_qp]);
102 : }
103 :
104 1634 : PressureBase::PressureBase(const InputParameters & parameters)
105 : : PressureBaseTempl<false>(parameters),
106 1634 : _q_dxi(nullptr),
107 1634 : _q_deta(nullptr),
108 1634 : _phi_dxi(nullptr),
109 1634 : _phi_deta(nullptr),
110 1634 : _use_displaced_mesh(this->template getParam<bool>("use_displaced_mesh")),
111 3268 : _fe(libMesh::n_threads())
112 : {
113 1634 : }
114 :
115 : Real
116 10705730 : 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 10705730 : const Real phi_dxi = (*_phi_dxi)[local_j][_qp];
124 10705730 : const Real phi_deta = _phi_deta ? (*_phi_deta)[local_j][_qp] : 0;
125 :
126 10705730 : const RealGradient & dqdxi = (*_q_dxi)[_qp];
127 : const RealGradient out_of_plane(0, 0, 1);
128 10705730 : 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 10705730 : RealGradient b(dqdxi(1) * dqdeta(2) - dqdxi(2) * dqdeta(1),
132 10705730 : dqdxi(2) * dqdeta(0) - dqdxi(0) * dqdeta(2),
133 10705730 : dqdxi(0) * dqdeta(1) - dqdxi(1) * dqdeta(0));
134 10705730 : const Real inv_length = 1 / (b * _normals[_qp]);
135 :
136 10705730 : 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 10705730 : 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 10705730 : const unsigned int index = 2 - (j + (i + 2) % 3) % 3;
144 :
145 10705730 : const Real variation_b = posneg * (phi_deta * dqdxi(index) - phi_dxi * dqdeta(index));
146 :
147 : Real rz_term = 0;
148 10705730 : if (_coord_type == Moose::COORD_RZ && j == _subproblem.getAxisymmetricRadialCoord())
149 : {
150 28296 : rz_term = _normals[_qp](i) * _phi[_j][_qp] / _q_point[_qp](0);
151 : }
152 :
153 10705730 : return computePressure() * _test[_i][_qp] * (inv_length * variation_b + rz_term);
154 : }
155 :
156 : Real
157 44513038 : PressureBase::computeStiffness(const unsigned int coupled_component)
158 : {
159 44513038 : if (_ndisp > 1)
160 : {
161 44512958 : const std::map<unsigned int, unsigned int>::iterator j_it = _node_map.find(_j);
162 44512958 : if (_test[_i][_qp] == 0 || j_it == _node_map.end())
163 : return 0;
164 :
165 10705730 : return computeFaceStiffness(j_it->second, coupled_component);
166 : }
167 :
168 80 : else if (_coord_type == Moose::COORD_RSPHERICAL)
169 : {
170 80 : return computePressure() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] *
171 80 : (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 22408094 : PressureBase::computeQpJacobian()
185 : {
186 22408094 : if (_use_displaced_mesh)
187 21501566 : return computeStiffness(_component);
188 :
189 : return 0;
190 : }
191 :
192 : Real
193 24052112 : PressureBase::computeQpOffDiagJacobian(const unsigned int jvar_num)
194 : {
195 24052112 : if (_use_displaced_mesh)
196 46105946 : for (unsigned int j = 0; j < _ndisp; ++j)
197 46105178 : if (jvar_num == _disp_var[j])
198 23011472 : return computeStiffness(j);
199 :
200 : return 0;
201 : }
202 :
203 : void
204 697084 : PressureBase::precalculateQpJacobian()
205 : {
206 697084 : if (_ndisp == 1)
207 : return;
208 :
209 697064 : if (_fe[_tid] == nullptr)
210 : {
211 1192 : const unsigned int dim = _sys.mesh().dimension() - 1;
212 1192 : QBase * const & qrule = _assembly.writeableQRuleFace();
213 1192 : _fe[_tid] = FEBase::build(dim, _var.feType());
214 1192 : _fe[_tid]->attach_quadrature_rule(qrule);
215 : }
216 :
217 697064 : _q_dxi = &_fe[_tid]->get_dxyzdxi();
218 697064 : _phi_dxi = &_fe[_tid]->get_dphidxi();
219 697064 : if (_coord_type == Moose::COORD_XYZ)
220 : {
221 687452 : _q_deta = &_fe[_tid]->get_dxyzdeta();
222 687452 : _phi_deta = &_fe[_tid]->get_dphideta();
223 : }
224 :
225 697064 : _fe[_tid]->reinit(_current_side_elem);
226 :
227 697064 : if (_coord_type == Moose::COORD_XYZ)
228 : {
229 687452 : if (_q_deta->empty())
230 37660 : _q_deta = nullptr;
231 687452 : if (_phi_deta->empty())
232 37660 : _phi_deta = nullptr;
233 : }
234 :
235 : // Compute node map (given elem node, supply face node)
236 : _node_map.clear();
237 697064 : const unsigned int num_node_elem = _current_elem->n_nodes();
238 697064 : const Node * const * elem_nodes = _current_elem->get_nodes();
239 697064 : const unsigned int num_node_face = _current_side_elem->n_nodes();
240 697064 : const Node * const * face_nodes = _current_side_elem->get_nodes();
241 : unsigned int num_found = 0;
242 4993500 : for (unsigned i = 0; i < num_node_elem; ++i)
243 : {
244 18461496 : for (unsigned j = 0; j < num_node_face; ++j)
245 16200774 : if (**(elem_nodes + i) == **(face_nodes + j))
246 : {
247 2732778 : _node_map[i] = j;
248 2732778 : ++num_found;
249 2732778 : break;
250 : }
251 4993500 : if (num_found == num_node_face)
252 : break;
253 : }
254 : }
255 :
256 : void
257 342764 : PressureBase::precalculateQpOffDiagJacobian(const MooseVariableFEBase & /*jvar*/)
258 : {
259 342764 : precalculateQpJacobian();
260 342764 : }
261 :
262 : template class PressureBaseTempl<false>;
263 : template class PressureBaseTempl<true>;
|
