24 #ifndef MULTIMESH_ELEMENTDATACONTAINER_H_
25 #define MULTIMESH_ELEMENTDATACONTAINER_H_
32 #include <dofs/dof_handler.h>
33 #include <hp/dof_handler.h>
35 using namespace dealii;
54 template<
template<
int,
int>
class DH,
typename VECTOR,
int dim>
61 "Dummy class, this constructor should never get called.",
62 "Multimesh_ElementDataContainer<dealii::DoFHandler , VECTOR, dim>::Multimesh_ElementDataContainer"));
76 template<
typename VECTOR,
int dim>
101 template<
template<
int,
int>
class FE,
typename SPARSITYPATTERN,
int dopedim,
104 const Quadrature<dim>& quad,
105 UpdateFlags update_flags,
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>*> ¶m_values,
111 const std::map<std::string, const VECTOR*> &domain_values) :
112 _param_values(param_values),
113 _domain_values(domain_values),
115 _tria_element(tria_element),
116 _fine_state_fe_values((sth.GetFESystem(
"state")), quad,
118 _fine_control_fe_values((sth.GetFESystem(
"control")), quad,
121 _state_index = sth.GetStateIndex();
122 if (_state_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);
142 ReInit(
unsigned int coarse_index,
unsigned int fine_index,
const FullMatrix<double>& prolongation_matrix);
150 GetNDoFsPerElement()
const;
154 GetMaterialId()
const;
156 GetNbrMaterialId(
unsigned int face)
const;
158 GetIsAtBoundary()
const;
160 GetElementDiameter()
const;
162 GetFEValuesState()
const;
164 GetFEValuesControl()
const;
172 GetParamValues(std::string name, Vector<double>& value)
const;
183 GetValuesState(std::string name, std::vector<double>& values)
const;
189 GetValuesState(std::string name, std::vector<Vector<double> >& values)
const;
196 GetValuesControl(std::string name, std::vector<double>& values)
const;
202 GetValuesControl(std::string name, std::vector<Vector<double> >& values)
const;
208 template<
int targetdim>
210 GetGradsState(std::string name, std::vector<Tensor<1, targetdim> >& values)
const;
216 template<
int targetdim>
218 GetGradsState(std::string name,
219 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const;
225 template<
int targetdim>
227 GetGradsControl(std::string name,
228 std::vector<Tensor<1, targetdim> >& values)
const;
234 template<
int targetdim>
236 GetGradsControl(std::string name,
237 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const;
244 GetStateIndex()
const;
246 GetControlIndex()
const;
248 GetCoarseIndex()
const {
return _coarse_index; }
250 GetFineIndex()
const {
return _fine_index; }
256 GetValues(
typename dealii::DoFHandler<dim>::cell_iterator element,
257 const FullMatrix<double>& prolongation,
259 std::vector<double>& values)
const;
265 GetValues(
typename dealii::DoFHandler<dim>::cell_iterator element,
266 const FullMatrix<double>& prolongation,
268 std::vector<Vector<double> >& values)
const;
273 template<
int targetdim>
275 GetGrads(
typename dealii::DoFHandler<dim>::cell_iterator element,
276 const FullMatrix<double>& prolongation,
278 std::vector<Tensor<1, targetdim> >& values)
const;
283 template<
int targetdim>
285 GetGrads(
typename dealii::DoFHandler<dim>::cell_iterator element,
286 const FullMatrix<double>& prolongation,
288 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const;
291 inline const std::map<std::string, const VECTOR*> &
292 GetDomainValues()
const
294 return _domain_values;
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;
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;
309 FullMatrix<double> _control_prolongation;
310 FullMatrix<double> _state_prolongation;
311 unsigned int _coarse_index, _fine_index;
313 unsigned int _n_q_points_per_element;
314 unsigned int _n_dofs_per_element;
321 template<
typename VECTOR,
int dim>
325 _coarse_index = coarse_index;
326 _fine_index = fine_index;
327 assert(this->GetControlIndex() < _element.size());
328 if(coarse_index == this->GetStateIndex())
330 _state_prolongation = prolongation_matrix;
331 _control_prolongation = IdentityMatrix(_element[this->GetControlIndex()]->get_fe().dofs_per_cell);
335 if(coarse_index == this->GetControlIndex())
337 _control_prolongation = prolongation_matrix;
338 _state_prolongation = IdentityMatrix(_element[this->GetStateIndex()]->get_fe().dofs_per_cell);
342 _control_prolongation = IdentityMatrix(_element[this->GetControlIndex()]->get_fe().dofs_per_cell);
343 _state_prolongation = IdentityMatrix(_element[this->GetStateIndex()]->get_fe().dofs_per_cell);
347 _fine_state_fe_values.reinit(_tria_element[GetFineIndex()]);
348 _fine_control_fe_values.reinit(_tria_element[GetFineIndex()]);
352 template<
typename VECTOR,
int dim>
356 return _n_dofs_per_element;
360 template<
typename VECTOR,
int dim>
364 return _n_q_points_per_element;
368 template<
typename VECTOR,
int dim>
372 return _tria_element[GetFineIndex()]->material_id();
376 template<
typename VECTOR,
int dim>
379 unsigned int face)
const
381 if (_tria_element[GetFineIndex()]->neighbor_index(face) != -1)
382 return _tria_element[GetFineIndex()]->neighbor(face)->material_id();
384 throw DOpEException(
"There is no neighbor with number " + face,
385 "Multimesh_ElementDataContainer::GetNbrMaterialId");
389 template<
typename VECTOR,
int dim>
393 return _tria_element[GetFineIndex()]->at_boundary();
396 template<
typename VECTOR,
int dim>
402 return _element[0]->diameter();
406 template<
typename VECTOR,
int dim>
410 return _fine_state_fe_values;
414 template<
typename VECTOR,
int dim>
418 return _fine_control_fe_values;
423 template<
typename VECTOR,
int dim>
426 std::string name, Vector<double>& value)
const
428 typename std::map<std::string, const Vector<double>*>::const_iterator it =
429 _param_values.find(name);
430 if (it == _param_values.end())
433 "Multimesh_ElementDataContainer::GetParamValues");
435 value = *(it->second);
439 template<
typename VECTOR,
int dim>
442 std::string name, std::vector<double>& values)
const
444 this->GetValues(_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
447 template<
typename VECTOR,
int dim>
450 std::string name, std::vector<Vector<double> >& values)
const
452 this->GetValues(_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
457 template<
typename VECTOR,
int dim>
460 std::string name, std::vector<double>& values)
const
462 this->GetValues(_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
466 template<
typename VECTOR,
int dim>
469 std::string name, std::vector<Vector<double> >& values)
const
471 this->GetValues(_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
475 template<
typename VECTOR,
int dim>
476 template<
int targetdim>
479 std::string name, std::vector<Tensor<1, targetdim> >& values)
const
481 this->GetGrads<targetdim> (_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
485 template<
typename VECTOR,
int dim>
486 template<
int targetdim>
489 std::string name, std::vector<std::vector<Tensor<1, targetdim> > >& values)
const
491 this->GetGrads<targetdim> (_element[this->GetStateIndex()],_state_prolongation,this->GetFEValuesState(), name, values);
496 template<
typename VECTOR,
int dim>
497 template<
int targetdim>
500 std::string name, std::vector<Tensor<1, targetdim> >& values)
const
502 this->GetGrads<targetdim> (_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
506 template<
typename VECTOR,
int dim>
507 template<
int targetdim>
510 std::string name, std::vector<std::vector<Tensor<1, targetdim> > >& values)
const
512 this->GetGrads<targetdim> (_element[this->GetControlIndex()],_control_prolongation,this->GetFEValuesControl(), name, values);
517 template<
typename VECTOR,
int dim>
526 template<
typename VECTOR,
int dim>
528 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetControlIndex()
const
530 return _control_index;
534 template<
typename VECTOR,
int dim>
536 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetValues(
537 typename dealii::DoFHandler<dim>::cell_iterator element,
538 const FullMatrix<double>& prolongation,
540 std::vector<double>& values)
const
542 typename std::map<std::string, const VECTOR*>::const_iterator it =
543 this->GetDomainValues().find(name);
544 if (it == this->GetDomainValues().end())
546 throw DOpEException(
"Did not find " + name,
547 "Multimesh_ElementDataContainer::GetValues");
550 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
556 prolongation.vmult(dof_values_transformed,dof_values);
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)
564 const double value = dof_values_transformed(shape_func);
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++;
575 template<
typename VECTOR,
int dim>
577 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetValues(
578 typename dealii::DoFHandler<dim>::cell_iterator element,
579 const FullMatrix<double>& prolongation,
581 std::vector<Vector<double> >& values)
const
583 typename std::map<std::string, const VECTOR*>::const_iterator it =
584 this->GetDomainValues().find(name);
585 if (it == this->GetDomainValues().end())
587 throw DOpEException(
"Did not find " + name,
588 "Multimesh_ElementDataContainer::GetValues");
591 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
597 prolongation.vmult(dof_values_transformed,dof_values);
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);
606 for (
unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
608 const double value = dof_values_transformed(shape_func);
612 if (element->get_fe().is_primitive(shape_func))
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);
619 for (
unsigned int c=0; c<n_components; ++c)
621 if (element->get_fe().get_nonzero_components(shape_func)[c] ==
false)
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);
631 template<
typename VECTOR,
int dim>
632 template<
int targetdim>
634 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetGrads(
635 typename dealii::DoFHandler<dim>::cell_iterator element,
636 const FullMatrix<double>& prolongation,
638 std::vector<Tensor<1, targetdim> >& values)
const
640 typename std::map<std::string, const VECTOR*>::const_iterator it =
641 this->GetDomainValues().find(name);
642 if (it == this->GetDomainValues().end())
644 throw DOpEException(
"Did not find " + name,
645 "Multimesh_ElementDataContainerBase::GetGrads");
648 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
654 prolongation.vmult(dof_values_transformed,dof_values);
657 unsigned int n_quadrature_points = GetNQPoints();
658 std::fill_n (values.begin(), n_quadrature_points, Tensor<1,targetdim>());
660 for (
unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
662 const double value = dof_values_transformed(shape_func);
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++;
676 template<
typename VECTOR,
int dim>
677 template<
int targetdim>
679 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetGrads(
680 typename dealii::DoFHandler<dim>::cell_iterator element,
681 const FullMatrix<double>& prolongation,
683 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const
685 typename std::map<std::string, const VECTOR*>::const_iterator it =
686 this->GetDomainValues().find(name);
687 if (it == this->GetDomainValues().end())
689 throw DOpEException(
"Did not find " + name,
690 "Multimesh_ElementDataContainerBase::GetGrads");
693 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
699 prolongation.vmult(dof_values_transformed,dof_values);
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>());
707 for (
unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
709 const double value = dof_values_transformed(shape_func);
713 if (element->get_fe().is_primitive(shape_func))
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);
720 for (
unsigned int c=0; c<n_components; ++c)
722 if (element->get_fe().get_nonzero_components(shape_func)[c] ==
false)
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);
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 > * > ¶m_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