libMesh
fem_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #include "libmesh/dof_map.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/equation_systems.h"
23 #include "libmesh/fe_base.h"
24 #include "libmesh/fem_context.h"
25 #include "libmesh/fem_system.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/numeric_vector.h"
29 #include "libmesh/parallel.h"
30 #include "libmesh/parallel_algebra.h"
31 #include "libmesh/parallel_ghost_sync.h"
32 #include "libmesh/quadrature.h"
33 #include "libmesh/sparse_matrix.h"
34 #include "libmesh/time_solver.h"
35 #include "libmesh/unsteady_solver.h" // For eulerian_residual
36 #include "libmesh/fe_interface.h"
37 
38 namespace {
39 using namespace libMesh;
40 
41 // give this guy some scope since there
42 // is underlying vector allocation upon
43 // creation/deletion
44 ConstElemRange elem_range;
45 
46 typedef Threads::spin_mutex femsystem_mutex;
47 femsystem_mutex assembly_mutex;
48 
49 void assemble_unconstrained_element_system(const FEMSystem & _sys,
50  const bool _get_jacobian,
51  const bool _constrain_heterogeneously,
52  FEMContext & _femcontext)
53 {
54  if (_sys.print_element_solutions)
55  {
56  std::streamsize old_precision = libMesh::out.precision();
58  if (_femcontext.has_elem())
59  libMesh::out << "U_elem " << _femcontext.get_elem().id();
60  else
61  libMesh::out << "U_scalar ";
62  libMesh::out << " = " << _femcontext.get_elem_solution() << std::endl;
63 
64  if (_sys.use_fixed_solution)
65  {
66  if (_femcontext.has_elem())
67  libMesh::out << "Ufixed_elem " << _femcontext.get_elem().id();
68  else
69  libMesh::out << "Ufixed_scalar ";
70  libMesh::out << " = " << _femcontext.get_elem_fixed_solution() << std::endl;
71  libMesh::out.precision(old_precision);
72  }
73  }
74 
75  // We need jacobians to do heterogeneous residual constraints
76  const bool need_jacobian =
77  (_get_jacobian || _constrain_heterogeneously);
78 
79  bool jacobian_computed =
80  _sys.time_solver->element_residual(need_jacobian, _femcontext);
81 
82  // Compute a numeric jacobian if we have to
83  if (need_jacobian && !jacobian_computed)
84  {
85  // Make sure we didn't compute a jacobian and lie about it
86  libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
87  // Logging of numerical jacobians is done separately
88  _sys.numerical_elem_jacobian(_femcontext);
89  }
90 
91  // Compute a numeric jacobian if we're asked to verify the
92  // analytic jacobian we got
93  if (need_jacobian && jacobian_computed &&
94  _sys.verify_analytic_jacobians != 0.0)
95  {
96  DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
97 
98  _femcontext.get_elem_jacobian().zero();
99  // Logging of numerical jacobians is done separately
100  _sys.numerical_elem_jacobian(_femcontext);
101 
102  Real analytic_norm = analytic_jacobian.l1_norm();
103  Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
104 
105  // If we can continue, we'll probably prefer the analytic jacobian
106  analytic_jacobian.swap(_femcontext.get_elem_jacobian());
107 
108  // The matrix "analytic_jacobian" will now hold the error matrix
109  analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
110  Real error_norm = analytic_jacobian.l1_norm();
111 
112  Real relative_error = error_norm /
113  std::max(analytic_norm, numerical_norm);
114 
115  if (relative_error > _sys.verify_analytic_jacobians)
116  {
117  libMesh::err << "Relative error " << relative_error
118  << " detected in analytic jacobian on element "
119  << _femcontext.get_elem().id() << '!' << std::endl;
120 
121  std::streamsize old_precision = libMesh::out.precision();
123  libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
124  << _femcontext.get_elem_jacobian() << std::endl;
125  analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
126  libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
127  << analytic_jacobian << std::endl;
128 
129  libMesh::out.precision(old_precision);
130 
131  libmesh_error_msg("Relative error too large, exiting!");
132  }
133  }
134 
135  for (_femcontext.side = 0;
136  _femcontext.side != _femcontext.get_elem().n_sides();
137  ++_femcontext.side)
138  {
139  // Don't compute on non-boundary sides unless requested
140  if (!_sys.get_physics()->compute_internal_sides &&
141  _femcontext.get_elem().neighbor_ptr(_femcontext.side) != libmesh_nullptr)
142  continue;
143 
144  // Any mesh movement has already been done (and restored,
145  // if the TimeSolver isn't broken), but
146  // reinitializing the side FE objects is still necessary
147  _femcontext.side_fe_reinit();
148 
149  DenseMatrix<Number> old_jacobian;
150  // If we're in DEBUG mode, we should always verify that the
151  // user's side_residual function doesn't alter our existing
152  // jacobian and then lie about it
153 #ifndef DEBUG
154  // Even if we're not in DEBUG mode, when we're verifying
155  // analytic jacobians we'll want to verify each side's
156  // jacobian contribution separately.
157  if (_sys.verify_analytic_jacobians != 0.0 && need_jacobian)
158 #endif // ifndef DEBUG
159  {
160  old_jacobian = _femcontext.get_elem_jacobian();
161  _femcontext.get_elem_jacobian().zero();
162  }
163 
164  jacobian_computed =
165  _sys.time_solver->side_residual(need_jacobian, _femcontext);
166 
167  // Compute a numeric jacobian if we have to
168  if (need_jacobian && !jacobian_computed)
169  {
170  // If we have already backed up old_jacobian,
171  // we can make sure side_residual didn't compute a
172  // jacobian and lie about it.
173  //
174  // If we haven't, then we need to, to let
175  // numerical_side_jacobian work.
176  if (old_jacobian.m())
177  libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
178  else
179  {
180  old_jacobian = _femcontext.get_elem_jacobian();
181  _femcontext.get_elem_jacobian().zero();
182  }
183 
184  // Logging of numerical jacobians is done separately
185  _sys.numerical_side_jacobian(_femcontext);
186 
187  // Add back in element interior numerical Jacobian
188  _femcontext.get_elem_jacobian() += old_jacobian;
189  }
190 
191  // Compute a numeric jacobian if we're asked to verify the
192  // analytic jacobian we got
193  else if (need_jacobian && jacobian_computed &&
194  _sys.verify_analytic_jacobians != 0.0)
195  {
196  DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
197 
198  _femcontext.get_elem_jacobian().zero();
199  // Logging of numerical jacobians is done separately
200  _sys.numerical_side_jacobian(_femcontext);
201 
202  Real analytic_norm = analytic_jacobian.l1_norm();
203  Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
204 
205  // If we can continue, we'll probably prefer the analytic jacobian
206  analytic_jacobian.swap(_femcontext.get_elem_jacobian());
207 
208  // The matrix "analytic_jacobian" will now hold the error matrix
209  analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
210  Real error_norm = analytic_jacobian.l1_norm();
211 
212  Real relative_error = error_norm /
213  std::max(analytic_norm, numerical_norm);
214 
215  if (relative_error > _sys.verify_analytic_jacobians)
216  {
217  libMesh::err << "Relative error " << relative_error
218  << " detected in analytic jacobian on element "
219  << _femcontext.get_elem().id()
220  << ", side "
221  << static_cast<unsigned int>(_femcontext.side) << '!' << std::endl;
222 
223  std::streamsize old_precision = libMesh::out.precision();
225  libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
226  << _femcontext.get_elem_jacobian() << std::endl;
227  analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
228  libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
229  << analytic_jacobian << std::endl;
230  libMesh::out.precision(old_precision);
231 
232  libmesh_error_msg("Relative error too large, exiting!");
233  }
234  // Once we've verified a side, we'll want to add back the
235  // rest of the accumulated jacobian
236  _femcontext.get_elem_jacobian() += old_jacobian;
237  }
238 
239  // In DEBUG mode, we've set elem_jacobian == 0, and we
240  // may have yet to add the old jacobian back
241 #ifdef DEBUG
242  else
243  {
244  _femcontext.get_elem_jacobian() += old_jacobian;
245  }
246 #endif // ifdef DEBUG
247  }
248 }
249 
250 void add_element_system(const FEMSystem & _sys,
251  const bool _get_residual,
252  const bool _get_jacobian,
253  const bool _constrain_heterogeneously,
254  const bool _no_constraints,
255  FEMContext & _femcontext)
256 {
257 #ifdef LIBMESH_ENABLE_CONSTRAINTS
258  if (_get_residual && _sys.print_element_residuals)
259  {
260  std::streamsize old_precision = libMesh::out.precision();
262  if (_femcontext.has_elem())
263  libMesh::out << "Rraw_elem " << _femcontext.get_elem().id();
264  else
265  libMesh::out << "Rraw_scalar ";
266  libMesh::out << " = " << _femcontext.get_elem_residual() << std::endl;
267  libMesh::out.precision(old_precision);
268  }
269 
270  // We turn off the asymmetric constraint application;
271  // enforce_constraints_exactly() should be called in the solver
272  if (_get_residual && _get_jacobian)
273  {
274  if (_constrain_heterogeneously)
276  (_femcontext.get_elem_jacobian(),
277  _femcontext.get_elem_residual(),
278  _femcontext.get_dof_indices(), false);
279  else if (!_no_constraints)
281  (_femcontext.get_elem_jacobian(),
282  _femcontext.get_elem_residual(),
283  _femcontext.get_dof_indices(), false);
284  // Do nothing if (_no_constraints)
285  }
286  else if (_get_residual)
287  {
288  if (_constrain_heterogeneously)
290  (_femcontext.get_elem_jacobian(),
291  _femcontext.get_elem_residual(),
292  _femcontext.get_dof_indices(), false);
293  else if (!_no_constraints)
295  (_femcontext.get_elem_residual(), _femcontext.get_dof_indices(), false);
296  // Do nothing if (_no_constraints)
297  }
298  else if (_get_jacobian)
299  {
300  // Heterogeneous and homogeneous constraints are the same on the
301  // matrix
302  // Only get these contribs if we are applying some constraints
303  if (!_no_constraints)
305  _femcontext.get_dof_indices(),
306  false);
307  }
308 #endif // #ifdef LIBMESH_ENABLE_CONSTRAINTS
309 
310  if (_get_residual && _sys.print_element_residuals)
311  {
312  std::streamsize old_precision = libMesh::out.precision();
314  if (_femcontext.has_elem())
315  libMesh::out << "R_elem " << _femcontext.get_elem().id();
316  else
317  libMesh::out << "R_scalar ";
318  libMesh::out << " = " << _femcontext.get_elem_residual() << std::endl;
319  libMesh::out.precision(old_precision);
320  }
321 
322  if (_get_jacobian && _sys.print_element_jacobians)
323  {
324  std::streamsize old_precision = libMesh::out.precision();
326  if (_femcontext.has_elem())
327  libMesh::out << "J_elem " << _femcontext.get_elem().id();
328  else
329  libMesh::out << "J_scalar ";
330  libMesh::out << " = " << _femcontext.get_elem_jacobian() << std::endl;
331  libMesh::out.precision(old_precision);
332  }
333 
334  { // A lock is necessary around access to the global system
335  femsystem_mutex::scoped_lock lock(assembly_mutex);
336 
337  if (_get_jacobian)
338  _sys.matrix->add_matrix (_femcontext.get_elem_jacobian(),
339  _femcontext.get_dof_indices());
340  if (_get_residual)
341  _sys.rhs->add_vector (_femcontext.get_elem_residual(),
342  _femcontext.get_dof_indices());
343  } // Scope for assembly mutex
344 }
345 
346 
347 
348 class AssemblyContributions
349 {
350 public:
354  AssemblyContributions(FEMSystem & sys,
355  bool get_residual,
356  bool get_jacobian,
357  bool constrain_heterogeneously,
358  bool no_constraints) :
359  _sys(sys),
360  _get_residual(get_residual),
361  _get_jacobian(get_jacobian),
362  _constrain_heterogeneously(constrain_heterogeneously),
363  _no_constraints(no_constraints) {}
364 
368  void operator()(const ConstElemRange & range) const
369  {
371  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
372  _sys.init_context(_femcontext);
373 
374  for (ConstElemRange::const_iterator elem_it = range.begin();
375  elem_it != range.end(); ++elem_it)
376  {
377  Elem * el = const_cast<Elem *>(*elem_it);
378 
379  _femcontext.pre_fe_reinit(_sys, el);
380  _femcontext.elem_fe_reinit();
381 
382  assemble_unconstrained_element_system
383  (_sys, _get_jacobian, _constrain_heterogeneously, _femcontext);
384 
385  add_element_system
386  (_sys, _get_residual, _get_jacobian,
387  _constrain_heterogeneously, _no_constraints, _femcontext);
388  }
389  }
390 
391 private:
392 
393  FEMSystem & _sys;
394 
395  const bool _get_residual, _get_jacobian, _constrain_heterogeneously, _no_constraints;
396 };
397 
398 class PostprocessContributions
399 {
400 public:
404  explicit
405  PostprocessContributions(FEMSystem & sys) : _sys(sys) {}
406 
410  void operator()(const ConstElemRange & range) const
411  {
413  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
414  _sys.init_context(_femcontext);
415 
416  for (ConstElemRange::const_iterator elem_it = range.begin();
417  elem_it != range.end(); ++elem_it)
418  {
419  Elem * el = const_cast<Elem *>(*elem_it);
420  _femcontext.pre_fe_reinit(_sys, el);
421 
422  // Optionally initialize all the interior FE objects on elem.
424  _femcontext.elem_fe_reinit();
425 
426  _sys.element_postprocess(_femcontext);
427 
428  for (_femcontext.side = 0;
429  _femcontext.side != _femcontext.get_elem().n_sides();
430  ++_femcontext.side)
431  {
432  // Don't compute on non-boundary sides unless requested
433  if (!_sys.postprocess_sides ||
435  _femcontext.get_elem().neighbor_ptr(_femcontext.side) != libmesh_nullptr))
436  continue;
437 
438  // Optionally initialize all the FE objects on this side.
440  _femcontext.side_fe_reinit();
441 
442  _sys.side_postprocess(_femcontext);
443  }
444  }
445  }
446 
447 private:
448 
449  FEMSystem & _sys;
450 };
451 
452 class QoIContributions
453 {
454 public:
458  explicit
459  QoIContributions(FEMSystem & sys,
460  DifferentiableQoI & diff_qoi,
461  const QoISet & qoi_indices) :
462  qoi(sys.qoi.size(), 0.), _sys(sys), _diff_qoi(diff_qoi),_qoi_indices(qoi_indices) {}
463 
467  QoIContributions(const QoIContributions & other,
468  Threads::split) :
469  qoi(other._sys.qoi.size(), 0.), _sys(other._sys), _diff_qoi(other._diff_qoi) {}
470 
474  void operator()(const ConstElemRange & range)
475  {
477  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
478  _diff_qoi.init_context(_femcontext);
479 
480  bool have_some_heterogenous_qoi_bc = false;
481 #ifdef LIBMESH_ENABLE_CONSTRAINTS
482  std::vector<bool> have_heterogenous_qoi_bc(_sys.qoi.size(), false);
483  for (std::size_t q=0; q != _sys.qoi.size(); ++q)
484  if (_qoi_indices.has_index(q) &&
486  {
487  have_heterogenous_qoi_bc[q] = true;
488  have_some_heterogenous_qoi_bc = true;
489  }
490 #endif
491 
492  if (have_some_heterogenous_qoi_bc)
493  _sys.init_context(_femcontext);
494 
495  for (ConstElemRange::const_iterator elem_it = range.begin();
496  elem_it != range.end(); ++elem_it)
497  {
498  Elem * el = const_cast<Elem *>(*elem_it);
499 
500  _femcontext.pre_fe_reinit(_sys, el);
501 
502  // We might have some heterogenous dofs here; let's see for
503  // certain
504 #ifdef LIBMESH_ENABLE_CONSTRAINTS
505  bool elem_has_some_heterogenous_qoi_bc = false;
506  std::vector<bool> elem_has_heterogenous_qoi_bc(_sys.qoi.size(), false);
507  if (have_some_heterogenous_qoi_bc)
508  {
509  for (std::size_t q=0; q != _sys.qoi.size(); ++q)
510  {
511  if (have_heterogenous_qoi_bc[q])
512  {
513  for (std::size_t d=0; d != _femcontext.get_dof_indices().size(); ++d)
515  (q, _femcontext.get_dof_indices()[d]) != Number(0))
516  {
517  elem_has_some_heterogenous_qoi_bc = true;
518  elem_has_heterogenous_qoi_bc[q] = true;
519  break;
520  }
521  }
522  }
523  }
524 #endif
525 
526  if (_diff_qoi.assemble_qoi_elements ||
527  elem_has_some_heterogenous_qoi_bc)
528  _femcontext.elem_fe_reinit();
529 
530  if (_diff_qoi.assemble_qoi_elements)
531  _diff_qoi.element_qoi(_femcontext, _qoi_indices);
532 
533  // If we have some heterogenous dofs here, those are
534  // themselves part of a regularized flux QoI which the library
535  // handles integrating
536 #ifdef LIBMESH_ENABLE_CONSTRAINTS
537  if (elem_has_some_heterogenous_qoi_bc)
538  {
539  _sys.time_solver->element_residual(false, _femcontext);
540 
541  for (std::size_t q=0; q != _sys.qoi.size(); ++q)
542  {
543  if (elem_has_heterogenous_qoi_bc[q])
544  {
545  for (std::size_t d=0; d != _femcontext.get_dof_indices().size(); ++d)
546  this->qoi[q] -= _femcontext.get_elem_residual()(d) *
548 
549  }
550  }
551  }
552 #endif
553 
554  for (_femcontext.side = 0;
555  _femcontext.side != _femcontext.get_elem().n_sides();
556  ++_femcontext.side)
557  {
558  // Don't compute on non-boundary sides unless requested
559  if (!_diff_qoi.assemble_qoi_sides ||
560  (!_diff_qoi.assemble_qoi_internal_sides &&
561  _femcontext.get_elem().neighbor_ptr(_femcontext.side) != libmesh_nullptr))
562  continue;
563 
564  _femcontext.side_fe_reinit();
565 
566  _diff_qoi.side_qoi(_femcontext, _qoi_indices);
567  }
568  }
569 
570  this->_diff_qoi.thread_join( this->qoi, _femcontext.get_qois(), _qoi_indices );
571  }
572 
573  void join (const QoIContributions & other)
574  {
575  libmesh_assert_equal_to (this->qoi.size(), other.qoi.size());
576  this->_diff_qoi.thread_join( this->qoi, other.qoi, _qoi_indices );
577  }
578 
579  std::vector<Number> qoi;
580 
581 private:
582 
583  FEMSystem & _sys;
584  DifferentiableQoI & _diff_qoi;
585 
586  const QoISet _qoi_indices;
587 };
588 
589 class QoIDerivativeContributions
590 {
591 public:
595  QoIDerivativeContributions(FEMSystem & sys,
596  const QoISet & qoi_indices,
597  DifferentiableQoI & qoi,
598  bool include_liftfunc,
599  bool apply_constraints) :
600  _sys(sys),
601  _qoi_indices(qoi_indices),
602  _qoi(qoi),
603  _include_liftfunc(include_liftfunc),
604  _apply_constraints(apply_constraints) {}
605 
609  void operator()(const ConstElemRange & range) const
610  {
612  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
613  _qoi.init_context(_femcontext);
614 
615  bool have_some_heterogenous_qoi_bc = false;
616 #ifdef LIBMESH_ENABLE_CONSTRAINTS
617  std::vector<bool> have_heterogenous_qoi_bc(_sys.qoi.size(), false);
618  if (_include_liftfunc || _apply_constraints)
619  for (std::size_t q=0; q != _sys.qoi.size(); ++q)
620  if (_qoi_indices.has_index(q) &&
622  {
623  have_heterogenous_qoi_bc[q] = true;
624  have_some_heterogenous_qoi_bc = true;
625  }
626 #endif
627 
628  if (have_some_heterogenous_qoi_bc)
629  _sys.init_context(_femcontext);
630 
631  for (ConstElemRange::const_iterator elem_it = range.begin();
632  elem_it != range.end(); ++elem_it)
633  {
634  Elem * el = const_cast<Elem *>(*elem_it);
635 
636  _femcontext.pre_fe_reinit(_sys, el);
637 
638  // We might have some heterogenous dofs here; let's see for
639  // certain
640 #ifdef LIBMESH_ENABLE_CONSTRAINTS
641  bool elem_has_some_heterogenous_qoi_bc = false;
642  std::vector<bool> elem_has_heterogenous_qoi_bc(_sys.qoi.size(), false);
643  if (have_some_heterogenous_qoi_bc)
644  {
645  for (std::size_t q=0; q != _sys.qoi.size(); ++q)
646  {
647  if (have_heterogenous_qoi_bc[q])
648  {
649  for (std::size_t d=0; d != _femcontext.get_dof_indices().size(); ++d)
651  (q, _femcontext.get_dof_indices()[d]) != Number(0))
652  {
653  elem_has_some_heterogenous_qoi_bc = true;
654  elem_has_heterogenous_qoi_bc[q] = true;
655  break;
656  }
657  }
658  }
659  }
660 #endif
661 
662  // If we're going to call a user integral, then we need FE
663  // information to call element_qoi.
664  // If we're going to evaluate lift-function-based components
665  // of a QoI, then we need FE information to assemble the
666  // element residual.
667  if (_qoi.assemble_qoi_elements ||
668  ((_include_liftfunc || _apply_constraints) &&
669  elem_has_some_heterogenous_qoi_bc))
670  _femcontext.elem_fe_reinit();
671 
672  if (_qoi.assemble_qoi_elements)
673  _qoi.element_qoi_derivative(_femcontext, _qoi_indices);
674 
675 #ifdef LIBMESH_ENABLE_CONSTRAINTS
676  // If we need to use heterogenous dofs here, we need the
677  // Jacobian either for the regularized flux QoI integration
678  // and/or for constraint application.
679  if ((_include_liftfunc || _apply_constraints) &&
680  elem_has_some_heterogenous_qoi_bc)
681  {
682  bool jacobian_computed = _sys.time_solver->element_residual(true, _femcontext);
683 
684  // If we're using numerical jacobians, above wont compute them
685  if (!jacobian_computed)
686  {
687  // Make sure we didn't compute a jacobian and lie about it
688  libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
689  // Logging of numerical jacobians is done separately
690  _sys.numerical_elem_jacobian(_femcontext);
691  }
692  }
693 
694  // If we have some heterogenous dofs here, those are
695  // themselves part of a regularized flux QoI which the library
696  // may handle integrating
697  if (_include_liftfunc && elem_has_some_heterogenous_qoi_bc)
698  {
699  for (std::size_t q=0; q != _sys.qoi.size(); ++q)
700  {
701  if (elem_has_heterogenous_qoi_bc[q])
702  {
703  for (std::size_t i=0; i != _femcontext.get_dof_indices().size(); ++i)
704  {
705  Number liftfunc_val =
707 
708  if (liftfunc_val != Number(0))
709  {
710  for (std::size_t j=0; j != _femcontext.get_dof_indices().size(); ++j)
711  _femcontext.get_qoi_derivatives()[q](j) -=
712  _femcontext.get_elem_jacobian()(i,j) *
713  liftfunc_val;
714  }
715  }
716  }
717  }
718  }
719 #endif
720 
721 
722  for (_femcontext.side = 0;
723  _femcontext.side != _femcontext.get_elem().n_sides();
724  ++_femcontext.side)
725  {
726  // Don't compute on non-boundary sides unless requested
727  if (!_qoi.assemble_qoi_sides ||
728  (!_qoi.assemble_qoi_internal_sides &&
729  _femcontext.get_elem().neighbor_ptr(_femcontext.side) != libmesh_nullptr))
730  continue;
731 
732  _femcontext.side_fe_reinit();
733 
734  _qoi.side_qoi_derivative(_femcontext, _qoi_indices);
735  }
736 
737  // We need some unmodified indices to use for constraining
738  // multiple vector
739  // FIXME - there should be a DofMap::constrain_element_vectors
740  // to do this more efficiently
741 #ifdef LIBMESH_ENABLE_CONSTRAINTS
742  std::vector<dof_id_type> original_dofs = _femcontext.get_dof_indices();
743 #endif
744 
745  { // A lock is necessary around access to the global system
746  femsystem_mutex::scoped_lock lock(assembly_mutex);
747 
748 #ifdef LIBMESH_ENABLE_CONSTRAINTS
749  // We'll need to see if any heterogenous constraints apply
750  // to the QoI dofs on this element *or* to any of the dofs
751  // they depend on, so let's get those dependencies
752  if (_apply_constraints)
753  _sys.get_dof_map().constrain_nothing(_femcontext.get_dof_indices());
754 #endif
755 
756  for (std::size_t i=0; i != _sys.qoi.size(); ++i)
757  if (_qoi_indices.has_index(i))
758  {
759 #ifdef LIBMESH_ENABLE_CONSTRAINTS
760  if (_apply_constraints)
761  {
762 #ifndef NDEBUG
763  bool has_heterogenous_constraint = false;
764  for (std::size_t d=0; d != _femcontext.get_dof_indices().size(); ++d)
766  (i, _femcontext.get_dof_indices()[d]) != Number(0))
767  {
768  has_heterogenous_constraint = true;
769  libmesh_assert(elem_has_heterogenous_qoi_bc[i]);
770  libmesh_assert(elem_has_some_heterogenous_qoi_bc);
771  break;
772  }
773 #else
774  bool has_heterogenous_constraint =
775  elem_has_heterogenous_qoi_bc[i];
776 #endif
777 
778  _femcontext.get_dof_indices() = original_dofs;
779 
780  if (has_heterogenous_constraint)
781  {
782  // Q_u gets used for *adjoint* solves, so we
783  // need K^T here.
784  DenseMatrix<Number> elem_jacobian_transpose;
785  _femcontext.get_elem_jacobian().get_transpose
786  (elem_jacobian_transpose);
787 
789  (elem_jacobian_transpose,
790  _femcontext.get_qoi_derivatives()[i],
791  _femcontext.get_dof_indices(), false, i);
792  }
793  else
794  {
796  (_femcontext.get_qoi_derivatives()[i],
797  _femcontext.get_dof_indices(), false);
798  }
799  }
800 #endif
801 
803  (_femcontext.get_qoi_derivatives()[i], _femcontext.get_dof_indices());
804  }
805  }
806  }
807  }
808 
809 private:
810 
811  FEMSystem & _sys;
812  const QoISet & _qoi_indices;
813  DifferentiableQoI & _qoi;
814  bool _include_liftfunc, _apply_constraints;
815 };
816 
817 
818 }
819 
820 
821 namespace libMesh
822 {
823 
824 
825 
826 
827 
829  const std::string & name_in,
830  const unsigned int number_in)
831  : Parent(es, name_in, number_in),
832  fe_reinit_during_postprocess(true),
833  numerical_jacobian_h(TOLERANCE),
834  verify_analytic_jacobians(0.0)
835 {
836 }
837 
838 
840 {
841 }
842 
843 
844 
846 {
847  // First initialize LinearImplicitSystem data
849 }
850 
851 
852 void FEMSystem::assembly (bool get_residual, bool get_jacobian,
853  bool apply_heterogeneous_constraints,
854  bool apply_no_constraints)
855 {
856  libmesh_assert(get_residual || get_jacobian);
857 
858 // Log residual and jacobian and combined performance separately
859 #ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
860  const char * log_name;
861  if (get_residual && get_jacobian)
862  log_name = "assembly()";
863  else if (get_residual)
864  log_name = "assembly(get_residual)";
865  else
866  log_name = "assembly(get_jacobian)";
867 
868  LOG_SCOPE(log_name, "FEMSystem");
869 #endif
870 
871  const MeshBase & mesh = this->get_mesh();
872 
873  // this->get_vector("_nonlinear_solution").localize
874  // (*current_local_nonlinear_solution,
875  // dof_map.get_send_list());
876  this->update();
877 
879  {
880  // this->get_vector("_nonlinear_solution").close();
881  this->solution->close();
882 
883  std::streamsize old_precision = libMesh::out.precision();
885  libMesh::out << "|U| = "
886  // << this->get_vector("_nonlinear_solution").l1_norm()
887  << this->solution->l1_norm()
888  << std::endl;
889  libMesh::out.precision(old_precision);
890  }
891  if (print_solutions)
892  {
893  std::streamsize old_precision = libMesh::out.precision();
895  // libMesh::out << "U = [" << this->get_vector("_nonlinear_solution")
896  libMesh::out << "U = [" << *(this->solution)
897  << "];" << std::endl;
898  libMesh::out.precision(old_precision);
899  }
900 
901  // Is this definitely necessary? [RHS]
902  // Yes. [RHS 2012]
903  if (get_jacobian)
904  matrix->zero();
905  if (get_residual)
906  rhs->zero();
907 
908  // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
909  if (verify_analytic_jacobians > 0.5)
910  {
911  libMesh::err << "WARNING! verify_analytic_jacobians was set "
912  << "to absurdly large value of "
913  << verify_analytic_jacobians << std::endl;
914  libMesh::err << "Resetting to 1e-6!" << std::endl;
916  }
917 
918  // In time-dependent problems, the nonlinear function we're trying
919  // to solve at each timestep may depend on the particular solver
920  // we're using
922 
923  // Build the residual and jacobian contributions on every active
924  // mesh element on this processor
926  (elem_range.reset(mesh.active_local_elements_begin(),
928  AssemblyContributions(*this, get_residual, get_jacobian,
929  apply_heterogeneous_constraints,
930  apply_no_constraints));
931 
932  // Check and see if we have SCALAR variables
933  bool have_scalar = false;
934  for (unsigned int i=0; i != this->n_variable_groups(); ++i)
935  {
936  if (this->variable_group(i).type().family == SCALAR)
937  {
938  have_scalar = true;
939  break;
940  }
941  }
942 
943  // SCALAR dofs are stored on the last processor, so we'll evaluate
944  // their equation terms there and only if we have a SCALAR variable
945  if (this->processor_id() == (this->n_processors()-1) && have_scalar)
946  {
948  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
949  this->init_context(_femcontext);
950  _femcontext.pre_fe_reinit(*this, libmesh_nullptr);
951 
952  bool jacobian_computed =
953  this->time_solver->nonlocal_residual(get_jacobian, _femcontext);
954 
955  // Nonlocal residuals are likely to be length 0, in which case we
956  // don't need to do any more. And we shouldn't try to do any
957  // more; lots of DenseVector/DenseMatrix code assumes rank>0.
958  if (_femcontext.get_elem_residual().size())
959  {
960  // Compute a numeric jacobian if we have to
961  if (get_jacobian && !jacobian_computed)
962  {
963  // Make sure we didn't compute a jacobian and lie about it
964  libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
965  // Logging of numerical jacobians is done separately
966  this->numerical_nonlocal_jacobian(_femcontext);
967  }
968 
969  // Compute a numeric jacobian if we're asked to verify the
970  // analytic jacobian we got
971  if (get_jacobian && jacobian_computed &&
972  this->verify_analytic_jacobians != 0.0)
973  {
974  DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
975 
976  _femcontext.get_elem_jacobian().zero();
977  // Logging of numerical jacobians is done separately
978  this->numerical_nonlocal_jacobian(_femcontext);
979 
980  Real analytic_norm = analytic_jacobian.l1_norm();
981  Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
982 
983  // If we can continue, we'll probably prefer the analytic jacobian
984  analytic_jacobian.swap(_femcontext.get_elem_jacobian());
985 
986  // The matrix "analytic_jacobian" will now hold the error matrix
987  analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
988  Real error_norm = analytic_jacobian.l1_norm();
989 
990  Real relative_error = error_norm /
991  std::max(analytic_norm, numerical_norm);
992 
993  if (relative_error > this->verify_analytic_jacobians)
994  {
995  libMesh::err << "Relative error " << relative_error
996  << " detected in analytic jacobian on nonlocal dofs!"
997  << std::endl;
998 
999  std::streamsize old_precision = libMesh::out.precision();
1000  libMesh::out.precision(16);
1001  libMesh::out << "J_analytic nonlocal = "
1002  << _femcontext.get_elem_jacobian() << std::endl;
1003  analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
1004  libMesh::out << "J_numeric nonlocal = "
1005  << analytic_jacobian << std::endl;
1006 
1007  libMesh::out.precision(old_precision);
1008 
1009  libmesh_error_msg("Relative error too large, exiting!");
1010  }
1011  }
1012 
1013  add_element_system
1014  (*this, get_residual, get_jacobian,
1015  apply_heterogeneous_constraints, apply_no_constraints, _femcontext);
1016  }
1017  }
1018 
1019  if (get_residual && (print_residual_norms || print_residuals))
1020  this->rhs->close();
1021  if (get_residual && print_residual_norms)
1022  {
1023  std::streamsize old_precision = libMesh::out.precision();
1024  libMesh::out.precision(16);
1025  libMesh::out << "|F| = " << this->rhs->l1_norm() << std::endl;
1026  libMesh::out.precision(old_precision);
1027  }
1028  if (get_residual && print_residuals)
1029  {
1030  std::streamsize old_precision = libMesh::out.precision();
1031  libMesh::out.precision(16);
1032  libMesh::out << "F = [" << *(this->rhs) << "];" << std::endl;
1033  libMesh::out.precision(old_precision);
1034  }
1035 
1036  if (get_jacobian && (print_jacobian_norms || print_jacobians))
1037  this->matrix->close();
1038  if (get_jacobian && print_jacobian_norms)
1039  {
1040  std::streamsize old_precision = libMesh::out.precision();
1041  libMesh::out.precision(16);
1042  libMesh::out << "|J| = " << this->matrix->l1_norm() << std::endl;
1043  libMesh::out.precision(old_precision);
1044  }
1045  if (get_jacobian && print_jacobians)
1046  {
1047  std::streamsize old_precision = libMesh::out.precision();
1048  libMesh::out.precision(16);
1049  libMesh::out << "J = [" << *(this->matrix) << "];" << std::endl;
1050  libMesh::out.precision(old_precision);
1051  }
1052 }
1053 
1054 
1055 
1057 {
1058  // We are solving the primal problem
1059  Parent::solve();
1060 
1061  // On a moving mesh we want the mesh to reflect the new solution
1062  this->mesh_position_set();
1063 }
1064 
1065 
1066 
1068 {
1069  // If we don't need to move the mesh, we're done
1070  if (_mesh_sys != this)
1071  return;
1072 
1073  MeshBase & mesh = this->get_mesh();
1074 
1075  UniquePtr<DiffContext> con = this->build_context();
1076  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
1077  this->init_context(_femcontext);
1078 
1079  // Move every mesh element we can
1080  for (const auto & elem : mesh.active_local_element_ptr_range())
1081  {
1082  // We need the algebraic data
1083  _femcontext.pre_fe_reinit(*this, elem);
1084  // And when asserts are on, we also need the FE so
1085  // we can assert that the mesh data is of the right type.
1086 #ifndef NDEBUG
1087  _femcontext.elem_fe_reinit();
1088 #endif
1089 
1090  // This code won't handle moving subactive elements
1091  libmesh_assert(!_femcontext.get_elem().has_children());
1092 
1093  _femcontext.elem_position_set(1.);
1094  }
1095 
1096  // We've now got positions set on all local nodes (and some
1097  // semilocal nodes); let's request positions for non-local nodes
1098  // from their processors.
1099 
1100  SyncNodalPositions sync_object(mesh);
1102  (this->comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
1103 }
1104 
1105 
1106 
1108 {
1109  LOG_SCOPE("postprocess()", "FEMSystem");
1110 
1111  const MeshBase & mesh = this->get_mesh();
1112 
1113  this->update();
1114 
1115  // Get the time solver object associated with the system, and tell it that
1116  // we are not solving the adjoint problem
1117  this->get_time_solver().set_is_adjoint(false);
1118 
1119  // Loop over every active mesh element on this processor
1121  mesh.active_local_elements_end()),
1122  PostprocessContributions(*this));
1123 }
1124 
1125 
1126 
1127 void FEMSystem::assemble_qoi (const QoISet & qoi_indices)
1128 {
1129  LOG_SCOPE("assemble_qoi()", "FEMSystem");
1130 
1131  const MeshBase & mesh = this->get_mesh();
1132 
1133  this->update();
1134 
1135  const unsigned int Nq = cast_int<unsigned int>(qoi.size());
1136 
1137  // the quantity of interest is assumed to be a sum of element and
1138  // side terms
1139  for (unsigned int i=0; i != Nq; ++i)
1140  if (qoi_indices.has_index(i))
1141  qoi[i] = 0;
1142 
1143  // Create a non-temporary qoi_contributions object, so we can query
1144  // its results after the reduction
1145  QoIContributions qoi_contributions(*this, *(this->diff_qoi), qoi_indices);
1146 
1147  // Loop over every active mesh element on this processor
1149  mesh.active_local_elements_end()),
1150  qoi_contributions);
1151 
1152  this->diff_qoi->parallel_op( this->comm(), this->qoi, qoi_contributions.qoi, qoi_indices );
1153 }
1154 
1155 
1156 
1157 void FEMSystem::assemble_qoi_derivative (const QoISet & qoi_indices,
1158  bool include_liftfunc,
1159  bool apply_constraints)
1160 {
1161  LOG_SCOPE("assemble_qoi_derivative()", "FEMSystem");
1162 
1163  const MeshBase & mesh = this->get_mesh();
1164 
1165  this->update();
1166 
1167  // The quantity of interest derivative assembly accumulates on
1168  // initially zero vectors
1169  for (std::size_t i=0; i != qoi.size(); ++i)
1170  if (qoi_indices.has_index(i))
1171  this->add_adjoint_rhs(i).zero();
1172 
1173  // Loop over every active mesh element on this processor
1175  mesh.active_local_elements_end()),
1176  QoIDerivativeContributions(*this, qoi_indices,
1177  *(this->diff_qoi),
1178  include_liftfunc,
1179  apply_constraints));
1180 
1181  for (std::size_t i=0; i != qoi.size(); ++i)
1182  if (qoi_indices.has_index(i))
1183  this->diff_qoi->finalize_derivative(this->get_adjoint_rhs(i),i);
1184 }
1185 
1186 
1187 
1189  FEMContext & context) const
1190 {
1191  // Logging is done by numerical_elem_jacobian
1192  // or numerical_side_jacobian
1193 
1194  DenseVector<Number> original_residual(context.get_elem_residual());
1195  DenseVector<Number> backwards_residual(context.get_elem_residual());
1196  DenseMatrix<Number> numeric_jacobian(context.get_elem_jacobian());
1197 #ifdef DEBUG
1198  DenseMatrix<Number> old_jacobian(context.get_elem_jacobian());
1199 #endif
1200 
1201  Real numerical_point_h = 0.;
1202  if (_mesh_sys == this)
1203  numerical_point_h = numerical_jacobian_h * context.get_elem().hmin();
1204 
1205  for (unsigned int v = 0; v != context.n_vars(); ++v)
1206  {
1207  const Real my_h = this->numerical_jacobian_h_for_var(v);
1208 
1209  unsigned int j_offset = libMesh::invalid_uint;
1210 
1211  if (!context.get_dof_indices(v).empty())
1212  {
1213  for (std::size_t i = 0; i != context.get_dof_indices().size(); ++i)
1214  if (context.get_dof_indices()[i] ==
1215  context.get_dof_indices(v)[0])
1216  j_offset = i;
1217 
1218  libmesh_assert_not_equal_to(j_offset, libMesh::invalid_uint);
1219  }
1220 
1221  for (std::size_t j = 0; j != context.get_dof_indices(v).size(); ++j)
1222  {
1223  const unsigned int total_j = j + j_offset;
1224 
1225  // Take the "minus" side of a central differenced first derivative
1226  Number original_solution = context.get_elem_solution(v)(j);
1227  context.get_elem_solution(v)(j) -= my_h;
1228 
1229  // Make sure to catch any moving mesh terms
1230  Real * coord = libmesh_nullptr;
1231  if (_mesh_sys == this)
1232  {
1233  if (_mesh_x_var == v)
1234  coord = &(context.get_elem().point(j)(0));
1235  else if (_mesh_y_var == v)
1236  coord = &(context.get_elem().point(j)(1));
1237  else if (_mesh_z_var == v)
1238  coord = &(context.get_elem().point(j)(2));
1239  }
1240  if (coord)
1241  {
1242  // We have enough information to scale the perturbations
1243  // here appropriately
1244  context.get_elem_solution(v)(j) = original_solution - numerical_point_h;
1245  *coord = libmesh_real(context.get_elem_solution(v)(j));
1246  }
1247 
1248  context.get_elem_residual().zero();
1249  ((*time_solver).*(res))(false, context);
1250 #ifdef DEBUG
1251  libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
1252 #endif
1253  backwards_residual = context.get_elem_residual();
1254 
1255  // Take the "plus" side of a central differenced first derivative
1256  context.get_elem_solution(v)(j) = original_solution + my_h;
1257  if (coord)
1258  {
1259  context.get_elem_solution()(j) = original_solution + numerical_point_h;
1260  *coord = libmesh_real(context.get_elem_solution(v)(j));
1261  }
1262  context.get_elem_residual().zero();
1263  ((*time_solver).*(res))(false, context);
1264 #ifdef DEBUG
1265  libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
1266 #endif
1267 
1268  context.get_elem_solution(v)(j) = original_solution;
1269  if (coord)
1270  {
1271  *coord = libmesh_real(context.get_elem_solution(v)(j));
1272  for (std::size_t i = 0; i != context.get_dof_indices().size(); ++i)
1273  {
1274  numeric_jacobian(i,total_j) =
1275  (context.get_elem_residual()(i) - backwards_residual(i)) /
1276  2. / numerical_point_h;
1277  }
1278  }
1279  else
1280  {
1281  for (std::size_t i = 0; i != context.get_dof_indices().size(); ++i)
1282  {
1283  numeric_jacobian(i,total_j) =
1284  (context.get_elem_residual()(i) - backwards_residual(i)) /
1285  2. / my_h;
1286  }
1287  }
1288  }
1289  }
1290 
1291  context.get_elem_residual() = original_residual;
1292  context.get_elem_jacobian() = numeric_jacobian;
1293 }
1294 
1295 
1296 
1298 {
1299  LOG_SCOPE("numerical_elem_jacobian()", "FEMSystem");
1301 }
1302 
1303 
1304 
1306 {
1307  LOG_SCOPE("numerical_side_jacobian()", "FEMSystem");
1309 }
1310 
1311 
1312 
1314 {
1315  LOG_SCOPE("numerical_nonlocal_jacobian()", "FEMSystem");
1317 }
1318 
1319 
1320 
1322 {
1323  FEMContext * fc = new FEMContext(*this);
1324 
1325  DifferentiablePhysics * phys = this->get_physics();
1326 
1327  libmesh_assert (phys);
1328 
1329  // If we are solving a moving mesh problem, tell that to the Context
1330  fc->set_mesh_system(phys->get_mesh_system());
1331  fc->set_mesh_x_var(phys->get_mesh_x_var());
1332  fc->set_mesh_y_var(phys->get_mesh_y_var());
1333  fc->set_mesh_z_var(phys->get_mesh_z_var());
1334 
1335  fc->set_deltat_pointer( &deltat );
1336 
1337  // If we are solving the adjoint problem, tell that to the Context
1338  fc->is_adjoint() = this->get_time_solver().is_adjoint();
1339 
1340  return UniquePtr<DiffContext>(fc);
1341 }
1342 
1343 
1344 
1346 {
1347  // Parent::init_context(c); // may be a good idea in derived classes
1348 
1349  // Although we do this in DiffSystem::build_context() and
1350  // FEMSystem::build_context() as well, we do it here just to be
1351  // extra sure that the deltat pointer gets set. Since the
1352  // intended behavior is for classes derived from FEMSystem to
1353  // call Parent::init_context() in their own init_context()
1354  // overloads, we can ensure that those classes get the correct
1355  // deltat pointers even if they have different build_context()
1356  // overloads.
1357  c.set_deltat_pointer ( &deltat );
1358 
1359  FEMContext & context = cast_ref<FEMContext &>(c);
1360 
1361  // Make sure we're prepared to do mass integration
1362  for (unsigned int var = 0; var != this->n_vars(); ++var)
1363  if (this->get_physics()->is_time_evolving(var))
1364  {
1365  // Request shape functions based on FEType
1366  switch( FEInterface::field_type( this->variable_type(var) ) )
1367  {
1368  case( TYPE_SCALAR ):
1369  {
1370  FEBase * elem_fe = libmesh_nullptr;
1371  context.get_element_fe(var, elem_fe);
1372  elem_fe->get_JxW();
1373  elem_fe->get_phi();
1374  }
1375  break;
1376  case( TYPE_VECTOR ):
1377  {
1379  context.get_element_fe(var, elem_fe);
1380  elem_fe->get_JxW();
1381  elem_fe->get_phi();
1382  }
1383  break;
1384  default:
1385  libmesh_error_msg("Unrecognized field type!");
1386  }
1387  }
1388 }
1389 
1390 
1391 
1393 {
1394  // This function makes no sense unless we've already picked out some
1395  // variable(s) to reflect mesh position coordinates
1396  if (!_mesh_sys)
1397  libmesh_error_msg("_mesh_sys was NULL!");
1398 
1399  // We currently assume mesh variables are in our own system
1400  if (_mesh_sys != this)
1401  libmesh_not_implemented();
1402 
1403  // Loop over every active mesh element on this processor
1404  const MeshBase & mesh = this->get_mesh();
1405 
1406  UniquePtr<DiffContext> con = this->build_context();
1407  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
1408  this->init_context(_femcontext);
1409 
1410  // Get the solution's mesh variables from every element
1411  for (const auto & elem : mesh.active_local_element_ptr_range())
1412  {
1413  _femcontext.pre_fe_reinit(*this, elem);
1414 
1415  _femcontext.elem_position_get();
1416 
1418  this->solution->insert(_femcontext.get_elem_solution(_mesh_x_var),
1419  _femcontext.get_dof_indices(_mesh_x_var) );
1421  this->solution->insert(_femcontext.get_elem_solution(_mesh_y_var),
1422  _femcontext.get_dof_indices(_mesh_y_var));
1424  this->solution->insert(_femcontext.get_elem_solution(_mesh_z_var),
1425  _femcontext.get_dof_indices(_mesh_z_var));
1426  }
1427 
1428  this->solution->close();
1429 
1430  // And make sure the current_local_solution is up to date too
1431  this->System::update();
1432 }
1433 
1434 } // namespace libMesh
T libmesh_real(T a)
OStreamProxy err
FEFamily family
The type of finite element.
Definition: fe_type.h:203
bool is_adjoint() const
Accessor for querying whether we need to do a primal or adjoint solve.
Definition: time_solver.h:232
virtual bool side_residual(bool request_jacobian, DiffContext &)=0
This method uses the DifferentiablePhysics side_time_derivative(), side_constraint(), and side_mass_residual() to build a full residual on an element&#39;s side.
bool has_children() const
Definition: elem.h:2295
const FEType & type() const
Definition: variable.h:119
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
This is the EquationSystems class.
virtual void parallel_op(const Parallel::Communicator &communicator, std::vector< Number > &sys_qoi, std::vector< Number > &local_qoi, const QoISet &qoi_indices)
Method to populate system qoi data structure with process-local qoi.
Definition: diff_qoi.C:39
void numerical_side_jacobian(FEMContext &context) const
Uses the results of multiple side_residual() calls to numerically differentiate the corresponding jac...
Definition: fem_system.C:1305
unsigned int _mesh_x_var
Variables from which to acquire moving mesh information.
Definition: diff_physics.h:546
void set_mesh_z_var(unsigned int z_var)
Accessor for z-variable of moving mesh System.
Definition: fem_context.h:859
Real verify_analytic_jacobians
If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calc...
Definition: fem_system.h:203
const DenseVector< Number > & get_elem_solution() const
Accessor for element solution.
Definition: diff_context.h:110
void elem_position_set(Real theta)
Uses the coordinate data specified by mesh_*_position configuration to set the geometry of elem to th...
Definition: fem_context.h:1195
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1582
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
Real numerical_jacobian_h_for_var(unsigned int var_num) const
If numerical_jacobian_h_for_var(var_num) is changed from its default value (numerical_jacobian_h), the FEMSystem will perturb solution vector entries for variable var_num by that amount when calculating finite differences with respect to that variable.
Definition: fem_system.h:253
void set_mesh_x_var(unsigned int x_var)
Accessor for x-variable of moving mesh System.
Definition: fem_context.h:831
bool is_adjoint() const
Accessor for querying whether we need to do a primal or adjoint solve.
Definition: diff_context.h:453
virtual void zero() libmesh_override
Set every element in the matrix to 0.
Definition: dense_matrix.h:792
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
Definition: threads_none.h:73
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Definition: diff_system.h:329
ImplicitSystem & sys
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: diff_system.C:108
bool has_elem() const
Test for current Elem object.
Definition: fem_context.h:865
Real numerical_jacobian_h
If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry b...
Definition: fem_system.h:175
Dummy "splitting object" used to distinguish splitting constructors from copy constructors.
Definition: threads_none.h:63
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
Constrains the element vector.
Definition: dof_map.h:1802
processor_id_type n_processors() const
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
virtual Real hmin() const
Definition: elem.C:458
MeshBase & mesh
Real l1_norm() const
Definition: dense_matrix.h:990
const DenseVector< Number > & get_elem_fixed_solution() const
Accessor for element fixed solution.
Definition: diff_context.h:214
static FEFieldType field_type(const FEType &fe_type)
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
void get_transpose(DenseMatrix< T > &dest) const
Put the tranposed matrix into dest.
Definition: dense_matrix.C:576
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
bool postprocess_sides
If postprocess_sides is true (it is false by default), the postprocessing loop will loop over all sid...
Definition: diff_system.h:302
static const Real TOLERANCE
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:52
bool(TimeSolver::* TimeSolverResPtr)(bool, DiffContext &)
Syntax sugar to make numerical_jacobian() declaration easier.
Definition: fem_system.h:208
bool is_time_evolving(unsigned int var) const
Definition: diff_physics.h:276
void elem_position_get()
Uses the geometry of elem to set the coordinate data specified by mesh_*_position configuration...
Definition: fem_context.C:1441
The libMesh namespace provides an interface to certain functionality in the library.
const_iterator begin() const
Beginning of the range.
Definition: stored_range.h:261
const VariableGroup & variable_group(unsigned int vg) const
Return a constant reference to VariableGroup vg.
Definition: system.h:2124
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:334
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
virtual Real l1_norm() const =0
long double max(long double a, double b)
This class provides a specific system class.
Definition: fem_system.h:53
virtual bool nonlocal_residual(bool request_jacobian, DiffContext &)=0
This method uses the DifferentiablePhysics nonlocal_time_derivative(), nonlocal_constraint(), and nonlocal_mass_residual() to build a full residual of non-local terms.
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
virtual UniquePtr< DiffContext > build_context() libmesh_override
Builds a FEMContext object with enough information to do evaluations on each element.
Definition: fem_system.C:1321
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: diff_system.C:152
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
bool print_element_residuals
Set print_element_residuals to true to print each R_elem contribution.
Definition: diff_system.h:344
libmesh_assert(j)
DifferentiableQoI * diff_qoi
Pointer to object to use for quantity of interest assembly evaluations.
Definition: diff_system.h:365
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
const System * get_mesh_system() const
Definition: diff_physics.h:623
void numerical_nonlocal_jacobian(FEMContext &context) const
Uses the results of multiple side_residual() calls to numerically differentiate the corresponding jac...
Definition: fem_system.C:1313
unsigned int get_mesh_x_var() const
Definition: diff_physics.h:635
virtual void set_mesh_system(System *sys)
Tells the FEMContext that system sys contains the isoparametric Lagrangian variables which correspond...
Definition: fem_context.h:805
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
bool print_element_solutions
Set print_element_solutions to true to print each U_elem input.
Definition: diff_system.h:339
virtual node_iterator nodes_begin()=0
Iterate over all the nodes in the Mesh.
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
virtual bool element_residual(bool request_jacobian, DiffContext &)=0
This method uses the DifferentiablePhysics element_time_derivative(), element_constraint(), and mass_residual() to build a full residual on an element.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
bool print_solution_norms
Set print_residual_norms to true to print |U| whenever it is used in an assembly() call...
Definition: diff_system.h:308
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
void set_is_adjoint(bool _is_adjoint_value)
Accessor for setting whether we need to do a primal or adjoint solve.
Definition: time_solver.h:239
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
virtual void assemble_qoi_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true) libmesh_override
Runs a qoi derivative assembly loop over all elements, and if assemble_qoi_sides is true over all sid...
Definition: fem_system.C:1157
unsigned int n_variable_groups() const
Definition: system.h:2094
const MeshBase & get_mesh() const
Definition: system.h:2014
const std::vector< DenseVector< Number > > & get_qoi_derivatives() const
Const accessor for QoI derivatives.
Definition: diff_context.h:330
const DofMap & get_dof_map() const
Definition: system.h:2030
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
virtual void zero()=0
Set all entries to 0.
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
Definition: diff_system.h:319
virtual void side_fe_reinit()
Reinitializes side FE objects on the current geometric element.
Definition: fem_context.C:1405
virtual void side_postprocess(DiffContext &)
Does any work that needs to be done on side of elem in a postprocessing loop.
Definition: diff_system.h:275
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
unsigned int m() const
virtual element_iterator active_local_elements_begin()=0
bool fe_reinit_during_postprocess
If fe_reinit_during_postprocess is true (it is true by default), FE objects will be reinit()ed with t...
Definition: fem_system.h:164
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
bool use_fixed_solution
A boolean to be set to true by systems using elem_fixed_solution, for optional use by e...
Definition: system.h:1493
const std::vector< Number > & get_qois() const
Const accessor for QoI vector.
Definition: diff_context.h:318
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Adds factor times mat to this matrix.
Definition: dense_matrix.h:883
System * _mesh_sys
System from which to acquire moving mesh information.
Definition: diff_physics.h:541
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
Definition: diff_system.h:324
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
void constrain_nothing(std::vector< dof_id_type > &dofs) const
Does not actually constrain anything, but modifies dofs in the same way as any of the constrain funct...
FEMSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
Definition: fem_system.C:828
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
Number has_heterogenous_adjoint_constraint(const unsigned int qoi_num, const dof_id_type dof) const
Definition: dof_map.h:1758
virtual unsigned int n_sides() const =0
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
const_iterator end() const
End of the range.
Definition: stored_range.h:266
virtual void elem_fe_reinit(const std::vector< Point > *const pts=libmesh_nullptr)
Reinitializes interior FE objects on the current geometric element.
Definition: fem_context.C:1382
const DifferentiablePhysics * get_physics() const
Definition: diff_system.h:169
virtual void element_postprocess(DiffContext &)
Does any work that needs to be done on elem in a postprocessing loop.
Definition: diff_system.h:269
NumericVector< Number > & add_adjoint_rhs(unsigned int i=0)
Definition: system.C:1041
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
bool print_solutions
Set print_solutions to true to print U whenever it is used in an assembly() call. ...
Definition: diff_system.h:314
std::vector< object_type >::const_iterator const_iterator
Allows an StoredRange to behave like an STL container.
Definition: stored_range.h:58
void heterogenously_constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
Constrains the element matrix and vector.
This class provides a specific system class.
Definition: diff_physics.h:74
This class provides a specific system class.
Definition: diff_qoi.h:49
bool compute_internal_sides
compute_internal_sides is false by default, indicating that side_* computations will only be done on ...
Definition: diff_physics.h:152
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
const Point & point(const unsigned int i) const
Definition: elem.h:1809
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix.
Definition: dof_map.h:1793
StoredRange< iterator_type, object_type > & reset(const iterator_type &first, const iterator_type &last)
Resets the StoredRange to contain [first,last).
Definition: stored_range.h:222
unsigned char side
Current side for side_* to examine.
Definition: fem_context.h:977
const Parallel::Communicator & comm() const
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int get_mesh_z_var() const
Definition: diff_physics.h:647
unsigned int get_mesh_y_var() const
Definition: diff_physics.h:641
void mesh_position_set()
Tells the FEMSystem to set the mesh nodal coordinates which should correspond to degree of freedom co...
Definition: fem_system.C:1067
void numerical_elem_jacobian(FEMContext &context) const
Uses the results of multiple element_residual() calls to numerically differentiate the corresponding ...
Definition: fem_system.C:1297
void set_mesh_y_var(unsigned int y_var)
Accessor for y-variable of moving mesh System.
Definition: fem_context.h:845
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override
Prepares matrix or rhs for matrix assembly.
Definition: fem_system.C:852
void heterogenously_constrain_element_vector(const DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
Constrains the element vector.
bool print_element_jacobians
Set print_element_jacobians to true to print each J_elem contribution.
Definition: diff_system.h:349
void parallel_reduce(const Range &range, Body &body)
Execute the provided reduction operation in parallel on the specified range.
Definition: threads_none.h:101
void set_deltat_pointer(Real *dt)
Points the _deltat member of this class at a timestep value stored in the creating System...
Definition: diff_context.C:138
bool has_heterogenous_adjoint_constraints(const unsigned int qoi_num) const
Definition: dof_map.h:1744
void numerical_jacobian(TimeSolverResPtr res, FEMContext &context) const
Uses the results of multiple res calls to numerically differentiate the corresponding jacobian...
Definition: fem_system.C:1188
virtual ~FEMSystem()
Destructor.
Definition: fem_system.C:839
unsigned int n_vars() const
Definition: system.h:2086
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
dof_id_type id() const
Definition: dof_object.h:632
void mesh_position_get()
Tells the FEMSystem to set the degree of freedom coefficients which should correspond to mesh nodal c...
Definition: fem_system.C:1392
virtual element_iterator active_local_elements_end()=0
virtual void assemble_qoi(const QoISet &indices=QoISet()) libmesh_override
Runs a qoi assembly loop over all elements, and if assemble_qoi_sides is true over all sides...
Definition: fem_system.C:1127
virtual Real l1_norm() const =0
std::streamsize precision() const
Get the associated write precision.
bool has_index(unsigned int) const
Return whether or not this index is in the set to be calculated.
Definition: qoi_set.h:221
This class forms the foundation from which generic finite elements may be derived.
virtual void finalize_derivative(NumericVector< Number > &derivatives, std::size_t qoi_index)
Method to finalize qoi derivatives which require more than just a simple sum of element contributions...
Definition: diff_qoi.C:51
virtual void postprocess() libmesh_override
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: fem_system.C:1107
processor_id_type processor_id() const
TimeSolver & get_time_solver()
Definition: diff_system.h:402
unsigned int n_vars() const
Number of variables in solution.
Definition: diff_context.h:98
virtual void init_context(DiffContext &) libmesh_override
Definition: fem_system.C:1345
NumericVector< Number > & get_adjoint_rhs(unsigned int i=0)
Definition: system.C:1051
This class provides a specific system class.