fem_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
23 #include "libmesh/fe_base.h"
24 #include "libmesh/fem_context.h"
25 #include "libmesh/fem_system.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/numeric_vector.h"
29 #include "libmesh/parallel.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
50  (const FEMSystem& _sys,
51  const bool _get_jacobian,
52  FEMContext &_femcontext)
53  {
54  bool jacobian_computed =
55  _sys.time_solver->element_residual(_get_jacobian, _femcontext);
56 
57  // Compute a numeric jacobian if we have to
58  if (_get_jacobian && !jacobian_computed)
59  {
60  // Make sure we didn't compute a jacobian and lie about it
61  libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
62  // Logging of numerical jacobians is done separately
63  _sys.numerical_elem_jacobian(_femcontext);
64  }
65 
66  // Compute a numeric jacobian if we're asked to verify the
67  // analytic jacobian we got
68  if (_get_jacobian && jacobian_computed &&
69  _sys.verify_analytic_jacobians != 0.0)
70  {
71  DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
72 
73  _femcontext.get_elem_jacobian().zero();
74  // Logging of numerical jacobians is done separately
75  _sys.numerical_elem_jacobian(_femcontext);
76 
77  Real analytic_norm = analytic_jacobian.l1_norm();
78  Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
79 
80  // If we can continue, we'll probably prefer the analytic jacobian
81  analytic_jacobian.swap(_femcontext.get_elem_jacobian());
82 
83  // The matrix "analytic_jacobian" will now hold the error matrix
84  analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
85  Real error_norm = analytic_jacobian.l1_norm();
86 
87  Real relative_error = error_norm /
88  std::max(analytic_norm, numerical_norm);
89 
90  if (relative_error > _sys.verify_analytic_jacobians)
91  {
92  libMesh::err << "Relative error " << relative_error
93  << " detected in analytic jacobian on element "
94  << _femcontext.get_elem().id() << '!' << std::endl;
95 
96  std::streamsize old_precision = libMesh::out.precision();
98  libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
99  << _femcontext.get_elem_jacobian() << std::endl;
100  analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
101  libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
102  << analytic_jacobian << std::endl;
103 
104  libMesh::out.precision(old_precision);
105 
106  libmesh_error();
107  }
108  }
109 
110  for (_femcontext.side = 0;
111  _femcontext.side != _femcontext.get_elem().n_sides();
112  ++_femcontext.side)
113  {
114  // Don't compute on non-boundary sides unless requested
115  if (!_sys.get_physics()->compute_internal_sides &&
116  _femcontext.get_elem().neighbor(_femcontext.side) != NULL)
117  continue;
118 
119  // Any mesh movement has already been done (and restored,
120  // if the TimeSolver isn't broken), but
121  // reinitializing the side FE objects is still necessary
122  _femcontext.side_fe_reinit();
123 
124  DenseMatrix<Number> old_jacobian;
125  // If we're in DEBUG mode, we should always verify that the
126  // user's side_residual function doesn't alter our existing
127  // jacobian and then lie about it
128 #ifndef DEBUG
129  // Even if we're not in DEBUG mode, when we're verifying
130  // analytic jacobians we'll want to verify each side's
131  // jacobian contribution separately.
132  /* PB: We also need to account for the case when the user wants to
133  use numerical Jacobians and not analytic Jacobians */
134  if ( (_sys.verify_analytic_jacobians != 0.0 && _get_jacobian) ||
135  (!jacobian_computed && _get_jacobian) )
136 #endif // ifndef DEBUG
137  {
138  old_jacobian = _femcontext.get_elem_jacobian();
139  _femcontext.get_elem_jacobian().zero();
140  }
141  jacobian_computed =
142  _sys.time_solver->side_residual(_get_jacobian, _femcontext);
143 
144  // Compute a numeric jacobian if we have to
145  if (_get_jacobian && !jacobian_computed)
146  {
147  // In DEBUG mode, we've already set elem_jacobian == 0,
148  // so we can make sure side_residual didn't compute a
149  // jacobian and lie about it
150 #ifdef DEBUG
151  libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
152 #endif
153  // Logging of numerical jacobians is done separately
154  _sys.numerical_side_jacobian(_femcontext);
155 
156  // Add back in element interior numerical Jacobian
157  _femcontext.get_elem_jacobian() += old_jacobian;
158  }
159 
160  // Compute a numeric jacobian if we're asked to verify the
161  // analytic jacobian we got
162  if (_get_jacobian && jacobian_computed &&
163  _sys.verify_analytic_jacobians != 0.0)
164  {
165  DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
166 
167  _femcontext.get_elem_jacobian().zero();
168  // Logging of numerical jacobians is done separately
169  _sys.numerical_side_jacobian(_femcontext);
170 
171  Real analytic_norm = analytic_jacobian.l1_norm();
172  Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
173 
174  // If we can continue, we'll probably prefer the analytic jacobian
175  analytic_jacobian.swap(_femcontext.get_elem_jacobian());
176 
177  // The matrix "analytic_jacobian" will now hold the error matrix
178  analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
179  Real error_norm = analytic_jacobian.l1_norm();
180 
181  Real relative_error = error_norm /
182  std::max(analytic_norm, numerical_norm);
183 
184  if (relative_error > _sys.verify_analytic_jacobians)
185  {
186  libMesh::err << "Relative error " << relative_error
187  << " detected in analytic jacobian on element "
188  << _femcontext.get_elem().id()
189  << ", side "
190  << static_cast<unsigned int>(_femcontext.side) << '!' << std::endl;
191 
192  std::streamsize old_precision = libMesh::out.precision();
194  libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
195  << _femcontext.get_elem_jacobian() << std::endl;
196  analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
197  libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
198  << analytic_jacobian << std::endl;
199  libMesh::out.precision(old_precision);
200 
201  libmesh_error();
202  }
203  // Once we've verified a side, we'll want to add back the
204  // rest of the accumulated jacobian
205  _femcontext.get_elem_jacobian() += old_jacobian;
206  }
207  // In DEBUG mode, we've set elem_jacobian == 0, and we
208  // may still need to add the old jacobian back
209 #ifdef DEBUG
210  if (_get_jacobian && jacobian_computed &&
211  _sys.verify_analytic_jacobians == 0.0)
212  {
213  _femcontext.get_elem_jacobian() += old_jacobian;
214  }
215 #endif // ifdef DEBUG
216  }
217  }
218 
219  class AssemblyContributions
220  {
221  public:
225  AssemblyContributions(FEMSystem &sys,
226  bool get_residual,
227  bool get_jacobian) :
228  _sys(sys),
229  _get_residual(get_residual),
230  _get_jacobian(get_jacobian) {}
231 
235  void operator()(const ConstElemRange &range) const
236  {
237  AutoPtr<DiffContext> con = _sys.build_context();
238  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
239  _sys.init_context(_femcontext);
240 
241  for (ConstElemRange::const_iterator elem_it = range.begin();
242  elem_it != range.end(); ++elem_it)
243  {
244  Elem *el = const_cast<Elem *>(*elem_it);
245 
246  _femcontext.pre_fe_reinit(_sys, el);
247  _femcontext.elem_fe_reinit();
248 
249  assemble_unconstrained_element_system
250  (_sys, _get_jacobian, _femcontext);
251 
252 #ifdef LIBMESH_ENABLE_CONSTRAINTS
253  // We turn off the asymmetric constraint application;
254  // enforce_constraints_exactly() should be called in the solver
255  if (_get_residual && _get_jacobian)
257  (_femcontext.get_elem_jacobian(), _femcontext.get_elem_residual(),
258  _femcontext.get_dof_indices(), false);
259  else if (_get_residual)
261  (_femcontext.get_elem_residual(), _femcontext.get_dof_indices(), false);
262  else if (_get_jacobian)
264  (_femcontext.get_elem_jacobian(), _femcontext.get_dof_indices(), false);
265 #endif // #ifdef LIBMESH_ENABLE_CONSTRAINTS
266 
267  if (_get_jacobian && _sys.print_element_jacobians)
268  {
269  std::streamsize old_precision = libMesh::out.precision();
271  libMesh::out << "J_elem " << _femcontext.get_elem().id() << " = "
272  << _femcontext.get_elem_jacobian() << std::endl;
273  libMesh::out.precision(old_precision);
274  }
275 
276  { // A lock is necessary around access to the global system
277  femsystem_mutex::scoped_lock lock(assembly_mutex);
278 
279  if (_get_jacobian)
280  _sys.matrix->add_matrix (_femcontext.get_elem_jacobian(),
281  _femcontext.get_dof_indices());
282  if (_get_residual)
283  _sys.rhs->add_vector (_femcontext.get_elem_residual(),
284  _femcontext.get_dof_indices());
285  } // Scope for assembly mutex
286 
287  }
288  }
289 
290  private:
291 
292  FEMSystem& _sys;
293 
294  const bool _get_residual, _get_jacobian;
295  };
296 
297  class PostprocessContributions
298  {
299  public:
303  explicit
304  PostprocessContributions(FEMSystem &sys) : _sys(sys) {}
305 
309  void operator()(const ConstElemRange &range) const
310  {
311  AutoPtr<DiffContext> con = _sys.build_context();
312  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
313  _sys.init_context(_femcontext);
314 
315  for (ConstElemRange::const_iterator elem_it = range.begin();
316  elem_it != range.end(); ++elem_it)
317  {
318  Elem *el = const_cast<Elem *>(*elem_it);
319  _femcontext.pre_fe_reinit(_sys, el);
320 
321  // Optionally initialize all the interior FE objects on elem.
323  _femcontext.elem_fe_reinit();
324 
325  _sys.element_postprocess(_femcontext);
326 
327  for (_femcontext.side = 0;
328  _femcontext.side != _femcontext.get_elem().n_sides();
329  ++_femcontext.side)
330  {
331  // Don't compute on non-boundary sides unless requested
332  if (!_sys.postprocess_sides ||
334  _femcontext.get_elem().neighbor(_femcontext.side) != NULL))
335  continue;
336 
337  // Optionally initialize all the FE objects on this side.
339  _femcontext.side_fe_reinit();
340 
341  _sys.side_postprocess(_femcontext);
342  }
343  }
344  }
345 
346  private:
347 
348  FEMSystem& _sys;
349  };
350 
351  class QoIContributions
352  {
353  public:
357  explicit
358  QoIContributions(FEMSystem &sys, DifferentiableQoI &diff_qoi,
359  const QoISet &qoi_indices) :
360  qoi(sys.qoi.size(), 0.), _sys(sys), _diff_qoi(diff_qoi),_qoi_indices(qoi_indices) {}
361 
365  QoIContributions(const QoIContributions &other,
366  Threads::split) :
367  qoi(other._sys.qoi.size(), 0.), _sys(other._sys), _diff_qoi(other._diff_qoi) {}
368 
372  void operator()(const ConstElemRange &range)
373  {
374  AutoPtr<DiffContext> con = _sys.build_context();
375  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
376  _diff_qoi.init_context(_femcontext);
377 
378  for (ConstElemRange::const_iterator elem_it = range.begin();
379  elem_it != range.end(); ++elem_it)
380  {
381  Elem *el = const_cast<Elem *>(*elem_it);
382 
383  _femcontext.pre_fe_reinit(_sys, el);
384 
385  if (_diff_qoi.assemble_qoi_elements)
386  {
387  _femcontext.elem_fe_reinit();
388 
389  _diff_qoi.element_qoi(_femcontext, _qoi_indices);
390  }
391 
392  for (_femcontext.side = 0;
393  _femcontext.side != _femcontext.get_elem().n_sides();
394  ++_femcontext.side)
395  {
396  // Don't compute on non-boundary sides unless requested
397  if (!_diff_qoi.assemble_qoi_sides ||
398  (!_diff_qoi.assemble_qoi_internal_sides &&
399  _femcontext.get_elem().neighbor(_femcontext.side) != NULL))
400  continue;
401 
402  _femcontext.side_fe_reinit();
403 
404  _diff_qoi.side_qoi(_femcontext, _qoi_indices);
405  }
406  }
407 
408  this->_diff_qoi.thread_join( this->qoi, _femcontext.get_qois(), _qoi_indices );
409  }
410 
411  void join (const QoIContributions& other)
412  {
413  libmesh_assert_equal_to (this->qoi.size(), other.qoi.size());
414  this->_diff_qoi.thread_join( this->qoi, other.qoi, _qoi_indices );
415  }
416 
417  std::vector<Number> qoi;
418 
419  private:
420 
421  FEMSystem& _sys;
422  DifferentiableQoI& _diff_qoi;
423 
424  const QoISet _qoi_indices;
425  };
426 
427  class QoIDerivativeContributions
428  {
429  public:
433  QoIDerivativeContributions(FEMSystem &sys, const QoISet& qoi_indices,
434  DifferentiableQoI &qoi ) :
435  _sys(sys), _qoi_indices(qoi_indices), _qoi(qoi) {}
436 
440  void operator()(const ConstElemRange &range) const
441  {
442  AutoPtr<DiffContext> con = _sys.build_context();
443  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
444  _qoi.init_context(_femcontext);
445 
446  for (ConstElemRange::const_iterator elem_it = range.begin();
447  elem_it != range.end(); ++elem_it)
448  {
449  Elem *el = const_cast<Elem *>(*elem_it);
450 
451  _femcontext.pre_fe_reinit(_sys, el);
452 
453  if (_qoi.assemble_qoi_elements)
454  {
455  _femcontext.elem_fe_reinit();
456 
457  _qoi.element_qoi_derivative(_femcontext, _qoi_indices);
458  }
459 
460  for (_femcontext.side = 0;
461  _femcontext.side != _femcontext.get_elem().n_sides();
462  ++_femcontext.side)
463  {
464  // Don't compute on non-boundary sides unless requested
465  if (!_qoi.assemble_qoi_sides ||
466  (!_qoi.assemble_qoi_internal_sides &&
467  _femcontext.get_elem().neighbor(_femcontext.side) != NULL))
468  continue;
469 
470  _femcontext.side_fe_reinit();
471 
472  _qoi.side_qoi_derivative(_femcontext, _qoi_indices);
473  }
474 
475  // We need some unmodified indices to use for constraining
476  // multiple vector
477  // FIXME - there should be a DofMap::constrain_element_vectors
478  // to do this more efficiently
479  std::vector<dof_id_type> original_dofs = _femcontext.get_dof_indices();
480 
481  { // A lock is necessary around access to the global system
482  femsystem_mutex::scoped_lock lock(assembly_mutex);
483 
484 #ifdef LIBMESH_ENABLE_CONSTRAINTS
485  // We'll need to see if any heterogenous constraints apply
486  // to the QoI dofs on this element *or* to any of the dofs
487  // they depend on, so let's get those dependencies
488  _sys.get_dof_map().constrain_nothing(_femcontext.get_dof_indices());
489 #endif
490 
491  for (unsigned int i=0; i != _sys.qoi.size(); ++i)
492  if (_qoi_indices.has_index(i))
493  {
494 #ifdef LIBMESH_ENABLE_CONSTRAINTS
495  bool has_heterogenous_constraint = false;
496  for (unsigned int d=0;
497  d != _femcontext.get_dof_indices().size(); ++d)
499  (i, _femcontext.get_dof_indices()[d]))
500  {
501  has_heterogenous_constraint = true;
502  break;
503  }
504 
505  _femcontext.get_dof_indices() = original_dofs;
506 
507  // If we're going to need K to impose a heterogenous
508  // constraint then we either already have it or we
509  // need to compute it
510  if (has_heterogenous_constraint)
511  {
512  assemble_unconstrained_element_system
513  (_sys, true, _femcontext);
514 
516  (_femcontext.get_elem_jacobian(),
517  _femcontext.get_qoi_derivatives()[i],
518  _femcontext.get_dof_indices(), false, i);
519  }
520  else
521  {
523  (_femcontext.get_qoi_derivatives()[i],
524  _femcontext.get_dof_indices(), false);
525  }
526 #endif
527 
529  (_femcontext.get_qoi_derivatives()[i], _femcontext.get_dof_indices());
530  }
531  }
532  }
533  }
534 
535  private:
536 
537  FEMSystem& _sys;
538  const QoISet& _qoi_indices;
539  DifferentiableQoI& _qoi;
540  };
541 
542 
543 }
544 
545 
546 namespace libMesh
547 {
548 
549 
550 
551 
552 
554  const std::string& name_in,
555  const unsigned int number_in)
556  : Parent(es, name_in, number_in),
557  fe_reinit_during_postprocess(true),
558  numerical_jacobian_h(TOLERANCE),
559  verify_analytic_jacobians(0.0)
560 {
561 }
562 
563 
565 {
566  this->clear();
567 }
568 
569 
570 
572 {
573  Parent::clear();
574 }
575 
576 
577 
579 {
580  // First initialize LinearImplicitSystem data
582 }
583 
584 
585 void FEMSystem::assembly (bool get_residual, bool get_jacobian)
586 {
587  libmesh_assert(get_residual || get_jacobian);
588  std::string log_name;
589  if (get_residual && get_jacobian)
590  log_name = "assembly()";
591  else if (get_residual)
592  log_name = "assembly(get_residual)";
593  else
594  log_name = "assembly(get_jacobian)";
595 
596  START_LOG(log_name, "FEMSystem");
597 
598  const MeshBase& mesh = this->get_mesh();
599 
600 // this->get_vector("_nonlinear_solution").localize
601 // (*current_local_nonlinear_solution,
602 // dof_map.get_send_list());
603  this->update();
604 
606  {
607 // this->get_vector("_nonlinear_solution").close();
608  this->solution->close();
609 
610  std::streamsize old_precision = libMesh::out.precision();
612  libMesh::out << "|U| = "
613 // << this->get_vector("_nonlinear_solution").l1_norm()
614  << this->solution->l1_norm()
615  << std::endl;
616  libMesh::out.precision(old_precision);
617  }
618  if (print_solutions)
619  {
620  std::streamsize old_precision = libMesh::out.precision();
622 // libMesh::out << "U = [" << this->get_vector("_nonlinear_solution")
623  libMesh::out << "U = [" << *(this->solution)
624  << "];" << std::endl;
625  libMesh::out.precision(old_precision);
626  }
627 
628  // Is this definitely necessary? [RHS]
629  // Yes. [RHS 2012]
630  if (get_jacobian)
631  matrix->zero();
632  if (get_residual)
633  rhs->zero();
634 
635  // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
636  if (verify_analytic_jacobians > 0.5)
637  {
638  libMesh::err << "WARNING! verify_analytic_jacobians was set "
639  << "to absurdly large value of "
640  << verify_analytic_jacobians << std::endl;
641  libMesh::err << "Resetting to 1e-6!" << std::endl;
643  }
644 
645  // In time-dependent problems, the nonlinear function we're trying
646  // to solve at each timestep may depend on the particular solver
647  // we're using
649 
650  // Build the residual and jacobian contributions on every active
651  // mesh element on this processor
652  Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(),
654  AssemblyContributions(*this, get_residual, get_jacobian));
655 
656 
657  if (get_residual && (print_residual_norms || print_residuals))
658  this->rhs->close();
659  if (get_residual && print_residual_norms)
660  {
661  std::streamsize old_precision = libMesh::out.precision();
663  libMesh::out << "|F| = " << this->rhs->l1_norm() << std::endl;
664  libMesh::out.precision(old_precision);
665  }
666  if (get_residual && print_residuals)
667  {
668  std::streamsize old_precision = libMesh::out.precision();
670  libMesh::out << "F = [" << *(this->rhs) << "];" << std::endl;
671  libMesh::out.precision(old_precision);
672  }
673 
674  if (get_jacobian && (print_jacobian_norms || print_jacobians))
675  this->matrix->close();
676  if (get_jacobian && print_jacobian_norms)
677  {
678  std::streamsize old_precision = libMesh::out.precision();
680  libMesh::out << "|J| = " << this->matrix->l1_norm() << std::endl;
681  libMesh::out.precision(old_precision);
682  }
683  if (get_jacobian && print_jacobians)
684  {
685  std::streamsize old_precision = libMesh::out.precision();
687  libMesh::out << "J = [" << *(this->matrix) << "];" << std::endl;
688  libMesh::out.precision(old_precision);
689  }
690  STOP_LOG(log_name, "FEMSystem");
691 }
692 
693 
694 
696 {
697  // We are solving the primal problem
698  Parent::solve();
699 
700  // On a moving mesh we want the mesh to reflect the new solution
701  this->mesh_position_set();
702 }
703 
704 
705 
707 {
708  // If we don't need to move the mesh, we're done
709  if (_mesh_sys != this)
710  return;
711 
712  MeshBase& mesh = this->get_mesh();
713 
714  AutoPtr<DiffContext> con = this->build_context();
715  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
716  this->init_context(_femcontext);
717 
718  // Move every mesh element we can
721  const MeshBase::const_element_iterator end_el =
723 
724  for ( ; el != end_el; ++el)
725  {
726  // We need the algebraic data
727  _femcontext.pre_fe_reinit(*this, *el);
728  // And when asserts are on, we also need the FE so
729  // we can assert that the mesh data is of the right type.
730 #ifndef NDEBUG
731  _femcontext.elem_fe_reinit();
732 #endif
733 
734  // This code won't handle moving subactive elements
735  libmesh_assert(!_femcontext.get_elem().has_children());
736 
737  _femcontext.elem_position_set(1.);
738  }
739 
740  // We've now got positions set on all local nodes (and some
741  // semilocal nodes); let's request positions for non-local nodes
742  // from their processors.
743 
744  SyncNodalPositions sync_object(mesh);
746  (this->comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
747 }
748 
749 
750 
752 {
753  START_LOG("postprocess()", "FEMSystem");
754 
755  const MeshBase& mesh = this->get_mesh();
756 
757  this->update();
758 
759  // Get the time solver object associated with the system, and tell it that
760  // we are not solving the adjoint problem
761  this->get_time_solver().set_is_adjoint(false);
762 
763  // Loop over every active mesh element on this processor
764  Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(),
766  PostprocessContributions(*this));
767 
768  STOP_LOG("postprocess()", "FEMSystem");
769 }
770 
771 
772 
773 void FEMSystem::assemble_qoi (const QoISet &qoi_indices)
774 {
775  START_LOG("assemble_qoi()", "FEMSystem");
776 
777  const MeshBase& mesh = this->get_mesh();
778 
779  this->update();
780 
781  const unsigned int Nq = libmesh_cast_int<unsigned int>(qoi.size());
782 
783  // the quantity of interest is assumed to be a sum of element and
784  // side terms
785  for (unsigned int i=0; i != Nq; ++i)
786  if (qoi_indices.has_index(i))
787  qoi[i] = 0;
788 
789  // Create a non-temporary qoi_contributions object, so we can query
790  // its results after the reduction
791  QoIContributions qoi_contributions(*this, *(this->diff_qoi), qoi_indices);
792 
793  // Loop over every active mesh element on this processor
796  qoi_contributions);
797 
798  this->diff_qoi->parallel_op( this->comm(), this->qoi, qoi_contributions.qoi, qoi_indices );
799 
800  STOP_LOG("assemble_qoi()", "FEMSystem");
801 }
802 
803 
804 
806 {
807  START_LOG("assemble_qoi_derivative()", "FEMSystem");
808 
809  const MeshBase& mesh = this->get_mesh();
810 
811  this->update();
812 
813  // The quantity of interest derivative assembly accumulates on
814  // initially zero vectors
815  for (unsigned int i=0; i != qoi.size(); ++i)
816  if (qoi_indices.has_index(i))
817  this->add_adjoint_rhs(i).zero();
818 
819  // Loop over every active mesh element on this processor
820  Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(),
822  QoIDerivativeContributions(*this, qoi_indices,
823  *(this->diff_qoi)));
824 
825  STOP_LOG("assemble_qoi_derivative()", "FEMSystem");
826 }
827 
828 
829 
830 void FEMSystem::numerical_jacobian (TimeSolverResPtr res,
831  FEMContext &context) const
832 {
833  // Logging is done by numerical_elem_jacobian
834  // or numerical_side_jacobian
835 
836  DenseVector<Number> original_residual(context.get_elem_residual());
837  DenseVector<Number> backwards_residual(context.get_elem_residual());
838  DenseMatrix<Number> numeric_jacobian(context.get_elem_jacobian());
839 #ifdef DEBUG
840  DenseMatrix<Number> old_jacobian(context.get_elem_jacobian());
841 #endif
842 
843  Real numerical_point_h = 0.;
844  if (_mesh_sys == this)
845  numerical_point_h = numerical_jacobian_h * context.get_elem().hmin();
846 
847  for (unsigned int j = 0; j != context.get_dof_indices().size(); ++j)
848  {
849  // Take the "minus" side of a central differenced first derivative
850  Number original_solution = context.get_elem_solution()(j);
851  context.get_elem_solution()(j) -= numerical_jacobian_h;
852 
853  // Make sure to catch any moving mesh terms
854  // FIXME - this could be less ugly
855  Real *coord = NULL;
856  if (_mesh_sys == this)
857  {
859  for (unsigned int k = 0;
860  k != context.get_dof_indices( _mesh_x_var ).size(); ++k)
861  if (context.get_dof_indices( _mesh_x_var )[k] ==
862  context.get_dof_indices()[j])
863  coord = &(context.get_elem().point(k)(0));
865  for (unsigned int k = 0;
866  k != context.get_dof_indices( _mesh_y_var ).size(); ++k)
867  if (context.get_dof_indices( _mesh_y_var )[k] ==
868  context.get_dof_indices()[j])
869  coord = &(context.get_elem().point(k)(1));
871  for (unsigned int k = 0;
872  k != context.get_dof_indices( _mesh_z_var ).size(); ++k)
873  if (context.get_dof_indices( _mesh_z_var )[k] ==
874  context.get_dof_indices()[j])
875  coord = &(context.get_elem().point(k)(2));
876  }
877  if (coord)
878  {
879  // We have enough information to scale the perturbations
880  // here appropriately
881  context.get_elem_solution()(j) = original_solution - numerical_point_h;
882  *coord = libmesh_real(context.get_elem_solution()(j));
883  }
884 
885  context.get_elem_residual().zero();
886  ((*time_solver).*(res))(false, context);
887 #ifdef DEBUG
888  libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
889 #endif
890  backwards_residual = context.get_elem_residual();
891 
892  // Take the "plus" side of a central differenced first derivative
893  context.get_elem_solution()(j) = original_solution + numerical_jacobian_h;
894  if (coord)
895  {
896  context.get_elem_solution()(j) = original_solution + numerical_point_h;
897  *coord = libmesh_real(context.get_elem_solution()(j));
898  }
899  context.get_elem_residual().zero();
900  ((*time_solver).*(res))(false, context);
901 #ifdef DEBUG
902  libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
903 #endif
904 
905  context.get_elem_solution()(j) = original_solution;
906  if (coord)
907  {
908  *coord = libmesh_real(context.get_elem_solution()(j));
909  for (unsigned int i = 0; i != context.get_dof_indices().size(); ++i)
910  {
911  numeric_jacobian(i,j) =
912  (context.get_elem_residual()(i) - backwards_residual(i)) /
913  2. / numerical_point_h;
914  }
915  }
916  else
917  {
918  for (unsigned int i = 0; i != context.get_dof_indices().size(); ++i)
919  {
920  numeric_jacobian(i,j) =
921  (context.get_elem_residual()(i) - backwards_residual(i)) /
923  }
924  }
925  }
926 
927  context.get_elem_residual() = original_residual;
928  context.get_elem_jacobian() = numeric_jacobian;
929 }
930 
931 
932 
934 {
935  START_LOG("numerical_elem_jacobian()", "FEMSystem");
937  STOP_LOG("numerical_elem_jacobian()", "FEMSystem");
938 }
939 
940 
941 
943 {
944  START_LOG("numerical_side_jacobian()", "FEMSystem");
946  STOP_LOG("numerical_side_jacobian()", "FEMSystem");
947 }
948 
949 
950 
952 {
953  FEMContext* fc = new FEMContext(*this);
954 
955  AutoPtr<DiffContext> ap(fc);
956 
957  DifferentiablePhysics* phys = this->get_physics();
958 
959  libmesh_assert (phys);
960 
961  // If we are solving a moving mesh problem, tell that to the Context
962  fc->set_mesh_system(phys->get_mesh_system());
963  fc->set_mesh_x_var(phys->get_mesh_x_var());
964  fc->set_mesh_y_var(phys->get_mesh_y_var());
965  fc->set_mesh_z_var(phys->get_mesh_z_var());
966 
967  ap->set_deltat_pointer( &deltat );
968 
969  // If we are solving the adjoint problem, tell that to the Context
970  ap->is_adjoint() = this->get_time_solver().is_adjoint();
971 
972  return ap;
973 }
974 
975 
976 
978 {
979  // Parent::init_context(c); // may be a good idea in derived classes
980 
981  // Although we do this in DiffSystem::build_context() and
982  // FEMSystem::build_context() as well, we do it here just to be
983  // extra sure that the deltat pointer gets set. Since the
984  // intended behavior is for classes derived from FEMSystem to
985  // call Parent::init_context() in their own init_context()
986  // overloads, we can ensure that those classes get the correct
987  // deltat pointers even if they have different build_context()
988  // overloads.
989  c.set_deltat_pointer ( &deltat );
990 
991  FEMContext &context = libmesh_cast_ref<FEMContext&>(c);
992 
993  // Make sure we're prepared to do mass integration
994  for (unsigned int var = 0; var != this->n_vars(); ++var)
995  if (this->get_physics()->is_time_evolving(var))
996  {
997  // Request shape functions based on FEType
998  switch( FEInterface::field_type( this->variable_type(var) ) )
999  {
1000  case( TYPE_SCALAR ):
1001  {
1002  FEBase* elem_fe = NULL;
1003  context.get_element_fe(var, elem_fe);
1004  elem_fe->get_JxW();
1005  elem_fe->get_phi();
1006  }
1007  break;
1008  case( TYPE_VECTOR ):
1009  {
1010  FEGenericBase<RealGradient>* elem_fe = NULL;
1011  context.get_element_fe(var, elem_fe);
1012  elem_fe->get_JxW();
1013  elem_fe->get_phi();
1014  }
1015  break;
1016  default:
1017  {
1018  libmesh_error();
1019  }
1020  }
1021  }
1022 }
1023 
1024 
1025 
1027 {
1028  // This function makes no sense unless we've already picked out some
1029  // variable(s) to reflect mesh position coordinates
1030  if (!_mesh_sys)
1031  libmesh_error();
1032 
1033  // We currently assume mesh variables are in our own system
1034  if (_mesh_sys != this)
1035  libmesh_not_implemented();
1036 
1037  // Loop over every active mesh element on this processor
1038  const MeshBase& mesh = this->get_mesh();
1039 
1042  const MeshBase::const_element_iterator end_el =
1044 
1045  AutoPtr<DiffContext> con = this->build_context();
1046  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
1047  this->init_context(_femcontext);
1048 
1049  // Get the solution's mesh variables from every element
1050  for ( ; el != end_el; ++el)
1051  {
1052  _femcontext.pre_fe_reinit(*this, *el);
1053 
1054  _femcontext.elem_position_get();
1055 
1057  this->solution->insert(_femcontext.get_elem_solution(_mesh_x_var),
1058  _femcontext.get_dof_indices(_mesh_x_var) );
1060  this->solution->insert(_femcontext.get_elem_solution(_mesh_y_var),
1061  _femcontext.get_dof_indices(_mesh_y_var));
1063  this->solution->insert(_femcontext.get_elem_solution(_mesh_z_var),
1064  _femcontext.get_dof_indices(_mesh_z_var));
1065  }
1066 
1067  this->solution->close();
1068 
1069  // And make sure the current_local_solution is up to date too
1070  this->System::update();
1071 }
1072 
1073 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo