www.mooseframework.org
Assembly.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "Assembly.h"
16 
17 // MOOSE includes
18 #include "SubProblem.h"
19 #include "ArbitraryQuadrature.h"
20 #include "SystemBase.h"
21 #include "MooseTypes.h"
22 #include "MooseMesh.h"
23 #include "MooseVariable.h"
24 #include "MooseVariableScalar.h"
25 #include "XFEMInterface.h"
26 
27 // libMesh
28 #include "libmesh/coupling_matrix.h"
29 #include "libmesh/dof_map.h"
30 #include "libmesh/elem.h"
31 #include "libmesh/equation_systems.h"
32 #include "libmesh/fe_interface.h"
33 #include "libmesh/node.h"
34 #include "libmesh/quadrature_gauss.h"
35 #include "libmesh/sparse_matrix.h"
36 #include "libmesh/tensor_value.h"
37 #include "libmesh/vector_value.h"
38 
40  : _sys(sys),
41  _nonlocal_cm(_sys.subproblem().nonlocalCouplingMatrix()),
42  _dof_map(_sys.dofMap()),
43  _tid(tid),
44  _mesh(sys.mesh()),
45  _mesh_dimension(_mesh.dimension()),
46  _current_qrule(NULL),
47  _current_qrule_volume(NULL),
48  _current_qrule_arbitrary(NULL),
49  _coord_type(Moose::COORD_XYZ),
50  _current_qrule_face(NULL),
51  _current_qface_arbitrary(NULL),
52  _current_qrule_neighbor(NULL),
53 
54  _current_elem(NULL),
55  _current_elem_volume(0),
56  _current_side(0),
57  _current_side_elem(NULL),
58  _current_side_volume(0),
59  _current_neighbor_elem(NULL),
60  _current_neighbor_side(0),
61  _current_neighbor_side_elem(NULL),
62  _need_neighbor_elem_volume(false),
63  _current_neighbor_volume(0),
64  _current_node(NULL),
65  _current_neighbor_node(NULL),
66  _current_elem_volume_computed(false),
67  _current_side_volume_computed(false),
68 
69  _should_use_fe_cache(false),
70  _currently_fe_caching(true),
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 {
79  // Build fe's for the helpers
80  buildFE(FEType(FIRST, LAGRANGE));
81  buildFaceFE(FEType(FIRST, LAGRANGE));
82  buildNeighborFE(FEType(FIRST, LAGRANGE));
83  buildFaceNeighborFE(FEType(FIRST, LAGRANGE));
84 
85  // Build an FE helper object for this type for each dimension up to the dimension of the current
86  // mesh
87  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
88  {
89  _holder_fe_helper[dim] = &_fe[dim][FEType(FIRST, LAGRANGE)];
90  (*_holder_fe_helper[dim])->get_phi();
91  (*_holder_fe_helper[dim])->get_dphi();
92  (*_holder_fe_helper[dim])->get_xyz();
93  (*_holder_fe_helper[dim])->get_JxW();
94 
95  _holder_fe_face_helper[dim] = &_fe_face[dim][FEType(FIRST, LAGRANGE)];
96  (*_holder_fe_face_helper[dim])->get_phi();
97  (*_holder_fe_face_helper[dim])->get_dphi();
98  (*_holder_fe_face_helper[dim])->get_xyz();
99  (*_holder_fe_face_helper[dim])->get_JxW();
100  (*_holder_fe_face_helper[dim])->get_normals();
101 
102  _holder_fe_face_neighbor_helper[dim] = &_fe_face_neighbor[dim][FEType(FIRST, LAGRANGE)];
103  (*_holder_fe_face_neighbor_helper[dim])->get_xyz();
104  (*_holder_fe_face_neighbor_helper[dim])->get_JxW();
105  (*_holder_fe_face_neighbor_helper[dim])->get_normals();
106 
107  _holder_fe_neighbor_helper[dim] = &_fe_neighbor[dim][FEType(FIRST, LAGRANGE)];
108  (*_holder_fe_neighbor_helper[dim])->get_xyz();
109  (*_holder_fe_neighbor_helper[dim])->get_JxW();
110  }
111 }
112 
114 {
115  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
116  for (auto & it : _fe[dim])
117  delete it.second;
118 
119  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
120  for (auto & it : _fe_face[dim])
121  delete it.second;
122 
123  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
124  for (auto & it : _fe_neighbor[dim])
125  delete it.second;
126 
127  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
128  for (auto & it : _fe_face_neighbor[dim])
129  delete it.second;
130 
131  for (auto & it : _holder_qrule_volume)
132  delete it.second;
133 
134  for (auto & it : _holder_qrule_arbitrary)
135  delete it.second;
136 
137  for (auto & it : _holder_qface_arbitrary)
138  delete it.second;
139 
140  for (auto & it : _holder_qrule_face)
141  delete it.second;
142 
143  for (auto & it : _holder_qrule_neighbor)
144  delete it.second;
145 
146  for (auto & it : _fe_shape_data)
147  delete it.second;
148 
149  for (auto & it : _fe_shape_data_face)
150  delete it.second;
151 
152  for (auto & it : _fe_shape_data_neighbor)
153  delete it.second;
154 
155  for (auto & it : _fe_shape_data_face_neighbor)
156  delete it.second;
157 
158  delete _current_side_elem;
160 
162 
163  _coord.release();
165 }
166 
167 void
169 {
170  if (!_fe_shape_data[type])
172 
173  // Build an FE object for this type for each dimension up to the dimension of the current mesh
174  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
175  {
176  if (!_fe[dim][type])
177  _fe[dim][type] = FEBase::build(dim, type).release();
178  _fe[dim][type]->get_phi();
179  _fe[dim][type]->get_dphi();
180  // Pre-request xyz. We have always computed xyz, but due to
181  // recent optimizations in libmesh, we now need to explicity
182  // request it, since apps (Yak) may rely on it being computed.
183  _fe[dim][type]->get_xyz();
184  if (_need_second_derivative.find(type) != _need_second_derivative.end())
185  _fe[dim][type]->get_d2phi();
186  }
187 }
188 
189 void
191 {
192  if (!_fe_shape_data_face[type])
194 
195  // Build an FE object for this type for each dimension up to the dimension of the current mesh
196  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
197  {
198  if (!_fe_face[dim][type])
199  _fe_face[dim][type] = FEBase::build(dim, type).release();
200  _fe_face[dim][type]->get_phi();
201  _fe_face[dim][type]->get_dphi();
202  if (_need_second_derivative.find(type) != _need_second_derivative.end())
203  _fe_face[dim][type]->get_d2phi();
204  }
205 }
206 
207 void
209 {
210  if (!_fe_shape_data_neighbor[type])
212 
213  // Build an FE object for this type for each dimension up to the dimension of the current mesh
214  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
215  {
216  if (!_fe_neighbor[dim][type])
217  _fe_neighbor[dim][type] = FEBase::build(dim, type).release();
218  _fe_neighbor[dim][type]->get_phi();
219  _fe_neighbor[dim][type]->get_dphi();
220  if (_need_second_derivative.find(type) != _need_second_derivative.end())
221  _fe_neighbor[dim][type]->get_d2phi();
222  }
223 }
224 
225 void
227 {
228  if (!_fe_shape_data_face_neighbor[type])
230 
231  // Build an FE object for this type for each dimension up to the dimension of the current mesh
232  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
233  {
234  if (!_fe_face_neighbor[dim][type])
235  _fe_face_neighbor[dim][type] = FEBase::build(dim, type).release();
236  _fe_face_neighbor[dim][type]->get_phi();
237  _fe_face_neighbor[dim][type]->get_dphi();
238  if (_need_second_derivative.find(type) != _need_second_derivative.end())
239  _fe_face_neighbor[dim][type]->get_d2phi();
240  }
241 }
242 
243 FEBase *&
244 Assembly::getFE(FEType type, unsigned int dim)
245 {
246  buildFE(type);
247  return _fe[dim][type];
248 }
249 
250 FEBase *&
251 Assembly::getFEFace(FEType type, unsigned int dim)
252 {
253  buildFaceFE(type);
254  return _fe_face[dim][type];
255 }
256 
257 FEBase *&
258 Assembly::getFEFaceNeighbor(FEType type, unsigned int dim)
259 {
260  buildFaceNeighborFE(type);
261  return _fe_neighbor[dim][type];
262 }
263 
264 const VariablePhiValue &
266 {
267  buildFE(type);
268  return _fe_shape_data[type]->_phi;
269 }
270 
271 const VariablePhiGradient &
273 {
274  buildFE(type);
275  return _fe_shape_data[type]->_grad_phi;
276 }
277 
278 const VariablePhiSecond &
280 {
282  buildFE(type);
283  return _fe_shape_data[type]->_second_phi;
284 }
285 
286 const VariablePhiValue &
288 {
289  buildFaceFE(type);
290  return _fe_shape_data_face[type]->_phi;
291 }
292 
293 const VariablePhiGradient &
295 {
296  buildFaceFE(type);
297  return _fe_shape_data_face[type]->_grad_phi;
298 }
299 
300 const VariablePhiSecond &
302 {
304  buildFaceFE(type);
305  return _fe_shape_data_face[type]->_second_phi;
306 }
307 
308 const VariablePhiValue &
310 {
311  buildNeighborFE(type);
312  return _fe_shape_data_neighbor[type]->_phi;
313 }
314 
315 const VariablePhiGradient &
317 {
318  buildNeighborFE(type);
319  return _fe_shape_data_neighbor[type]->_grad_phi;
320 }
321 
322 const VariablePhiSecond &
324 {
326  buildNeighborFE(type);
327  return _fe_shape_data_neighbor[type]->_second_phi;
328 }
329 
330 const VariablePhiValue &
332 {
333  buildFaceNeighborFE(type);
334  return _fe_shape_data_face_neighbor[type]->_phi;
335 }
336 
337 const VariablePhiGradient &
339 {
340  buildFaceNeighborFE(type);
341  return _fe_shape_data_face_neighbor[type]->_grad_phi;
342 }
343 
344 const VariablePhiSecond &
346 {
348  buildFaceNeighborFE(type);
349  return _fe_shape_data_face_neighbor[type]->_second_phi;
350 }
351 
352 const Real &
354 {
357 }
358 
359 void
360 Assembly::createQRules(QuadratureType type, Order order, Order volume_order, Order face_order)
361 {
362  _holder_qrule_volume.clear();
363  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
364  _holder_qrule_volume[dim] = QBase::build(type, dim, volume_order).release();
365 
366  _holder_qrule_face.clear();
367  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
368  _holder_qrule_face[dim] = QBase::build(type, dim - 1, face_order).release();
369 
370  _holder_qrule_neighbor.clear();
371  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
372  _holder_qrule_neighbor[dim] = new ArbitraryQuadrature(dim, face_order);
373 
374  _holder_qrule_arbitrary.clear();
375  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
376  _holder_qrule_arbitrary[dim] = new ArbitraryQuadrature(dim, order);
377 }
378 
379 void
380 Assembly::setVolumeQRule(QBase * qrule, unsigned int dim)
381 {
382  _current_qrule = qrule;
383 
384  if (qrule) // Don't set a NULL qrule
385  {
386  for (auto & it : _fe[dim])
387  it.second->attach_quadrature_rule(_current_qrule);
388  }
389 }
390 
391 void
392 Assembly::setFaceQRule(QBase * qrule, unsigned int dim)
393 {
394  _current_qrule_face = qrule;
395 
396  for (auto & it : _fe_face[dim])
397  it.second->attach_quadrature_rule(_current_qrule_face);
398 }
399 
400 void
401 Assembly::setNeighborQRule(QBase * qrule, unsigned int dim)
402 {
403  _current_qrule_neighbor = qrule;
404 
405  for (auto & it : _fe_face_neighbor[dim])
406  it.second->attach_quadrature_rule(_current_qrule_neighbor);
407 }
408 
409 void
411 {
412  for (auto & it : _element_fe_shape_data_cache)
413  it.second->_invalidated = true;
414 }
415 
416 void
418 {
419  unsigned int dim = elem->dim();
420  ElementFEShapeData * efesd = NULL;
421 
422  // Whether or not we're going to do FE caching this time through
423  bool do_caching = _should_use_fe_cache && _currently_fe_caching;
424 
425  if (do_caching)
426  {
427  efesd = _element_fe_shape_data_cache[elem->id()];
428 
429  if (!efesd)
430  {
431  efesd = new ElementFEShapeData;
432  _element_fe_shape_data_cache[elem->id()] = efesd;
433  efesd->_invalidated = true;
434  }
435  }
436 
437  for (const auto & it : _fe[dim])
438  {
439  FEBase * fe = it.second;
440  const FEType & fe_type = it.first;
441 
442  _current_fe[fe_type] = fe;
443 
444  FEShapeData * fesd = _fe_shape_data[fe_type];
445 
446  FEShapeData * cached_fesd = NULL;
447  if (do_caching)
448  cached_fesd = efesd->_shape_data[fe_type];
449 
450  if (!cached_fesd || efesd->_invalidated)
451  {
452  fe->reinit(elem);
453 
454  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe->get_phi()));
455  fesd->_grad_phi.shallowCopy(
456  const_cast<std::vector<std::vector<RealGradient>> &>(fe->get_dphi()));
457  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
458  fesd->_second_phi.shallowCopy(
459  const_cast<std::vector<std::vector<RealTensor>> &>(fe->get_d2phi()));
460 
461  if (do_caching)
462  {
463  if (!cached_fesd)
464  {
465  cached_fesd = new FEShapeData;
466  efesd->_shape_data[fe_type] = cached_fesd;
467  }
468 
469  *cached_fesd = *fesd;
470  }
471  }
472  else // This means we have valid cached shape function values for this element / fe_type combo
473  {
474  fesd->_phi.shallowCopy(cached_fesd->_phi);
475  fesd->_grad_phi.shallowCopy(cached_fesd->_grad_phi);
476  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
477  fesd->_second_phi.shallowCopy(cached_fesd->_second_phi);
478  }
479  }
480 
481  // During that last loop the helper objects will have been reinitialized as well
482  // We need to dig out the q_points and JxW from it.
483  if (!do_caching || efesd->_invalidated)
484  {
486  const_cast<std::vector<Point> &>((*_holder_fe_helper[dim])->get_xyz()));
487  _current_JxW.shallowCopy(const_cast<std::vector<Real> &>((*_holder_fe_helper[dim])->get_JxW()));
488 
489  if (do_caching)
490  {
491  efesd->_q_points = _current_q_points;
492  efesd->_JxW = _current_JxW;
493  }
494  }
495  else // Use cached values
496  {
498  _current_JxW.shallowCopy(efesd->_JxW);
499  }
500 
501  if (do_caching)
502  efesd->_invalidated = false;
503 
504  if (_xfem != NULL)
506 }
507 
508 void
509 Assembly::reinitFEFace(const Elem * elem, unsigned int side)
510 {
511  unsigned int dim = elem->dim();
512 
513  for (const auto & it : _fe_face[dim])
514  {
515  FEBase * fe_face = it.second;
516  const FEType & fe_type = it.first;
517  FEShapeData * fesd = _fe_shape_data_face[fe_type];
518  fe_face->reinit(elem, side);
519  _current_fe_face[fe_type] = fe_face;
520 
521  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_face->get_phi()));
522  fesd->_grad_phi.shallowCopy(
523  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face->get_dphi()));
524  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
525  fesd->_second_phi.shallowCopy(
526  const_cast<std::vector<std::vector<RealTensor>> &>(fe_face->get_d2phi()));
527  }
528 
529  // During that last loop the helper objects will have been reinitialized as well
530  // We need to dig out the q_points and JxW from it.
532  const_cast<std::vector<Point> &>((*_holder_fe_face_helper[dim])->get_xyz()));
534  const_cast<std::vector<Real> &>((*_holder_fe_face_helper[dim])->get_JxW()));
536  const_cast<std::vector<Point> &>((*_holder_fe_face_helper[dim])->get_normals()));
537 }
538 
539 void
540 Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
541 {
542  unsigned int neighbor_dim = neighbor->dim();
543 
544  // reinit neighbor face
545  for (const auto & it : _fe_face_neighbor[neighbor_dim])
546  {
547  FEBase * fe_face_neighbor = it.second;
548  FEType fe_type = it.first;
549  FEShapeData * fesd = _fe_shape_data_face_neighbor[fe_type];
550 
551  fe_face_neighbor->reinit(neighbor, &reference_points);
552 
553  _current_fe_face_neighbor[fe_type] = fe_face_neighbor;
554 
555  fesd->_phi.shallowCopy(
556  const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor->get_phi()));
557  fesd->_grad_phi.shallowCopy(
558  const_cast<std::vector<std::vector<RealGradient>> &>(fe_face_neighbor->get_dphi()));
559  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
560  fesd->_second_phi.shallowCopy(
561  const_cast<std::vector<std::vector<RealTensor>> &>(fe_face_neighbor->get_d2phi()));
562  }
563 }
564 
565 void
566 Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
567 {
568  unsigned int neighbor_dim = neighbor->dim();
569 
570  // reinit neighbor face
571  for (const auto & it : _fe_neighbor[neighbor_dim])
572  {
573  FEBase * fe_neighbor = it.second;
574  FEType fe_type = it.first;
575  FEShapeData * fesd = _fe_shape_data_neighbor[fe_type];
576 
577  fe_neighbor->reinit(neighbor, &reference_points);
578 
579  _current_fe_neighbor[fe_type] = fe_neighbor;
580 
581  fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_neighbor->get_phi()));
582  fesd->_grad_phi.shallowCopy(
583  const_cast<std::vector<std::vector<RealGradient>> &>(fe_neighbor->get_dphi()));
584  if (_need_second_derivative.find(fe_type) != _need_second_derivative.end())
585  fesd->_second_phi.shallowCopy(
586  const_cast<std::vector<std::vector<RealTensor>> &>(fe_neighbor->get_d2phi()));
587  }
588 }
589 
590 void
591 Assembly::reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points)
592 {
593  unsigned int neighbor_dim = neighbor->dim();
594 
595  ArbitraryQuadrature * neighbor_rule = _holder_qrule_neighbor[neighbor_dim];
596  neighbor_rule->setPoints(reference_points);
597  setNeighborQRule(neighbor_rule, neighbor_dim);
598 
600  mooseAssert(_current_neighbor_subdomain_id == _current_neighbor_elem->subdomain_id(),
601  "current neighbor subdomain has been set incorrectly");
602 
603  // Calculate the volume of the neighbor
605  {
606  unsigned int dim = neighbor->dim();
607  FEBase * fe = *_holder_fe_neighbor_helper[dim];
608  QBase * qrule = _holder_qrule_volume[dim];
609 
610  fe->attach_quadrature_rule(qrule);
611  fe->reinit(neighbor);
612 
613  // set the coord transformation
614  _coord_neighbor.resize(qrule->n_points());
615  Moose::CoordinateSystemType coord_type =
617  unsigned int rz_radial_coord = _sys.subproblem().getAxisymmetricRadialCoord();
618  const std::vector<Real> & JxW = fe->get_JxW();
619  const std::vector<Point> & q_points = fe->get_xyz();
620 
621  switch (coord_type)
622  {
623  case Moose::COORD_XYZ:
624  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
625  _coord_neighbor[qp] = 1.;
626  break;
627 
628  case Moose::COORD_RZ:
629  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
630  _coord_neighbor[qp] = 2 * M_PI * q_points[qp](rz_radial_coord);
631  break;
632 
634  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
635  _coord_neighbor[qp] = 4 * M_PI * q_points[qp](0) * q_points[qp](0);
636  break;
637 
638  default:
639  mooseError("Unknown coordinate system");
640  break;
641  }
642 
644  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
645  _current_neighbor_volume += JxW[qp] * _coord_neighbor[qp];
646  }
647 }
648 
649 void
650 Assembly::reinit(const Elem * elem)
651 {
653  _current_neighbor_elem = nullptr;
654  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
655  "current subdomain has been set incorrectly");
657 
658  unsigned int elem_dimension = elem->dim();
659 
661 
662  // Make sure the qrule is the right one
664  setVolumeQRule(_current_qrule_volume, elem_dimension);
665 
666  _currently_fe_caching = true;
667 
668  reinitFE(elem);
669 
671 }
672 
673 void
674 Assembly::setCoordinateTransformation(const QBase * qrule, const MooseArray<Point> & q_points)
675 {
676  _coord.resize(qrule->n_points());
678  unsigned int rz_radial_coord = _sys.subproblem().getAxisymmetricRadialCoord();
679 
680  switch (_coord_type)
681  {
682  case Moose::COORD_XYZ:
683  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
684  _coord[qp] = 1.;
685  break;
686 
687  case Moose::COORD_RZ:
688  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
689  _coord[qp] = 2 * M_PI * q_points[qp](rz_radial_coord);
690  break;
691 
693  for (unsigned int qp = 0; qp < qrule->n_points(); qp++)
694  _coord[qp] = 4 * M_PI * q_points[qp](0) * q_points[qp](0);
695  break;
696 
697  default:
698  mooseError("Unknown coordinate system");
699  break;
700  }
701 }
702 
703 void
705 {
707  return;
708 
710 
712  for (unsigned int qp = 0; qp < _current_qrule->n_points(); qp++)
714 
716 }
717 
718 void
720 {
722  return;
723 
725 
727  for (unsigned int qp = 0; qp < _current_qrule_face->n_points(); qp++)
729 
731 }
732 
733 void
734 Assembly::reinitAtPhysical(const Elem * elem, const std::vector<Point> & physical_points)
735 {
737  _current_neighbor_elem = nullptr;
738  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
739  "current subdomain has been set incorrectly");
741 
742  FEInterface::inverse_map(elem->dim(),
743  (*_holder_fe_helper[elem->dim()])->get_fe_type(),
744  elem,
745  physical_points,
747 
748  _currently_fe_caching = false;
749 
751 
752  // Save off the physical points
753  _current_physical_points = physical_points;
754 }
755 
756 void
757 Assembly::reinit(const Elem * elem, const std::vector<Point> & reference_points)
758 {
760  _current_neighbor_elem = nullptr;
761  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
762  "current subdomain has been set incorrectly");
764 
765  unsigned int elem_dimension = _current_elem->dim();
766 
768 
769  // Make sure the qrule is the right one
771  setVolumeQRule(_current_qrule_arbitrary, elem_dimension);
772 
773  _current_qrule_arbitrary->setPoints(reference_points);
774 
775  _currently_fe_caching = false;
776 
777  reinitFE(elem);
778 
780 }
781 
782 void
783 Assembly::reinit(const Elem * elem, unsigned int side)
784 {
786  _current_neighbor_elem = nullptr;
787  mooseAssert(_current_subdomain_id == _current_elem->subdomain_id(),
788  "current subdomain has been set incorrectly");
792 
793  unsigned int elem_dimension = _current_elem->dim();
794 
795  if (_current_qrule_face != _holder_qrule_face[elem_dimension])
796  {
797  _current_qrule_face = _holder_qrule_face[elem_dimension];
798  setFaceQRule(_current_qrule_face, elem_dimension);
799  }
800 
801  if (_current_side_elem)
802  delete _current_side_elem;
803  _current_side_elem = elem->build_side_ptr(side).release();
804 
805  reinitFEFace(elem, side);
806 
808 }
809 
810 void
811 Assembly::reinit(const Node * node)
812 {
814  _current_neighbor_node = NULL;
815 }
816 
817 void
819 {
821 }
822 
823 void
825  unsigned int side,
826  const Elem * neighbor,
827  unsigned int neighbor_side)
828 {
829  _current_neighbor_side = neighbor_side;
830 
831  reinit(elem, side);
832 
833  unsigned int neighbor_dim = neighbor->dim();
834 
835  std::vector<Point> reference_points;
836  FEInterface::inverse_map(
837  neighbor_dim, FEType(), neighbor, _current_q_points_face.stdVector(), reference_points);
838 
839  reinitFEFaceNeighbor(neighbor, reference_points);
840  reinitNeighbor(neighbor, reference_points);
841 }
842 
843 void
845  unsigned int neighbor_side,
846  const std::vector<Point> & physical_points)
847 {
849  _current_neighbor_side_elem = neighbor->build_side_ptr(neighbor_side).release();
850 
851  std::vector<Point> reference_points;
852 
853  unsigned int neighbor_dim = neighbor->dim();
854  FEInterface::inverse_map(neighbor_dim, FEType(), neighbor, physical_points, reference_points);
855 
856  // first do the side element
858  reinitNeighbor(_current_neighbor_side_elem, reference_points);
859  // compute JxW on the neighbor's face
860  unsigned int neighbor_side_dim = _current_neighbor_side_elem->dim();
861  _current_JxW_neighbor.shallowCopy(const_cast<std::vector<Real> &>(
862  (*_holder_fe_face_neighbor_helper[neighbor_side_dim])->get_JxW()));
863 
864  reinitFEFaceNeighbor(neighbor, reference_points);
865  reinitNeighbor(neighbor, reference_points);
866 
867  // Save off the physical points
868  _current_physical_points = physical_points;
869 }
870 
871 void
873  const std::vector<Point> & physical_points)
874 {
875  std::vector<Point> reference_points;
876 
877  unsigned int neighbor_dim = neighbor->dim();
878  FEInterface::inverse_map(neighbor_dim, FEType(), neighbor, physical_points, reference_points);
879 
880  reinitFENeighbor(neighbor, reference_points);
881  reinitNeighbor(neighbor, reference_points);
882  // Save off the physical points
883  _current_physical_points = physical_points;
884 }
885 
886 DenseMatrix<Number> &
887 Assembly::jacobianBlock(unsigned int ivar, unsigned int jvar)
888 {
889  _jacobian_block_used[ivar][jvar] = 1;
890  return _sub_Kee[ivar][_block_diagonal_matrix ? 0 : jvar];
891 }
892 
893 DenseMatrix<Number> &
894 Assembly::jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar)
895 {
896  _jacobian_block_nonlocal_used[ivar][jvar] = 1;
897  return _sub_Keg[ivar][_block_diagonal_matrix ? 0 : jvar];
898 }
899 
900 DenseMatrix<Number> &
901 Assembly::jacobianBlockNeighbor(Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar)
902 {
903  _jacobian_block_neighbor_used[ivar][jvar] = 1;
905  {
906  switch (type)
907  {
908  default:
910  return _sub_Kee[ivar][0];
912  return _sub_Ken[ivar][0];
914  return _sub_Kne[ivar][0];
916  return _sub_Knn[ivar][0];
917  }
918  }
919  else
920  {
921  switch (type)
922  {
923  default:
925  return _sub_Kee[ivar][jvar];
927  return _sub_Ken[ivar][jvar];
929  return _sub_Kne[ivar][jvar];
931  return _sub_Knn[ivar][jvar];
932  }
933  }
934 }
935 
936 void
937 Assembly::init(const CouplingMatrix * cm)
938 {
939  _cm = cm;
940 
941  unsigned int n_vars = _sys.nVariables();
942 
943  // I want the blocks to go by columns first to reduce copying of shape function in assembling
944  // "full" Jacobian
945  _cm_entry.clear();
946  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
947  for (const auto & jvar : vars)
948  {
949  unsigned int j = jvar->number();
950  for (const auto & ivar : vars)
951  {
952  unsigned int i = ivar->number();
953  if ((*_cm)(i, j) != 0)
954  _cm_entry.push_back(std::make_pair(ivar, jvar));
955  }
956  }
957 
958  unsigned int max_rows_per_column = 0; // If this is 1 then we are using a block-diagonal
959  // preconditioner and we can save a bunch of memory.
960  for (unsigned int i = 0; i < n_vars; i++)
961  {
962  unsigned int max_rows_per_this_column = 0;
963  for (unsigned int j = 0; j < n_vars; j++)
964  {
965  if ((*_cm)(i, j) != 0)
966  max_rows_per_this_column++;
967  }
968  max_rows_per_column = std::max(max_rows_per_column, max_rows_per_this_column);
969  }
970 
971  if (max_rows_per_column == 1 && _sys.getScalarVariables(_tid).size() == 0)
972  _block_diagonal_matrix = true;
973 
974  // two vectors: one for time residual contributions and one for non-time residual contributions
975  _sub_Re.resize(2);
976  _sub_Rn.resize(2);
977  for (unsigned int i = 0; i < _sub_Re.size(); i++)
978  {
979  _sub_Re[i].resize(n_vars);
980  _sub_Rn[i].resize(n_vars);
981  }
982 
983  _sub_Kee.resize(n_vars);
984  _sub_Keg.resize(n_vars);
985  _sub_Ken.resize(n_vars);
986  _sub_Kne.resize(n_vars);
987  _sub_Knn.resize(n_vars);
988  _jacobian_block_used.resize(n_vars);
989  _jacobian_block_nonlocal_used.resize(n_vars);
990  _jacobian_block_neighbor_used.resize(n_vars);
991 
992  for (unsigned int i = 0; i < n_vars; ++i)
993  {
995  {
996  _sub_Kee[i].resize(n_vars);
997  _sub_Keg[i].resize(n_vars);
998  _sub_Ken[i].resize(n_vars);
999  _sub_Kne[i].resize(n_vars);
1000  _sub_Knn[i].resize(n_vars);
1001  }
1002  else
1003  {
1004  _sub_Kee[i].resize(1);
1005  _sub_Keg[i].resize(1);
1006  _sub_Ken[i].resize(1);
1007  _sub_Kne[i].resize(1);
1008  _sub_Knn[i].resize(1);
1009  }
1010  _jacobian_block_used[i].resize(n_vars);
1011  _jacobian_block_nonlocal_used[i].resize(n_vars);
1012  _jacobian_block_neighbor_used[i].resize(n_vars);
1013  }
1014 }
1015 
1016 void
1018 {
1019  _cm_nonlocal_entry.clear();
1020  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1021  for (const auto & jvar : vars)
1022  {
1023  unsigned int j = jvar->number();
1024  for (const auto & ivar : vars)
1025  {
1026  unsigned int i = ivar->number();
1027  if (_nonlocal_cm(i, j) != 0)
1028  _cm_nonlocal_entry.push_back(std::make_pair(ivar, jvar));
1029  }
1030  }
1031 }
1032 
1033 void
1035 {
1036  for (const auto & it : _cm_entry)
1037  {
1038  MooseVariable & ivar = *(it.first);
1039  MooseVariable & jvar = *(it.second);
1040 
1041  unsigned int vi = ivar.number();
1042  unsigned int vj = jvar.number();
1043 
1044  jacobianBlock(vi, vj).resize(ivar.dofIndices().size(), jvar.dofIndices().size());
1045  jacobianBlock(vi, vj).zero();
1046  _jacobian_block_used[vi][vj] = 0;
1047  }
1048 
1049  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1050  for (const auto & var : vars)
1051  for (unsigned int i = 0; i < _sub_Re.size(); i++)
1052  {
1053  _sub_Re[i][var->number()].resize(var->dofIndices().size());
1054  _sub_Re[i][var->number()].zero();
1055  }
1056 }
1057 
1058 void
1060 {
1061  for (const auto & it : _cm_nonlocal_entry)
1062  {
1063  MooseVariable & ivar = *(it.first);
1064  MooseVariable & jvar = *(it.second);
1065 
1066  unsigned int vi = ivar.number();
1067  unsigned int vj = jvar.number();
1068 
1069  jacobianBlockNonlocal(vi, vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
1070  jacobianBlockNonlocal(vi, vj).zero();
1071  _jacobian_block_nonlocal_used[vi][vj] = 0;
1072  }
1073 }
1074 
1075 void
1077 {
1078  for (const auto & it : _cm_entry)
1079  {
1080  MooseVariable & ivar = *(it.first);
1081  MooseVariable & jvar = *(it.second);
1082 
1083  unsigned int vi = ivar.number();
1084  unsigned int vj = jvar.number();
1085 
1086  if (vi == var->number() || vj == var->number())
1087  jacobianBlock(vi, vj).resize(ivar.dofIndices().size(), jvar.dofIndices().size());
1088  }
1089 
1090  for (unsigned int i = 0; i < _sub_Re.size(); i++)
1091  {
1092  _sub_Re[i][var->number()].resize(var->dofIndices().size());
1093  _sub_Re[i][var->number()].zero();
1094  }
1095 }
1096 
1097 void
1099 {
1100  for (const auto & it : _cm_nonlocal_entry)
1101  {
1102  MooseVariable & ivar = *(it.first);
1103  MooseVariable & jvar = *(it.second);
1104 
1105  unsigned int vi = ivar.number();
1106  unsigned int vj = jvar.number();
1107 
1108  if (vi == var->number() || vj == var->number())
1109  jacobianBlockNonlocal(vi, vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
1110  }
1111 }
1112 
1113 void
1115 {
1116  for (const auto & it : _cm_entry)
1117  {
1118  MooseVariable & ivar = *(it.first);
1119  MooseVariable & jvar = *(it.second);
1120 
1121  unsigned int vi = ivar.number();
1122  unsigned int vj = jvar.number();
1123 
1125  .resize(ivar.dofIndices().size(), jvar.dofIndicesNeighbor().size());
1127 
1129  .resize(ivar.dofIndicesNeighbor().size(), jvar.dofIndices().size());
1131 
1133  .resize(ivar.dofIndicesNeighbor().size(), jvar.dofIndicesNeighbor().size());
1135 
1136  _jacobian_block_neighbor_used[vi][vj] = 0;
1137  }
1138 
1139  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1140  for (const auto & var : vars)
1141  for (unsigned int i = 0; i < _sub_Rn.size(); i++)
1142  {
1143  _sub_Rn[i][var->number()].resize(var->dofIndicesNeighbor().size());
1144  _sub_Rn[i][var->number()].zero();
1145  }
1146 }
1147 
1148 void
1149 Assembly::prepareBlock(unsigned int ivar,
1150  unsigned int jvar,
1151  const std::vector<dof_id_type> & dof_indices)
1152 {
1153  jacobianBlock(ivar, jvar).resize(dof_indices.size(), dof_indices.size());
1154  jacobianBlock(ivar, jvar).zero();
1155  _jacobian_block_used[ivar][jvar] = 0;
1156 
1157  for (unsigned int i = 0; i < _sub_Re.size(); i++)
1158  {
1159  _sub_Re[i][ivar].resize(dof_indices.size());
1160  _sub_Re[i][ivar].zero();
1161  }
1162 }
1163 
1164 void
1166  unsigned int jvar,
1167  const std::vector<dof_id_type> & idof_indices,
1168  const std::vector<dof_id_type> & jdof_indices)
1169 {
1170  jacobianBlockNonlocal(ivar, jvar).resize(idof_indices.size(), jdof_indices.size());
1171  jacobianBlockNonlocal(ivar, jvar).zero();
1172  _jacobian_block_nonlocal_used[ivar][jvar] = 0;
1173 }
1174 
1175 void
1177 {
1178  const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
1179  for (const auto & ivar : vars)
1180  {
1181  unsigned int idofs = ivar->dofIndices().size();
1182 
1183  for (unsigned int i = 0; i < _sub_Re.size(); i++)
1184  {
1185  _sub_Re[i][ivar->number()].resize(idofs);
1186  _sub_Re[i][ivar->number()].zero();
1187  }
1188 
1189  for (const auto & jvar : vars)
1190  {
1191  unsigned int jdofs = jvar->dofIndices().size();
1192 
1193  jacobianBlock(ivar->number(), jvar->number()).resize(idofs, jdofs);
1194  jacobianBlock(ivar->number(), jvar->number()).zero();
1195  _jacobian_block_used[ivar->number()][jvar->number()] = 0;
1196  }
1197  }
1198 }
1199 
1200 void
1202 {
1203  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1204  const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
1205 
1206  for (const auto & ivar : scalar_vars)
1207  {
1208  unsigned int idofs = ivar->dofIndices().size();
1209 
1210  for (const auto & jvar : vars)
1211  {
1212  unsigned int jdofs = jvar->dofIndices().size();
1213 
1214  jacobianBlock(ivar->number(), jvar->number()).resize(idofs, jdofs);
1215  jacobianBlock(ivar->number(), jvar->number()).zero();
1216  _jacobian_block_used[ivar->number()][jvar->number()] = 0;
1217 
1218  jacobianBlock(jvar->number(), ivar->number()).resize(jdofs, idofs);
1219  jacobianBlock(jvar->number(), ivar->number()).zero();
1220  _jacobian_block_used[jvar->number()][ivar->number()] = 0;
1221  }
1222  }
1223 }
1224 
1225 void
1226 Assembly::copyShapes(unsigned int var)
1227 {
1228  MooseVariable & v = _sys.getVariable(_tid, var);
1229 
1230  _phi.shallowCopy(v.phi());
1232 
1233  if (v.computingSecond())
1235 }
1236 
1237 void
1238 Assembly::copyFaceShapes(unsigned int var)
1239 {
1240  MooseVariable & v = _sys.getVariable(_tid, var);
1241 
1244 
1245  if (v.computingSecond())
1247 }
1248 
1249 void
1251 {
1252  MooseVariable & v = _sys.getVariable(_tid, var);
1253 
1254  if (v.usesPhi())
1256  if (v.usesGradPhi())
1258  if (v.usesSecondPhi())
1260 
1261  if (v.usesPhi())
1263  if (v.usesGradPhi())
1265  if (v.usesSecondPhi())
1267 }
1268 
1269 void
1270 Assembly::addResidualBlock(NumericVector<Number> & residual,
1271  DenseVector<Number> & res_block,
1272  const std::vector<dof_id_type> & dof_indices,
1273  Real scaling_factor)
1274 {
1275  if (dof_indices.size() > 0 && res_block.size())
1276  {
1277  _temp_dof_indices = dof_indices;
1278  _dof_map.constrain_element_vector(res_block, _temp_dof_indices, false);
1279 
1280  if (scaling_factor != 1.0)
1281  {
1282  _tmp_Re = res_block;
1283  _tmp_Re *= scaling_factor;
1284  residual.add_vector(_tmp_Re, _temp_dof_indices);
1285  }
1286  else
1287  {
1288  residual.add_vector(res_block, _temp_dof_indices);
1289  }
1290  }
1291 }
1292 
1293 void
1294 Assembly::cacheResidualBlock(std::vector<Real> & cached_residual_values,
1295  std::vector<dof_id_type> & cached_residual_rows,
1296  DenseVector<Number> & res_block,
1297  std::vector<dof_id_type> & dof_indices,
1298  Real scaling_factor)
1299 {
1300  if (dof_indices.size() > 0 && res_block.size())
1301  {
1302  _temp_dof_indices = dof_indices;
1303  _dof_map.constrain_element_vector(res_block, _temp_dof_indices, false);
1304 
1305  if (scaling_factor != 1.0)
1306  {
1307  _tmp_Re = res_block;
1308  _tmp_Re *= scaling_factor;
1309 
1310  for (unsigned int i = 0; i < _tmp_Re.size(); i++)
1311  {
1312  cached_residual_values.push_back(_tmp_Re(i));
1313  cached_residual_rows.push_back(_temp_dof_indices[i]);
1314  }
1315  }
1316  else
1317  {
1318  for (unsigned int i = 0; i < res_block.size(); i++)
1319  {
1320  cached_residual_values.push_back(res_block(i));
1321  cached_residual_rows.push_back(_temp_dof_indices[i]);
1322  }
1323  }
1324  }
1325 
1326  res_block.zero();
1327 }
1328 
1329 void
1330 Assembly::addResidual(NumericVector<Number> & residual,
1331  Moose::KernelType type /* = Moose::KT_NONTIME*/)
1332 {
1333  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1334  for (const auto & var : vars)
1336  residual, _sub_Re[type][var->number()], var->dofIndices(), var->scalingFactor());
1337 }
1338 
1339 void
1340 Assembly::addResidualNeighbor(NumericVector<Number> & residual,
1341  Moose::KernelType type /* = Moose::KT_NONTIME*/)
1342 {
1343  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1344  for (const auto & var : vars)
1346  residual, _sub_Rn[type][var->number()], var->dofIndicesNeighbor(), var->scalingFactor());
1347 }
1348 
1349 void
1350 Assembly::addResidualScalar(NumericVector<Number> & residual,
1351  Moose::KernelType type /* = Moose::KT_NONTIME*/)
1352 {
1353  // add the scalar variables residuals
1354  const std::vector<MooseVariableScalar *> & vars = _sys.getScalarVariables(_tid);
1355  for (const auto & var : vars)
1357  residual, _sub_Re[type][var->number()], var->dofIndices(), var->scalingFactor());
1358 }
1359 
1360 void
1362 {
1363  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1364  for (const auto & var : vars)
1365  for (unsigned int i = 0; i < _sub_Re.size(); i++)
1369  _sub_Re[i][var->number()],
1370  var->dofIndices(),
1371  var->scalingFactor());
1372 }
1373 
1374 void
1376 {
1377  _cached_residual_values[type].push_back(value);
1378  _cached_residual_rows[type].push_back(dof);
1379 }
1380 
1381 void
1383 {
1384  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1385  for (const auto & var : vars)
1386  for (unsigned int i = 0; i < _sub_Re.size(); i++)
1390  _sub_Rn[i][var->number()],
1391  var->dofIndicesNeighbor(),
1392  var->scalingFactor());
1393 }
1394 
1395 void
1396 Assembly::cacheResidualNodes(DenseVector<Number> & res, std::vector<dof_id_type> & dof_index)
1397 {
1398  // Add the residual value and dof_index to cached_residual_values and cached_residual_rows
1399  // respectively.
1400  // This is used by NodalConstraint.C to cache the residual calculated for master and slave node.
1402  for (unsigned int i = 0; i < dof_index.size(); ++i)
1403  {
1404  _cached_residual_values[type].push_back(res(i));
1405  _cached_residual_rows[type].push_back(dof_index[i]);
1406  }
1407 }
1408 
1409 void
1411 {
1412  for (unsigned int i = 0; i < _sub_Re.size(); i++)
1413  {
1415  if (!_sys.hasResidualVector(type))
1416  {
1417  _cached_residual_values[i].clear();
1418  _cached_residual_rows[i].clear();
1419  continue;
1420  }
1422  }
1423 }
1424 
1425 void
1426 Assembly::addCachedResidual(NumericVector<Number> & residual, Moose::KernelType type)
1427 {
1428  if (!_sys.hasResidualVector(type))
1429  {
1430  _cached_residual_values[type].clear();
1431  _cached_residual_rows[type].clear();
1432  return;
1433  }
1434 
1435  std::vector<Real> & cached_residual_values = _cached_residual_values[type];
1436  std::vector<dof_id_type> & cached_residual_rows = _cached_residual_rows[type];
1437 
1438  mooseAssert(cached_residual_values.size() == cached_residual_rows.size(),
1439  "Number of cached residuals and number of rows must match!");
1440 
1441  residual.add_vector(cached_residual_values, cached_residual_rows);
1442 
1443  if (_max_cached_residuals < cached_residual_values.size())
1444  _max_cached_residuals = cached_residual_values.size();
1445 
1446  // Try to be more efficient from now on
1447  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
1448  cached_residual_values.clear();
1449  cached_residual_values.reserve(_max_cached_residuals * 2);
1450 
1451  cached_residual_rows.clear();
1452  cached_residual_rows.reserve(_max_cached_residuals * 2);
1453 }
1454 
1455 void
1456 Assembly::setResidualBlock(NumericVector<Number> & residual,
1457  DenseVector<Number> & res_block,
1458  std::vector<dof_id_type> & dof_indices,
1459  Real scaling_factor)
1460 {
1461  if (dof_indices.size() > 0)
1462  {
1463  std::vector<dof_id_type> di(dof_indices);
1464  _dof_map.constrain_element_vector(res_block, di, false);
1465 
1466  if (scaling_factor != 1.0)
1467  {
1468  _tmp_Re = res_block;
1469  _tmp_Re *= scaling_factor;
1470  residual.insert(_tmp_Re, di);
1471  }
1472  else
1473  residual.insert(res_block, di);
1474  }
1475 }
1476 
1477 void
1478 Assembly::setResidual(NumericVector<Number> & residual,
1479  Moose::KernelType type /* = Moose::KT_NONTIME*/)
1480 {
1481  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1482  for (const auto & var : vars)
1484  residual, _sub_Re[type][var->number()], var->dofIndices(), var->scalingFactor());
1485 }
1486 
1487 void
1488 Assembly::setResidualNeighbor(NumericVector<Number> & residual,
1489  Moose::KernelType type /* = Moose::KT_NONTIME*/)
1490 {
1491  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1492  for (const auto & var : vars)
1494  residual, _sub_Rn[type][var->number()], var->dofIndicesNeighbor(), var->scalingFactor());
1495 }
1496 
1497 void
1498 Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
1499  DenseMatrix<Number> & jac_block,
1500  const std::vector<dof_id_type> & idof_indices,
1501  const std::vector<dof_id_type> & jdof_indices,
1502  Real scaling_factor)
1503 {
1504  if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m())
1505  {
1506  std::vector<dof_id_type> di(idof_indices);
1507  std::vector<dof_id_type> dj(jdof_indices);
1508  _dof_map.constrain_element_matrix(jac_block, di, dj, false);
1509 
1510  if (scaling_factor != 1.0)
1511  {
1512  _tmp_Ke = jac_block;
1513  _tmp_Ke *= scaling_factor;
1514  jacobian.add_matrix(_tmp_Ke, di, dj);
1515  }
1516  else
1517  jacobian.add_matrix(jac_block, di, dj);
1518  }
1519 }
1520 
1521 void
1522 Assembly::cacheJacobianBlock(DenseMatrix<Number> & jac_block,
1523  std::vector<dof_id_type> & idof_indices,
1524  std::vector<dof_id_type> & jdof_indices,
1525  Real scaling_factor)
1526 {
1527  if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m())
1528  {
1529  std::vector<dof_id_type> di(idof_indices);
1530  std::vector<dof_id_type> dj(jdof_indices);
1531  _dof_map.constrain_element_matrix(jac_block, di, dj, false);
1532 
1533  if (scaling_factor != 1.0)
1534  jac_block *= scaling_factor;
1535 
1536  for (unsigned int i = 0; i < di.size(); i++)
1537  for (unsigned int j = 0; j < dj.size(); j++)
1538  {
1539  _cached_jacobian_values.push_back(jac_block(i, j));
1540  _cached_jacobian_rows.push_back(di[i]);
1541  _cached_jacobian_cols.push_back(dj[j]);
1542  }
1543  }
1544  jac_block.zero();
1545 }
1546 
1547 void
1548 Assembly::cacheJacobianBlockNonlocal(DenseMatrix<Number> & jac_block,
1549  const std::vector<dof_id_type> & idof_indices,
1550  const std::vector<dof_id_type> & jdof_indices,
1551  Real scaling_factor)
1552 {
1553  if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.n() && jac_block.m())
1554  {
1555  std::vector<dof_id_type> di(idof_indices);
1556  std::vector<dof_id_type> dj(jdof_indices);
1557  _dof_map.constrain_element_matrix(jac_block, di, dj, false);
1558 
1559  if (scaling_factor != 1.0)
1560  jac_block *= scaling_factor;
1561 
1562  for (unsigned int i = 0; i < di.size(); i++)
1563  for (unsigned int j = 0; j < dj.size(); j++)
1564  if (jac_block(i, j) != 0.0) // no storage allocated for unimplemented jacobian terms,
1565  // maintaining maximum sparsity possible
1566  {
1567  _cached_jacobian_values.push_back(jac_block(i, j));
1568  _cached_jacobian_rows.push_back(di[i]);
1569  _cached_jacobian_cols.push_back(dj[j]);
1570  }
1571  }
1572  jac_block.zero();
1573 }
1574 
1575 void
1576 Assembly::addCachedJacobian(SparseMatrix<Number> & jacobian)
1577 {
1579  mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
1580  "Error: Cached data sizes MUST be the same!");
1581 
1582  for (unsigned int i = 0; i < _cached_jacobian_rows.size(); i++)
1584 
1587 
1588  // Try to be more efficient from now on
1589  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
1590  _cached_jacobian_values.clear();
1592 
1593  _cached_jacobian_rows.clear();
1595 
1596  _cached_jacobian_cols.clear();
1598 }
1599 
1600 void
1601 Assembly::addJacobian(SparseMatrix<Number> & jacobian)
1602 {
1603  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1604  for (const auto & ivar : vars)
1605  for (const auto & jvar : vars)
1606  if ((*_cm)(ivar->number(), jvar->number()) != 0 &&
1607  _jacobian_block_used[ivar->number()][jvar->number()])
1608  addJacobianBlock(jacobian,
1609  jacobianBlock(ivar->number(), jvar->number()),
1610  ivar->dofIndices(),
1611  jvar->dofIndices(),
1612  ivar->scalingFactor());
1613 
1614  // Possibly add jacobian contributions from off-diagonal blocks coming from the scalar variables
1615  if (_sys.getScalarVariables(_tid).size() > 0)
1616  {
1617  const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
1618  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1619  for (const auto & ivar : scalar_vars)
1620  {
1621  for (const auto & jvar : vars)
1622  {
1623  // We only add jacobian blocks for d(nl-var)/d(scalar-var) now (these we generated by
1624  // kernels and BCs)
1625  // jacobian blocks d(scalar-var)/d(nl-var) are added later with addJacobianScalar
1626  if ((*_cm)(jvar->number(), ivar->number()) != 0 &&
1627  _jacobian_block_used[jvar->number()][ivar->number()])
1628  addJacobianBlock(jacobian,
1629  jacobianBlock(jvar->number(), ivar->number()),
1630  jvar->dofIndices(),
1631  ivar->dofIndices(),
1632  jvar->scalingFactor());
1633  if ((*_cm)(ivar->number(), jvar->number()) != 0 &&
1634  _jacobian_block_used[ivar->number()][jvar->number()])
1635  addJacobianBlock(jacobian,
1636  jacobianBlock(ivar->number(), jvar->number()),
1637  ivar->dofIndices(),
1638  jvar->dofIndices(),
1639  ivar->scalingFactor());
1640  }
1641  }
1642  }
1643 }
1644 
1645 void
1646 Assembly::addJacobianNonlocal(SparseMatrix<Number> & jacobian)
1647 {
1648  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1649  for (const auto & ivar : vars)
1650  for (const auto & jvar : vars)
1651  if (_nonlocal_cm(ivar->number(), jvar->number()) != 0 &&
1652  _jacobian_block_nonlocal_used[ivar->number()][jvar->number()])
1653  addJacobianBlock(jacobian,
1654  jacobianBlockNonlocal(ivar->number(), jvar->number()),
1655  ivar->dofIndices(),
1656  jvar->allDofIndices(),
1657  ivar->scalingFactor());
1658 }
1659 
1660 void
1661 Assembly::addJacobianNeighbor(SparseMatrix<Number> & jacobian)
1662 {
1663  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1664  for (const auto & ivar : vars)
1665  for (const auto & jvar : vars)
1666  if ((*_cm)(ivar->number(), jvar->number()) != 0 &&
1667  _jacobian_block_neighbor_used[ivar->number()][jvar->number()])
1668  {
1670  jacobian,
1671  jacobianBlockNeighbor(Moose::ElementNeighbor, ivar->number(), jvar->number()),
1672  ivar->dofIndices(),
1673  jvar->dofIndicesNeighbor(),
1674  ivar->scalingFactor());
1675 
1677  jacobian,
1678  jacobianBlockNeighbor(Moose::NeighborElement, ivar->number(), jvar->number()),
1679  ivar->dofIndicesNeighbor(),
1680  jvar->dofIndices(),
1681  ivar->scalingFactor());
1682 
1684  jacobian,
1685  jacobianBlockNeighbor(Moose::NeighborNeighbor, ivar->number(), jvar->number()),
1686  ivar->dofIndicesNeighbor(),
1687  jvar->dofIndicesNeighbor(),
1688  ivar->scalingFactor());
1689  }
1690 }
1691 
1692 void
1694 {
1695  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1696  for (const auto & ivar : vars)
1697  for (const auto & jvar : vars)
1698  if ((*_cm)(ivar->number(), jvar->number()) != 0 &&
1699  _jacobian_block_used[ivar->number()][jvar->number()])
1700  cacheJacobianBlock(jacobianBlock(ivar->number(), jvar->number()),
1701  ivar->dofIndices(),
1702  jvar->dofIndices(),
1703  ivar->scalingFactor());
1704 
1705  // Possibly add jacobian contributions from off-diagonal blocks coming from the scalar variables
1706  if (_sys.getScalarVariables(_tid).size() > 0)
1707  {
1708  const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
1709  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1710  for (std::vector<MooseVariableScalar *>::const_iterator it = scalar_vars.begin();
1711  it != scalar_vars.end();
1712  ++it)
1713  {
1714  MooseVariableScalar & ivar = *(*it);
1715  for (std::vector<MooseVariable *>::const_iterator jt = vars.begin(); jt != vars.end(); ++jt)
1716  {
1717  MooseVariable & jvar = *(*jt);
1718  // We only add jacobian blocks for d(nl-var)/d(scalar-var) now (these we generated by
1719  // kernels and BCs)
1720  // jacobian blocks d(scalar-var)/d(nl-var) are added later with addJacobianScalar
1721  if ((*_cm)(jvar.number(), ivar.number()) != 0 &&
1722  _jacobian_block_used[jvar.number()][ivar.number()])
1724  jvar.dofIndices(),
1725  ivar.dofIndices(),
1726  jvar.scalingFactor());
1727  if ((*_cm)(ivar.number(), jvar.number()) != 0 &&
1728  _jacobian_block_used[ivar.number()][jvar.number()])
1730  ivar.dofIndices(),
1731  jvar.dofIndices(),
1732  ivar.scalingFactor());
1733  }
1734  }
1735  }
1736 }
1737 
1738 void
1740 {
1741  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1742  for (const auto & ivar : vars)
1743  for (const auto & jvar : vars)
1744  if (_nonlocal_cm(ivar->number(), jvar->number()) != 0 &&
1745  _jacobian_block_nonlocal_used[ivar->number()][jvar->number()])
1746  cacheJacobianBlockNonlocal(jacobianBlockNonlocal(ivar->number(), jvar->number()),
1747  ivar->dofIndices(),
1748  jvar->allDofIndices(),
1749  ivar->scalingFactor());
1750 }
1751 
1752 void
1754 {
1755  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1756  for (const auto & ivar : vars)
1757  for (const auto & jvar : vars)
1758  if ((*_cm)(ivar->number(), jvar->number()) != 0 &&
1759  _jacobian_block_neighbor_used[ivar->number()][jvar->number()])
1760  {
1762  jacobianBlockNeighbor(Moose::ElementNeighbor, ivar->number(), jvar->number()),
1763  ivar->dofIndices(),
1764  jvar->dofIndicesNeighbor(),
1765  ivar->scalingFactor());
1767  jacobianBlockNeighbor(Moose::NeighborElement, ivar->number(), jvar->number()),
1768  ivar->dofIndicesNeighbor(),
1769  jvar->dofIndices(),
1770  ivar->scalingFactor());
1772  jacobianBlockNeighbor(Moose::NeighborNeighbor, ivar->number(), jvar->number()),
1773  ivar->dofIndicesNeighbor(),
1774  jvar->dofIndicesNeighbor(),
1775  ivar->scalingFactor());
1776  }
1777 }
1778 
1779 void
1780 Assembly::addJacobianBlock(SparseMatrix<Number> & jacobian,
1781  unsigned int ivar,
1782  unsigned int jvar,
1783  const DofMap & dof_map,
1784  std::vector<dof_id_type> & dof_indices)
1785 {
1786  DenseMatrix<Number> & ke = jacobianBlock(ivar, jvar);
1787 
1788  // stick it into the matrix
1789  std::vector<dof_id_type> di(dof_indices);
1790  dof_map.constrain_element_matrix(ke, di, false);
1791 
1792  Real scaling_factor = _sys.getVariable(_tid, ivar).scalingFactor();
1793  if (scaling_factor != 1.0)
1794  {
1795  _tmp_Ke = ke;
1796  _tmp_Ke *= scaling_factor;
1797  jacobian.add_matrix(_tmp_Ke, di);
1798  }
1799  else
1800  jacobian.add_matrix(ke, di);
1801 }
1802 
1803 void
1804 Assembly::addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
1805  unsigned int ivar,
1806  unsigned int jvar,
1807  const DofMap & dof_map,
1808  const std::vector<dof_id_type> & idof_indices,
1809  const std::vector<dof_id_type> & jdof_indices)
1810 {
1811  DenseMatrix<Number> & keg = jacobianBlockNonlocal(ivar, jvar);
1812 
1813  std::vector<dof_id_type> di(idof_indices);
1814  std::vector<dof_id_type> dg(jdof_indices);
1815  dof_map.constrain_element_matrix(keg, di, dg, false);
1816 
1817  Real scaling_factor = _sys.getVariable(_tid, ivar).scalingFactor();
1818  if (scaling_factor != 1.0)
1819  {
1820  _tmp_Ke = keg;
1821  _tmp_Ke *= scaling_factor;
1822  jacobian.add_matrix(_tmp_Ke, di, dg);
1823  }
1824  else
1825  jacobian.add_matrix(keg, di, dg);
1826 }
1827 
1828 void
1829 Assembly::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
1830  unsigned int ivar,
1831  unsigned int jvar,
1832  const DofMap & dof_map,
1833  std::vector<dof_id_type> & dof_indices,
1834  std::vector<dof_id_type> & neighbor_dof_indices)
1835 {
1836  DenseMatrix<Number> & kee = jacobianBlock(ivar, jvar);
1837  DenseMatrix<Number> & ken = jacobianBlockNeighbor(Moose::ElementNeighbor, ivar, jvar);
1838  DenseMatrix<Number> & kne = jacobianBlockNeighbor(Moose::NeighborElement, ivar, jvar);
1839  DenseMatrix<Number> & knn = jacobianBlockNeighbor(Moose::NeighborNeighbor, ivar, jvar);
1840 
1841  std::vector<dof_id_type> di(dof_indices);
1842  std::vector<dof_id_type> dn(neighbor_dof_indices);
1843  // stick it into the matrix
1844  dof_map.constrain_element_matrix(kee, di, false);
1845  dof_map.constrain_element_matrix(ken, di, dn, false);
1846  dof_map.constrain_element_matrix(kne, dn, di, false);
1847  dof_map.constrain_element_matrix(knn, dn, false);
1848 
1849  Real scaling_factor = _sys.getVariable(_tid, ivar).scalingFactor();
1850  if (scaling_factor != 1.0)
1851  {
1852  _tmp_Ke = ken;
1853  _tmp_Ke *= scaling_factor;
1854  jacobian.add_matrix(_tmp_Ke, di, dn);
1855 
1856  _tmp_Ke = kne;
1857  _tmp_Ke *= scaling_factor;
1858  jacobian.add_matrix(_tmp_Ke, dn, di);
1859 
1860  _tmp_Ke = knn;
1861  _tmp_Ke *= scaling_factor;
1862  jacobian.add_matrix(_tmp_Ke, dn);
1863  }
1864  else
1865  {
1866  jacobian.add_matrix(ken, di, dn);
1867  jacobian.add_matrix(kne, dn, di);
1868  jacobian.add_matrix(knn, dn);
1869  }
1870 }
1871 
1872 void
1873 Assembly::addJacobianScalar(SparseMatrix<Number> & jacobian)
1874 {
1875  const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
1876  for (const auto & ivar : scalar_vars)
1877  for (const auto & jvar : scalar_vars)
1878  if ((*_cm)(ivar->number(), jvar->number()) != 0 &&
1879  _jacobian_block_used[ivar->number()][jvar->number()])
1880  addJacobianBlock(jacobian,
1881  jacobianBlock(ivar->number(), jvar->number()),
1882  ivar->dofIndices(),
1883  jvar->dofIndices(),
1884  ivar->scalingFactor());
1885 }
1886 
1887 void
1888 Assembly::addJacobianOffDiagScalar(SparseMatrix<Number> & jacobian, unsigned int ivar)
1889 {
1890  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
1892  for (const auto & var_j : vars)
1893  if ((*_cm)(var_i.number(), var_j->number()) != 0 &&
1894  _jacobian_block_used[var_i.number()][var_j->number()])
1895  addJacobianBlock(jacobian,
1896  jacobianBlock(var_i.number(), var_j->number()),
1897  var_i.dofIndices(),
1898  var_j->dofIndices(),
1899  var_i.scalingFactor());
1900 }
1901 
1902 void
1903 Assembly::cacheJacobianContribution(numeric_index_type i, numeric_index_type j, Real value)
1904 {
1907  _cached_jacobian_contribution_vals.push_back(value);
1908 }
1909 
1910 void
1911 Assembly::setCachedJacobianContributions(SparseMatrix<Number> & jacobian)
1912 {
1913  // First zero the rows (including the diagonals) to prepare for
1914  // setting the cached values.
1915  jacobian.zero_rows(_cached_jacobian_contribution_rows, 0.0);
1916 
1917  // TODO: Use SparseMatrix::set_values() for efficiency
1918  for (unsigned int i = 0; i < _cached_jacobian_contribution_vals.size(); ++i)
1919  jacobian.set(_cached_jacobian_contribution_rows[i],
1922 
1924 }
1925 
1926 void
1927 Assembly::zeroCachedJacobianContributions(SparseMatrix<Number> & jacobian)
1928 {
1929  // First zero the rows (including the diagonals) to prepare for
1930  // setting the cached values.
1931  jacobian.zero_rows(_cached_jacobian_contribution_rows, 0.0);
1932 
1934 }
1935 
1936 void
1937 Assembly::addCachedJacobianContributions(SparseMatrix<Number> & jacobian)
1938 {
1939  // TODO: Use SparseMatrix::add_values() for efficiency
1940  for (unsigned int i = 0; i < _cached_jacobian_contribution_vals.size(); ++i)
1941  jacobian.add(_cached_jacobian_contribution_rows[i],
1944 
1946 }
1947 
1948 void
1950 {
1951  unsigned int orig_size = _cached_jacobian_contribution_rows.size();
1952 
1956 
1957  // It's possible (though massively unlikely) that clear() will
1958  // change the capacity of the vectors, so let's be paranoid and
1959  // explicitly reserve() the same amount of memory to avoid multiple
1960  // push_back() induced allocations. We reserve 20% more than the
1961  // original size that was cached to account for variations in the
1962  // number of BCs assigned to each thread (for when the Jacobian
1963  // contributions are computed threaded).
1964  _cached_jacobian_contribution_rows.reserve(1.2 * orig_size);
1965  _cached_jacobian_contribution_cols.reserve(1.2 * orig_size);
1966  _cached_jacobian_contribution_vals.reserve(1.2 * orig_size);
1967 }
1968 
1969 void
1971 {
1972  mooseAssert(_xfem != NULL, "This function should not be called if xfem is inactive");
1973 
1975  return;
1976 
1977  MooseArray<Real> xfem_weight_multipliers;
1978  if (_xfem->getXFEMWeights(xfem_weight_multipliers, elem, _current_qrule, _current_q_points))
1979  {
1980  mooseAssert(xfem_weight_multipliers.size() == _current_JxW.size(),
1981  "Size of weight multipliers in xfem doesn't match number of quadrature points");
1982  for (unsigned i = 0; i < xfem_weight_multipliers.size(); i++)
1983  {
1984  _current_JxW[i] = _current_JxW[i] * xfem_weight_multipliers[i];
1985  }
1986  xfem_weight_multipliers.release();
1987  }
1988 }
VariablePhiSecond _second_phi_face
Definition: Assembly.h:843
const VariablePhiValue & phiFace()
VariablePhiSecond _second_phi_face_neighbor
Definition: Assembly.h:851
Ok - here&#39;s the design.
Definition: Assembly.h:870
const std::vector< dof_id_type > & allDofIndices() const
Get all global dofindices for the variable.
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face
types of finite elements
Definition: Assembly.h:735
bool _need_neighbor_elem_volume
true is apps need to compute neighbor element volume
Definition: Assembly.h:800
ArbitraryQuadrature * _current_qrule_arbitrary
The current arbitrary quadrature rule used within the element interior.
Definition: Assembly.h:714
std::map< FEType, FEBase * > _current_fe
The "volume" fe object that matches the current elem.
Definition: Assembly.h:693
void addCachedResidual(NumericVector< Number > &residual, Moose::KernelType type)
Adds the values that have been cached by calling cacheResidual() and or cacheResidualNeighbor() to th...
Definition: Assembly.C:1426
void buildNeighborFE(FEType type)
Build FEs for a neighbor with a type.
Definition: Assembly.C:208
bool usesPhi()
Whether or not this variable is actually using the shape function value.
const VariablePhiSecond & feSecondPhiFace(FEType type)
Definition: Assembly.C:301
std::map< unsigned int, FEBase ** > _holder_fe_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:706
SystemBase & _sys
Definition: Assembly.h:667
const VariablePhiSecond & secondPhiFace()
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:905
unsigned int _max_cached_residuals
Definition: Assembly.h:907
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:1780
MooseArray< Point > _current_physical_points
This will be filled up with the physical points passed into reinitAtPhysical() if it is called...
Definition: Assembly.h:813
void addResidual(NumericVector< Number > &residual, Moose::KernelType type=Moose::KT_NONTIME)
Definition: Assembly.C:1330
std::vector< std::vector< unsigned char > > _jacobian_block_nonlocal_used
Definition: Assembly.h:677
void setResidualNeighbor(NumericVector< Number > &residual, Moose::KernelType type=Moose::KT_NONTIME)
Definition: Assembly.C:1488
void cacheJacobianBlock(DenseMatrix< Number > &jac_block, std::vector< dof_id_type > &idof_indices, std::vector< dof_id_type > &jdof_indices, Real scaling_factor)
Definition: Assembly.C:1522
MooseArray< Real > _coord_neighbor
The current coordinate transformation coefficients.
Definition: Assembly.h:777
const std::vector< MooseVariableScalar * > & getScalarVariables(THREAD_ID tid)
Definition: SystemBase.h:422
QBase * _current_qrule
The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac k...
Definition: Assembly.h:710
std::map< unsigned int, FEBase ** > _holder_fe_face_neighbor_helper
Definition: Assembly.h:768
virtual Moose::CoordinateSystemType getCoordSystem(SubdomainID sid)=0
Assembly(SystemBase &sys, THREAD_ID tid)
Definition: Assembly.C:39
Class for stuff related to variables.
Definition: MooseVariable.h:43
FEGenericBase< Real > FEBase
Definition: Assembly.h:34
void init()
Deprecated init method.
MooseArray< Real > _current_JxW_neighbor
The current transformed jacobian weights on a neighbor&#39;s face.
Definition: Assembly.h:775
std::shared_ptr< XFEMInterface > _xfem
The XFEM controller.
Definition: Assembly.h:690
virtual ~Assembly()
Definition: Assembly.C:113
void cacheResidualBlock(std::vector< Real > &cached_residual_values, std::vector< dof_id_type > &cached_residual_rows, DenseVector< Number > &res_block, std::vector< dof_id_type > &dof_indices, Real scaling_factor)
Definition: Assembly.C:1294
MooseArray< std::vector< Real > > _phi
Definition: Assembly.h:856
void cacheResidualContribution(dof_id_type dof, Real value, Moose::KernelType type)
Cache individual residual contributions.
Definition: Assembly.C:1375
const VariablePhiGradient & feGradPhiFaceNeighbor(FEType type)
Definition: Assembly.C:338
MooseArray< Real > _coord
The current coordinate transformation coefficients.
Definition: Assembly.h:722
void resize(const unsigned int size)
Change the number of elements the array can store.
Definition: MooseArray.h:205
const VariablePhiValue & phiNeighbor()
void reinitFE(const Elem *elem)
Just an internal helper function to reinit the volume FE objects.
Definition: Assembly.C:417
void prepareNonlocal()
Definition: Assembly.C:1059
VariablePhiGradient _grad_phi
Definition: Assembly.h:838
VariablePhiSecond _second_phi_neighbor
Definition: Assembly.h:847
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:734
std::map< unsigned int, ArbitraryQuadrature * > _holder_qrule_neighbor
Holds arbitrary qrules for each dimension.
Definition: Assembly.h:773
void addJacobianNonlocal(SparseMatrix< Number > &jacobian)
Definition: Assembly.C:1646
bool _current_elem_volume_computed
Boolean to indicate whether current element volumes has been computed.
Definition: Assembly.h:808
const VariablePhiValue & fePhiFace(FEType type)
Definition: Assembly.C:287
bool _currently_fe_caching
Whether or not fe should currently be cached - This will be false if something funky is going on with...
Definition: Assembly.h:893
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
Real _current_neighbor_volume
Volume of the current neighbor.
Definition: Assembly.h:802
std::map< FEType, FEBase * > _current_fe_face
The "face" fe object that matches the current elem.
Definition: Assembly.h:695
std::vector< std::pair< MooseVariable *, MooseVariable * > > _cm_entry
Entries in the coupling matrix (only for field variables)
Definition: Assembly.h:673
const VariablePhiValue & fePhiFaceNeighbor(FEType type)
Definition: Assembly.C:331
const Elem * _current_neighbor_elem
The current neighbor "element".
Definition: Assembly.h:792
std::vector< std::vector< DenseVector< Number > > > _sub_Rn
residual contributions for each variable from the neighbor
Definition: Assembly.h:818
const VariablePhiSecond & feSecondPhiFaceNeighbor(FEType type)
Definition: Assembly.C:345
const VariablePhiSecond & feSecondPhi(FEType type)
Definition: Assembly.C:279
void addJacobian(SparseMatrix< Number > &jacobian)
Definition: Assembly.C:1601
void setPoints(const std::vector< Point > &points)
MooseArray< Real > _current_JxW_face
The current transformed jacobian weights on a face.
Definition: Assembly.h:747
void modifyWeightsDueToXFEM(const Elem *elem)
Update the integration weights for XFEM partial elements.
Definition: Assembly.C:1970
bool _invalidated
Whether or not this data is invalid (needs to be recached) note that there is no constructor so the v...
Definition: Assembly.h:877
DenseMatrix< Number > _tmp_Ke
auxiliary matrix for scaling jacobians (optimization to avoid expensive construction/destruction) ...
Definition: Assembly.h:834
const Elem * _current_elem
The current "element" we are currently on.
Definition: Assembly.h:780
THREAD_ID _tid
Thread number (id)
Definition: Assembly.h:683
const Elem *& neighbor()
Return the neighbor element.
Definition: Assembly.h:235
const VariablePhiGradient & feGradPhiNeighbor(FEType type)
Definition: Assembly.C:316
void reinitFEFaceNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:540
const MooseArray< Real > & JxW()
Returns the reference to the transformed jacobian weights.
Definition: Assembly.h:147
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:844
const VariablePhiGradient & feGradPhi(FEType type)
Definition: Assembly.C:272
void prepareNeighbor()
Definition: Assembly.C:1114
const Real & neighborVolume()
Returns the reference to the current neighbor volume.
Definition: Assembly.C:353
std::map< FEType, bool > _need_second_derivative
Definition: Assembly.h:571
std::vector< numeric_index_type > _cached_jacobian_contribution_cols
Definition: Assembly.h:932
Real _current_elem_volume
Volume of the current element.
Definition: Assembly.h:784
std::map< FEType, FEBase * > _current_fe_face_neighbor
The "neighbor face" fe object that matches the current elem.
Definition: Assembly.h:699
void shallowCopy(const MooseArray &rhs)
Doesn&#39;t actually make a copy of the data.
Definition: MooseArray.h:284
std::map< FEType, FEShapeData * > _fe_shape_data_face_neighbor
Definition: Assembly.h:899
void reinitNodeNeighbor(const Node *node)
Reinitialize assembly data for a neighbor node.
Definition: Assembly.C:818
bool usesSecondPhi()
Whether or not this variable is actually using the shape function second derivative.
Base class for a system (of equations)
Definition: SystemBase.h:91
void setVolumeQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for volume integration.
Definition: Assembly.C:380
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:919
void buildFaceFE(FEType type)
Build FEs for a face with a type.
Definition: Assembly.C:190
VariablePhiValue _phi
Definition: Assembly.h:837
MooseArray< Point > _q_points
Cached xyz positions of quadrature points.
Definition: Assembly.h:883
void prepareScalar()
Definition: Assembly.C:1176
const CouplingMatrix * _cm
Coupling matrices.
Definition: Assembly.h:670
std::vector< std::vector< DenseMatrix< Number > > > _sub_Kee
jacobian contributions
Definition: Assembly.h:823
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:360
DenseVector< Number > _tmp_Re
auxiliary vector for scaling residuals (optimization to avoid expensive construction/destruction) ...
Definition: Assembly.h:820
unsigned int _mesh_dimension
Definition: Assembly.h:687
FEBase *& getFE(FEType type, unsigned int dim)
Get a reference to a pointer that will contain the current volume FE.
Definition: Assembly.C:244
void reinit(const Elem *elem)
Reinitialize objects (JxW, q_points, ...) for an elements.
Definition: Assembly.C:650
void prepareVariable(MooseVariable *var)
Used for preparing the dense residual and jacobian blocks for one particular variable.
Definition: Assembly.C:1076
QBase * _current_qrule_volume
The current volumetric quadrature for the element.
Definition: Assembly.h:712
const DofMap & _dof_map
DOF map.
Definition: Assembly.h:681
MooseArray< std::vector< RealGradient > > _grad_phi
Definition: Assembly.h:857
const VariablePhiSecond & secondPhiFaceNeighbor()
void copyNeighborShapes(unsigned int var)
Definition: Assembly.C:1250
unsigned int _max_cached_jacobians
Definition: Assembly.h:916
void cacheResidual()
Takes the values that are currently in _sub_Re and appends them to the cached values.
Definition: Assembly.C:1361
void copyShapes(unsigned int var)
Definition: Assembly.C:1226
std::vector< std::vector< DenseMatrix< Number > > > _sub_Kne
jacobian contributions from the neighbor and element
Definition: Assembly.h:829
std::map< dof_id_type, ElementFEShapeData * > _element_fe_shape_data_cache
Cached shape function values stored by element.
Definition: Assembly.h:887
void setResidual(NumericVector< Number > &residual, Moose::KernelType type=Moose::KT_NONTIME)
Definition: Assembly.C:1478
const CouplingMatrix & _nonlocal_cm
Definition: Assembly.h:671
std::vector< Point > _temp_reference_points
Temporary work data for reinitAtPhysical()
Definition: Assembly.h:925
unsigned int _current_neighbor_side
The current side of the selected neighboring element (valid only when working with sides) ...
Definition: Assembly.h:796
unsigned int & side()
Returns the current side.
Definition: Assembly.h:211
void computeCurrentFaceVolume()
Definition: Assembly.C:719
const Node *& node()
Returns the reference to the node.
Definition: Assembly.h:269
std::map< unsigned int, QBase * > _holder_qrule_face
Holds face qrules for each dimension.
Definition: Assembly.h:751
void prepareOffDiagScalar()
Definition: Assembly.C:1201
Implements a fake quadrature rule where you can specify the locations (in the reference domain) of th...
virtual MooseVariable & getVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:103
QBase * _current_qrule_neighbor
quadrature rule used on neighbors
Definition: Assembly.h:771
std::map< FEType, FEShapeData * > _fe_shape_data_neighbor
Definition: Assembly.h:898
const VariablePhiSecond & secondPhi()
std::vector< dof_id_type > _cached_jacobian_cols
Column where the corresponding cached value should go.
Definition: Assembly.h:914
std::map< FEType, FEBase * > _current_fe_neighbor
The "neighbor" fe object that matches the current elem.
Definition: Assembly.h:697
void addCachedJacobian(SparseMatrix< Number > &jacobian)
Adds the values that have been cached by calling cacheJacobian() and or cacheJacobianNeighbor() to th...
Definition: Assembly.C:1576
SubdomainID _current_neighbor_subdomain_id
The current neighbor subdomain ID.
Definition: Assembly.h:794
std::map< unsigned int, FEBase ** > _holder_fe_neighbor_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:767
const VariablePhiValue & phi()
std::vector< std::vector< DenseVector< Number > > > _sub_Re
residual contributions for each variable from the element
Definition: Assembly.h:816
void reinitNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:591
void addResidualScalar(NumericVector< Number > &residual, Moose::KernelType type=Moose::KT_NONTIME)
Definition: Assembly.C:1350
const VariablePhiGradient & gradPhi()
void prepare()
Definition: Assembly.C:1034
std::vector< T > stdVector()
Extremely inefficient way to produce a std::vector from a MooseArray!
Definition: MooseArray.h:330
std::vector< dof_id_type > & dofIndicesNeighbor()
Get DOF indices for currently selected element.
DofMap & dof_map
void setFaceQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for face integration.
Definition: Assembly.C:392
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:824
std::vector< std::vector< DenseMatrix< Number > > > _sub_Knn
jacobian contributions from the neighbor
Definition: Assembly.h:831
std::vector< std::vector< DenseMatrix< Number > > > _sub_Ken
jacobian contributions from the element and neighbor
Definition: Assembly.h:827
std::vector< numeric_index_type > _cached_jacobian_contribution_rows
Definition: Assembly.h:931
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_neighbor
types of finite elements
Definition: Assembly.h:764
VariablePhiGradient _grad_phi_face_neighbor
Definition: Assembly.h:850
void buildFaceNeighborFE(FEType type)
Build FEs for a neighbor face with a type.
Definition: Assembly.C:226
MooseArray< Real > _current_JxW
The current list of transformed jacobian weights.
Definition: Assembly.h:718
virtual bool checkNonlocalCouplingRequirement()
Definition: SubProblem.h:62
std::vector< dof_id_type > _temp_dof_indices
Temporary work vector to keep from reallocating it.
Definition: Assembly.h:922
void buildFE(FEType type)
Build FEs with a type.
Definition: Assembly.C:168
void setResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, std::vector< dof_id_type > &dof_indices, Real scaling_factor)
Definition: Assembly.C:1456
void cacheResidualNeighbor()
Takes the values that are currently in _sub_Ke and appends them to the cached values.
Definition: Assembly.C:1382
std::map< unsigned int, ArbitraryQuadrature * > _holder_qrule_arbitrary
Holds arbitrary qrules for each dimension.
Definition: Assembly.h:726
const VariablePhiValue & fePhi(FEType type)
Definition: Assembly.C:265
virtual SubProblem & subproblem()
Definition: SystemBase.h:103
VariablePhiGradient _grad_phi_neighbor
Definition: Assembly.h:846
const std::vector< MooseVariable * > & getVariables(THREAD_ID tid)
Definition: SystemBase.h:418
VariablePhiValue _phi_face
Definition: Assembly.h:841
void addJacobianOffDiagScalar(SparseMatrix< Number > &jacobian, unsigned int ivar)
Definition: Assembly.C:1888
std::vector< Real > _cached_jacobian_values
Values cached by calling cacheJacobian()
Definition: Assembly.h:910
void cacheJacobianContribution(numeric_index_type i, numeric_index_type j, Real value)
Caches the Jacobian entry &#39;value&#39;, to eventually be added/set in the (i,j) location of the matrix...
Definition: Assembly.C:1903
KernelType
Definition: MooseTypes.h:162
std::map< unsigned int, FEBase ** > _holder_fe_face_helper
Each dimension&#39;s helper objects.
Definition: Assembly.h:737
VariablePhiValue _phi_face_neighbor
Definition: Assembly.h:849
std::map< unsigned int, std::map< FEType, FEBase * > > _fe
Each dimension&#39;s actual fe objects indexed on type.
Definition: Assembly.h:704
std::vector< dof_id_type > _cached_jacobian_rows
Row where the corresponding cached value should go.
Definition: Assembly.h:912
VariablePhiValue _phi_neighbor
Definition: Assembly.h:845
std::map< unsigned int, ArbitraryQuadrature * > _holder_qface_arbitrary
Holds arbitrary face qrules for each dimension.
Definition: Assembly.h:753
const VariablePhiGradient & feGradPhiFace(FEType type)
Definition: Assembly.C:294
const VariablePhiGradient & gradPhiFaceNeighbor()
void zeroCachedJacobianContributions(SparseMatrix< Number > &jacobian)
Zero out previously-cached Jacobian rows.
Definition: Assembly.C:1927
Real _current_side_volume
Volume of the current side element.
Definition: Assembly.h:790
CoordinateSystemType
Definition: MooseTypes.h:212
std::vector< std::vector< Real > > _cached_residual_values
Values cached by calling cacheResidual() (the first vector is for TIME vs NONTIME) ...
Definition: Assembly.h:902
const Node * _current_node
The current node we are working with.
Definition: Assembly.h:804
void addCachedJacobianContributions(SparseMatrix< Number > &jacobian)
Adds previously-cached Jacobian values via SparseMatrix::add() calls.
Definition: Assembly.C:1937
DGJacobianType
Definition: MooseTypes.h:190
std::vector< Real > _cached_jacobian_contribution_vals
Storage for cached Jacobian entries.
Definition: Assembly.h:930
FEBase *& getFEFaceNeighbor(FEType type, unsigned int dim)
Get a reference to a pointer that will contain the current "neighbor" FE.
Definition: Assembly.C:258
MooseArray< Point > _current_normals
The current Normal vectors at the quadrature points.
Definition: Assembly.h:749
MatType type
unsigned int number() const
Get variable number coming from libMesh.
void cacheJacobian()
Takes the values that are currently in _sub_Kee and appends them to the cached values.
Definition: Assembly.C:1693
const VariablePhiValue & fePhiNeighbor(FEType type)
Definition: Assembly.C:309
virtual NumericVector< Number > & residualVector(Moose::KernelType)
Definition: SystemBase.h:158
void invalidateCache()
Invalidate any currently cached data.
Definition: Assembly.C:410
std::vector< std::pair< MooseVariable *, MooseVariable * > > _cm_nonlocal_entry
Definition: Assembly.h:674
const Elem *& elem()
Return the current element.
Definition: Assembly.h:189
MooseArray< std::vector< RealTensor > > _second_phi
Definition: Assembly.h:858
bool computingSecond()
Whether or not this variable is computing any second derivatives.
std::vector< dof_id_type > & dofIndices()
void setNeighborQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for neighbor integration.
Definition: Assembly.C:401
MooseArray< Real > _JxW
Cached JxW.
Definition: Assembly.h:880
const VariablePhiSecond & secondPhiNeighbor()
void prepareVariableNonlocal(MooseVariable *var)
Definition: Assembly.C:1098
unsigned int size() const
The number of elements that can currently be stored in the array.
Definition: MooseArray.h:250
void addResidualNeighbor(NumericVector< Number > &residual, Moose::KernelType type=Moose::KT_NONTIME)
Definition: Assembly.C:1340
void release()
Manually deallocates the data pointer.
Definition: MooseArray.h:56
const VariablePhiSecond & feSecondPhiNeighbor(FEType type)
Definition: Assembly.C:323
void reinitFENeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Definition: Assembly.C:566
void initNonlocalCoupling()
Create pair of variables requiring nonlocal jacobian contributions.
Definition: Assembly.C:1017
void setCoordinateTransformation(const QBase *qrule, const MooseArray< Point > &q_points)
Definition: Assembly.C:674
void cacheJacobianNeighbor()
Takes the values that are currently in the neighbor Dense Matrices and appends them to the cached val...
Definition: Assembly.C:1753
std::map< FEType, FEShapeData * > _shape_data
This is where the cached shape functions will be held.
Definition: Assembly.h:874
Class for scalar variables (they are different).
const VariablePhiGradient & gradPhiFace()
void cacheResidualNodes(DenseVector< Number > &res, std::vector< dof_id_type > &dof_index)
Lets an external class cache residual at a set of nodes.
Definition: Assembly.C:1396
bool _should_use_fe_cache
Whether or not fe cache should be built at all.
Definition: Assembly.h:890
unsigned int _current_side
The current side of the selected element (valid only when working with sides)
Definition: Assembly.h:786
const Elem * _current_neighbor_side_elem
The current side element of the ncurrent neighbor element.
Definition: Assembly.h:798
std::map< unsigned int, QBase * > _holder_qrule_volume
Holds volume qrules for each dimension.
Definition: Assembly.h:724
std::vector< std::vector< DenseMatrix< Number > > > _sub_Keg
Definition: Assembly.h:824
bool usesGradPhi()
Whether or not this variable is actually using the shape function gradient.
bool _current_side_volume_computed
Boolean to indicate whether current element side volumes has been computed.
Definition: Assembly.h:810
FEBase *& getFEFace(FEType type, unsigned int dim)
Get a reference to a pointer that will contain the current "face" FE.
Definition: Assembly.C:251
std::vector< std::vector< unsigned char > > _jacobian_block_used
Flag that indicates if the jacobian block was used.
Definition: Assembly.h:676
virtual bool hasResidualVector(Moose::KernelType) const
Definition: SystemBase.h:159
void prepareBlock(unsigned int ivar, unsigned jvar, const std::vector< dof_id_type > &dof_indices)
Definition: Assembly.C:1149
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face_neighbor
Definition: Assembly.h:765
Definition: Moose.h:84
void addJacobianScalar(SparseMatrix< Number > &jacobian)
Definition: Assembly.C:1873
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)
Definition: Assembly.C:1548
void computeCurrentElemVolume()
Definition: Assembly.C:704
std::map< FEType, FEShapeData * > _fe_shape_data
Definition: Assembly.h:896
DenseMatrix< Number > & jacobianBlock(unsigned int ivar, unsigned int jvar)
Definition: Assembly.C:887
MooseArray< Point > _current_q_points
The current list of quadrature points.
Definition: Assembly.h:716
Moose::CoordinateSystemType _coord_type
The coordinate system.
Definition: Assembly.h:720
const VariablePhiValue & phiFaceNeighbor()
VariablePhiSecond _second_phi
Definition: Assembly.h:839
std::map< FEType, FEShapeData * > _fe_shape_data_face
Definition: Assembly.h:897
void clearCachedJacobianContributions()
Clear any currently cached jacobian contributions.
Definition: Assembly.C:1949
MooseArray< Point > _current_q_points_face
The current quadrature points on a face.
Definition: Assembly.h:745
std::vector< std::vector< unsigned char > > _jacobian_block_neighbor_used
Flag that indicates if the jacobian block for neighbor was used.
Definition: Assembly.h:679
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:1165
const Elem * _current_side_elem
The current "element" making up the side we are currently on.
Definition: Assembly.h:788
void copyFaceShapes(unsigned int var)
Definition: Assembly.C:1238
DenseMatrix< Number > & jacobianBlockNeighbor(Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar)
Definition: Assembly.C:901
void cacheJacobianNonlocal()
Takes the values that are currently in _sub_Keg and appends them to the cached values.
Definition: Assembly.C:1739
void addJacobianNeighbor(SparseMatrix< Number > &jacobian)
Definition: Assembly.C:1661
VariablePhiGradient _grad_phi_face
Definition: Assembly.h:842
virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a scalar variable with specified number.
Definition: SystemBase.C:121
void addCachedResiduals()
Definition: Assembly.C:1410
DenseMatrix< Number > & jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar)
Definition: Assembly.C:894
virtual unsigned int nVariables()
Get the number of variables in this system.
Definition: SystemBase.C:580
SubdomainID _current_subdomain_id
The current subdomain ID.
Definition: Assembly.h:782
void addResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, Real scaling_factor)
Definition: Assembly.C:1270
QBase * _current_qrule_face
quadrature rule used on faces
Definition: Assembly.h:741
unsigned int THREAD_ID
Definition: MooseTypes.h:79
const Node * _current_neighbor_node
The current neighboring node we are working with.
Definition: Assembly.h:806
void setCachedJacobianContributions(SparseMatrix< Number > &jacobian)
Sets previously-cached Jacobian values via SparseMatrix::set() calls.
Definition: Assembly.C:1911
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:1804
void reinitFEFace(const Elem *elem, unsigned int side)
Just an internal helper function to reinit the face FE objects.
Definition: Assembly.C:509
const VariablePhiGradient & gradPhiNeighbor()
void scalingFactor(Real factor)
Set the scaling factor for this variable.
unsigned int getAxisymmetricRadialCoord()
Returns the desired radial direction for RZ coordinate transformation.
Definition: SubProblem.C:402