24 #ifndef MULTIMESH_ELEMENTDATACONTAINER_H_
25 #define MULTIMESH_ELEMENTDATACONTAINER_H_
34 #include <deal.II/dofs/dof_handler.h>
35 #include <deal.II/hp/dof_handler.h>
37 using namespace dealii;
56 template<
template<
int,
int>
class DH,
typename VECTOR,
int dim>
63 "Dummy class, this constructor should never get called.",
64 "Multimesh_ElementDataContainer<dealii::DoFHandler , VECTOR, dim>::Multimesh_ElementDataContainer"));
78 template<
typename VECTOR,
int dim>
103 template<
template<
int,
int>
class FE,
typename SPARSITYPATTERN,
int dopedim,
106 const Quadrature<dim>& quad,
107 UpdateFlags update_flags,
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>*> ¶m_values,
113 const std::map<std::string, const VECTOR*> &domain_values) :
114 param_values_(param_values),
115 domain_values_(domain_values),
117 tria_element_(tria_element),
118 fine_state_fe_values_((sth.GetFESystem(
"state")), quad,
120 fine_control_fe_values_((sth.GetFESystem(
"control")), quad,
123 state_index_ = sth.GetStateIndex();
124 if (state_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);
144 ReInit(
unsigned int coarse_index,
unsigned int fine_index,
const FullMatrix<double>& prolongation_matrix);
152 GetNDoFsPerElement()
const;
156 GetMaterialId()
const;
158 GetNbrMaterialId(
unsigned int face)
const;
160 GetIsAtBoundary()
const;
162 GetElementDiameter()
const;
164 GetFEValuesState()
const;
166 GetFEValuesControl()
const;
174 GetParamValues(std::string name, Vector<double>& value)
const;
185 GetValuesState(std::string name, std::vector<double>& values)
const;
191 GetValuesState(std::string name, std::vector<Vector<double> >& values)
const;
198 GetValuesControl(std::string name, std::vector<double>& values)
const;
204 GetValuesControl(std::string name, std::vector<Vector<double> >& values)
const;
210 template<
int targetdim>
212 GetGradsState(std::string name, std::vector<Tensor<1, targetdim> >& values)
const;
218 template<
int targetdim>
220 GetGradsState(std::string name,
221 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const;
227 template<
int targetdim>
229 GetGradsControl(std::string name,
230 std::vector<Tensor<1, targetdim> >& values)
const;
236 template<
int targetdim>
238 GetGradsControl(std::string name,
239 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const;
246 GetStateIndex()
const;
248 GetControlIndex()
const;
250 GetCoarseIndex()
const {
return coarse_index_; }
252 GetFineIndex()
const {
return fine_index_; }
258 GetValues(
typename dealii::DoFHandler<dim>::cell_iterator element,
259 const FullMatrix<double>& prolongation,
261 std::vector<double>& values)
const;
267 GetValues(
typename dealii::DoFHandler<dim>::cell_iterator element,
268 const FullMatrix<double>& prolongation,
270 std::vector<Vector<double> >& values)
const;
275 template<
int targetdim>
277 GetGrads(
typename dealii::DoFHandler<dim>::cell_iterator element,
278 const FullMatrix<double>& prolongation,
280 std::vector<Tensor<1, targetdim> >& values)
const;
285 template<
int targetdim>
287 GetGrads(
typename dealii::DoFHandler<dim>::cell_iterator element,
288 const FullMatrix<double>& prolongation,
290 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const;
293 inline const std::map<std::string, const VECTOR*> &
294 GetDomainValues()
const
296 return domain_values_;
301 const std::map<std::string, const Vector<double>*> ¶m_values_;
302 const std::map<std::string, const VECTOR*> &domain_values_;
303 unsigned int state_index_;
304 unsigned int control_index_;
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_;
311 FullMatrix<double> control_prolongation_;
312 FullMatrix<double> state_prolongation_;
313 unsigned int coarse_index_, fine_index_;
315 unsigned int n_q_points_per_element_;
316 unsigned int n_dofs_per_element_;
323 template<
typename VECTOR,
int dim>
327 coarse_index_ = coarse_index;
328 fine_index_ = fine_index;
329 assert(this->GetControlIndex() < element_.size());
330 if(coarse_index == this->GetStateIndex())
332 state_prolongation_ = prolongation_matrix;
333 control_prolongation_ = IdentityMatrix(element_[this->GetControlIndex()]->get_fe().dofs_per_cell);
337 if(coarse_index == this->GetControlIndex())
339 control_prolongation_ = prolongation_matrix;
340 state_prolongation_ = IdentityMatrix(element_[this->GetStateIndex()]->get_fe().dofs_per_cell);
344 control_prolongation_ = IdentityMatrix(element_[this->GetControlIndex()]->get_fe().dofs_per_cell);
345 state_prolongation_ = IdentityMatrix(element_[this->GetStateIndex()]->get_fe().dofs_per_cell);
349 fine_state_fe_values_.reinit(tria_element_[GetFineIndex()]);
350 fine_control_fe_values_.reinit(tria_element_[GetFineIndex()]);
354 template<
typename VECTOR,
int dim>
358 return n_dofs_per_element_;
362 template<
typename VECTOR,
int dim>
366 return n_q_points_per_element_;
370 template<
typename VECTOR,
int dim>
374 return tria_element_[GetFineIndex()]->material_id();
378 template<
typename VECTOR,
int dim>
381 unsigned int face)
const
383 if (tria_element_[GetFineIndex()]->neighbor_index(face) != -1)
384 return tria_element_[GetFineIndex()]->neighbor(face)->material_id();
387 std::stringstream out;
388 out <<
"There is no neighbor with number " << face;
390 "Multimesh_ElementDataContainer::GetNbrMaterialId");
395 template<
typename VECTOR,
int dim>
399 return tria_element_[GetFineIndex()]->at_boundary();
402 template<
typename VECTOR,
int dim>
408 return element_[0]->diameter();
412 template<
typename VECTOR,
int dim>
416 return fine_state_fe_values_;
420 template<
typename VECTOR,
int dim>
424 return fine_control_fe_values_;
429 template<
typename VECTOR,
int dim>
432 std::string name, Vector<double>& value)
const
434 typename std::map<std::string, const Vector<double>*>::const_iterator it =
435 param_values_.find(name);
436 if (it == param_values_.end())
439 "Multimesh_ElementDataContainer::GetParamValues");
441 value = *(it->second);
445 template<
typename VECTOR,
int dim>
448 std::string name, std::vector<double>& values)
const
450 this->GetValues(element_[this->GetStateIndex()],state_prolongation_,this->GetFEValuesState(), name, values);
453 template<
typename VECTOR,
int dim>
456 std::string name, std::vector<Vector<double> >& values)
const
458 this->GetValues(element_[this->GetStateIndex()],state_prolongation_,this->GetFEValuesState(), name, values);
463 template<
typename VECTOR,
int dim>
466 std::string name, std::vector<double>& values)
const
468 this->GetValues(element_[this->GetControlIndex()],control_prolongation_,this->GetFEValuesControl(), name, values);
472 template<
typename VECTOR,
int dim>
475 std::string name, std::vector<Vector<double> >& values)
const
477 this->GetValues(element_[this->GetControlIndex()],control_prolongation_,this->GetFEValuesControl(), name, values);
481 template<
typename VECTOR,
int dim>
482 template<
int targetdim>
485 std::string name, std::vector<Tensor<1, targetdim> >& values)
const
487 this->GetGrads<targetdim> (element_[this->GetStateIndex()],state_prolongation_,this->GetFEValuesState(), name, values);
491 template<
typename VECTOR,
int dim>
492 template<
int targetdim>
495 std::string name, std::vector<std::vector<Tensor<1, targetdim> > >& values)
const
497 this->GetGrads<targetdim> (element_[this->GetStateIndex()],state_prolongation_,this->GetFEValuesState(), name, values);
502 template<
typename VECTOR,
int dim>
503 template<
int targetdim>
506 std::string name, std::vector<Tensor<1, targetdim> >& values)
const
508 this->GetGrads<targetdim> (element_[this->GetControlIndex()],control_prolongation_,this->GetFEValuesControl(), name, values);
512 template<
typename VECTOR,
int dim>
513 template<
int targetdim>
516 std::string name, std::vector<std::vector<Tensor<1, targetdim> > >& values)
const
518 this->GetGrads<targetdim> (element_[this->GetControlIndex()],control_prolongation_,this->GetFEValuesControl(), name, values);
523 template<
typename VECTOR,
int dim>
532 template<
typename VECTOR,
int dim>
534 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetControlIndex()
const
536 return control_index_;
540 template<
typename VECTOR,
int dim>
542 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetValues(
543 typename dealii::DoFHandler<dim>::cell_iterator element,
544 const FullMatrix<double>& prolongation,
546 std::vector<double>& values)
const
548 typename std::map<std::string, const VECTOR*>::const_iterator it =
549 this->GetDomainValues().find(name);
550 if (it == this->GetDomainValues().end())
552 throw DOpEException(
"Did not find " + name,
553 "Multimesh_ElementDataContainer::GetValues");
556 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
562 prolongation.vmult(dof_values_transformed,dof_values);
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)
570 const double value = dof_values_transformed(shape_func);
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++;
581 template<
typename VECTOR,
int dim>
583 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetValues(
584 typename dealii::DoFHandler<dim>::cell_iterator element,
585 const FullMatrix<double>& prolongation,
587 std::vector<Vector<double> >& values)
const
589 typename std::map<std::string, const VECTOR*>::const_iterator it =
590 this->GetDomainValues().find(name);
591 if (it == this->GetDomainValues().end())
593 throw DOpEException(
"Did not find " + name,
594 "Multimesh_ElementDataContainer::GetValues");
597 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
603 prolongation.vmult(dof_values_transformed,dof_values);
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);
612 for (
unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
614 const double value = dof_values_transformed(shape_func);
618 if (element->get_fe().is_primitive(shape_func))
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);
625 for (
unsigned int c=0; c<n_components; ++c)
627 if (element->get_fe().get_nonzero_components(shape_func)[c] ==
false)
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);
637 template<
typename VECTOR,
int dim>
638 template<
int targetdim>
640 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetGrads(
641 typename dealii::DoFHandler<dim>::cell_iterator element,
642 const FullMatrix<double>& prolongation,
644 std::vector<Tensor<1, targetdim> >& values)
const
646 typename std::map<std::string, const VECTOR*>::const_iterator it =
647 this->GetDomainValues().find(name);
648 if (it == this->GetDomainValues().end())
650 throw DOpEException(
"Did not find " + name,
651 "Multimesh_ElementDataContainerBase::GetGrads");
654 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
660 prolongation.vmult(dof_values_transformed,dof_values);
663 unsigned int n_quadrature_points = GetNQPoints();
664 std::fill_n (values.begin(), n_quadrature_points, Tensor<1,targetdim>());
666 for (
unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
668 const double value = dof_values_transformed(shape_func);
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++;
682 template<
typename VECTOR,
int dim>
683 template<
int targetdim>
685 Multimesh_ElementDataContainer<dealii::DoFHandler, VECTOR, dim>::GetGrads(
686 typename dealii::DoFHandler<dim>::cell_iterator element,
687 const FullMatrix<double>& prolongation,
689 std::vector<std::vector<Tensor<1, targetdim> > >& values)
const
691 typename std::map<std::string, const VECTOR*>::const_iterator it =
692 this->GetDomainValues().find(name);
693 if (it == this->GetDomainValues().end())
695 throw DOpEException(
"Did not find " + name,
696 "Multimesh_ElementDataContainerBase::GetGrads");
699 unsigned int dofs_per_element = element->get_fe().dofs_per_cell;
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);
705 prolongation.vmult(dof_values_transformed,dof_values);
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>());
713 for (
unsigned int shape_func=0; shape_func<dofs_per_element; ++shape_func)
715 const double value = dof_values_transformed(shape_func);
719 if (element->get_fe().is_primitive(shape_func))
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);
726 for (
unsigned int c=0; c<n_components; ++c)
728 if (element->get_fe().get_nonzero_components(shape_func)[c] ==
false)
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);
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: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