https://mooseframework.inl.gov
Assembly.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
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 "Assembly.h"
11 
12 // MOOSE includes
13 #include "SubProblem.h"
14 #include "ArbitraryQuadrature.h"
15 #include "SystemBase.h"
16 #include "MooseTypes.h"
17 #include "MooseMesh.h"
18 #include "MooseVariableFE.h"
19 #include "MooseVariableScalar.h"
20 #include "XFEMInterface.h"
21 #include "DisplacedSystem.h"
22 #include "MooseMeshUtils.h"
23 
24 // libMesh
25 #include "libmesh/coupling_matrix.h"
26 #include "libmesh/dof_map.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/equation_systems.h"
29 #include "libmesh/fe_interface.h"
30 #include "libmesh/node.h"
31 #include "libmesh/quadrature_gauss.h"
32 #include "libmesh/sparse_matrix.h"
33 #include "libmesh/tensor_value.h"
34 #include "libmesh/vector_value.h"
35 #include "libmesh/fe.h"
36 #include "libmesh/static_condensation.h"
37 
38 using namespace libMesh;
39 
40 template <typename P, typename C>
41 void
43  const SubdomainID sub_id,
44  const P & point,
45  C & factor,
46  const SubdomainID neighbor_sub_id)
47 {
48  coordTransformFactor(s.mesh(), sub_id, point, factor, neighbor_sub_id);
49 }
50 
51 template <typename P, typename C>
52 void
54  const SubdomainID sub_id,
55  const P & point,
56  C & factor,
57  const SubdomainID libmesh_dbg_var(neighbor_sub_id))
58 {
59  mooseAssert(neighbor_sub_id != libMesh::Elem::invalid_subdomain_id
60  ? mesh.getCoordSystem(sub_id) == mesh.getCoordSystem(neighbor_sub_id)
61  : true,
62  "Coordinate systems must be the same between element and neighbor");
63  const auto coord_type = mesh.getCoordSystem(sub_id);
64 
65  if (coord_type == Moose::COORD_RZ)
66  {
67  if (mesh.usingGeneralAxisymmetricCoordAxes())
68  {
69  const auto & axis = mesh.getGeneralAxisymmetricCoordAxis(sub_id);
71  }
72  else
74  point, factor, coord_type, mesh.getAxisymmetricRadialCoord());
75  }
76  else
78 }
79 
81  : _sys(sys),
82  _subproblem(_sys.subproblem()),
83  _displaced(dynamic_cast<DisplacedSystem *>(&sys) ? true : false),
84  _nonlocal_cm(_subproblem.nonlocalCouplingMatrix(_sys.number())),
85  _computing_residual(_subproblem.currentlyComputingResidual()),
86  _computing_jacobian(_subproblem.currentlyComputingJacobian()),
87  _computing_residual_and_jacobian(_subproblem.currentlyComputingResidualAndJacobian()),
88  _dof_map(_sys.dofMap()),
89  _tid(tid),
90  _mesh(sys.mesh()),
91  _mesh_dimension(_mesh.dimension()),
92  _helper_type(_mesh.hasSecondOrderElements() ? SECOND : FIRST, LAGRANGE),
93  _user_added_fe_of_helper_type(false),
94  _user_added_fe_face_of_helper_type(false),
95  _user_added_fe_face_neighbor_of_helper_type(false),
96  _user_added_fe_neighbor_of_helper_type(false),
97  _user_added_fe_lower_of_helper_type(false),
98  _building_helpers(false),
99  _current_qrule(nullptr),
100  _current_qrule_volume(nullptr),
101  _current_qrule_arbitrary(nullptr),
102  _coord_type(Moose::COORD_XYZ),
103  _current_qrule_face(nullptr),
104  _current_qface_arbitrary(nullptr),
105  _current_qrule_neighbor(nullptr),
106  _need_JxW_neighbor(false),
107  _qrule_msm(nullptr),
108  _custom_mortar_qrule(false),
109  _current_qrule_lower(nullptr),
110 
111  _current_elem(nullptr),
112  _current_elem_volume(0),
113  _current_side(0),
114  _current_side_elem(nullptr),
115  _current_side_volume(0),
116  _current_neighbor_elem(nullptr),
117  _current_neighbor_side(0),
118  _current_neighbor_side_elem(nullptr),
119  _need_neighbor_elem_volume(false),
120  _current_neighbor_volume(0),
121  _current_node(nullptr),
122  _current_neighbor_node(nullptr),
123  _current_elem_volume_computed(false),
124  _current_side_volume_computed(false),
125 
126  _current_lower_d_elem(nullptr),
127  _current_neighbor_lower_d_elem(nullptr),
128  _need_lower_d_elem_volume(false),
129  _need_neighbor_lower_d_elem_volume(false),
130  _need_dual(false),
131 
132  _residual_vector_tags(_subproblem.getVectorTags(Moose::VECTOR_TAG_RESIDUAL)),
133  _cached_residual_values(2), // The 2 is for TIME and NONTIME
134  _cached_residual_rows(2), // The 2 is for TIME and NONTIME
135  _max_cached_residuals(0),
136  _max_cached_jacobians(0),
137 
138  _block_diagonal_matrix(false),
139  _calculate_xyz(false),
140  _calculate_face_xyz(false),
141  _calculate_curvatures(false),
142  _calculate_ad_coord(false),
143  _have_p_refinement(false),
144  _sc(nullptr)
145 {
146  const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
147  _building_helpers = true;
148  // Build fe's for the helpers
149  buildFE(FEType(helper_order, LAGRANGE));
150  buildFaceFE(FEType(helper_order, LAGRANGE));
151  buildNeighborFE(FEType(helper_order, LAGRANGE));
152  buildFaceNeighborFE(FEType(helper_order, LAGRANGE));
153  buildLowerDFE(FEType(helper_order, LAGRANGE));
154  _building_helpers = false;
155 
156  // Build an FE helper object for this type for each dimension up to the dimension of the current
157  // mesh
158  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
159  {
160  _holder_fe_helper[dim] = _fe[dim][FEType(helper_order, LAGRANGE)];
164  }
165 
166  for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
168 
169  // request phi, dphi, xyz, JxW, etc. data
171 
172  // For 3D mortar, mortar segments are always TRI3 elements so we want FIRST LAGRANGE regardless
173  // of discretization
174  _fe_msm = (_mesh_dimension == 2)
177  // This FE object should not take part in p-refinement
178  _fe_msm->add_p_level_in_reinit(false);
179  _JxW_msm = &_fe_msm->get_JxW();
180  // Prerequest xyz so that it is computed for _fe_msm so that it can be used for calculating
181  // _coord_msm
182  _fe_msm->get_xyz();
183 
186 }
187 
189 {
190  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
191  for (auto & it : _fe[dim])
192  delete it.second;
193 
194  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
195  for (auto & it : _fe_face[dim])
196  delete it.second;
197 
198  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
199  for (auto & it : _fe_neighbor[dim])
200  delete it.second;
201 
202  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
203  for (auto & it : _fe_face_neighbor[dim])
204  delete it.second;
205 
206  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
207  for (auto & it : _fe_lower[dim])
208  delete it.second;
209 
210  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
211  for (auto & it : _vector_fe[dim])
212  delete it.second;
213 
214  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
215  for (auto & it : _vector_fe_face[dim])
216  delete it.second;
217 
218  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
219  for (auto & it : _vector_fe_neighbor[dim])
220  delete it.second;
221 
222  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
223  for (auto & it : _vector_fe_face_neighbor[dim])
224  delete it.second;
225 
226  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
227  for (auto & it : _vector_fe_lower[dim])
228  delete it.second;
229 
230  for (auto & it : _ad_grad_phi_data)
231  it.second.release();
232 
233  for (auto & it : _ad_vector_grad_phi_data)
234  it.second.release();
235 
236  for (auto & it : _ad_grad_phi_data_face)
237  it.second.release();
238 
239  for (auto & it : _ad_vector_grad_phi_data_face)
240  it.second.release();
241 
243 
244  _coord.release();
247 
248  _ad_JxW.release();
255  _ad_coord.release();
256 
257  delete _qrule_msm;
258 }
259 
260 const MooseArray<Real> &
262 {
263  _need_JxW_neighbor = true;
264  return _current_JxW_neighbor;
265 }
266 
267 void
269 {
270  if (!_building_helpers && type == _helper_type)
272 
273  if (!_fe_shape_data[type])
274  _fe_shape_data[type] = std::make_unique<FEShapeData>();
275 
276  // Build an FE object for this type for each dimension up to the dimension of the current mesh
277  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
278  {
279  if (!_fe[dim][type])
280  _fe[dim][type] = FEGenericBase<Real>::build(dim, type).release();
281 
282  _fe[dim][type]->get_phi();
283  _fe[dim][type]->get_dphi();
284  // Pre-request xyz. We have always computed xyz, but due to
285  // recent optimizations in libmesh, we now need to explicity
286  // request it, since apps (Yak) may rely on it being computed.
287  _fe[dim][type]->get_xyz();
288  if (_need_second_derivative.count(type))
289  _fe[dim][type]->get_d2phi();
290  }
291 }
292 
293 void
295 {
296  if (!_building_helpers && type == _helper_type)
298 
299  if (!_fe_shape_data_face[type])
300  _fe_shape_data_face[type] = std::make_unique<FEShapeData>();
301 
302  // Build an FE object for this type for each dimension up to the dimension of the current mesh
303  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
304  {
305  if (!_fe_face[dim][type])
306  _fe_face[dim][type] = FEGenericBase<Real>::build(dim, type).release();
307 
308  _fe_face[dim][type]->get_phi();
309  _fe_face[dim][type]->get_dphi();
310  if (_need_second_derivative.count(type))
311  _fe_face[dim][type]->get_d2phi();
312  }
313 }
314 
315 void
317 {
318  if (!_building_helpers && type == _helper_type)
320 
321  if (!_fe_shape_data_neighbor[type])
322  _fe_shape_data_neighbor[type] = std::make_unique<FEShapeData>();
323 
324  // Build an FE object for this type for each dimension up to the dimension of the current mesh
325  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
326  {
327  if (!_fe_neighbor[dim][type])
328  _fe_neighbor[dim][type] = FEGenericBase<Real>::build(dim, type).release();
329 
330  _fe_neighbor[dim][type]->get_phi();
331  _fe_neighbor[dim][type]->get_dphi();
332  if (_need_second_derivative_neighbor.count(type))
333  _fe_neighbor[dim][type]->get_d2phi();
334  }
335 }
336 
337 void
339 {
340  if (!_building_helpers && type == _helper_type)
342 
343  if (!_fe_shape_data_face_neighbor[type])
344  _fe_shape_data_face_neighbor[type] = std::make_unique<FEShapeData>();
345 
346  // Build an FE object for this type for each dimension up to the dimension of the current mesh
347  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
348  {
349  if (!_fe_face_neighbor[dim][type])
350  _fe_face_neighbor[dim][type] = FEGenericBase<Real>::build(dim, type).release();
351 
352  _fe_face_neighbor[dim][type]->get_phi();
353  _fe_face_neighbor[dim][type]->get_dphi();
354  if (_need_second_derivative_neighbor.count(type))
355  _fe_face_neighbor[dim][type]->get_d2phi();
356  }
357 }
358 
359 void
361 {
362  if (!_building_helpers && type == _helper_type)
364 
365  if (!_fe_shape_data_lower[type])
366  _fe_shape_data_lower[type] = std::make_unique<FEShapeData>();
367 
368  // Build an FE object for this type for each dimension up to the dimension of
369  // the current mesh minus one (because this is for lower-dimensional
370  // elements!)
371  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
372  {
373  if (!_fe_lower[dim][type])
374  _fe_lower[dim][type] = FEGenericBase<Real>::build(dim, type).release();
375 
376  _fe_lower[dim][type]->get_phi();
377  _fe_lower[dim][type]->get_dphi();
378  if (_need_second_derivative.count(type))
379  _fe_lower[dim][type]->get_d2phi();
380  }
381 }
382 
383 void
385 {
386  if (!_fe_shape_data_dual_lower[type])
387  _fe_shape_data_dual_lower[type] = std::make_unique<FEShapeData>();
388 
389  // Build an FE object for this type for each dimension up to the dimension of
390  // the current mesh minus one (because this is for lower-dimensional
391  // elements!)
392  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
393  {
394  if (!_fe_lower[dim][type])
395  _fe_lower[dim][type] = FEGenericBase<Real>::build(dim, type).release();
396 
397  _fe_lower[dim][type]->get_dual_phi();
398  _fe_lower[dim][type]->get_dual_dphi();
399  if (_need_second_derivative.count(type))
400  _fe_lower[dim][type]->get_dual_d2phi();
401  }
402 }
403 
404 void
406 {
407  if (!_vector_fe_shape_data_lower[type])
408  _vector_fe_shape_data_lower[type] = std::make_unique<VectorFEShapeData>();
409 
410  // Build an FE object for this type for each dimension up to the dimension of
411  // the current mesh minus one (because this is for lower-dimensional
412  // elements!)
413  unsigned int dim = ((type.family == LAGRANGE_VEC) || (type.family == MONOMIAL_VEC)) ? 0 : 2;
414  const auto ending_dim = cast_int<unsigned int>(_mesh_dimension - 1);
415  if (ending_dim < dim)
416  return;
417  for (; dim <= ending_dim; dim++)
418  {
419  if (!_vector_fe_lower[dim][type])
420  _vector_fe_lower[dim][type] = FEVectorBase::build(dim, type).release();
421 
422  _vector_fe_lower[dim][type]->get_phi();
423  _vector_fe_lower[dim][type]->get_dphi();
424  if (_need_second_derivative.count(type))
425  _vector_fe_lower[dim][type]->get_d2phi();
426  }
427 }
428 
429 void
431 {
433  _vector_fe_shape_data_dual_lower[type] = std::make_unique<VectorFEShapeData>();
434 
435  // Build an FE object for this type for each dimension up to the dimension of
436  // the current mesh minus one (because this is for lower-dimensional
437  // elements!)
438  unsigned int dim = ((type.family == LAGRANGE_VEC) || (type.family == MONOMIAL_VEC)) ? 0 : 2;
439  const auto ending_dim = cast_int<unsigned int>(_mesh_dimension - 1);
440  if (ending_dim < dim)
441  return;
442  for (; dim <= ending_dim; dim++)
443  {
444  if (!_vector_fe_lower[dim][type])
445  _vector_fe_lower[dim][type] = FEVectorBase::build(dim, type).release();
446 
447  _vector_fe_lower[dim][type]->get_dual_phi();
448  _vector_fe_lower[dim][type]->get_dual_dphi();
449  if (_need_second_derivative.count(type))
450  _vector_fe_lower[dim][type]->get_dual_d2phi();
451  }
452 }
453 
454 void
456 {
457  if (!_vector_fe_shape_data[type])
458  _vector_fe_shape_data[type] = std::make_unique<VectorFEShapeData>();
459 
460  // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
461  unsigned int min_dim;
462  if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
463  type.family == L2_RAVIART_THOMAS)
464  min_dim = 2;
465  else
466  min_dim = 0;
467 
468  // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
469  // current mesh
470  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
471  {
472  if (!_vector_fe[dim][type])
473  _vector_fe[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
474 
475  _vector_fe[dim][type]->get_phi();
476  _vector_fe[dim][type]->get_dphi();
477  if (_need_curl.count(type))
478  _vector_fe[dim][type]->get_curl_phi();
479  if (_need_div.count(type))
480  _vector_fe[dim][type]->get_div_phi();
481  _vector_fe[dim][type]->get_xyz();
482  }
483 }
484 
485 void
487 {
488  if (!_vector_fe_shape_data_face[type])
489  _vector_fe_shape_data_face[type] = std::make_unique<VectorFEShapeData>();
490 
491  // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
492  unsigned int min_dim;
493  if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
494  type.family == L2_RAVIART_THOMAS)
495  min_dim = 2;
496  else
497  min_dim = 0;
498 
499  // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
500  // current mesh
501  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
502  {
503  if (!_vector_fe_face[dim][type])
504  _vector_fe_face[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
505 
506  _vector_fe_face[dim][type]->get_phi();
507  _vector_fe_face[dim][type]->get_dphi();
508  if (_need_curl.count(type))
509  _vector_fe_face[dim][type]->get_curl_phi();
510  if (_need_face_div.count(type))
511  _vector_fe_face[dim][type]->get_div_phi();
512  }
513 }
514 
515 void
517 {
519  _vector_fe_shape_data_neighbor[type] = std::make_unique<VectorFEShapeData>();
520 
521  // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
522  unsigned int min_dim;
523  if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
524  type.family == L2_RAVIART_THOMAS)
525  min_dim = 2;
526  else
527  min_dim = 0;
528 
529  // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
530  // current mesh
531  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
532  {
533  if (!_vector_fe_neighbor[dim][type])
534  _vector_fe_neighbor[dim][type] = FEGenericBase<VectorValue<Real>>::build(dim, type).release();
535 
536  _vector_fe_neighbor[dim][type]->get_phi();
537  _vector_fe_neighbor[dim][type]->get_dphi();
538  if (_need_curl.count(type))
539  _vector_fe_neighbor[dim][type]->get_curl_phi();
540  if (_need_neighbor_div.count(type))
541  _vector_fe_neighbor[dim][type]->get_div_phi();
542  }
543 }
544 
545 void
547 {
549  _vector_fe_shape_data_face_neighbor[type] = std::make_unique<VectorFEShapeData>();
550 
551  // Note that NEDELEC_ONE and RAVIART_THOMAS elements can only be built for dimension > 2
552  unsigned int min_dim;
553  if (type.family == NEDELEC_ONE || type.family == RAVIART_THOMAS ||
554  type.family == L2_RAVIART_THOMAS)
555  min_dim = 2;
556  else
557  min_dim = 0;
558 
559  // Build an FE object for this type for each dimension from the min_dim up to the dimension of the
560  // current mesh
561  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
562  {
563  if (!_vector_fe_face_neighbor[dim][type])
565  FEGenericBase<VectorValue<Real>>::build(dim, type).release();
566 
567  _vector_fe_face_neighbor[dim][type]->get_phi();
568  _vector_fe_face_neighbor[dim][type]->get_dphi();
569  if (_need_curl.count(type))
570  _vector_fe_face_neighbor[dim][type]->get_curl_phi();
571  if (_need_face_neighbor_div.count(type))
572  _vector_fe_face_neighbor[dim][type]->get_div_phi();
573  }
574 }
575 
576 void
578 {
579  auto & qdefault = _qrules[Moose::ANY_BLOCK_ID];
580  mooseAssert(qdefault.size() > 0, "default quadrature must be initialized before order bumps");
581 
582  unsigned int ndims = _mesh_dimension + 1; // must account for 0-dimensional quadrature.
583  auto & qvec = _qrules[block];
584  if (qvec.size() != ndims || !qvec[0].vol)
585  createQRules(qdefault[0].vol->type(),
586  qdefault[0].arbitrary_vol->get_order(),
587  volume_order,
588  qdefault[0].face->get_order(),
589  block);
590  else if (qvec[0].vol->get_order() < volume_order)
591  createQRules(qvec[0].vol->type(),
592  qvec[0].arbitrary_vol->get_order(),
593  volume_order,
594  qvec[0].face->get_order(),
595  block);
596  // otherwise do nothing - quadrature order is already as high as requested
597 }
598 
599 void
601 {
602  auto & qdefault = _qrules[Moose::ANY_BLOCK_ID];
603  mooseAssert(qdefault.size() > 0, "default quadrature must be initialized before order bumps");
604 
605  unsigned int ndims = _mesh_dimension + 1; // must account for 0-dimensional quadrature.
606  auto & qvec = _qrules[block];
607  if (qvec.size() != ndims || !qvec[0].vol)
608  createQRules(qdefault[0].vol->type(), order, order, order, block);
609  else if (qvec[0].vol->get_order() < order || qvec[0].face->get_order() < order)
610  createQRules(qvec[0].vol->type(),
611  std::max(order, qvec[0].arbitrary_vol->get_order()),
612  std::max(order, qvec[0].vol->get_order()),
613  std::max(order, qvec[0].face->get_order()),
614  block);
615  // otherwise do nothing - quadrature order is already as high as requested
616 }
617 
618 void
620  Order order,
621  Order volume_order,
622  Order face_order,
623  SubdomainID block,
624  bool allow_negative_qweights)
625 {
626  auto & qvec = _qrules[block];
627  unsigned int ndims = _mesh_dimension + 1; // must account for 0-dimensional quadrature.
628  if (qvec.size() != ndims)
629  qvec.resize(ndims);
630 
631  for (unsigned int i = 0; i < qvec.size(); i++)
632  {
633  int dim = i;
634  auto & q = qvec[dim];
635  q.vol = QBase::build(type, dim, volume_order);
636  q.vol->allow_rules_with_negative_weights = allow_negative_qweights;
637  q.face = QBase::build(type, dim - 1, face_order);
638  q.face->allow_rules_with_negative_weights = allow_negative_qweights;
639  q.fv_face = QBase::build(QMONOMIAL, dim - 1, CONSTANT);
640  q.fv_face->allow_rules_with_negative_weights = allow_negative_qweights;
641  q.neighbor = std::make_unique<ArbitraryQuadrature>(dim - 1, face_order);
642  q.neighbor->allow_rules_with_negative_weights = allow_negative_qweights;
643  q.arbitrary_vol = std::make_unique<ArbitraryQuadrature>(dim, order);
644  q.arbitrary_vol->allow_rules_with_negative_weights = allow_negative_qweights;
645  q.arbitrary_face = std::make_unique<ArbitraryQuadrature>(dim - 1, face_order);
646  q.arbitrary_face->allow_rules_with_negative_weights = allow_negative_qweights;
647  }
648 
649  delete _qrule_msm;
650  _custom_mortar_qrule = false;
651  _qrule_msm = QBase::build(type, _mesh_dimension - 1, face_order).release();
652  _qrule_msm->allow_rules_with_negative_weights = allow_negative_qweights;
653  _fe_msm->attach_quadrature_rule(_qrule_msm);
654 }
655 
656 void
657 Assembly::setVolumeQRule(QBase * qrule, unsigned int dim)
658 {
659  _current_qrule = qrule;
660 
661  if (qrule) // Don't set a NULL qrule
662  {
663  for (auto & it : _fe[dim])
664  it.second->attach_quadrature_rule(qrule);
665  for (auto & it : _vector_fe[dim])
666  it.second->attach_quadrature_rule(qrule);
667  if (!_unique_fe_helper.empty())
668  {
669  mooseAssert(dim < _unique_fe_helper.size(), "We should not be indexing out of bounds");
670  _unique_fe_helper[dim]->attach_quadrature_rule(qrule);
671  }
672  }
673 }
674 
675 void
676 Assembly::setFaceQRule(QBase * qrule, unsigned int dim)
677 {
678  _current_qrule_face = qrule;
679 
680  for (auto & it : _fe_face[dim])
681  it.second->attach_quadrature_rule(qrule);
682  for (auto & it : _vector_fe_face[dim])
683  it.second->attach_quadrature_rule(qrule);
684  if (!_unique_fe_face_helper.empty())
685  {
686  mooseAssert(dim < _unique_fe_face_helper.size(), "We should not be indexing out of bounds");
687  _unique_fe_face_helper[dim]->attach_quadrature_rule(qrule);
688  }
689 }
690 
691 void
692 Assembly::setLowerQRule(QBase * qrule, unsigned int dim)
693 {
694  // The lower-dimensional quadrature rule matches the face quadrature rule
695  setFaceQRule(qrule, dim);
696 
697  _current_qrule_lower = qrule;
698 
699  for (auto & it : _fe_lower[dim])
700  it.second->attach_quadrature_rule(qrule);
701  for (auto & it : _vector_fe_lower[dim])
702  it.second->attach_quadrature_rule(qrule);
703  if (!_unique_fe_lower_helper.empty())
704  {
705  mooseAssert(dim < _unique_fe_lower_helper.size(), "We should not be indexing out of bounds");
706  _unique_fe_lower_helper[dim]->attach_quadrature_rule(qrule);
707  }
708 }
709 
710 void
711 Assembly::setNeighborQRule(QBase * qrule, unsigned int dim)
712 {
713  _current_qrule_neighbor = qrule;
714 
715  for (auto & it : _fe_face_neighbor[dim])
716  it.second->attach_quadrature_rule(qrule);
717  for (auto & it : _vector_fe_face_neighbor[dim])
718  it.second->attach_quadrature_rule(qrule);
719  if (!_unique_fe_face_neighbor_helper.empty())
720  {
721  mooseAssert(dim < _unique_fe_face_neighbor_helper.size(),
722  "We should not be indexing out of bounds");
723  _unique_fe_face_neighbor_helper[dim]->attach_quadrature_rule(qrule);
724  }
725 }
726 
727 void
729 {
730  _current_qrule = nullptr;
731  _current_qrule_face = nullptr;
732  _current_qrule_lower = nullptr;
733  _current_qrule_neighbor = nullptr;
734 }
735 
736 void
738 {
739  if (order != _qrule_msm->get_order())
740  {
741  // If custom mortar qrule has not yet been specified
743  {
744  _custom_mortar_qrule = true;
745  const unsigned int dim = _qrule_msm->get_dim();
746  const QuadratureType type = _qrule_msm->type();
747  delete _qrule_msm;
748 
749  _qrule_msm = QBase::build(type, dim, order).release();
750  _fe_msm->attach_quadrature_rule(_qrule_msm);
751  }
752  else
753  mooseError("Mortar quadrature_order: ",
754  order,
755  " does not match previously specified quadrature_order: ",
757  ". Quadrature_order (when specified) must match for all mortar constraints.");
758  }
759 }
760 
761 void
763 {
764  unsigned int dim = elem->dim();
765 
766  for (const auto & it : _fe[dim])
767  {
768  FEBase & fe = *it.second;
769  const FEType & fe_type = it.first;
770 
771  _current_fe[fe_type] = &fe;
772 
773  FEShapeData & fesd = *_fe_shape_data[fe_type];
774 
775  fe.reinit(elem);
776 
777  fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_phi()));
778  fesd._grad_phi.shallowCopy(
779  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_dphi()));
780  if (_need_second_derivative.count(fe_type))
781  fesd._second_phi.shallowCopy(
782  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe.get_d2phi()));
783  }
784  for (const auto & it : _vector_fe[dim])
785  {
786  FEVectorBase & fe = *it.second;
787  const FEType & fe_type = it.first;
788 
789  _current_vector_fe[fe_type] = &fe;
790 
791  VectorFEShapeData & fesd = *_vector_fe_shape_data[fe_type];
792 
793  fe.reinit(elem);
794 
795  fesd._phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_phi()));
796  fesd._grad_phi.shallowCopy(
797  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe.get_dphi()));
798  if (_need_second_derivative.count(fe_type))
799  fesd._second_phi.shallowCopy(
800  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe.get_d2phi()));
801  if (_need_curl.count(fe_type))
802  fesd._curl_phi.shallowCopy(
803  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_curl_phi()));
804  if (_need_div.count(fe_type))
805  fesd._div_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_div_phi()));
806  }
807  if (!_unique_fe_helper.empty())
808  {
809  mooseAssert(dim < _unique_fe_helper.size(), "We should be in bounds here");
810  _unique_fe_helper[dim]->reinit(elem);
811  }
812 
813  // During that last loop the helper objects will have been reinitialized as well
814  // We need to dig out the q_points and JxW from it.
816  const_cast<std::vector<Point> &>(_holder_fe_helper[dim]->get_xyz()));
817  _current_JxW.shallowCopy(const_cast<std::vector<Real> &>(_holder_fe_helper[dim]->get_JxW()));
818 
820  {
821  auto n_qp = _current_qrule->n_points();
823  if (_displaced)
824  {
825  const auto & qw = _current_qrule->get_weights();
826  for (unsigned int qp = 0; qp != n_qp; qp++)
828  }
829  else
830  {
831  for (unsigned qp = 0; qp < n_qp; ++qp)
832  _ad_JxW[qp] = _current_JxW[qp];
833  if (_calculate_xyz)
834  for (unsigned qp = 0; qp < n_qp; ++qp)
836  }
837 
838  for (const auto & it : _fe[dim])
839  {
840  FEBase & fe = *it.second;
841  auto fe_type = it.first;
842  auto num_shapes = FEInterface::n_shape_functions(fe_type, elem);
843  auto & grad_phi = _ad_grad_phi_data[fe_type];
844 
845  grad_phi.resize(num_shapes);
846  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
847  grad_phi[i].resize(n_qp);
848 
849  if (_displaced)
850  computeGradPhiAD(elem, n_qp, grad_phi, &fe);
851  else
852  {
853  const auto & regular_grad_phi = _fe_shape_data[fe_type]->_grad_phi;
854  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
855  for (unsigned qp = 0; qp < n_qp; ++qp)
856  grad_phi[i][qp] = regular_grad_phi[i][qp];
857  }
858  }
859  for (const auto & it : _vector_fe[dim])
860  {
861  FEVectorBase & fe = *it.second;
862  auto fe_type = it.first;
863  auto num_shapes = FEInterface::n_shape_functions(fe_type, elem);
864  auto & grad_phi = _ad_vector_grad_phi_data[fe_type];
865 
866  grad_phi.resize(num_shapes);
867  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
868  grad_phi[i].resize(n_qp);
869 
870  if (_displaced)
871  computeGradPhiAD(elem, n_qp, grad_phi, &fe);
872  else
873  {
874  const auto & regular_grad_phi = _vector_fe_shape_data[fe_type]->_grad_phi;
875  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
876  for (unsigned qp = 0; qp < n_qp; ++qp)
877  grad_phi[i][qp] = regular_grad_phi[i][qp];
878  }
879  }
880  }
881 
882  auto n = numExtraElemIntegers();
883  for (auto i : make_range(n))
886 
887  if (_xfem != nullptr)
889 }
890 
891 template <typename OutputType>
892 void
894  unsigned int n_qp,
897 {
898  // This function relies on the fact that FE::reinit has already been called. FE::reinit will
899  // importantly have already called FEMap::init_shape_functions which will have computed
900  // these quantities at the integration/quadrature points: dphidxi,
901  // dphideta, and dphidzeta (e.g. \nabla phi w.r.t. reference coordinates). These *phi* quantities
902  // are independent of mesh displacements when using a quadrature rule.
903  //
904  // Note that a user could have specified custom integration points (e.g. independent of a
905  // quadrature rule) which could very well depend on displacements. In that case even the *phi*
906  // quantities from the above paragraph would be a function of the displacements and we would be
907  // missing that derivative information in the calculations below
908 
909  auto dim = elem->dim();
910  const auto & dphidxi = fe->get_dphidxi();
911  const auto & dphideta = fe->get_dphideta();
912  const auto & dphidzeta = fe->get_dphidzeta();
913  auto num_shapes = grad_phi.size();
914 
915  switch (dim)
916  {
917  case 0:
918  {
919  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
920  for (unsigned qp = 0; qp < n_qp; ++qp)
921  grad_phi[i][qp] = 0;
922  break;
923  }
924 
925  case 1:
926  {
927  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
928  for (unsigned qp = 0; qp < n_qp; ++qp)
929  {
930  grad_phi[i][qp].slice(0) = dphidxi[i][qp] * _ad_dxidx_map[qp];
931  grad_phi[i][qp].slice(1) = dphidxi[i][qp] * _ad_dxidy_map[qp];
932  grad_phi[i][qp].slice(2) = dphidxi[i][qp] * _ad_dxidz_map[qp];
933  }
934  break;
935  }
936 
937  case 2:
938  {
939  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
940  for (unsigned qp = 0; qp < n_qp; ++qp)
941  {
942  grad_phi[i][qp].slice(0) =
943  dphidxi[i][qp] * _ad_dxidx_map[qp] + dphideta[i][qp] * _ad_detadx_map[qp];
944  grad_phi[i][qp].slice(1) =
945  dphidxi[i][qp] * _ad_dxidy_map[qp] + dphideta[i][qp] * _ad_detady_map[qp];
946  grad_phi[i][qp].slice(2) =
947  dphidxi[i][qp] * _ad_dxidz_map[qp] + dphideta[i][qp] * _ad_detadz_map[qp];
948  }
949  break;
950  }
951 
952  case 3:
953  {
954  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
955  for (unsigned qp = 0; qp < n_qp; ++qp)
956  {
957  grad_phi[i][qp].slice(0) = dphidxi[i][qp] * _ad_dxidx_map[qp] +
958  dphideta[i][qp] * _ad_detadx_map[qp] +
959  dphidzeta[i][qp] * _ad_dzetadx_map[qp];
960  grad_phi[i][qp].slice(1) = dphidxi[i][qp] * _ad_dxidy_map[qp] +
961  dphideta[i][qp] * _ad_detady_map[qp] +
962  dphidzeta[i][qp] * _ad_dzetady_map[qp];
963  grad_phi[i][qp].slice(2) = dphidxi[i][qp] * _ad_dxidz_map[qp] +
964  dphideta[i][qp] * _ad_detadz_map[qp] +
965  dphidzeta[i][qp] * _ad_dzetadz_map[qp];
966  }
967  break;
968  }
969  }
970 }
971 
972 void
973 Assembly::resizeADMappingObjects(unsigned int n_qp, unsigned int dim)
974 {
975  _ad_dxyzdxi_map.resize(n_qp);
976  _ad_dxidx_map.resize(n_qp);
977  _ad_dxidy_map.resize(n_qp); // 1D element may live in 2D ...
978  _ad_dxidz_map.resize(n_qp); // ... or 3D
979 
980  if (dim > 1)
981  {
982  _ad_dxyzdeta_map.resize(n_qp);
983  _ad_detadx_map.resize(n_qp);
984  _ad_detady_map.resize(n_qp);
985  _ad_detadz_map.resize(n_qp);
986 
987  if (dim > 2)
988  {
989  _ad_dxyzdzeta_map.resize(n_qp);
990  _ad_dzetadx_map.resize(n_qp);
991  _ad_dzetady_map.resize(n_qp);
992  _ad_dzetadz_map.resize(n_qp);
993  }
994  }
995 
996  _ad_jac.resize(n_qp);
997  _ad_JxW.resize(n_qp);
998  if (_calculate_xyz)
999  _ad_q_points.resize(n_qp);
1000 }
1001 
1002 void
1004  const std::vector<Real> & qw,
1005  unsigned p,
1006  FEBase * fe)
1007 {
1008  // This function relies on the fact that FE::reinit has already been called. FE::reinit will
1009  // importantly have already called FEMap::init_reference_to_physical_map which will have computed
1010  // these quantities at the integration/quadrature points: phi_map, dphidxi_map,
1011  // dphideta_map, and dphidzeta_map (e.g. phi and \nabla phi w.r.t reference coordinates). *_map is
1012  // used to denote that quantities are in reference to a mapping Lagrange FE object. The FE<Dim,
1013  // LAGRANGE> objects used for mapping will in general have an order matching the order of the
1014  // mesh. These *phi*_map quantities are independent of mesh displacements when using a quadrature
1015  // rule.
1016  //
1017  // Note that a user could have specified custom integration points (e.g. independent of a
1018  // quadrature rule) which could very well depend on displacements. In that case even the *phi*_map
1019  // quantities from the above paragraph would be a function of the displacements and we would be
1020  // missing that derivative information in the calculations below
1021  //
1022  // Important quantities calculated by this method:
1023  // - _ad_JxW;
1024  // - _ad_q_points;
1025  // And the following quantities are important because they are used in the computeGradPhiAD method
1026  // to calculate the shape function gradients with respect to the physical coordinates
1027  // dphi/dphys = dphi/dref * dref/dphys:
1028  // - _ad_dxidx_map;
1029  // - _ad_dxidy_map;
1030  // - _ad_dxidz_map;
1031  // - _ad_detadx_map;
1032  // - _ad_detady_map;
1033  // - _ad_detadz_map;
1034  // - _ad_dzetadx_map;
1035  // - _ad_dzetady_map;
1036  // - _ad_dzetadz_map;
1037  //
1038  // Some final notes. This method will be called both when we are reinit'ing in the volume and on
1039  // faces. When reinit'ing on faces, computation of _ad_JxW will be garbage because we will be
1040  // using dummy quadrature weights. _ad_q_points computation is also currently extraneous during
1041  // face reinit because we compute _ad_q_points_face in the computeFaceMap method. However,
1042  // computation of dref/dphys is absolutely necessary (and the reason we call this method for the
1043  // face case) for both volume and face reinit
1044 
1045  auto dim = elem->dim();
1046  const auto & elem_nodes = elem->get_nodes();
1047  auto num_shapes = FEInterface::n_shape_functions(fe->get_fe_type(), elem);
1048  const auto & phi_map = fe->get_fe_map().get_phi_map();
1049  const auto & dphidxi_map = fe->get_fe_map().get_dphidxi_map();
1050  const auto & dphideta_map = fe->get_fe_map().get_dphideta_map();
1051  const auto & dphidzeta_map = fe->get_fe_map().get_dphidzeta_map();
1052  const auto sys_num = _sys.number();
1053  const bool do_derivatives =
1054  ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
1055 
1056  switch (dim)
1057  {
1058  case 0:
1059  {
1060  _ad_jac[p] = 1.0;
1061  _ad_JxW[p] = qw[p];
1062  if (_calculate_xyz)
1063  _ad_q_points[p] = *elem_nodes[0];
1064  break;
1065  }
1066 
1067  case 1:
1068  {
1069  if (_calculate_xyz)
1070  _ad_q_points[p].zero();
1071 
1072  _ad_dxyzdxi_map[p].zero();
1073 
1074  for (std::size_t i = 0; i < num_shapes; i++)
1075  {
1076  libmesh_assert(elem_nodes[i]);
1077  const Node & node = *elem_nodes[i];
1078  libMesh::VectorValue<ADReal> elem_point = node;
1079  if (do_derivatives)
1080  for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1081  if (node.n_dofs(sys_num, disp_num))
1083  elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1084 
1085  _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1086 
1087  if (_calculate_xyz)
1088  _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1089  }
1090 
1091  _ad_jac[p] = _ad_dxyzdxi_map[p].norm();
1092 
1093  if (_ad_jac[p].value() <= -TOLERANCE * TOLERANCE)
1094  {
1095  static bool failing = false;
1096  if (!failing)
1097  {
1098  failing = true;
1100  libmesh_error_msg("ERROR: negative Jacobian " << _ad_jac[p].value() << " at point index "
1101  << p << " in element " << elem->id());
1102  }
1103  else
1104  return;
1105  }
1106 
1107  const auto jacm2 = 1. / _ad_jac[p] / _ad_jac[p];
1108  _ad_dxidx_map[p] = jacm2 * _ad_dxyzdxi_map[p](0);
1109  _ad_dxidy_map[p] = jacm2 * _ad_dxyzdxi_map[p](1);
1110  _ad_dxidz_map[p] = jacm2 * _ad_dxyzdxi_map[p](2);
1111 
1112  _ad_JxW[p] = _ad_jac[p] * qw[p];
1113 
1114  break;
1115  }
1116 
1117  case 2:
1118  {
1119  if (_calculate_xyz)
1120  _ad_q_points[p].zero();
1121  _ad_dxyzdxi_map[p].zero();
1122  _ad_dxyzdeta_map[p].zero();
1123 
1124  for (std::size_t i = 0; i < num_shapes; i++)
1125  {
1126  libmesh_assert(elem_nodes[i]);
1127  const Node & node = *elem_nodes[i];
1128  libMesh::VectorValue<ADReal> elem_point = node;
1129  if (do_derivatives)
1130  for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1131  if (node.n_dofs(sys_num, disp_num))
1133  elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1134 
1135  _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1136  _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
1137 
1138  if (_calculate_xyz)
1139  _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1140  }
1141 
1142  const auto &dx_dxi = _ad_dxyzdxi_map[p](0), &dx_deta = _ad_dxyzdeta_map[p](0),
1143  &dy_dxi = _ad_dxyzdxi_map[p](1), &dy_deta = _ad_dxyzdeta_map[p](1),
1144  &dz_dxi = _ad_dxyzdxi_map[p](2), &dz_deta = _ad_dxyzdeta_map[p](2);
1145 
1146  const auto g11 = (dx_dxi * dx_dxi + dy_dxi * dy_dxi + dz_dxi * dz_dxi);
1147 
1148  const auto g12 = (dx_dxi * dx_deta + dy_dxi * dy_deta + dz_dxi * dz_deta);
1149 
1150  const auto & g21 = g12;
1151 
1152  const auto g22 = (dx_deta * dx_deta + dy_deta * dy_deta + dz_deta * dz_deta);
1153 
1154  auto det = (g11 * g22 - g12 * g21);
1155 
1156  if (det.value() <= -TOLERANCE * TOLERANCE)
1157  {
1158  static bool failing = false;
1159  if (!failing)
1160  {
1161  failing = true;
1163  libmesh_error_msg("ERROR: negative Jacobian " << det << " at point index " << p
1164  << " in element " << elem->id());
1165  }
1166  else
1167  return;
1168  }
1169  else if (det.value() <= 0.)
1170  det.value() = TOLERANCE * TOLERANCE;
1171 
1172  const auto inv_det = 1. / det;
1173  _ad_jac[p] = std::sqrt(det);
1174 
1175  _ad_JxW[p] = _ad_jac[p] * qw[p];
1176 
1177  const auto g11inv = g22 * inv_det;
1178  const auto g12inv = -g12 * inv_det;
1179  const auto g21inv = -g21 * inv_det;
1180  const auto g22inv = g11 * inv_det;
1181 
1182  _ad_dxidx_map[p] = g11inv * dx_dxi + g12inv * dx_deta;
1183  _ad_dxidy_map[p] = g11inv * dy_dxi + g12inv * dy_deta;
1184  _ad_dxidz_map[p] = g11inv * dz_dxi + g12inv * dz_deta;
1185 
1186  _ad_detadx_map[p] = g21inv * dx_dxi + g22inv * dx_deta;
1187  _ad_detady_map[p] = g21inv * dy_dxi + g22inv * dy_deta;
1188  _ad_detadz_map[p] = g21inv * dz_dxi + g22inv * dz_deta;
1189 
1190  break;
1191  }
1192 
1193  case 3:
1194  {
1195  if (_calculate_xyz)
1196  _ad_q_points[p].zero();
1197  _ad_dxyzdxi_map[p].zero();
1198  _ad_dxyzdeta_map[p].zero();
1199  _ad_dxyzdzeta_map[p].zero();
1200 
1201  for (std::size_t i = 0; i < num_shapes; i++)
1202  {
1203  libmesh_assert(elem_nodes[i]);
1204  const Node & node = *elem_nodes[i];
1205  libMesh::VectorValue<ADReal> elem_point = node;
1206  if (do_derivatives)
1207  for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1208  if (node.n_dofs(sys_num, disp_num))
1210  elem_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1211 
1212  _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
1213  _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
1214  _ad_dxyzdzeta_map[p].add_scaled(elem_point, dphidzeta_map[i][p]);
1215 
1216  if (_calculate_xyz)
1217  _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
1218  }
1219 
1220  const auto &dx_dxi = _ad_dxyzdxi_map[p](0), &dy_dxi = _ad_dxyzdxi_map[p](1),
1221  &dz_dxi = _ad_dxyzdxi_map[p](2), &dx_deta = _ad_dxyzdeta_map[p](0),
1222  &dy_deta = _ad_dxyzdeta_map[p](1), &dz_deta = _ad_dxyzdeta_map[p](2),
1223  &dx_dzeta = _ad_dxyzdzeta_map[p](0), &dy_dzeta = _ad_dxyzdzeta_map[p](1),
1224  &dz_dzeta = _ad_dxyzdzeta_map[p](2);
1225 
1226  _ad_jac[p] = (dx_dxi * (dy_deta * dz_dzeta - dz_deta * dy_dzeta) +
1227  dy_dxi * (dz_deta * dx_dzeta - dx_deta * dz_dzeta) +
1228  dz_dxi * (dx_deta * dy_dzeta - dy_deta * dx_dzeta));
1229 
1230  if (_ad_jac[p].value() <= -TOLERANCE * TOLERANCE)
1231  {
1232  static bool failing = false;
1233  if (!failing)
1234  {
1235  failing = true;
1237  libmesh_error_msg("ERROR: negative Jacobian " << _ad_jac[p].value() << " at point index "
1238  << p << " in element " << elem->id());
1239  }
1240  else
1241  return;
1242  }
1243 
1244  _ad_JxW[p] = _ad_jac[p] * qw[p];
1245 
1246  const auto inv_jac = 1. / _ad_jac[p];
1247 
1248  _ad_dxidx_map[p] = (dy_deta * dz_dzeta - dz_deta * dy_dzeta) * inv_jac;
1249  _ad_dxidy_map[p] = (dz_deta * dx_dzeta - dx_deta * dz_dzeta) * inv_jac;
1250  _ad_dxidz_map[p] = (dx_deta * dy_dzeta - dy_deta * dx_dzeta) * inv_jac;
1251 
1252  _ad_detadx_map[p] = (dz_dxi * dy_dzeta - dy_dxi * dz_dzeta) * inv_jac;
1253  _ad_detady_map[p] = (dx_dxi * dz_dzeta - dz_dxi * dx_dzeta) * inv_jac;
1254  _ad_detadz_map[p] = (dy_dxi * dx_dzeta - dx_dxi * dy_dzeta) * inv_jac;
1255 
1256  _ad_dzetadx_map[p] = (dy_dxi * dz_deta - dz_dxi * dy_deta) * inv_jac;
1257  _ad_dzetady_map[p] = (dz_dxi * dx_deta - dx_dxi * dz_deta) * inv_jac;
1258  _ad_dzetadz_map[p] = (dx_dxi * dy_deta - dy_dxi * dx_deta) * inv_jac;
1259 
1260  break;
1261  }
1262 
1263  default:
1264  libmesh_error_msg("Invalid dim = " << dim);
1265  }
1266 }
1267 
1268 void
1269 Assembly::reinitFEFace(const Elem * elem, unsigned int side)
1270 {
1271  unsigned int dim = elem->dim();
1272 
1273  for (const auto & it : _fe_face[dim])
1274  {
1275  FEBase & fe_face = *it.second;
1276  const FEType & fe_type = it.first;
1277  FEShapeData & fesd = *_fe_shape_data_face[fe_type];
1278  fe_face.reinit(elem, side);
1279  _current_fe_face[fe_type] = &fe_face;
1280 
1281  fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face.get_phi()));
1282  fesd._grad_phi.shallowCopy(
1283  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_dphi()));
1284  if (_need_second_derivative.count(fe_type))
1285  fesd._second_phi.shallowCopy(
1286  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
1287  }
1288  for (const auto & it : _vector_fe_face[dim])
1289  {
1290  FEVectorBase & fe_face = *it.second;
1291  const FEType & fe_type = it.first;
1292 
1293  _current_vector_fe_face[fe_type] = &fe_face;
1294 
1295  VectorFEShapeData & fesd = *_vector_fe_shape_data_face[fe_type];
1296 
1297  fe_face.reinit(elem, side);
1298 
1299  fesd._phi.shallowCopy(
1300  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_phi()));
1301  fesd._grad_phi.shallowCopy(
1302  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_dphi()));
1303  if (_need_second_derivative.count(fe_type))
1304  fesd._second_phi.shallowCopy(
1305  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
1306  if (_need_curl.count(fe_type))
1307  fesd._curl_phi.shallowCopy(
1308  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
1309  if (_need_face_div.count(fe_type))
1310  fesd._div_phi.shallowCopy(
1311  const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
1312  }
1313  if (!_unique_fe_face_helper.empty())
1314  {
1315  mooseAssert(dim < _unique_fe_face_helper.size(), "We should be in bounds here");
1316  _unique_fe_face_helper[dim]->reinit(elem, side);
1317  }
1318 
1319  // During that last loop the helper objects will have been reinitialized as well
1320  // We need to dig out the q_points and JxW from it.
1322  const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_xyz()));
1324  const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_JxW()));
1326  const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_normals()));
1327 
1328  _mapped_normals.resize(_current_normals.size(), Eigen::Map<RealDIMValue>(nullptr));
1329  for (unsigned int i = 0; i < _current_normals.size(); i++)
1330  // Note: this does NOT do any allocation. It is "reconstructing" the object in place
1331  new (&_mapped_normals[i]) Eigen::Map<RealDIMValue>(const_cast<Real *>(&_current_normals[i](0)));
1332 
1335  const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_curvatures()));
1336 
1337  computeADFace(*elem, side);
1338 
1339  if (_xfem != nullptr)
1341 
1342  auto n = numExtraElemIntegers();
1343  for (auto i : make_range(n))
1346 }
1347 
1348 void
1349 Assembly::computeFaceMap(const Elem & elem, const unsigned int side, const std::vector<Real> & qw)
1350 {
1351  // Important quantities calculated by this method:
1352  // - _ad_JxW_face
1353  // - _ad_q_points_face
1354  // - _ad_normals
1355  // - _ad_curvatures
1356 
1357  const Elem & side_elem = _compute_face_map_side_elem_builder(elem, side);
1358  const auto dim = elem.dim();
1359  const auto n_qp = qw.size();
1360  const auto & dpsidxi_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsidxi();
1361  const auto & dpsideta_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsideta();
1362  const auto & psi_map = _holder_fe_face_helper[dim]->get_fe_map().get_psi();
1363  std::vector<std::vector<Real>> const * d2psidxi2_map = nullptr;
1364  std::vector<std::vector<Real>> const * d2psidxideta_map = nullptr;
1365  std::vector<std::vector<Real>> const * d2psideta2_map = nullptr;
1366  const auto sys_num = _sys.number();
1367  const bool do_derivatives = ADReal::do_derivatives && sys_num == _subproblem.currentNlSysNum();
1368 
1370  {
1371  d2psidxi2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxi2();
1372  d2psidxideta_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxideta();
1373  d2psideta2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psideta2();
1374  }
1375 
1376  switch (dim)
1377  {
1378  case 1:
1379  {
1380  if (!n_qp)
1381  break;
1382 
1383  if (side_elem.node_id(0) == elem.node_id(0))
1384  _ad_normals[0] = Point(-1.);
1385  else
1386  _ad_normals[0] = Point(1.);
1387 
1388  VectorValue<ADReal> side_point;
1389  if (_calculate_face_xyz)
1390  {
1391  const Node & node = side_elem.node_ref(0);
1392  side_point = node;
1393 
1394  if (do_derivatives)
1395  for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1397  side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1398  }
1399 
1400  for (const auto p : make_range(n_qp))
1401  {
1402  if (_calculate_face_xyz)
1403  {
1404  _ad_q_points_face[p].zero();
1405  _ad_q_points_face[p].add_scaled(side_point, psi_map[0][p]);
1406  }
1407 
1408  _ad_normals[p] = _ad_normals[0];
1409  _ad_JxW_face[p] = 1.0 * qw[p];
1410  }
1411 
1412  break;
1413  }
1414 
1415  case 2:
1416  {
1417  _ad_dxyzdxi_map.resize(n_qp);
1419  _ad_d2xyzdxi2_map.resize(n_qp);
1420 
1421  for (const auto p : make_range(n_qp))
1422  _ad_dxyzdxi_map[p].zero();
1423  if (_calculate_face_xyz)
1424  for (const auto p : make_range(n_qp))
1425  _ad_q_points_face[p].zero();
1427  for (const auto p : make_range(n_qp))
1428  _ad_d2xyzdxi2_map[p].zero();
1429 
1430  const auto n_mapping_shape_functions =
1431  FE<2, LAGRANGE>::n_dofs(&side_elem, side_elem.default_order());
1432 
1433  for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1434  {
1435  const Node & node = side_elem.node_ref(i);
1436  VectorValue<ADReal> side_point = node;
1437 
1438  if (do_derivatives)
1439  for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1441  side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1442 
1443  for (const auto p : make_range(n_qp))
1444  _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1445  if (_calculate_face_xyz)
1446  for (const auto p : make_range(n_qp))
1447  _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1449  for (const auto p : make_range(n_qp))
1450  _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1451  }
1452 
1453  for (const auto p : make_range(n_qp))
1454  {
1455  _ad_normals[p] =
1456  (VectorValue<ADReal>(_ad_dxyzdxi_map[p](1), -_ad_dxyzdxi_map[p](0), 0.)).unit();
1457  const auto the_jac = _ad_dxyzdxi_map[p].norm();
1458  _ad_JxW_face[p] = the_jac * qw[p];
1460  {
1461  const auto numerator = _ad_d2xyzdxi2_map[p] * _ad_normals[p];
1462  const auto denominator = _ad_dxyzdxi_map[p].norm_sq();
1463  libmesh_assert_not_equal_to(denominator, 0);
1464  _ad_curvatures[p] = numerator / denominator;
1465  }
1466  }
1467 
1468  break;
1469  }
1470 
1471  case 3:
1472  {
1473  _ad_dxyzdxi_map.resize(n_qp);
1474  _ad_dxyzdeta_map.resize(n_qp);
1476  {
1477  _ad_d2xyzdxi2_map.resize(n_qp);
1478  _ad_d2xyzdxideta_map.resize(n_qp);
1479  _ad_d2xyzdeta2_map.resize(n_qp);
1480  }
1481 
1482  for (const auto p : make_range(n_qp))
1483  {
1484  _ad_dxyzdxi_map[p].zero();
1485  _ad_dxyzdeta_map[p].zero();
1486  }
1487  if (_calculate_face_xyz)
1488  for (const auto p : make_range(n_qp))
1489  _ad_q_points_face[p].zero();
1491  for (const auto p : make_range(n_qp))
1492  {
1493  _ad_d2xyzdxi2_map[p].zero();
1494  _ad_d2xyzdxideta_map[p].zero();
1495  _ad_d2xyzdeta2_map[p].zero();
1496  }
1497 
1498  const unsigned int n_mapping_shape_functions =
1499  FE<3, LAGRANGE>::n_dofs(&side_elem, side_elem.default_order());
1500 
1501  for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1502  {
1503  const Node & node = side_elem.node_ref(i);
1504  VectorValue<ADReal> side_point = node;
1505 
1506  if (do_derivatives)
1507  for (const auto & [disp_num, direction] : _disp_numbers_and_directions)
1509  side_point(direction).derivatives(), node.dof_number(sys_num, disp_num, 0), 1.);
1510 
1511  for (const auto p : make_range(n_qp))
1512  {
1513  _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1514  _ad_dxyzdeta_map[p].add_scaled(side_point, dpsideta_map[i][p]);
1515  }
1516  if (_calculate_face_xyz)
1517  for (const auto p : make_range(n_qp))
1518  _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1520  for (const auto p : make_range(n_qp))
1521  {
1522  _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1523  _ad_d2xyzdxideta_map[p].add_scaled(side_point, (*d2psidxideta_map)[i][p]);
1524  _ad_d2xyzdeta2_map[p].add_scaled(side_point, (*d2psideta2_map)[i][p]);
1525  }
1526  }
1527 
1528  for (const auto p : make_range(n_qp))
1529  {
1530  _ad_normals[p] = _ad_dxyzdxi_map[p].cross(_ad_dxyzdeta_map[p]).unit();
1531 
1532  const auto &dxdxi = _ad_dxyzdxi_map[p](0), &dxdeta = _ad_dxyzdeta_map[p](0),
1533  &dydxi = _ad_dxyzdxi_map[p](1), &dydeta = _ad_dxyzdeta_map[p](1),
1534  &dzdxi = _ad_dxyzdxi_map[p](2), &dzdeta = _ad_dxyzdeta_map[p](2);
1535 
1536  const auto g11 = (dxdxi * dxdxi + dydxi * dydxi + dzdxi * dzdxi);
1537 
1538  const auto g12 = (dxdxi * dxdeta + dydxi * dydeta + dzdxi * dzdeta);
1539 
1540  const auto & g21 = g12;
1541 
1542  const auto g22 = (dxdeta * dxdeta + dydeta * dydeta + dzdeta * dzdeta);
1543 
1544  const auto the_jac = std::sqrt(g11 * g22 - g12 * g21);
1545 
1546  _ad_JxW_face[p] = the_jac * qw[p];
1547 
1549  {
1550  const auto L = -_ad_d2xyzdxi2_map[p] * _ad_normals[p];
1551  const auto M = -_ad_d2xyzdxideta_map[p] * _ad_normals[p];
1552  const auto N = -_ad_d2xyzdeta2_map[p] * _ad_normals[p];
1553  const auto E = _ad_dxyzdxi_map[p].norm_sq();
1554  const auto F = _ad_dxyzdxi_map[p] * _ad_dxyzdeta_map[p];
1555  const auto G = _ad_dxyzdeta_map[p].norm_sq();
1556 
1557  const auto numerator = E * N - 2. * F * M + G * L;
1558  const auto denominator = E * G - F * F;
1559  libmesh_assert_not_equal_to(denominator, 0.);
1560  _ad_curvatures[p] = 0.5 * numerator / denominator;
1561  }
1562  }
1563 
1564  break;
1565  }
1566 
1567  default:
1568  mooseError("Invalid dimension dim = ", dim);
1569  }
1570 }
1571 
1572 void
1573 Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1574 {
1575  unsigned int neighbor_dim = neighbor->dim();
1576 
1577  // reinit neighbor face
1578  for (const auto & it : _fe_face_neighbor[neighbor_dim])
1579  {
1580  FEBase & fe_face_neighbor = *it.second;
1581  FEType fe_type = it.first;
1582  FEShapeData & fesd = *_fe_shape_data_face_neighbor[fe_type];
1583 
1584  fe_face_neighbor.reinit(neighbor, &reference_points);
1585 
1586  _current_fe_face_neighbor[fe_type] = &fe_face_neighbor;
1587 
1588  fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_phi()));
1589  fesd._grad_phi.shallowCopy(
1590  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor.get_dphi()));
1591  if (_need_second_derivative_neighbor.count(fe_type))
1592  fesd._second_phi.shallowCopy(
1593  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
1594  }
1595  for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
1596  {
1597  FEVectorBase & fe_face_neighbor = *it.second;
1598  const FEType & fe_type = it.first;
1599 
1600  _current_vector_fe_face_neighbor[fe_type] = &fe_face_neighbor;
1601 
1603 
1604  fe_face_neighbor.reinit(neighbor, &reference_points);
1605 
1606  fesd._phi.shallowCopy(
1607  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor.get_phi()));
1608  fesd._grad_phi.shallowCopy(
1609  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_dphi()));
1610  if (_need_second_derivative.count(fe_type))
1611  fesd._second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
1612  fe_face_neighbor.get_d2phi()));
1613  if (_need_curl.count(fe_type))
1614  fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
1615  fe_face_neighbor.get_curl_phi()));
1616  if (_need_face_neighbor_div.count(fe_type))
1617  fesd._div_phi.shallowCopy(
1618  const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
1619  }
1620  if (!_unique_fe_face_neighbor_helper.empty())
1621  {
1622  mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
1623  "We should be in bounds here");
1624  _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
1625  }
1626 
1628  const_cast<std::vector<Point> &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz()));
1629 }
1630 
1631 void
1632 Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1633 {
1634  unsigned int neighbor_dim = neighbor->dim();
1635 
1636  // reinit neighbor face
1637  for (const auto & it : _fe_neighbor[neighbor_dim])
1638  {
1639  FEBase & fe_neighbor = *it.second;
1640  FEType fe_type = it.first;
1641  FEShapeData & fesd = *_fe_shape_data_neighbor[fe_type];
1642 
1643  fe_neighbor.reinit(neighbor, &reference_points);
1644 
1645  _current_fe_neighbor[fe_type] = &fe_neighbor;
1646 
1647  fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_phi()));
1648  fesd._grad_phi.shallowCopy(
1649  const_cast<std::vector<std::vector<RealGradient>> &>(fe_neighbor.get_dphi()));
1650  if (_need_second_derivative_neighbor.count(fe_type))
1651  fesd._second_phi.shallowCopy(
1652  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor.get_d2phi()));
1653  }
1654  for (const auto & it : _vector_fe_neighbor[neighbor_dim])
1655  {
1656  FEVectorBase & fe_neighbor = *it.second;
1657  const FEType & fe_type = it.first;
1658 
1659  _current_vector_fe_neighbor[fe_type] = &fe_neighbor;
1660 
1662 
1663  fe_neighbor.reinit(neighbor, &reference_points);
1664 
1665  fesd._phi.shallowCopy(
1666  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_phi()));
1667  fesd._grad_phi.shallowCopy(
1668  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor.get_dphi()));
1669  if (_need_second_derivative.count(fe_type))
1670  fesd._second_phi.shallowCopy(
1671  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_neighbor.get_d2phi()));
1672  if (_need_curl.count(fe_type))
1673  fesd._curl_phi.shallowCopy(
1674  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_curl_phi()));
1675  if (_need_neighbor_div.count(fe_type))
1676  fesd._div_phi.shallowCopy(
1677  const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_div_phi()));
1678  }
1679  if (!_unique_fe_neighbor_helper.empty())
1680  {
1681  mooseAssert(neighbor_dim < _unique_fe_neighbor_helper.size(), "We should be in bounds here");
1682  _unique_fe_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
1683  }
1684 }
1685 
1686 void
1687 Assembly::reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1688 {
1689  unsigned int neighbor_dim = neighbor->dim();
1691  "Neighbor subdomain ID has not been correctly set");
1692 
1693  ArbitraryQuadrature * neighbor_rule =
1694  qrules(neighbor_dim, _current_neighbor_subdomain_id).neighbor.get();
1695  neighbor_rule->setPoints(reference_points);
1696  setNeighborQRule(neighbor_rule, neighbor_dim);
1697 
1700  "current neighbor subdomain has been set incorrectly");
1701 
1702  // Calculate the volume of the neighbor
1704  {
1705  unsigned int dim = neighbor->dim();
1707  QBase * qrule = qrules(dim).vol.get();
1708 
1709  fe.attach_quadrature_rule(qrule);
1710  fe.reinit(neighbor);
1711 
1712  const std::vector<Real> & JxW = fe.get_JxW();
1713  MooseArray<Point> q_points;
1714  q_points.shallowCopy(const_cast<std::vector<Point> &>(fe.get_xyz()));
1715 
1717 
1719  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1721  }
1722 
1723  auto n = numExtraElemIntegers();
1724  for (auto i : make_range(n))
1727 }
1728 
1729 template <typename Points, typename Coords>
1730 void
1732  const Points & q_points,
1733  Coords & coord,
1734  SubdomainID sub_id)
1735 {
1736 
1737  mooseAssert(qrule, "The quadrature rule is null in Assembly::setCoordinateTransformation");
1738  auto n_points = qrule->n_points();
1739  mooseAssert(n_points == q_points.size(),
1740  "The number of points in the quadrature rule doesn't match the number of passed-in "
1741  "points in Assembly::setCoordinateTransformation");
1742 
1743  // Make sure to honor the name of this method and set the _coord_type member because users may
1744  // make use of the const Moose::CoordinateSystem & coordTransformation() { return _coord_type; }
1745  // API. MaterialBase for example uses it
1747 
1748  coord.resize(n_points);
1749  for (unsigned int qp = 0; qp < n_points; qp++)
1750  coordTransformFactor(_subproblem, sub_id, q_points[qp], coord[qp]);
1751 }
1752 
1753 void
1755 {
1757  return;
1758 
1761  if (_calculate_ad_coord)
1764 
1765  _current_elem_volume = 0.;
1766  for (unsigned int qp = 0; qp < _current_qrule->n_points(); qp++)
1768 
1770 }
1771 
1772 void
1774 {
1776  return;
1777 
1780  if (_calculate_ad_coord)
1783 
1784  _current_side_volume = 0.;
1785  for (unsigned int qp = 0; qp < _current_qrule_face->n_points(); qp++)
1787 
1789 }
1790 
1791 void
1792 Assembly::reinitAtPhysical(const Elem * elem, const std::vector<Point> & physical_points)
1793 {
1794  _current_elem = elem;
1795  _current_neighbor_elem = nullptr;
1797  "current subdomain has been set incorrectly");
1799 
1800  FEMap::inverse_map(elem->dim(), elem, physical_points, _temp_reference_points);
1801 
1803 
1804  // Save off the physical points
1805  _current_physical_points = physical_points;
1806 }
1807 
1808 void
1809 Assembly::setVolumeQRule(const Elem * const elem)
1810 {
1811  unsigned int elem_dimension = elem->dim();
1812  _current_qrule_volume = qrules(elem_dimension).vol.get();
1813  // Make sure the qrule is the right one
1815  setVolumeQRule(_current_qrule_volume, elem_dimension);
1816 }
1817 
1818 void
1819 Assembly::reinit(const Elem * elem)
1820 {
1821  _current_elem = elem;
1822  _current_neighbor_elem = nullptr;
1824  "current subdomain has been set incorrectly");
1826  if (_sc)
1828 
1830  reinitFE(elem);
1831 
1833 }
1834 
1835 void
1836 Assembly::reinit(const Elem * elem, const std::vector<Point> & reference_points)
1837 {
1838  _current_elem = elem;
1839  _current_neighbor_elem = nullptr;
1841  "current subdomain has been set incorrectly");
1843 
1844  unsigned int elem_dimension = _current_elem->dim();
1845 
1846  _current_qrule_arbitrary = qrules(elem_dimension).arbitrary_vol.get();
1847 
1848  // Make sure the qrule is the right one
1850  setVolumeQRule(_current_qrule_arbitrary, elem_dimension);
1851 
1852  _current_qrule_arbitrary->setPoints(reference_points);
1853 
1854  reinitFE(elem);
1855 
1857 }
1858 
1859 void
1861 {
1862  _current_elem = &fi.elem();
1864  _current_side = fi.elemSideID();
1867  "current subdomain has been set incorrectly");
1868 
1871 
1872  prepareResidual();
1873  prepareNeighbor();
1875 
1876  unsigned int dim = _current_elem->dim();
1877  if (_current_qrule_face != qrules(dim).fv_face.get())
1878  {
1879  setFaceQRule(qrules(dim).fv_face.get(), dim);
1880  // The order of the element that is used for initing here doesn't matter since this will just
1881  // be used for constant monomials (which only need a single integration point)
1882  if (dim == 3)
1883  _current_qrule_face->init(QUAD4, /* p_level = */ 0, /* simple_type_only = */ true);
1884  else
1885  _current_qrule_face->init(EDGE2, /* p_level = */ 0, /* simple_type_only = */ true);
1886  }
1887 
1889 
1890  mooseAssert(_current_qrule_face->n_points() == 1,
1891  "Our finite volume quadrature rule should always yield a single point");
1892 
1893  // We've initialized the reference points. Now we need to compute the physical location of the
1894  // quadrature points. We do not do any FE initialization so we cannot simply copy over FE
1895  // results like we do in reinitFEFace. Instead we handle the computation of the physical
1896  // locations manually
1898  const auto & ref_points = _current_qrule_face->get_points();
1899  const auto & ref_point = ref_points[0];
1900  auto physical_point = FEMap::map(_current_side_elem->dim(), _current_side_elem, ref_point);
1901  _current_q_points_face[0] = physical_point;
1902 
1904  {
1906  "current neighbor subdomain has been set incorrectly");
1907  // Now handle the neighbor qrule/qpoints
1908  ArbitraryQuadrature * const neighbor_rule =
1910  // Here we are setting a reference point that is correct for the neighbor *side* element. It
1911  // would be wrong if this reference point is used for a volumetric FE reinit with the neighbor
1912  neighbor_rule->setPoints(ref_points);
1913  setNeighborQRule(neighbor_rule, _current_neighbor_elem->dim());
1915  _current_q_points_face_neighbor[0] = std::move(physical_point);
1916  }
1917 }
1918 
1919 QBase *
1920 Assembly::qruleFace(const Elem * elem, unsigned int side)
1921 {
1922  return qruleFaceHelper<QBase>(elem, side, [](QRules & q) { return q.face.get(); });
1923 }
1924 
1926 Assembly::qruleArbitraryFace(const Elem * elem, unsigned int side)
1927 {
1928  return qruleFaceHelper<ArbitraryQuadrature>(
1929  elem, side, [](QRules & q) { return q.arbitrary_face.get(); });
1930 }
1931 
1932 void
1933 Assembly::setFaceQRule(const Elem * const elem, const unsigned int side)
1934 {
1935  const auto elem_dimension = elem->dim();
1937  auto rule = qruleFace(elem, side);
1938  if (_current_qrule_face != rule)
1939  setFaceQRule(rule, elem_dimension);
1940 }
1941 
1942 void
1943 Assembly::reinit(const Elem * const elem, const unsigned int side)
1944 {
1945  _current_elem = elem;
1946  _current_neighbor_elem = nullptr;
1948  "current subdomain has been set incorrectly");
1949  _current_side = side;
1952 
1954 
1957 
1959 }
1960 
1961 void
1962 Assembly::reinit(const Elem * elem, unsigned int side, const std::vector<Point> & reference_points)
1963 {
1964  _current_elem = elem;
1965  _current_neighbor_elem = nullptr;
1967  "current subdomain has been set incorrectly");
1968  _current_side = side;
1971 
1972  unsigned int elem_dimension = _current_elem->dim();
1973 
1975 
1976  // Make sure the qrule is the right one
1979 
1980  _current_qrule_arbitrary->setPoints(reference_points);
1981 
1983 
1985 
1987 }
1988 
1989 void
1990 Assembly::reinit(const Node * node)
1991 {
1992  _current_node = node;
1993  _current_neighbor_node = NULL;
1994 }
1995 
1996 void
1998  unsigned int side,
1999  const Elem * neighbor,
2000  unsigned int neighbor_side,
2001  const std::vector<Point> * neighbor_reference_points)
2002 {
2003  _current_neighbor_side = neighbor_side;
2004 
2005  reinit(elem, side);
2006 
2007  unsigned int neighbor_dim = neighbor->dim();
2008 
2009  if (neighbor_reference_points)
2010  _current_neighbor_ref_points = *neighbor_reference_points;
2011  else
2012  FEMap::inverse_map(
2014 
2016 
2019 }
2020 
2021 void
2023  unsigned int elem_side,
2024  Real tolerance,
2025  const std::vector<Point> * const pts,
2026  const std::vector<Real> * const weights)
2027 {
2028  _current_elem = elem;
2029 
2030  unsigned int elem_dim = elem->dim();
2031 
2032  // Attach the quadrature rules
2033  if (pts)
2034  {
2035  auto face_rule = qruleArbitraryFace(elem, elem_side);
2036  face_rule->setPoints(*pts);
2037  setFaceQRule(face_rule, elem_dim);
2038  }
2039  else
2040  {
2041  auto rule = qruleFace(elem, elem_side);
2042  if (_current_qrule_face != rule)
2043  setFaceQRule(rule, elem_dim);
2044  }
2045 
2046  // reinit face
2047  for (const auto & it : _fe_face[elem_dim])
2048  {
2049  FEBase & fe_face = *it.second;
2050  FEType fe_type = it.first;
2051  FEShapeData & fesd = *_fe_shape_data_face[fe_type];
2052 
2053  fe_face.reinit(elem, elem_side, tolerance, pts, weights);
2054 
2055  _current_fe_face[fe_type] = &fe_face;
2056 
2057  fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face.get_phi()));
2058  fesd._grad_phi.shallowCopy(
2059  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face.get_dphi()));
2060  if (_need_second_derivative_neighbor.count(fe_type))
2061  fesd._second_phi.shallowCopy(
2062  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
2063  }
2064  for (const auto & it : _vector_fe_face[elem_dim])
2065  {
2066  FEVectorBase & fe_face = *it.second;
2067  const FEType & fe_type = it.first;
2068 
2069  _current_vector_fe_face[fe_type] = &fe_face;
2070 
2071  VectorFEShapeData & fesd = *_vector_fe_shape_data_face[fe_type];
2072 
2073  fe_face.reinit(elem, elem_side, tolerance, pts, weights);
2074 
2075  fesd._phi.shallowCopy(
2076  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_phi()));
2077  fesd._grad_phi.shallowCopy(
2078  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_dphi()));
2079  if (_need_second_derivative.count(fe_type))
2080  fesd._second_phi.shallowCopy(
2081  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
2082  if (_need_curl.count(fe_type))
2083  fesd._curl_phi.shallowCopy(
2084  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
2085  if (_need_face_div.count(fe_type))
2086  fesd._div_phi.shallowCopy(
2087  const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
2088  }
2089  if (!_unique_fe_face_helper.empty())
2090  {
2091  mooseAssert(elem_dim < _unique_fe_face_helper.size(), "We should be in bounds here");
2092  _unique_fe_face_helper[elem_dim]->reinit(elem, elem_side, tolerance, pts, weights);
2093  }
2094 
2095  // During that last loop the helper objects will have been reinitialized
2097  const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_xyz()));
2099  const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_normals()));
2100  _current_tangents.shallowCopy(const_cast<std::vector<std::vector<Point>> &>(
2101  _holder_fe_face_helper[elem_dim]->get_tangents()));
2102  // Note that if the user did pass in points and not weights to this method, JxW will be garbage
2103  // and should not be used
2105  const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_JxW()));
2108  const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_curvatures()));
2109 
2110  computeADFace(*elem, elem_side);
2111 }
2112 
2113 void
2114 Assembly::computeADFace(const Elem & elem, const unsigned int side)
2115 {
2116  const auto dim = elem.dim();
2117 
2118  if (_subproblem.haveADObjects())
2119  {
2120  auto n_qp = _current_qrule_face->n_points();
2121  resizeADMappingObjects(n_qp, dim);
2122  _ad_normals.resize(n_qp);
2123  _ad_JxW_face.resize(n_qp);
2124  if (_calculate_face_xyz)
2125  _ad_q_points_face.resize(n_qp);
2127  _ad_curvatures.resize(n_qp);
2128 
2129  if (_displaced)
2130  {
2131  const auto & qw = _current_qrule_face->get_weights();
2132  computeFaceMap(elem, side, qw);
2133  const std::vector<Real> dummy_qw(n_qp, 1.);
2134 
2135  for (unsigned int qp = 0; qp != n_qp; qp++)
2137  }
2138  else
2139  {
2140  for (unsigned qp = 0; qp < n_qp; ++qp)
2141  {
2142  _ad_JxW_face[qp] = _current_JxW_face[qp];
2143  _ad_normals[qp] = _current_normals[qp];
2144  }
2145  if (_calculate_face_xyz)
2146  for (unsigned qp = 0; qp < n_qp; ++qp)
2149  for (unsigned qp = 0; qp < n_qp; ++qp)
2150  _ad_curvatures[qp] = _curvatures[qp];
2151  }
2152 
2153  for (const auto & it : _fe_face[dim])
2154  {
2155  FEBase & fe = *it.second;
2156  auto fe_type = it.first;
2157  auto num_shapes = FEInterface::n_shape_functions(fe_type, &elem);
2158  auto & grad_phi = _ad_grad_phi_data_face[fe_type];
2159 
2160  grad_phi.resize(num_shapes);
2161  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2162  grad_phi[i].resize(n_qp);
2163 
2164  const auto & regular_grad_phi = _fe_shape_data_face[fe_type]->_grad_phi;
2165 
2166  if (_displaced)
2167  computeGradPhiAD(&elem, n_qp, grad_phi, &fe);
2168  else
2169  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2170  for (unsigned qp = 0; qp < n_qp; ++qp)
2171  grad_phi[i][qp] = regular_grad_phi[i][qp];
2172  }
2173  for (const auto & it : _vector_fe_face[dim])
2174  {
2175  FEVectorBase & fe = *it.second;
2176  auto fe_type = it.first;
2177  auto num_shapes = FEInterface::n_shape_functions(fe_type, &elem);
2178  auto & grad_phi = _ad_vector_grad_phi_data_face[fe_type];
2179 
2180  grad_phi.resize(num_shapes);
2181  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2182  grad_phi[i].resize(n_qp);
2183 
2184  const auto & regular_grad_phi = _vector_fe_shape_data_face[fe_type]->_grad_phi;
2185 
2186  if (_displaced)
2187  computeGradPhiAD(&elem, n_qp, grad_phi, &fe);
2188  else
2189  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2190  for (unsigned qp = 0; qp < n_qp; ++qp)
2191  grad_phi[i][qp] = regular_grad_phi[i][qp];
2192  }
2193  }
2194 }
2195 
2196 void
2198  unsigned int neighbor_side,
2199  Real tolerance,
2200  const std::vector<Point> * const pts,
2201  const std::vector<Real> * const weights)
2202 {
2204 
2205  unsigned int neighbor_dim = neighbor->dim();
2206 
2207  ArbitraryQuadrature * neighbor_rule =
2208  qrules(neighbor_dim, neighbor->subdomain_id()).neighbor.get();
2209  neighbor_rule->setPoints(*pts);
2210 
2211  // Attach this quadrature rule to all the _fe_face_neighbor FE objects. This
2212  // has to have garbage quadrature weights but that's ok because we never
2213  // actually use the JxW coming from these FE reinit'd objects, e.g. we use the
2214  // JxW coming from the element face reinit for DGKernels or we use the JxW
2215  // coming from reinit of the mortar segment element in the case of mortar
2216  setNeighborQRule(neighbor_rule, neighbor_dim);
2217 
2218  // reinit neighbor face
2219  for (const auto & it : _fe_face_neighbor[neighbor_dim])
2220  {
2221  FEBase & fe_face_neighbor = *it.second;
2222  FEType fe_type = it.first;
2223  FEShapeData & fesd = *_fe_shape_data_face_neighbor[fe_type];
2224 
2225  fe_face_neighbor.reinit(neighbor, neighbor_side, tolerance, pts, weights);
2226 
2227  _current_fe_face_neighbor[fe_type] = &fe_face_neighbor;
2228 
2229  fesd._phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_phi()));
2230  fesd._grad_phi.shallowCopy(
2231  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor.get_dphi()));
2232  if (_need_second_derivative_neighbor.count(fe_type))
2233  fesd._second_phi.shallowCopy(
2234  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
2235  }
2236  for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
2237  {
2238  FEVectorBase & fe_face_neighbor = *it.second;
2239  const FEType & fe_type = it.first;
2240 
2241  _current_vector_fe_face_neighbor[fe_type] = &fe_face_neighbor;
2242 
2244 
2245  fe_face_neighbor.reinit(neighbor, neighbor_side, tolerance, pts, weights);
2246 
2247  fesd._phi.shallowCopy(
2248  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor.get_phi()));
2249  fesd._grad_phi.shallowCopy(
2250  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_dphi()));
2251  if (_need_second_derivative.count(fe_type))
2252  fesd._second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
2253  fe_face_neighbor.get_d2phi()));
2254  if (_need_curl.count(fe_type))
2255  fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
2256  fe_face_neighbor.get_curl_phi()));
2257  if (_need_face_neighbor_div.count(fe_type))
2258  fesd._div_phi.shallowCopy(
2259  const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
2260  }
2261  if (!_unique_fe_face_neighbor_helper.empty())
2262  {
2263  mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
2264  "We should be in bounds here");
2265  _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(
2266  neighbor, neighbor_side, tolerance, pts, weights);
2267  }
2268  // During that last loop the helper objects will have been reinitialized as well
2269  // We need to dig out the q_points from it
2271  const_cast<std::vector<Point> &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz()));
2272 }
2273 
2274 void
2276  const std::vector<Point> & pts,
2277  const std::vector<Real> & JxW)
2278 {
2279  const unsigned int elem_dim = elem->dim();
2280  mooseAssert(elem_dim == _mesh_dimension - 1,
2281  "Dual shape functions should only be computed on lower dimensional face elements");
2282 
2283  for (const auto & it : _fe_lower[elem_dim])
2284  {
2285  FEBase & fe_lower = *it.second;
2286  // We use customized quadrature rule for integration along the mortar segment elements
2287  fe_lower.set_calculate_default_dual_coeff(false);
2288  fe_lower.reinit_dual_shape_coeffs(elem, pts, JxW);
2289  }
2290 }
2291 
2292 void
2294  const std::vector<Point> * const pts,
2295  const std::vector<Real> * const weights)
2296 {
2298 
2299  const unsigned int elem_dim = elem->dim();
2300  mooseAssert(elem_dim < _mesh_dimension,
2301  "The lower dimensional element should truly be a lower dimensional element");
2302 
2303  if (pts)
2304  {
2305  // Lower rule matches the face rule for the higher dimensional element
2306  ArbitraryQuadrature * lower_rule = qrules(elem_dim + 1).arbitrary_face.get();
2307 
2308  // This also sets the quadrature weights to unity
2309  lower_rule->setPoints(*pts);
2310 
2311  if (weights)
2312  lower_rule->setWeights(*weights);
2313 
2314  setLowerQRule(lower_rule, elem_dim);
2315  }
2316  else if (_current_qrule_lower != qrules(elem_dim + 1).face.get())
2317  setLowerQRule(qrules(elem_dim + 1).face.get(), elem_dim);
2318 
2319  for (const auto & it : _fe_lower[elem_dim])
2320  {
2321  FEBase & fe_lower = *it.second;
2322  FEType fe_type = it.first;
2323 
2324  fe_lower.reinit(elem);
2325 
2326  if (FEShapeData * fesd = _fe_shape_data_lower[fe_type].get())
2327  {
2328  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower.get_phi()));
2329  fesd->_grad_phi.shallowCopy(
2330  const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower.get_dphi()));
2331  if (_need_second_derivative_neighbor.count(fe_type))
2332  fesd->_second_phi.shallowCopy(
2333  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_d2phi()));
2334  }
2335 
2336  // Dual shape functions need to be computed after primal basis being initialized
2337  if (FEShapeData * fesd = _fe_shape_data_dual_lower[fe_type].get())
2338  {
2339  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower.get_dual_phi()));
2340  fesd->_grad_phi.shallowCopy(
2341  const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower.get_dual_dphi()));
2342  if (_need_second_derivative_neighbor.count(fe_type))
2343  fesd->_second_phi.shallowCopy(
2344  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_dual_d2phi()));
2345  }
2346  }
2347  if (!_unique_fe_lower_helper.empty())
2348  {
2349  mooseAssert(elem_dim < _unique_fe_lower_helper.size(), "We should be in bounds here");
2350  _unique_fe_lower_helper[elem_dim]->reinit(elem);
2351  }
2352 
2354  return;
2355 
2356  if (pts && !weights)
2357  {
2358  // We only have dummy weights so the JxWs computed during our FE reinits are meaningless and
2359  // we cannot use them
2360 
2362  // We are in a Cartesian coordinate system and we can just use the element volume method
2363  // which has fast computation for certain element types
2365  else
2366  // We manually compute the volume taking the curvilinear coordinate transformations into
2367  // account
2369  }
2370  else
2371  {
2372  // During that last loop the helper objects will have been reinitialized as well
2373  FEBase & helper_fe = *_holder_fe_lower_helper[elem_dim];
2374  const auto & physical_q_points = helper_fe.get_xyz();
2375  const auto & JxW = helper_fe.get_JxW();
2376  MooseArray<Real> coord;
2378  _current_qrule_lower, physical_q_points, coord, elem->subdomain_id());
2380  for (const auto qp : make_range(_current_qrule_lower->n_points()))
2381  _current_lower_d_elem_volume += JxW[qp] * coord[qp];
2382  }
2383 }
2384 
2385 void
2387 {
2388  mooseAssert(elem->dim() < _mesh_dimension,
2389  "You should be calling reinitNeighborLowerDElem on a lower dimensional element");
2390 
2392 
2394  return;
2395 
2397  // We are in a Cartesian coordinate system and we can just use the element volume method which
2398  // has fast computation for certain element types
2400  else
2401  // We manually compute the volume taking the curvilinear coordinate transformations into
2402  // account
2404 }
2405 
2406 void
2408 {
2409  mooseAssert(elem->dim() == _mesh_dimension - 1,
2410  "You should be calling reinitMortarElem on a lower dimensional element");
2411 
2412  _fe_msm->reinit(elem);
2413  _msm_elem = elem;
2414 
2415  MooseArray<Point> array_q_points;
2416  array_q_points.shallowCopy(const_cast<std::vector<Point> &>(_fe_msm->get_xyz()));
2418 }
2419 
2420 void
2421 Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
2422  unsigned int neighbor_side,
2423  const std::vector<Point> & physical_points)
2424 {
2425  unsigned int neighbor_dim = neighbor->dim();
2426  FEMap::inverse_map(neighbor_dim, neighbor, physical_points, _current_neighbor_ref_points);
2427 
2428  if (_need_JxW_neighbor)
2429  {
2430  mooseAssert(
2431  physical_points.size() == 1,
2432  "If reinitializing with more than one point, then I am dubious of your use case. Perhaps "
2433  "you are performing a DG type method and you are reinitializing using points from the "
2434  "element face. In such a case your neighbor JxW must have its index order 'match' the "
2435  "element JxW index order, e.g. imagining a vertical 1D face with two quadrature points, "
2436  "if "
2437  "index 0 for elem JxW corresponds to the 'top' quadrature point, then index 0 for "
2438  "neighbor "
2439  "JxW must also correspond to the 'top' quadrature point. And libMesh/MOOSE has no way to "
2440  "guarantee that with multiple quadrature points.");
2441 
2443 
2444  // With a single point our size-1 JxW should just be the element volume
2447  }
2448 
2451 
2452  // Save off the physical points
2453  _current_physical_points = physical_points;
2454 }
2455 
2456 void
2457 Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
2458  const std::vector<Point> & physical_points)
2459 {
2460  unsigned int neighbor_dim = neighbor->dim();
2461  FEMap::inverse_map(neighbor_dim, neighbor, physical_points, _current_neighbor_ref_points);
2462 
2465  // Save off the physical points
2466  _current_physical_points = physical_points;
2467 }
2468 
2469 void
2471 {
2472  _cm = cm;
2473 
2474  unsigned int n_vars = _sys.nVariables();
2475 
2476  _cm_ss_entry.clear();
2477  _cm_sf_entry.clear();
2478  _cm_fs_entry.clear();
2479  _cm_ff_entry.clear();
2480 
2481  auto & vars = _sys.getVariables(_tid);
2482 
2483  _block_diagonal_matrix = true;
2484  for (auto & ivar : vars)
2485  {
2486  auto i = ivar->number();
2487  if (i >= _component_block_diagonal.size())
2488  _component_block_diagonal.resize(i + 1, true);
2489 
2490  auto ivar_start = _cm_ff_entry.size();
2491  for (unsigned int k = 0; k < ivar->count(); ++k)
2492  {
2493  unsigned int iv = i + k;
2494  for (const auto & j : ConstCouplingRow(iv, *_cm))
2495  {
2496  if (_sys.isScalarVariable(j))
2497  {
2498  auto & jvar = _sys.getScalarVariable(_tid, j);
2499  _cm_fs_entry.push_back(std::make_pair(ivar, &jvar));
2500  _block_diagonal_matrix = false;
2501  }
2502  else
2503  {
2504  auto & jvar = _sys.getVariable(_tid, j);
2505  auto pair = std::make_pair(ivar, &jvar);
2506  auto c = ivar_start;
2507  // check if the pair has been pushed or not
2508  bool has_pair = false;
2509  for (; c < _cm_ff_entry.size(); ++c)
2510  if (_cm_ff_entry[c] == pair)
2511  {
2512  has_pair = true;
2513  break;
2514  }
2515  if (!has_pair)
2516  _cm_ff_entry.push_back(pair);
2517  // only set having diagonal matrix to false when ivar and jvar numbers are different
2518  // Note: for array variables, since we save the entire local Jacobian of all components,
2519  // even there are couplings among components of the same array variable, we still
2520  // do not set the flag to false.
2521  if (i != jvar.number())
2522  _block_diagonal_matrix = false;
2523  else if (iv != j)
2524  _component_block_diagonal[i] = false;
2525  }
2526  }
2527  }
2528  }
2529 
2530  auto & scalar_vars = _sys.getScalarVariables(_tid);
2531 
2532  for (auto & ivar : scalar_vars)
2533  {
2534  auto i = ivar->number();
2535  if (i >= _component_block_diagonal.size())
2536  _component_block_diagonal.resize(i + 1, true);
2537 
2538  for (const auto & j : ConstCouplingRow(i, *_cm))
2539  if (_sys.isScalarVariable(j))
2540  {
2541  auto & jvar = _sys.getScalarVariable(_tid, j);
2542  _cm_ss_entry.push_back(std::make_pair(ivar, &jvar));
2543  }
2544  else
2545  {
2546  auto & jvar = _sys.getVariable(_tid, j);
2547  _cm_sf_entry.push_back(std::make_pair(ivar, &jvar));
2548  }
2549  }
2550 
2551  if (_block_diagonal_matrix && scalar_vars.size() != 0)
2552  _block_diagonal_matrix = false;
2553 
2554  auto num_vector_tags = _residual_vector_tags.size();
2555 
2556  _sub_Re.resize(num_vector_tags);
2557  _sub_Rn.resize(num_vector_tags);
2558  _sub_Rl.resize(num_vector_tags);
2559  for (MooseIndex(_sub_Re) i = 0; i < _sub_Re.size(); i++)
2560  {
2561  _sub_Re[i].resize(n_vars);
2562  _sub_Rn[i].resize(n_vars);
2563  _sub_Rl[i].resize(n_vars);
2564  }
2565 
2566  _cached_residual_values.resize(num_vector_tags);
2567  _cached_residual_rows.resize(num_vector_tags);
2568 
2569  auto num_matrix_tags = _subproblem.numMatrixTags();
2570 
2571  _cached_jacobian_values.resize(num_matrix_tags);
2572  _cached_jacobian_rows.resize(num_matrix_tags);
2573  _cached_jacobian_cols.resize(num_matrix_tags);
2574 
2575  // Element matrices
2576  _sub_Kee.resize(num_matrix_tags);
2577  _sub_Keg.resize(num_matrix_tags);
2578  _sub_Ken.resize(num_matrix_tags);
2579  _sub_Kne.resize(num_matrix_tags);
2580  _sub_Knn.resize(num_matrix_tags);
2581  _sub_Kll.resize(num_matrix_tags);
2582  _sub_Kle.resize(num_matrix_tags);
2583  _sub_Kln.resize(num_matrix_tags);
2584  _sub_Kel.resize(num_matrix_tags);
2585  _sub_Knl.resize(num_matrix_tags);
2586 
2587  _jacobian_block_used.resize(num_matrix_tags);
2588  _jacobian_block_neighbor_used.resize(num_matrix_tags);
2589  _jacobian_block_lower_used.resize(num_matrix_tags);
2590  _jacobian_block_nonlocal_used.resize(num_matrix_tags);
2591 
2592  for (MooseIndex(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
2593  {
2594  _sub_Keg[tag].resize(n_vars);
2595  _sub_Ken[tag].resize(n_vars);
2596  _sub_Kne[tag].resize(n_vars);
2597  _sub_Knn[tag].resize(n_vars);
2598  _sub_Kee[tag].resize(n_vars);
2599  _sub_Kll[tag].resize(n_vars);
2600  _sub_Kle[tag].resize(n_vars);
2601  _sub_Kln[tag].resize(n_vars);
2602  _sub_Kel[tag].resize(n_vars);
2603  _sub_Knl[tag].resize(n_vars);
2604 
2605  _jacobian_block_used[tag].resize(n_vars);
2607  _jacobian_block_lower_used[tag].resize(n_vars);
2609  for (MooseIndex(n_vars) i = 0; i < n_vars; ++i)
2610  {
2612  {
2613  _sub_Kee[tag][i].resize(n_vars);
2614  _sub_Keg[tag][i].resize(n_vars);
2615  _sub_Ken[tag][i].resize(n_vars);
2616  _sub_Kne[tag][i].resize(n_vars);
2617  _sub_Knn[tag][i].resize(n_vars);
2618  _sub_Kll[tag][i].resize(n_vars);
2619  _sub_Kle[tag][i].resize(n_vars);
2620  _sub_Kln[tag][i].resize(n_vars);
2621  _sub_Kel[tag][i].resize(n_vars);
2622  _sub_Knl[tag][i].resize(n_vars);
2623 
2624  _jacobian_block_used[tag][i].resize(n_vars);
2625  _jacobian_block_neighbor_used[tag][i].resize(n_vars);
2626  _jacobian_block_lower_used[tag][i].resize(n_vars);
2627  _jacobian_block_nonlocal_used[tag][i].resize(n_vars);
2628  }
2629  else
2630  {
2631  _sub_Kee[tag][i].resize(1);
2632  _sub_Keg[tag][i].resize(1);
2633  _sub_Ken[tag][i].resize(1);
2634  _sub_Kne[tag][i].resize(1);
2635  _sub_Knn[tag][i].resize(1);
2636  _sub_Kll[tag][i].resize(1);
2637  _sub_Kle[tag][i].resize(1);
2638  _sub_Kln[tag][i].resize(1);
2639  _sub_Kel[tag][i].resize(1);
2640  _sub_Knl[tag][i].resize(1);
2641 
2642  _jacobian_block_used[tag][i].resize(1);
2643  _jacobian_block_neighbor_used[tag][i].resize(1);
2644  _jacobian_block_lower_used[tag][i].resize(1);
2645  _jacobian_block_nonlocal_used[tag][i].resize(1);
2646  }
2647  }
2648  }
2649 }
2650 
2651 void
2653 {
2654  _cm_nonlocal_entry.clear();
2655 
2656  auto & vars = _sys.getVariables(_tid);
2657 
2658  for (auto & ivar : vars)
2659  {
2660  auto i = ivar->number();
2661  auto ivar_start = _cm_nonlocal_entry.size();
2662  for (unsigned int k = 0; k < ivar->count(); ++k)
2663  {
2664  unsigned int iv = i + k;
2665  for (const auto & j : ConstCouplingRow(iv, _nonlocal_cm))
2666  if (!_sys.isScalarVariable(j))
2667  {
2668  auto & jvar = _sys.getVariable(_tid, j);
2669  auto pair = std::make_pair(ivar, &jvar);
2670  auto c = ivar_start;
2671  // check if the pair has been pushed or not
2672  bool has_pair = false;
2673  for (; c < _cm_nonlocal_entry.size(); ++c)
2674  if (_cm_nonlocal_entry[c] == pair)
2675  {
2676  has_pair = true;
2677  break;
2678  }
2679  if (!has_pair)
2680  _cm_nonlocal_entry.push_back(pair);
2681  }
2682  }
2683  }
2684 }
2685 
2686 void
2688 {
2689  for (const auto & it : _cm_ff_entry)
2690  {
2691  MooseVariableFEBase & ivar = *(it.first);
2692  MooseVariableFEBase & jvar = *(it.second);
2693 
2694  unsigned int vi = ivar.number();
2695  unsigned int vj = jvar.number();
2696 
2697  unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2698 
2699  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2700  {
2701  jacobianBlock(vi, vj, LocalDataKey{}, tag)
2702  .resize(ivar.dofIndices().size() * ivar.count(), jvar.dofIndices().size() * jcount);
2703  jacobianBlockUsed(tag, vi, vj, false);
2704  }
2705  }
2706 }
2707 
2708 void
2710 {
2711  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2712  for (const auto & var : vars)
2713  for (auto & tag_Re : _sub_Re)
2714  tag_Re[var->number()].resize(var->dofIndices().size() * var->count());
2715 }
2716 
2717 void
2719 {
2721  prepareResidual();
2722 }
2723 
2724 void
2726 {
2727  for (const auto & it : _cm_nonlocal_entry)
2728  {
2729  MooseVariableFEBase & ivar = *(it.first);
2730  MooseVariableFEBase & jvar = *(it.second);
2731 
2732  unsigned int vi = ivar.number();
2733  unsigned int vj = jvar.number();
2734 
2735  unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2736 
2737  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2738  tag < _jacobian_block_nonlocal_used.size();
2739  tag++)
2740  {
2741  jacobianBlockNonlocal(vi, vj, LocalDataKey{}, tag)
2742  .resize(ivar.dofIndices().size() * ivar.count(), jvar.allDofIndices().size() * jcount);
2743  jacobianBlockNonlocalUsed(tag, vi, vj, false);
2744  }
2745  }
2746 }
2747 
2748 void
2750 {
2751  for (const auto & it : _cm_ff_entry)
2752  {
2753  MooseVariableFEBase & ivar = *(it.first);
2754  MooseVariableFEBase & jvar = *(it.second);
2755 
2756  unsigned int vi = ivar.number();
2757  unsigned int vj = jvar.number();
2758 
2759  unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2760 
2761  if (vi == var->number() || vj == var->number())
2762  {
2763  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2764  {
2765  jacobianBlock(vi, vj, LocalDataKey{}, tag)
2766  .resize(ivar.dofIndices().size() * ivar.count(), jvar.dofIndices().size() * jcount);
2767  jacobianBlockUsed(tag, vi, vj, false);
2768  }
2769  }
2770  }
2771 
2772  for (auto & tag_Re : _sub_Re)
2773  tag_Re[var->number()].resize(var->dofIndices().size() * var->count());
2774 }
2775 
2776 void
2778 {
2779  for (const auto & it : _cm_nonlocal_entry)
2780  {
2781  MooseVariableFEBase & ivar = *(it.first);
2782  MooseVariableFEBase & jvar = *(it.second);
2783 
2784  unsigned int vi = ivar.number();
2785  unsigned int vj = jvar.number();
2786 
2787  unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2788 
2789  if (vi == var->number() || vj == var->number())
2790  {
2791  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2792  tag < _jacobian_block_nonlocal_used.size();
2793  tag++)
2794  {
2795  jacobianBlockNonlocal(vi, vj, LocalDataKey{}, tag)
2796  .resize(ivar.dofIndices().size() * ivar.count(), jvar.allDofIndices().size() * jcount);
2797  jacobianBlockNonlocalUsed(tag, vi, vj);
2798  }
2799  }
2800  }
2801 }
2802 
2803 void
2805 {
2806  for (const auto & it : _cm_ff_entry)
2807  {
2808  MooseVariableFEBase & ivar = *(it.first);
2809  MooseVariableFEBase & jvar = *(it.second);
2810 
2811  unsigned int vi = ivar.number();
2812  unsigned int vj = jvar.number();
2813 
2814  unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2815 
2816  for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
2817  tag < _jacobian_block_neighbor_used.size();
2818  tag++)
2819  {
2821  .resize(ivar.dofIndices().size() * ivar.count(),
2822  jvar.dofIndicesNeighbor().size() * jcount);
2823 
2825  .resize(ivar.dofIndicesNeighbor().size() * ivar.count(),
2826  jvar.dofIndices().size() * jcount);
2827 
2829  .resize(ivar.dofIndicesNeighbor().size() * ivar.count(),
2830  jvar.dofIndicesNeighbor().size() * jcount);
2831 
2832  jacobianBlockNeighborUsed(tag, vi, vj, false);
2833  }
2834  }
2835 
2836  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2837  for (const auto & var : vars)
2838  for (auto & tag_Rn : _sub_Rn)
2839  tag_Rn[var->number()].resize(var->dofIndicesNeighbor().size() * var->count());
2840 }
2841 
2842 void
2844 {
2845  for (const auto & it : _cm_ff_entry)
2846  {
2847  MooseVariableFEBase & ivar = *(it.first);
2848  MooseVariableFEBase & jvar = *(it.second);
2849 
2850  unsigned int vi = ivar.number();
2851  unsigned int vj = jvar.number();
2852 
2853  unsigned int jcount = (vi == vj && _component_block_diagonal[vi]) ? 1 : jvar.count();
2854 
2855  for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
2856  tag++)
2857  {
2858  // To cover all possible cases we should have 9 combinations below for every 2-permutation
2859  // of Lower,Secondary,Primary. However, 4 cases will in general be covered by calls to
2860  // prepare() and prepareNeighbor(). These calls will cover SecondarySecondary
2861  // (ElementElement), SecondaryPrimary (ElementNeighbor), PrimarySecondary (NeighborElement),
2862  // and PrimaryPrimary (NeighborNeighbor). With these covered we only need to prepare the 5
2863  // remaining below
2864 
2865  // derivatives w.r.t. lower dimensional residuals
2867  .resize(ivar.dofIndicesLower().size() * ivar.count(),
2868  jvar.dofIndicesLower().size() * jcount);
2869 
2871  .resize(ivar.dofIndicesLower().size() * ivar.count(),
2872  jvar.dofIndices().size() * jvar.count());
2873 
2875  .resize(ivar.dofIndicesLower().size() * ivar.count(),
2876  jvar.dofIndicesNeighbor().size() * jvar.count());
2877 
2878  // derivatives w.r.t. interior secondary residuals
2880  .resize(ivar.dofIndices().size() * ivar.count(),
2881  jvar.dofIndicesLower().size() * jvar.count());
2882 
2883  // derivatives w.r.t. interior primary residuals
2885  .resize(ivar.dofIndicesNeighbor().size() * ivar.count(),
2886  jvar.dofIndicesLower().size() * jvar.count());
2887 
2888  jacobianBlockLowerUsed(tag, vi, vj, false);
2889  }
2890  }
2891 
2892  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2893  for (const auto & var : vars)
2894  for (auto & tag_Rl : _sub_Rl)
2895  tag_Rl[var->number()].resize(var->dofIndicesLower().size() * var->count());
2896 }
2897 
2898 void
2899 Assembly::prepareBlock(unsigned int ivar,
2900  unsigned int jvar,
2901  const std::vector<dof_id_type> & dof_indices)
2902 {
2903  const auto & iv = _sys.getVariable(_tid, ivar);
2904  const auto & jv = _sys.getVariable(_tid, jvar);
2905  const unsigned int ivn = iv.number();
2906  const unsigned int jvn = jv.number();
2907  const unsigned int icount = iv.count();
2908  unsigned int jcount = jv.count();
2909  if (ivn == jvn && _component_block_diagonal[ivn])
2910  jcount = 1;
2911 
2912  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2913  {
2914  jacobianBlock(ivn, jvn, LocalDataKey{}, tag)
2915  .resize(dof_indices.size() * icount, dof_indices.size() * jcount);
2916  jacobianBlockUsed(tag, ivn, jvn, false);
2917  }
2918 
2919  for (auto & tag_Re : _sub_Re)
2920  tag_Re[ivn].resize(dof_indices.size() * icount);
2921 }
2922 
2923 void
2925  unsigned int jvar,
2926  const std::vector<dof_id_type> & idof_indices,
2927  const std::vector<dof_id_type> & jdof_indices)
2928 {
2929  const auto & iv = _sys.getVariable(_tid, ivar);
2930  const auto & jv = _sys.getVariable(_tid, jvar);
2931  const unsigned int ivn = iv.number();
2932  const unsigned int jvn = jv.number();
2933  const unsigned int icount = iv.count();
2934  unsigned int jcount = jv.count();
2935  if (ivn == jvn && _component_block_diagonal[ivn])
2936  jcount = 1;
2937 
2938  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2939  tag < _jacobian_block_nonlocal_used.size();
2940  tag++)
2941  {
2942  jacobianBlockNonlocal(ivn, jvn, LocalDataKey{}, tag)
2943  .resize(idof_indices.size() * icount, jdof_indices.size() * jcount);
2944 
2945  jacobianBlockNonlocalUsed(tag, ivn, jvn, false);
2946  }
2947 }
2948 
2949 void
2951 {
2952  const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
2953  for (const auto & ivar : vars)
2954  {
2955  auto idofs = ivar->dofIndices().size();
2956 
2957  for (auto & tag_Re : _sub_Re)
2958  tag_Re[ivar->number()].resize(idofs);
2959 
2960  for (const auto & jvar : vars)
2961  {
2962  auto jdofs = jvar->dofIndices().size();
2963 
2964  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2965  {
2966  jacobianBlock(ivar->number(), jvar->number(), LocalDataKey{}, tag).resize(idofs, jdofs);
2967  jacobianBlockUsed(tag, ivar->number(), jvar->number(), false);
2968  }
2969  }
2970  }
2971 }
2972 
2973 void
2975 {
2976  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2977  const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
2978 
2979  for (const auto & ivar : scalar_vars)
2980  {
2981  auto idofs = ivar->dofIndices().size();
2982 
2983  for (const auto & jvar : vars)
2984  {
2985  auto jdofs = jvar->dofIndices().size() * jvar->count();
2986  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2987  {
2988  jacobianBlock(ivar->number(), jvar->number(), LocalDataKey{}, tag).resize(idofs, jdofs);
2989  jacobianBlockUsed(tag, ivar->number(), jvar->number(), false);
2990 
2991  jacobianBlock(jvar->number(), ivar->number(), LocalDataKey{}, tag).resize(jdofs, idofs);
2992  jacobianBlockUsed(tag, jvar->number(), ivar->number(), false);
2993  }
2994  }
2995  }
2996 }
2997 
2998 template <typename T>
2999 void
3001 {
3002  phi(v).shallowCopy(v.phi());
3003  gradPhi(v).shallowCopy(v.gradPhi());
3004  if (v.computingSecond())
3005  secondPhi(v).shallowCopy(v.secondPhi());
3006 }
3007 
3008 void
3009 Assembly::copyShapes(unsigned int var)
3010 {
3011  auto & v = _sys.getVariable(_tid, var);
3012  if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3013  {
3014  auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3015  copyShapes(v);
3016  }
3017  else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3018  {
3019  auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3020  copyShapes(v);
3021  }
3022  else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3023  {
3024  auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3025  copyShapes(v);
3026  if (v.computingCurl())
3027  curlPhi(v).shallowCopy(v.curlPhi());
3028  if (v.computingDiv())
3029  divPhi(v).shallowCopy(v.divPhi());
3030  }
3031  else
3032  mooseError("Unsupported variable field type!");
3033 }
3034 
3035 template <typename T>
3036 void
3038 {
3039  phiFace(v).shallowCopy(v.phiFace());
3040  gradPhiFace(v).shallowCopy(v.gradPhiFace());
3041  if (v.computingSecond())
3042  secondPhiFace(v).shallowCopy(v.secondPhiFace());
3043 }
3044 
3045 void
3046 Assembly::copyFaceShapes(unsigned int var)
3047 {
3048  auto & v = _sys.getVariable(_tid, var);
3049  if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3050  {
3051  auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3052  copyFaceShapes(v);
3053  }
3054  else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3055  {
3056  auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3057  copyFaceShapes(v);
3058  }
3059  else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3060  {
3061  auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3062  copyFaceShapes(v);
3063  if (v.computingCurl())
3064  _vector_curl_phi_face.shallowCopy(v.curlPhi());
3065  if (v.computingDiv())
3066  _vector_div_phi_face.shallowCopy(v.divPhi());
3067  }
3068  else
3069  mooseError("Unsupported variable field type!");
3070 }
3071 
3072 template <typename T>
3073 void
3075 {
3076  if (v.usesPhiNeighbor())
3077  {
3078  phiFaceNeighbor(v).shallowCopy(v.phiFaceNeighbor());
3079  phiNeighbor(v).shallowCopy(v.phiNeighbor());
3080  }
3081  if (v.usesGradPhiNeighbor())
3082  {
3083  gradPhiFaceNeighbor(v).shallowCopy(v.gradPhiFaceNeighbor());
3084  gradPhiNeighbor(v).shallowCopy(v.gradPhiNeighbor());
3085  }
3086  if (v.usesSecondPhiNeighbor())
3087  {
3088  secondPhiFaceNeighbor(v).shallowCopy(v.secondPhiFaceNeighbor());
3089  secondPhiNeighbor(v).shallowCopy(v.secondPhiNeighbor());
3090  }
3091 }
3092 
3093 void
3095 {
3096  auto & v = _sys.getVariable(_tid, var);
3097  if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
3098  {
3099  auto & v = _sys.getActualFieldVariable<Real>(_tid, var);
3100  copyNeighborShapes(v);
3101  }
3102  else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
3103  {
3104  auto & v = _sys.getActualFieldVariable<RealEigenVector>(_tid, var);
3105  copyNeighborShapes(v);
3106  }
3107  else if (v.fieldType() == Moose::VarFieldType::VAR_FIELD_VECTOR)
3108  {
3109  auto & v = _sys.getActualFieldVariable<RealVectorValue>(_tid, var);
3110  copyNeighborShapes(v);
3111  }
3112  else
3113  mooseError("Unsupported variable field type!");
3114 }
3115 
3118  Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
3119 {
3120  if (type == Moose::ElementElement)
3121  jacobianBlockUsed(tag, ivar, jvar, true);
3122  else
3123  jacobianBlockNeighborUsed(tag, ivar, jvar, true);
3124 
3126  {
3127  switch (type)
3128  {
3129  default:
3130  case Moose::ElementElement:
3131  return _sub_Kee[tag][ivar][0];
3133  return _sub_Ken[tag][ivar][0];
3135  return _sub_Kne[tag][ivar][0];
3137  return _sub_Knn[tag][ivar][0];
3138  }
3139  }
3140  else
3141  {
3142  switch (type)
3143  {
3144  default:
3145  case Moose::ElementElement:
3146  return _sub_Kee[tag][ivar][jvar];
3148  return _sub_Ken[tag][ivar][jvar];
3150  return _sub_Kne[tag][ivar][jvar];
3152  return _sub_Knn[tag][ivar][jvar];
3153  }
3154  }
3155 }
3156 
3159  unsigned int ivar,
3160  unsigned int jvar,
3161  LocalDataKey,
3162  TagID tag)
3163 {
3164  jacobianBlockLowerUsed(tag, ivar, jvar, true);
3166  {
3167  switch (type)
3168  {
3169  default:
3170  case Moose::LowerLower:
3171  return _sub_Kll[tag][ivar][0];
3172  case Moose::LowerSecondary:
3173  return _sub_Kle[tag][ivar][0];
3174  case Moose::LowerPrimary:
3175  return _sub_Kln[tag][ivar][0];
3176  case Moose::SecondaryLower:
3177  return _sub_Kel[tag][ivar][0];
3179  return _sub_Kee[tag][ivar][0];
3181  return _sub_Ken[tag][ivar][0];
3182  case Moose::PrimaryLower:
3183  return _sub_Knl[tag][ivar][0];
3185  return _sub_Kne[tag][ivar][0];
3186  case Moose::PrimaryPrimary:
3187  return _sub_Knn[tag][ivar][0];
3188  }
3189  }
3190  else
3191  {
3192  switch (type)
3193  {
3194  default:
3195  case Moose::LowerLower:
3196  return _sub_Kll[tag][ivar][jvar];
3197  case Moose::LowerSecondary:
3198  return _sub_Kle[tag][ivar][jvar];
3199  case Moose::LowerPrimary:
3200  return _sub_Kln[tag][ivar][jvar];
3201  case Moose::SecondaryLower:
3202  return _sub_Kel[tag][ivar][jvar];
3204  return _sub_Kee[tag][ivar][jvar];
3206  return _sub_Ken[tag][ivar][jvar];
3207  case Moose::PrimaryLower:
3208  return _sub_Knl[tag][ivar][jvar];
3210  return _sub_Kne[tag][ivar][jvar];
3211  case Moose::PrimaryPrimary:
3212  return _sub_Knn[tag][ivar][jvar];
3213  }
3214  }
3215 }
3216 
3217 void
3219  std::vector<dof_id_type> & dof_indices,
3220  const std::vector<Real> & scaling_factor,
3221  bool is_nodal)
3222 {
3223  // For an array variable, ndof is the number of dofs of the zero-th component and
3224  // ntdof is the number of dofs of all components.
3225  // For standard or vector variables, ndof will be the same as ntdof.
3226  auto ndof = dof_indices.size();
3227  auto ntdof = res_block.size();
3228  auto count = ntdof / ndof;
3229  mooseAssert(count == scaling_factor.size(), "Inconsistent of number of components");
3230  mooseAssert(count * ndof == ntdof, "Inconsistent of number of components");
3231  if (count > 1)
3232  {
3233  // expanding dof indices
3234  dof_indices.resize(ntdof);
3235  unsigned int p = 0;
3236  for (MooseIndex(count) j = 0; j < count; ++j)
3237  for (MooseIndex(ndof) i = 0; i < ndof; ++i)
3238  {
3239  dof_indices[p] = dof_indices[i] + (is_nodal ? j : j * ndof);
3240  res_block(p) *= scaling_factor[j];
3241  ++p;
3242  }
3243  }
3244  else
3245  {
3246  if (scaling_factor[0] != 1.0)
3247  res_block *= scaling_factor[0];
3248  }
3249 
3250  _dof_map.constrain_element_vector(res_block, dof_indices, false);
3251 }
3252 
3253 void
3255  DenseVector<Number> & res_block,
3256  const std::vector<dof_id_type> & dof_indices,
3257  const std::vector<Real> & scaling_factor,
3258  bool is_nodal)
3259 {
3260  if (dof_indices.size() > 0 && res_block.size())
3261  {
3262  _temp_dof_indices = dof_indices;
3263  _tmp_Re = res_block;
3264  processLocalResidual(_tmp_Re, _temp_dof_indices, scaling_factor, is_nodal);
3266  }
3267 }
3268 
3269 void
3270 Assembly::cacheResidualBlock(std::vector<Real> & cached_residual_values,
3271  std::vector<dof_id_type> & cached_residual_rows,
3272  DenseVector<Number> & res_block,
3273  const std::vector<dof_id_type> & dof_indices,
3274  const std::vector<Real> & scaling_factor,
3275  bool is_nodal)
3276 {
3277  if (dof_indices.size() > 0 && res_block.size())
3278  {
3279  _temp_dof_indices = dof_indices;
3280  _tmp_Re = res_block;
3281  processLocalResidual(_tmp_Re, _temp_dof_indices, scaling_factor, is_nodal);
3282 
3283  for (MooseIndex(_tmp_Re) i = 0; i < _tmp_Re.size(); i++)
3284  {
3285  cached_residual_values.push_back(_tmp_Re(i));
3286  cached_residual_rows.push_back(_temp_dof_indices[i]);
3287  }
3288  }
3289 
3290  res_block.zero();
3291 }
3292 
3293 void
3295  DenseVector<Number> & res_block,
3296  const std::vector<dof_id_type> & dof_indices,
3297  const std::vector<Real> & scaling_factor,
3298  bool is_nodal)
3299 {
3300  if (dof_indices.size() > 0)
3301  {
3302  std::vector<dof_id_type> di(dof_indices);
3303  _tmp_Re = res_block;
3304  processLocalResidual(_tmp_Re, di, scaling_factor, is_nodal);
3305  residual.insert(_tmp_Re, di);
3306  }
3307 }
3308 
3309 void
3311 {
3312  mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3313  "Non-residual tag in Assembly::addResidual");
3314 
3315  auto & tag_Re = _sub_Re[vector_tag._type_id];
3316  NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3317  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3318  for (const auto & var : vars)
3319  addResidualBlock(residual,
3320  tag_Re[var->number()],
3321  var->dofIndices(),
3322  var->arrayScalingFactor(),
3323  var->isNodal());
3324 }
3325 
3326 void
3327 Assembly::addResidual(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3328 {
3329  for (const auto & vector_tag : vector_tags)
3330  if (_sys.hasVector(vector_tag._id))
3331  addResidual(vector_tag);
3332 }
3333 
3334 void
3336 {
3337  mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3338  "Non-residual tag in Assembly::addResidualNeighbor");
3339 
3340  auto & tag_Rn = _sub_Rn[vector_tag._type_id];
3341  NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3342  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3343  for (const auto & var : vars)
3344  addResidualBlock(residual,
3345  tag_Rn[var->number()],
3346  var->dofIndicesNeighbor(),
3347  var->arrayScalingFactor(),
3348  var->isNodal());
3349 }
3350 
3351 void
3352 Assembly::addResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3353 {
3354  for (const auto & vector_tag : vector_tags)
3355  if (_sys.hasVector(vector_tag._id))
3356  addResidualNeighbor(vector_tag);
3357 }
3358 
3359 void
3361 {
3362  mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3363  "Non-residual tag in Assembly::addResidualLower");
3364 
3365  auto & tag_Rl = _sub_Rl[vector_tag._type_id];
3366  NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3367  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3368  for (const auto & var : vars)
3369  addResidualBlock(residual,
3370  tag_Rl[var->number()],
3371  var->dofIndicesLower(),
3372  var->arrayScalingFactor(),
3373  var->isNodal());
3374 }
3375 
3376 void
3377 Assembly::addResidualLower(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3378 {
3379  for (const auto & vector_tag : vector_tags)
3380  if (_sys.hasVector(vector_tag._id))
3381  addResidualLower(vector_tag);
3382 }
3383 
3384 // private method, so no key required
3385 void
3387 {
3388  mooseAssert(vector_tag._type == Moose::VECTOR_TAG_RESIDUAL,
3389  "Non-residual tag in Assembly::addResidualScalar");
3390 
3391  // add the scalar variables residuals
3392  auto & tag_Re = _sub_Re[vector_tag._type_id];
3393  NumericVector<Number> & residual = _sys.getVector(vector_tag._id);
3394  const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
3395  for (const auto & var : vars)
3397  residual, tag_Re[var->number()], var->dofIndices(), var->arrayScalingFactor(), false);
3398 }
3399 
3400 void
3401 Assembly::addResidualScalar(GlobalDataKey, const std::vector<VectorTag> & vector_tags)
3402 {
3403  for (const auto & vector_tag : vector_tags)
3404  if (_sys.hasVector(vector_tag._id))
3405  addResidualScalar(vector_tag);
3406 }
3407 
3408 void
3409 Assembly::cacheResidual(GlobalDataKey, const std::vector<VectorTag> & tags)
3410 {
3411  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3412  for (const auto & var : vars)
3413  for (const auto & vector_tag : tags)
3414  if (_sys.hasVector(vector_tag._id))
3415  cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3416  _cached_residual_rows[vector_tag._type_id],
3417  _sub_Re[vector_tag._type_id][var->number()],
3418  var->dofIndices(),
3419  var->arrayScalingFactor(),
3420  var->isNodal());
3421 }
3422 
3423 // private method, so no key required
3424 void
3426 {
3427  const VectorTag & tag = _subproblem.getVectorTag(tag_id);
3428 
3429  _cached_residual_values[tag._type_id].push_back(value);
3430  _cached_residual_rows[tag._type_id].push_back(dof);
3431 }
3432 
3433 // private method, so no key required
3434 void
3435 Assembly::cacheResidual(dof_id_type dof, Real value, const std::set<TagID> & tags)
3436 {
3437  for (auto & tag : tags)
3438  cacheResidual(dof, value, tag);
3439 }
3440 
3441 void
3443  const std::vector<dof_id_type> & dof_index,
3444  LocalDataKey,
3445  TagID tag)
3446 {
3447  // Add the residual value and dof_index to cached_residual_values and cached_residual_rows
3448  // respectively.
3449  // This is used by NodalConstraint.C to cache the residual calculated for primary and secondary
3450  // node.
3451  const VectorTag & vector_tag = _subproblem.getVectorTag(tag);
3452  for (MooseIndex(dof_index) i = 0; i < dof_index.size(); ++i)
3453  {
3454  _cached_residual_values[vector_tag._type_id].push_back(res(i));
3455  _cached_residual_rows[vector_tag._type_id].push_back(dof_index[i]);
3456  }
3457 }
3458 
3459 void
3460 Assembly::cacheResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & tags)
3461 {
3462  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3463  for (const auto & var : vars)
3464  for (const auto & vector_tag : tags)
3465  if (_sys.hasVector(vector_tag._id))
3466  cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3467  _cached_residual_rows[vector_tag._type_id],
3468  _sub_Rn[vector_tag._type_id][var->number()],
3469  var->dofIndicesNeighbor(),
3470  var->arrayScalingFactor(),
3471  var->isNodal());
3472 }
3473 
3474 void
3475 Assembly::cacheResidualLower(GlobalDataKey, const std::vector<VectorTag> & tags)
3476 {
3477  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3478  for (const auto & var : vars)
3479  for (const auto & vector_tag : tags)
3480  if (_sys.hasVector(vector_tag._id))
3481  cacheResidualBlock(_cached_residual_values[vector_tag._type_id],
3482  _cached_residual_rows[vector_tag._type_id],
3483  _sub_Rl[vector_tag._type_id][var->number()],
3484  var->dofIndicesLower(),
3485  var->arrayScalingFactor(),
3486  var->isNodal());
3487 }
3488 
3489 void
3490 Assembly::addCachedResiduals(GlobalDataKey, const std::vector<VectorTag> & tags)
3491 {
3492  for (const auto & vector_tag : tags)
3493  {
3494  if (!_sys.hasVector(vector_tag._id))
3495  {
3496  _cached_residual_values[vector_tag._type_id].clear();
3497  _cached_residual_rows[vector_tag._type_id].clear();
3498  continue;
3499  }
3500  addCachedResidualDirectly(_sys.getVector(vector_tag._id), GlobalDataKey{}, vector_tag);
3501  }
3502 }
3503 
3504 void
3506 {
3507  for (const auto & vector_tag : _residual_vector_tags)
3508  clearCachedResiduals(vector_tag);
3509 }
3510 
3511 // private method, so no key required
3512 void
3514 {
3515  auto & values = _cached_residual_values[vector_tag._type_id];
3516  auto & rows = _cached_residual_rows[vector_tag._type_id];
3517 
3518  mooseAssert(values.size() == rows.size(),
3519  "Number of cached residuals and number of rows must match!");
3520 
3521  // Keep track of the largest size so we can use it to reserve and avoid
3522  // as much dynamic allocation as possible
3523  if (_max_cached_residuals < values.size())
3524  _max_cached_residuals = values.size();
3525 
3526  // Clear both vectors (keeps the capacity the same)
3527  values.clear();
3528  rows.clear();
3529  // And then reserve: use 2 as a fudge factor to *really* avoid dynamic allocation!
3530  values.reserve(_max_cached_residuals * 2);
3531  rows.reserve(_max_cached_residuals * 2);
3532 }
3533 
3534 void
3536  GlobalDataKey,
3537  const VectorTag & vector_tag)
3538 {
3539  const auto & values = _cached_residual_values[vector_tag._type_id];
3540  const auto & rows = _cached_residual_rows[vector_tag._type_id];
3541 
3542  mooseAssert(values.size() == rows.size(),
3543  "Number of cached residuals and number of rows must match!");
3544 
3545  if (!values.empty())
3546  {
3547  residual.add_vector(values, rows);
3548  clearCachedResiduals(vector_tag);
3549  }
3550 }
3551 
3552 void
3554 {
3555  auto & tag_Re = _sub_Re[vector_tag._type_id];
3556  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3557  for (const auto & var : vars)
3558  setResidualBlock(residual,
3559  tag_Re[var->number()],
3560  var->dofIndices(),
3561  var->arrayScalingFactor(),
3562  var->isNodal());
3563 }
3564 
3565 void
3567  GlobalDataKey,
3568  const VectorTag & vector_tag)
3569 {
3570  auto & tag_Rn = _sub_Rn[vector_tag._type_id];
3571  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3572  for (const auto & var : vars)
3573  setResidualBlock(residual,
3574  tag_Rn[var->number()],
3575  var->dofIndicesNeighbor(),
3576  var->arrayScalingFactor(),
3577  var->isNodal());
3578 }
3579 
3580 // private method, so no key required
3581 void
3583  DenseMatrix<Number> & jac_block,
3584  const MooseVariableBase & ivar,
3585  const MooseVariableBase & jvar,
3586  const std::vector<dof_id_type> & idof_indices,
3587  const std::vector<dof_id_type> & jdof_indices)
3588 {
3589  if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3590  return;
3591  if (jac_block.n() == 0 || jac_block.m() == 0)
3592  return;
3593 
3594  auto & scaling_factor = ivar.arrayScalingFactor();
3595 
3596  for (unsigned int i = 0; i < ivar.count(); ++i)
3597  {
3598  unsigned int iv = ivar.number();
3599  for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3600  {
3601  unsigned int jv = jvar.number();
3602  if (jt < jv || jt >= jv + jvar.count())
3603  continue;
3604  unsigned int j = jt - jv;
3605 
3606  auto di = ivar.componentDofIndices(idof_indices, i);
3607  auto dj = jvar.componentDofIndices(jdof_indices, j);
3608  auto indof = di.size();
3609  auto jndof = dj.size();
3610 
3611  unsigned int jj = j;
3612  if (iv == jv && _component_block_diagonal[iv])
3613  // here i must be equal to j
3614  jj = 0;
3615 
3616  auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3617  if (scaling_factor[i] != 1.0)
3618  sub *= scaling_factor[i];
3619 
3620  // If we're computing the jacobian for automatically scaling variables we do not want
3621  // to constrain the element matrix because it introduces 1s on the diagonal for the
3622  // constrained dofs
3624  _dof_map.constrain_element_matrix(sub, di, dj, false);
3625 
3626  jacobian.add_matrix(sub, di, dj);
3627  }
3628  }
3629 }
3630 
3631 // private method, so no key required
3632 void
3634  const MooseVariableBase & ivar,
3635  const MooseVariableBase & jvar,
3636  const std::vector<dof_id_type> & idof_indices,
3637  const std::vector<dof_id_type> & jdof_indices,
3638  TagID tag)
3639 {
3640  if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3641  return;
3642  if (jac_block.n() == 0 || jac_block.m() == 0)
3643  return;
3644  if (!_sys.hasMatrix(tag))
3645  return;
3646 
3647  auto & scaling_factor = ivar.arrayScalingFactor();
3648 
3649  for (unsigned int i = 0; i < ivar.count(); ++i)
3650  {
3651  unsigned int iv = ivar.number();
3652  for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3653  {
3654  unsigned int jv = jvar.number();
3655  if (jt < jv || jt >= jv + jvar.count())
3656  continue;
3657  unsigned int j = jt - jv;
3658 
3659  auto di = ivar.componentDofIndices(idof_indices, i);
3660  auto dj = jvar.componentDofIndices(jdof_indices, j);
3661  auto indof = di.size();
3662  auto jndof = dj.size();
3663 
3664  unsigned int jj = j;
3665  if (iv == jv && _component_block_diagonal[iv])
3666  // here i must be equal to j
3667  jj = 0;
3668 
3669  auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3670  if (scaling_factor[i] != 1.0)
3671  sub *= scaling_factor[i];
3672 
3673  // If we're computing the jacobian for automatically scaling variables we do not want
3674  // to constrain the element matrix because it introduces 1s on the diagonal for the
3675  // constrained dofs
3677  _dof_map.constrain_element_matrix(sub, di, dj, false);
3678 
3679  for (MooseIndex(di) i = 0; i < di.size(); i++)
3680  for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3681  {
3682  _cached_jacobian_values[tag].push_back(sub(i, j));
3683  _cached_jacobian_rows[tag].push_back(di[i]);
3684  _cached_jacobian_cols[tag].push_back(dj[j]);
3685  }
3686  }
3687  }
3688 
3689  jac_block.zero();
3690 }
3691 
3692 // private method, so no key required
3693 void
3695  const MooseVariableBase & ivar,
3696  const MooseVariableBase & jvar,
3697  const std::vector<dof_id_type> & idof_indices,
3698  const std::vector<dof_id_type> & jdof_indices,
3699  TagID tag)
3700 {
3701  if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3702  return;
3703  if (jac_block.n() == 0 || jac_block.m() == 0)
3704  return;
3705  if (!_sys.hasMatrix(tag))
3706  return;
3707 
3708  auto & scaling_factor = ivar.arrayScalingFactor();
3709 
3710  for (unsigned int i = 0; i < ivar.count(); ++i)
3711  {
3712  unsigned int iv = ivar.number();
3713  for (const auto & jt : ConstCouplingRow(iv + i, *_cm))
3714  {
3715  unsigned int jv = jvar.number();
3716  if (jt < jv || jt >= jv + jvar.count())
3717  continue;
3718  unsigned int j = jt - jv;
3719 
3720  auto di = ivar.componentDofIndices(idof_indices, i);
3721  auto dj = jvar.componentDofIndices(jdof_indices, j);
3722  auto indof = di.size();
3723  auto jndof = dj.size();
3724 
3725  unsigned int jj = j;
3726  if (iv == jv && _component_block_diagonal[iv])
3727  // here i must be equal to j
3728  jj = 0;
3729 
3730  auto sub = jac_block.sub_matrix(i * indof, indof, jj * jndof, jndof);
3731  if (scaling_factor[i] != 1.0)
3732  sub *= scaling_factor[i];
3733 
3734  _dof_map.constrain_element_matrix(sub, di, dj, false);
3735 
3736  for (MooseIndex(di) i = 0; i < di.size(); i++)
3737  for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3738  if (sub(i, j) != 0.0) // no storage allocated for unimplemented jacobian terms,
3739  // maintaining maximum sparsity possible
3740  {
3741  _cached_jacobian_values[tag].push_back(sub(i, j));
3742  _cached_jacobian_rows[tag].push_back(di[i]);
3743  _cached_jacobian_cols[tag].push_back(dj[j]);
3744  }
3745  }
3746  }
3747 
3748  jac_block.zero();
3749 }
3750 
3751 void
3753  const std::vector<dof_id_type> & idof_indices,
3754  const std::vector<dof_id_type> & jdof_indices,
3755  Real scaling_factor,
3756  LocalDataKey,
3757  TagID tag)
3758 {
3759  // Only cache data when the matrix exists
3760  if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m() &&
3761  _sys.hasMatrix(tag))
3762  {
3763  std::vector<dof_id_type> di(idof_indices);
3764  std::vector<dof_id_type> dj(jdof_indices);
3765 
3766  // If we're computing the jacobian for automatically scaling variables we do not want to
3767  // constrain the element matrix because it introduces 1s on the diagonal for the constrained
3768  // dofs
3770  _dof_map.constrain_element_matrix(jac_block, di, dj, false);
3771 
3772  if (scaling_factor != 1.0)
3773  jac_block *= scaling_factor;
3774 
3775  for (MooseIndex(di) i = 0; i < di.size(); i++)
3776  for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3777  {
3778  _cached_jacobian_values[tag].push_back(jac_block(i, j));
3779  _cached_jacobian_rows[tag].push_back(di[i]);
3780  _cached_jacobian_cols[tag].push_back(dj[j]);
3781  }
3782  }
3783  jac_block.zero();
3784 }
3785 
3786 Real
3787 Assembly::elementVolume(const Elem * elem) const
3788 {
3789  FEType fe_type(elem->default_order(), LAGRANGE);
3790  std::unique_ptr<FEBase> fe(FEBase::build(elem->dim(), fe_type));
3791 
3792  // references to the quadrature points and weights
3793  const std::vector<Real> & JxW = fe->get_JxW();
3794  const std::vector<Point> & q_points = fe->get_xyz();
3795 
3796  // The default quadrature rule should integrate the mass matrix,
3797  // thus it should be plenty to compute the volume
3798  QGauss qrule(elem->dim(), fe_type.default_quadrature_order());
3799  fe->attach_quadrature_rule(&qrule);
3800  fe->reinit(elem);
3801 
3802  // perform a sanity check to ensure that size of quad rule and size of q_points is
3803  // identical
3804  mooseAssert(qrule.n_points() == q_points.size(),
3805  "The number of points in the quadrature rule doesn't match the number of passed-in "
3806  "points in Assembly::setCoordinateTransformation");
3807 
3808  // compute the coordinate transformation
3809  Real vol = 0;
3810  for (unsigned int qp = 0; qp < qrule.n_points(); ++qp)
3811  {
3812  Real coord;
3813  coordTransformFactor(_subproblem, elem->subdomain_id(), q_points[qp], coord);
3814  vol += JxW[qp] * coord;
3815  }
3816  return vol;
3817 }
3818 
3819 void
3821 {
3822 #ifndef NDEBUG
3824  {
3825  mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
3826  "Error: Cached data sizes MUST be the same!");
3827  for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3828  mooseAssert(_cached_jacobian_rows[i].size() == _cached_jacobian_cols[i].size(),
3829  "Error: Cached data sizes MUST be the same for a given tag!");
3830  }
3831 #endif
3832 
3833  for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3834  if (_sys.hasMatrix(i))
3835  for (MooseIndex(_cached_jacobian_rows[i]) j = 0; j < _cached_jacobian_rows[i].size(); j++)
3837  _cached_jacobian_cols[i][j],
3838  _cached_jacobian_values[i][j]);
3839 
3840  for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3841  {
3842  if (!_sys.hasMatrix(i))
3843  continue;
3844 
3847 
3848  // Try to be more efficient from now on
3849  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
3850  _cached_jacobian_values[i].clear();
3852 
3853  _cached_jacobian_rows[i].clear();
3855 
3856  _cached_jacobian_cols[i].clear();
3858  }
3859 }
3860 
3861 inline void
3863 {
3864  auto i = ivar.number();
3865  auto j = jvar.number();
3866  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
3867  if (jacobianBlockUsed(tag, i, j) && _sys.hasMatrix(tag))
3869  jacobianBlock(i, j, LocalDataKey{}, tag),
3870  ivar,
3871  jvar,
3872  ivar.dofIndices(),
3873  jvar.dofIndices());
3874 }
3875 
3876 void
3878 {
3879  for (const auto & it : _cm_ff_entry)
3880  addJacobianCoupledVarPair(*it.first, *it.second);
3881 
3882  for (const auto & it : _cm_sf_entry)
3883  addJacobianCoupledVarPair(*it.first, *it.second);
3884 
3885  for (const auto & it : _cm_fs_entry)
3886  addJacobianCoupledVarPair(*it.first, *it.second);
3887 }
3888 
3889 void
3891 {
3892  for (const auto & it : _cm_nonlocal_entry)
3893  {
3894  auto ivar = it.first;
3895  auto jvar = it.second;
3896  auto i = ivar->number();
3897  auto j = jvar->number();
3898  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
3899  tag < _jacobian_block_nonlocal_used.size();
3900  tag++)
3901  if (jacobianBlockNonlocalUsed(tag, i, j) && _sys.hasMatrix(tag))
3903  jacobianBlockNonlocal(i, j, LocalDataKey{}, tag),
3904  *ivar,
3905  *jvar,
3906  ivar->dofIndices(),
3907  jvar->allDofIndices());
3908  }
3909 }
3910 
3911 void
3913 {
3914  for (const auto & it : _cm_ff_entry)
3915  {
3916  auto ivar = it.first;
3917  auto jvar = it.second;
3918  auto i = ivar->number();
3919  auto j = jvar->number();
3920  for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3921  tag < _jacobian_block_neighbor_used.size();
3922  tag++)
3923  if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
3924  {
3927  *ivar,
3928  *jvar,
3929  ivar->dofIndices(),
3930  jvar->dofIndicesNeighbor());
3931 
3934  *ivar,
3935  *jvar,
3936  ivar->dofIndicesNeighbor(),
3937  jvar->dofIndices());
3938 
3941  *ivar,
3942  *jvar,
3943  ivar->dofIndicesNeighbor(),
3944  jvar->dofIndicesNeighbor());
3945  }
3946  }
3947 }
3948 
3949 void
3951 {
3952  for (const auto & it : _cm_ff_entry)
3953  {
3954  auto ivar = it.first;
3955  auto jvar = it.second;
3956  auto i = ivar->number();
3957  auto j = jvar->number();
3958  for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
3959  tag++)
3960  if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
3961  {
3964  *ivar,
3965  *jvar,
3966  ivar->dofIndicesLower(),
3967  jvar->dofIndicesLower());
3968 
3971  *ivar,
3972  *jvar,
3973  ivar->dofIndicesLower(),
3974  jvar->dofIndicesNeighbor());
3975 
3978  *ivar,
3979  *jvar,
3980  ivar->dofIndicesLower(),
3981  jvar->dofIndices());
3982 
3985  *ivar,
3986  *jvar,
3987  ivar->dofIndicesNeighbor(),
3988  jvar->dofIndicesLower());
3989 
3992  *ivar,
3993  *jvar,
3994  ivar->dofIndices(),
3995  jvar->dofIndicesLower());
3996  }
3997 
3998  for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3999  tag < _jacobian_block_neighbor_used.size();
4000  tag++)
4001  if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
4002  {
4005  *ivar,
4006  *jvar,
4007  ivar->dofIndices(),
4008  jvar->dofIndicesNeighbor());
4009 
4012  *ivar,
4013  *jvar,
4014  ivar->dofIndicesNeighbor(),
4015  jvar->dofIndices());
4016 
4019  *ivar,
4020  *jvar,
4021  ivar->dofIndicesNeighbor(),
4022  jvar->dofIndicesNeighbor());
4023  }
4024  }
4025 }
4026 
4027 void
4029 {
4030  for (const auto & it : _cm_ff_entry)
4031  {
4032  auto ivar = it.first;
4033  auto jvar = it.second;
4034  auto i = ivar->number();
4035  auto j = jvar->number();
4036  for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
4037  tag++)
4038  if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
4039  {
4042  *ivar,
4043  *jvar,
4044  ivar->dofIndicesLower(),
4045  jvar->dofIndicesLower());
4046 
4049  *ivar,
4050  *jvar,
4051  ivar->dofIndicesLower(),
4052  jvar->dofIndices());
4053 
4056  *ivar,
4057  *jvar,
4058  ivar->dofIndices(),
4059  jvar->dofIndicesLower());
4060  }
4061  }
4062 }
4063 
4064 void
4066 {
4067  for (const auto & it : _cm_ff_entry)
4068  cacheJacobianCoupledVarPair(*it.first, *it.second);
4069 
4070  for (const auto & it : _cm_fs_entry)
4071  cacheJacobianCoupledVarPair(*it.first, *it.second);
4072 
4073  for (const auto & it : _cm_sf_entry)
4074  cacheJacobianCoupledVarPair(*it.first, *it.second);
4075 }
4076 
4077 // private method, so no key required
4078 void
4080  const MooseVariableBase & jvar)
4081 {
4082  auto i = ivar.number();
4083  auto j = jvar.number();
4084  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
4085  if (jacobianBlockUsed(tag, i, j) && _sys.hasMatrix(tag))
4087  ivar,
4088  jvar,
4089  ivar.dofIndices(),
4090  jvar.dofIndices(),
4091  tag);
4092 }
4093 
4094 void
4096 {
4097  for (const auto & it : _cm_nonlocal_entry)
4098  {
4099  auto ivar = it.first;
4100  auto jvar = it.second;
4101  auto i = ivar->number();
4102  auto j = jvar->number();
4103  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
4104  tag < _jacobian_block_nonlocal_used.size();
4105  tag++)
4106  if (jacobianBlockNonlocalUsed(tag, i, j) && _sys.hasMatrix(tag))
4108  *ivar,
4109  *jvar,
4110  ivar->dofIndices(),
4111  jvar->allDofIndices(),
4112  tag);
4113  }
4114 }
4115 
4116 void
4118 {
4119  for (const auto & it : _cm_ff_entry)
4120  {
4121  auto ivar = it.first;
4122  auto jvar = it.second;
4123  auto i = ivar->number();
4124  auto j = jvar->number();
4125 
4126  for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
4127  tag < _jacobian_block_neighbor_used.size();
4128  tag++)
4129  if (jacobianBlockNeighborUsed(tag, i, j) && _sys.hasMatrix(tag))
4130  {
4132  *ivar,
4133  *jvar,
4134  ivar->dofIndices(),
4135  jvar->dofIndicesNeighbor(),
4136  tag);
4138  *ivar,
4139  *jvar,
4140  ivar->dofIndicesNeighbor(),
4141  jvar->dofIndices(),
4142  tag);
4145  *ivar,
4146  *jvar,
4147  ivar->dofIndicesNeighbor(),
4148  jvar->dofIndicesNeighbor(),
4149  tag);
4150  }
4151  }
4152 }
4153 
4154 void
4156 {
4157  for (const auto & it : _cm_ff_entry)
4158  {
4159  auto ivar = it.first;
4160  auto jvar = it.second;
4161  auto i = ivar->number();
4162  auto j = jvar->number();
4163  for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
4164  tag++)
4165  if (jacobianBlockLowerUsed(tag, i, j) && _sys.hasMatrix(tag))
4166  {
4168  *ivar,
4169  *jvar,
4170  ivar->dofIndicesLower(),
4171  jvar->dofIndicesLower(),
4172  tag);
4173 
4175  *ivar,
4176  *jvar,
4177  ivar->dofIndicesLower(),
4178  jvar->dofIndices(),
4179  tag);
4180 
4182  *ivar,
4183  *jvar,
4184  ivar->dofIndicesLower(),
4185  jvar->dofIndicesNeighbor(),
4186  tag);
4187 
4189  *ivar,
4190  *jvar,
4191  ivar->dofIndices(),
4192  jvar->dofIndicesLower(),
4193  tag);
4194 
4197  *ivar,
4198  *jvar,
4199  ivar->dofIndices(),
4200  jvar->dofIndices(),
4201  tag);
4202 
4204  *ivar,
4205  *jvar,
4206  ivar->dofIndices(),
4207  jvar->dofIndicesNeighbor(),
4208  tag);
4209 
4211  *ivar,
4212  *jvar,
4213  ivar->dofIndicesNeighbor(),
4214  jvar->dofIndicesLower(),
4215  tag);
4216 
4218  *ivar,
4219  *jvar,
4220  ivar->dofIndicesNeighbor(),
4221  jvar->dofIndices(),
4222  tag);
4223 
4225  *ivar,
4226  *jvar,
4227  ivar->dofIndicesNeighbor(),
4228  jvar->dofIndicesNeighbor(),
4229  tag);
4230  }
4231  }
4232 }
4233 
4234 void
4236  unsigned int ivar,
4237  unsigned int jvar,
4238  const DofMap & dof_map,
4239  std::vector<dof_id_type> & dof_indices,
4240  GlobalDataKey,
4241  const std::set<TagID> & tags)
4242 {
4243  for (auto tag : tags)
4244  addJacobianBlock(jacobian, ivar, jvar, dof_map, dof_indices, GlobalDataKey{}, tag);
4245 }
4246 
4247 void
4249  unsigned int ivar,
4250  unsigned int jvar,
4251  const DofMap & dof_map,
4252  std::vector<dof_id_type> & dof_indices,
4253  GlobalDataKey,
4254  TagID tag)
4255 {
4256  if (dof_indices.size() == 0)
4257  return;
4258  if (!(*_cm)(ivar, jvar))
4259  return;
4260 
4261  auto & iv = _sys.getVariable(_tid, ivar);
4262  auto & jv = _sys.getVariable(_tid, jvar);
4263  auto & scaling_factor = iv.arrayScalingFactor();
4264 
4265  const unsigned int ivn = iv.number();
4266  const unsigned int jvn = jv.number();
4267  auto & ke = jacobianBlock(ivn, jvn, LocalDataKey{}, tag);
4268 
4269  // It is guaranteed by design iv.number <= ivar since iv is obtained
4270  // through SystemBase::getVariable with ivar.
4271  // Most of times ivar will just be equal to iv.number except for array variables,
4272  // where ivar could be a number for a component of an array variable but calling
4273  // getVariable will return the array variable that has the number of the 0th component.
4274  // It is the same for jvar.
4275  const unsigned int i = ivar - ivn;
4276  const unsigned int j = jvar - jvn;
4277 
4278  // DoF indices are independently given
4279  auto di = dof_indices;
4280  auto dj = dof_indices;
4281 
4282  auto indof = di.size();
4283  auto jndof = dj.size();
4284 
4285  unsigned int jj = j;
4286  if (ivar == jvar && _component_block_diagonal[ivn])
4287  jj = 0;
4288 
4289  auto sub = ke.sub_matrix(i * indof, indof, jj * jndof, jndof);
4290  // If we're computing the jacobian for automatically scaling variables we do not want to
4291  // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4292  // dofs
4294  dof_map.constrain_element_matrix(sub, di, dj, false);
4295 
4296  if (scaling_factor[i] != 1.0)
4297  sub *= scaling_factor[i];
4298 
4299  jacobian.add_matrix(sub, di, dj);
4300 }
4301 
4302 void
4304  const unsigned int ivar,
4305  const unsigned int jvar,
4306  const DofMap & dof_map,
4307  const std::vector<dof_id_type> & idof_indices,
4308  const std::vector<dof_id_type> & jdof_indices,
4309  GlobalDataKey,
4310  const TagID tag)
4311 {
4312  if (idof_indices.size() == 0 || jdof_indices.size() == 0)
4313  return;
4314  if (jacobian.n() == 0 || jacobian.m() == 0)
4315  return;
4316  if (!(*_cm)(ivar, jvar))
4317  return;
4318 
4319  auto & iv = _sys.getVariable(_tid, ivar);
4320  auto & jv = _sys.getVariable(_tid, jvar);
4321  auto & scaling_factor = iv.arrayScalingFactor();
4322 
4323  const unsigned int ivn = iv.number();
4324  const unsigned int jvn = jv.number();
4325  auto & keg = jacobianBlockNonlocal(ivn, jvn, LocalDataKey{}, tag);
4326 
4327  // It is guaranteed by design iv.number <= ivar since iv is obtained
4328  // through SystemBase::getVariable with ivar.
4329  // Most of times ivar will just be equal to iv.number except for array variables,
4330  // where ivar could be a number for a component of an array variable but calling
4331  // getVariable will return the array variable that has the number of the 0th component.
4332  // It is the same for jvar.
4333  const unsigned int i = ivar - ivn;
4334  const unsigned int j = jvar - jvn;
4335 
4336  // DoF indices are independently given
4337  auto di = idof_indices;
4338  auto dj = jdof_indices;
4339 
4340  auto indof = di.size();
4341  auto jndof = dj.size();
4342 
4343  unsigned int jj = j;
4344  if (ivar == jvar && _component_block_diagonal[ivn])
4345  jj = 0;
4346 
4347  auto sub = keg.sub_matrix(i * indof, indof, jj * jndof, jndof);
4348  // If we're computing the jacobian for automatically scaling variables we do not want to
4349  // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4350  // dofs
4352  dof_map.constrain_element_matrix(sub, di, dj, false);
4353 
4354  if (scaling_factor[i] != 1.0)
4355  sub *= scaling_factor[i];
4356 
4357  jacobian.add_matrix(sub, di, dj);
4358 }
4359 
4360 void
4362  const unsigned int ivar,
4363  const unsigned int jvar,
4364  const DofMap & dof_map,
4365  const std::vector<dof_id_type> & idof_indices,
4366  const std::vector<dof_id_type> & jdof_indices,
4367  GlobalDataKey,
4368  const std::set<TagID> & tags)
4369 {
4370  for (auto tag : tags)
4372  jacobian, ivar, jvar, dof_map, idof_indices, jdof_indices, GlobalDataKey{}, tag);
4373 }
4374 
4375 void
4377  const unsigned int ivar,
4378  const unsigned int jvar,
4379  const DofMap & dof_map,
4380  std::vector<dof_id_type> & dof_indices,
4381  std::vector<dof_id_type> & neighbor_dof_indices,
4382  GlobalDataKey,
4383  const TagID tag)
4384 {
4385  if (dof_indices.size() == 0 && neighbor_dof_indices.size() == 0)
4386  return;
4387  if (!(*_cm)(ivar, jvar))
4388  return;
4389 
4390  auto & iv = _sys.getVariable(_tid, ivar);
4391  auto & jv = _sys.getVariable(_tid, jvar);
4392  auto & scaling_factor = iv.arrayScalingFactor();
4393 
4394  const unsigned int ivn = iv.number();
4395  const unsigned int jvn = jv.number();
4396  auto & ken = jacobianBlockNeighbor(Moose::ElementNeighbor, ivn, jvn, LocalDataKey{}, tag);
4397  auto & kne = jacobianBlockNeighbor(Moose::NeighborElement, ivn, jvn, LocalDataKey{}, tag);
4398  auto & knn = jacobianBlockNeighbor(Moose::NeighborNeighbor, ivn, jvn, LocalDataKey{}, tag);
4399 
4400  // It is guaranteed by design iv.number <= ivar since iv is obtained
4401  // through SystemBase::getVariable with ivar.
4402  // Most of times ivar will just be equal to iv.number except for array variables,
4403  // where ivar could be a number for a component of an array variable but calling
4404  // getVariable will return the array variable that has the number of the 0th component.
4405  // It is the same for jvar.
4406  const unsigned int i = ivar - ivn;
4407  const unsigned int j = jvar - jvn;
4408  // DoF indices are independently given
4409  auto dc = dof_indices;
4410  auto dn = neighbor_dof_indices;
4411  auto cndof = dc.size();
4412  auto nndof = dn.size();
4413 
4414  unsigned int jj = j;
4415  if (ivar == jvar && _component_block_diagonal[ivn])
4416  jj = 0;
4417 
4418  auto suben = ken.sub_matrix(i * cndof, cndof, jj * nndof, nndof);
4419  auto subne = kne.sub_matrix(i * nndof, nndof, jj * cndof, cndof);
4420  auto subnn = knn.sub_matrix(i * nndof, nndof, jj * nndof, nndof);
4421 
4422  // If we're computing the jacobian for automatically scaling variables we do not want to
4423  // constrain the element matrix because it introduces 1s on the diagonal for the constrained
4424  // dofs
4426  {
4427  dof_map.constrain_element_matrix(suben, dc, dn, false);
4428  dof_map.constrain_element_matrix(subne, dn, dc, false);
4429  dof_map.constrain_element_matrix(subnn, dn, dn, false);
4430  }
4431 
4432  if (scaling_factor[i] != 1.0)
4433  {
4434  suben *= scaling_factor[i];
4435  subne *= scaling_factor[i];
4436  subnn *= scaling_factor[i];
4437  }
4438 
4439  jacobian.add_matrix(suben, dc, dn);
4440  jacobian.add_matrix(subne, dn, dc);
4441  jacobian.add_matrix(subnn, dn, dn);
4442 }
4443 
4444 void
4446  const unsigned int ivar,
4447  const unsigned int jvar,
4448  const DofMap & dof_map,
4449  std::vector<dof_id_type> & dof_indices,
4450  std::vector<dof_id_type> & neighbor_dof_indices,
4451  GlobalDataKey,
4452  const std::set<TagID> & tags)
4453 {
4454  for (const auto tag : tags)
4456  jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices, GlobalDataKey{}, tag);
4457 }
4458 
4459 void
4461 {
4462  for (const auto & it : _cm_ss_entry)
4463  addJacobianCoupledVarPair(*it.first, *it.second);
4464 }
4465 
4466 void
4468 {
4469  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
4471  for (const auto & var_j : vars)
4472  addJacobianCoupledVarPair(var_i, *var_j);
4473 }
4474 
4475 void
4478 {
4479  _cached_jacobian_rows[tag].push_back(i);
4480  _cached_jacobian_cols[tag].push_back(j);
4481  _cached_jacobian_values[tag].push_back(value);
4482 }
4483 
4484 void
4487  Real value,
4488  LocalDataKey,
4489  const std::set<TagID> & tags)
4490 {
4491  for (auto tag : tags)
4492  if (_sys.hasMatrix(tag))
4493  cacheJacobian(i, j, value, LocalDataKey{}, tag);
4494 }
4495 
4496 void
4498 {
4499  for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4500  if (_sys.hasMatrix(tag))
4501  {
4502  // First zero the rows (including the diagonals) to prepare for
4503  // setting the cached values.
4505 
4506  // TODO: Use SparseMatrix::set_values() for efficiency
4507  for (MooseIndex(_cached_jacobian_values) i = 0; i < _cached_jacobian_values[tag].size(); ++i)
4508  _sys.getMatrix(tag).set(_cached_jacobian_rows[tag][i],
4509  _cached_jacobian_cols[tag][i],
4510  _cached_jacobian_values[tag][i]);
4511  }
4512 
4514 }
4515 
4516 void
4518 {
4519  for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4520  if (_sys.hasMatrix(tag))
4522 
4524 }
4525 
4526 void
4528 {
4529  for (MooseIndex(_cached_jacobian_rows) tag = 0; tag < _cached_jacobian_rows.size(); tag++)
4530  {
4531  _cached_jacobian_rows[tag].clear();
4532  _cached_jacobian_cols[tag].clear();
4533  _cached_jacobian_values[tag].clear();
4534  }
4535 }
4536 
4537 void
4539 {
4540  mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
4541 
4543  return;
4544 
4545  MooseArray<Real> xfem_weight_multipliers;
4546  if (_xfem->getXFEMWeights(xfem_weight_multipliers, elem, _current_qrule, _current_q_points))
4547  {
4548  mooseAssert(xfem_weight_multipliers.size() == _current_JxW.size(),
4549  "Size of weight multipliers in xfem doesn't match number of quadrature points");
4550  for (unsigned i = 0; i < xfem_weight_multipliers.size(); i++)
4551  _current_JxW[i] = _current_JxW[i] * xfem_weight_multipliers[i];
4552 
4553  xfem_weight_multipliers.release();
4554  }
4555 }
4556 
4557 void
4558 Assembly::modifyFaceWeightsDueToXFEM(const Elem * elem, unsigned int side)
4559 {
4560  mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
4561 
4563  return;
4564 
4565  MooseArray<Real> xfem_face_weight_multipliers;
4566  if (_xfem->getXFEMFaceWeights(
4567  xfem_face_weight_multipliers, elem, _current_qrule_face, _current_q_points_face, side))
4568  {
4569  mooseAssert(xfem_face_weight_multipliers.size() == _current_JxW_face.size(),
4570  "Size of weight multipliers in xfem doesn't match number of quadrature points");
4571  for (unsigned i = 0; i < xfem_face_weight_multipliers.size(); i++)
4572  _current_JxW_face[i] = _current_JxW_face[i] * xfem_face_weight_multipliers[i];
4573 
4574  xfem_face_weight_multipliers.release();
4575  }
4576 }
4577 
4578 void
4580 {
4581  _scaling_vector = &_sys.getVector("scaling_factors");
4582 }
4583 
4584 void
4585 Assembly::modifyArbitraryWeights(const std::vector<Real> & weights)
4586 {
4587  mooseAssert(_current_qrule == _current_qrule_arbitrary, "Rule should be arbitrary");
4588  mooseAssert(weights.size() == _current_physical_points.size(), "Size mismatch");
4589 
4590  for (MooseIndex(weights.size()) i = 0; i < weights.size(); ++i)
4591  _current_JxW[i] = weights[i];
4592 }
4593 
4594 template <>
4596 Assembly::fePhi<VectorValue<Real>>(FEType type) const
4597 {
4598  buildVectorFE(type);
4599  return _vector_fe_shape_data[type]->_phi;
4600 }
4601 
4602 template <>
4604 Assembly::feGradPhi<VectorValue<Real>>(FEType type) const
4605 {
4606  buildVectorFE(type);
4607  return _vector_fe_shape_data[type]->_grad_phi;
4608 }
4609 
4610 template <>
4612 Assembly::feSecondPhi<VectorValue<Real>>(FEType type) const
4613 {
4614  _need_second_derivative.insert(type);
4615  buildVectorFE(type);
4616  return _vector_fe_shape_data[type]->_second_phi;
4617 }
4618 
4619 template <>
4621 Assembly::fePhiLower<VectorValue<Real>>(FEType type) const
4622 {
4623  buildVectorLowerDFE(type);
4624  return _vector_fe_shape_data_lower[type]->_phi;
4625 }
4626 
4627 template <>
4629 Assembly::feDualPhiLower<VectorValue<Real>>(FEType type) const
4630 {
4631  buildVectorDualLowerDFE(type);
4632  return _vector_fe_shape_data_dual_lower[type]->_phi;
4633 }
4634 
4635 template <>
4637 Assembly::feGradPhiLower<VectorValue<Real>>(FEType type) const
4638 {
4639  buildVectorLowerDFE(type);
4640  return _vector_fe_shape_data_lower[type]->_grad_phi;
4641 }
4642 
4643 template <>
4645 Assembly::feGradDualPhiLower<VectorValue<Real>>(FEType type) const
4646 {
4647  buildVectorDualLowerDFE(type);
4648  return _vector_fe_shape_data_dual_lower[type]->_grad_phi;
4649 }
4650 
4651 template <>
4653 Assembly::fePhiFace<VectorValue<Real>>(FEType type) const
4654 {
4655  buildVectorFaceFE(type);
4656  return _vector_fe_shape_data_face[type]->_phi;
4657 }
4658 
4659 template <>
4661 Assembly::feGradPhiFace<VectorValue<Real>>(FEType type) const
4662 {
4663  buildVectorFaceFE(type);
4664  return _vector_fe_shape_data_face[type]->_grad_phi;
4665 }
4666 
4667 template <>
4669 Assembly::feSecondPhiFace<VectorValue<Real>>(FEType type) const
4670 {
4671  _need_second_derivative.insert(type);
4672  buildVectorFaceFE(type);
4673 
4674  // If we're building for a face we probably need to build for a
4675  // neighbor while _need_second_derivative is set;
4676  // onInterface/reinitNeighbor/etc don't distinguish
4677  buildVectorFaceNeighborFE(type);
4678 
4679  return _vector_fe_shape_data_face[type]->_second_phi;
4680 }
4681 
4682 template <>
4684 Assembly::fePhiNeighbor<VectorValue<Real>>(FEType type) const
4685 {
4686  buildVectorNeighborFE(type);
4687  return _vector_fe_shape_data_neighbor[type]->_phi;
4688 }
4689 
4690 template <>
4692 Assembly::feGradPhiNeighbor<VectorValue<Real>>(FEType type) const
4693 {
4694  buildVectorNeighborFE(type);
4695  return _vector_fe_shape_data_neighbor[type]->_grad_phi;
4696 }
4697 
4698 template <>
4700 Assembly::feSecondPhiNeighbor<VectorValue<Real>>(FEType type) const
4701 {
4702  _need_second_derivative_neighbor.insert(type);
4703  buildVectorNeighborFE(type);
4704  return _vector_fe_shape_data_neighbor[type]->_second_phi;
4705 }
4706 
4707 template <>
4709 Assembly::fePhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4710 {
4711  buildVectorFaceNeighborFE(type);
4712  return _vector_fe_shape_data_face_neighbor[type]->_phi;
4713 }
4714 
4715 template <>
4717 Assembly::feGradPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4718 {
4719  buildVectorFaceNeighborFE(type);
4720  return _vector_fe_shape_data_face_neighbor[type]->_grad_phi;
4721 }
4722 
4723 template <>
4725 Assembly::feSecondPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4726 {
4727  _need_second_derivative_neighbor.insert(type);
4728  buildVectorFaceNeighborFE(type);
4729  return _vector_fe_shape_data_face_neighbor[type]->_second_phi;
4730 }
4731 
4732 template <>
4734 Assembly::feCurlPhi<VectorValue<Real>>(FEType type) const
4735 {
4736  _need_curl.insert(type);
4737  buildVectorFE(type);
4738  return _vector_fe_shape_data[type]->_curl_phi;
4739 }
4740 
4741 template <>
4743 Assembly::feCurlPhiFace<VectorValue<Real>>(FEType type) const
4744 {
4745  _need_curl.insert(type);
4746  buildVectorFaceFE(type);
4747 
4748  // If we're building for a face we probably need to build for a
4749  // neighbor while _need_curl is set;
4750  // onInterface/reinitNeighbor/etc don't distinguish
4751  buildVectorFaceNeighborFE(type);
4752 
4753  return _vector_fe_shape_data_face[type]->_curl_phi;
4754 }
4755 
4756 template <>
4758 Assembly::feCurlPhiNeighbor<VectorValue<Real>>(FEType type) const
4759 {
4760  _need_curl.insert(type);
4761  buildVectorNeighborFE(type);
4762  return _vector_fe_shape_data_neighbor[type]->_curl_phi;
4763 }
4764 
4765 template <>
4767 Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4768 {
4769  _need_curl.insert(type);
4770  buildVectorFaceNeighborFE(type);
4771 
4772  return _vector_fe_shape_data_face_neighbor[type]->_curl_phi;
4773 }
4774 
4775 template <>
4777 Assembly::feDivPhi<VectorValue<Real>>(FEType type) const
4778 {
4779  _need_div.insert(type);
4780  buildVectorFE(type);
4781  return _vector_fe_shape_data[type]->_div_phi;
4782 }
4783 
4784 template <>
4786 Assembly::feDivPhiFace<VectorValue<Real>>(FEType type) const
4787 {
4788  _need_face_div.insert(type);
4789  buildVectorFaceFE(type);
4790 
4791  // If we're building for a face we probably need to build for a
4792  // neighbor while _need_face_div is set;
4793  // onInterface/reinitNeighbor/etc don't distinguish
4794  buildVectorFaceNeighborFE(type);
4795 
4796  return _vector_fe_shape_data_face[type]->_div_phi;
4797 }
4798 
4799 template <>
4801 Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type) const
4802 {
4803  _need_neighbor_div.insert(type);
4804  buildVectorNeighborFE(type);
4805  return _vector_fe_shape_data_neighbor[type]->_div_phi;
4806 }
4807 
4808 template <>
4810 Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
4811 {
4812  _need_face_neighbor_div.insert(type);
4813  buildVectorFaceNeighborFE(type);
4814  return _vector_fe_shape_data_face_neighbor[type]->_div_phi;
4815 }
4816 
4817 const MooseArray<ADReal> &
4819 {
4820  _calculate_curvatures = true;
4821  const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
4822  const FEType helper_type(helper_order, LAGRANGE);
4823  // Must prerequest the second derivatives. Sadly because there is only one
4824  // _need_second_derivative map for both volumetric and face FE objects we must request both here
4825  feSecondPhi<Real>(helper_type);
4826  feSecondPhiFace<Real>(helper_type);
4827  return _ad_curvatures;
4828 }
4829 
4830 void
4832 {
4833  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
4834  {
4835  _holder_fe_helper[dim]->get_phi();
4836  _holder_fe_helper[dim]->get_dphi();
4837  _holder_fe_helper[dim]->get_xyz();
4838  _holder_fe_helper[dim]->get_JxW();
4839 
4840  _holder_fe_face_helper[dim]->get_phi();
4841  _holder_fe_face_helper[dim]->get_dphi();
4842  _holder_fe_face_helper[dim]->get_xyz();
4843  _holder_fe_face_helper[dim]->get_JxW();
4844  _holder_fe_face_helper[dim]->get_normals();
4845 
4848  _holder_fe_face_neighbor_helper[dim]->get_normals();
4849 
4850  _holder_fe_neighbor_helper[dim]->get_xyz();
4851  _holder_fe_neighbor_helper[dim]->get_JxW();
4852  }
4853 
4854  for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
4855  {
4856  // We need these computations in order to compute correct lower-d element volumes in
4857  // curvilinear coordinates
4858  _holder_fe_lower_helper[dim]->get_xyz();
4859  _holder_fe_lower_helper[dim]->get_JxW();
4860  }
4861 }
4862 
4863 void
4864 Assembly::havePRefinement(const std::unordered_set<FEFamily> & disable_families)
4865 {
4866  if (_have_p_refinement)
4867  // Already performed tasks for p-refinement
4868  return;
4869 
4870  const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
4871  const FEType helper_type(helper_order, LAGRANGE);
4872  auto process_fe =
4873  [&disable_families](const unsigned int num_dimensionalities, auto & fe_container)
4874  {
4875  if (!disable_families.empty())
4876  for (const auto dim : make_range(num_dimensionalities))
4877  {
4878  auto fe_container_it = fe_container.find(dim);
4879  if (fe_container_it != fe_container.end())
4880  for (auto & [fe_type, fe_ptr] : fe_container_it->second)
4881  if (disable_families.count(fe_type.family))
4882  fe_ptr->add_p_level_in_reinit(false);
4883  }
4884  };
4885  auto process_fe_and_helpers = [process_fe, &helper_type](auto & unique_helper_container,
4886  auto & helper_container,
4887  const unsigned int num_dimensionalities,
4888  const bool user_added_helper_type,
4889  auto & fe_container)
4890  {
4891  unique_helper_container.resize(num_dimensionalities);
4892  for (const auto dim : make_range(num_dimensionalities))
4893  {
4894  auto & unique_helper = unique_helper_container[dim];
4895  unique_helper = FEGenericBase<Real>::build(dim, helper_type);
4896  // don't participate in p-refinement
4897  unique_helper->add_p_level_in_reinit(false);
4898  helper_container[dim] = unique_helper.get();
4899 
4900  // If the user did not request the helper type then we should erase it from our FE container
4901  // so that they're not penalized (in the "we should be able to do p-refinement sense") for
4902  // our perhaps silly helpers
4903  if (!user_added_helper_type)
4904  {
4905  auto & fe_container_dim = libmesh_map_find(fe_container, dim);
4906  auto fe_it = fe_container_dim.find(helper_type);
4907  mooseAssert(fe_it != fe_container_dim.end(), "We should have the helper type");
4908  delete fe_it->second;
4909  fe_container_dim.erase(fe_it);
4910  }
4911  }
4912 
4913  process_fe(num_dimensionalities, fe_container);
4914  };
4915 
4916  // Handle scalar field families
4917  process_fe_and_helpers(_unique_fe_helper,
4919  _mesh_dimension + 1,
4921  _fe);
4922  process_fe_and_helpers(_unique_fe_face_helper,
4924  _mesh_dimension + 1,
4926  _fe_face);
4927  process_fe_and_helpers(_unique_fe_face_neighbor_helper,
4929  _mesh_dimension + 1,
4932  process_fe_and_helpers(_unique_fe_neighbor_helper,
4934  _mesh_dimension + 1,
4936  _fe_neighbor);
4937  process_fe_and_helpers(_unique_fe_lower_helper,
4941  _fe_lower);
4942  // Handle vector field families
4943  process_fe(_mesh_dimension + 1, _vector_fe);
4944  process_fe(_mesh_dimension + 1, _vector_fe_face);
4945  process_fe(_mesh_dimension + 1, _vector_fe_neighbor);
4946  process_fe(_mesh_dimension + 1, _vector_fe_face_neighbor);
4947  process_fe(_mesh_dimension, _vector_fe_lower);
4948 
4950 
4951  _have_p_refinement = true;
4952 }
4953 
4954 template void coordTransformFactor<Point, Real>(const SubProblem & s,
4955  SubdomainID sub_id,
4956  const Point & point,
4957  Real & factor,
4958  SubdomainID neighbor_sub_id);
4959 template void coordTransformFactor<ADPoint, ADReal>(const SubProblem & s,
4960  SubdomainID sub_id,
4961  const ADPoint & point,
4962  ADReal & factor,
4963  SubdomainID neighbor_sub_id);
4964 template void coordTransformFactor<Point, Real>(const MooseMesh & mesh,
4965  SubdomainID sub_id,
4966  const Point & point,
4967  Real & factor,
4968  SubdomainID neighbor_sub_id);
4970  SubdomainID sub_id,
4971  const ADPoint & point,
4972  ADReal & factor,
4973  SubdomainID neighbor_sub_id);
4974 
4975 template <>
4977 Assembly::genericQPoints<false>() const
4978 {
4979  return qPoints();
4980 }
4981 
4982 template <>
4984 Assembly::genericQPoints<true>() const
4985 {
4986  return adQPoints();
4987 }
const Elem *const & elem() const
Return the current element.
Definition: Assembly.h:375
virtual MooseMesh & mesh()=0
MooseArray< VectorValue< ADReal > > _ad_normals
Definition: Assembly.h:2817
void copyShapes(MooseVariableField< T > &v)
Definition: Assembly.C:3000
libMesh::ElemSideBuilder _compute_face_map_side_elem_builder
In place side element builder for computeFaceMap()
Definition: Assembly.h:2851
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face
types of finite elements
Definition: Assembly.h:2484
bool _need_neighbor_elem_volume
true is apps need to compute neighbor element volume
Definition: Assembly.h:2588
void cacheResidualNeighbor(GlobalDataKey, const std::vector< VectorTag > &tags)
Takes the values that are currently in _sub_Rn of all field variables and appends them to the cached ...
Definition: Assembly.C:3460
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
ArbitraryQuadrature * _current_qrule_arbitrary
The current arbitrary quadrature rule used within the element interior.
Definition: Assembly.h:2384
std::map< FEType, FEBase * > _current_fe
The "volume" fe object that matches the current elem.
Definition: Assembly.h:2352
MooseArray< Real > _curvatures
Definition: Assembly.h:2819
const std::vector< MooseVariableFieldBase * > & getVariables(THREAD_ID tid)
Definition: SystemBase.h:752
std::vector< ADReal > _ad_detadz_map
Definition: Assembly.h:2811
VectorVariablePhiValue _phi
Definition: Assembly.h:2727
std::unique_ptr< FEGenericBase< Real > > build(const unsigned int dim, const FEType &fet)
SystemBase & _sys
Definition: Assembly.h:2282
std::vector< std::vector< dof_id_type > > _cached_residual_rows
Where the cached values should go (the first vector is for TIME vs NONTIME)
Definition: Assembly.h:2771
std::vector< std::vector< dof_id_type > > _cached_jacobian_rows
Row where the corresponding cached value should go.
Definition: Assembly.h:2778
unsigned int _max_cached_residuals
Definition: Assembly.h:2773
std::map< FEType, ADTemplateVariablePhiGradient< Real > > _ad_grad_phi_data_face
Definition: Assembly.h:2752
std::vector< dof_id_type > componentDofIndices(const std::vector< dof_id_type > &dof_indices, unsigned int component) const
Obtain DoF indices of a component with the indices of the 0th component.
virtual_for_inffe const std::vector< std::vector< OutputDivergence > > & get_div_phi() const
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
bool _user_added_fe_lower_of_helper_type
Definition: Assembly.h:2335
MooseArray< Point > _current_physical_points
This will be filled up with the physical points passed into reinitAtPhysical() if it is called...
Definition: Assembly.h:2616
const VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField< Real > &) const
Definition: Assembly.h:1334
std::vector< std::unique_ptr< FEBase > > _unique_fe_lower_helper
Definition: Assembly.h:2343
Order
void print_info(std::ostream &os=libMesh::out) const
MooseArray< ADReal > _ad_coord
The AD version of the current coordinate transformation coefficients.
Definition: Assembly.h:2396
virtual const FieldVariablePhiValue & phi() const =0
Return the variable&#39;s elemental shape functions.
void setResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Set a local residual block to a global residual vector with proper scaling.
Definition: Assembly.C:3294
void buildNeighborFE(FEType type) const
Build FEs for a neighbor with a type.
Definition: Assembly.C:316
const VariablePhiSecond & secondPhi() const
Definition: Assembly.h:1299
void setMortarQRule(Order order)
Specifies a custom qrule for integration on mortar segment mesh.
Definition: Assembly.C:737
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe
Each dimension&#39;s actual vector fe objects indexed on type.
Definition: Assembly.h:2374
typename OutputTools< typename Moose::ADType< T >::type >::VariablePhiGradient ADTemplateVariablePhiGradient
Definition: MooseTypes.h:637
const std::vector< std::vector< OutputShape > > & get_dual_phi() const
void buildVectorLowerDFE(FEType type) const
Build Vector FEs for a lower dimensional element with a type.
Definition: Assembly.C:405
MooseArray< Real > _coord_neighbor
The current coordinate transformation coefficients.
Definition: Assembly.h:2541
virtual void haveADObjects(bool have_ad_objects)
Method for setting whether we have any ad objects.
Definition: SubProblem.h:767
void buildFE(FEType type) const
Build FEs with a type.
Definition: Assembly.C:268
const std::vector< MooseVariableScalar * > & getScalarVariables(THREAD_ID tid)
Definition: SystemBase.h:757
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kle
dlower/dsecondary (or dlower/delement)
Definition: Assembly.h:2661
void set_current_elem(const Elem &elem)
virtual void zero() override final
void cacheResidualLower(GlobalDataKey, const std::vector< VectorTag > &tags)
Takes the values that are currently in _sub_Rl and appends them to the cached values.
Definition: Assembly.C:3475
const unsigned int invalid_uint
void setFaceQRule(libMesh::QBase *qrule, unsigned int dim)
Set the qrule to be used for face integration.
Definition: Assembly.C:676
bool hasVector(const std::string &tag_name) const
Check if the named vector exists in the system.
Definition: SystemBase.C:907
LAGRANGE_VEC
Assembly(SystemBase &sys, THREAD_ID tid)
Definition: Assembly.C:80
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data
Shape function values, gradients, second derivatives for each FE type.
Definition: Assembly.h:2735
void buildLowerDDualFE(FEType type) const
Definition: Assembly.C:384
unsigned int TagID
Definition: MooseTypes.h:210
virtual bool checkNonlocalCouplingRequirement() const =0
MooseArray< Real > _current_JxW_neighbor
The current transformed jacobian weights on a neighbor&#39;s face.
Definition: Assembly.h:2539
std::shared_ptr< XFEMInterface > _xfem
The XFEM controller.
Definition: Assembly.h:2349
void reinitNeighborLowerDElem(const Elem *elem)
reinitialize a neighboring lower dimensional element
Definition: Assembly.C:2386
virtual ~Assembly()
Definition: Assembly.C:188
void prepareJacobianBlock()
Sizes and zeroes the Jacobian blocks used for the current element.
Definition: Assembly.C:2687
DenseMatrix< Number > & jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Get local Jacobian block from non-local contribution for a pair of variables and a tag...
Definition: Assembly.h:1124
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
virtual void zero() override final
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kel
dsecondary/dlower (or delement/dlower)
Definition: Assembly.h:2665
const std::vector< Real > & arrayScalingFactor() const
MooseArray< Real > _coord
The current coordinate transformation coefficients.
Definition: Assembly.h:2394
unsigned int number() const
Get variable number coming from libMesh.
bool allow_rules_with_negative_weights
DenseMatrix< Number > & jacobianBlockNeighbor(Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Get local Jacobian block of a DG Jacobian type for a pair of variables and a tag. ...
Definition: Assembly.C:3117
void reinitFE(const Elem *elem)
Just an internal helper function to reinit the volume FE objects.
Definition: Assembly.C:762
TagID _id
The id associated with the vector tag.
Definition: VectorTag.h:30
virtual unsigned int currentNlSysNum() const =0
libMesh::QBase * _current_qrule_neighbor
quadrature rule used on neighbors
Definition: Assembly.h:2533
void jacobianBlockUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not Jacobian coupling between ivar and jvar is used to the value used.
Definition: Assembly.h:2209
void coordTransformFactor(const P &point, C &factor, const Moose::CoordinateSystemType coord_type, const unsigned int rz_radial_coord=libMesh::invalid_uint)
compute a coordinate transformation factor
void prepareNonlocal()
Definition: Assembly.C:2725
std::map< unsigned int, FEBase * > _holder_fe_neighbor_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:2522
std::unique_ptr< FEBase > _fe_msm
A FE object for working on mortar segement elements.
Definition: Assembly.h:2551
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_face_neighbor
Definition: Assembly.h:2738
void createQRules(QuadratureType type, Order order, Order volume_order, Order face_order, SubdomainID block, bool allow_negative_qweights=true)
Creates block-specific volume, face and arbitrary qrules based on the orders and the flag of whether ...
Definition: Assembly.C:619
const std::vector< std::vector< Real > > & get_dphidzeta_map() const
const Elem & elem() const
Definition: FaceInfo.h:81
void reinitAtPhysical(const Elem *elem, const std::vector< Point > &physical_points)
Reinitialize the assembly data at specific physical point in the given element.
Definition: Assembly.C:1792
void addCachedResidualDirectly(NumericVector< Number > &residual, GlobalDataKey, const VectorTag &vector_tag)
Adds the values that have been cached by calling cacheResidual(), cacheResidualNeighbor(), and/or cacheResidualLower() to a user-defined residual (that is, not necessarily the vector that vector_tag points to)
Definition: Assembly.C:3535
bool _current_elem_volume_computed
Boolean to indicate whether current element volumes has been computed.
Definition: Assembly.h:2596
std::vector< ADReal > _ad_detady_map
Definition: Assembly.h:2810
void jacobianBlockLowerUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not lower Jacobian coupling between ivar and jvar is used to the value used...
Definition: Assembly.h:2245
const VectorVariablePhiDivergence & divPhi(const MooseVariableField< RealVectorValue > &) const
Definition: Assembly.h:1359
std::vector< ADReal > _ad_dzetadz_map
Definition: Assembly.h:2814
std::vector< std::unique_ptr< FEBase > > _unique_fe_neighbor_helper
Definition: Assembly.h:2342
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kll
dlower/dlower
Definition: Assembly.h:2659
bool _user_added_fe_face_neighbor_of_helper_type
Definition: Assembly.h:2333
void setResidualNeighbor(NumericVector< Number > &residual, GlobalDataKey, const VectorTag &vector_tag)
Sets local neighbor residuals of all field variables to the global residual vector for a tag...
Definition: Assembly.C:3566
bool _have_p_refinement
Whether we have ever conducted p-refinement.
Definition: Assembly.h:2871
VariablePhiValue _phi
Definition: Assembly.h:2717
char ** vars
Real _current_neighbor_volume
Volume of the current neighbor.
Definition: Assembly.h:2590
void jacobianBlockNonlocalUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not nonlocal Jacobian coupling between ivar and jvar is used to the value used...
Definition: Assembly.h:2263
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
ArbitraryQuadrature * qruleArbitraryFace(const Elem *elem, unsigned int side)
Definition: Assembly.C:1926
std::map< FEType, FEBase * > _current_fe_face
The "face" fe object that matches the current elem.
Definition: Assembly.h:2354
const VectorVariablePhiCurl & curlPhi(const MooseVariableField< RealVectorValue > &) const
Definition: Assembly.h:1355
std::vector< Point > _current_neighbor_ref_points
The current reference points on the neighbor element.
Definition: Assembly.h:2874
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kln
dlower/dprimary (or dlower/dneighbor)
Definition: Assembly.h:2663
void setVolumeQRule(libMesh::QBase *qrule, unsigned int dim)
Set the qrule to be used for volume integration.
Definition: Assembly.C:657
const Elem * _current_neighbor_elem
The current neighbor "element".
Definition: Assembly.h:2580
void addJacobianNonlocal(GlobalDataKey)
Adds non-local Jacobian to the global Jacobian matrices.
Definition: Assembly.C:3890
std::vector< std::vector< DenseVector< Number > > > _sub_Rn
Definition: Assembly.h:2631
unsigned int count() const
Get the number of components Note: For standard and vector variables, the number is one...
virtual const FieldVariablePhiGradient & gradPhiNeighbor() const =0
Return the gradients of the variable&#39;s shape functions on a neighboring element.
const VariablePhiValue & phi() const
Definition: Assembly.h:1290
MooseMesh & _mesh
Definition: Assembly.h:2322
TagTypeID _type_id
The index for this tag into a vector that contains tags of only its type ordered by ID...
Definition: VectorTag.h:47
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_neighbor
Definition: Assembly.h:2737
std::map< unsigned int, FEBase * > _holder_fe_face_neighbor_helper
Definition: Assembly.h:2523
unsigned int n_elem_integers() const
MeshBase & mesh
const VariablePhiGradient & gradPhi() const
Definition: Assembly.h:1297
MooseArray< Real > _current_JxW_face
The current transformed jacobian weights on a face.
Definition: Assembly.h:2498
void modifyWeightsDueToXFEM(const Elem *elem)
Update the integration weights for XFEM partial elements.
Definition: Assembly.C:4538
const std::vector< Real > & get_weights() const
void coordTransformFactorRZGeneral(const P &point, const std::pair< Point, RealVectorValue > &axis, C &factor)
Computes a coordinate transformation factor for a general axisymmetric axis.
const Elem * _current_elem
The current "element" we are currently on.
Definition: Assembly.h:2566
const VariablePhiSecond & secondPhiNeighbor(const MooseVariableField< Real > &) const
Definition: Assembly.h:1325
bool computingScalingJacobian() const
Whether we are computing an initial Jacobian for automatic variable scaling.
Definition: SystemBase.C:1519
THREAD_ID _tid
Thread number (id)
Definition: Assembly.h:2320
void reinitFEFaceNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:1573
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:153
virtual const FieldVariablePhiSecond & secondPhi() const =0
Return the rank-2 tensor of second derivatives of the variable&#39;s elemental shape functions.
void reinitNeighborAtPhysical(const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > &physical_points)
Reinitializes the neighbor at the physical coordinates on neighbor side given.
std::vector< std::unique_ptr< FEBase > > _unique_fe_face_helper
Definition: Assembly.h:2340
QuadratureType
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_nonlocal_used
Definition: Assembly.h:2312
std::vector< std::pair< MooseVariableScalar *, MooseVariableFieldBase * > > _cm_sf_entry
Entries in the coupling matrix for scalar variables vs field variables.
Definition: Assembly.h:2305
void prepareNeighbor()
Definition: Assembly.C:2804
void resizeADMappingObjects(unsigned int n_qp, unsigned int dim)
resize any objects that contribute to automatic differentiation-related mapping calculations ...
Definition: Assembly.C:973
OutputTools< Real >::VariablePhiValue VariablePhiValue
Definition: MooseTypes.h:319
unsigned int m() const
virtual bool hasMatrix(TagID tag) const
Check if the tagged matrix exists in the system.
Definition: SystemBase.h:351
const std::vector< std::vector< Real > > & get_phi_map() const
libMesh::StaticCondensation * _sc
A pointer to the static condensation class. Null if not present.
Definition: Assembly.h:2877
void setCoordinateTransformation(const libMesh::QBase *qrule, const Points &q_points, Coords &coord, SubdomainID sub_id)
Definition: Assembly.C:1731
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_face_neighbor
Definition: Assembly.h:2519
const VariablePhiValue & phiFaceNeighbor(const MooseVariableField< Real > &) const
Definition: Assembly.h:1330
unsigned int elemSideID() const
Definition: FaceInfo.h:109
std::vector< bool > _component_block_diagonal
An flag array Indiced by variable index to show if there is no component-wise coupling for the variab...
Definition: Assembly.h:2788
Real _current_elem_volume
Volume of the current element.
Definition: Assembly.h:2572
std::map< FEType, FEBase * > _current_fe_face_neighbor
The "neighbor face" fe object that matches the current elem.
Definition: Assembly.h:2358
const MooseArray< ADReal > & adCurvatures() const
Definition: Assembly.C:4818
This class provides an interface for common operations on field variables of both FE and FV types wit...
std::map< FEType, FEVectorBase * > _current_vector_fe_face
The "face" vector fe object that matches the current elem.
Definition: Assembly.h:2363
void shallowCopy(const MooseArray &rhs)
Doesn&#39;t actually make a copy of the data.
Definition: MooseArray.h:296
std::vector< std::unique_ptr< FEBase > > _unique_fe_helper
Containers for holding unique FE helper types if we are doing p-refinement.
Definition: Assembly.h:2339
void processLocalResidual(DenseVector< Number > &res_block, std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Appling scaling, constraints to the local residual block and populate the full DoF indices for array ...
Definition: Assembly.C:3218
const FEType _helper_type
The finite element type of the FE helper classes.
Definition: Assembly.h:2328
void addResidualLower(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add local neighbor residuals of all field variables for a set of tags onto the global residual vector...
Definition: Assembly.C:3377
DenseMatrix< Number > & jacobianBlock(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Get local Jacobian block for a pair of variables and a tag.
Definition: Assembly.h:1113
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
void addJacobianScalar(GlobalDataKey)
Add Jacobians for pairs of scalar variables into the global Jacobian matrices.
Definition: Assembly.C:4460
const std::vector< std::vector< OutputShape > > & get_dphideta() const
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Knn
jacobian contributions from the neighbor <Tag, ivar, jvar>
Definition: Assembly.h:2657
Base class for a system (of equations)
Definition: SystemBase.h:84
libMesh::QBase * _current_qrule_face
quadrature rule used on faces
Definition: Assembly.h:2492
const VariablePhiValue & phiNeighbor(const MooseVariableField< Real > &) const
Definition: Assembly.h:1317
unsigned int numExtraElemIntegers() const
Number of extra element integers Assembly tracked.
Definition: Assembly.h:334
std::set< FEType > _need_second_derivative_neighbor
Definition: Assembly.h:2836
void setWeights(const std::vector< libMesh::Real > &weights)
Set the quadrature weights.
DualNumber< Real, DNDerivativeType, true > ADReal
Definition: ADRealForward.h:46
std::vector< std::vector< Real > > _cached_jacobian_values
Values cached by calling cacheJacobian()
Definition: Assembly.h:2776
std::vector< std::pair< MooseVariableFieldBase *, MooseVariableFieldBase * > > _cm_ff_entry
Entries in the coupling matrix for field variables.
Definition: Assembly.h:2301
void addJacobianNeighborLowerD(GlobalDataKey)
Add all portions of the Jacobian except PrimaryPrimary, e.g.
Definition: Assembly.C:3950
std::map< FEType, FEVectorBase * > _current_vector_fe
The "volume" vector fe object that matches the current elem.
Definition: Assembly.h:2361
VariablePhiGradient _grad_phi
Definition: Assembly.h:2718
template void coordTransformFactor< ADPoint, ADReal >(const SubProblem &s, SubdomainID sub_id, const ADPoint &point, ADReal &factor, SubdomainID neighbor_sub_id)
std::vector< VectorValue< ADReal > > _ad_dxyzdeta_map
Definition: Assembly.h:2798
void prepareScalar()
Definition: Assembly.C:2950
virtual const FieldVariablePhiValue & phiNeighbor() const =0
Return the variable&#39;s shape functions on a neighboring element.
void modifyFaceWeightsDueToXFEM(const Elem *elem, unsigned int side=0)
Update the face integration weights for XFEM partial elements.
Definition: Assembly.C:4558
void reinitElemAndNeighbor(const Elem *elem, unsigned int side, const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > *neighbor_reference_points=nullptr)
Reinitialize an element and its neighbor along a particular side.
Definition: Assembly.C:1997
OutputTools< Real >::VariablePhiSecond VariablePhiSecond
Definition: MooseTypes.h:321
std::unique_ptr< libMesh::QBase > face
area/face (meshdim-1) quadrature rule
Definition: Assembly.h:2414
unsigned int neighborSideID() const
Definition: FaceInfo.h:110
DenseVector< Number > _tmp_Re
auxiliary vector for scaling residuals (optimization to avoid expensive construction/destruction) ...
Definition: Assembly.h:2636
unsigned int _mesh_dimension
Definition: Assembly.h:2324
virtual QuadratureType type() const=0
std::vector< Eigen::Map< RealDIMValue > > _mapped_normals
Mapped normals.
Definition: Assembly.h:2502
VectorVariablePhiDivergence _vector_div_phi_face
Definition: Assembly.h:2700
void cacheJacobian(GlobalDataKey)
Takes the values that are currently in _sub_Kee and appends them to the cached values.
Definition: Assembly.C:4065
void reinit(const Elem *elem)
Reinitialize objects (JxW, q_points, ...) for an elements.
const NumericVector< Real > * _scaling_vector
The map from global index to variable scaling factor.
Definition: Assembly.h:2844
auto max(const L &left, const R &right)
const std::vector< std::vector< OutputShape > > & get_dphidzeta() const
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value)=0
void addCachedResiduals(GlobalDataKey, const std::vector< VectorTag > &tags)
Pushes all cached residuals to the global residual vectors associated with each tag.
Definition: Assembly.C:3490
Data structure for tracking/grouping a set of quadrature rules for a particular dimensionality of mes...
Definition: Assembly.h:2400
MONOMIAL_VEC
const std::vector< std::vector< Real > > & get_dphideta_map() const
unsigned int _max_cached_jacobians
Definition: Assembly.h:2782
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
DenseMatrix< Number > & jacobianBlockMortar(Moose::ConstraintJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Returns the jacobian block for the given mortar Jacobian type.
Definition: Assembly.C:3158
QUAD4
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Ken
jacobian contributions from the element and neighbor <Tag, ivar, jvar>
Definition: Assembly.h:2653
FEType get_fe_type() const
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_dual_lower
Definition: Assembly.h:2748
virtual unsigned int nVariables() const
Get the number of variables in this system.
Definition: SystemBase.C:874
virtual const FieldVariablePhiValue & phiFaceNeighbor() const =0
Return the variable&#39;s shape functions on a neighboring element face.
unsigned int n_dofs(const unsigned int s, const unsigned int var=libMesh::invalid_uint) const
std::vector< VectorValue< ADReal > > _ad_dxyzdzeta_map
Definition: Assembly.h:2799
std::vector< Point > _temp_reference_points
Temporary work data for reinitAtPhysical()
Definition: Assembly.h:2794
std::vector< std::pair< MooseVariableScalar *, MooseVariableScalar * > > _cm_ss_entry
Entries in the coupling matrix for scalar variables.
Definition: Assembly.h:2307
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kne
jacobian contributions from the neighbor and element <Tag, ivar, jvar>
Definition: Assembly.h:2655
const Elem * _current_neighbor_lower_d_elem
The current neighboring lower dimensional element.
Definition: Assembly.h:2603
unsigned int _current_neighbor_side
The current side of the selected neighboring element (valid only when working with sides) ...
Definition: Assembly.h:2584
void cacheResidualNodes(const DenseVector< Number > &res, const std::vector< dof_id_type > &dof_index, LocalDataKey, TagID tag)
Lets an external class cache residual at a set of nodes.
Definition: Assembly.C:3442
void init(const libMesh::CouplingMatrix *cm)
Initialize the Assembly object and set the CouplingMatrix for use throughout.
Definition: Assembly.C:2470
void prepareLowerD()
Prepare the Jacobians and residuals for a lower dimensional element.
Definition: Assembly.C:2843
bool _user_added_fe_neighbor_of_helper_type
Definition: Assembly.h:2334
const std::vector< std::vector< OutputGradient > > & get_dphi() const
void buildFaceNeighborFE(FEType type) const
Build FEs for a neighbor face with a type.
Definition: Assembly.C:338
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_face
Definition: Assembly.h:2744
void reinitLowerDElem(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)
Reinitialize FE data for a lower dimenesional element with a given set of reference points...
Definition: Assembly.C:2293
void computeCurrentFaceVolume()
Definition: Assembly.C:1773
std::vector< std::vector< DenseVector< Number > > > _sub_Rl
residual contributions for each variable from the lower dimensional element
Definition: Assembly.h:2633
This data structure is used to store geometric and variable related metadata about each cell face in ...
Definition: FaceInfo.h:36
QRules & qrules(unsigned int dim)
Definition: Assembly.h:2460
virtual const FieldVariablePhiSecond & secondPhiFaceNeighbor() const =0
Return the rank-2 tensor of second derivatives of the variable&#39;s shape functions on a neighboring ele...
virtual const std::vector< dof_id_type > & dofIndicesNeighbor() const =0
Get neighbor DOF indices for currently selected element.
static const subdomain_id_type invalid_subdomain_id
CONSTANT
std::unique_ptr< libMesh::QBase > vol
volume/elem (meshdim) quadrature rule
Definition: Assembly.h:2412
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
void prepareOffDiagScalar()
Definition: Assembly.C:2974
VectorVariablePhiSecond _second_phi
Definition: Assembly.h:2729
unsigned int size() const
The number of elements that can currently be stored in the array.
Definition: MooseArray.h:259
Implements a fake quadrature rule where you can specify the locations (in the reference domain) of th...
void addJacobianCoupledVarPair(const MooseVariableBase &ivar, const MooseVariableBase &jvar)
Adds element matrices for ivar rows and jvar columns to the global Jacobian matrices.
Definition: Assembly.C:3862
std::vector< ADReal > _ad_dzetadx_map
Definition: Assembly.h:2812
void addJacobianLowerD(GlobalDataKey)
Add portions of the Jacobian of LowerLower, LowerSecondary, and SecondaryLower for boundary condition...
Definition: Assembly.C:4028
virtual const FieldVariablePhiSecond & secondPhiFace() const =0
Return the rank-2 tensor of second derivatives of the variable&#39;s shape functions on an element face...
void setPoints(const std::vector< libMesh::Point > &points)
Set the quadrature points.
const Elem * neighborPtr() const
Definition: FaceInfo.h:84
VectorVariablePhiGradient _grad_phi
Definition: Assembly.h:2728
MooseArray< ADReal > _ad_curvatures
Definition: Assembly.h:2820
void computeGradPhiAD(const Elem *elem, unsigned int n_qp, ADTemplateVariablePhiGradient< OutputType > &grad_phi, libMesh::FEGenericBase< OutputType > *fe)
compute gradient of phi possibly with derivative information with respect to nonlinear displacement v...
ArbitraryQuadrature * _current_qrule_arbitrary_face
The current arbitrary quadrature rule used on the element face.
Definition: Assembly.h:2386
void modifyArbitraryWeights(const std::vector< Real > &weights)
Modify the weights when using the arbitrary quadrature rule.
Definition: Assembly.C:4585
std::vector< std::pair< MooseVariableFieldBase *, MooseVariableScalar * > > _cm_fs_entry
Entries in the coupling matrix for field variables vs scalar variables.
Definition: Assembly.h:2303
const Node & node_ref(const unsigned int i) const
const std::vector< Real > * _JxW_msm
A JxW for working on mortar segement elements.
Definition: Assembly.h:2549
dof_id_type id() const
std::vector< ADReal > _ad_dzetady_map
Definition: Assembly.h:2813
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3443
void coordTransformFactor(const SubProblem &s, const SubdomainID sub_id, const P &point, C &factor, const SubdomainID neighbor_sub_id)
Computes a conversion multiplier for use when computing integraals for the current coordinate system ...
Definition: Assembly.C:42
dof_id_type numeric_index_type
std::map< FEType, ADTemplateVariablePhiGradient< Real > > _ad_grad_phi_data
Definition: Assembly.h:2750
template void coordTransformFactor< Point, Real >(const SubProblem &s, SubdomainID sub_id, const Point &point, Real &factor, SubdomainID neighbor_sub_id)
std::unordered_map< SubdomainID, std::vector< QRules > > _qrules
Holds quadrature rules for each dimension.
Definition: Assembly.h:2429
unsigned int n_vars
void prepareResidual()
Sizes and zeroes the residual for the current element.
Definition: Assembly.C:2709
VectorVariablePhiCurl _curl_phi
Definition: Assembly.h:2730
std::map< FEType, FEBase * > _current_fe_neighbor
The "neighbor" fe object that matches the current elem.
Definition: Assembly.h:2356
Real elementVolume(const Elem *elem) const
On-demand computation of volume element accounting for RZ/RSpherical.
Definition: Assembly.C:3787
void cacheJacobianNeighbor(GlobalDataKey)
Takes the values that are currently in the neighbor Dense Matrices and appends them to the cached val...
Definition: Assembly.C:4117
SubdomainID _current_neighbor_subdomain_id
The current neighbor subdomain ID.
Definition: Assembly.h:2582
const std::vector< std::vector< OutputGradient > > & get_dual_dphi() const
std::vector< std::vector< DenseVector< Number > > > _sub_Re
Definition: Assembly.h:2630
void reinitNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Reinitializes the neighbor side using reference coordinates.
Definition: Assembly.C:1687
std::vector< ADReal > _ad_dxidz_map
Definition: Assembly.h:2808
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_used
Flag that indicates if the jacobian block was used.
Definition: Assembly.h:2311
void prepare()
Definition: Assembly.C:2718
virtual numeric_index_type m() const=0
const Node *const * get_nodes() const
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:88
libMesh::QBase * _current_qrule_lower
quadrature rule used on lower dimensional elements.
Definition: Assembly.h:2562
SubProblem & _subproblem
Definition: Assembly.h:2283
libmesh_assert(ctx)
std::vector< ADReal > _ad_detadx_map
Definition: Assembly.h:2809
void buildFaceFE(FEType type) const
Build FEs for a face with a type.
Definition: Assembly.C:294
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_neighbor
types of finite elements
Definition: Assembly.h:2516
void cacheResidual(GlobalDataKey, const std::vector< VectorTag > &tags)
Takes the values that are currently in _sub_Re of all field variables and appends them to the cached ...
Definition: Assembly.C:3409
bool _calculate_xyz
Definition: Assembly.h:2827
void addJacobianBlockNonlocalTags(libMesh::SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const libMesh::DofMap &dof_map, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, GlobalDataKey, const std::set< TagID > &tags)
Adds non-local element matrix for ivar rows and jvar columns to the global Jacobian matrix...
Definition: Assembly.C:4361
MooseArray< std::vector< Point > > _current_tangents
The current tangent vectors at the quadrature points.
Definition: Assembly.h:2504
virtual const std::vector< dof_id_type > & dofIndices() const
Get local DoF indices.
bool _calculate_curvatures
Definition: Assembly.h:2829
const std::vector< VectorTag > & _residual_vector_tags
The residual vector tags that Assembly could possibly contribute to.
Definition: Assembly.h:2765
virtual_for_inffe const std::vector< Real > & get_JxW() const
std::map< FEType, FEVectorBase * > _current_vector_fe_face_neighbor
The "neighbor face" vector fe object that matches the current elem.
Definition: Assembly.h:2367
std::vector< dof_id_type > _neighbor_extra_elem_ids
Extra element IDs of neighbor.
Definition: Assembly.h:2509
MooseArray< Real > _current_JxW
The current list of transformed jacobian weights.
Definition: Assembly.h:2390
Order get_order() const
std::vector< VectorValue< ADReal > > _ad_dxyzdxi_map
AD quantities.
Definition: Assembly.h:2797
virtual void zero_rows(std::vector< numeric_index_type > &rows, T diag_value=0.0)
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
std::vector< std::pair< unsigned int, unsigned short > > _disp_numbers_and_directions
Container of displacement numbers and directions.
Definition: Assembly.h:2825
std::vector< dof_id_type > _temp_dof_indices
Temporary work vector to keep from reallocating it.
Definition: Assembly.h:2791
std::set< FEType > _need_neighbor_div
Definition: Assembly.h:2840
L2_RAVIART_THOMAS
OutputTools< Real >::VariablePhiCurl VariablePhiCurl
Definition: MooseTypes.h:322
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kee
Definition: Assembly.h:2649
std::vector< VectorValue< ADReal > > _ad_d2xyzdxi2_map
Definition: Assembly.h:2800
unsigned int get_dim() const
std::vector< ADReal > _ad_jac
Definition: Assembly.h:2803
std::set< FEType > _need_curl
Definition: Assembly.h:2837
void addJacobianBlockNonlocal(libMesh::SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const libMesh::DofMap &dof_map, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, GlobalDataKey, TagID tag)
Adds non-local element matrix for ivar rows and jvar columns to the global Jacobian matrix...
Definition: Assembly.C:4303
unsigned int n_points() const
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_lower
Definition: Assembly.h:2747
OStreamProxy err(std::cerr)
const VariablePhiSecond & secondPhiFace(const MooseVariableField< Real > &) const
Definition: Assembly.h:1312
std::map< FEType, ADTemplateVariablePhiGradient< RealVectorValue > > _ad_vector_grad_phi_data
Definition: Assembly.h:2751
virtual_for_inffe const std::vector< Point > & get_xyz() const
bool usesPhiNeighbor() const
Whether or not this variable is actually using the shape function value.
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_face
types of vector finite elements
Definition: Assembly.h:2486
libMesh::QBase * _current_qrule_volume
The current volumetric quadrature for the element.
Definition: Assembly.h:2382
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
virtual bool computingSecond() const =0
Whether or not this variable is computing any second derivatives.
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1159
const Elem * _msm_elem
Definition: Assembly.h:2853
std::vector< T > stdVector() const
Extremely inefficient way to produce a std::vector from a MooseArray!
Definition: MooseArray.h:344
libMesh::QBase * qruleFace(const Elem *elem, unsigned int side)
This is an abstraction over the internal qrules function.
Definition: Assembly.C:1920
bool _calculate_ad_coord
Whether to calculate coord with AD.
Definition: Assembly.h:2833
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_lower
Definition: Assembly.h:2739
const Elem * _current_lower_d_elem
The current lower dimensional element.
Definition: Assembly.h:2601
void addResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Add a local residual block to a global residual vector with proper scaling.
Definition: Assembly.C:3254
void addResidual(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add local residuals of all field variables for a set of tags onto the global residual vectors associa...
Definition: Assembly.C:3327
const VariablePhiValue & phiFace() const
Definition: Assembly.h:1305
std::set< FEType > _need_second_derivative
Definition: Assembly.h:2835
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_lower
FE objects for lower dimensional elements.
Definition: Assembly.h:2526
bool hasSecondOrderElements()
check if the mesh has SECOND order elements
Definition: MooseMesh.C:3706
std::set< FEType > _need_face_neighbor_div
Definition: Assembly.h:2841
libMesh::QBase * _qrule_msm
A qrule object for working on mortar segement elements.
Definition: Assembly.h:2556
libMesh::ElemSideBuilder _current_neighbor_side_elem_builder
In place side element builder for _current_neighbor_side_elem.
Definition: Assembly.h:2849
OutputTools< Real >::VariablePhiDivergence VariablePhiDivergence
Definition: MooseTypes.h:323
Real _current_lower_d_elem_volume
The current lower dimensional element volume.
Definition: Assembly.h:2607
libMesh::QBase * _current_qrule
The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac k...
Definition: Assembly.h:2380
std::map< unsigned int, std::map< FEType, FEBase * > > _fe
Each dimension&#39;s actual fe objects indexed on type.
Definition: Assembly.h:2372
const MooseArray< Real > & JxW() const
Returns the reference to the transformed jacobian weights.
Definition: Assembly.h:248
void reinitMortarElem(const Elem *elem)
reinitialize a mortar segment mesh element in order to get a proper JxW
Definition: Assembly.C:2407
Real _current_side_volume
Volume of the current side element.
Definition: Assembly.h:2578
std::set< FEType > _need_div
Definition: Assembly.h:2838
std::vector< std::vector< Real > > _cached_residual_values
Values cached by calling cacheResidual() (the first vector is for TIME vs NONTIME) ...
Definition: Assembly.h:2768
const bool _displaced
Definition: Assembly.h:2285
virtual bool usesSecondPhiNeighbor() const =0
Whether or not this variable is actually using the shape function second derivatives.
virtual const FieldVariablePhiSecond & secondPhiNeighbor() const =0
Return the rank-2 tensor of second derivatives of the variable&#39;s shape functions on a neighboring ele...
void buildVectorFE(FEType type) const
Build Vector FEs with a type.
Definition: Assembly.C:455
virtual const FieldVariablePhiGradient & gradPhiFaceNeighbor() const =0
Return the gradients of the variable&#39;s shape functions on a neighboring element face.
void addJacobianNeighborTags(libMesh::SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const libMesh::DofMap &dof_map, std::vector< dof_id_type > &dof_indices, std::vector< dof_id_type > &neighbor_dof_indices, GlobalDataKey, const std::set< TagID > &tags)
Adds three neighboring element matrices for ivar rows and jvar columns to the global Jacobian matrix...
Definition: Assembly.C:4445
void reinitFVFace(const FaceInfo &fi)
Definition: Assembly.C:1860
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_lower
Vector FE objects for lower dimensional elements.
Definition: Assembly.h:2528
const MooseArray< Real > & JxWNeighbor() const
Returns the reference to the transformed jacobian weights on a current face.
Definition: Assembly.C:261
const Node * _current_node
The current node we are working with.
Definition: Assembly.h:2592
bool _building_helpers
Whether we are currently building the FE classes for the helpers.
Definition: Assembly.h:2346
bool _calculate_face_xyz
Definition: Assembly.h:2828
virtual unsigned int numMatrixTags() const
The total number of tags.
Definition: SubProblem.h:248
DGJacobianType
Definition: MooseTypes.h:749
OutputTools< Real >::VariablePhiGradient VariablePhiGradient
Definition: MooseTypes.h:320
std::vector< ADReal > _ad_dxidx_map
Definition: Assembly.h:2806
virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a scalar variable with specified number.
Definition: SystemBase.C:144
std::vector< ADReal > _ad_dxidy_map
Definition: Assembly.h:2807
std::vector< VectorValue< ADReal > > _ad_d2xyzdxideta_map
Definition: Assembly.h:2801
void set_calculate_default_dual_coeff(const bool val)
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_face_neighbor
Definition: Assembly.h:2746
void addJacobian(GlobalDataKey)
Adds all local Jacobian to the global Jacobian matrices.
Definition: Assembly.C:3877
MooseArray< Point > _current_normals
The current Normal vectors at the quadrature points.
Definition: Assembly.h:2500
virtual const FieldVariablePhiValue & phiFace() const =0
Return the variable&#39;s shape functions on an element face.
void buildVectorFaceFE(FEType type) const
Build Vector FEs for a face with a type.
Definition: Assembly.C:486
void havePRefinement(const std::unordered_set< FEFamily > &disable_p_refinement_for_families)
Indicate that we have p-refinement.
Definition: Assembly.C:4864
const SubdomainID ANY_BLOCK_ID
Definition: MooseTypes.C:19
MooseArray< VectorValue< ADReal > > _ad_q_points
Definition: Assembly.h:2805
const VariablePhiGradient & gradPhiNeighbor(const MooseVariableField< Real > &) const
Definition: Assembly.h:1321
MooseArray< VectorValue< ADReal > > _ad_q_points_face
Definition: Assembly.h:2818
virtual const FieldVariablePhiGradient & gradPhi() const =0
Return the gradients of the variable&#39;s elemental shape functions.
RAVIART_THOMAS
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const libMesh::CouplingMatrix * _cm
Coupling matrices.
Definition: Assembly.h:2288
void hasScalingVector()
signals this object that a vector containing variable scaling factors should be used when doing resid...
Definition: Assembly.C:4579
void addJacobianBlock(libMesh::SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const libMesh::DofMap &dof_map, std::vector< dof_id_type > &dof_indices, GlobalDataKey, TagID tag)
Adds element matrix for ivar rows and jvar columns to the global Jacobian matrix. ...
Generic class for solving transient nonlinear problems.
Definition: SubProblem.h:78
EDGE2
subdomain_id_type subdomain_id() const
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_dual_lower
Definition: Assembly.h:2740
void zeroCachedJacobian(GlobalDataKey)
Zero out previously-cached Jacobian rows.
Definition: Assembly.C:4517
void reinitNeighborFaceRef(const Elem *neighbor_elem, unsigned int neighbor_side, Real tolerance, const std::vector< Point > *const pts, const std::vector< Real > *const weights=nullptr)
Reinitialize FE data for the given neighbor_element on the given side with a given set of reference p...
Definition: Assembly.C:2197
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sqrt(_arg)) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tanh
void setNeighborQRule(libMesh::QBase *qrule, unsigned int dim)
Set the qrule to be used for neighbor integration.
Definition: Assembly.C:711
const std::vector< Point > & get_points() const
virtual unsigned short dim() const=0
bool usesGradPhiNeighbor() const
Whether or not this variable is actually using the shape function gradient.
void prepareVariableNonlocal(MooseVariableFieldBase *var)
Definition: Assembly.C:2777
void release()
Manually deallocates the data pointer.
Definition: MooseArray.h:66
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_face
Definition: Assembly.h:2736
virtual libMesh::SparseMatrix< Number > & getMatrix(TagID tag)
Get a raw SparseMatrix.
Definition: SystemBase.C:1007
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data
Shape function values, gradients, second derivatives for each vector FE type.
Definition: Assembly.h:2743
void addResidualScalar(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add residuals of all scalar variables for a set of tags onto the global residual vectors associated w...
Definition: Assembly.C:3401
MooseArray< Real > _coord_msm
The coordinate transformation coefficients evaluated on the quadrature points of the mortar segment m...
Definition: Assembly.h:2544
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Keg
Definition: Assembly.h:2650
bool _need_lower_d_elem_volume
Whether we need to compute the lower dimensional element volume.
Definition: Assembly.h:2605
bool _user_added_fe_face_of_helper_type
Definition: Assembly.h:2332
const std::vector< dof_id_type > & allDofIndices() const
Get all global dofindices for the variable.
void buildVectorNeighborFE(FEType type) const
Build Vector FEs for a neighbor with a type.
Definition: Assembly.C:516
std::unique_ptr< ArbitraryQuadrature > arbitrary_vol
volume/elem (meshdim) custom points quadrature rule
Definition: Assembly.h:2418
void reinitFENeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:1632
void initNonlocalCoupling()
Create pair of variables requiring nonlocal jacobian contributions.
Definition: Assembly.C:2652
ConstraintJacobianType
Definition: MooseTypes.h:796
std::map< unsigned int, FEBase * > _holder_fe_face_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:2488
const std::vector< std::vector< Real > > & get_dphidxi_map() const
const Elem *const & neighbor() const
Return the neighbor element.
Definition: Assembly.h:431
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_neighbor
Definition: Assembly.h:2518
virtual Real volume() const
void cacheJacobianNonlocal(GlobalDataKey)
Takes the values that are currently in _sub_Keg and appends them to the cached values.
Definition: Assembly.C:4095
Class for scalar variables (they are different).
IntRange< T > make_range(T beg, T end)
std::vector< std::unique_ptr< FEBase > > _unique_fe_face_neighbor_helper
Definition: Assembly.h:2341
void resize(unsigned int size)
Change the number of elements the array can store.
Definition: MooseArray.h:216
void computeSinglePointMapAD(const Elem *elem, const std::vector< Real > &qw, unsigned p, FEBase *fe)
compute the finite element reference-physical mapping quantities (such as JxW) with possible dependen...
Definition: Assembly.C:1003
void helpersRequestData()
request phi, dphi, xyz, JxW, etc.
Definition: Assembly.C:4831
std::vector< std::pair< MooseVariableFieldBase *, MooseVariableFieldBase * > > _cm_nonlocal_entry
Entries in the coupling matrix for field variables for nonlocal calculations.
Definition: Assembly.h:2309
void addJacobianOffDiagScalar(unsigned int ivar, GlobalDataKey)
Add Jacobians for a scalar variables with all other field variables into the global Jacobian matrices...
Definition: Assembly.C:4467
virtual void init(const Elem &e, unsigned int p_level=invalid_uint)
void buildLowerDFE(FEType type) const
Build FEs for a lower dimensional element with a type.
Definition: Assembly.C:360
virtual unsigned int size() const override final
unsigned int _current_side
The current side of the selected element (valid only when working with sides)
Definition: Assembly.h:2574
void setCachedJacobian(GlobalDataKey)
Sets previously-cached Jacobian values via SparseMatrix::set() calls.
Definition: Assembly.C:4497
VariablePhiSecond _second_phi
Definition: Assembly.h:2719
const libMesh::CouplingMatrix & _nonlocal_cm
Definition: Assembly.h:2289
std::map< unsigned int, FEBase * > _holder_fe_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:2376
void prepareVariable(MooseVariableFieldBase *var)
Used for preparing the dense residual and jacobian blocks for one particular variable.
Definition: Assembly.C:2749
const Elem * _current_neighbor_side_elem
The current side element of the ncurrent neighbor element.
Definition: Assembly.h:2586
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
Definition: ADReal.h:21
void reinitElemFaceRef(const Elem *elem, unsigned int elem_side, Real tolerance, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)
Reinitialize FE data for the given element on the given side, optionally with a given set of referenc...
Definition: Assembly.C:2022
void copyFaceShapes(MooseVariableField< T > &v)
Definition: Assembly.C:3037
std::vector< dof_id_type > _extra_elem_ids
Extra element IDs.
Definition: Assembly.h:2507
Moose::CoordinateSystemType getCoordSystem(SubdomainID sid) const
Definition: SubProblem.C:1278
void clearCachedQRules()
Set the cached quadrature rules to nullptr.
Definition: Assembly.C:728
void addJacobianBlockTags(libMesh::SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const libMesh::DofMap &dof_map, std::vector< dof_id_type > &dof_indices, GlobalDataKey, const std::set< TagID > &tags)
Add element matrix for ivar rows and jvar columns to the global Jacobian matrix for given tags...
Definition: Assembly.C:4235
void addResidualNeighbor(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add local neighbor residuals of all field variables for a set of tags onto the global residual vector...
Definition: Assembly.C:3352
bool _current_side_volume_computed
Boolean to indicate whether current element side volumes has been computed.
Definition: Assembly.h:2598
virtual bool isScalarVariable(unsigned int var_name) const
Definition: SystemBase.C:868
virtual void reinit_dual_shape_coeffs(const Elem *, const std::vector< Point > &, const std::vector< Real > &)
std::map< FEType, ADTemplateVariablePhiGradient< RealVectorValue > > _ad_vector_grad_phi_data_face
Definition: Assembly.h:2754
Moose::VectorTagType _type
The type of the vector tag.
Definition: VectorTag.h:53
void jacobianBlockNeighborUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not neighbor Jacobian coupling between ivar and jvar is used to the value used...
Definition: Assembly.h:2227
std::unique_ptr< ArbitraryQuadrature > arbitrary_face
area/face (meshdim-1) custom points quadrature rule
Definition: Assembly.h:2420
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
Definition: MooseTypes.h:146
virtual const std::vector< dof_id_type > & dofIndicesLower() const =0
Get dof indices for the current lower dimensional element (this is meaningful when performing mortar ...
void prepareBlock(unsigned int ivar, unsigned jvar, const std::vector< dof_id_type > &dof_indices)
Definition: Assembly.C:2899
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face_neighbor
Definition: Assembly.h:2517
void addJacobianNeighbor(GlobalDataKey)
Add ElementNeighbor, NeighborElement, and NeighborNeighbor portions of the Jacobian for compute objec...
Definition: Assembly.C:3912
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_lower_used
Flag that indicates if the jacobian block for the lower dimensional element was used.
Definition: Assembly.h:2316
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
VectorVariablePhiDivergence _div_phi
Definition: Assembly.h:2731
void cacheJacobianBlock(DenseMatrix< Number > &jac_block, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, Real scaling_factor, LocalDataKey, TagID tag)
Cache a local Jacobian block with the provided rows (idof_indices) and columns (jdof_indices) for eve...
const std::vector< std::vector< OutputShape > > & get_dphidxi() const
MooseVariableFieldBase & getVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:89
NEDELEC_ONE
bool _custom_mortar_qrule
Flag specifying whether a custom quadrature rule has been specified for mortar segment mesh...
Definition: Assembly.h:2558
const unsigned int & side() const
Returns the current side.
Definition: Assembly.h:407
void computeCurrentElemVolume()
Definition: Assembly.C:1754
void cacheJacobianBlockNonzero(DenseMatrix< Number > &jac_block, const MooseVariableBase &ivar, const MooseVariableBase &jvar, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, TagID tag)
Push non-zeros of a local Jacobian block with proper scaling into cache for a certain tag...
Definition: Assembly.C:3694
Real _current_neighbor_lower_d_elem_volume
The current neighboring lower dimensional element volume.
Definition: Assembly.h:2611
const FEMap & get_fe_map() const
bool _need_neighbor_lower_d_elem_volume
Whether we need to compute the neighboring lower dimensional element volume.
Definition: Assembly.h:2609
MooseArray< Point > _current_q_points
The current list of quadrature points.
Definition: Assembly.h:2388
bool _user_added_fe_of_helper_type
Whether user code requested a FEType the same as our _helper_type.
Definition: Assembly.h:2331
Moose::CoordinateSystemType _coord_type
The coordinate system.
Definition: Assembly.h:2392
void setLowerQRule(libMesh::QBase *qrule, unsigned int dim)
Set the qrule to be used for lower dimensional integration.
Definition: Assembly.C:692
void cacheJacobianMortar(GlobalDataKey)
Cache all portions of the Jacobian, e.g.
Definition: Assembly.C:4155
std::unique_ptr< ArbitraryQuadrature > neighbor
area/face (meshdim-1) custom points quadrature rule for DG
Definition: Assembly.h:2422
unsigned int n() const
void buildVectorFaceNeighborFE(FEType type) const
Build Vector FEs for a neighbor face with a type.
Definition: Assembly.C:546
void buildVectorDualLowerDFE(FEType type) const
Definition: Assembly.C:430
void bumpVolumeQRuleOrder(Order volume_order, SubdomainID block)
Increases the element/volume quadrature order for the specified mesh block if and only if the current...
Definition: Assembly.C:577
Storage for all of the information pretaining to a vector tag.
Definition: VectorTag.h:17
const Node *const & node() const
Returns the reference to the node.
Definition: Assembly.h:506
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_neighbor_used
Flag that indicates if the jacobian block for neighbor was used.
Definition: Assembly.h:2314
std::map< unsigned int, FEBase * > _holder_fe_lower_helper
helper object for transforming coordinates for lower dimensional element quadrature points ...
Definition: Assembly.h:2530
const libMesh::DofMap & _dof_map
DOF map.
Definition: Assembly.h:2318
virtual Order default_order() const=0
void computeADFace(const Elem &elem, const unsigned int side)
compute AD things on an element face
Definition: Assembly.C:2114
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_neighbor
Definition: Assembly.h:2745
void setResidual(NumericVector< Number > &residual, GlobalDataKey, const VectorTag &vector_tag)
Sets local residuals of all field variables to the global residual vector for a tag.
Definition: Assembly.C:3553
std::vector< std::vector< dof_id_type > > _cached_jacobian_cols
Column where the corresponding cached value should go.
Definition: Assembly.h:2780
void addCachedJacobian(GlobalDataKey)
Adds the values that have been cached by calling cacheJacobian() and or cacheJacobianNeighbor() to th...
Definition: Assembly.C:3820
DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
MooseArray< Point > _current_q_points_face
The current quadrature points on a face.
Definition: Assembly.h:2496
void cacheResidualBlock(std::vector< Real > &cached_residual_values, std::vector< dof_id_type > &cached_residual_rows, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Push a local residual block with proper scaling into cache.
Definition: Assembly.C:3270
void prepareBlockNonlocal(unsigned int ivar, unsigned jvar, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices)
Definition: Assembly.C:2924
virtual NumericVector< Number > & getVector(const std::string &name)
Get a raw NumericVector by name.
Definition: SystemBase.C:916
void bumpAllQRuleOrder(Order order, SubdomainID block)
Increases the element/volume and face/area quadrature orders for the specified mesh block if and only...
Definition: Assembly.C:600
const VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField< Real > &) const
Definition: Assembly.h:1338
const Elem * _current_side_elem
The current "element" making up the side we are currently on.
Definition: Assembly.h:2576
bool _block_diagonal_matrix
Will be true if our preconditioning matrix is a block-diagonal matrix. Which means that we can take s...
Definition: Assembly.h:2785
bool _need_JxW_neighbor
Flag to indicate that JxW_neighbor is needed.
Definition: Assembly.h:2537
virtual const VectorTag & getVectorTag(const TagID tag_id) const
Get a VectorTag from a TagID.
Definition: SubProblem.C:161
void copyNeighborShapes(MooseVariableField< T > &v)
Definition: Assembly.C:3074
void cacheJacobianCoupledVarPair(const MooseVariableBase &ivar, const MooseVariableBase &jvar)
Caches element matrix for ivar rows and jvar columns.
Definition: Assembly.C:4079
QMONOMIAL
dof_id_type node_id(const unsigned int i) const
MooseArray< Point > _current_q_points_face_neighbor
The current quadrature points on the neighbor face.
Definition: Assembly.h:2535
virtual void attach_quadrature_rule(QBase *q)=0
libMesh::ElemSideBuilder _current_side_elem_builder
In place side element builder for _current_side_elem.
Definition: Assembly.h:2847
VectorVariablePhiCurl _vector_curl_phi_face
Definition: Assembly.h:2699
Key structure for APIs adding/caching local element residuals/Jacobians.
Definition: Assembly.h:833
std::set< FEType > _need_face_div
Definition: Assembly.h:2839
SubdomainID _current_subdomain_id
The current subdomain ID.
Definition: Assembly.h:2568
dof_id_type get_extra_integer(const unsigned int index) const
virtual numeric_index_type n() const=0
std::map< FEType, FEVectorBase * > _current_vector_fe_neighbor
The "neighbor" vector fe object that matches the current elem.
Definition: Assembly.h:2365
void clearCachedResiduals(GlobalDataKey)
Clears all of the residuals in _cached_residual_rows and _cached_residual_values. ...
Definition: Assembly.C:3505
virtual const FieldVariablePhiGradient & gradPhiFace() const =0
Return the gradients of the variable&#39;s shape functions on an element face.
unsigned int THREAD_ID
Definition: MooseTypes.h:209
const Node * _current_neighbor_node
The current neighboring node we are working with.
Definition: Assembly.h:2594
std::vector< VectorValue< ADReal > > _ad_d2xyzdeta2_map
Definition: Assembly.h:2802
void clearCachedJacobian()
Clear any currently cached jacobians.
Definition: Assembly.C:4527
uint8_t dof_id_type
virtual_for_inffe const std::vector< std::vector< OutputShape > > & get_curl_phi() const
const std::vector< std::vector< OutputShape > > & get_phi() const
const VariablePhiGradient & gradPhiFace() const
Definition: Assembly.h:1307
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Knl
dprimary/dlower (or dneighbor/dlower)
Definition: Assembly.h:2667
const std::vector< std::vector< OutputTensor > > & get_dual_d2phi() const
MooseVariableField< T > & getActualFieldVariable(THREAD_ID tid, const std::string &var_name)
Returns a field variable pointer - this includes finite volume variables.
Definition: SystemBase.C:117
void reinitFEFace(const Elem *elem, unsigned int side)
Just an internal helper function to reinit the face FE objects.
Definition: Assembly.C:1269
MooseArray< ADReal > _ad_JxW
Definition: Assembly.h:2804
MooseArray< ADReal > _ad_JxW_face
Definition: Assembly.h:2816
void computeFaceMap(const Elem &elem, const unsigned int side, const std::vector< Real > &qw)
Definition: Assembly.C:1349
void reinitDual(const Elem *elem, const std::vector< Point > &pts, const std::vector< Real > &JxW)
Reintialize dual basis coefficients based on a customized quadrature rule.
Definition: Assembly.C:2275
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Key structure for APIs manipulating global vectors/matrices.
Definition: Assembly.h:815