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