www.mooseframework.org
Assembly.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "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 
23 // libMesh
24 #include "libmesh/coupling_matrix.h"
25 #include "libmesh/dof_map.h"
26 #include "libmesh/elem.h"
27 #include "libmesh/equation_systems.h"
28 #include "libmesh/fe_interface.h"
29 #include "libmesh/node.h"
30 #include "libmesh/quadrature_gauss.h"
31 #include "libmesh/sparse_matrix.h"
32 #include "libmesh/tensor_value.h"
33 #include "libmesh/vector_value.h"
34 #include "libmesh/fe.h"
35 
37  : _sys(sys),
38  _subproblem(_sys.subproblem()),
39  _displaced(dynamic_cast<DisplacedSystem *>(&sys) ? true : false),
40  _nonlocal_cm(_subproblem.nonlocalCouplingMatrix()),
41  _computing_jacobian(_subproblem.currentlyComputingJacobian()),
42  _dof_map(_sys.dofMap()),
43  _tid(tid),
44  _mesh(sys.mesh()),
45  _mesh_dimension(_mesh.dimension()),
46  _current_qrule(nullptr),
47  _current_qrule_volume(nullptr),
48  _current_qrule_arbitrary(nullptr),
49  _coord_type(Moose::COORD_XYZ),
50  _current_qrule_face(nullptr),
51  _current_qface_arbitrary(nullptr),
52  _current_qrule_neighbor(nullptr),
53  _qrule_msm(nullptr),
54 
55  _current_elem(nullptr),
56  _current_elem_volume(0),
57  _current_side(0),
58  _current_side_elem(nullptr),
59  _current_side_volume(0),
60  _current_neighbor_elem(nullptr),
61  _current_neighbor_side(0),
62  _current_neighbor_side_elem(nullptr),
63  _need_neighbor_elem_volume(false),
64  _current_neighbor_volume(0),
65  _current_node(nullptr),
66  _current_neighbor_node(nullptr),
67  _current_elem_volume_computed(false),
68  _current_side_volume_computed(false),
69 
70  _current_lower_d_elem(nullptr),
71 
72  _cached_residual_values(2), // The 2 is for TIME and NONTIME
73  _cached_residual_rows(2), // The 2 is for TIME and NONTIME
74 
75  _max_cached_residuals(0),
76  _max_cached_jacobians(0),
77  _block_diagonal_matrix(false),
78  _calculate_xyz(false),
79  _calculate_face_xyz(false),
80  _calculate_curvatures(false)
81 {
82  Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
83  // Build fe's for the helpers
84  buildFE(FEType(helper_order, LAGRANGE));
85  buildFaceFE(FEType(helper_order, LAGRANGE));
86  buildNeighborFE(FEType(helper_order, LAGRANGE));
87  buildFaceNeighborFE(FEType(helper_order, LAGRANGE));
88 
89  // Build an FE helper object for this type for each dimension up to the dimension of the current
90  // mesh
91  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
92  {
93  _holder_fe_helper[dim] = &_fe[dim][FEType(helper_order, LAGRANGE)];
94  (*_holder_fe_helper[dim])->get_phi();
95  (*_holder_fe_helper[dim])->get_dphi();
96  (*_holder_fe_helper[dim])->get_xyz();
97  (*_holder_fe_helper[dim])->get_JxW();
98 
99  _holder_fe_face_helper[dim] = &_fe_face[dim][FEType(helper_order, LAGRANGE)];
100  (*_holder_fe_face_helper[dim])->get_phi();
101  (*_holder_fe_face_helper[dim])->get_dphi();
102  (*_holder_fe_face_helper[dim])->get_xyz();
103  (*_holder_fe_face_helper[dim])->get_JxW();
104  (*_holder_fe_face_helper[dim])->get_normals();
105 
106  _holder_fe_face_neighbor_helper[dim] = &_fe_face_neighbor[dim][FEType(helper_order, LAGRANGE)];
107  (*_holder_fe_face_neighbor_helper[dim])->get_xyz();
108  (*_holder_fe_face_neighbor_helper[dim])->get_JxW();
109  (*_holder_fe_face_neighbor_helper[dim])->get_normals();
110 
111  _holder_fe_neighbor_helper[dim] = &_fe_neighbor[dim][FEType(helper_order, LAGRANGE)];
112  (*_holder_fe_neighbor_helper[dim])->get_xyz();
113  (*_holder_fe_neighbor_helper[dim])->get_JxW();
114  }
115 
116  _fe_msm = FEGenericBase<Real>::build(_mesh_dimension - 1, FEType(helper_order, LAGRANGE));
117  _JxW_msm = &_fe_msm->get_JxW();
118 }
119 
121 {
122  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
123  for (auto & it : _fe[dim])
124  delete it.second;
125 
126  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
127  for (auto & it : _fe_face[dim])
128  delete it.second;
129 
130  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
131  for (auto & it : _fe_neighbor[dim])
132  delete it.second;
133 
134  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
135  for (auto & it : _fe_face_neighbor[dim])
136  delete it.second;
137 
138  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
139  for (auto & it : _fe_lower[dim])
140  delete it.second;
141 
142  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
143  for (auto & it : _vector_fe[dim])
144  delete it.second;
145 
146  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
147  for (auto & it : _vector_fe_face[dim])
148  delete it.second;
149 
150  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
151  for (auto & it : _vector_fe_neighbor[dim])
152  delete it.second;
153 
154  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
155  for (auto & it : _vector_fe_face_neighbor[dim])
156  delete it.second;
157 
158  for (auto & it : _holder_qrule_volume)
159  delete it.second;
160 
161  for (auto & it : _holder_qrule_arbitrary)
162  delete it.second;
163 
164  for (auto & it : _holder_qrule_arbitrary_face)
165  delete it.second;
166 
167  for (auto & it : _holder_qrule_face)
168  delete it.second;
169 
170  for (auto & it : _holder_qrule_neighbor)
171  delete it.second;
172 
173  for (auto & it : _fe_shape_data)
174  delete it.second;
175 
176  for (auto & it : _fe_shape_data_face)
177  delete it.second;
178 
179  for (auto & it : _fe_shape_data_neighbor)
180  delete it.second;
181 
182  for (auto & it : _fe_shape_data_face_neighbor)
183  delete it.second;
184 
185  for (auto & it : _vector_fe_shape_data)
186  delete it.second;
187 
188  for (auto & it : _vector_fe_shape_data_face)
189  delete it.second;
190 
191  for (auto & it : _vector_fe_shape_data_neighbor)
192  delete it.second;
193 
194  for (auto & it : _vector_fe_shape_data_face_neighbor)
195  delete it.second;
196 
197  for (auto & it : _fe_shape_data_lower)
198  delete it.second;
199 
200  for (auto & it : _ad_grad_phi_data)
201  it.second.release();
202 
203  for (auto & it : _ad_vector_grad_phi_data)
204  it.second.release();
205 
206  for (auto & it : _ad_grad_phi_data_face)
207  it.second.release();
208 
209  for (auto & it : _ad_vector_grad_phi_data_face)
210  it.second.release();
211 
212  delete _current_side_elem;
214 
216 
217  _coord.release();
219 
220  _ad_JxW.release();
227  _ad_coord.release();
228 
229  if (_qrule_msm)
230  delete _qrule_msm;
231 }
232 
233 void
234 Assembly::buildFE(FEType type) const
235 {
236  if (!_fe_shape_data[type])
238 
239  // Build an FE object for this type for each dimension up to the dimension of the current mesh
240  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
241  {
242  if (!_fe[dim][type])
243  _const_fe[dim][type] = _fe[dim][type] = FEGenericBase<Real>::build(dim, type).release();
244 
245  _fe[dim][type]->get_phi();
246  _fe[dim][type]->get_dphi();
247  // Pre-request xyz. We have always computed xyz, but due to
248  // recent optimizations in libmesh, we now need to explicity
249  // request it, since apps (Yak) may rely on it being computed.
250  _fe[dim][type]->get_xyz();
252  _fe[dim][type]->get_d2phi();
253  }
254 }
255 
256 void
258 {
261 
262  // Build an FE object for this type for each dimension up to the dimension of the current mesh
263  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
264  {
265  if (!_fe_face[dim][type])
266  _const_fe_face[dim][type] = _fe_face[dim][type] =
267  FEGenericBase<Real>::build(dim, type).release();
268 
269  _fe_face[dim][type]->get_phi();
270  _fe_face[dim][type]->get_dphi();
272  _fe_face[dim][type]->get_d2phi();
273  }
274 }
275 
276 void
278 {
281 
282  // Build an FE object for this type for each dimension up to the dimension of the current mesh
283  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
284  {
285  if (!_fe_neighbor[dim][type])
286  _const_fe_neighbor[dim][type] = _fe_neighbor[dim][type] =
287  FEGenericBase<Real>::build(dim, type).release();
288 
289  _fe_neighbor[dim][type]->get_phi();
290  _fe_neighbor[dim][type]->get_dphi();
292  _fe_neighbor[dim][type]->get_d2phi();
293  }
294 }
295 
296 void
298 {
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_neighbor[dim][type])
307  FEGenericBase<Real>::build(dim, type).release();
308 
309  _fe_face_neighbor[dim][type]->get_phi();
310  _fe_face_neighbor[dim][type]->get_dphi();
312  _fe_face_neighbor[dim][type]->get_d2phi();
313  }
314 }
315 
316 void
318 {
321 
322  // Build an FE object for this type for each dimension up to the dimension of
323  // the current mesh minus one (because this is for lower-dimensional
324  // elements!)
325  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
326  {
327  if (!_fe_lower[dim][type])
328  _const_fe_lower[dim][type] = _fe_lower[dim][type] =
329  FEGenericBase<Real>::build(dim, type).release();
330 
331  _fe_lower[dim][type]->get_phi();
332  _fe_lower[dim][type]->get_dphi();
334  _fe_lower[dim][type]->get_d2phi();
335  }
336 }
337 
338 void
340 {
343 
344  // Build an FE object for this type for each dimension up to the dimension of
345  // the current mesh minus one (because this is for lower-dimensional
346  // elements!)
347  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
348  {
349  if (!_vector_fe_lower[dim][type])
351  FEVectorBase::build(dim, type).release();
352 
353  _vector_fe_lower[dim][type]->get_phi();
354  _vector_fe_lower[dim][type]->get_dphi();
356  _vector_fe_lower[dim][type]->get_d2phi();
357  }
358 }
359 
360 void
362 {
365 
366  unsigned int min_dim;
367  if (type.family == LAGRANGE_VEC)
368  min_dim = 0;
369  else
370  min_dim = 2;
371 
372  // Build an FE object for this type for each dimension up to the dimension of the current mesh
373  // Note that NEDELEC_ONE elements can only be built for dimension > 2. The for loop logic should
374  // be modified for LAGRANGE_VEC
375  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
376  {
377  if (!_vector_fe[dim][type])
378  _const_vector_fe[dim][type] = _vector_fe[dim][type] =
379  FEGenericBase<VectorValue<Real>>::build(dim, type).release();
380 
381  _vector_fe[dim][type]->get_phi();
382  _vector_fe[dim][type]->get_dphi();
383  if (type.family == NEDELEC_ONE)
384  _vector_fe[dim][type]->get_curl_phi();
385  // Pre-request xyz. We have always computed xyz, but due to
386  // recent optimizations in libmesh, we now need to explicity
387  // request it, since apps (Yak) may rely on it being computed.
388  _vector_fe[dim][type]->get_xyz();
389  }
390 }
391 
392 void
394 {
397 
398  unsigned int min_dim;
399  if (type.family == LAGRANGE_VEC)
400  min_dim = 0;
401  else
402  min_dim = 2;
403 
404  // Build an VectorFE object for this type for each dimension up to the dimension of the current
405  // mesh
406  // Note that NEDELEC_ONE elements can only be built for dimension > 2. The for loop logic should
407  // be modified for LAGRANGE_VEC
408  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
409  {
410  if (!_vector_fe_face[dim][type])
412  FEGenericBase<VectorValue<Real>>::build(dim, type).release();
413 
414  _vector_fe_face[dim][type]->get_phi();
415  _vector_fe_face[dim][type]->get_dphi();
416  if (type.family == NEDELEC_ONE)
417  _vector_fe_face[dim][type]->get_curl_phi();
418  }
419 }
420 
421 void
423 {
426 
427  unsigned int min_dim;
428  if (type.family == LAGRANGE_VEC)
429  min_dim = 0;
430  else
431  min_dim = 2;
432 
433  // Build an VectorFE object for this type for each dimension up to the dimension of the current
434  // mesh
435  // Note that NEDELEC_ONE elements can only be built for dimension > 2. The for loop logic should
436  // be modified for LAGRANGE_VEC
437  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
438  {
439  if (!_vector_fe_neighbor[dim][type])
441  FEGenericBase<VectorValue<Real>>::build(dim, type).release();
442 
443  _vector_fe_neighbor[dim][type]->get_phi();
444  _vector_fe_neighbor[dim][type]->get_dphi();
445  if (type.family == NEDELEC_ONE)
446  _vector_fe_neighbor[dim][type]->get_curl_phi();
447  }
448 }
449 
450 void
452 {
455 
456  unsigned int min_dim;
457  if (type.family == LAGRANGE_VEC)
458  min_dim = 0;
459  else
460  min_dim = 2;
461 
462  // Build an VectorFE object for this type for each dimension up to the dimension of the current
463  // mesh
464  // Note that NEDELEC_ONE elements can only be built for dimension > 2. The for loop logic should
465  // be modified for LAGRANGE_VEC
466  for (unsigned int dim = min_dim; dim <= _mesh_dimension; dim++)
467  {
468  if (!_vector_fe_face_neighbor[dim][type])
470  FEGenericBase<VectorValue<Real>>::build(dim, type).release();
471 
472  _vector_fe_face_neighbor[dim][type]->get_phi();
473  _vector_fe_face_neighbor[dim][type]->get_dphi();
474  if (type.family == NEDELEC_ONE)
475  _vector_fe_face_neighbor[dim][type]->get_curl_phi();
476  }
477 }
478 
479 const Real &
481 {
484 }
485 
486 void
487 Assembly::createQRules(QuadratureType type, Order order, Order volume_order, Order face_order)
488 {
489  _holder_qrule_volume.clear();
490  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
491  _holder_qrule_volume[dim] = QBase::build(type, dim, volume_order).release();
492 
493  _holder_qrule_face.clear();
494  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
495  _holder_qrule_face[dim] = QBase::build(type, dim - 1, face_order).release();
496 
497  _holder_qrule_neighbor.clear();
498  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
499  _holder_qrule_neighbor[dim] = new ArbitraryQuadrature(dim - 1, face_order);
500 
501  _holder_qrule_arbitrary.clear();
502  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
503  _holder_qrule_arbitrary[dim] = new ArbitraryQuadrature(dim, order);
504 
506  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
507  _holder_qrule_arbitrary_face[dim] = new ArbitraryQuadrature(dim - 1, face_order);
508 
509  _const_qrule_msm = _qrule_msm = QBase::build(type, _mesh_dimension - 1, face_order).release();
510  _fe_msm->attach_quadrature_rule(_qrule_msm);
511 }
512 
513 void
514 Assembly::setVolumeQRule(QBase * qrule, unsigned int dim)
515 {
517 
518  if (qrule) // Don't set a NULL qrule
519  {
520  for (auto & it : _fe[dim])
521  it.second->attach_quadrature_rule(qrule);
522  for (auto & it : _vector_fe[dim])
523  it.second->attach_quadrature_rule(qrule);
524  }
525 }
526 
527 void
528 Assembly::setFaceQRule(QBase * qrule, unsigned int dim)
529 {
531 
532  for (auto & it : _fe_face[dim])
533  it.second->attach_quadrature_rule(qrule);
534  for (auto & it : _vector_fe_face[dim])
535  it.second->attach_quadrature_rule(qrule);
536 }
537 
538 void
539 Assembly::setNeighborQRule(QBase * qrule, unsigned int dim)
540 {
542 
543  for (auto & it : _fe_face_neighbor[dim])
544  it.second->attach_quadrature_rule(qrule);
545  for (auto & it : _vector_fe_face_neighbor[dim])
546  it.second->attach_quadrature_rule(qrule);
547 }
548 
549 void
550 Assembly::reinitFE(const Elem * elem)
551 {
552  unsigned int dim = elem->dim();
553 
554  for (const auto & it : _fe[dim])
555  {
556  FEBase * fe = it.second;
557  const FEType & fe_type = it.first;
558 
559  _current_fe[fe_type] = fe;
560 
561  FEShapeData * fesd = _fe_shape_data[fe_type];
562 
563  fe->reinit(elem);
564 
565  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe->get_phi()));
566  fesd->_grad_phi.shallowCopy(
567  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe->get_dphi()));
568  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
569  fesd->_second_phi.shallowCopy(
570  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe->get_d2phi()));
571  }
572  for (const auto & it : _vector_fe[dim])
573  {
574  FEVectorBase * fe = it.second;
575  const FEType & fe_type = it.first;
576 
577  _current_vector_fe[fe_type] = fe;
578 
579  VectorFEShapeData * fesd = _vector_fe_shape_data[fe_type];
580 
581  fe->reinit(elem);
582 
583  fesd->_phi.shallowCopy(
584  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe->get_phi()));
585  fesd->_grad_phi.shallowCopy(
586  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe->get_dphi()));
587  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
588  fesd->_second_phi.shallowCopy(
589  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe->get_d2phi()));
590  if (_need_curl.find(fe_type) != _need_curl.end())
591  fesd->_curl_phi.shallowCopy(
592  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe->get_curl_phi()));
593  }
594 
595  // During that last loop the helper objects will have been reinitialized as well
596  // We need to dig out the q_points and JxW from it.
598  const_cast<std::vector<Point> &>((*_holder_fe_helper[dim])->get_xyz()));
599  _current_JxW.shallowCopy(const_cast<std::vector<Real> &>((*_holder_fe_helper[dim])->get_JxW()));
600 
602  {
603  auto n_qp = _current_qrule->n_points();
604  resizeMappingObjects(n_qp, dim);
605  _ad_JxW.resize(n_qp);
606  if (_calculate_xyz)
607  _ad_q_points.resize(n_qp);
608  if (_displaced)
609  {
610  const auto & qw = _current_qrule->get_weights();
611  if (elem->has_affine_map())
612  computeAffineMapAD(elem, qw, n_qp, *_holder_fe_helper[dim]);
613  else
614  for (unsigned int qp = 0; qp != n_qp; qp++)
616  }
617  else
618  for (unsigned qp = 0; qp < n_qp; ++qp)
619  {
620  _ad_JxW[qp] = _current_JxW[qp];
621  if (_calculate_xyz)
623  }
624 
625  for (const auto & it : _fe[dim])
626  {
627  FEBase * fe = it.second;
628  auto fe_type = it.first;
629  auto num_shapes = fe->n_shape_functions();
630  auto & grad_phi = _ad_grad_phi_data[fe_type];
631 
632  grad_phi.resize(num_shapes);
633  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
634  grad_phi[i].resize(n_qp);
635 
636  if (_displaced)
637  computeGradPhiAD(elem, n_qp, grad_phi, fe);
638  else
639  {
640  const auto & regular_grad_phi = _fe_shape_data[fe_type]->_grad_phi;
641  for (unsigned qp = 0; qp < n_qp; ++qp)
642  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
643  grad_phi[i][qp] = regular_grad_phi[i][qp];
644  }
645  }
646  for (const auto & it : _vector_fe[dim])
647  {
648  FEVectorBase * fe = it.second;
649  auto fe_type = it.first;
650  auto num_shapes = fe->n_shape_functions();
651  auto & grad_phi = _ad_vector_grad_phi_data[fe_type];
652 
653  grad_phi.resize(num_shapes);
654  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
655  grad_phi[i].resize(n_qp);
656 
657  if (_displaced)
658  computeGradPhiAD(elem, n_qp, grad_phi, fe);
659  else
660  {
661  const auto & regular_grad_phi = _vector_fe_shape_data[fe_type]->_grad_phi;
662  for (unsigned qp = 0; qp < n_qp; ++qp)
663  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
664  grad_phi[i][qp] = regular_grad_phi[i][qp];
665  }
666  }
667  }
668 
669  if (_xfem != nullptr)
671 }
672 
673 template <typename OutputType>
674 void
676  const Elem * elem,
677  unsigned int n_qp,
679  FEGenericBase<OutputType> * fe)
680 {
681  auto dim = elem->dim();
682  const auto & dphidxi = fe->get_dphidxi();
683  const auto & dphideta = fe->get_dphideta();
684  const auto & dphidzeta = fe->get_dphidzeta();
685  auto num_shapes = grad_phi.size();
686 
687  switch (dim)
688  {
689  case 0:
690  {
691  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
692  for (unsigned qp = 0; qp < n_qp; ++qp)
693  grad_phi[i][qp] = 0;
694  break;
695  }
696 
697  case 1:
698  {
699  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
700  for (unsigned qp = 0; qp < n_qp; ++qp)
701  {
702  grad_phi[i][qp](0) = dphidxi[i][qp] * _ad_dxidx_map[qp];
703  grad_phi[i][qp](1) = dphidxi[i][qp] * _ad_dxidy_map[qp];
704  grad_phi[i][qp](2) = dphidxi[i][qp] * _ad_dxidz_map[qp];
705  }
706  break;
707  }
708 
709  case 2:
710  {
711  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
712  for (unsigned qp = 0; qp < n_qp; ++qp)
713  {
714  grad_phi[i][qp](0) =
715  dphidxi[i][qp] * _ad_dxidx_map[qp] + dphideta[i][qp] * _ad_detadx_map[qp];
716  grad_phi[i][qp](1) =
717  dphidxi[i][qp] * _ad_dxidy_map[qp] + dphideta[i][qp] * _ad_detady_map[qp];
718  grad_phi[i][qp](2) =
719  dphidxi[i][qp] * _ad_dxidz_map[qp] + dphideta[i][qp] * _ad_detadz_map[qp];
720  }
721  break;
722  }
723 
724  case 3:
725  {
726  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
727  for (unsigned qp = 0; qp < n_qp; ++qp)
728  {
729  grad_phi[i][qp](0) = dphidxi[i][qp] * _ad_dxidx_map[qp] +
730  dphideta[i][qp] * _ad_detadx_map[qp] +
731  dphidzeta[i][qp] * _ad_dzetadx_map[qp];
732  grad_phi[i][qp](1) = dphidxi[i][qp] * _ad_dxidy_map[qp] +
733  dphideta[i][qp] * _ad_detady_map[qp] +
734  dphidzeta[i][qp] * _ad_dzetady_map[qp];
735  grad_phi[i][qp](2) = dphidxi[i][qp] * _ad_dxidz_map[qp] +
736  dphideta[i][qp] * _ad_detadz_map[qp] +
737  dphidzeta[i][qp] * _ad_dzetadz_map[qp];
738  }
739  break;
740  }
741  }
742 }
743 
744 void
745 Assembly::resizeMappingObjects(unsigned int n_qp, unsigned int dim)
746 {
747  _ad_dxyzdxi_map.resize(n_qp);
748  _ad_dxidx_map.resize(n_qp);
749  _ad_dxidy_map.resize(n_qp); // 1D element may live in 2D ...
750  _ad_dxidz_map.resize(n_qp); // ... or 3D
751 
752  if (dim > 1)
753  {
754  _ad_dxyzdeta_map.resize(n_qp);
755  _ad_detadx_map.resize(n_qp);
756  _ad_detady_map.resize(n_qp);
757  _ad_detadz_map.resize(n_qp);
758 
759  if (dim > 2)
760  {
761  _ad_dxyzdzeta_map.resize(n_qp);
762  _ad_dzetadx_map.resize(n_qp);
763  _ad_dzetady_map.resize(n_qp);
764  _ad_dzetadz_map.resize(n_qp);
765  }
766  }
767 
768  _ad_jac.resize(n_qp);
769 }
770 
771 void
773  const std::vector<Real> & qw,
774  unsigned int n_qp,
775  FEBase * fe)
776 {
777  computeSinglePointMapAD(elem, qw, 0, fe);
778 
779  for (unsigned int p = 1; p < n_qp; p++)
780  {
781  // Compute xyz at all other quadrature points. Note that this map method call is only really
782  // referring to the volumetric element...face quadrature points are calculated in computeFaceMap
783  if (_calculate_xyz)
784  {
785  auto num_shapes = fe->n_shape_functions();
786  const auto & elem_nodes = elem->get_nodes();
787  const auto & phi_map = fe->get_fe_map().get_phi_map();
788  _ad_q_points[p].zero();
789  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
790  {
791  mooseAssert(elem_nodes[i], "The node is null!");
792  VectorValue<DualReal> elem_point = *elem_nodes[i];
793  unsigned dimension = 0;
794  for (const auto & disp_num : _displacements)
795  elem_point(dimension++).derivatives()[disp_num * _sys.getMaxVarNDofsPerElem() + i] = 1.;
796 
797  _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
798  }
799  }
800 
801  // Now copy over other map data for each extra quadrature point
802 
807 
808  if (elem->dim() > 1)
809  {
814 
815  if (elem->dim() > 2)
816  {
821  }
822  }
823  _ad_jac[p] = _ad_jac[0];
824  _ad_JxW[p] = _ad_JxW[0] / qw[0] * qw[p];
825  }
826 }
827 
828 void
830  const std::vector<Real> & qw,
831  unsigned p,
832  FEBase * fe)
833 {
834  auto dim = elem->dim();
835  const auto & elem_nodes = elem->get_nodes();
836  auto num_shapes = fe->n_shape_functions();
837  const auto & phi_map = fe->get_fe_map().get_phi_map();
838  const auto & dphidxi_map = fe->get_fe_map().get_dphidxi_map();
839  const auto & dphideta_map = fe->get_fe_map().get_dphideta_map();
840  const auto & dphidzeta_map = fe->get_fe_map().get_dphidzeta_map();
841 
842  switch (dim)
843  {
844  case 0:
845  {
846  _ad_jac[p] = 1.0;
847  _ad_JxW[p] = qw[p];
848  if (_calculate_xyz)
849  _ad_q_points[p] = *elem_nodes[0];
850  break;
851  }
852 
853  case 1:
854  {
855  if (_calculate_xyz)
856  _ad_q_points[p].zero();
857 
858  _ad_dxyzdxi_map[p].zero();
859 
860  for (std::size_t i = 0; i < num_shapes; i++)
861  {
862  libmesh_assert(elem_nodes[i]);
863  libMesh::VectorValue<DualReal> elem_point = *elem_nodes[i];
864  unsigned dimension = 0;
865  for (const auto & disp_num : _displacements)
866  elem_point(dimension++).derivatives()[disp_num * _sys.getMaxVarNDofsPerElem() + i] = 1.;
867 
868  _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
869 
870  if (_calculate_xyz)
871  _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
872  }
873 
874  _ad_jac[p] = _ad_dxyzdxi_map[p].norm();
875 
876  if (_ad_jac[p].value() <= -TOLERANCE * TOLERANCE)
877  {
878  static bool failing = false;
879  if (!failing)
880  {
881  failing = true;
882  elem->print_info(libMesh::err);
883  libmesh_error_msg("ERROR: negative Jacobian " << _ad_jac[p].value() << " at point index "
884  << p << " in element " << elem->id());
885  }
886  else
887  return;
888  }
889 
890  const auto jacm2 = 1. / _ad_jac[p] / _ad_jac[p];
891  _ad_dxidx_map[p] = jacm2 * _ad_dxyzdxi_map[p](0);
892  _ad_dxidy_map[p] = jacm2 * _ad_dxyzdxi_map[p](1);
893  _ad_dxidz_map[p] = jacm2 * _ad_dxyzdxi_map[p](2);
894 
895  _ad_JxW[p] = _ad_jac[p] * qw[p];
896 
897  break;
898  }
899 
900  case 2:
901  {
902  if (_calculate_xyz)
903  _ad_q_points[p].zero();
904  _ad_dxyzdxi_map[p].zero();
905  _ad_dxyzdeta_map[p].zero();
906 
907  for (std::size_t i = 0; i < num_shapes; i++)
908  {
909  libmesh_assert(elem_nodes[i]);
910  libMesh::VectorValue<DualReal> elem_point = *elem_nodes[i];
911  unsigned dimension = 0;
912  for (const auto & disp_num : _displacements)
913  elem_point(dimension++).derivatives()[disp_num * _sys.getMaxVarNDofsPerElem() + i] = 1.;
914 
915  _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
916  _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
917 
918  if (_calculate_xyz)
919  _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
920  }
921 
922  const auto &dx_dxi = _ad_dxyzdxi_map[p](0), dx_deta = _ad_dxyzdeta_map[p](0),
923  dy_dxi = _ad_dxyzdxi_map[p](1), dy_deta = _ad_dxyzdeta_map[p](1),
924  dz_dxi = _ad_dxyzdxi_map[p](2), dz_deta = _ad_dxyzdeta_map[p](2);
925 
926  const auto g11 = (dx_dxi * dx_dxi + dy_dxi * dy_dxi + dz_dxi * dz_dxi);
927 
928  const auto g12 = (dx_dxi * dx_deta + dy_dxi * dy_deta + dz_dxi * dz_deta);
929 
930  const auto g21 = g12;
931 
932  const auto g22 = (dx_deta * dx_deta + dy_deta * dy_deta + dz_deta * dz_deta);
933 
934  auto det = (g11 * g22 - g12 * g21);
935 
936  if (det.value() <= -TOLERANCE * TOLERANCE)
937  {
938  static bool failing = false;
939  if (!failing)
940  {
941  failing = true;
942  elem->print_info(libMesh::err);
943  libmesh_error_msg("ERROR: negative Jacobian " << det << " at point index " << p
944  << " in element " << elem->id());
945  }
946  else
947  return;
948  }
949  else if (det.value() <= 0.)
950  det.value() = TOLERANCE * TOLERANCE;
951 
952  const auto inv_det = 1. / det;
953  _ad_jac[p] = std::sqrt(det);
954 
955  _ad_JxW[p] = _ad_jac[p] * qw[p];
956 
957  const auto g11inv = g22 * inv_det;
958  const auto g12inv = -g12 * inv_det;
959  const auto g21inv = -g21 * inv_det;
960  const auto g22inv = g11 * inv_det;
961 
962  _ad_dxidx_map[p] = g11inv * dx_dxi + g12inv * dx_deta;
963  _ad_dxidy_map[p] = g11inv * dy_dxi + g12inv * dy_deta;
964  _ad_dxidz_map[p] = g11inv * dz_dxi + g12inv * dz_deta;
965 
966  _ad_detadx_map[p] = g21inv * dx_dxi + g22inv * dx_deta;
967  _ad_detady_map[p] = g21inv * dy_dxi + g22inv * dy_deta;
968  _ad_detadz_map[p] = g21inv * dz_dxi + g22inv * dz_deta;
969 
970  break;
971  }
972 
973  case 3:
974  {
975  if (_calculate_xyz)
976  _ad_q_points[p].zero();
977  _ad_dxyzdxi_map[p].zero();
978  _ad_dxyzdeta_map[p].zero();
979  _ad_dxyzdzeta_map[p].zero();
980 
981  for (std::size_t i = 0; i < num_shapes; i++)
982  {
983  libmesh_assert(elem_nodes[i]);
984  libMesh::VectorValue<DualReal> elem_point = *elem_nodes[i];
985  unsigned dimension = 0;
986  for (const auto & disp_num : _displacements)
987  elem_point(dimension++).derivatives()[disp_num * _sys.getMaxVarNDofsPerElem() + i] = 1.;
988 
989  _ad_dxyzdxi_map[p].add_scaled(elem_point, dphidxi_map[i][p]);
990  _ad_dxyzdeta_map[p].add_scaled(elem_point, dphideta_map[i][p]);
991  _ad_dxyzdzeta_map[p].add_scaled(elem_point, dphidzeta_map[i][p]);
992 
993  if (_calculate_xyz)
994  _ad_q_points[p].add_scaled(elem_point, phi_map[i][p]);
995  }
996 
997  const auto dx_dxi = _ad_dxyzdxi_map[p](0), dy_dxi = _ad_dxyzdxi_map[p](1),
998  dz_dxi = _ad_dxyzdxi_map[p](2), dx_deta = _ad_dxyzdeta_map[p](0),
999  dy_deta = _ad_dxyzdeta_map[p](1), dz_deta = _ad_dxyzdeta_map[p](2),
1000  dx_dzeta = _ad_dxyzdzeta_map[p](0), dy_dzeta = _ad_dxyzdzeta_map[p](1),
1001  dz_dzeta = _ad_dxyzdzeta_map[p](2);
1002 
1003  _ad_jac[p] = (dx_dxi * (dy_deta * dz_dzeta - dz_deta * dy_dzeta) +
1004  dy_dxi * (dz_deta * dx_dzeta - dx_deta * dz_dzeta) +
1005  dz_dxi * (dx_deta * dy_dzeta - dy_deta * dx_dzeta));
1006 
1007  if (_ad_jac[p].value() <= -TOLERANCE * TOLERANCE)
1008  {
1009  static bool failing = false;
1010  if (!failing)
1011  {
1012  failing = true;
1013  elem->print_info(libMesh::err);
1014  libmesh_error_msg("ERROR: negative Jacobian " << _ad_jac[p].value() << " at point index "
1015  << p << " in element " << elem->id());
1016  }
1017  else
1018  return;
1019  }
1020 
1021  _ad_JxW[p] = _ad_jac[p] * qw[p];
1022 
1023  const auto inv_jac = 1. / _ad_jac[p];
1024 
1025  _ad_dxidx_map[p] = (dy_deta * dz_dzeta - dz_deta * dy_dzeta) * inv_jac;
1026  _ad_dxidy_map[p] = (dz_deta * dx_dzeta - dx_deta * dz_dzeta) * inv_jac;
1027  _ad_dxidz_map[p] = (dx_deta * dy_dzeta - dy_deta * dx_dzeta) * inv_jac;
1028 
1029  _ad_detadx_map[p] = (dz_dxi * dy_dzeta - dy_dxi * dz_dzeta) * inv_jac;
1030  _ad_detady_map[p] = (dx_dxi * dz_dzeta - dz_dxi * dx_dzeta) * inv_jac;
1031  _ad_detadz_map[p] = (dy_dxi * dx_dzeta - dx_dxi * dy_dzeta) * inv_jac;
1032 
1033  _ad_dzetadx_map[p] = (dy_dxi * dz_deta - dz_dxi * dy_deta) * inv_jac;
1034  _ad_dzetady_map[p] = (dz_dxi * dx_deta - dx_dxi * dz_deta) * inv_jac;
1035  _ad_dzetadz_map[p] = (dx_dxi * dy_deta - dy_dxi * dx_deta) * inv_jac;
1036 
1037  break;
1038  }
1039 
1040  default:
1041  libmesh_error_msg("Invalid dim = " << dim);
1042  }
1043 }
1044 
1045 void
1046 Assembly::reinitFEFace(const Elem * elem, unsigned int side)
1047 {
1048  unsigned int dim = elem->dim();
1049 
1050  for (const auto & it : _fe_face[dim])
1051  {
1052  FEBase * fe_face = it.second;
1053  const FEType & fe_type = it.first;
1054  FEShapeData * fesd = _fe_shape_data_face[fe_type];
1055  fe_face->reinit(elem, side);
1056  _current_fe_face[fe_type] = fe_face;
1057 
1058  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face->get_phi()));
1059  fesd->_grad_phi.shallowCopy(
1060  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face->get_dphi()));
1061  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
1062  fesd->_second_phi.shallowCopy(
1063  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face->get_d2phi()));
1064  }
1065  for (const auto & it : _vector_fe_face[dim])
1066  {
1067  FEVectorBase * fe_face = it.second;
1068  const FEType & fe_type = it.first;
1069 
1070  _current_vector_fe_face[fe_type] = fe_face;
1071 
1073 
1074  fe_face->reinit(elem, side);
1075 
1076  fesd->_phi.shallowCopy(
1077  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face->get_phi()));
1078  fesd->_grad_phi.shallowCopy(
1079  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face->get_dphi()));
1080  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
1081  fesd->_second_phi.shallowCopy(
1082  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face->get_d2phi()));
1083  if (_need_curl.find(fe_type) != _need_curl.end())
1084  fesd->_curl_phi.shallowCopy(
1085  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face->get_curl_phi()));
1086  }
1087 
1088  // During that last loop the helper objects will have been reinitialized as well
1089  // We need to dig out the q_points and JxW from it.
1091  const_cast<std::vector<Point> &>((*_holder_fe_face_helper[dim])->get_xyz()));
1093  const_cast<std::vector<Real> &>((*_holder_fe_face_helper[dim])->get_JxW()));
1095  const_cast<std::vector<Point> &>((*_holder_fe_face_helper[dim])->get_normals()));
1098  const_cast<std::vector<Real> &>((*_holder_fe_face_helper[dim])->get_curvatures()));
1099 
1101 
1102  if (_xfem != nullptr)
1104 }
1105 
1106 void
1107 Assembly::computeFaceMap(unsigned dim, const std::vector<Real> & qw, const Elem * side)
1108 {
1109  const auto n_qp = qw.size();
1110  const Elem * elem = side->parent();
1111  auto side_number = elem->which_side_am_i(side);
1112  const auto & dpsidxi_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_dpsidxi();
1113  const auto & dpsideta_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_dpsideta();
1114  const auto & psi_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_psi();
1115  std::vector<std::vector<Real>> const * d2psidxi2_map = nullptr;
1116  std::vector<std::vector<Real>> const * d2psidxideta_map = nullptr;
1117  std::vector<std::vector<Real>> const * d2psideta2_map = nullptr;
1119  {
1120  d2psidxi2_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psidxi2();
1121  d2psidxideta_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psidxideta();
1122  d2psideta2_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psideta2();
1123  }
1124 
1125  switch (dim)
1126  {
1127  case 1:
1128  {
1129  if (!n_qp)
1130  break;
1131 
1132  if (side->node_id(0) == elem->node_id(0))
1133  _ad_normals[0] = Point(-1.);
1134  else
1135  _ad_normals[0] = Point(1.);
1136 
1137  VectorValue<DualReal> side_point;
1138  if (_calculate_face_xyz)
1139  {
1140  side_point = side->point(0);
1141  auto element_node_number = elem->which_node_am_i(side_number, 0);
1142 
1143  unsigned dimension = 0;
1144  for (const auto & disp_num : _displacements)
1145  side_point(dimension++)
1146  .derivatives()[disp_num * _sys.getMaxVarNDofsPerElem() + element_node_number] = 1.;
1147  }
1148 
1149  for (unsigned int p = 0; p < n_qp; p++)
1150  {
1151  if (_calculate_face_xyz)
1152  {
1153  _ad_q_points_face[p].zero();
1154  _ad_q_points_face[p].add_scaled(side_point, psi_map[0][p]);
1155  }
1156 
1157  _ad_normals[p] = _ad_normals[0];
1158  _ad_JxW_face[p] = 1.0 * qw[p];
1159  }
1160 
1161  break;
1162  }
1163 
1164  case 2:
1165  {
1166  _ad_dxyzdxi_map.resize(n_qp);
1168  _ad_d2xyzdxi2_map.resize(n_qp);
1169 
1170  for (unsigned int p = 0; p < n_qp; p++)
1171  {
1172  _ad_dxyzdxi_map[p].zero();
1173  if (_calculate_face_xyz)
1174  _ad_q_points_face[p].zero();
1176  _ad_d2xyzdxi2_map[p].zero();
1177  }
1178 
1179  const auto n_mapping_shape_functions =
1180  FE<2, LAGRANGE>::n_shape_functions(side->type(), side->default_order());
1181 
1182  for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1183  {
1184  VectorValue<DualReal> side_point = side->point(i);
1185  auto element_node_number = elem->which_node_am_i(side_number, i);
1186 
1187  unsigned dimension = 0;
1188  for (const auto & disp_num : _displacements)
1189  side_point(dimension++)
1190  .derivatives()[disp_num * _sys.getMaxVarNDofsPerElem() + element_node_number] = 1.;
1191 
1192  for (unsigned int p = 0; p < n_qp; p++)
1193  {
1194  _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1195  if (_calculate_face_xyz)
1196  _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1198  _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1199  }
1200  }
1201 
1202  for (unsigned int p = 0; p < n_qp; p++)
1203  {
1204  _ad_normals[p] =
1205  (VectorValue<DualReal>(_ad_dxyzdxi_map[p](1), -_ad_dxyzdxi_map[p](0), 0.)).unit();
1206  const auto the_jac = _ad_dxyzdxi_map[p].norm();
1207  _ad_JxW_face[p] = the_jac * qw[p];
1209  {
1210  const auto numerator = _ad_d2xyzdxi2_map[p] * _ad_normals[p];
1211  const auto denominator = _ad_dxyzdxi_map[p].norm_sq();
1212  libmesh_assert_not_equal_to(denominator, 0);
1213  _ad_curvatures[p] = numerator / denominator;
1214  }
1215  }
1216 
1217  break;
1218  }
1219 
1220  case 3:
1221  {
1222  _ad_dxyzdxi_map.resize(n_qp);
1223  _ad_dxyzdeta_map.resize(n_qp);
1225  {
1226  _ad_d2xyzdxi2_map.resize(n_qp);
1227  _ad_d2xyzdxideta_map.resize(n_qp);
1228  _ad_d2xyzdeta2_map.resize(n_qp);
1229  }
1230 
1231  for (unsigned int p = 0; p < n_qp; p++)
1232  {
1233  _ad_dxyzdxi_map[p].zero();
1234  _ad_dxyzdeta_map[p].zero();
1235  if (_calculate_face_xyz)
1236  _ad_q_points_face[p].zero();
1238  {
1239  _ad_d2xyzdxi2_map[p].zero();
1240  _ad_d2xyzdxideta_map[p].zero();
1241  _ad_d2xyzdeta2_map[p].zero();
1242  }
1243  }
1244 
1245  const unsigned int n_mapping_shape_functions =
1246  FE<3, LAGRANGE>::n_shape_functions(side->type(), side->default_order());
1247 
1248  for (unsigned int i = 0; i < n_mapping_shape_functions; i++)
1249  {
1250  VectorValue<DualReal> side_point = side->point(i);
1251  auto element_node_number = elem->which_node_am_i(side_number, i);
1252 
1253  unsigned dimension = 0;
1254  for (const auto & disp_num : _displacements)
1255  side_point(dimension++)
1256  .derivatives()[disp_num * _sys.getMaxVarNDofsPerElem() + element_node_number] = 1.;
1257 
1258  for (unsigned int p = 0; p < n_qp; p++)
1259  {
1260  _ad_dxyzdxi_map[p].add_scaled(side_point, dpsidxi_map[i][p]);
1261  _ad_dxyzdeta_map[p].add_scaled(side_point, dpsideta_map[i][p]);
1262  if (_calculate_face_xyz)
1263  _ad_q_points_face[p].add_scaled(side_point, psi_map[i][p]);
1265  {
1266  _ad_d2xyzdxi2_map[p].add_scaled(side_point, (*d2psidxi2_map)[i][p]);
1267  _ad_d2xyzdxideta_map[p].add_scaled(side_point, (*d2psidxideta_map)[i][p]);
1268  _ad_d2xyzdeta2_map[p].add_scaled(side_point, (*d2psideta2_map)[i][p]);
1269  }
1270  }
1271  }
1272 
1273  for (unsigned int p = 0; p < n_qp; p++)
1274  {
1275  _ad_normals[p] = _ad_dxyzdxi_map[p].cross(_ad_dxyzdeta_map[p]).unit();
1276 
1277  const auto &dxdxi = _ad_dxyzdxi_map[p](0), dxdeta = _ad_dxyzdeta_map[p](0),
1278  dydxi = _ad_dxyzdxi_map[p](1), dydeta = _ad_dxyzdeta_map[p](1),
1279  dzdxi = _ad_dxyzdxi_map[p](2), dzdeta = _ad_dxyzdeta_map[p](2);
1280 
1281  const auto g11 = (dxdxi * dxdxi + dydxi * dydxi + dzdxi * dzdxi);
1282 
1283  const auto g12 = (dxdxi * dxdeta + dydxi * dydeta + dzdxi * dzdeta);
1284 
1285  const auto g21 = g12;
1286 
1287  const auto g22 = (dxdeta * dxdeta + dydeta * dydeta + dzdeta * dzdeta);
1288 
1289  const auto the_jac = std::sqrt(g11 * g22 - g12 * g21);
1290 
1291  _ad_JxW_face[p] = the_jac * qw[p];
1292 
1294  {
1295  const auto L = -_ad_d2xyzdxi2_map[p] * _ad_normals[p];
1296  const auto M = -_ad_d2xyzdxideta_map[p] * _ad_normals[p];
1297  const auto N = -_ad_d2xyzdeta2_map[p] * _ad_normals[p];
1298  const auto E = _ad_dxyzdxi_map[p].norm_sq();
1299  const auto F = _ad_dxyzdxi_map[p] * _ad_dxyzdeta_map[p];
1300  const auto G = _ad_dxyzdeta_map[p].norm_sq();
1301 
1302  const auto numerator = E * N - 2. * F * M + G * L;
1303  const auto denominator = E * G - F * F;
1304  libmesh_assert_not_equal_to(denominator, 0.);
1305  _ad_curvatures[p] = 0.5 * numerator / denominator;
1306  }
1307  }
1308 
1309  break;
1310  }
1311 
1312  default:
1313  mooseError("Invalid dimension dim = ", dim);
1314  }
1315 }
1316 
1317 void
1318 Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1319 {
1320  unsigned int neighbor_dim = neighbor->dim();
1321 
1322  // reinit neighbor face
1323  for (const auto & it : _fe_face_neighbor[neighbor_dim])
1324  {
1325  FEBase * fe_face_neighbor = it.second;
1326  FEType fe_type = it.first;
1327  FEShapeData * fesd = _fe_shape_data_face_neighbor[fe_type];
1328 
1329  fe_face_neighbor->reinit(neighbor, &reference_points);
1330 
1331  _current_fe_face_neighbor[fe_type] = fe_face_neighbor;
1332 
1333  fesd->_phi.shallowCopy(
1334  const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor->get_phi()));
1335  fesd->_grad_phi.shallowCopy(
1336  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor->get_dphi()));
1338  fesd->_second_phi.shallowCopy(
1339  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor->get_d2phi()));
1340  }
1341  for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
1342  {
1343  FEVectorBase * fe_face_neighbor = it.second;
1344  const FEType & fe_type = it.first;
1345 
1346  _current_vector_fe_face_neighbor[fe_type] = fe_face_neighbor;
1347 
1349 
1350  fe_face_neighbor->reinit(neighbor, &reference_points);
1351 
1352  fesd->_phi.shallowCopy(
1353  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor->get_phi()));
1354  fesd->_grad_phi.shallowCopy(
1355  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor->get_dphi()));
1356  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
1357  fesd->_second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
1358  fe_face_neighbor->get_d2phi()));
1359  if (_need_curl.find(fe_type) != _need_curl.end())
1360  fesd->_curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
1361  fe_face_neighbor->get_curl_phi()));
1362  }
1363 }
1364 
1365 void
1366 Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1367 {
1368  unsigned int neighbor_dim = neighbor->dim();
1369 
1370  // reinit neighbor face
1371  for (const auto & it : _fe_neighbor[neighbor_dim])
1372  {
1373  FEBase * fe_neighbor = it.second;
1374  FEType fe_type = it.first;
1375  FEShapeData * fesd = _fe_shape_data_neighbor[fe_type];
1376 
1377  fe_neighbor->reinit(neighbor, &reference_points);
1378 
1379  _current_fe_neighbor[fe_type] = fe_neighbor;
1380 
1381  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_neighbor->get_phi()));
1382  fesd->_grad_phi.shallowCopy(
1383  const_cast<std::vector<std::vector<RealGradient>> &>(fe_neighbor->get_dphi()));
1385  fesd->_second_phi.shallowCopy(
1386  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor->get_d2phi()));
1387  }
1388  for (const auto & it : _vector_fe_neighbor[neighbor_dim])
1389  {
1390  FEVectorBase * fe_neighbor = it.second;
1391  const FEType & fe_type = it.first;
1392 
1393  _current_vector_fe_neighbor[fe_type] = fe_neighbor;
1394 
1396 
1397  fe_neighbor->reinit(neighbor, &reference_points);
1398 
1399  fesd->_phi.shallowCopy(
1400  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor->get_phi()));
1401  fesd->_grad_phi.shallowCopy(
1402  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor->get_dphi()));
1403  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
1404  fesd->_second_phi.shallowCopy(
1405  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_neighbor->get_d2phi()));
1406  if (_need_curl.find(fe_type) != _need_curl.end())
1407  fesd->_curl_phi.shallowCopy(
1408  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor->get_curl_phi()));
1409  }
1410 }
1411 
1412 void
1413 Assembly::reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
1414 {
1415  unsigned int neighbor_dim = neighbor->dim();
1416 
1417  ArbitraryQuadrature * neighbor_rule = _holder_qrule_neighbor[neighbor_dim];
1418  neighbor_rule->setPoints(reference_points);
1419  setNeighborQRule(neighbor_rule, neighbor_dim);
1420 
1422  mooseAssert(_current_neighbor_subdomain_id == _current_neighbor_elem->subdomain_id(),
1423  "current neighbor subdomain has been set incorrectly");
1424 
1425  // Calculate the volume of the neighbor
1427  {
1428  unsigned int dim = neighbor->dim();
1429  FEBase * fe = *_holder_fe_neighbor_helper[dim];
1430  QBase * qrule = _holder_qrule_volume[dim];
1431 
1432  fe->attach_quadrature_rule(qrule);
1433  fe->reinit(neighbor);
1434 
1435  // set the coord transformation
1436  _coord_neighbor.resize(qrule->n_points());
1437  Moose::CoordinateSystemType coord_type =
1439  unsigned int rz_radial_coord = _subproblem.getAxisymmetricRadialCoord();
1440  const std::vector<Real> & JxW = fe->get_JxW();
1441  const std::vector<Point> & q_points = fe->get_xyz();
1442 
1443  switch (coord_type)
1444  {
1445  case Moose::COORD_XYZ:
1446  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1447  _coord_neighbor[qp] = 1.;
1448  break;
1449 
1450  case Moose::COORD_RZ:
1451  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1452  _coord_neighbor[qp] = 2 * M_PI * q_points[qp](rz_radial_coord);
1453  break;
1454 
1456  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1457  _coord_neighbor[qp] = 4 * M_PI * q_points[qp](0) * q_points[qp](0);
1458  break;
1459 
1460  default:
1461  mooseError("Unknown coordinate system");
1462  break;
1463  }
1464 
1466  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1468  }
1469 }
1470 
1471 template <ComputeStage compute_stage>
1472 void
1474  const ADPoint & q_points,
1475  MooseArray<ADReal> & coord)
1476 {
1477  coord.resize(qrule->n_points());
1479  unsigned int rz_radial_coord = _subproblem.getAxisymmetricRadialCoord();
1480 
1481  switch (_coord_type)
1482  {
1483  case Moose::COORD_XYZ:
1484  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1485  coord[qp] = 1.;
1486  break;
1487 
1488  case Moose::COORD_RZ:
1489  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1490  coord[qp] = 2 * M_PI * q_points[qp](rz_radial_coord);
1491  break;
1492 
1494  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
1495  coord[qp] = 4 * M_PI * q_points[qp](0) * q_points[qp](0);
1496  break;
1497 
1498  default:
1499  mooseError("Unknown coordinate system");
1500  break;
1501  }
1502 }
1503 
1504 template void Assembly::setCoordinateTransformation<ComputeStage::RESIDUAL>(
1505  const QBase *, const MooseArray<Point> &, MooseArray<Real> &);
1506 template void Assembly::setCoordinateTransformation<ComputeStage::JACOBIAN>(
1507  const QBase *, const MooseArray<VectorValue<DualReal>> &, MooseArray<DualReal> &);
1508 
1509 void
1511 {
1513  return;
1514 
1515  setCoordinateTransformation<ComputeStage::RESIDUAL>(_current_qrule, _current_q_points, _coord);
1517  setCoordinateTransformation<ComputeStage::JACOBIAN>(_current_qrule, _ad_q_points, _ad_coord);
1518 
1519  _current_elem_volume = 0.;
1520  for (unsigned int qp = 0; qp < _current_qrule->n_points(); qp++)
1522 
1524 }
1525 
1526 void
1528 {
1530  return;
1531 
1532  setCoordinateTransformation<ComputeStage::RESIDUAL>(
1535  setCoordinateTransformation<ComputeStage::JACOBIAN>(
1537 
1538  _current_side_volume = 0.;
1539  for (unsigned int qp = 0; qp < _current_qrule_face->n_points(); qp++)
1541 
1543 }
1544 
1545 void
1546 Assembly::reinitAtPhysical(const Elem * elem, const std::vector<Point> & physical_points)
1547 {
1548  _current_elem = elem;
1549  _current_neighbor_elem = nullptr;
1550  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1551  "current subdomain has been set incorrectly");
1553 
1554  FEInterface::inverse_map(elem->dim(),
1555  (*_holder_fe_helper[elem->dim()])->get_fe_type(),
1556  elem,
1557  physical_points,
1559 
1561 
1562  // Save off the physical points
1563  _current_physical_points = physical_points;
1564 }
1565 
1566 void
1567 Assembly::reinit(const Elem * elem)
1568 {
1569  _current_elem = elem;
1570  _current_neighbor_elem = nullptr;
1571  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1572  "current subdomain has been set incorrectly");
1574 
1575  unsigned int elem_dimension = elem->dim();
1576 
1577  _current_qrule_volume = _holder_qrule_volume[elem_dimension];
1578 
1579  // Make sure the qrule is the right one
1581  setVolumeQRule(_current_qrule_volume, elem_dimension);
1582 
1583  reinitFE(elem);
1584 
1586 }
1587 
1588 void
1589 Assembly::reinit(const Elem * elem, const std::vector<Point> & reference_points)
1590 {
1591  _current_elem = elem;
1592  _current_neighbor_elem = nullptr;
1593  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1594  "current subdomain has been set incorrectly");
1596 
1597  unsigned int elem_dimension = _current_elem->dim();
1598 
1600 
1601  // Make sure the qrule is the right one
1603  setVolumeQRule(_current_qrule_arbitrary, elem_dimension);
1604 
1605  _current_qrule_arbitrary->setPoints(reference_points);
1606 
1607  reinitFE(elem);
1608 
1610 }
1611 
1612 void
1613 Assembly::reinit(const Elem * elem, unsigned int side)
1614 {
1615  _current_elem = elem;
1616  _current_neighbor_elem = nullptr;
1617  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1618  "current subdomain has been set incorrectly");
1619  _current_side = side;
1622 
1623  unsigned int elem_dimension = _current_elem->dim();
1624 
1625  // Make sure the qrule is the right one
1626  if (_current_qrule_face != _holder_qrule_face[elem_dimension])
1627  {
1628  _current_qrule_face = _holder_qrule_face[elem_dimension];
1629  setFaceQRule(_current_qrule_face, elem_dimension);
1630  }
1631 
1632  if (_current_side_elem)
1633  delete _current_side_elem;
1634  _current_side_elem = elem->build_side_ptr(side).release();
1635 
1637 
1639 }
1640 
1641 void
1642 Assembly::reinit(const Elem * elem, unsigned int side, const std::vector<Point> & reference_points)
1643 {
1644  _current_elem = elem;
1645  _current_neighbor_elem = nullptr;
1646  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
1647  "current subdomain has been set incorrectly");
1648  _current_side = side;
1651 
1652  unsigned int elem_dimension = _current_elem->dim();
1653 
1655 
1656  // Make sure the qrule is the right one
1659 
1660  _current_qrule_arbitrary->setPoints(reference_points);
1661 
1662  if (_current_side_elem)
1663  delete _current_side_elem;
1664  _current_side_elem = elem->build_side_ptr(side).release();
1665 
1667 
1669 }
1670 
1671 void
1672 Assembly::reinit(const Node * node)
1673 {
1674  _current_node = node;
1675  _current_neighbor_node = NULL;
1676 }
1677 
1678 void
1680  unsigned int side,
1681  const Elem * neighbor,
1682  unsigned int neighbor_side)
1683 {
1684  _current_neighbor_side = neighbor_side;
1685 
1686  reinit(elem, side);
1687 
1688  unsigned int neighbor_dim = neighbor->dim();
1689 
1690  std::vector<Point> reference_points;
1691  FEInterface::inverse_map(
1692  neighbor_dim, FEType(), neighbor, _current_q_points_face.stdVector(), reference_points);
1693 
1694  reinitFEFaceNeighbor(neighbor, reference_points);
1695  reinitNeighbor(neighbor, reference_points);
1696 }
1697 
1698 void
1700  unsigned int elem_side,
1701  Real tolerance,
1702  const std::vector<Point> * const pts,
1703  const std::vector<Real> * const weights)
1704 {
1705  _current_elem = elem;
1706 
1707  unsigned int elem_dim = elem->dim();
1708 
1709  // Attach the quadrature rules
1710  if (pts)
1711  {
1712  ArbitraryQuadrature * face_rule = _holder_qrule_arbitrary_face[elem_dim];
1713  face_rule->setPoints(*pts);
1714  setFaceQRule(face_rule, elem_dim);
1715  }
1716  else
1717  {
1718  if (_current_qrule_face != _holder_qrule_face[elem_dim])
1719  {
1721  setFaceQRule(_current_qrule_face, elem_dim);
1722  }
1723  }
1724 
1725  // reinit face
1726  for (const auto & it : _fe_face[elem_dim])
1727  {
1728  FEBase * fe_face = it.second;
1729  FEType fe_type = it.first;
1730  FEShapeData * fesd = _fe_shape_data_face[fe_type];
1731 
1732  fe_face->reinit(elem, elem_side, tolerance, pts, weights);
1733 
1734  _current_fe_face[fe_type] = fe_face;
1735 
1736  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face->get_phi()));
1737  fesd->_grad_phi.shallowCopy(
1738  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face->get_dphi()));
1740  fesd->_second_phi.shallowCopy(
1741  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face->get_d2phi()));
1742  }
1743  for (const auto & it : _vector_fe_face[elem_dim])
1744  {
1745  FEVectorBase * fe_face = it.second;
1746  const FEType & fe_type = it.first;
1747 
1748  _current_vector_fe_face[fe_type] = fe_face;
1749 
1751 
1752  fe_face->reinit(elem, elem_side, tolerance, pts, weights);
1753 
1754  fesd->_phi.shallowCopy(
1755  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face->get_phi()));
1756  fesd->_grad_phi.shallowCopy(
1757  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face->get_dphi()));
1758  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
1759  fesd->_second_phi.shallowCopy(
1760  const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face->get_d2phi()));
1761  if (_need_curl.find(fe_type) != _need_curl.end())
1762  fesd->_curl_phi.shallowCopy(
1763  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face->get_curl_phi()));
1764  }
1765  // During that last loop the helper objects will have been reinitialized
1767  const_cast<std::vector<Point> &>((*_holder_fe_face_helper[elem_dim])->get_xyz()));
1769  const_cast<std::vector<Point> &>((*_holder_fe_face_helper[elem_dim])->get_normals()));
1770  // Note that if the user did pass in points and not weights to this method, JxW will be garbage
1771  // and should not be used
1773  const_cast<std::vector<Real> &>((*_holder_fe_face_helper[elem_dim])->get_JxW()));
1776  const_cast<std::vector<Real> &>((*_holder_fe_face_helper[elem_dim])->get_curvatures()));
1777 
1778  computeADFace(elem, elem_side);
1779 }
1780 
1781 void
1782 Assembly::computeADFace(const Elem * elem, unsigned int side)
1783 {
1784  auto dim = elem->dim();
1785 
1787  {
1788  const std::unique_ptr<const Elem> side_elem(elem->build_side_ptr(side));
1789 
1790  auto n_qp = _current_qrule_face->n_points();
1791  resizeMappingObjects(n_qp, dim);
1792  _ad_normals.resize(n_qp);
1793  _ad_JxW_face.resize(n_qp);
1794  if (_calculate_face_xyz)
1795  _ad_q_points_face.resize(n_qp);
1797  _ad_curvatures.resize(n_qp);
1798 
1799  if (_displaced)
1800  {
1801  const auto & qw = _current_qrule_face->get_weights();
1802  computeFaceMap(dim, qw, side_elem.get());
1803  std::vector<Real> dummy_qw(n_qp, 1.);
1804 
1805  if (elem->has_affine_map())
1806  computeAffineMapAD(elem, dummy_qw, n_qp, *_holder_fe_face_helper[dim]);
1807  else
1808  for (unsigned int qp = 0; qp != n_qp; qp++)
1809  computeSinglePointMapAD(elem, dummy_qw, qp, *_holder_fe_face_helper[dim]);
1810  }
1811  else
1812  for (unsigned qp = 0; qp < n_qp; ++qp)
1813  {
1814  _ad_JxW_face[qp] = _current_JxW_face[qp];
1815  if (_calculate_face_xyz)
1817  _ad_normals[qp] = _current_normals[qp];
1819  _ad_curvatures[qp] = _curvatures[qp];
1820  }
1821 
1822  for (const auto & it : _fe_face[dim])
1823  {
1824  FEBase * fe = it.second;
1825  auto fe_type = it.first;
1826  auto num_shapes = fe->n_shape_functions();
1827  auto & grad_phi = _ad_grad_phi_data_face[fe_type];
1828 
1829  grad_phi.resize(num_shapes);
1830  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
1831  grad_phi[i].resize(n_qp);
1832 
1833  const auto & regular_grad_phi = _fe_shape_data_face[fe_type]->_grad_phi;
1834 
1835  if (_displaced)
1836  computeGradPhiAD(elem, n_qp, grad_phi, fe);
1837  else
1838  for (unsigned qp = 0; qp < n_qp; ++qp)
1839  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
1840  grad_phi[i][qp] = regular_grad_phi[i][qp];
1841  }
1842  for (const auto & it : _vector_fe_face[dim])
1843  {
1844  FEVectorBase * fe = it.second;
1845  auto fe_type = it.first;
1846  auto num_shapes = fe->n_shape_functions();
1847  auto & grad_phi = _ad_vector_grad_phi_data_face[fe_type];
1848 
1849  grad_phi.resize(num_shapes);
1850  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
1851  grad_phi[i].resize(n_qp);
1852 
1853  const auto & regular_grad_phi = _vector_fe_shape_data_face[fe_type]->_grad_phi;
1854 
1855  if (_displaced)
1856  computeGradPhiAD(elem, n_qp, grad_phi, fe);
1857  else
1858  for (unsigned qp = 0; qp < n_qp; ++qp)
1859  for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
1860  grad_phi[i][qp] = regular_grad_phi[i][qp];
1861  }
1862  }
1863 }
1864 
1865 void
1866 Assembly::reinitNeighborFaceRef(const Elem * neighbor,
1867  unsigned int neighbor_side,
1868  Real tolerance,
1869  const std::vector<Point> * const pts,
1870  const std::vector<Real> * const weights)
1871 {
1873 
1874  unsigned int neighbor_dim = neighbor->dim();
1875 
1876  // _holder_qrule_neighbor really does contain pointers to ArbitraryQuadrature
1877  ArbitraryQuadrature * neighbor_rule = _holder_qrule_neighbor[neighbor_dim];
1878  neighbor_rule->setPoints(*pts);
1879 
1880  // Attach this quadrature rule to all the _fe_face_neighbor FE objects. This
1881  // has to have garbage quadrature weights but that's ok because we never
1882  // actually use the JxW coming from these FE reinit'd objects, e.g. we use the
1883  // JxW coming from the element face reinit for DGKernels or we use the JxW
1884  // coming from reinit of the mortar segment element in the case of mortar
1885  setNeighborQRule(neighbor_rule, neighbor_dim);
1886 
1887  // reinit neighbor face
1888  for (const auto & it : _fe_face_neighbor[neighbor_dim])
1889  {
1890  FEBase * fe_face_neighbor = it.second;
1891  FEType fe_type = it.first;
1892  FEShapeData * fesd = _fe_shape_data_face_neighbor[fe_type];
1893 
1894  fe_face_neighbor->reinit(neighbor, neighbor_side, tolerance, pts, weights);
1895 
1896  _current_fe_face_neighbor[fe_type] = fe_face_neighbor;
1897 
1898  fesd->_phi.shallowCopy(
1899  const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor->get_phi()));
1900  fesd->_grad_phi.shallowCopy(
1901  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor->get_dphi()));
1903  fesd->_second_phi.shallowCopy(
1904  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor->get_d2phi()));
1905  }
1906  for (const auto & it : _vector_fe_face_neighbor[neighbor_dim])
1907  {
1908  FEVectorBase * fe_face_neighbor = it.second;
1909  const FEType & fe_type = it.first;
1910 
1911  _current_vector_fe_face_neighbor[fe_type] = fe_face_neighbor;
1912 
1914 
1915  fe_face_neighbor->reinit(neighbor, neighbor_side, tolerance, pts, weights);
1916 
1917  fesd->_phi.shallowCopy(
1918  const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face_neighbor->get_phi()));
1919  fesd->_grad_phi.shallowCopy(
1920  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor->get_dphi()));
1921  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
1922  fesd->_second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
1923  fe_face_neighbor->get_d2phi()));
1924  if (_need_curl.find(fe_type) != _need_curl.end())
1925  fesd->_curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
1926  fe_face_neighbor->get_curl_phi()));
1927  }
1928  // During that last loop the helper objects will have been reinitialized as well
1929  // We need to dig out the q_points from it
1930  _current_q_points_face_neighbor.shallowCopy(const_cast<std::vector<Point> &>(
1931  (*_holder_fe_face_neighbor_helper[neighbor_dim])->get_xyz()));
1932 }
1933 
1934 void
1936  const std::vector<Point> * const pts,
1937  const std::vector<Real> * const weights)
1938 {
1939  mooseAssert(pts->size(),
1940  "Currently reinitialization of lower d elements is only supported with custom "
1941  "quadrature points; there is no fall-back quadrature rule. Consequently make sure "
1942  "you never try to use JxW coming from a fe_lower object unless you are also passing "
1943  "a weights argument");
1944 
1946 
1947  unsigned int elem_dim = elem->dim();
1948 
1949  for (const auto & it : _fe_lower[elem_dim])
1950  {
1951  FEBase * fe_lower = it.second;
1952  FEType fe_type = it.first;
1953  FEShapeData * fesd = _fe_shape_data_lower[fe_type];
1954 
1955  fe_lower->reinit(elem, pts, weights);
1956 
1957  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower->get_phi()));
1958  fesd->_grad_phi.shallowCopy(
1959  const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower->get_dphi()));
1961  fesd->_second_phi.shallowCopy(
1962  const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower->get_d2phi()));
1963  }
1964 }
1965 
1966 void
1967 Assembly::reinitMortarElem(const Elem * elem)
1968 {
1969  mooseAssert(elem->dim() == _mesh_dimension - 1,
1970  "You should be calling reinitMortarElem on a lower dimensional element");
1971 
1972  _fe_msm->reinit(elem);
1973 }
1974 
1975 void
1977  unsigned int neighbor_side,
1978  const std::vector<Point> & physical_points)
1979 {
1981  _current_neighbor_side_elem = neighbor->build_side_ptr(neighbor_side).release();
1982 
1983  std::vector<Point> reference_points;
1984 
1985  unsigned int neighbor_dim = neighbor->dim();
1986  FEInterface::inverse_map(neighbor_dim, FEType(), neighbor, physical_points, reference_points);
1987 
1988  // first do the side element
1990  reinitNeighbor(_current_neighbor_side_elem, reference_points);
1991  // compute JxW on the neighbor's face
1992  unsigned int neighbor_side_dim = _current_neighbor_side_elem->dim();
1993  _current_JxW_neighbor.shallowCopy(const_cast<std::vector<Real> &>(
1994  (*_holder_fe_face_neighbor_helper[neighbor_side_dim])->get_JxW()));
1995 
1996  reinitFEFaceNeighbor(neighbor, reference_points);
1997  reinitNeighbor(neighbor, reference_points);
1998 
1999  // Save off the physical points
2000  _current_physical_points = physical_points;
2001 }
2002 
2003 void
2005  const std::vector<Point> & physical_points)
2006 {
2007  std::vector<Point> reference_points;
2008 
2009  unsigned int neighbor_dim = neighbor->dim();
2010  FEInterface::inverse_map(neighbor_dim, FEType(), neighbor, physical_points, reference_points);
2011 
2012  reinitFENeighbor(neighbor, reference_points);
2013  reinitNeighbor(neighbor, reference_points);
2014  // Save off the physical points
2015  _current_physical_points = physical_points;
2016 }
2017 
2018 DenseMatrix<Number> &
2019 Assembly::jacobianBlock(unsigned int ivar, unsigned int jvar, TagID tag /* = 0 */)
2020 {
2021  _jacobian_block_used[tag][ivar][jvar] = 1;
2022  return _sub_Kee[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
2023 }
2024 
2025 DenseMatrix<Number> &
2026 Assembly::jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar, TagID tag /* = 0*/)
2027 {
2028  _jacobian_block_nonlocal_used[tag][ivar][jvar] = 1;
2029  return _sub_Keg[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
2030 }
2031 
2032 DenseMatrix<Number> &
2034  unsigned int ivar,
2035  unsigned int jvar,
2036  TagID tag /*=0*/)
2037 {
2038  _jacobian_block_neighbor_used[tag][ivar][jvar] = 1;
2040  {
2041  switch (type)
2042  {
2043  default:
2044  case Moose::ElementElement:
2045  return _sub_Kee[tag][ivar][0];
2047  return _sub_Ken[tag][ivar][0];
2049  return _sub_Kne[tag][ivar][0];
2051  return _sub_Knn[tag][ivar][0];
2052  }
2053  }
2054  else
2055  {
2056  switch (type)
2057  {
2058  default:
2059  case Moose::ElementElement:
2060  return _sub_Kee[tag][ivar][jvar];
2062  return _sub_Ken[tag][ivar][jvar];
2064  return _sub_Kne[tag][ivar][jvar];
2066  return _sub_Knn[tag][ivar][jvar];
2067  }
2068  }
2069 }
2070 
2071 DenseMatrix<Number> &
2073  unsigned int ivar,
2074  unsigned int jvar,
2075  TagID tag /*=0*/)
2076 {
2077  _jacobian_block_lower_used[tag][ivar][jvar] = 1;
2079  {
2080  switch (type)
2081  {
2082  default:
2083  case Moose::LowerLower:
2084  return _sub_Kll[tag][ivar][0];
2085  case Moose::LowerSlave:
2086  return _sub_Kle[tag][ivar][0];
2087  case Moose::LowerMaster:
2088  return _sub_Kln[tag][ivar][0];
2089  case Moose::SlaveLower:
2090  return _sub_Kel[tag][ivar][0];
2091  case Moose::SlaveSlave:
2092  return _sub_Kee[tag][ivar][0];
2093  case Moose::SlaveMaster:
2094  return _sub_Ken[tag][ivar][0];
2095  case Moose::MasterLower:
2096  return _sub_Knl[tag][ivar][0];
2097  case Moose::MasterSlave:
2098  return _sub_Kne[tag][ivar][0];
2099  case Moose::MasterMaster:
2100  return _sub_Knn[tag][ivar][0];
2101  }
2102  }
2103  else
2104  {
2105  switch (type)
2106  {
2107  default:
2108  case Moose::LowerLower:
2109  return _sub_Kll[tag][ivar][jvar];
2110  case Moose::LowerSlave:
2111  return _sub_Kle[tag][ivar][jvar];
2112  case Moose::LowerMaster:
2113  return _sub_Kln[tag][ivar][jvar];
2114  case Moose::SlaveLower:
2115  return _sub_Kel[tag][ivar][jvar];
2116  case Moose::SlaveSlave:
2117  return _sub_Kee[tag][ivar][jvar];
2118  case Moose::SlaveMaster:
2119  return _sub_Ken[tag][ivar][jvar];
2120  case Moose::MasterLower:
2121  return _sub_Knl[tag][ivar][jvar];
2122  case Moose::MasterSlave:
2123  return _sub_Kne[tag][ivar][jvar];
2124  case Moose::MasterMaster:
2125  return _sub_Knn[tag][ivar][jvar];
2126  }
2127  }
2128 }
2129 
2130 void
2131 Assembly::init(const CouplingMatrix * cm)
2132 {
2133  _cm = cm;
2134 
2135  unsigned int n_vars = _sys.nVariables();
2136 
2137  _cm_ss_entry.clear();
2138  _cm_sf_entry.clear();
2139  _cm_fs_entry.clear();
2140  _cm_ff_entry.clear();
2141 
2142  auto & vars = _sys.getVariables(_tid);
2143 
2144  _block_diagonal_matrix = true;
2145  for (auto & ivar : vars)
2146  {
2147  auto i = ivar->number();
2148  for (const auto & j : ConstCouplingRow(i, *_cm))
2149  {
2150  if (_sys.isScalarVariable(j))
2151  {
2152  auto & jvar = _sys.getScalarVariable(_tid, j);
2153  _cm_fs_entry.push_back(std::make_pair(ivar, &jvar));
2154  _block_diagonal_matrix = false;
2155  }
2156  else
2157  {
2158  auto & jvar = _sys.getVariable(_tid, j);
2159  _cm_ff_entry.push_back(std::make_pair(ivar, &jvar));
2160  if (i != j)
2161  _block_diagonal_matrix = false;
2162  }
2163  }
2164  }
2165 
2166  auto & scalar_vars = _sys.getScalarVariables(_tid);
2167 
2168  for (auto & ivar : scalar_vars)
2169  {
2170  auto i = ivar->number();
2171  for (const auto & j : ConstCouplingRow(i, *_cm))
2172  if (_sys.isScalarVariable(j))
2173  {
2174  auto & jvar = _sys.getScalarVariable(_tid, j);
2175  _cm_ss_entry.push_back(std::make_pair(ivar, &jvar));
2176  }
2177  else
2178  {
2179  auto & jvar = _sys.getVariable(_tid, j);
2180  _cm_sf_entry.push_back(std::make_pair(ivar, &jvar));
2181  }
2182  }
2183 
2184  if (_block_diagonal_matrix && scalar_vars.size() != 0)
2185  _block_diagonal_matrix = false;
2186 
2187  auto num_vector_tags = _subproblem.numVectorTags();
2188  auto num_matrix_tags = _subproblem.numMatrixTags();
2189 
2190  _sub_Re.resize(num_vector_tags);
2191  _sub_Rn.resize(num_vector_tags);
2192  _sub_Rl.resize(num_vector_tags);
2193  for (MooseIndex(_sub_Re) i = 0; i < _sub_Re.size(); i++)
2194  {
2195  _sub_Re[i].resize(n_vars);
2196  _sub_Rn[i].resize(n_vars);
2197  _sub_Rl[i].resize(n_vars);
2198  }
2199 
2200  _cached_residual_values.resize(num_vector_tags);
2201  _cached_residual_rows.resize(num_vector_tags);
2202 
2203  _cached_jacobian_values.resize(num_matrix_tags);
2204  _cached_jacobian_rows.resize(num_matrix_tags);
2205  _cached_jacobian_cols.resize(num_matrix_tags);
2206 
2207  // Element matrices
2208  _sub_Kee.resize(num_matrix_tags);
2209  _sub_Keg.resize(num_matrix_tags);
2210  _sub_Ken.resize(num_matrix_tags);
2211  _sub_Kne.resize(num_matrix_tags);
2212  _sub_Knn.resize(num_matrix_tags);
2213  _sub_Kll.resize(num_matrix_tags);
2214  _sub_Kle.resize(num_matrix_tags);
2215  _sub_Kln.resize(num_matrix_tags);
2216  _sub_Kel.resize(num_matrix_tags);
2217  _sub_Knl.resize(num_matrix_tags);
2218 
2219  _jacobian_block_used.resize(num_matrix_tags);
2220  _jacobian_block_neighbor_used.resize(num_matrix_tags);
2221  _jacobian_block_lower_used.resize(num_matrix_tags);
2222  _jacobian_block_nonlocal_used.resize(num_matrix_tags);
2223 
2224  for (MooseIndex(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
2225  {
2226  _sub_Keg[tag].resize(n_vars);
2227  _sub_Ken[tag].resize(n_vars);
2228  _sub_Kne[tag].resize(n_vars);
2229  _sub_Knn[tag].resize(n_vars);
2230  _sub_Kee[tag].resize(n_vars);
2231  _sub_Kll[tag].resize(n_vars);
2232  _sub_Kle[tag].resize(n_vars);
2233  _sub_Kln[tag].resize(n_vars);
2234  _sub_Kel[tag].resize(n_vars);
2235  _sub_Knl[tag].resize(n_vars);
2236 
2237  _jacobian_block_used[tag].resize(n_vars);
2238  _jacobian_block_neighbor_used[tag].resize(n_vars);
2239  _jacobian_block_lower_used[tag].resize(n_vars);
2240  _jacobian_block_nonlocal_used[tag].resize(n_vars);
2241  for (MooseIndex(n_vars) i = 0; i < n_vars; ++i)
2242  {
2244  {
2245  _sub_Kee[tag][i].resize(n_vars);
2246  _sub_Keg[tag][i].resize(n_vars);
2247  _sub_Ken[tag][i].resize(n_vars);
2248  _sub_Kne[tag][i].resize(n_vars);
2249  _sub_Knn[tag][i].resize(n_vars);
2250  _sub_Kll[tag][i].resize(n_vars);
2251  _sub_Kle[tag][i].resize(n_vars);
2252  _sub_Kln[tag][i].resize(n_vars);
2253  _sub_Kel[tag][i].resize(n_vars);
2254  _sub_Knl[tag][i].resize(n_vars);
2255  }
2256  else
2257  {
2258  _sub_Kee[tag][i].resize(1);
2259  _sub_Keg[tag][i].resize(1);
2260  _sub_Ken[tag][i].resize(1);
2261  _sub_Kne[tag][i].resize(1);
2262  _sub_Knn[tag][i].resize(1);
2263  _sub_Kll[tag][i].resize(1);
2264  _sub_Kle[tag][i].resize(1);
2265  _sub_Kln[tag][i].resize(1);
2266  _sub_Kel[tag][i].resize(1);
2267  _sub_Knl[tag][i].resize(1);
2268  }
2269  _jacobian_block_used[tag][i].resize(n_vars);
2270  _jacobian_block_neighbor_used[tag][i].resize(n_vars);
2271  _jacobian_block_lower_used[tag][i].resize(n_vars);
2272  _jacobian_block_nonlocal_used[tag][i].resize(n_vars);
2273  }
2274  }
2275 
2276  // Cached Jacobian contributions
2277  _cached_jacobian_contribution_vals.resize(num_matrix_tags);
2278  _cached_jacobian_contribution_rows.resize(num_matrix_tags);
2279  _cached_jacobian_contribution_cols.resize(num_matrix_tags);
2280 }
2281 
2282 void
2284 {
2285  _cm_nonlocal_entry.clear();
2286 
2287  auto & vars = _sys.getVariables(_tid);
2288 
2289  for (auto & ivar : vars)
2290  {
2291  auto i = ivar->number();
2292  for (const auto & j : ConstCouplingRow(i, _nonlocal_cm))
2293  if (!_sys.isScalarVariable(j))
2294  {
2295  auto & jvar = _sys.getVariable(_tid, j);
2296  _cm_nonlocal_entry.push_back(std::make_pair(ivar, &jvar));
2297  }
2298  }
2299 }
2300 
2301 void
2303 {
2304  for (const auto & it : _cm_ff_entry)
2305  {
2306  MooseVariableFEBase & ivar = *(it.first);
2307  MooseVariableFEBase & jvar = *(it.second);
2308 
2309  unsigned int vi = ivar.number();
2310  unsigned int vj = jvar.number();
2311 
2312  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2313  {
2314  jacobianBlock(vi, vj, tag).resize(ivar.dofIndices().size(), jvar.dofIndices().size());
2315  jacobianBlock(vi, vj, tag).zero();
2316  _jacobian_block_used[tag][vi][vj] = 0;
2317  }
2318  }
2319 }
2320 
2321 void
2323 {
2324  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2325  for (const auto & var : vars)
2326  for (MooseIndex(_sub_Re) tag = 0; tag < _sub_Re.size(); tag++)
2327  {
2328  _sub_Re[tag][var->number()].resize(var->dofIndices().size());
2329  _sub_Re[tag][var->number()].zero();
2330  }
2331 }
2332 
2333 void
2335 {
2337  prepareResidual();
2338 }
2339 
2340 void
2342 {
2343  for (const auto & it : _cm_nonlocal_entry)
2344  {
2345  MooseVariableFEBase & ivar = *(it.first);
2346  MooseVariableFEBase & jvar = *(it.second);
2347 
2348  unsigned int vi = ivar.number();
2349  unsigned int vj = jvar.number();
2350 
2351  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2352  tag < _jacobian_block_nonlocal_used.size();
2353  tag++)
2354  {
2355  jacobianBlockNonlocal(vi, vj, tag)
2356  .resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
2357  jacobianBlockNonlocal(vi, vj, tag).zero();
2358  _jacobian_block_nonlocal_used[tag][vi][vj] = 0;
2359  }
2360  }
2361 }
2362 
2363 void
2365 {
2366  for (const auto & it : _cm_ff_entry)
2367  {
2368  MooseVariableFEBase & ivar = *(it.first);
2369  MooseVariableFEBase & jvar = *(it.second);
2370 
2371  unsigned int vi = ivar.number();
2372  unsigned int vj = jvar.number();
2373 
2374  if (vi == var->number() || vj == var->number())
2375  {
2376  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2377  {
2378  jacobianBlock(vi, vj, tag).resize(ivar.dofIndices().size(), jvar.dofIndices().size());
2379  jacobianBlock(vi, vj, tag).zero();
2380  _jacobian_block_used[tag][vi][vj] = 0;
2381  }
2382  }
2383  }
2384 
2385  for (MooseIndex(_sub_Re) tag = 0; tag < _sub_Re.size(); tag++)
2386  {
2387  _sub_Re[tag][var->number()].resize(var->dofIndices().size());
2388  _sub_Re[tag][var->number()].zero();
2389  }
2390 }
2391 
2392 void
2394 {
2395  for (const auto & it : _cm_nonlocal_entry)
2396  {
2397  MooseVariableFEBase & ivar = *(it.first);
2398  MooseVariableFEBase & jvar = *(it.second);
2399 
2400  unsigned int vi = ivar.number();
2401  unsigned int vj = jvar.number();
2402 
2403  if (vi == var->number() || vj == var->number())
2404  {
2405  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2406  tag < _jacobian_block_nonlocal_used.size();
2407  tag++)
2408  {
2409  jacobianBlockNonlocal(vi, vj, tag)
2410  .resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
2411  jacobianBlockNonlocal(vi, vj, tag).zero();
2412  _jacobian_block_nonlocal_used[tag][vi][vj] = 0;
2413  }
2414  }
2415  }
2416 }
2417 
2418 void
2420 {
2421  for (const auto & it : _cm_ff_entry)
2422  {
2423  MooseVariableFEBase & ivar = *(it.first);
2424  MooseVariableFEBase & jvar = *(it.second);
2425 
2426  unsigned int vi = ivar.number();
2427  unsigned int vj = jvar.number();
2428 
2429  for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
2430  tag < _jacobian_block_neighbor_used.size();
2431  tag++)
2432  {
2434  .resize(ivar.dofIndices().size(), jvar.dofIndicesNeighbor().size());
2435  jacobianBlockNeighbor(Moose::ElementNeighbor, vi, vj, tag).zero();
2436 
2438  .resize(ivar.dofIndicesNeighbor().size(), jvar.dofIndices().size());
2439  jacobianBlockNeighbor(Moose::NeighborElement, vi, vj, tag).zero();
2440 
2442  .resize(ivar.dofIndicesNeighbor().size(), jvar.dofIndicesNeighbor().size());
2443  jacobianBlockNeighbor(Moose::NeighborNeighbor, vi, vj, tag).zero();
2444 
2445  _jacobian_block_neighbor_used[tag][vi][vj] = 0;
2446  }
2447  }
2448 
2449  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2450  for (const auto & var : vars)
2451  for (MooseIndex(_sub_Rn) tag = 0; tag < _sub_Rn.size(); tag++)
2452  {
2453  _sub_Rn[tag][var->number()].resize(var->dofIndicesNeighbor().size());
2454  _sub_Rn[tag][var->number()].zero();
2455  }
2456 }
2457 
2458 void
2460 {
2461  for (const auto & it : _cm_ff_entry)
2462  {
2463  MooseVariableFEBase & ivar = *(it.first);
2464  MooseVariableFEBase & jvar = *(it.second);
2465 
2466  unsigned int vi = ivar.number();
2467  unsigned int vj = jvar.number();
2468 
2469  for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
2470  tag++)
2471  {
2472  // To cover all possible cases we should have 9 combinations below for every 2-permutation of
2473  // Lower,Slave,Master. However, 4 cases will in general be covered by calls to prepare() and
2474  // prepareNeighbor(). These calls will cover SlaveSlave (ElementElement), SlaveMaster
2475  // (ElementNeighbor), MasterSlave (NeighborElement), and MasterMaster (NeighborNeighbor). With
2476  // these covered we only need to prepare the 5 remaining below
2477 
2478  // derivatives w.r.t. lower dimensional residuals
2480  .resize(ivar.dofIndicesLower().size(), jvar.dofIndicesLower().size());
2481  jacobianBlockLower(Moose::LowerLower, vi, vj, tag).zero();
2482 
2484  .resize(ivar.dofIndicesLower().size(), jvar.dofIndices().size());
2485  jacobianBlockLower(Moose::LowerSlave, vi, vj, tag).zero();
2486 
2488  .resize(ivar.dofIndicesLower().size(), jvar.dofIndicesNeighbor().size());
2489  jacobianBlockLower(Moose::LowerMaster, vi, vj, tag).zero();
2490 
2491  // derivatives w.r.t. interior slave residuals
2493  .resize(ivar.dofIndices().size(), jvar.dofIndicesLower().size());
2494  jacobianBlockLower(Moose::SlaveLower, vi, vj, tag).zero();
2495 
2496  // derivatives w.r.t. interior master residuals
2498  .resize(ivar.dofIndicesNeighbor().size(), jvar.dofIndicesLower().size());
2499  jacobianBlockLower(Moose::MasterLower, vi, vj, tag).zero();
2500 
2501  _jacobian_block_lower_used[tag][vi][vj] = 0;
2502  }
2503  }
2504 
2505  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2506  for (const auto & var : vars)
2507  for (MooseIndex(_sub_Rl) tag = 0; tag < _sub_Rl.size(); tag++)
2508  {
2509  _sub_Rl[tag][var->number()].resize(var->dofIndicesLower().size());
2510  _sub_Rl[tag][var->number()].zero();
2511  }
2512 }
2513 
2514 void
2515 Assembly::prepareBlock(unsigned int ivar,
2516  unsigned int jvar,
2517  const std::vector<dof_id_type> & dof_indices)
2518 {
2519  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2520  {
2521  jacobianBlock(ivar, jvar, tag).resize(dof_indices.size(), dof_indices.size());
2522  jacobianBlock(ivar, jvar, tag).zero();
2523  _jacobian_block_used[tag][ivar][jvar] = 0;
2524  }
2525 
2526  for (MooseIndex(_sub_Re) tag = 0; tag < _sub_Re.size(); tag++)
2527  {
2528  _sub_Re[tag][ivar].resize(dof_indices.size());
2529  _sub_Re[tag][ivar].zero();
2530  }
2531 }
2532 
2533 void
2535  unsigned int jvar,
2536  const std::vector<dof_id_type> & idof_indices,
2537  const std::vector<dof_id_type> & jdof_indices)
2538 {
2539  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
2540  tag < _jacobian_block_nonlocal_used.size();
2541  tag++)
2542  {
2543  jacobianBlockNonlocal(ivar, jvar, tag).resize(idof_indices.size(), jdof_indices.size());
2544  jacobianBlockNonlocal(ivar, jvar, tag).zero();
2545  _jacobian_block_nonlocal_used[tag][ivar][jvar] = 0;
2546  }
2547 }
2548 
2549 void
2551 {
2552  const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
2553  for (const auto & ivar : vars)
2554  {
2555  auto idofs = ivar->dofIndices().size();
2556 
2557  for (MooseIndex(_sub_Re) tag = 0; tag < _sub_Re.size(); tag++)
2558  {
2559  _sub_Re[tag][ivar->number()].resize(idofs);
2560  _sub_Re[tag][ivar->number()].zero();
2561  }
2562 
2563  for (const auto & jvar : vars)
2564  {
2565  auto jdofs = jvar->dofIndices().size();
2566 
2567  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2568  {
2569  jacobianBlock(ivar->number(), jvar->number(), tag).resize(idofs, jdofs);
2570  jacobianBlock(ivar->number(), jvar->number(), tag).zero();
2571  _jacobian_block_used[tag][ivar->number()][jvar->number()] = 0;
2572  }
2573  }
2574  }
2575 }
2576 
2577 void
2579 {
2580  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2581  const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
2582 
2583  for (const auto & ivar : scalar_vars)
2584  {
2585  auto idofs = ivar->dofIndices().size();
2586 
2587  for (const auto & jvar : vars)
2588  {
2589  auto jdofs = jvar->dofIndices().size();
2590  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
2591  {
2592  jacobianBlock(ivar->number(), jvar->number(), tag).resize(idofs, jdofs);
2593  jacobianBlock(ivar->number(), jvar->number(), tag).zero();
2594  _jacobian_block_used[tag][ivar->number()][jvar->number()] = 0;
2595 
2596  jacobianBlock(jvar->number(), ivar->number(), tag).resize(jdofs, idofs);
2597  jacobianBlock(jvar->number(), ivar->number(), tag).zero();
2598  _jacobian_block_used[tag][jvar->number()][ivar->number()] = 0;
2599  }
2600  }
2601  }
2602 }
2603 
2604 template <typename T>
2605 void
2607 {
2608  phi(v).shallowCopy(v.phi());
2609  gradPhi(v).shallowCopy(v.gradPhi());
2610  if (v.computingSecond())
2611  secondPhi(v).shallowCopy(v.secondPhi());
2612 }
2613 
2614 void
2615 Assembly::copyShapes(unsigned int var)
2616 {
2617  try
2618  {
2619  MooseVariable & v = _sys.getFieldVariable<Real>(_tid, var);
2620  copyShapes(v);
2621  }
2622  catch (std::out_of_range & e)
2623  {
2625  copyShapes(v);
2626  if (v.computingCurl())
2627  curlPhi(v).shallowCopy(v.curlPhi());
2628  }
2629 }
2630 
2631 template <typename T>
2632 void
2634 {
2635  phiFace(v).shallowCopy(v.phiFace());
2636  gradPhiFace(v).shallowCopy(v.gradPhiFace());
2637  if (v.computingSecond())
2638  secondPhiFace(v).shallowCopy(v.secondPhiFace());
2639 }
2640 
2641 void
2642 Assembly::copyFaceShapes(unsigned int var)
2643 {
2644  try
2645  {
2646  MooseVariable & v = _sys.getFieldVariable<Real>(_tid, var);
2647  copyFaceShapes(v);
2648  }
2649  catch (std::out_of_range & e)
2650  {
2652  copyFaceShapes(v);
2653  if (v.computingCurl())
2654  _vector_curl_phi_face.shallowCopy(v.curlPhi());
2655  }
2656 }
2657 
2658 template <typename T>
2659 void
2661 {
2662  if (v.usesPhiNeighbor())
2663  {
2664  phiFaceNeighbor(v).shallowCopy(v.phiFaceNeighbor());
2665  phiNeighbor(v).shallowCopy(v.phiNeighbor());
2666  }
2667  if (v.usesGradPhiNeighbor())
2668  {
2669  gradPhiFaceNeighbor(v).shallowCopy(v.gradPhiFaceNeighbor());
2670  gradPhiNeighbor(v).shallowCopy(v.gradPhiNeighbor());
2671  }
2672  if (v.usesSecondPhiNeighbor())
2673  {
2674  secondPhiFaceNeighbor(v).shallowCopy(v.secondPhiFaceNeighbor());
2675  secondPhiNeighbor(v).shallowCopy(v.secondPhiNeighbor());
2676  }
2677 }
2678 
2679 void
2681 {
2682  try
2683  {
2684  MooseVariable & v = _sys.getFieldVariable<Real>(_tid, var);
2685  copyNeighborShapes(v);
2686  }
2687  catch (std::out_of_range & e)
2688  {
2690  copyNeighborShapes(v);
2691  }
2692 }
2693 
2694 void
2695 Assembly::addResidualBlock(NumericVector<Number> & residual,
2696  DenseVector<Number> & res_block,
2697  const std::vector<dof_id_type> & dof_indices,
2698  Real scaling_factor)
2699 {
2700  if (dof_indices.size() > 0 && res_block.size())
2701  {
2702  _temp_dof_indices = dof_indices;
2703  _dof_map.constrain_element_vector(res_block, _temp_dof_indices, false);
2704 
2705  if (scaling_factor != 1.0)
2706  {
2707  _tmp_Re = res_block;
2708  _tmp_Re *= scaling_factor;
2709  residual.add_vector(_tmp_Re, _temp_dof_indices);
2710  }
2711  else
2712  {
2713  residual.add_vector(res_block, _temp_dof_indices);
2714  }
2715  }
2716 }
2717 
2718 void
2719 Assembly::cacheResidualBlock(std::vector<Real> & cached_residual_values,
2720  std::vector<dof_id_type> & cached_residual_rows,
2721  DenseVector<Number> & res_block,
2722  const std::vector<dof_id_type> & dof_indices,
2723  Real scaling_factor)
2724 {
2725  if (dof_indices.size() > 0 && res_block.size())
2726  {
2727  _temp_dof_indices = dof_indices;
2728  _dof_map.constrain_element_vector(res_block, _temp_dof_indices, false);
2729 
2730  if (scaling_factor != 1.0)
2731  {
2732  _tmp_Re = res_block;
2733  _tmp_Re *= scaling_factor;
2734 
2735  for (MooseIndex(_tmp_Re) i = 0; i < _tmp_Re.size(); i++)
2736  {
2737  cached_residual_values.push_back(_tmp_Re(i));
2738  cached_residual_rows.push_back(_temp_dof_indices[i]);
2739  }
2740  }
2741  else
2742  {
2743  for (MooseIndex(res_block) i = 0; i < res_block.size(); i++)
2744  {
2745  cached_residual_values.push_back(res_block(i));
2746  cached_residual_rows.push_back(_temp_dof_indices[i]);
2747  }
2748  }
2749  }
2750 
2751  res_block.zero();
2752 }
2753 
2754 void
2755 Assembly::addResidual(NumericVector<Number> & residual, TagID tag_id /* = 0 */)
2756 {
2757  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2758  for (const auto & var : vars)
2760  residual, _sub_Re[tag_id][var->number()], var->dofIndices(), var->scalingFactor());
2761 }
2762 
2763 void
2764 Assembly::addResidual(const std::map<TagName, TagID> & tags)
2765 {
2766  for (auto & tag : tags)
2767  if (_sys.hasVector(tag.second))
2768  addResidual(_sys.getVector(tag.second), tag.second);
2769 }
2770 
2771 void
2772 Assembly::addResidualNeighbor(NumericVector<Number> & residual, TagID tag_id /* = 0 */)
2773 {
2774  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2775  for (const auto & var : vars)
2777  residual, _sub_Rn[tag_id][var->number()], var->dofIndicesNeighbor(), var->scalingFactor());
2778 }
2779 
2780 void
2781 Assembly::addResidualNeighbor(const std::map<TagName, TagID> & tags)
2782 {
2783  for (auto & tag : tags)
2784  if (_sys.hasVector(tag.second))
2785  addResidualNeighbor(_sys.getVector(tag.second), tag.second);
2786 }
2787 
2788 void
2790 {
2791  // add the scalar variables residuals
2792  const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
2793  for (const auto & var : vars)
2794  if (_sys.hasVector(tag_id))
2796  _sub_Re[tag_id][var->number()],
2797  var->dofIndices(),
2798  var->scalingFactor());
2799 }
2800 
2801 void
2802 Assembly::addResidualScalar(const std::map<TagName, TagID> & tags)
2803 {
2804  for (auto & tag : tags)
2805  addResidualScalar(tag.second);
2806 }
2807 
2808 void
2810 {
2811  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2812  for (const auto & var : vars)
2813  {
2814  for (MooseIndex(_cached_residual_values) tag = 0; tag < _cached_residual_values.size(); tag++)
2815  if (_sys.hasVector(tag))
2817  _cached_residual_rows[tag],
2818  _sub_Re[tag][var->number()],
2819  var->dofIndices(),
2820  var->scalingFactor());
2821  }
2822 }
2823 
2824 void
2825 Assembly::cacheResidualContribution(dof_id_type dof, Real value, TagID tag_id)
2826 {
2827  _cached_residual_values[tag_id].push_back(value);
2828  _cached_residual_rows[tag_id].push_back(dof);
2829 }
2830 
2831 void
2832 Assembly::cacheResidualContribution(dof_id_type dof, Real value, const std::set<TagID> & tags)
2833 {
2834  for (auto & tag : tags)
2835  cacheResidualContribution(dof, value, tag);
2836 }
2837 
2838 void
2840 {
2841  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2842  for (const auto & var : vars)
2843  {
2844  for (MooseIndex(_cached_residual_values) tag = 0; tag < _cached_residual_values.size(); tag++)
2845  {
2846  if (_sys.hasVector(tag))
2848  _cached_residual_rows[tag],
2849  _sub_Rn[tag][var->number()],
2850  var->dofIndicesNeighbor(),
2851  var->scalingFactor());
2852  }
2853  }
2854 }
2855 
2856 void
2858 {
2859  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2860  for (const auto & var : vars)
2861  {
2862  for (MooseIndex(_cached_residual_values) tag = 0; tag < _cached_residual_values.size(); tag++)
2863  {
2864  if (_sys.hasVector(tag))
2866  _cached_residual_rows[tag],
2867  _sub_Rl[tag][var->number()],
2868  var->dofIndicesLower(),
2869  var->scalingFactor());
2870  }
2871  }
2872 }
2873 
2874 void
2875 Assembly::cacheResidualNodes(const DenseVector<Number> & res,
2876  const std::vector<dof_id_type> & dof_index,
2877  TagID tag /* = 0*/)
2878 {
2879  // Add the residual value and dof_index to cached_residual_values and cached_residual_rows
2880  // respectively.
2881  // This is used by NodalConstraint.C to cache the residual calculated for master and slave node.
2882  for (MooseIndex(dof_index) i = 0; i < dof_index.size(); ++i)
2883  {
2884  _cached_residual_values[tag].push_back(res(i));
2885  _cached_residual_rows[tag].push_back(dof_index[i]);
2886  }
2887 }
2888 
2889 void
2891 {
2892  for (MooseIndex(_cached_residual_values) tag = 0; tag < _cached_residual_values.size(); tag++)
2893  {
2894  if (!_sys.hasVector(tag))
2895  {
2896  _cached_residual_values[tag].clear();
2897  _cached_residual_rows[tag].clear();
2898  continue;
2899  }
2900  addCachedResidual(_sys.getVector(tag), tag);
2901  }
2902 }
2903 
2904 void
2905 Assembly::addCachedResidual(NumericVector<Number> & residual, TagID tag_id)
2906 {
2907  if (!_sys.hasVector(tag_id))
2908  {
2909  // Only clean up things when tag exists
2910  if (_subproblem.vectorTagExists(tag_id))
2911  {
2912  _cached_residual_values[tag_id].clear();
2913  _cached_residual_rows[tag_id].clear();
2914  }
2915  return;
2916  }
2917 
2918  std::vector<Real> & cached_residual_values = _cached_residual_values[tag_id];
2919  std::vector<dof_id_type> & cached_residual_rows = _cached_residual_rows[tag_id];
2920 
2921  mooseAssert(cached_residual_values.size() == cached_residual_rows.size(),
2922  "Number of cached residuals and number of rows must match!");
2923  if (cached_residual_values.size())
2924  residual.add_vector(cached_residual_values, cached_residual_rows);
2925 
2926  if (_max_cached_residuals < cached_residual_values.size())
2927  _max_cached_residuals = cached_residual_values.size();
2928 
2929  // Try to be more efficient from now on
2930  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
2931  cached_residual_values.clear();
2932  cached_residual_values.reserve(_max_cached_residuals * 2);
2933 
2934  cached_residual_rows.clear();
2935  cached_residual_rows.reserve(_max_cached_residuals * 2);
2936 }
2937 
2938 void
2939 Assembly::setResidualBlock(NumericVector<Number> & residual,
2940  DenseVector<Number> & res_block,
2941  const std::vector<dof_id_type> & dof_indices,
2942  Real scaling_factor)
2943 {
2944  if (dof_indices.size() > 0)
2945  {
2946  std::vector<dof_id_type> di(dof_indices);
2947  _dof_map.constrain_element_vector(res_block, di, false);
2948 
2949  if (scaling_factor != 1.0)
2950  {
2951  _tmp_Re = res_block;
2952  _tmp_Re *= scaling_factor;
2953  residual.insert(_tmp_Re, di);
2954  }
2955  else
2956  residual.insert(res_block, di);
2957  }
2958 }
2959 
2960 void
2961 Assembly::setResidual(NumericVector<Number> & residual, TagID tag_id /* = 0 */)
2962 {
2963  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2964  for (const auto & var : vars)
2966  residual, _sub_Re[tag_id][var->number()], var->dofIndices(), var->scalingFactor());
2967 }
2968 
2969 void
2970 Assembly::setResidualNeighbor(NumericVector<Number> & residual, TagID tag_id /* = 0 */)
2971 {
2972  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
2973  for (const auto & var : vars)
2975  residual, _sub_Rn[tag_id][var->number()], var->dofIndicesNeighbor(), var->scalingFactor());
2976 }
2977 
2978 void
2979 Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
2980  DenseMatrix<Number> & jac_block,
2981  const std::vector<dof_id_type> & idof_indices,
2982  const std::vector<dof_id_type> & jdof_indices,
2983  Real scaling_factor)
2984 {
2985  if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m())
2986  {
2987  std::vector<dof_id_type> di(idof_indices);
2988  std::vector<dof_id_type> dj(jdof_indices);
2989  _dof_map.constrain_element_matrix(jac_block, di, dj, false);
2990 
2991  if (scaling_factor != 1.0)
2992  {
2993  _tmp_Ke = jac_block;
2994  _tmp_Ke *= scaling_factor;
2995  jacobian.add_matrix(_tmp_Ke, di, dj);
2996  }
2997  else
2998  jacobian.add_matrix(jac_block, di, dj);
2999  }
3000 }
3001 
3002 void
3003 Assembly::cacheJacobianBlock(DenseMatrix<Number> & jac_block,
3004  const std::vector<dof_id_type> & idof_indices,
3005  const std::vector<dof_id_type> & jdof_indices,
3006  Real scaling_factor,
3007  TagID tag /*=0*/)
3008 {
3009  // Only cache data when the matrix exists
3010  if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m() &&
3011  _sys.hasMatrix(tag))
3012  {
3013  std::vector<dof_id_type> di(idof_indices);
3014  std::vector<dof_id_type> dj(jdof_indices);
3015  _dof_map.constrain_element_matrix(jac_block, di, dj, false);
3016 
3017  if (scaling_factor != 1.0)
3018  jac_block *= scaling_factor;
3019 
3020  for (MooseIndex(di) i = 0; i < di.size(); i++)
3021  for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3022  {
3023  _cached_jacobian_values[tag].push_back(jac_block(i, j));
3024  _cached_jacobian_rows[tag].push_back(di[i]);
3025  _cached_jacobian_cols[tag].push_back(dj[j]);
3026  }
3027  }
3028  jac_block.zero();
3029 }
3030 
3031 void
3032 Assembly::cacheJacobianBlockNonlocal(DenseMatrix<Number> & jac_block,
3033  const std::vector<dof_id_type> & idof_indices,
3034  const std::vector<dof_id_type> & jdof_indices,
3035  Real scaling_factor,
3036  TagID tag /*= 0*/)
3037 {
3038  if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m() &&
3039  _sys.hasMatrix(tag))
3040  {
3041  std::vector<dof_id_type> di(idof_indices);
3042  std::vector<dof_id_type> dj(jdof_indices);
3043  _dof_map.constrain_element_matrix(jac_block, di, dj, false);
3044 
3045  if (scaling_factor != 1.0)
3046  jac_block *= scaling_factor;
3047 
3048  for (MooseIndex(di) i = 0; i < di.size(); i++)
3049  for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3050  if (jac_block(i, j) != 0.0) // no storage allocated for unimplemented jacobian terms,
3051  // maintaining maximum sparsity possible
3052  {
3053  _cached_jacobian_values[tag].push_back(jac_block(i, j));
3054  _cached_jacobian_rows[tag].push_back(di[i]);
3055  _cached_jacobian_cols[tag].push_back(dj[j]);
3056  }
3057  }
3058  jac_block.zero();
3059 }
3060 
3061 void
3062 Assembly::addCachedJacobian(SparseMatrix<Number> & /*jacobian*/)
3063 {
3064  mooseDeprecated(" Please use addCachedJacobian() ");
3065 
3067 }
3068 
3069 void
3071 {
3073  {
3074  mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
3075  "Error: Cached data sizes MUST be the same!");
3076  for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3077  mooseAssert(_cached_jacobian_rows[i].size() == _cached_jacobian_cols[i].size(),
3078  "Error: Cached data sizes MUST be the same for a given tag!");
3079  }
3080 
3081  for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3082  if (_sys.hasMatrix(i))
3083  for (MooseIndex(_cached_jacobian_rows[i]) j = 0; j < _cached_jacobian_rows[i].size(); j++)
3084  _sys.getMatrix(i).add(_cached_jacobian_rows[i][j],
3085  _cached_jacobian_cols[i][j],
3086  _cached_jacobian_values[i][j]);
3087 
3088  for (MooseIndex(_cached_jacobian_rows) i = 0; i < _cached_jacobian_rows.size(); i++)
3089  {
3090  if (!_sys.hasMatrix(i))
3091  continue;
3092 
3095 
3096  // Try to be more efficient from now on
3097  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
3098  _cached_jacobian_values[i].clear();
3100 
3101  _cached_jacobian_rows[i].clear();
3103 
3104  _cached_jacobian_cols[i].clear();
3106  }
3107 }
3108 
3109 inline void
3111 {
3112  auto i = ivar->number();
3113  auto j = jvar->number();
3114  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
3115  if (_jacobian_block_used[tag][i][j] && _sys.hasMatrix(tag))
3117  jacobianBlock(i, j, tag),
3118  ivar->dofIndices(),
3119  jvar->dofIndices(),
3120  ivar->scalingFactor());
3121 }
3122 
3123 void
3125 {
3126  for (const auto & it : _cm_ff_entry)
3127  addJacobianCoupledVarPair(it.first, it.second);
3128 
3129  for (const auto & it : _cm_sf_entry)
3130  addJacobianCoupledVarPair(it.first, it.second);
3131 
3132  for (const auto & it : _cm_fs_entry)
3133  addJacobianCoupledVarPair(it.first, it.second);
3134 }
3135 
3136 void
3138 {
3139  for (const auto & it : _cm_nonlocal_entry)
3140  {
3141  auto ivar = it.first;
3142  auto jvar = it.second;
3143  auto i = ivar->number();
3144  auto j = jvar->number();
3145  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
3146  tag < _jacobian_block_nonlocal_used.size();
3147  tag++)
3148  if (_jacobian_block_nonlocal_used[tag][i][j] && _sys.hasMatrix(tag))
3150  jacobianBlockNonlocal(i, j, tag),
3151  ivar->dofIndices(),
3152  jvar->allDofIndices(),
3153  ivar->scalingFactor());
3154  }
3155 }
3156 
3157 void
3159 {
3160  for (const auto & it : _cm_ff_entry)
3161  {
3162  auto ivar = it.first;
3163  auto jvar = it.second;
3164  auto i = ivar->number();
3165  auto j = jvar->number();
3166  for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3167  tag < _jacobian_block_neighbor_used.size();
3168  tag++)
3169  if (_jacobian_block_neighbor_used[tag][i][j] && _sys.hasMatrix(tag))
3170  {
3173  ivar->dofIndices(),
3174  jvar->dofIndicesNeighbor(),
3175  ivar->scalingFactor());
3176 
3179  ivar->dofIndicesNeighbor(),
3180  jvar->dofIndices(),
3181  ivar->scalingFactor());
3182 
3185  ivar->dofIndicesNeighbor(),
3186  jvar->dofIndicesNeighbor(),
3187  ivar->scalingFactor());
3188  }
3189  }
3190 }
3191 
3192 void
3194 {
3195  for (const auto & it : _cm_ff_entry)
3196  {
3197  auto ivar = it.first;
3198  auto jvar = it.second;
3199  auto i = ivar->number();
3200  auto j = jvar->number();
3201  for (MooseIndex(_jacobian_block_lower_used) tag = 0; tag < _jacobian_block_lower_used.size();
3202  tag++)
3203  if (_jacobian_block_lower_used[tag][i][j] && _sys.hasMatrix(tag))
3204  {
3207  ivar->dofIndicesLower(),
3208  jvar->dofIndicesLower(),
3209  ivar->scalingFactor());
3210 
3213  ivar->dofIndicesLower(),
3214  jvar->dofIndices(),
3215  ivar->scalingFactor());
3216 
3219  ivar->dofIndicesLower(),
3220  jvar->dofIndicesNeighbor(),
3221  ivar->scalingFactor());
3222 
3225  ivar->dofIndices(),
3226  jvar->dofIndicesLower(),
3227  ivar->scalingFactor());
3228 
3231  ivar->dofIndicesNeighbor(),
3232  jvar->dofIndicesLower(),
3233  ivar->scalingFactor());
3234  }
3235  }
3236 }
3237 
3238 inline void
3240 {
3241  auto i = ivar->number();
3242  auto j = jvar->number();
3243  for (MooseIndex(_jacobian_block_used) tag = 0; tag < _jacobian_block_used.size(); tag++)
3244  if (_jacobian_block_used[tag][i][j] && _sys.hasMatrix(tag))
3245  cacheJacobianBlock(jacobianBlock(i, j, tag),
3246  ivar->dofIndices(),
3247  jvar->dofIndices(),
3248  ivar->scalingFactor(),
3249  tag);
3250 }
3251 
3252 void
3254 {
3255  for (const auto & it : _cm_ff_entry)
3256  cacheJacobianCoupledVarPair(it.first, it.second);
3257 
3258  for (const auto & it : _cm_fs_entry)
3259  cacheJacobianCoupledVarPair(it.first, it.second);
3260 
3261  for (const auto & it : _cm_sf_entry)
3262  cacheJacobianCoupledVarPair(it.first, it.second);
3263 }
3264 
3265 void
3267 {
3268  for (const auto & it : _cm_nonlocal_entry)
3269  {
3270  auto ivar = it.first;
3271  auto jvar = it.second;
3272  auto i = ivar->number();
3273  auto j = jvar->number();
3274  for (MooseIndex(_jacobian_block_nonlocal_used) tag = 0;
3275  tag < _jacobian_block_nonlocal_used.size();
3276  tag++)
3277  if (_jacobian_block_nonlocal_used[tag][i][j] && _sys.hasMatrix(tag))
3279  ivar->dofIndices(),
3280  jvar->allDofIndices(),
3281  ivar->scalingFactor(),
3282  tag);
3283  }
3284 }
3285 
3286 void
3288 {
3289  for (const auto & it : _cm_ff_entry)
3290  {
3291  auto ivar = it.first;
3292  auto jvar = it.second;
3293  auto i = ivar->number();
3294  auto j = jvar->number();
3295 
3296  for (MooseIndex(_jacobian_block_neighbor_used) tag = 0;
3297  tag < _jacobian_block_neighbor_used.size();
3298  tag++)
3299  if (_jacobian_block_neighbor_used[tag][i][j] && _sys.hasMatrix(tag))
3300  {
3302  ivar->dofIndices(),
3303  jvar->dofIndicesNeighbor(),
3304  ivar->scalingFactor(),
3305  tag);
3307  ivar->dofIndicesNeighbor(),
3308  jvar->dofIndices(),
3309  ivar->scalingFactor(),
3310  tag);
3312  ivar->dofIndicesNeighbor(),
3313  jvar->dofIndicesNeighbor(),
3314  ivar->scalingFactor(),
3315  tag);
3316  }
3317  }
3318 }
3319 
3320 void
3321 Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
3322  unsigned int ivar,
3323  unsigned int jvar,
3324  const DofMap & dof_map,
3325  std::vector<dof_id_type> & dof_indices)
3326 {
3327  DenseMatrix<Number> & ke = jacobianBlock(ivar, jvar);
3328 
3329  // stick it into the matrix
3330  std::vector<dof_id_type> di(dof_indices);
3331  dof_map.constrain_element_matrix(ke, di, false);
3332 
3333  Real scaling_factor = _sys.getVariable(_tid, ivar).scalingFactor();
3334  if (scaling_factor != 1.0)
3335  {
3336  _tmp_Ke = ke;
3337  _tmp_Ke *= scaling_factor;
3338  jacobian.add_matrix(_tmp_Ke, di);
3339  }
3340  else
3341  jacobian.add_matrix(ke, di);
3342 }
3343 
3344 void
3345 Assembly::addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
3346  unsigned int ivar,
3347  unsigned int jvar,
3348  const DofMap & dof_map,
3349  const std::vector<dof_id_type> & idof_indices,
3350  const std::vector<dof_id_type> & jdof_indices)
3351 {
3352  DenseMatrix<Number> & keg = jacobianBlockNonlocal(ivar, jvar);
3353 
3354  std::vector<dof_id_type> di(idof_indices);
3355  std::vector<dof_id_type> dg(jdof_indices);
3356  dof_map.constrain_element_matrix(keg, di, dg, false);
3357 
3358  Real scaling_factor = _sys.getVariable(_tid, ivar).scalingFactor();
3359  if (scaling_factor != 1.0)
3360  {
3361  _tmp_Ke = keg;
3362  _tmp_Ke *= scaling_factor;
3363  jacobian.add_matrix(_tmp_Ke, di, dg);
3364  }
3365  else
3366  jacobian.add_matrix(keg, di, dg);
3367 }
3368 
3369 void
3370 Assembly::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
3371  unsigned int ivar,
3372  unsigned int jvar,
3373  const DofMap & dof_map,
3374  std::vector<dof_id_type> & dof_indices,
3375  std::vector<dof_id_type> & neighbor_dof_indices)
3376 {
3377  DenseMatrix<Number> & kee = jacobianBlock(ivar, jvar);
3378  DenseMatrix<Number> & ken = jacobianBlockNeighbor(Moose::ElementNeighbor, ivar, jvar);
3379  DenseMatrix<Number> & kne = jacobianBlockNeighbor(Moose::NeighborElement, ivar, jvar);
3380  DenseMatrix<Number> & knn = jacobianBlockNeighbor(Moose::NeighborNeighbor, ivar, jvar);
3381 
3382  std::vector<dof_id_type> di(dof_indices);
3383  std::vector<dof_id_type> dn(neighbor_dof_indices);
3384  // stick it into the matrix
3385  dof_map.constrain_element_matrix(kee, di, false);
3386  dof_map.constrain_element_matrix(ken, di, dn, false);
3387  dof_map.constrain_element_matrix(kne, dn, di, false);
3388  dof_map.constrain_element_matrix(knn, dn, false);
3389 
3390  Real scaling_factor = _sys.getVariable(_tid, ivar).scalingFactor();
3391  if (scaling_factor != 1.0)
3392  {
3393  _tmp_Ke = ken;
3394  _tmp_Ke *= scaling_factor;
3395  jacobian.add_matrix(_tmp_Ke, di, dn);
3396 
3397  _tmp_Ke = kne;
3398  _tmp_Ke *= scaling_factor;
3399  jacobian.add_matrix(_tmp_Ke, dn, di);
3400 
3401  _tmp_Ke = knn;
3402  _tmp_Ke *= scaling_factor;
3403  jacobian.add_matrix(_tmp_Ke, dn);
3404  }
3405  else
3406  {
3407  jacobian.add_matrix(ken, di, dn);
3408  jacobian.add_matrix(kne, dn, di);
3409  jacobian.add_matrix(knn, dn);
3410  }
3411 }
3412 
3413 void
3415 {
3416  for (const auto & it : _cm_ss_entry)
3417  addJacobianCoupledVarPair(it.first, it.second);
3418 }
3419 
3420 void
3422 {
3423  const std::vector<MooseVariableFEBase *> & vars = _sys.getVariables(_tid);
3425  for (const auto & var_j : vars)
3426  addJacobianCoupledVarPair(&var_i, var_j);
3427 }
3428 
3429 void
3431  numeric_index_type j,
3432  Real value,
3433  TagID tag /* = 0*/)
3434 {
3435  _cached_jacobian_contribution_rows[tag].push_back(i);
3436  _cached_jacobian_contribution_cols[tag].push_back(j);
3437  _cached_jacobian_contribution_vals[tag].push_back(value);
3438 }
3439 
3440 void
3442  numeric_index_type j,
3443  Real value,
3444  const std::set<TagID> & tags)
3445 {
3446  for (auto tag : tags)
3447  if (_sys.hasMatrix(tag))
3448  cacheJacobianContribution(i, j, value, tag);
3449 }
3450 
3451 void
3453 {
3454  for (MooseIndex(_cached_jacobian_contribution_rows) tag = 0;
3456  tag++)
3457  if (_sys.hasMatrix(tag))
3458  {
3459  // First zero the rows (including the diagonals) to prepare for
3460  // setting the cached values.
3461  _sys.getMatrix(tag).zero_rows(_cached_jacobian_contribution_rows[tag], 0.0);
3462 
3463  // TODO: Use SparseMatrix::set_values() for efficiency
3464  for (MooseIndex(_cached_jacobian_contribution_vals) i = 0;
3465  i < _cached_jacobian_contribution_vals[tag].size();
3466  ++i)
3470  }
3471 
3473 }
3474 
3475 void
3477 {
3478  for (MooseIndex(_cached_jacobian_contribution_rows) tag = 0;
3480  tag++)
3481  if (_sys.hasMatrix(tag))
3482  _sys.getMatrix(tag).zero_rows(_cached_jacobian_contribution_rows[tag], 0.0);
3483 
3485 }
3486 
3487 void
3489 {
3490  for (MooseIndex(_cached_jacobian_contribution_rows) tag = 0;
3492  tag++)
3493  if (_sys.hasMatrix(tag))
3494  {
3495  // TODO: Use SparseMatrix::set_values() for efficiency
3496  for (MooseIndex(_cached_jacobian_contribution_vals[tag]) i = 0;
3497  i < _cached_jacobian_contribution_vals[tag].size();
3498  ++i)
3502  }
3503 
3505 }
3506 
3507 void
3509 {
3510  for (MooseIndex(_cached_jacobian_contribution_rows) tag = 0;
3512  tag++)
3513  {
3514  auto orig_size = _cached_jacobian_contribution_rows[tag].size();
3515 
3519 
3520  // It's possible (though massively unlikely) that clear() will
3521  // change the capacity of the vectors, so let's be paranoid and
3522  // explicitly reserve() the same amount of memory to avoid multiple
3523  // push_back() induced allocations. We reserve 20% more than the
3524  // original size that was cached to account for variations in the
3525  // number of BCs assigned to each thread (for when the Jacobian
3526  // contributions are computed threaded).
3527  _cached_jacobian_contribution_rows[tag].reserve(1.2 * orig_size);
3528  _cached_jacobian_contribution_cols[tag].reserve(1.2 * orig_size);
3529  _cached_jacobian_contribution_vals[tag].reserve(1.2 * orig_size);
3530  }
3531 }
3532 
3533 void
3535 {
3536  mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
3537 
3539  return;
3540 
3541  MooseArray<Real> xfem_weight_multipliers;
3542  if (_xfem->getXFEMWeights(xfem_weight_multipliers, elem, _current_qrule, _current_q_points))
3543  {
3544  mooseAssert(xfem_weight_multipliers.size() == _current_JxW.size(),
3545  "Size of weight multipliers in xfem doesn't match number of quadrature points");
3546  for (unsigned i = 0; i < xfem_weight_multipliers.size(); i++)
3547  _current_JxW[i] = _current_JxW[i] * xfem_weight_multipliers[i];
3548 
3549  xfem_weight_multipliers.release();
3550  }
3551 }
3552 
3553 void
3554 Assembly::modifyFaceWeightsDueToXFEM(const Elem * elem, unsigned int side)
3555 {
3556  mooseAssert(_xfem != nullptr, "This function should not be called if xfem is inactive");
3557 
3559  return;
3560 
3561  MooseArray<Real> xfem_face_weight_multipliers;
3562  if (_xfem->getXFEMFaceWeights(
3563  xfem_face_weight_multipliers, elem, _current_qrule_face, _current_q_points_face, side))
3564  {
3565  mooseAssert(xfem_face_weight_multipliers.size() == _current_JxW_face.size(),
3566  "Size of weight multipliers in xfem doesn't match number of quadrature points");
3567  for (unsigned i = 0; i < xfem_face_weight_multipliers.size(); i++)
3568  _current_JxW_face[i] = _current_JxW_face[i] * xfem_face_weight_multipliers[i];
3569 
3570  xfem_face_weight_multipliers.release();
3571  }
3572 }
3573 
3574 template <>
3576 Assembly::fePhi<VectorValue<Real>>(FEType type) const
3577 {
3578  buildVectorFE(type);
3579  return _vector_fe_shape_data[type]->_phi;
3580 }
3581 
3582 template <>
3584 Assembly::feGradPhi<VectorValue<Real>>(FEType type) const
3585 {
3586  buildVectorFE(type);
3587  return _vector_fe_shape_data[type]->_grad_phi;
3588 }
3589 
3590 template <>
3592 Assembly::feSecondPhi<VectorValue<Real>>(FEType type) const
3593 {
3594  _need_second_derivative[type] = true;
3595  buildVectorFE(type);
3596  return _vector_fe_shape_data[type]->_second_phi;
3597 }
3598 
3599 template <>
3601 Assembly::fePhiLower<VectorValue<Real>>(FEType type) const
3602 {
3603  buildVectorFE(type);
3604  return _vector_fe_shape_data_lower[type]->_phi;
3605 }
3606 
3607 template <>
3609 Assembly::feGradPhiLower<VectorValue<Real>>(FEType type) const
3610 {
3611  buildVectorFE(type);
3612  return _vector_fe_shape_data_lower[type]->_grad_phi;
3613 }
3614 
3615 template <>
3617 Assembly::fePhiFace<VectorValue<Real>>(FEType type) const
3618 {
3619  buildVectorFaceFE(type);
3620  return _vector_fe_shape_data_face[type]->_phi;
3621 }
3622 
3623 template <>
3625 Assembly::feGradPhiFace<VectorValue<Real>>(FEType type) const
3626 {
3627  buildVectorFaceFE(type);
3628  return _vector_fe_shape_data_face[type]->_grad_phi;
3629 }
3630 
3631 template <>
3633 Assembly::feSecondPhiFace<VectorValue<Real>>(FEType type) const
3634 {
3635  _need_second_derivative[type] = true;
3636  buildVectorFaceFE(type);
3637  return _vector_fe_shape_data_face[type]->_second_phi;
3638 }
3639 
3640 template <>
3642 Assembly::fePhiNeighbor<VectorValue<Real>>(FEType type) const
3643 {
3644  buildVectorNeighborFE(type);
3645  return _vector_fe_shape_data_neighbor[type]->_phi;
3646 }
3647 
3648 template <>
3650 Assembly::feGradPhiNeighbor<VectorValue<Real>>(FEType type) const
3651 {
3652  buildVectorNeighborFE(type);
3653  return _vector_fe_shape_data_neighbor[type]->_grad_phi;
3654 }
3655 
3656 template <>
3658 Assembly::feSecondPhiNeighbor<VectorValue<Real>>(FEType type) const
3659 {
3660  _need_second_derivative_neighbor[type] = true;
3661  buildVectorNeighborFE(type);
3662  return _vector_fe_shape_data_neighbor[type]->_second_phi;
3663 }
3664 
3665 template <>
3667 Assembly::fePhiFaceNeighbor<VectorValue<Real>>(FEType type) const
3668 {
3669  buildVectorFaceNeighborFE(type);
3670  return _vector_fe_shape_data_face_neighbor[type]->_phi;
3671 }
3672 
3673 template <>
3675 Assembly::feGradPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
3676 {
3677  buildVectorFaceNeighborFE(type);
3678  return _vector_fe_shape_data_face_neighbor[type]->_grad_phi;
3679 }
3680 
3681 template <>
3683 Assembly::feSecondPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
3684 {
3685  _need_second_derivative_neighbor[type] = true;
3686  buildVectorFaceNeighborFE(type);
3687  return _vector_fe_shape_data_face_neighbor[type]->_second_phi;
3688 }
3689 
3690 template <>
3692 Assembly::feCurlPhi<VectorValue<Real>>(FEType type) const
3693 {
3694  _need_curl[type] = true;
3695  buildVectorFE(type);
3696  return _vector_fe_shape_data[type]->_curl_phi;
3697 }
3698 
3699 template <>
3701 Assembly::feCurlPhiFace<VectorValue<Real>>(FEType type) const
3702 {
3703  _need_curl[type] = true;
3704  buildVectorFaceFE(type);
3705  return _vector_fe_shape_data_face[type]->_curl_phi;
3706 }
3707 
3708 template <>
3710 Assembly::feCurlPhiNeighbor<VectorValue<Real>>(FEType type) const
3711 {
3712  _need_curl[type] = true;
3713  buildVectorNeighborFE(type);
3714  return _vector_fe_shape_data_neighbor[type]->_curl_phi;
3715 }
3716 
3717 template <>
3719 Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const
3720 {
3721  _need_curl[type] = true;
3722  buildVectorFaceNeighborFE(type);
3723  return _vector_fe_shape_data_face_neighbor[type]->_curl_phi;
3724 }
3725 
3726 template <>
3728 Assembly::adGradPhi<Real, ComputeStage::JACOBIAN>(const MooseVariableFE<Real> & v) const
3729 {
3730  return _ad_grad_phi_data.at(v.feType());
3731 }
3732 
3733 template <>
3735 Assembly::adGradPhi<RealVectorValue, ComputeStage::JACOBIAN>(
3736  const MooseVariableFE<RealVectorValue> & v) const
3737 {
3738  return _ad_vector_grad_phi_data.at(v.feType());
3739 }
3740 
3741 template <>
3743 Assembly::adJxW<RESIDUAL>() const
3744 {
3745  return _current_JxW;
3746 }
3747 
3748 template <ComputeStage compute_stage>
3749 const MooseArray<ADReal> &
3751 {
3752  _calculate_curvatures = true;
3753  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
3754  (*_holder_fe_face_helper.at(dim))->get_curvatures();
3755  return _curvatures;
3756 }
3757 
3759 Assembly::adCurvatures<RESIDUAL>() const;
3760 
3761 template void Assembly::computeGradPhiAD<Real>(
3762  const Elem * elem,
3763  unsigned int n_qp,
3765  FEGenericBase<Real> * fe);
std::vector< std::pair< MooseVariableFEBase *, MooseVariableScalar * > > _cm_fs_entry
Entries in the coupling matrix for field variables vs scalar variables.
Definition: Assembly.h:1296
const Elem *const & elem() const
Return the current element.
Definition: Assembly.h:283
MooseArray< DualReal > _ad_JxW
Definition: Assembly.h:1668
unsigned int getAxisymmetricRadialCoord() const
Returns the desired radial direction for RZ coordinate transformation.
Definition: SubProblem.C:605
virtual bool hasMatrix(TagID tag) const
Check if the tagged matrix exists in the system.
Definition: SystemBase.C:805
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face
types of finite elements
Definition: Assembly.h:1389
std::map< unsigned int, std::map< FEType, const FEBase * > > _const_fe_lower
FE objects for lower dimensional elements.
Definition: Assembly.h:1441
bool _need_neighbor_elem_volume
true is apps need to compute neighbor element volume
Definition: Assembly.h:1495
ArbitraryQuadrature * _current_qrule_arbitrary
The current arbitrary quadrature rule used within the element interior.
Definition: Assembly.h:1361
std::map< FEType, FEBase * > _current_fe
The "volume" fe object that matches the current elem.
Definition: Assembly.h:1323
MooseArray< Real > _curvatures
Definition: Assembly.h:1683
VectorVariablePhiValue _phi
Definition: Assembly.h:1596
std::map< unsigned int, FEBase ** > _holder_fe_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:1351
SystemBase & _sys
Definition: Assembly.h:1282
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:1631
const FieldVariablePhiValue & phi() const
std::vector< std::vector< dof_id_type > > _cached_jacobian_rows
Row where the corresponding cached value should go.
Definition: Assembly.h:1638
unsigned int _max_cached_residuals
Definition: Assembly.h:1633
MooseArray< VectorValue< DualReal > > _ad_q_points_face
Definition: Assembly.h:1682
void addJacobianBlock(SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const DofMap &dof_map, std::vector< dof_id_type > &dof_indices)
Definition: Assembly.C:3321
std::vector< DualReal > _ad_detadz_map
Definition: Assembly.h:1675
MooseArray< Point > _current_physical_points
This will be filled up with the physical points passed into reinitAtPhysical() if it is called...
Definition: Assembly.h:1511
std::vector< DualReal > _ad_dxidy_map
Definition: Assembly.h:1671
const FieldVariablePhiSecond & secondPhiFace() const
const FieldVariablePhiGradient & gradPhiFace() const
std::map< unsigned int, ArbitraryQuadrature * > _holder_qrule_arbitrary_face
Holds arbitrary qrules for each dimension for faces.
Definition: Assembly.h:1380
MooseArray< VectorValue< DualReal > > _ad_normals
Definition: Assembly.h:1681
const QBase * _const_current_qrule
The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac k...
Definition: Assembly.h:1355
std::map< unsigned int, std::map< FEType, const FEVectorBase * > > _const_vector_fe_face_neighbor
Definition: Assembly.h:1432
const FieldVariablePhiCurl & curlPhi() const
void buildNeighborFE(FEType type) const
Build FEs for a neighbor with a type.
Definition: Assembly.C:277
void cacheResidualNodes(const DenseVector< Number > &res, const std::vector< dof_id_type > &dof_index, TagID tag=0)
Lets an external class cache residual at a set of nodes.
Definition: Assembly.C:2875
const VariablePhiSecond & secondPhi() const
Definition: Assembly.h:782
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe
Each dimension&#39;s actual vector fe objects indexed on type.
Definition: Assembly.h:1347
void buildVectorLowerDFE(FEType type) const
Build Vector FEs for a lower dimensional element with a type.
Definition: Assembly.C:339
std::vector< VectorValue< DualReal > > _ad_d2xyzdxi2_map
Definition: Assembly.h:1664
MooseArray< Real > _coord_neighbor
The current coordinate transformation coefficients.
Definition: Assembly.h:1458
virtual void haveADObjects(bool have_ad_objects)
Method for setting whether we have any ad objects.
Definition: SubProblem.h:618
void buildFE(FEType type) const
Build FEs with a type.
Definition: Assembly.C:234
QBase * _qrule_msm
A qrule object for working on mortar segement elements.
Definition: Assembly.h:1470
const std::vector< MooseVariableScalar * > & getScalarVariables(THREAD_ID tid)
Definition: SystemBase.h:614
QBase * _current_qrule
The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac k...
Definition: Assembly.h:1357
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kle
dlower/dslave (or dlower/delement)
Definition: Assembly.h:1535
std::map< unsigned int, FEBase ** > _holder_fe_face_neighbor_helper
Definition: Assembly.h:1436
virtual Moose::CoordinateSystemType getCoordSystem(SubdomainID sid)=0
bool hasVector(const std::string &tag_name) const
Check if the named vector exists in the system.
Definition: SystemBase.C:701
MooseArray< DualReal > _ad_JxW_face
Definition: Assembly.h:1680
void computeGradPhiAD(const Elem *elem, unsigned int n_qp, typename VariableTestGradientType< OutputType, ComputeStage::JACOBIAN >::type &grad_phi, FEGenericBase< OutputType > *fe)
Definition: Assembly.C:675
Assembly(SystemBase &sys, THREAD_ID tid)
Definition: Assembly.C:36
unsigned int TagID
Definition: MooseTypes.h:162
void init()
Deprecated init method.
const VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariable &) const
Definition: Assembly.h:809
PetscInt N
MooseArray< Real > _current_JxW_neighbor
The current transformed jacobian weights on a neighbor&#39;s face.
Definition: Assembly.h:1456
std::shared_ptr< XFEMInterface > _xfem
The XFEM controller.
Definition: Assembly.h:1320
const FieldVariablePhiSecond & secondPhi() const
virtual ~Assembly()
Definition: Assembly.C:120
void prepareJacobianBlock()
Sizes and zeroes the Jacobian blocks used for the current element.
Definition: Assembly.C:2302
void setResidualNeighbor(NumericVector< Number > &residual, TagID tag_id=0)
Definition: Assembly.C:2970
const FieldVariablePhiSecond & secondPhiFaceNeighbor() const
VectorValue< Real > RealVectorValue
Definition: Assembly.h:31
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:207
std::vector< DualReal > _ad_dzetady_map
Definition: Assembly.h:1677
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kel
dslave/dlower (or delement/dlower)
Definition: Assembly.h:1539
MooseArray< Real > _coord
The current coordinate transformation coefficients.
Definition: Assembly.h:1371
unsigned int number() const
Get variable number coming from libMesh.
void resize(unsigned int size)
Change the number of elements the array can store.
Definition: MooseArray.h:218
MooseArray< VectorValue< DualReal > > _ad_q_points
Definition: Assembly.h:1669
void reinitFE(const Elem *elem)
Just an internal helper function to reinit the volume FE objects.
Definition: Assembly.C:550
std::map< FEType, VectorFEShapeData * > _vector_fe_shape_data_neighbor
Definition: Assembly.h:1612
void prepareNonlocal()
Definition: Assembly.C:2341
size_t getMaxVarNDofsPerElem() const
Gets the maximum number of dofs used by any one variable on any one element.
Definition: SystemBase.h:442
std::unique_ptr< FEBase > _fe_msm
A FE object for working on mortar segement elements.
Definition: Assembly.h:1465
const VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariable &) const
Definition: Assembly.h:805
std::vector< DualReal > _ad_dzetadz_map
Definition: Assembly.h:1678
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:1546
std::map< unsigned int, ArbitraryQuadrature * > _holder_qrule_neighbor
Holds arbitrary qrules for each dimension.
Definition: Assembly.h:1454
DenseMatrix< Number > & jacobianBlock(unsigned int ivar, unsigned int jvar, TagID tag=0)
Definition: Assembly.C:2019
bool _current_elem_volume_computed
Boolean to indicate whether current element volumes has been computed.
Definition: Assembly.h:1503
std::map< unsigned int, std::map< FEType, const FEVectorBase * > > _const_vector_fe_lower
Vector FE objects for lower dimensional elements.
Definition: Assembly.h:1445
bool usesPhiNeighbor() const
Whether or not this variable is actually using the shape function value.
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kll
dlower/dlower
Definition: Assembly.h:1533
VariablePhiValue _phi
Definition: Assembly.h:1587
const FieldVariablePhiValue & phiFace() const
Real _current_neighbor_volume
Volume of the current neighbor.
Definition: Assembly.h:1497
std::map< FEType, FEBase * > _current_fe_face
The "face" fe object that matches the current elem.
Definition: Assembly.h:1325
void addResidualNeighbor(NumericVector< Number > &residual, TagID tag_id=0)
Definition: Assembly.C:2772
const FieldVariablePhiSecond & secondPhiNeighbor() const
void computeAffineMapAD(const Elem *elem, const std::vector< Real > &qw, unsigned int n_qp, FEBase *fe)
Definition: Assembly.C:772
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kln
dlower/dmaster (or dlower/dneighbor)
Definition: Assembly.h:1537
const Elem * _current_neighbor_elem
The current neighbor "element".
Definition: Assembly.h:1487
std::vector< std::vector< DenseVector< Number > > > _sub_Rn
residual contributions for each variable from the neighbor
Definition: Assembly.h:1516
const VariablePhiValue & phi() const
Definition: Assembly.h:772
MooseMesh & _mesh
Definition: Assembly.h:1315
void cacheJacobianCoupledVarPair(MooseVariableBase *ivar, MooseVariableBase *jvar)
Caches element matrix for ivar rows and jvar columns.
Definition: Assembly.C:3239
void setPoints(const std::vector< Point > &points)
const VariablePhiGradient & gradPhi() const
Definition: Assembly.h:780
MooseArray< Real > _current_JxW_face
The current transformed jacobian weights on a face.
Definition: Assembly.h:1409
void modifyWeightsDueToXFEM(const Elem *elem)
Update the integration weights for XFEM partial elements.
Definition: Assembly.C:3534
DenseMatrix< Number > _tmp_Ke
auxiliary matrix for scaling jacobians (optimization to avoid expensive construction/destruction) ...
Definition: Assembly.h:1544
const VariablePhiSecond & secondPhiFace(const MooseVariable &) const
Definition: Assembly.h:789
const Elem * _current_elem
The current "element" we are currently on.
Definition: Assembly.h:1475
const std::vector< MooseVariableFEBase * > & getVariables(THREAD_ID tid)
Definition: SystemBase.h:609
THREAD_ID _tid
Thread number (id)
Definition: Assembly.h:1313
void reinitFEFaceNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:1318
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.
Definition: Assembly.C:1976
const MooseArray< ADReal > & adCurvatures() const
Definition: Assembly.C:3750
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_nonlocal_used
Definition: Assembly.h:1305
void prepareNeighbor()
Definition: Assembly.C:2419
OutputTools< Real >::VariablePhiValue VariablePhiValue
Definition: MooseTypes.h:201
std::map< FEType, typename VariableTestGradientType< Real, ComputeStage::JACOBIAN >::type > _ad_grad_phi_data
Definition: Assembly.h:1617
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_face_neighbor
Definition: Assembly.h:1427
std::map< FEType, bool > _need_second_derivative
Definition: Assembly.h:1692
std::vector< std::pair< MooseVariableFEBase *, MooseVariableFEBase * > > _cm_ff_entry
Entries in the coupling matrix for field variables.
Definition: Assembly.h:1294
Real _current_elem_volume
Volume of the current element.
Definition: Assembly.h:1479
std::map< FEType, FEBase * > _current_fe_face_neighbor
The "neighbor face" fe object that matches the current elem.
Definition: Assembly.h:1329
std::map< FEType, FEVectorBase * > _current_vector_fe_face
The "face" vector fe object that matches the current elem.
Definition: Assembly.h:1334
void shallowCopy(const MooseArray &rhs)
Doesn&#39;t actually make a copy of the data.
Definition: MooseArray.h:294
std::vector< DualReal > _ad_dzetadx_map
Definition: Assembly.h:1676
std::map< FEType, FEShapeData * > _fe_shape_data_face_neighbor
Definition: Assembly.h:1606
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Knn
jacobian contributions from the neighbor <Tag, ivar, jvar>
Definition: Assembly.h:1531
Base class for a system (of equations)
Definition: SystemBase.h:92
void setVolumeQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for volume integration.
Definition: Assembly.C:514
void zeroCachedJacobianContributions()
Zero out previously-cached Jacobian rows.
Definition: Assembly.C:3476
unsigned int _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:1645
void addJacobianNonlocal()
Definition: Assembly.C:3137
std::vector< std::vector< Real > > _cached_jacobian_values
Values cached by calling cacheJacobian()
Definition: Assembly.h:1636
std::map< FEType, FEVectorBase * > _current_vector_fe
The "volume" vector fe object that matches the current elem.
Definition: Assembly.h:1332
VariablePhiGradient _grad_phi
Definition: Assembly.h:1588
void prepareScalar()
Definition: Assembly.C:2550
void copyNeighborShapes(MooseVariableFE< T > &v)
Definition: Assembly.C:2660
void modifyFaceWeightsDueToXFEM(const Elem *elem, unsigned int side=0)
Update the face integration weights for XFEM partial elements.
Definition: Assembly.C:3554
const CouplingMatrix * _cm
Coupling matrices.
Definition: Assembly.h:1288
std::map< FEType, bool > _need_second_derivative_neighbor
Definition: Assembly.h:1693
void createQRules(QuadratureType type, Order order, Order volume_order, Order face_order)
Creates the volume, face and arbitrary qrules based on the orders passed in.
Definition: Assembly.C:487
OutputTools< Real >::VariablePhiSecond VariablePhiSecond
Definition: MooseTypes.h:203
DenseVector< Number > _tmp_Re
auxiliary vector for scaling residuals (optimization to avoid expensive construction/destruction) ...
Definition: Assembly.h:1520
unsigned int _mesh_dimension
Definition: Assembly.h:1317
void reinit(const Elem *elem)
Reinitialize objects (JxW, q_points, ...) for an elements.
Definition: Assembly.C:1567
QBase * _current_qrule_volume
The current volumetric quadrature for the element.
Definition: Assembly.h:1359
const DofMap & _dof_map
DOF map.
Definition: Assembly.h:1311
std::vector< std::vector< Real > > _cached_jacobian_contribution_vals
Storage for cached Jacobian entries.
Definition: Assembly.h:1656
std::vector< DualReal > _ad_dxidx_map
Definition: Assembly.h:1670
unsigned int _max_cached_jacobians
Definition: Assembly.h:1642
const VariablePhiValue & phiFaceNeighbor(const MooseVariable &) const
Definition: Assembly.h:801
void cacheResidual()
Takes the values that are currently in _sub_Re and appends them to the cached values.
Definition: Assembly.C:2809
std::map< FEType, VectorFEShapeData * > _vector_fe_shape_data_face_neighbor
Definition: Assembly.h:1613
bool usesGradPhiNeighbor() const
Whether or not this variable is actually using the shape function gradient.
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Ken
jacobian contributions from the element and neighbor <Tag, ivar, jvar>
Definition: Assembly.h:1527
const QBase * _const_qrule_msm
A pointer to const qrule_msm.
Definition: Assembly.h:1472
virtual unsigned int nVariables() const
Get the number of variables in this system.
Definition: SystemBase.C:692
const CouplingMatrix & _nonlocal_cm
Definition: Assembly.h:1289
MooseVariableFEBase & getVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:105
std::vector< Point > _temp_reference_points
Temporary work data for reinitAtPhysical()
Definition: Assembly.h:1651
std::vector< std::pair< MooseVariableScalar *, MooseVariableScalar * > > _cm_ss_entry
Entries in the coupling matrix for scalar variables.
Definition: Assembly.h:1300
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kne
jacobian contributions from the neighbor and element <Tag, ivar, jvar>
Definition: Assembly.h:1529
unsigned int _current_neighbor_side
The current side of the selected neighboring element (valid only when working with sides) ...
Definition: Assembly.h:1491
void prepareLowerD()
Prepare the Jacobians and residuals for a lower dimensional element.
Definition: Assembly.C:2459
void buildFaceNeighborFE(FEType type) const
Build FEs for a neighbor face with a type.
Definition: Assembly.C:297
void computeCurrentFaceVolume()
Definition: Assembly.C:1527
std::vector< unsigned > _displacements
Definition: Assembly.h:1686
std::vector< std::vector< DenseVector< Number > > > _sub_Rl
residual contributions for each variable from the lower dimensional element
Definition: Assembly.h:1518
void computeADFace(const Elem *elem, unsigned int side)
compute AD things on an element face
Definition: Assembly.C:1782
const FieldVariablePhiGradient & gradPhiNeighbor() const
DenseMatrix< Number > & jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar, TagID tag=0)
Definition: Assembly.C:2026
std::map< unsigned int, QBase * > _holder_qrule_face
Holds face qrules for each dimension.
Definition: Assembly.h:1413
void cacheJacobianBlockNonlocal(DenseMatrix< Number > &jac_block, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, Real scaling_factor, TagID tag=0)
Definition: Assembly.C:3032
void prepareOffDiagScalar()
Definition: Assembly.C:2578
std::map< unsigned int, std::map< FEType, const FEBase * > > _const_fe_neighbor
Definition: Assembly.h:1428
VectorVariablePhiSecond _second_phi
Definition: Assembly.h:1598
std::map< unsigned int, std::map< FEType, const FEBase * > > _const_fe_face_neighbor
Definition: Assembly.h:1429
unsigned int size() const
The number of elements that can currently be stored in the array.
Definition: MooseArray.h:259
std::map< unsigned int, std::map< FEType, const FEBase * > > _const_fe_face
types of finite elements
Definition: Assembly.h:1391
Implements a fake quadrature rule where you can specify the locations (in the reference domain) of th...
void addCachedResidual(NumericVector< Number > &residual, TagID tag_id)
Adds the values that have been cached by calling cacheResidual(), cacheResidualNeighbor(), and/or cacheResidualLower() to the residual.
Definition: Assembly.C:2905
void copyShapes(MooseVariableFE< T > &v)
Definition: Assembly.C:2606
void addJacobian()
Definition: Assembly.C:3124
void setCoordinateTransformation(const QBase *qrule, const ADPoint &q_points, MooseArray< ADReal > &coord)
Definition: Assembly.C:1473
void setResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, Real scaling_factor)
Definition: Assembly.C:2939
const FieldVariablePhiValue & phiFaceNeighbor() const
QBase * _current_qrule_neighbor
quadrature rule used on neighbors
Definition: Assembly.h:1450
std::map< FEType, FEShapeData * > _fe_shape_data_neighbor
Definition: Assembly.h:1605
std::vector< DualReal > _ad_dxidz_map
Definition: Assembly.h:1672
std::map< unsigned int, std::map< FEType, const FEBase * > > _const_fe
Each dimension&#39;s actual fe objects indexed on type.
Definition: Assembly.h:1345
std::vector< std::vector< numeric_index_type > > _cached_jacobian_contribution_cols
Definition: Assembly.h:1658
VectorVariablePhiGradient _grad_phi
Definition: Assembly.h:1597
ArbitraryQuadrature * _current_qrule_arbitrary_face
The current arbitrary quadrature rule used on the element face.
Definition: Assembly.h:1363
const std::vector< Real > * _JxW_msm
A JxW for working on mortar segement elements.
Definition: Assembly.h:1463
void prepareResidual()
Sizes and zeroes the residual for the current element.
Definition: Assembly.C:2322
void prepareVariableNonlocal(MooseVariableFEBase *var)
Definition: Assembly.C:2393
void addJacobianOffDiagScalar(unsigned int ivar)
Definition: Assembly.C:3421
VectorVariablePhiCurl _curl_phi
Definition: Assembly.h:1599
std::map< FEType, FEBase * > _current_fe_neighbor
The "neighbor" fe object that matches the current elem.
Definition: Assembly.h:1327
void cacheJacobianContribution(numeric_index_type i, numeric_index_type j, Real value, TagID tag=0)
Caches the Jacobian entry &#39;value&#39;, to eventually be added/set in the (i,j) location of the matrix...
Definition: Assembly.C:3430
SubdomainID _current_neighbor_subdomain_id
The current neighbor subdomain ID.
Definition: Assembly.h:1489
std::map< unsigned int, FEBase ** > _holder_fe_neighbor_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:1435
std::vector< std::vector< DenseVector< Number > > > _sub_Re
residual contributions for each variable from the element
Definition: Assembly.h:1514
void reinitNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:1413
const FieldVariablePhiValue & phiNeighbor() const
DenseMatrix< Number > & jacobianBlockLower(Moose::ConstraintJacobianType type, unsigned int ivar, unsigned int jvar, TagID tag=0)
Returns the jacobian block for the given mortar Jacobian type.
Definition: Assembly.C:2072
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_used
Flag that indicates if the jacobian block was used.
Definition: Assembly.h:1304
void prepare()
Definition: Assembly.C:2334
DofMap & dof_map
void setFaceQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for face integration.
Definition: Assembly.C:528
void reinitElemAndNeighbor(const Elem *elem, unsigned int side, const Elem *neighbor, unsigned int neighbor_side)
Reinitialize an element and its neighbor along a particular side.
Definition: Assembly.C:1679
std::vector< DualReal > _ad_jac
Definition: Assembly.h:1667
void addJacobianLower()
Add LowerLower, LowerSlave (LowerElement), LowerMaster (LowerNeighbor), SlaveLower (ElementLower)...
Definition: Assembly.C:3193
SubProblem & _subproblem
Definition: Assembly.h:1283
void buildFaceFE(FEType type) const
Build FEs for a face with a type.
Definition: Assembly.C:257
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_neighbor
types of finite elements
Definition: Assembly.h:1424
bool _calculate_xyz
Definition: Assembly.h:1688
virtual const std::vector< dof_id_type > & dofIndices() const
Get local DoF indices.
bool _calculate_curvatures
Definition: Assembly.h:1690
std::map< FEType, FEVectorBase * > _current_vector_fe_face_neighbor
The "neighbor face" vector fe object that matches the current elem.
Definition: Assembly.h:1338
MooseArray< Real > _current_JxW
The current list of transformed jacobian weights.
Definition: Assembly.h:1367
virtual bool checkNonlocalCouplingRequirement()
Definition: SubProblem.h:69
void prepareVariable(MooseVariableFEBase *var)
Used for preparing the dense residual and jacobian blocks for one particular variable.
Definition: Assembly.C:2364
std::vector< dof_id_type > _temp_dof_indices
Temporary work vector to keep from reallocating it.
Definition: Assembly.h:1648
OutputTools< Real >::VariablePhiCurl VariablePhiCurl
Definition: MooseTypes.h:204
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kee
jacobian contributions <Tag, ivar, jvar>
Definition: Assembly.h:1523
void cacheResidualNeighbor()
Takes the values that are currently in _sub_Rn and appends them to the cached values.
Definition: Assembly.C:2839
std::map< unsigned int, ArbitraryQuadrature * > _holder_qrule_arbitrary
Holds arbitrary qrules for each dimension.
Definition: Assembly.h:1378
void mooseDeprecated(Args &&... args)
Emit a deprecated code/feature message with the given stringified, concatenated args.
Definition: MooseError.h:236
virtual bool vectorTagExists(TagID tag)
Check to see if a particular Tag exists.
Definition: SubProblem.h:112
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_face
types of vector finite elements
Definition: Assembly.h:1393
std::map< FEType, typename VariableTestGradientType< RealVectorValue, ComputeStage::JACOBIAN >::type > _ad_vector_grad_phi_data_face
Definition: Assembly.h:1625
std::vector< T > stdVector() const
Extremely inefficient way to produce a std::vector from a MooseArray!
Definition: MooseArray.h:340
const VariablePhiSecond & secondPhiNeighbor(const MooseVariable &) const
Definition: Assembly.h:796
const Elem * _current_lower_d_elem
The current lower dimensional element.
Definition: Assembly.h:1508
std::map< FEType, typename VariableTestGradientType< RealVectorValue, ComputeStage::JACOBIAN >::type > _ad_vector_grad_phi_data
Definition: Assembly.h:1620
std::vector< VectorValue< DualReal > > _ad_dxyzdeta_map
Definition: Assembly.h:1662
const VariablePhiValue & phiFace() const
Definition: Assembly.h:785
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 ...
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_lower
FE objects for lower dimensional elements.
Definition: Assembly.h:1439
bool hasSecondOrderElements()
check if the mesh has SECOND order elements
Definition: MooseMesh.C:2708
const VectorVariablePhiCurl & curlPhi(const VectorMooseVariable &) const
Definition: Assembly.h:823
FEGenericBase< Real > FEBase
Definition: Assembly.h:34
std::map< unsigned int, FEBase ** > _holder_fe_face_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:1397
virtual unsigned int numVectorTags() const
The total number of tags.
Definition: SubProblem.h:122
std::map< unsigned int, std::map< FEType, FEBase * > > _fe
Each dimension&#39;s actual fe objects indexed on type.
Definition: Assembly.h:1343
const MooseArray< Real > & JxW() const
Returns the reference to the transformed jacobian weights.
Definition: Assembly.h:192
void reinitMortarElem(const Elem *elem)
reinitialize a mortar segment mesh element in order to get a proper JxW
Definition: Assembly.C:1967
Real _current_side_volume
Volume of the current side element.
Definition: Assembly.h:1485
CoordinateSystemType
Definition: MooseTypes.h:556
std::vector< std::vector< Real > > _cached_residual_values
Values cached by calling cacheResidual() (the first vector is for TIME vs NONTIME) ...
Definition: Assembly.h:1628
void addJacobianCoupledVarPair(MooseVariableBase *ivar, MooseVariableBase *jvar)
Adds element matrix for ivar rows and jvar columns.
Definition: Assembly.C:3110
const bool _displaced
Definition: Assembly.h:1285
void addJacobianNeighbor()
Add ElementNeighbor, NeighborElement, and NeighborNeighbor portions of the Jacobian for compute objec...
Definition: Assembly.C:3158
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, TagID tag=0)
Definition: Assembly.C:3003
std::map< FEType, VectorFEShapeData * > _vector_fe_shape_data_lower
Definition: Assembly.h:1614
void buildVectorFE(FEType type) const
Build Vector FEs with a type.
Definition: Assembly.C:361
MooseArray< DualReal > _ad_curvatures
Definition: Assembly.h:1684
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_lower
Vector FE objects for lower dimensional elements.
Definition: Assembly.h:1443
std::map< unsigned int, std::map< FEType, const FEVectorBase * > > _const_vector_fe
Each dimension&#39;s actual vector fe objects indexed on type.
Definition: Assembly.h:1349
std::vector< VectorValue< DualReal > > _ad_d2xyzdxideta_map
Definition: Assembly.h:1665
const Node * _current_node
The current node we are working with.
Definition: Assembly.h:1499
bool _calculate_face_xyz
Definition: Assembly.h:1689
std::map< FEType, typename VariableTestGradientType< Real, ComputeStage::JACOBIAN >::type > _ad_grad_phi_data_face
Definition: Assembly.h:1622
std::vector< std::pair< MooseVariableScalar *, MooseVariableFEBase * > > _cm_sf_entry
Entries in the coupling matrix for scalar variables vs field variables.
Definition: Assembly.h:1298
virtual unsigned int numMatrixTags() const
The total number of tags.
Definition: SubProblem.h:157
DGJacobianType
Definition: MooseTypes.h:515
void addCachedJacobian()
Definition: Assembly.C:3070
OutputTools< Real >::VariablePhiGradient VariablePhiGradient
Definition: MooseTypes.h:202
MooseArray< Point > _current_normals
The current Normal vectors at the quadrature points.
Definition: Assembly.h:1411
MatType type
void reinitLowerDElemRef(const Elem *elem, const std::vector< Point > *const pts, 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:1935
void buildVectorFaceFE(FEType type) const
Build Vector FEs for a face with a type.
Definition: Assembly.C:393
void cacheJacobian()
Takes the values that are currently in _sub_Kee and appends them to the cached values.
Definition: Assembly.C:3253
forward declarations
Definition: MooseArray.h:16
std::vector< std::pair< MooseVariableFEBase *, MooseVariableFEBase * > > _cm_nonlocal_entry
Entries in the coupling matrix for field variables for nonlocal calculations.
Definition: Assembly.h:1302
const QBase * _const_current_qrule_face
quadrature rule used on faces
Definition: Assembly.h:1401
void setNeighborQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for neighbor integration.
Definition: Assembly.C:539
std::map< FEType, VectorFEShapeData * > _vector_fe_shape_data
Shape function values, gradients, second derivatives for each vector FE type.
Definition: Assembly.h:1610
std::vector< VectorValue< DualReal > > _ad_dxyzdzeta_map
Definition: Assembly.h:1663
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:1866
PetscInt M
void setResidual(NumericVector< Number > &residual, TagID tag_id=0)
Definition: Assembly.C:2961
MooseVariableFE< T > & getFieldVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:136
void release()
Manually deallocates the data pointer.
Definition: MooseArray.h:61
virtual SparseMatrix< Number > & getMatrix(TagID tag)
Get a raw SparseMatrix.
Definition: SystemBase.C:811
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Keg
Definition: Assembly.h:1524
void resizeMappingObjects(unsigned int n_qp, unsigned int dim)
Definition: Assembly.C:745
const std::vector< dof_id_type > & allDofIndices() const
Get all global dofindices for the variable.
std::map< unsigned int, std::map< FEType, const FEVectorBase * > > _const_vector_fe_neighbor
Definition: Assembly.h:1430
void buildVectorNeighborFE(FEType type) const
Build Vector FEs for a neighbor with a type.
Definition: Assembly.C:422
void reinitFENeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:1366
void initNonlocalCoupling()
Create pair of variables requiring nonlocal jacobian contributions.
Definition: Assembly.C:2283
void cacheJacobianNeighbor()
Takes the values that are currently in the neighbor Dense Matrices and appends them to the cached val...
Definition: Assembly.C:3287
ConstraintJacobianType
Definition: MooseTypes.h:543
void addResidual(NumericVector< Number > &residual, TagID tag_id=0)
Definition: Assembly.C:2755
const bool & _computing_jacobian
Definition: Assembly.h:1291
bool computingCurl() const
Whether or not this variable is computing the curl.
const Elem *const & neighbor() const
Return the neighbor element.
Definition: Assembly.h:329
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_neighbor
Definition: Assembly.h:1426
const FieldVariablePhiGradient & gradPhi() const
void addResidualScalar(TagID tag_id)
Definition: Assembly.C:2789
Class for scalar variables (they are different).
void computeFaceMap(unsigned dim, const std::vector< Real > &qw, const Elem *side)
Definition: Assembly.C:1107
void computeSinglePointMapAD(const Elem *elem, const std::vector< Real > &qw, unsigned p, FEBase *fe)
Definition: Assembly.C:829
void addJacobianScalar()
Definition: Assembly.C:3414
const QBase * _const_current_qrule_neighbor
quadrature rule used on neighbors
Definition: Assembly.h:1448
std::vector< DualReal > _ad_detadx_map
Definition: Assembly.h:1673
void buildLowerDFE(FEType type) const
Build FEs for a lower dimensional element with a type.
Definition: Assembly.C:317
unsigned int _current_side
The current side of the selected element (valid only when working with sides)
Definition: Assembly.h:1481
VariablePhiSecond _second_phi
Definition: Assembly.h:1589
const Elem * _current_neighbor_side_elem
The current side element of the ncurrent neighbor element.
Definition: Assembly.h:1493
std::map< unsigned int, QBase * > _holder_qrule_volume
Holds volume qrules for each dimension.
Definition: Assembly.h:1376
void addCachedJacobianContributions()
Adds previously-cached Jacobian values via SparseMatrix::add() calls.
Definition: Assembly.C:3488
DenseMatrix< Number > & jacobianBlockNeighbor(Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, TagID tag=0)
Definition: Assembly.C:2033
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:1699
const Real & neighborVolume() const
Returns the reference to the current neighbor volume.
Definition: Assembly.C:480
bool _current_side_volume_computed
Boolean to indicate whether current element side volumes has been computed.
Definition: Assembly.h:1505
virtual bool isScalarVariable(unsigned int var_name) const
Definition: SystemBase.C:686
void cacheResidualContribution(dof_id_type dof, Real value, TagID tag_id)
Cache individual residual contributions.
Definition: Assembly.C:2825
std::vector< VectorValue< DualReal > > _ad_d2xyzdeta2_map
Definition: Assembly.h:1666
const VariablePhiValue & phiNeighbor(const MooseVariable &) const
Definition: Assembly.h:791
void prepareBlock(unsigned int ivar, unsigned jvar, const std::vector< dof_id_type > &dof_indices)
Definition: Assembly.C:2515
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face_neighbor
Definition: Assembly.h:1425
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:1309
Definition: Moose.h:111
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, Real scaling_factor)
Definition: Assembly.C:2719
const unsigned int & side() const
Returns the current side.
Definition: Assembly.h:305
void computeCurrentElemVolume()
Definition: Assembly.C:1510
std::map< FEType, FEShapeData * > _fe_shape_data
Shape function values, gradients, second derivatives for each FE type.
Definition: Assembly.h:1603
std::map< FEType, FEShapeData * > _fe_shape_data_lower
Definition: Assembly.h:1607
bool usesSecondPhiNeighbor() const
Whether or not this variable is actually using the shape function second derivative on a neighbor...
virtual const std::vector< dof_id_type > & dofIndicesNeighbor() const =0
Get neighbor DOF indices for currently selected element.
MooseArray< Point > _current_q_points
The current list of quadrature points.
Definition: Assembly.h:1365
Moose::CoordinateSystemType _coord_type
The coordinate system.
Definition: Assembly.h:1369
void cacheResidualLower()
Takes the values that are currently in _sub_Rl and appends them to the cached values.
Definition: Assembly.C:2857
void buildVectorFaceNeighborFE(FEType type) const
Build Vector FEs for a neighbor face with a type.
Definition: Assembly.C:451
bool computingSecond() const
Whether or not this variable is computing any second derivatives.
const Node *const & node() const
Returns the reference to the node.
Definition: Assembly.h:381
std::vector< std::vector< numeric_index_type > > _cached_jacobian_contribution_rows
Definition: Assembly.h:1657
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:1307
std::map< FEType, FEShapeData * > _fe_shape_data_face
Definition: Assembly.h:1604
std::vector< DualReal > _ad_detady_map
Definition: Assembly.h:1674
std::vector< std::vector< dof_id_type > > _cached_jacobian_cols
Column where the corresponding cached value should go.
Definition: Assembly.h:1640
void clearCachedJacobianContributions()
Clear any currently cached jacobian contributions.
Definition: Assembly.C:3508
const VariablePhiGradient & gradPhiNeighbor(const MooseVariable &) const
Definition: Assembly.h:792
MooseArray< Point > _current_q_points_face
The current quadrature points on a face.
Definition: Assembly.h:1407
std::map< FEType, VectorFEShapeData * > _vector_fe_shape_data_face
Definition: Assembly.h:1611
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:2534
virtual NumericVector< Number > & getVector(const std::string &name)
Get a raw NumericVector.
Definition: SystemBase.C:751
const Elem * _current_side_elem
The current "element" making up the side we are currently on.
Definition: Assembly.h:1483
MooseArray< DualReal > _ad_coord
The AD version of the current coordinate transformation coefficients.
Definition: Assembly.h:1373
std::map< FEType, bool > _need_curl
Definition: Assembly.h:1694
void cacheJacobianNonlocal()
Takes the values that are currently in _sub_Keg and appends them to the cached values.
Definition: Assembly.C:3266
MooseArray< Point > _current_q_points_face_neighbor
The current quadrature points on the neighbor face.
Definition: Assembly.h:1452
FEGenericBase< VectorValue< Real > > FEVectorBase
Definition: Assembly.h:36
void copyFaceShapes(MooseVariableFE< T > &v)
Definition: Assembly.C:2633
virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a scalar variable with specified number.
Definition: SystemBase.C:149
void addCachedResiduals()
Definition: Assembly.C:2890
VectorVariablePhiCurl _vector_curl_phi_face
Definition: Assembly.h:1572
std::map< unsigned int, std::map< FEType, const FEVectorBase * > > _const_vector_fe_face
types of vector finite elements
Definition: Assembly.h:1395
SubdomainID _current_subdomain_id
The current subdomain ID.
Definition: Assembly.h:1477
void addResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, Real scaling_factor)
Definition: Assembly.C:2695
QBase * _current_qrule_face
quadrature rule used on faces
Definition: Assembly.h:1403
std::map< FEType, FEVectorBase * > _current_vector_fe_neighbor
The "neighbor" vector fe object that matches the current elem.
Definition: Assembly.h:1336
unsigned int THREAD_ID
Definition: MooseTypes.h:161
const Node * _current_neighbor_node
The current neighboring node we are working with.
Definition: Assembly.h:1501
const VariablePhiGradient & gradPhiFace() const
Definition: Assembly.h:787
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Knl
dmaster/dlower (or dneighbor/dlower)
Definition: Assembly.h:1541
void addJacobianBlockNonlocal(SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const DofMap &dof_map, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices)
Definition: Assembly.C:3345
void setCachedJacobianContributions()
Sets previously-cached Jacobian values via SparseMatrix::set() calls.
Definition: Assembly.C:3452
std::vector< VectorValue< DualReal > > _ad_dxyzdxi_map
AD quantities.
Definition: Assembly.h:1661
void reinitFEFace(const Elem *elem, unsigned int side)
Just an internal helper function to reinit the face FE objects.
Definition: Assembly.C:1046
const FieldVariablePhiGradient & gradPhiFaceNeighbor() const
void scalingFactor(Real factor)
Set the scaling factor for this variable.