DOpE
multimesh_elementdatacontainer.h
Go to the documentation of this file.
1 
24 #ifndef MULTIMESH_ELEMENTDATACONTAINER_H_
25 #define MULTIMESH_ELEMENTDATACONTAINER_H_
26 
27 #include "spacetimehandler.h"
28 #include "statespacetimehandler.h"
29 #include "fevalues_wrapper.h"
30 #include "dopeexception.h"
31 
32 #include <dofs/dof_handler.h>
33 #include <hp/dof_handler.h>
34 
35 using namespace dealii;
36 
37 namespace DOpE
38 {
39 
54  template<template<int, int> class DH, typename VECTOR, int dim>
56  {
57  public:
59  {
60  throw(DOpE::DOpEException(
61  "Dummy class, this constructor should never get called.",
62  "Multimesh_ElementDataContainer<dealii::DoFHandler , VECTOR, dim>::Multimesh_ElementDataContainer"));
63  }
64  ;
65  };
66 
76  template<typename VECTOR, int dim>
77  class Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>
78  {
79 
80  public:
101  template<template<int, int> class FE, typename SPARSITYPATTERN, int dopedim,
102  int dealdim>
104  const Quadrature<dim>& quad,
105  UpdateFlags update_flags,
106  SpaceTimeHandler<FE, dealii::DoFHandler, SPARSITYPATTERN,
107  VECTOR, dopedim, dealdim>& sth,
108  const typename std::vector<typename dealii::DoFHandler<dim>::cell_iterator>& element,
109  const typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& tria_element,
110  const std::map<std::string, const Vector<double>*> &param_values,
111  const std::map<std::string, const VECTOR*> &domain_values) :
112  _param_values(param_values),
113  _domain_values(domain_values),
114  _element(element),
115  _tria_element(tria_element),
116  _fine_state_fe_values((sth.GetFESystem("state")), quad,
117  update_flags),
118  _fine_control_fe_values((sth.GetFESystem("control")), quad,
119  update_flags)
120  {
121  _state_index = sth.GetStateIndex();
122  if (_state_index == 1)
123  _control_index = 0;
124  else
125  _control_index = 1;
126  _n_q_points_per_element = quad.size();
127  _n_dofs_per_element = element[0]->get_fe().dofs_per_cell;
128  _control_prolongation = IdentityMatrix(_element[this->GetControlIndex()]->get_fe().dofs_per_cell);
129  _state_prolongation = IdentityMatrix(_element[this->GetStateIndex()]->get_fe().dofs_per_cell);
130  }
131 
132 
134  {
135  }
136  /*********************************************/
137  /*
138  * This function reinits the FEValues on the actual element. Should
139  * be called prior to any of the get-functions.
140  */
141  inline void
142  ReInit(unsigned int coarse_index,unsigned int fine_index, const FullMatrix<double>& prolongation_matrix);
143 
144  /*********************************************/
149  inline unsigned int
150  GetNDoFsPerElement() const;
151  inline unsigned int
152  GetNQPoints() const;
153  inline unsigned int
154  GetMaterialId() const;
155  inline unsigned int
156  GetNbrMaterialId(unsigned int face) const;
157  inline bool
158  GetIsAtBoundary() const;
159  inline double
160  GetElementDiameter() const;
161  inline const DOpEWrapper::FEValues<dim>&
162  GetFEValuesState() const;
163  inline const DOpEWrapper::FEValues<dim>&
164  GetFEValuesControl() const;
165 
166  /**********************************************/
167  /*
168  * Looks up the given name in _parameter_data and returns the corresponding value
169  * through 'value'.
170  */
171  void
172  GetParamValues(std::string name, Vector<double>& value) const;
173 
174  /*********************************************/
179  /*
180  * Writes the values of the state variable at the quadrature points into values.
181  */
182  inline void
183  GetValuesState(std::string name, std::vector<double>& values) const;
184  /*********************************************/
185  /*
186  * Same as above for the Vector valued case.
187  */
188  inline void
189  GetValuesState(std::string name, std::vector<Vector<double> >& values) const;
190 
191  /*********************************************/
192  /*
193  * Writes the values of the control variable at the quadrature points into values
194  */
195  inline void
196  GetValuesControl(std::string name, std::vector<double>& values) const;
197  /*********************************************/
198  /*
199  * Same as above for the Vector valued case.
200  */
201  inline void
202  GetValuesControl(std::string name, std::vector<Vector<double> >& values) const;
203  /*********************************************/
204  /*
205  * Writes the values of the state gradient at the quadrature points into values.
206  */
207 
208  template<int targetdim>
209  inline void
210  GetGradsState(std::string name, std::vector<Tensor<1, targetdim> >& values) const;
211 
212  /*********************************************/
213  /*
214  * Same as above for the Vector valued case.
215  */
216  template<int targetdim>
217  inline void
218  GetGradsState(std::string name,
219  std::vector<std::vector<Tensor<1, targetdim> > >& values) const;
220 
221  /*********************************************/
222  /*
223  * Writes the values of the control gradient at the quadrature points into values.
224  */
225  template<int targetdim>
226  inline void
227  GetGradsControl(std::string name,
228  std::vector<Tensor<1, targetdim> >& values) const;
229 
230  /*********************************************/
231  /*
232  * Same as above for the Vector valued case.
233  */
234  template<int targetdim>
235  inline void
236  GetGradsControl(std::string name,
237  std::vector<std::vector<Tensor<1, targetdim> > >& values) const;
238 
239  private:
240  /*
241  * Helper Functions
242  */
243  unsigned int
244  GetStateIndex() const;
245  unsigned int
246  GetControlIndex() const;
247  unsigned int
248  GetCoarseIndex() const { return _coarse_index; }
249  unsigned int
250  GetFineIndex() const { return _fine_index; }
251  /***********************************************************/
255  inline void
256  GetValues(typename dealii::DoFHandler<dim>::cell_iterator element,
257  const FullMatrix<double>& prolongation,
258  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
259  std::vector<double>& values) const;
260  /***********************************************************/
264  inline void
265  GetValues(typename dealii::DoFHandler<dim>::cell_iterator element,
266  const FullMatrix<double>& prolongation,
267  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
268  std::vector<Vector<double> >& values) const;
269  /***********************************************************/
273  template<int targetdim>
274  inline void
275  GetGrads(typename dealii::DoFHandler<dim>::cell_iterator element,
276  const FullMatrix<double>& prolongation,
277  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
278  std::vector<Tensor<1, targetdim> >& values) const;
279  /***********************************************************/
283  template<int targetdim>
284  inline void
285  GetGrads(typename dealii::DoFHandler<dim>::cell_iterator element,
286  const FullMatrix<double>& prolongation,
287  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
288  std::vector<std::vector<Tensor<1, targetdim> > >& values) const;
289  /***********************************************************/
290 
291  inline const std::map<std::string, const VECTOR*> &
292  GetDomainValues() const
293  {
294  return _domain_values;
295  }
296 
297  /***********************************************************/
298  //"global" member data, part of every instantiation
299  const std::map<std::string, const Vector<double>*> &_param_values;
300  const std::map<std::string, const VECTOR*> &_domain_values;
301  unsigned int _state_index;
302  unsigned int _control_index;
303 
304  const typename std::vector<typename dealii::DoFHandler<dim>::cell_iterator>& _element;
305  const typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& _tria_element;
306  DOpEWrapper::FEValues<dim> _fine_state_fe_values;
307  DOpEWrapper::FEValues<dim> _fine_control_fe_values;
308 
309  FullMatrix<double> _control_prolongation;
310  FullMatrix<double> _state_prolongation;
311  unsigned int _coarse_index, _fine_index;
312 
313  unsigned int _n_q_points_per_element;
314  unsigned int _n_dofs_per_element;
315  };
316 
317  /***********************************************************************/
318  /************************IMPLEMENTATION for DoFHandler*********************************/
319  /***********************************************************************/
320 
321  template<typename VECTOR, int dim>
322  void
323  DOpE::Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::ReInit(unsigned int coarse_index,unsigned int fine_index, const FullMatrix<double>& prolongation_matrix)
324  {
325  _coarse_index = coarse_index;
326  _fine_index = fine_index;
327  assert(this->GetControlIndex() < _element.size());
328  if(coarse_index == this->GetStateIndex())
329  {
330  _state_prolongation = prolongation_matrix;
331  _control_prolongation = IdentityMatrix(_element[this->GetControlIndex()]->get_fe().dofs_per_cell);
332  }
333  else
334  {
335  if(coarse_index == this->GetControlIndex())
336  {
337  _control_prolongation = prolongation_matrix;
338  _state_prolongation = IdentityMatrix(_element[this->GetStateIndex()]->get_fe().dofs_per_cell);
339  }
340  else
341  {
342  _control_prolongation = IdentityMatrix(_element[this->GetControlIndex()]->get_fe().dofs_per_cell);
343  _state_prolongation = IdentityMatrix(_element[this->GetStateIndex()]->get_fe().dofs_per_cell);
344  _fine_index = 0;
345  }
346  }
347  _fine_state_fe_values.reinit(_tria_element[GetFineIndex()]);
348  _fine_control_fe_values.reinit(_tria_element[GetFineIndex()]);
349  }
350 
351  /***********************************************************************/
352  template<typename VECTOR, int dim>
353  unsigned int
355  {
356  return _n_dofs_per_element;
357  }
358 
359  /**********************************************/
360  template<typename VECTOR, int dim>
361  unsigned int
363  {
364  return _n_q_points_per_element;
365  }
366 
367  /**********************************************/
368  template<typename VECTOR, int dim>
369  unsigned int
371  {
372  return _tria_element[GetFineIndex()]->material_id();
373  }
374 
375  /**********************************************/
376  template<typename VECTOR, int dim>
377  unsigned int
379  unsigned int face) const
380  {
381  if (_tria_element[GetFineIndex()]->neighbor_index(face) != -1)
382  return _tria_element[GetFineIndex()]->neighbor(face)->material_id();
383  else
384  throw DOpEException("There is no neighbor with number " + face,
385  "Multimesh_ElementDataContainer::GetNbrMaterialId");
386  }
387 
388  /**********************************************/
389  template<typename VECTOR, int dim>
390  bool
392  {
393  return _tria_element[GetFineIndex()]->at_boundary();
394  }
395 
396  template<typename VECTOR, int dim>
397  double
399  {
400  //Note that we return the diameter of the _element[0], which may be the coarse one,
401  //But it is the one for which we do the computation
402  return _element[0]->diameter();
403  }
404 
405  /**********************************************/
406  template<typename VECTOR, int dim>
409  {
410  return _fine_state_fe_values;
411  }
412 
413  /**********************************************/
414  template<typename VECTOR, int dim>
417  {
418  return _fine_control_fe_values;
419  }
420 
421  /**********************************************/
422 
423  template<typename VECTOR, int dim>
424  void
426  std::string name, Vector<double>& value) const
427  {
428  typename std::map<std::string, const Vector<double>*>::const_iterator it =
429  _param_values.find(name);
430  if (it == _param_values.end())
431  {
432  throw DOpEException("Did not find " + name,
433  "Multimesh_ElementDataContainer::GetParamValues");
434  }
435  value = *(it->second);
436  }
437 
438  /*********************************************/
439  template<typename VECTOR, int dim>
440  void
442  std::string name, std::vector<double>& values) const
443  {
444  this->GetValues(_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
445  }
446  /*********************************************/
447  template<typename VECTOR, int dim>
448  void
450  std::string name, std::vector<Vector<double> >& values) const
451  {
452  this->GetValues(_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
453 
454  }
455 
456  /*********************************************/
457  template<typename VECTOR, int dim>
458  void
460  std::string name, std::vector<double>& values) const
461  {
462  this->GetValues(_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
463  }
464 
465  /*********************************************/
466  template<typename VECTOR, int dim>
467  void
469  std::string name, std::vector<Vector<double> >& values) const
470  {
471  this->GetValues(_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
472  }
473 
474  /*********************************************/
475  template<typename VECTOR, int dim>
476  template<int targetdim>
477  void
479  std::string name, std::vector<Tensor<1, targetdim> >& values) const
480  {
481  this->GetGrads<targetdim> (_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
482  }
483 
484  /*********************************************/
485  template<typename VECTOR, int dim>
486  template<int targetdim>
487  void
489  std::string name, std::vector<std::vector<Tensor<1, targetdim> > >& values) const
490  {
491  this->GetGrads<targetdim> (_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
492  }
493 
494  /***********************************************************************/
495 
496  template<typename VECTOR, int dim>
497  template<int targetdim>
498  void
500  std::string name, std::vector<Tensor<1, targetdim> >& values) const
501  {
502  this->GetGrads<targetdim> (_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
503  }
504  /***********************************************************************/
505 
506  template<typename VECTOR, int dim>
507  template<int targetdim>
508  void
510  std::string name, std::vector<std::vector<Tensor<1, targetdim> > >& values) const
511  {
512  this->GetGrads<targetdim> (_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
513  }
514 
515  /***********************************************************************/
516 
517  template<typename VECTOR, int dim>
518  unsigned int
520  {
521  return _state_index;
522  }
523 
524  /***********************************************************************/
525 
526  template<typename VECTOR, int dim>
527  unsigned int
528  Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetControlIndex() const
529  {
530  return _control_index;
531  }
532 
533  /***********************************************************************/
534  template<typename VECTOR, int dim>
535  void
536  Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetValues(
537  typename dealii::DoFHandler<dim>::cell_iterator element,
538  const FullMatrix<double>& prolongation,
539  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
540  std::vector<double>& values) const
541  {
542  typename std::map<std::string, const VECTOR*>::const_iterator it =
543  this->GetDomainValues().find(name);
544  if (it == this->GetDomainValues().end())
545  {
546  throw DOpEException("Did not find " + name,
547  "Multimesh_ElementDataContainer::GetValues");
548  }
549 
550  unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
551  //Now we get the values on the real element
552  dealii::Vector<double> dof_values(dofs_per_element);
553  dealii::Vector<double> dof_values_transformed(dofs_per_element);
554  element->get_dof_values (*(it->second), dof_values);
555  //Now compute the real values at the nodal points
556  prolongation.vmult(dof_values_transformed,dof_values);
557 
558  //Copied from deal FEValuesBase<dim,spacedim>::get_function_values
559  // see deal.II/source/fe/fe_values.cc
560  unsigned int n_quadrature_points = GetNQPoints();
561  std::fill_n (values.begin(), n_quadrature_points, 0);
562  for (unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
563  {
564  const double value = dof_values_transformed(shape_func);
565  if (value == 0.)
566  continue;
567 
568  const double *shape_value_ptr = &(fe_values.shape_value(shape_func, 0));
569  for (unsigned int point=0; point<n_quadrature_points; ++point)
570  values[point] += value * *shape_value_ptr++;
571  }
572  }
573 
574  /***********************************************************************/
575  template<typename VECTOR, int dim>
576  void
577  Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetValues(
578  typename dealii::DoFHandler<dim>::cell_iterator element,
579  const FullMatrix<double>& prolongation,
580  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
581  std::vector<Vector<double> >& values) const
582  {
583  typename std::map<std::string, const VECTOR*>::const_iterator it =
584  this->GetDomainValues().find(name);
585  if (it == this->GetDomainValues().end())
586  {
587  throw DOpEException("Did not find " + name,
588  "Multimesh_ElementDataContainer::GetValues");
589  }
590 
591  unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
592  //Now we get the values on the real element
593  dealii::Vector<double> dof_values(dofs_per_element);
594  dealii::Vector<double> dof_values_transformed(dofs_per_element);
595  element->get_dof_values (*(it->second), dof_values);
596  //Now compute the real values at the nodal points
597  prolongation.vmult(dof_values_transformed,dof_values);
598 
599  //Copied from deal FEValuesBase<dim,spacedim>::get_function_values
600  // see deal.II/source/fe/fe_values.cc
601  const unsigned int n_components = element->get_fe().n_components();
602  unsigned int n_quadrature_points = GetNQPoints();
603  for (unsigned i=0;i<values.size();++i)
604  std::fill_n (values[i].begin(), values[i].size(), 0);
605 
606  for (unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
607  {
608  const double value = dof_values_transformed(shape_func);
609  if (value == 0.)
610  continue;
611 
612  if (element->get_fe().is_primitive(shape_func))
613  {
614  const unsigned int comp = element->get_fe().system_to_component_index(shape_func).first;
615  for (unsigned int point=0; point<n_quadrature_points; ++point)
616  values[point](comp) += value * fe_values.shape_value(shape_func,point);
617  }
618  else
619  for (unsigned int c=0; c<n_components; ++c)
620  {
621  if (element->get_fe().get_nonzero_components(shape_func)[c] == false)
622  continue;
623  for (unsigned int point=0; point<n_quadrature_points; ++point)
624  values[point](c) += value * fe_values.shape_value_component(shape_func,point,c);
625  }
626  }
627  }
628 
629  /***********************************************************************/
630 
631  template<typename VECTOR, int dim>
632  template<int targetdim>
633  void
634  Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetGrads(
635  typename dealii::DoFHandler<dim>::cell_iterator element,
636  const FullMatrix<double>& prolongation,
637  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
638  std::vector<Tensor<1, targetdim> >& values) const
639  {
640  typename std::map<std::string, const VECTOR*>::const_iterator it =
641  this->GetDomainValues().find(name);
642  if (it == this->GetDomainValues().end())
643  {
644  throw DOpEException("Did not find " + name,
645  "Multimesh_ElementDataContainerBase::GetGrads");
646  }
647 
648  unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
649  //Now we get the values on the real element
650  dealii::Vector<double> dof_values(dofs_per_element);
651  dealii::Vector<double> dof_values_transformed(dofs_per_element);
652  element->get_dof_values (*(it->second), dof_values);
653  //Now compute the real values at the nodal points
654  prolongation.vmult(dof_values_transformed,dof_values);
655 
656  //Copied from deal FEValuesBase<dim,spacedim>::get_function_gradients
657  unsigned int n_quadrature_points = GetNQPoints();
658  std::fill_n (values.begin(), n_quadrature_points, Tensor<1,targetdim>());
659 
660  for (unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
661  {
662  const double value = dof_values_transformed(shape_func);
663  if (value == 0.)
664  continue;
665 
666  const Tensor<1,targetdim> *shape_gradient_ptr
667  = &(fe_values.shape_grad(shape_func,0));
668  for (unsigned int point=0; point<n_quadrature_points; ++point)
669  values[point] += value * *shape_gradient_ptr++;
670  }
671 
672  }
673 
674  /***********************************************************************/
675 
676  template<typename VECTOR, int dim>
677  template<int targetdim>
678  void
679  Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetGrads(
680  typename dealii::DoFHandler<dim>::cell_iterator element,
681  const FullMatrix<double>& prolongation,
682  const DOpEWrapper::FEValues<dim>& fe_values, std::string name,
683  std::vector<std::vector<Tensor<1, targetdim> > >& values) const
684  {
685  typename std::map<std::string, const VECTOR*>::const_iterator it =
686  this->GetDomainValues().find(name);
687  if (it == this->GetDomainValues().end())
688  {
689  throw DOpEException("Did not find " + name,
690  "Multimesh_ElementDataContainerBase::GetGrads");
691  }
692 
693  unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
694  //Now we get the values on the real element
695  dealii::Vector<double> dof_values(dofs_per_element);
696  dealii::Vector<double> dof_values_transformed(dofs_per_element);
697  element->get_dof_values (*(it->second), dof_values);
698  //Now compute the real values at the nodal points
699  prolongation.vmult(dof_values_transformed,dof_values);
700 
701  //Copied from deal FEValuesBase<dim,spacedim>::get_function_gradients
702  const unsigned int n_components = element->get_fe().n_components();
703  unsigned int n_quadrature_points = GetNQPoints();
704  for (unsigned i=0;i<values.size();++i)
705  std::fill_n (values[i].begin(), values[i].size(), Tensor<1,dim>());
706 
707  for (unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
708  {
709  const double value = dof_values_transformed(shape_func);
710  if (value == 0.)
711  continue;
712 
713  if (element->get_fe().is_primitive(shape_func))
714  {
715  const unsigned int comp = element->get_fe().system_to_component_index(shape_func).first;
716  for (unsigned int point=0; point<n_quadrature_points; ++point)
717  values[point][comp] += value * fe_values.shape_grad(shape_func,point);
718  }
719  else
720  for (unsigned int c=0; c<n_components; ++c)
721  {
722  if (element->get_fe().get_nonzero_components(shape_func)[c] == false)
723  continue;
724  for (unsigned int point=0; point<n_quadrature_points; ++point)
725  values[point][c] += value * fe_values.shape_grad_component(shape_func,point,c);
726  }
727  }
728  }
729 
730  /***********************************************************************/
731  /************************END*OF*IMPLEMENTATION**************************/
732  /***********************************************************************/
733 }//end of namespace
734 
735 #endif /* MULTIMESH_ELEMENTDATACONTAINER_H_ */
Multimesh_ElementDataContainer(const Quadrature< dim > &quad, UpdateFlags update_flags, SpaceTimeHandler< FE, dealii::DoFHandler, SPARSITYPATTERN, VECTOR, dopedim, dealdim > &sth, const typename std::vector< typename dealii::DoFHandler< dim >::cell_iterator > &element, const typename std::vector< typename dealii::Triangulation< dim >::cell_iterator > &tria_element, const std::map< std::string, const Vector< double > * > &param_values, const std::map< std::string, const VECTOR * > &domain_values)
Definition: multimesh_elementdatacontainer.h:103
Definition: multimesh_elementdatacontainer.h:55
~Multimesh_ElementDataContainer()
Definition: multimesh_elementdatacontainer.h:133
Definition: spacetimehandler.h:71
Definition: fevalues_wrapper.h:43
Definition: dopeexception.h:35
Multimesh_ElementDataContainer()
Definition: multimesh_elementdatacontainer.h:58