fem_context.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/boundary_info.h"
21 #include "libmesh/dof_map.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/fe_base.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/fem_context.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/system.h"
30 #include "libmesh/time_solver.h"
31 #include "libmesh/unsteady_solver.h" // For euler_residual
32 
33 namespace libMesh
34 {
35 
37  : DiffContext(sys),
38  _mesh_sys(NULL),
39  _mesh_x_var(0),
40  _mesh_y_var(0),
41  _mesh_z_var(0),
42  side(0), edge(0),
43  _boundary_info(sys.get_mesh().boundary_info.get()),
44  elem(NULL),
45  dim(sys.get_mesh().mesh_dimension()),
46  element_qrule(NULL), side_qrule(NULL),
47  edge_qrule(NULL)
48 {
49  // We need to know which of our variables has the hardest
50  // shape functions to numerically integrate.
51 
52  unsigned int nv = sys.n_vars();
53 
54  libmesh_assert (nv);
55  FEType hardest_fe_type = sys.variable_type(0);
56 
57  for (unsigned int i=0; i != nv; ++i)
58  {
59  FEType fe_type = sys.variable_type(i);
60 
61  // FIXME - we don't yet handle mixed finite elements from
62  // different families which require different quadrature rules
63  // libmesh_assert_equal_to (fe_type.family, hardest_fe_type.family);
64 
65  if (fe_type.order > hardest_fe_type.order)
66  hardest_fe_type = fe_type;
67  }
68 
69  // Create an adequate quadrature rule
70  element_qrule = hardest_fe_type.default_quadrature_rule
71  (dim, sys.extra_quadrature_order).release();
72  side_qrule = hardest_fe_type.default_quadrature_rule
73  (dim-1, sys.extra_quadrature_order).release();
74  if (dim == 3)
75  edge_qrule = hardest_fe_type.default_quadrature_rule
76  (1, sys.extra_quadrature_order).release();
77 
78  // Next, create finite element objects
79  // Preserving backward compatibility here for now
80  // Should move to the protected/FEAbstract interface
81  _element_fe_var.resize(nv);
82  _side_fe_var.resize(nv);
83  if (dim == 3)
84  _edge_fe_var.resize(nv);
85 
86  for (unsigned int i=0; i != nv; ++i)
87  {
88  FEType fe_type = sys.variable_type(i);
89 
90  if ( _element_fe[fe_type] == NULL )
91  {
92  _element_fe[fe_type] = FEAbstract::build(dim, fe_type).release();
93  _element_fe[fe_type]->attach_quadrature_rule(element_qrule);
94  _side_fe[fe_type] = FEAbstract::build(dim, fe_type).release();
95  _side_fe[fe_type]->attach_quadrature_rule(side_qrule);
96 
97  if (dim == 3)
98  {
99  _edge_fe[fe_type] = FEAbstract::build(dim, fe_type).release();
100  _edge_fe[fe_type]->attach_quadrature_rule(edge_qrule);
101  }
102  }
103  _element_fe_var[i] = _element_fe[fe_type];
104  _side_fe_var[i] = _side_fe[fe_type];
105  if (dim == 3)
106  _edge_fe_var[i] = _edge_fe[fe_type];
107 
108  }
109 }
110 
111 
113 {
114  // We don't want to store AutoPtrs in STL containers, but we don't
115  // want to leak memory either
116  for (std::map<FEType, FEAbstract *>::iterator i = _element_fe.begin();
117  i != _element_fe.end(); ++i)
118  delete i->second;
119  _element_fe.clear();
120 
121  for (std::map<FEType, FEAbstract *>::iterator i = _side_fe.begin();
122  i != _side_fe.end(); ++i)
123  delete i->second;
124  _side_fe.clear();
125 
126  for (std::map<FEType, FEAbstract *>::iterator i = _edge_fe.begin();
127  i != _edge_fe.end(); ++i)
128  delete i->second;
129  _edge_fe.clear();
130 
131  delete element_qrule;
132  element_qrule = NULL;
133 
134  delete side_qrule;
135  side_qrule = NULL;
136 
137  if (edge_qrule)
138  delete edge_qrule;
139  side_qrule = NULL;
140 }
141 
142 
143 
145 {
146  return _boundary_info->has_boundary_id(elem, side, id);
147 }
148 
149 
150 std::vector<boundary_id_type> FEMContext::side_boundary_ids() const
151 {
153 }
154 
155 
156 
157 Number FEMContext::interior_value(unsigned int var, unsigned int qp) const
158 {
159  Number u = 0.;
160 
161  this->interior_value( var, qp, u );
162 
163  return u;
164 }
165 
166 template<typename OutputType>
167 void FEMContext::interior_value(unsigned int var, unsigned int qp,
168  OutputType& u) const
169 {
170  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
171 
172  // Get local-to-global dof index lookup
173  libmesh_assert_greater (dof_indices.size(), var);
174  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
175  (dof_indices_var[var].size());
176 
177  // Get current local coefficients
181 
182  // Get finite element object
183  FEGenericBase<OutputShape>* fe = NULL;
184  this->get_element_fe<OutputShape>( var, fe );
185 
186  // Get shape function values at quadrature point
187  const std::vector<std::vector<OutputShape> > &phi = fe->get_phi();
188 
189  // Accumulate solution value
190  u = 0.;
191 
192  for (unsigned int l=0; l != n_dofs; l++)
193  u += phi[l][qp] * coef(l);
194 
195  return;
196 }
197 
198 
199 template<typename OutputType>
200 void FEMContext::interior_values (unsigned int var,
201  const NumericVector<Number> & _system_vector,
202  std::vector<OutputType>& u_vals) const
203 {
204  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
205 
206  // Get local-to-global dof index lookup
207  libmesh_assert_greater (dof_indices.size(), var);
208  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
209  (dof_indices_var[var].size());
210 
211  // Get current local coefficients
212  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);
213 
214  // Get the finite element object
215  FEGenericBase<OutputShape>* fe = NULL;
216  this->get_element_fe<OutputShape>( var, fe );
217 
218  // Get shape function values at quadrature point
219  const std::vector<std::vector<OutputShape> > &phi = fe->get_phi();
220 
221  // Loop over all the q_points on this element
222  for (unsigned int qp=0; qp != u_vals.size(); qp++)
223  {
224  OutputType &u = u_vals[qp];
225 
226  // Compute the value at this q_point
227  u = 0.;
228 
229  for (unsigned int l=0; l != n_dofs; l++)
230  u += phi[l][qp] * coef(l);
231  }
232 
233  return;
234 }
235 
236 Gradient FEMContext::interior_gradient(unsigned int var, unsigned int qp) const
237 {
238  Gradient du;
239 
240  this->interior_gradient( var, qp, du );
241 
242  return du;
243 }
244 
245 
246 
247 template<typename OutputType>
248 void FEMContext::interior_gradient(unsigned int var, unsigned int qp,
249  OutputType& du) const
250 {
251  typedef typename TensorTools::MakeReal
253  OutputShape;
254 
255  // Get local-to-global dof index lookup
256  libmesh_assert_greater (dof_indices.size(), var);
257  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
258  (dof_indices_var[var].size());
259 
260  // Get current local coefficients
264 
265  // Get finite element object
266  FEGenericBase<OutputShape>* fe = NULL;
267  this->get_element_fe<OutputShape>( var, fe );
268 
269  // Get shape function values at quadrature point
270  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = fe->get_dphi();
271 
272  // Accumulate solution derivatives
273  du = 0;
274 
275  for (unsigned int l=0; l != n_dofs; l++)
276  du.add_scaled(dphi[l][qp], coef(l));
277 
278  return;
279 }
280 
281 
282 
283 template<typename OutputType>
285  (unsigned int var,
286  const NumericVector<Number> & _system_vector,
287  std::vector<OutputType>& du_vals) const
288 {
289  typedef typename TensorTools::MakeReal
291  OutputShape;
292 
293  // Get local-to-global dof index lookup
294  libmesh_assert_greater (dof_indices.size(), var);
295  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
296  (dof_indices_var[var].size());
297 
298  // Get current local coefficients
299  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);
300 
301  // Get finite element object
302  FEGenericBase<OutputShape>* fe = NULL;
303  this->get_element_fe<OutputShape>( var, fe );
304 
305  // Get shape function values at quadrature point
306  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = fe->get_dphi();
307 
308  // Loop over all the q_points in this finite element
309  for (unsigned int qp=0; qp != du_vals.size(); qp++)
310  {
311  OutputType &du = du_vals[qp];
312 
313  // Compute the gradient at this q_point
314  du = 0;
315 
316  for (unsigned int l=0; l != n_dofs; l++)
317  du.add_scaled(dphi[l][qp], coef(l));
318  }
319 
320  return;
321 }
322 
323 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
324 Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) const
325 {
326  Tensor d2u;
327 
328  this->interior_hessian( var, qp, d2u );
329 
330  return d2u;
331 }
332 
333 template<typename OutputType>
334 void FEMContext::interior_hessian(unsigned int var, unsigned int qp,
335  OutputType& d2u) const
336 {
337  typedef typename TensorTools::MakeReal<
340  OutputType>::type>::type>::type
341  OutputShape;
342 
343  // Get local-to-global dof index lookup
344  libmesh_assert_greater (dof_indices.size(), var);
345  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
346  (dof_indices_var[var].size());
347 
348  // Get current local coefficients
352 
353  // Get finite element object
354  FEGenericBase<OutputShape>* fe = NULL;
355  this->get_element_fe<OutputShape>( var, fe );
356 
357  // Get shape function values at quadrature point
358  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = fe->get_d2phi();
359 
360  // Accumulate solution second derivatives
361  d2u = 0.0;
362 
363  for (unsigned int l=0; l != n_dofs; l++)
364  d2u.add_scaled(d2phi[l][qp], coef(l));
365 
366  return;
367 }
368 
369 
370 template<typename OutputType>
372  (unsigned int var,
373  const NumericVector<Number> & _system_vector,
374  std::vector<OutputType>& d2u_vals) const
375 {
376  typedef typename TensorTools::MakeReal<
379  OutputType>::type>::type>::type
380  OutputShape;
381 
382  // Get local-to-global dof index lookup
383  libmesh_assert_greater (dof_indices.size(), var);
384  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
385  (dof_indices_var[var].size());
386 
387  // Get current local coefficients
388  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);
389 
390  // Get finite element object
391  FEGenericBase<OutputShape>* fe = NULL;
392  this->get_element_fe<OutputShape>( var, fe );
393 
394  // Get shape function values at quadrature point
395  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = fe->get_d2phi();
396 
397  // Loop over all the q_points in this finite element
398  for (unsigned int qp=0; qp != d2u_vals.size(); qp++)
399  {
400  OutputType &d2u = d2u_vals[qp];
401 
402  // Compute the gradient at this q_point
403  d2u = 0;
404 
405  for (unsigned int l=0; l != n_dofs; l++)
406  d2u.add_scaled(d2phi[l][qp], coef(l));
407  }
408 
409  return;
410 }
411 
412 
413 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
414 
415 
416 template<typename OutputType>
417 void FEMContext::interior_curl(unsigned int var, unsigned int qp,
418  OutputType& curl_u) const
419 {
420  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
421 
422  // Get local-to-global dof index lookup
423  libmesh_assert_greater (dof_indices.size(), var);
424  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
425  (dof_indices_var[var].size());
426 
427  // Get current local coefficients
431 
432  // Get finite element object
433  FEGenericBase<OutputShape>* fe = NULL;
434  this->get_element_fe<OutputShape>( var, fe );
435 
436  // Get shape function values at quadrature point
437  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> > &curl_phi = fe->get_curl_phi();
438 
439  // Accumulate solution curl
440  curl_u = 0.;
441 
442  for (unsigned int l=0; l != n_dofs; l++)
443  curl_u.add_scaled(curl_phi[l][qp], coef(l));
444 
445  return;
446 }
447 
448 
449 template<typename OutputType>
450 void FEMContext::interior_div(unsigned int var, unsigned int qp,
451  OutputType& div_u) const
452 {
453  typedef typename
455  <typename TensorTools::MakeReal<OutputType>::type>::type OutputShape;
456 
457  // Get local-to-global dof index lookup
458  libmesh_assert_greater (dof_indices.size(), var);
459  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
460  (dof_indices_var[var].size());
461 
462  // Get current local coefficients
466 
467  // Get finite element object
468  FEGenericBase<OutputShape>* fe = NULL;
469  this->get_element_fe<OutputShape>( var, fe );
470 
471  // Get shape function values at quadrature point
472  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence> > &div_phi = fe->get_div_phi();
473 
474  // Accumulate solution curl
475  div_u = 0.;
476 
477  for (unsigned int l=0; l != n_dofs; l++)
478  div_u += div_phi[l][qp] * coef(l);
479 
480  return;
481 }
482 
483 
484 Number FEMContext::side_value(unsigned int var, unsigned int qp) const
485 {
486  Number u = 0.;
487 
488  this->side_value( var, qp, u );
489 
490  return u;
491 }
492 
493 
494 template<typename OutputType>
495 void FEMContext::side_value(unsigned int var, unsigned int qp,
496  OutputType& u) const
497 {
498  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
499 
500  // Get local-to-global dof index lookup
501  libmesh_assert_greater (dof_indices.size(), var);
502  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
503  (dof_indices_var[var].size());
504 
505  // Get current local coefficients
509 
510  // Get finite element object
511  FEGenericBase<OutputShape>* the_side_fe = NULL;
512  this->get_side_fe<OutputShape>( var, the_side_fe );
513 
514  // Get shape function values at quadrature point
515  const std::vector<std::vector<OutputShape> > &phi = the_side_fe->get_phi();
516 
517  // Accumulate solution value
518  u = 0.;
519 
520  for (unsigned int l=0; l != n_dofs; l++)
521  u += phi[l][qp] * coef(l);
522 
523  return;
524 }
525 
526 
527 template<typename OutputType>
529  (unsigned int var,
530  const NumericVector<Number> & _system_vector,
531  std::vector<OutputType>& u_vals) const
532 {
533  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
534 
535  // Get local-to-global dof index lookup
536  libmesh_assert_greater (dof_indices.size(), var);
537  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
538  (dof_indices_var[var].size());
539 
540  // Get current local coefficients
541  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);
542 
543  // Get the finite element object
544  FEGenericBase<OutputShape>* the_side_fe = NULL;
545  this->get_side_fe<OutputShape>( var, the_side_fe );
546 
547  // Get shape function values at quadrature point
548  const std::vector<std::vector<OutputShape> > &phi = the_side_fe->get_phi();
549 
550  // Loop over all the q_points on this element
551  for (unsigned int qp=0; qp != u_vals.size(); qp++)
552  {
553  OutputType &u = u_vals[qp];
554 
555  // Compute the value at this q_point
556  u = 0.;
557 
558  for (unsigned int l=0; l != n_dofs; l++)
559  u += phi[l][qp] * coef(l);
560  }
561 
562  return;
563 }
564 
565 Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) const
566 {
567  Gradient du;
568 
569  this->side_gradient( var, qp, du );
570 
571  return du;
572 }
573 
574 
575 template<typename OutputType>
576 void FEMContext::side_gradient(unsigned int var, unsigned int qp,
577  OutputType& du) const
578 {
579  typedef typename TensorTools::MakeReal
581  OutputShape;
582 
583  // Get local-to-global dof index lookup
584  libmesh_assert_greater (dof_indices.size(), var);
585  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
586  (dof_indices_var[var].size());
587 
588  // Get current local coefficients
592 
593  // Get finite element object
594  FEGenericBase<OutputShape>* the_side_fe = NULL;
595  this->get_side_fe<OutputShape>( var, the_side_fe );
596 
597  // Get shape function values at quadrature point
598  const std::vector<std::vector< typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = the_side_fe->get_dphi();
599 
600  // Accumulate solution derivatives
601  du = 0.;
602 
603  for (unsigned int l=0; l != n_dofs; l++)
604  du.add_scaled(dphi[l][qp], coef(l));
605 
606  return;
607 }
608 
609 
610 
611 template<typename OutputType>
613  (unsigned int var,
614  const NumericVector<Number> & _system_vector,
615  std::vector<OutputType>& du_vals) const
616 {
617  typedef typename TensorTools::MakeReal
619  OutputShape;
620 
621  // Get local-to-global dof index lookup
622  libmesh_assert_greater (dof_indices.size(), var);
623  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
624  (dof_indices_var[var].size());
625 
626  // Get current local coefficients
627  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);
628 
629  // Get finite element object
630  FEGenericBase<OutputShape>* the_side_fe = NULL;
631  this->get_side_fe<OutputShape>( var, the_side_fe );
632 
633  // Get shape function values at quadrature point
634  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = the_side_fe->get_dphi();
635 
636  // Loop over all the q_points in this finite element
637  for (unsigned int qp=0; qp != du_vals.size(); qp++)
638  {
639  OutputType &du = du_vals[qp];
640 
641  du = 0;
642 
643  // Compute the gradient at this q_point
644  for (unsigned int l=0; l != n_dofs; l++)
645  du.add_scaled(dphi[l][qp], coef(l));
646  }
647 
648  return;
649 }
650 
651 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
652 Tensor FEMContext::side_hessian(unsigned int var, unsigned int qp) const
653 {
654  Tensor d2u;
655 
656  this->side_hessian( var, qp, d2u );
657 
658  return d2u;
659 }
660 
661 template<typename OutputType>
662 void FEMContext::side_hessian(unsigned int var, unsigned int qp,
663  OutputType& d2u) const
664 {
665  typedef typename TensorTools::MakeReal<
668  OutputType>::type>::type>::type
669  OutputShape;
670 
671  // Get local-to-global dof index lookup
672  libmesh_assert_greater (dof_indices.size(), var);
673  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
674  (dof_indices_var[var].size());
675 
676  // Get current local coefficients
680 
681  // Get finite element object
682  FEGenericBase<OutputShape>* the_side_fe = NULL;
683  this->get_side_fe<OutputShape>( var, the_side_fe );
684 
685  // Get shape function values at quadrature point
686  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = the_side_fe->get_d2phi();
687 
688  // Accumulate solution second derivatives
689  d2u = 0.0;
690 
691  for (unsigned int l=0; l != n_dofs; l++)
692  d2u.add_scaled(d2phi[l][qp], coef(l));
693 
694  return;
695 }
696 
697 
698 template<typename OutputType>
700  (unsigned int var,
701  const NumericVector<Number> & _system_vector,
702  std::vector<OutputType>& d2u_vals) const
703 {
704  typedef typename TensorTools::MakeReal<
707  OutputType>::type>::type>::type
708  OutputShape;
709 
710  // Get local-to-global dof index lookup
711  libmesh_assert_greater (dof_indices.size(), var);
712  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
713  (dof_indices_var[var].size());
714 
715  // Get current local coefficients
716  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);
717 
718  // Get finite element object
719  FEGenericBase<OutputShape>* the_side_fe = NULL;
720  this->get_side_fe<OutputShape>( var, the_side_fe );
721 
722  // Get shape function values at quadrature point
723  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = the_side_fe->get_d2phi();
724 
725  // Loop over all the q_points in this finite element
726  for (unsigned int qp=0; qp != d2u_vals.size(); qp++)
727  {
728  OutputType &d2u = d2u_vals[qp];
729 
730  // Compute the gradient at this q_point
731  d2u = 0;
732 
733  for (unsigned int l=0; l != n_dofs; l++)
734  d2u.add_scaled(d2phi[l][qp], coef(l));
735  }
736 
737  return;
738 }
739 
740 
741 
742 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
743 
744 
745 
746 Number FEMContext::point_value(unsigned int var, const Point &p) const
747 {
748  Number u = 0.;
749 
750  this->point_value( var, p, u );
751 
752  return u;
753 }
754 
755 template<typename OutputType>
756 void FEMContext::point_value(unsigned int var, const Point &p,
757  OutputType& u) const
758 {
759  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
760 
761  // Get local-to-global dof index lookup
762  libmesh_assert_greater (dof_indices.size(), var);
763  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
764  (dof_indices_var[var].size());
765 
766  // Get current local coefficients
770 
771  // Get finite element object
772  FEGenericBase<OutputShape>* fe = NULL;
773  this->get_element_fe<OutputShape>( var, fe );
774 
775  // Build a FE for calculating u(p)
776  AutoPtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );
777 
778  // Get the values of the shape function derivatives
779  const std::vector<std::vector<OutputShape> >& phi = fe_new->get_phi();
780 
781  u = 0.;
782 
783  for (unsigned int l=0; l != n_dofs; l++)
784  u += phi[l][0] * coef(l);
785 
786  return;
787 }
788 
789 
790 
791 Gradient FEMContext::point_gradient(unsigned int var, const Point &p) const
792 {
793  Gradient grad_u;
794 
795  this->point_gradient( var, p, grad_u );
796 
797  return grad_u;
798 }
799 
800 
801 
802 template<typename OutputType>
803 void FEMContext::point_gradient(unsigned int var, const Point &p,
804  OutputType& grad_u) const
805 {
806  typedef typename TensorTools::MakeReal
808  OutputShape;
809 
810  // Get local-to-global dof index lookup
811  libmesh_assert_greater (dof_indices.size(), var);
812  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
813  (dof_indices_var[var].size());
814 
815  // Get current local coefficients
819 
820  // Get finite element object
821  FEGenericBase<OutputShape>* fe = NULL;
822  this->get_element_fe<OutputShape>( var, fe );
823 
824  // Build a FE for calculating u(p)
825  AutoPtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );
826 
827  // Get the values of the shape function derivatives
828  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >& dphi = fe_new->get_dphi();
829 
830  grad_u = 0.0;
831 
832  for (unsigned int l=0; l != n_dofs; l++)
833  grad_u.add_scaled(dphi[l][0], coef(l));
834 
835  return;
836 }
837 
838 
839 
840 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
841 
842 Tensor FEMContext::point_hessian(unsigned int var, const Point &p) const
843 {
844  Tensor hess_u;
845 
846  this->point_hessian( var, p, hess_u );
847 
848  return hess_u;
849 }
850 
851 
852 template<typename OutputType>
853 void FEMContext::point_hessian(unsigned int var, const Point &p,
854  OutputType& hess_u) const
855 {
856  typedef typename TensorTools::MakeReal<
859  OutputType>::type>::type>::type
860  OutputShape;
861 
862  // Get local-to-global dof index lookup
863  libmesh_assert_greater (dof_indices.size(), var);
864  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
865  (dof_indices_var[var].size());
866 
867  // Get current local coefficients
871 
872  // Get finite element object
873  FEGenericBase<OutputShape>* fe = NULL;
874  this->get_element_fe<OutputShape>( var, fe );
875 
876  // Build a FE for calculating u(p)
877  AutoPtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );
878 
879  // Get the values of the shape function derivatives
880  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >& d2phi = fe_new->get_d2phi();
881 
882  hess_u = 0.0;
883 
884  for (unsigned int l=0; l != n_dofs; l++)
885  hess_u.add_scaled(d2phi[l][0], coef(l));
886 
887  return;
888 }
889 
890 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
891 
892 
893 template<typename OutputType>
894 void FEMContext::point_curl(unsigned int var, const Point &p,
895  OutputType& curl_u) const
896 {
897  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
898 
899  // Get local-to-global dof index lookup
900  libmesh_assert_greater (dof_indices.size(), var);
901  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
902  (dof_indices_var[var].size());
903 
904  // Get current local coefficients
908 
909  // Get finite element object
910  FEGenericBase<OutputShape>* fe = NULL;
911  this->get_element_fe<OutputShape>( var, fe );
912 
913  // Build a FE for calculating u(p)
914  AutoPtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );
915 
916  // Get the values of the shape function derivatives
917  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> >& curl_phi = fe_new->get_curl_phi();
918 
919  curl_u = 0.0;
920 
921  for (unsigned int l=0; l != n_dofs; l++)
922  curl_u.add_scaled(curl_phi[l][0], coef(l));
923 
924  return;
925 }
926 
927 
928 
929 Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp) const
930 {
931  Number u = 0.;
932 
933  this->fixed_interior_value( var, qp, u );
934 
935  return u;
936 }
937 
938 
939 
940 template<typename OutputType>
941 void FEMContext::fixed_interior_value(unsigned int var, unsigned int qp,
942  OutputType& u) const
943 {
944  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
945 
946  // Get local-to-global dof index lookup
947  libmesh_assert_greater (dof_indices.size(), var);
948  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
949  (dof_indices_var[var].size());
950 
951  // Get current local coefficients
955 
956  // Get finite element object
957  FEGenericBase<OutputShape>* fe = NULL;
958  this->get_element_fe<OutputShape>( var, fe );
959 
960  // Get shape function values at quadrature point
961  const std::vector<std::vector<OutputShape> > &phi = fe->get_phi();
962 
963  // Accumulate solution value
964  u = 0.0;
965 
966  for (unsigned int l=0; l != n_dofs; l++)
967  u += phi[l][qp] * coef(l);
968 
969  return;
970 }
971 
972 
973 
974 Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp) const
975 {
976  Gradient du;
977 
978  this->fixed_interior_gradient( var, qp, du );
979 
980  return du;
981 }
982 
983 
984 template<typename OutputType>
985 void FEMContext::FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp,
986  OutputType& du) const
987 {
988  typedef typename TensorTools::MakeReal
990  OutputShape;
991 
992  // Get local-to-global dof index lookup
993  libmesh_assert_greater (dof_indices.size(), var);
994  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
995  (dof_indices_var[var].size());
996 
997  // Get current local coefficients
998  libmesh_assert_greater (elem_fixed_subsolutions.size(), var);
999  libmesh_assert(elem_fixed_subsolutions[var]);
1000  DenseSubVector<Number> &coef = *elem_fixed_subsolutions[var];
1001 
1002  // Get finite element object
1003  FEGenericBase<OutputShape>* fe = NULL;
1004  this->get_element_fe<OutputShape>( var, fe );
1005 
1006  // Get shape function values at quadrature point
1007  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = fe->get_dphi();
1008 
1009  // Accumulate solution derivatives
1010  du = 0.0;
1011 
1012  for (unsigned int l=0; l != n_dofs; l++)
1013  du.add_scaled(dphi[l][qp], coef(l));
1014 
1015  return;
1016 }
1017 
1018 
1019 
1020 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1021 Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp) const
1022 {
1023  Tensor d2u;
1024 
1025  this->fixed_interior_hessian( var, qp, d2u );
1026 
1027  return d2u;
1028 }
1029 
1030 
1031 template<typename OutputType>
1032 void FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp,
1033  OutputType& d2u) const
1034 {
1035  typedef typename TensorTools::MakeReal<
1036  typename TensorTools::DecrementRank<
1037  typename TensorTools::DecrementRank<
1038  OutputType>::type>::type>::type
1039  OutputShape;
1040 
1041  // Get local-to-global dof index lookup
1042  libmesh_assert_greater (dof_indices.size(), var);
1043  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1044  (dof_indices_var[var].size());
1045 
1046  // Get current local coefficients
1050 
1051  // Get finite element object
1052  FEGenericBase<OutputShape>* fe = NULL;
1053  this->get_element_fe<OutputShape>( var, fe );
1054 
1055  // Get shape function values at quadrature point
1056  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = fe->get_d2phi();
1057 
1058  // Accumulate solution second derivatives
1059  d2u = 0.0;
1060 
1061  for (unsigned int l=0; l != n_dofs; l++)
1062  d2u.add_scaled(d2phi[l][qp], coef(l));
1063 
1064  return;
1065 }
1066 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1067 
1068 
1069 
1070 Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp) const
1071 {
1072  Number u = 0.;
1073 
1074  this->fixed_side_value( var, qp, u );
1075 
1076  return u;
1077 }
1078 
1079 
1080 template<typename OutputType>
1081 void FEMContext::fixed_side_value(unsigned int var, unsigned int qp,
1082  OutputType& u) const
1083 {
1084  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
1085 
1086  // Get local-to-global dof index lookup
1087  libmesh_assert_greater (dof_indices.size(), var);
1088  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1089  (dof_indices_var[var].size());
1090 
1091  // Get current local coefficients
1095 
1096  // Get finite element object
1097  FEGenericBase<OutputShape>* the_side_fe = NULL;
1098  this->get_side_fe<OutputShape>( var, the_side_fe );
1099 
1100  // Get shape function values at quadrature point
1101  const std::vector<std::vector<OutputShape> > &phi = the_side_fe->get_phi();
1102 
1103  // Accumulate solution value
1104  u = 0.0;
1105 
1106  for (unsigned int l=0; l != n_dofs; l++)
1107  u += phi[l][qp] * coef(l);
1108 
1109  return;
1110 }
1111 
1112 
1113 
1114 Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp) const
1115 {
1116  Gradient du;
1117 
1118  this->fixed_side_gradient( var, qp, du );
1119 
1120  return du;
1121 }
1122 
1123 
1124 template<typename OutputType>
1125 void FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp,
1126  OutputType& du) const
1127 {
1128  typedef typename TensorTools::MakeReal
1130  OutputShape;
1131 
1132  // Get local-to-global dof index lookup
1133  libmesh_assert_greater (dof_indices.size(), var);
1134  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1135  (dof_indices_var[var].size());
1136 
1137  // Get current local coefficients
1141 
1142  // Get finite element object
1143  FEGenericBase<OutputShape>* the_side_fe = NULL;
1144  this->get_side_fe<OutputShape>( var, the_side_fe );
1145 
1146  // Get shape function values at quadrature point
1147  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = the_side_fe->get_dphi();
1148 
1149  // Accumulate solution derivatives
1150  du = 0.0;
1151 
1152  for (unsigned int l=0; l != n_dofs; l++)
1153  du.add_scaled(dphi[l][qp], coef(l));
1154 
1155  return;
1156 }
1157 
1158 
1159 
1160 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1161 Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp) const
1162 {
1163  Tensor d2u;
1164 
1165  this->fixed_side_hessian( var, qp, d2u );
1166 
1167  return d2u;
1168 }
1169 
1170 template<typename OutputType>
1171 void FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp,
1172  OutputType& d2u) const
1173 {
1174  typedef typename TensorTools::MakeReal<
1175  typename TensorTools::DecrementRank<
1176  typename TensorTools::DecrementRank<
1177  OutputType>::type>::type>::type
1178  OutputShape;
1179 
1180  // Get local-to-global dof index lookup
1181  libmesh_assert_greater (dof_indices.size(), var);
1182  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1183  (dof_indices_var[var].size());
1184 
1185  // Get current local coefficients
1189 
1190  // Get finite element object
1191  FEGenericBase<OutputShape>* the_side_fe = NULL;
1192  this->get_side_fe<OutputShape>( var, the_side_fe );
1193 
1194  // Get shape function values at quadrature point
1195  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = the_side_fe->get_d2phi();
1196 
1197  // Accumulate solution second derivatives
1198  d2u = 0.0;
1199 
1200  for (unsigned int l=0; l != n_dofs; l++)
1201  d2u.add_scaled(d2phi[l][qp], coef(l));
1202 
1203  return;
1204 }
1205 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1206 
1207 
1208 
1209 Number FEMContext::fixed_point_value(unsigned int var, const Point &p) const
1210 {
1211  Number u = 0.;
1212 
1213  this->fixed_point_value( var, p, u );
1214 
1215  return u;
1216 }
1217 
1218 template<typename OutputType>
1219 void FEMContext::fixed_point_value(unsigned int var, const Point &p,
1220  OutputType& u) const
1221 {
1222  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
1223 
1224  // Get local-to-global dof index lookup
1225  libmesh_assert_greater (dof_indices.size(), var);
1226  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1227  (dof_indices_var[var].size());
1228 
1229  // Get current local coefficients
1233 
1234  // Get finite element object
1235  FEGenericBase<OutputShape>* fe = NULL;
1236  this->get_element_fe<OutputShape>( var, fe );
1237 
1238  // Build a FE for calculating u(p)
1239  AutoPtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );
1240 
1241  // Get the values of the shape function derivatives
1242  const std::vector<std::vector<OutputShape> >& phi = fe_new->get_phi();
1243 
1244  u = 0.;
1245 
1246  for (unsigned int l=0; l != n_dofs; l++)
1247  u += phi[l][0] * coef(l);
1248 
1249  return;
1250 }
1251 
1252 
1253 
1254 Gradient FEMContext::fixed_point_gradient(unsigned int var, const Point &p) const
1255 {
1256  Gradient grad_u;
1257 
1258  this->fixed_point_gradient( var, p, grad_u );
1259 
1260  return grad_u;
1261 }
1262 
1263 
1264 
1265 template<typename OutputType>
1266 void FEMContext::fixed_point_gradient(unsigned int var, const Point &p,
1267  OutputType& grad_u) const
1268 {
1269  typedef typename TensorTools::MakeReal
1271  OutputShape;
1272 
1273  // Get local-to-global dof index lookup
1274  libmesh_assert_greater (dof_indices.size(), var);
1275  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1276  (dof_indices_var[var].size());
1277 
1278  // Get current local coefficients
1282 
1283  // Get finite element object
1284  FEGenericBase<OutputShape>* fe = NULL;
1285  this->get_element_fe<OutputShape>( var, fe );
1286 
1287  // Build a FE for calculating u(p)
1288  AutoPtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );
1289 
1290  // Get the values of the shape function derivatives
1291  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >& dphi = fe_new->get_dphi();
1292 
1293  grad_u = 0.0;
1294 
1295  for (unsigned int l=0; l != n_dofs; l++)
1296  grad_u.add_scaled(dphi[l][0], coef(l));
1297 
1298  return;
1299 }
1300 
1301 
1302 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1303 
1304 Tensor FEMContext::fixed_point_hessian(unsigned int var, const Point &p) const
1305 {
1306  Tensor hess_u;
1307 
1308  this->fixed_point_hessian( var, p, hess_u );
1309 
1310  return hess_u;
1311 }
1312 
1313 
1314 
1315 template<typename OutputType>
1316 void FEMContext::fixed_point_hessian(unsigned int var, const Point &p,
1317  OutputType& hess_u) const
1318 {
1319  typedef typename TensorTools::MakeReal<
1320  typename TensorTools::DecrementRank<
1321  typename TensorTools::DecrementRank<
1322  OutputType>::type>::type>::type
1323  OutputShape;
1324 
1325  // Get local-to-global dof index lookup
1326  libmesh_assert_greater (dof_indices.size(), var);
1327  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1328  (dof_indices_var[var].size());
1329 
1330  // Get current local coefficients
1334 
1335  // Get finite element object
1336  FEGenericBase<OutputShape>* fe = NULL;
1337  this->get_element_fe<OutputShape>( var, fe );
1338 
1339  // Build a FE for calculating u(p)
1340  AutoPtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );
1341 
1342  // Get the values of the shape function derivatives
1343  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >& d2phi = fe_new->get_d2phi();
1344 
1345  hess_u = 0.0;
1346 
1347  for (unsigned int l=0; l != n_dofs; l++)
1348  hess_u.add_scaled(d2phi[l][0], coef(l));
1349 
1350  return;
1351 }
1352 
1353 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1354 
1355 
1356 
1357 
1358 
1360 {
1361  // Update the "time" variable of this context object
1362  this->_update_time_from_system(theta);
1363 
1364  // Handle a moving element if necessary.
1365  if (_mesh_sys)
1366  // We assume that the ``default'' state
1367  // of the mesh is its final, theta=1.0
1368  // position, so we don't bother with
1369  // mesh motion in that case.
1370  if (theta != 1.0)
1371  {
1372  // FIXME - ALE is not threadsafe yet!
1373  libmesh_assert_equal_to (libMesh::n_threads(), 1);
1374 
1375  elem_position_set(theta);
1376  }
1377  elem_fe_reinit();
1378 }
1379 
1380 
1382 {
1383  // Update the "time" variable of this context object
1384  this->_update_time_from_system(theta);
1385 
1386  // Handle a moving element if necessary
1387  if (_mesh_sys)
1388  {
1389  // FIXME - not threadsafe yet!
1390  elem_position_set(theta);
1391  side_fe_reinit();
1392  }
1393 }
1394 
1395 
1397 {
1398  // Update the "time" variable of this context object
1399  this->_update_time_from_system(theta);
1400 
1401  // Handle a moving element if necessary
1402  if (_mesh_sys)
1403  {
1404  // FIXME - not threadsafe yet!
1405  elem_position_set(theta);
1406  edge_fe_reinit();
1407  }
1408 }
1409 
1410 
1412 {
1413  // Initialize all the interior FE objects on elem.
1414  // Logging of FE::reinit is done in the FE functions
1415  std::map<FEType, FEAbstract *>::iterator local_fe_end = _element_fe.end();
1416  for (std::map<FEType, FEAbstract *>::iterator i = _element_fe.begin();
1417  i != local_fe_end; ++i)
1418  {
1419  i->second->reinit(elem);
1420  }
1421 }
1422 
1423 
1425 {
1426  // Initialize all the side FE objects on elem/side.
1427  // Logging of FE::reinit is done in the FE functions
1428  std::map<FEType, FEAbstract *>::iterator local_fe_end = _side_fe.end();
1429  for (std::map<FEType, FEAbstract *>::iterator i = _side_fe.begin();
1430  i != local_fe_end; ++i)
1431  {
1432  i->second->reinit(elem, side);
1433  }
1434 }
1435 
1436 
1437 
1439 {
1440  libmesh_assert_equal_to (dim, 3);
1441 
1442  // Initialize all the interior FE objects on elem/edge.
1443  // Logging of FE::reinit is done in the FE functions
1444  std::map<FEType, FEAbstract *>::iterator local_fe_end = _edge_fe.end();
1445  for (std::map<FEType, FEAbstract *>::iterator i = _edge_fe.begin();
1446  i != local_fe_end; ++i)
1447  {
1448  i->second->edge_reinit(elem, edge);
1449  }
1450 }
1451 
1452 
1453 
1455 {
1456  // This is too expensive to call unless we've been asked to move the mesh
1458 
1459  // This will probably break with threading when two contexts are
1460  // operating on elements which share a node
1461  libmesh_assert_equal_to (libMesh::n_threads(), 1);
1462 
1463  // If the coordinate data is in our own system, it's already
1464  // been set up for us
1465 // if (_mesh_sys == this->number())
1466 // {
1467  unsigned int n_nodes = elem->n_nodes();
1468  // For simplicity we demand that mesh coordinates be stored
1469  // in a format that allows a direct copy
1471  (_element_fe_var[_mesh_x_var]->get_fe_type().family
1472  == LAGRANGE &&
1473  _element_fe_var[_mesh_x_var]->get_fe_type().order
1474  == elem->default_order()));
1476  (_element_fe_var[_mesh_y_var]->get_fe_type().family
1477  == LAGRANGE &&
1478  _element_fe_var[_mesh_y_var]->get_fe_type().order
1479  == elem->default_order()));
1481  (_element_fe_var[_mesh_z_var]->get_fe_type().family
1482  == LAGRANGE &&
1483  _element_fe_var[_mesh_z_var]->get_fe_type().order
1484  == elem->default_order()));
1485 
1486  // Get degree of freedom coefficients from point coordinates
1488  for (unsigned int i=0; i != n_nodes; ++i)
1489  (*elem_subsolutions[_mesh_x_var])(i) = elem->point(i)(0);
1490 
1492  for (unsigned int i=0; i != n_nodes; ++i)
1493  (*elem_subsolutions[_mesh_y_var])(i) = elem->point(i)(1);
1494 
1496  for (unsigned int i=0; i != n_nodes; ++i)
1497  (*elem_subsolutions[_mesh_z_var])(i) = elem->point(i)(2);
1498 // }
1499  // FIXME - If the coordinate data is not in our own system, someone
1500  // had better get around to implementing that... - RHS
1501 // else
1502 // {
1503 // libmesh_not_implemented();
1504 // }
1505 }
1506 
1507 
1508 
1509 // We can ignore the theta argument in the current use of this
1510 // function, because elem_subsolutions will already have been set to
1511 // the theta value.
1512 //
1513 // To enable loose mesh movement coupling things will need to change.
1515 {
1516  // This is too expensive to call unless we've been asked to move the mesh
1518 
1519  // This will probably break with threading when two contexts are
1520  // operating on elements which share a node
1521  libmesh_assert_equal_to (libMesh::n_threads(), 1);
1522 
1523  // If the coordinate data is in our own system, it's already
1524  // been set up for us, and we can ignore our input parameter theta
1525 // if (_mesh_sys == this->number())
1526 // {
1527  unsigned int n_nodes = elem->n_nodes();
1528  // For simplicity we demand that mesh coordinates be stored
1529  // in a format that allows a direct copy
1531  (_element_fe_var[_mesh_x_var]->get_fe_type().family
1532  == LAGRANGE &&
1533  elem_subsolutions[_mesh_x_var]->size() == n_nodes));
1535  (_element_fe_var[_mesh_y_var]->get_fe_type().family
1536  == LAGRANGE &&
1537  elem_subsolutions[_mesh_y_var]->size() == n_nodes));
1539  (_element_fe_var[_mesh_z_var]->get_fe_type().family
1540  == LAGRANGE &&
1541  elem_subsolutions[_mesh_z_var]->size() == n_nodes));
1542 
1543  // Set the new point coordinates
1545  for (unsigned int i=0; i != n_nodes; ++i)
1546  const_cast<Elem*>(elem)->point(i)(0) =
1548 
1550  for (unsigned int i=0; i != n_nodes; ++i)
1551  const_cast<Elem*>(elem)->point(i)(1) =
1553 
1555  for (unsigned int i=0; i != n_nodes; ++i)
1556  const_cast<Elem*>(elem)->point(i)(2) =
1558 // }
1559  // FIXME - If the coordinate data is not in our own system, someone
1560  // had better get around to implementing that... - RHS
1561 // else
1562 // {
1563 // libmesh_not_implemented();
1564 // }
1565 }
1566 
1567 
1568 
1569 
1570 
1571 /*
1572 void FEMContext::reinit(const FEMSystem &sys, Elem *e)
1573 {
1574  // Initialize our elem pointer, algebraic objects
1575  this->pre_fe_reinit(e);
1576 
1577  // Moving the mesh may be necessary
1578  // Reinitializing the FE objects is definitely necessary
1579  this->elem_reinit(1.);
1580 }
1581 */
1582 
1583 
1584 
1585 void FEMContext::pre_fe_reinit(const System &sys, const Elem *e)
1586 {
1587  elem = e;
1588 
1589  // Initialize the per-element data for elem.
1591  const unsigned int n_dofs = libmesh_cast_int<unsigned int>
1592  (dof_indices.size());
1593  const std::size_t n_qoi = sys.qoi.size();
1594 
1595  if (sys.use_fixed_solution)
1596  elem_fixed_solution.resize(n_dofs);
1597 
1599 
1600  // These resize calls also zero out the residual and jacobian
1601  elem_residual.resize(n_dofs);
1602  elem_jacobian.resize(n_dofs, n_dofs);
1603 
1604  elem_qoi_derivative.resize(n_qoi);
1605  elem_qoi_subderivatives.resize(n_qoi);
1606  for (std::size_t q=0; q != n_qoi; ++q)
1607  elem_qoi_derivative[q].resize(n_dofs);
1608 
1609  // Initialize the per-variable data for elem.
1610  {
1611  unsigned int sub_dofs = 0;
1612  for (unsigned int i=0; i != sys.n_vars(); ++i)
1613  {
1615 
1616  const unsigned int n_dofs_var = libmesh_cast_int<unsigned int>
1617  (dof_indices_var[i].size());
1618 
1619  elem_subsolutions[i]->reposition
1620  (sub_dofs, n_dofs_var);
1621 
1622  if (sys.use_fixed_solution)
1623  elem_fixed_subsolutions[i]->reposition
1624  (sub_dofs, n_dofs_var);
1625 
1626  elem_subresiduals[i]->reposition
1627  (sub_dofs, n_dofs_var);
1628 
1629  for (std::size_t q=0; q != n_qoi; ++q)
1630  elem_qoi_subderivatives[q][i]->reposition
1631  (sub_dofs, n_dofs_var);
1632 
1633  for (unsigned int j=0; j != i; ++j)
1634  {
1635  const unsigned int n_dofs_var_j =
1636  libmesh_cast_int<unsigned int>
1637  (dof_indices_var[j].size());
1638 
1639  elem_subjacobians[i][j]->reposition
1640  (sub_dofs, elem_subresiduals[j]->i_off(),
1641  n_dofs_var, n_dofs_var_j);
1642  elem_subjacobians[j][i]->reposition
1643  (elem_subresiduals[j]->i_off(), sub_dofs,
1644  n_dofs_var_j, n_dofs_var);
1645  }
1646  elem_subjacobians[i][i]->reposition
1647  (sub_dofs, sub_dofs,
1648  n_dofs_var,
1649  n_dofs_var);
1650  sub_dofs += n_dofs_var;
1651  }
1652  libmesh_assert_equal_to (sub_dofs, n_dofs);
1653  }
1654 
1655  // Now do the localization for the user requested vectors
1657  const DiffContext::localized_vectors_iterator localized_vec_end = localized_vectors.end();
1658 
1659  for(; localized_vec_it != localized_vec_end; ++localized_vec_it)
1660  {
1661  const NumericVector<Number>& current_localized_vector = *localized_vec_it->first;
1662  DenseVector<Number>& target_vector = localized_vec_it->second.first;
1663 
1664  current_localized_vector.get(dof_indices, target_vector.get_values());
1665 
1666  // Initialize the per-variable data for elem.
1667  unsigned int sub_dofs = 0;
1668  for (unsigned int i=0; i != sys.n_vars(); ++i)
1669  {
1670  const unsigned int n_dofs_var = libmesh_cast_int<unsigned int>
1671  (dof_indices_var[i].size());
1673 
1674  localized_vec_it->second.second[i]->reposition
1675  (sub_dofs, n_dofs_var);
1676 
1677  sub_dofs += n_dofs_var;
1678  }
1679  libmesh_assert_equal_to (sub_dofs, n_dofs);
1680  }
1681 }
1682 
1683 
1684 
1686 {
1687  // Update the "time" variable based on the value of theta. For this
1688  // to work, we need to know the value of deltat, a pointer to which is now
1689  // stored by our parent DiffContext class. Note: get_deltat_value() will
1690  // assert in debug mode if the requested pointer is NULL.
1691  const Real deltat = this->get_deltat_value();
1692 
1693  this->time = theta*(this->system_time + deltat) + (1.-theta)*this->system_time;
1694 }
1695 
1696 
1697 
1698 template<typename OutputShape>
1700  const Point &p ) const
1701 {
1702  FEType fe_type = fe->get_fe_type();
1704 
1705  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
1706  // Build a vector of point co-ordinates to send to reinit
1707  std::vector<Point> coor(1, FEInterface::inverse_map(dim, fe_type, elem, p));
1708 
1709  // Reinitialize the element and compute the shape function values at coor
1710  fe_new->reinit (elem, &coor);
1711 
1712  return fe_new;
1713 }
1714 
1715 
1716 
1717 
1718 
1719 // Instantiate member function templates
1720 template void FEMContext::interior_value<Number>(unsigned int, unsigned int, Number&) const;
1721 template void FEMContext::interior_values<Number>(unsigned int, const NumericVector<Number> &,
1722  std::vector<Number>&) const;
1723 template void FEMContext::interior_value<Gradient>(unsigned int, unsigned int, Gradient&) const;
1724 template void FEMContext::interior_values<Gradient>(unsigned int, const NumericVector<Number> &,
1725  std::vector<Gradient>&) const;
1726 
1727 template void FEMContext::interior_gradient<Gradient>(unsigned int, unsigned int, Gradient&) const;
1728 template void FEMContext::interior_gradients<Gradient>(unsigned int, const NumericVector<Number> &,
1729  std::vector<Gradient>&) const;
1730 template void FEMContext::interior_gradient<Tensor>(unsigned int, unsigned int, Tensor&) const;
1731 template void FEMContext::interior_gradients<Tensor>(unsigned int, const NumericVector<Number> &,
1732  std::vector<Tensor>&) const;
1733 
1734 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1735 template void FEMContext::interior_hessian<Tensor>(unsigned int, unsigned int, Tensor&) const;
1736 template void FEMContext::interior_hessians<Tensor>(unsigned int, const NumericVector<Number> &,
1737  std::vector<Tensor>&) const;
1738 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1739 //template void FEMContext::interior_hessian<??>(unsigned int, unsigned int, ??&) const;
1740 //template void FEMContext::interior_hessians<??>(unsigned int, const NumericVector<Number> &,
1741 // std::vector<??>&) const;
1742 #endif
1743 
1744 template void FEMContext::interior_curl<Gradient>(unsigned int, unsigned int, Gradient&) const;
1745 
1746 template void FEMContext::interior_div<Number>(unsigned int, unsigned int, Number&) const;
1747 
1748 template void FEMContext::side_value<Number>(unsigned int, unsigned int, Number&) const;
1749 template void FEMContext::side_value<Gradient>(unsigned int, unsigned int, Gradient&) const;
1750 template void FEMContext::side_values<Number>(unsigned int, const NumericVector<Number> &,
1751  std::vector<Number>&) const;
1752 template void FEMContext::side_values<Gradient>(unsigned int, const NumericVector<Number> &,
1753  std::vector<Gradient>&) const;
1754 
1755 template void FEMContext::side_gradient<Gradient>(unsigned int, unsigned int, Gradient&) const;
1756 template void FEMContext::side_gradients<Gradient>(unsigned int, const NumericVector<Number> &,
1757  std::vector<Gradient>&) const;
1758 template void FEMContext::side_gradient<Tensor>(unsigned int, unsigned int, Tensor&) const;
1759 template void FEMContext::side_gradients<Tensor>(unsigned int, const NumericVector<Number> &,
1760  std::vector<Tensor>&) const;
1761 
1762 
1763 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1764 template void FEMContext::side_hessian<Tensor>(unsigned int, unsigned int, Tensor&) const;
1765 template void FEMContext::side_hessians<Tensor>(unsigned int, const NumericVector<Number> &,
1766  std::vector<Tensor>&) const;
1767 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1768 //template void FEMContext::side_hessian<??>(unsigned int, unsigned int,
1769 // ??&) const;
1770 //template void FEMContext::side_hessians<??>(unsigned int, const NumericVector<Number> &,
1771 // std::vector<??>&) const;
1772 #endif
1773 
1774 template void FEMContext::point_value<Number>(unsigned int, const Point&, Number&) const;
1775 template void FEMContext::point_value<Gradient>(unsigned int, const Point&, Gradient&) const;
1776 
1777 template void FEMContext::point_gradient<Gradient>(unsigned int, const Point&, Gradient&) const;
1778 template void FEMContext::point_gradient<Tensor>(unsigned int, const Point&, Tensor&) const;
1779 
1780 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1781 template void FEMContext::point_hessian<Tensor>(unsigned int, const Point&, Tensor&) const;
1782 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1783 //template void FEMContext::point_hessian<??>(unsigned int, const Point&, ??&) const;
1784 #endif
1785 
1786 template void FEMContext::point_curl<Gradient>(unsigned int, const Point&, Gradient&) const;
1787 
1788 template void FEMContext::fixed_interior_value<Number>(unsigned int, unsigned int, Number&) const;
1789 template void FEMContext::fixed_interior_value<Gradient>(unsigned int, unsigned int, Gradient&) const;
1790 
1791 template void FEMContext::fixed_interior_gradient<Gradient>(unsigned int, unsigned int, Gradient&) const;
1792 template void FEMContext::fixed_interior_gradient<Tensor>(unsigned int, unsigned int, Tensor&) const;
1793 
1794 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1795 template void FEMContext::fixed_interior_hessian<Tensor>(unsigned int, unsigned int, Tensor&) const;
1796 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1797 //template void FEMContext::fixed_interior_hessian<??>(unsigned int, unsigned int, ??&) const;
1798 #endif
1799 
1800 template void FEMContext::fixed_side_value<Number>(unsigned int, unsigned int, Number&) const;
1801 template void FEMContext::fixed_side_value<Gradient>(unsigned int, unsigned int, Gradient&) const;
1802 
1803 template void FEMContext::fixed_side_gradient<Gradient>(unsigned int, unsigned int, Gradient&) const;
1804 template void FEMContext::fixed_side_gradient<Tensor>(unsigned int, unsigned int, Tensor&) const;
1805 
1806 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1807 template void FEMContext::fixed_side_hessian<Tensor>(unsigned int, unsigned int, Tensor&) const;
1808 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1809 //template void FEMContext::fixed_side_hessian<??>(unsigned int, unsigned int, ??&) const;
1810 #endif
1811 
1812 template void FEMContext::fixed_point_value<Number>(unsigned int, const Point&, Number&) const;
1813 template void FEMContext::fixed_point_value<Gradient>(unsigned int, const Point&, Gradient&) const;
1814 
1815 template void FEMContext::fixed_point_gradient<Gradient>(unsigned int, const Point&, Gradient&) const;
1816 template void FEMContext::fixed_point_gradient<Tensor>(unsigned int, const Point&, Tensor&) const;
1817 
1818 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1819 template void FEMContext::fixed_point_hessian<Tensor>(unsigned int, const Point&, Tensor&) const;
1820 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1821 //template void FEMContext::fixed_point_hessian<??>(unsigned int, const Point&, ??&) const;
1822 #endif
1823 
1826 
1827 
1828 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo