https://mooseframework.inl.gov
PressureBase.C
Go to the documentation of this file.
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>
20 {
22  params.addDeprecatedParam<unsigned int>(
23  "component", "The component for the pressure", "This parameter is no longer necessary");
24  params.addParam<bool>("use_displaced_mesh", true, "Whether to use the displaced mesh.");
25  return params;
26 }
27 
28 template <bool is_ad>
31 {
32  auto params = emptyInputParameters();
33  params.addRequiredCoupledVar("displacements",
34  "The string of displacements suitable for the problem statement");
35  return params;
36 }
37 
38 template <bool is_ad>
40  : PressureBaseParent<is_ad>(parameters),
41  _ndisp(this->coupledComponents("displacements")),
42  _component(
43  [this, &parameters]()
44  {
45  for (unsigned int i = 0; i < _ndisp; ++i)
46  _disp_var.push_back(this->coupled("displacements", i));
47 
48  for (unsigned int i = 0; i < _ndisp; ++i)
49  {
50  if (_var.number() == _disp_var[i])
51  {
52  if (parameters.isParamSetByUser("component") &&
53  i != this->template getParam<unsigned int>("component"))
54  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  return i;
62  }
63  }
64 
65  this->paramError("variable", "The BC variable should be a displacement variable.");
66  }())
67 {
68 }
69 
70 template <bool is_ad>
71 void
73 {
74  auto boundary_ids = this->boundaryIDs();
75  std::set<SubdomainID> block_ids;
76  for (auto bndry_id : boundary_ids)
77  {
78  auto bids = _mesh.getBoundaryConnectedBlocks(bndry_id);
79  block_ids.insert(bids.begin(), bids.end());
80  }
81 
82  if (block_ids.size())
83  _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  _coord_type = _fe_problem.mesh().getUniqueCoordSystem();
89  }
90  for (auto blk_id : block_ids)
91  {
92  if (_coord_type != _fe_problem.getCoordSystem(blk_id))
93  mooseError("The Pressure BC requires subdomains to have the same coordinate system.");
94  }
95 }
96 
97 template <bool is_ad>
100 {
101  return computePressure() * (_normals[_qp](_component) * _test[_i][_qp]);
102 }
103 
105  : PressureBaseTempl<false>(parameters),
106  _q_dxi(nullptr),
107  _q_deta(nullptr),
108  _phi_dxi(nullptr),
109  _phi_deta(nullptr),
110  _use_displaced_mesh(this->template getParam<bool>("use_displaced_mesh")),
111  _fe(libMesh::n_threads())
112 {
113 }
114 
115 Real
116 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  const Real phi_dxi = (*_phi_dxi)[local_j][_qp];
124  const Real phi_deta = _phi_deta ? (*_phi_deta)[local_j][_qp] : 0;
125 
126  const RealGradient & dqdxi = (*_q_dxi)[_qp];
127  const RealGradient out_of_plane(0, 0, 1);
128  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  RealGradient b(dqdxi(1) * dqdeta(2) - dqdxi(2) * dqdeta(1),
132  dqdxi(2) * dqdeta(0) - dqdxi(0) * dqdeta(2),
133  dqdxi(0) * dqdeta(1) - dqdxi(1) * dqdeta(0));
134  const Real inv_length = 1 / (b * _normals[_qp]);
135 
136  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  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  const unsigned int index = 2 - (j + (i + 2) % 3) % 3;
144 
145  const Real variation_b = posneg * (phi_deta * dqdxi(index) - phi_dxi * dqdeta(index));
146 
147  Real rz_term = 0;
148  if (_coord_type == Moose::COORD_RZ && j == _subproblem.getAxisymmetricRadialCoord())
149  {
150  rz_term = _normals[_qp](i) * _phi[_j][_qp] / _q_point[_qp](0);
151  }
152 
153  return computePressure() * _test[_i][_qp] * (inv_length * variation_b + rz_term);
154 }
155 
156 Real
157 PressureBase::computeStiffness(const unsigned int coupled_component)
158 {
159  if (_ndisp > 1)
160  {
161  const std::map<unsigned int, unsigned int>::iterator j_it = _node_map.find(_j);
162  if (_test[_i][_qp] == 0 || j_it == _node_map.end())
163  return 0;
164 
165  return computeFaceStiffness(j_it->second, coupled_component);
166  }
167 
169  {
170  return computePressure() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] *
171  (2 / _q_point[_qp](0));
172  }
173 
175  {
176  return computePressure() * _normals[_qp](_component) * _test[_i][_qp] * _phi[_j][_qp] /
177  _q_point[_qp](0);
178  }
179 
180  return 0;
181 }
182 
183 Real
185 {
188 
189  return 0;
190 }
191 
192 Real
193 PressureBase::computeQpOffDiagJacobian(const unsigned int jvar_num)
194 {
196  for (unsigned int j = 0; j < _ndisp; ++j)
197  if (jvar_num == _disp_var[j])
198  return computeStiffness(j);
199 
200  return 0;
201 }
202 
203 void
205 {
206  if (_ndisp == 1)
207  return;
208 
209  if (_fe[_tid] == nullptr)
210  {
211  const unsigned int dim = _sys.mesh().dimension() - 1;
212  QBase * const & qrule = _assembly.writeableQRuleFace();
213  _fe[_tid] = FEBase::build(dim, _var.feType());
214  _fe[_tid]->attach_quadrature_rule(qrule);
215  }
216 
217  _q_dxi = &_fe[_tid]->get_dxyzdxi();
218  _phi_dxi = &_fe[_tid]->get_dphidxi();
220  {
221  _q_deta = &_fe[_tid]->get_dxyzdeta();
222  _phi_deta = &_fe[_tid]->get_dphideta();
223  }
224 
225  _fe[_tid]->reinit(_current_side_elem);
226 
228  {
229  if (_q_deta->empty())
230  _q_deta = nullptr;
231  if (_phi_deta->empty())
232  _phi_deta = nullptr;
233  }
234 
235  // Compute node map (given elem node, supply face node)
236  _node_map.clear();
237  const unsigned int num_node_elem = _current_elem->n_nodes();
238  const Node * const * elem_nodes = _current_elem->get_nodes();
239  const unsigned int num_node_face = _current_side_elem->n_nodes();
240  const Node * const * face_nodes = _current_side_elem->get_nodes();
241  unsigned int num_found = 0;
242  for (unsigned i = 0; i < num_node_elem; ++i)
243  {
244  for (unsigned j = 0; j < num_node_face; ++j)
245  if (**(elem_nodes + i) == **(face_nodes + j))
246  {
247  _node_map[i] = j;
248  ++num_found;
249  break;
250  }
251  if (num_found == num_node_face)
252  break;
253  }
254 }
255 
256 void
258 {
260 }
261 
262 template class PressureBaseTempl<false>;
263 template class PressureBaseTempl<true>;
Moose::GenericType< Real, is_ad > GenericReal
virtual GenericReal< is_ad > computeQpResidual() override final
Definition: PressureBase.C:99
unsigned int n_threads()
void addDeprecatedParam(const std::string &name, const T &value, const std::string &doc_string, const std::string &deprecation_message)
typename std::conditional< is_ad, ADIntegratedBC, IntegratedBC >::type PressureBaseParent
Pressure applies a pressure on a given boundary in the direction defined by component.
Definition: PressureBase.h:26
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const std::vector< std::vector< Real > > * _phi_dxi
Definition: PressureBase.h:86
void mooseError(Args &&... args)
unsigned int dim
const bool _use_displaced_mesh
Definition: PressureBase.h:88
virtual Real computeQpOffDiagJacobian(const unsigned int jvar_num) override final
Definition: PressureBase.C:193
Moose::CoordinateSystemType _coord_type
Coordinate system type.
Definition: PressureBase.h:54
Real computeStiffness(const unsigned int coupled_component)
Definition: PressureBase.C:157
std::vector< unsigned int > _disp_var
Variable numbers of coupled displacement variables.
Definition: PressureBase.h:45
virtual GenericReal< is_ad > computePressure() const=0
COORD_RSPHERICAL
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
const std::vector< RealGradient > * _q_dxi
Definition: PressureBase.h:84
InputParameters emptyInputParameters()
static InputParameters actionParams()
Definition: PressureBase.C:30
std::vector< std::unique_ptr< FEBase > > _fe
Definition: PressureBase.h:91
void mooseInfo(Args &&... args)
InputParameters validParams()
const unsigned int _component
displacement component to apply the bc to
Definition: PressureBase.h:51
virtual void precalculateQpJacobian() override final
Definition: PressureBase.C:204
bool isParamSetByUser(const std::string &name) const
static InputParameters validParams()
Definition: PressureBase.C:19
virtual void initialSetup() override
Definition: PressureBase.C:72
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< RealGradient > * _q_deta
Definition: PressureBase.h:85
const unsigned int _ndisp
Number of displacement variables.
Definition: PressureBase.h:48
std::map< unsigned int, unsigned int > _node_map
Definition: PressureBase.h:93
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
virtual Real computeQpJacobian() override final
Definition: PressureBase.C:184
const std::vector< std::vector< Real > > * _phi_deta
Definition: PressureBase.h:87
PressureBaseTempl(const InputParameters &parameters)
Definition: PressureBase.C:39
Real computeFaceStiffness(const unsigned int local_j, const unsigned int coupled_component)
Definition: PressureBase.C:116
PressureBase(const InputParameters &parameters)
Definition: PressureBase.C:104
virtual void precalculateQpOffDiagJacobian(const MooseVariableFEBase &jvar) override final
Definition: PressureBase.C:257