24 #ifndef _IntegratorMixed_H_
25 #define _IntegratorMixed_H_
27 #include <lac/vector.h>
28 #include <lac/block_sparsity_pattern.h>
29 #include <lac/block_sparse_matrix.h>
31 #include <numerics/matrix_tools.h>
32 #include <numerics/vector_tools.h>
34 #include <base/function.h>
36 #include <dofs/dof_tools.h>
38 #include <fe/mapping_q1.h>
45 #if DEAL_II_VERSION_GTE(7,3)
46 #include <base/types.h>
68 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
79 template<
typename PROBLEM>
81 template<
typename PROBLEM,
typename MATRIX>
83 template<
typename PROBLEM>
87 template<
typename PROBLEM>
89 template<
typename PROBLEM>
91 template<
typename PROBLEM>
93 template<
typename PROBLEM>
95 template<
typename PROBLEM>
97 template<
typename PROBLEM>
100 template<
typename PROBLEM>
102 template<
typename PROBLEM>
104 template<
typename PROBLEM>
106 template<
typename PROBLEM,
typename MATRIX>
109 inline void AddDomainData(std::string name,
const VECTOR* new_data);
112 inline void AddParamData(std::string name,
const dealii::Vector<SCALAR>* new_data);
118 inline const std::map<std::string, const VECTOR*>&
GetDomainData()
const;
119 inline const std::map<std::string, const dealii::Vector<SCALAR>*>&
GetParamData()
const;
122 INTEGRATORDATACONT & _idc;
124 std::map<std::string, const VECTOR*> _domain_data;
125 std::map<std::string, const dealii::Vector<SCALAR>*> _param_data;
130 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
131 int dimlow,
int dimhigh>
140 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
149 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
158 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
160 template<
typename PROBLEM>
163 bool apply_boundary_values)
168 const unsigned int dofs =
169 pde.GetBaseProblem().GetSpaceTimeHandler()->GetControlNDoFs();
171 dealii::Vector<SCALAR> local_vector(dofs);
172 const auto& dof_handler =
173 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
174 auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
175 auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
178 _idc.InitializeEDC(pde.GetUpdateFlags(),
179 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
180 this->GetParamData(), this->GetDomainData());
181 auto& edc = _idc.GetElementDataContainer();
183 bool need_faces = pde.HasFaces();
184 std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
185 bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
186 _idc.InitializeFDC(pde.GetFaceUpdateFlags(),
187 *(pde.GetBaseProblem().GetSpaceTimeHandler()),
189 this->GetParamData(),
190 this->GetDomainData());
191 auto & fdc = _idc.GetFaceDataContainer();
193 for (; element[0] != endc[0]; element[0]++)
195 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
197 if (element[dh] == endc[dh])
199 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
200 "IntegratorMixedDimensions::ComputeNonlinearResidual");
207 pde.ElementRhs(edc,local_vector, -1.);
209 if(need_boundary_integrals)
211 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
213 if (element[0]->face(face)->at_boundary()
215 (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
218 pde.BoundaryRhs(fdc,local_vector,-1.);
224 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
226 if (element[0]->neighbor_index(face) != -1)
229 pde.FaceRhs(fdc,local_vector,-1.);
234 for (
unsigned int i = 0; i < dofs; ++i)
236 residual(i) += local_vector(i);
239 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
246 pde.ElementEquation(edc, local_vector, 1., 1.);
248 for (
unsigned int i = 0; i < dofs; ++i)
250 residual(i) += local_vector(i);
253 if (pde.HasControlInDirichletData())
255 ApplyTransposedInitialBoundaryValues(pde,residual, -1.);
258 if (apply_boundary_values)
260 ApplyNewtonBoundaryValues(pde,residual);
266 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
268 template<
typename PROBLEM>
271 bool apply_boundary_values)
276 const unsigned int dofs =
277 pde.GetBaseProblem().GetSpaceTimeHandler()->GetControlNDoFs();
279 dealii::Vector<SCALAR> local_vector(dofs);
280 const auto& dof_handler =
281 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
282 auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
283 auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
286 _idc.InitializeEDC(pde.GetUpdateFlags(),
287 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
288 this->GetParamData(), this->GetDomainData());
289 auto& edc = _idc.GetElementDataContainer();
291 bool need_faces = pde.HasFaces();
292 std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
293 bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
294 _idc.InitializeFDC(pde.GetFaceUpdateFlags(),
295 *(pde.GetBaseProblem().GetSpaceTimeHandler()),
297 this->GetParamData(),
298 this->GetDomainData());
299 auto & fdc = _idc.GetFaceDataContainer();
301 for (; element[0] != endc[0]; element[0]++)
303 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
305 if (element[dh] == endc[dh])
307 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
308 "IntegratorMixedDimensions::ComputeNonlinearResidual");
315 pde.ElementRhs(edc,local_vector, 1.);
317 if(need_boundary_integrals)
319 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
321 if (element[0]->face(face)->at_boundary()
323 (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
326 pde.BoundaryRhs(fdc,local_vector,1.);
332 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
334 if (element[0]->neighbor_index(face) != -1)
337 pde.FaceRhs(fdc,local_vector,1.);
342 for (
unsigned int i = 0; i < dofs; ++i)
344 residual(i) += local_vector(i);
347 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
353 if (pde.HasControlInDirichletData())
355 ApplyTransposedInitialBoundaryValues(pde,residual, -1.);
358 if (apply_boundary_values)
360 ApplyNewtonBoundaryValues(pde,residual);
367 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
369 template<
typename PROBLEM,
typename MATRIX>
372 throw DOpEException(
"You should not use this function, try VoidLinearSolver instead.",
373 "IntegratorMixedDimensions::ComputeMatrix");
378 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
380 template<
typename PROBLEM>
384 pde.ComputeLocalControlConstraints(constraints,this->GetParamData(),this->GetDomainData());
388 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
389 int dimlow,
int dimhigh>
390 template<
typename PROBLEM>
393 dimhigh>::ComputeDomainScalar(PROBLEM& pde)
397 const unsigned int n_q_points = this->GetQuadratureFormula()->size();
399 const auto& dof_handler =
400 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
401 auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
402 auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
404 _idc.InitializeEDC(pde.GetUpdateFlags(),
405 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
406 this->GetParamData(), this->GetDomainData());
407 auto& edc = _idc.GetElementDataContainer();
410 bool need_faces = pde.HasFaces();
412 for (; element[0] != endc[0]; element[0]++)
414 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
416 if (element[dh] == endc[dh])
418 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
419 "IntegratorMixedDimensions::ComputeDomainScalar");
423 ret += pde.ElementFunctional(edc);
428 "IntegratorMixedDimensions::ComputeDomainScalar");
431 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
441 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
443 template<
typename PROBLEM>
446 if (pde.GetFEValuesNeededToBeInitialized())
448 this->InitializeFEValues();
454 ret += pde.PointFunctional(this->GetParamData(), this->GetDomainData());
461 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
463 template<
typename PROBLEM>
471 const std::vector<const DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler >*>& dof_handler =
472 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
473 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler >::active_element_iterator>
474 element(dof_handler.size());
475 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler>::active_element_iterator>
476 endc(dof_handler.size());
478 for (
unsigned int dh = 0; dh < dof_handler.size(); dh++)
480 element[dh] = dof_handler[dh]->begin_active();
481 endc[dh] = dof_handler[dh]->end();
486 pde.GetFaceUpdateFlags(),
487 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
488 this->GetParamData(), this->GetDomainData());
490 std::vector<unsigned int> boundary_functional_colors =
491 pde.GetBoundaryFunctionalColors();
492 bool need_boundary_integrals = (boundary_functional_colors.size() > 0);
493 if (!need_boundary_integrals)
496 "IntegratorMixedDimensions::ComputeBoundaryScalar");
499 for (; element[0] != endc[0]; element[0]++)
501 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
503 if (element[dh] == endc[dh])
505 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
506 "IntegratorMixedDimensions::ComputeBoundaryScalar");
510 if(need_boundary_integrals)
512 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
514 if (element[0]->face(face)->at_boundary()
516 (find(boundary_functional_colors.begin(),boundary_functional_colors.end(),
517 element[0]->face(face)->boundary_indicator()) != boundary_functional_colors.end()))
521 ret += pde.BoundaryFunctional(fdc);
525 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
535 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
537 template<
typename PROBLEM>
544 const std::vector<const DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler >*>& dof_handler =
545 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
546 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler>::active_element_iterator>
547 element(dof_handler.size());
548 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler>::active_element_iterator>
549 endc(dof_handler.size());
551 for (
unsigned int dh = 0; dh < dof_handler.size(); dh++)
553 element[dh] = dof_handler[dh]->begin_active();
554 endc[dh] = dof_handler[dh]->end();
558 pde.GetFaceUpdateFlags(),
559 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
560 this->GetParamData(), this->GetDomainData());
562 bool need_faces = pde.HasFaces();
565 throw DOpEException(
"No faces required!",
"IntegratorMixedDimensions::ComputeFaceScalar");
568 for (; element[0] != endc[0]; element[0]++)
570 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
572 if (element[dh] == endc[dh])
574 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
575 "IntegratorMixedDimensions::ComputeFaceScalar");
581 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
584 ret +=pde.FaceFunctional(fdc);
587 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
598 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
600 template<
typename PROBLEM>
604 ret = pde.AlgebraicFunctional(this->GetParamData(), this->GetDomainData());
608 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
610 template<
typename PROBLEM>
617 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
619 template<
typename PROBLEM>
627 pde.GetBaseProblem().GetSpaceTimeHandler()->GetControlNDoFs();
628 dealii::Vector<SCALAR> local_vector(dofs);
630 std::vector<unsigned int> dirichlet_colors = pde.GetTransposedDirichletColors();
631 std::vector<bool> selected_components;
632 if (dirichlet_colors.size() > 0)
634 selected_components.resize(
635 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0]->n_dofs());
636 const std::vector<Point<dimhigh> >& support_points =
637 pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapDoFToSupportPoints();
639 for (
unsigned int i = 0; i < dirichlet_colors.size(); i++)
641 unsigned int color = dirichlet_colors[i];
642 std::vector<bool> comp_mask = pde.GetTransposedDirichletCompMask(color);
643 std::vector<bool> current_comp(comp_mask.size(),
false);
644 #if DEAL_II_VERSION_GTE(7,3)
645 std::set<types::boundary_id> boundary_indicators;
647 std::set<unsigned char> boundary_indicators;
649 boundary_indicators.insert(color);
650 for (
unsigned int j = 0; j < comp_mask.size(); j++)
653 current_comp[j - 1] =
false;
656 current_comp[j] =
true;
658 DoFTools::extract_boundary_dofs(pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0]->GetDEALDoFHandler(),
659 current_comp, selected_components, boundary_indicators);
662 pde.GetTransposedDirichletValues(color,
663 this->GetParamData(),
664 this->GetDomainData());
665 for (
unsigned int k = 0; k < selected_components.size(); k++)
667 if (selected_components[k])
670 DD.
value(support_points[k], j, k, local_vector);
671 for (
unsigned int l = 0; l < dofs; ++l)
673 u(l) += scale * local_vector(l);
685 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
687 template<
typename PROBLEM>
695 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
697 template<
typename PROBLEM,
typename MATRIX>
708 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
712 const VECTOR* new_data)
714 if (_domain_data.find(name) != _domain_data.end())
716 throw DOpEException(
"Adding multiple Data with name " + name +
" is prohibited!",
717 "IntegratorMixedDimensions::AddDomainData");
719 _domain_data.insert(std::pair<std::string, const VECTOR*>(name, new_data));
724 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
729 typename std::map<std::string, const VECTOR *>::iterator it = _domain_data.find(name);
730 if (it == _domain_data.end())
732 throw DOpEException(
"Deleting Data " + name +
" is impossible! Data not found",
733 "IntegratorMixedDimensions::DeleteDomainData");
735 _domain_data.erase(it);
740 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
743 SCALAR, dimlow, dimhigh>::GetDomainData()
const
750 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
754 const dealii::Vector<
757 if (_param_data.find(name) != _param_data.end())
759 throw DOpEException(
"Adding multiple Data with name " + name +
" is prohibited!",
760 "IntegratorMixedDimensions::AddParamData");
762 _param_data.insert(std::pair<std::string,
const dealii::Vector<SCALAR>*>(name, new_data));
767 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
772 typename std::map<std::string, const dealii::Vector<SCALAR>*>::iterator it =
773 _param_data.find(name);
774 if (it == _param_data.end())
776 throw DOpEException(
"Deleting Data " + name +
" is impossible! Data not found",
777 "IntegratorMixedDimensions::DeleteParamData");
779 _param_data.erase(it);
784 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
787 INTEGRATORDATACONT, VECTOR, SCALAR, dimlow, dimhigh>::GetParamData()
const
793 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
794 int dimlow,
int dimhigh>
795 const INTEGRATORDATACONT&
void ComputeMatrix(PROBLEM &pde, MATRIX &matrix)
Definition: integratormixeddims.h:370
const std::map< std::string, const VECTOR * > & GetDomainData() const
Definition: integratormixeddims.h:743
IntegratorMixedDimensions(INTEGRATORDATACONT &idc)
Definition: integratormixeddims.h:133
const INTEGRATORDATACONT & GetIntegratorDataContainer() const
Definition: integratormixeddims.h:796
~IntegratorMixedDimensions()
Definition: integratormixeddims.h:142
const std::map< std::string, const dealii::Vector< SCALAR > * > & GetParamData() const
Definition: integratormixeddims.h:787
void ApplyNewtonBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integratormixeddims.h:688
void ComputeNonlinearRhs(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integratormixeddims.h:269
virtual void value(const dealii::Point< dealdim > &p, const unsigned int component, const unsigned int dof_number, dealii::Vector< double > &local_vector) const =0
void ReInit()
Definition: integratormixeddims.h:151
void DeleteDomainData(std::string name)
Definition: integratormixeddims.h:726
void ComputeLocalControlConstraints(PROBLEM &pde, VECTOR &constraints)
Definition: integratormixeddims.h:381
Definition: facedatacontainer.h:48
SCALAR ComputeAlgebraicScalar(PROBLEM &pde)
Definition: integratormixeddims.h:601
SCALAR ComputeDomainScalar(PROBLEM &pde)
Definition: integratormixeddims.h:393
void DeleteParamData(std::string name)
Definition: integratormixeddims.h:769
void ApplyInitialBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integratormixeddims.h:611
Definition: integratormixeddims.h:70
void ApplyTransposedInitialBoundaryValues(PROBLEM &pde, VECTOR &u, SCALAR scale)
Definition: integratormixeddims.h:621
void ComputeNonlinearResidual(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integratormixeddims.h:161
SCALAR ComputeFaceScalar(PROBLEM &pde)
Definition: integratormixeddims.h:538
Definition: transposeddirichletdatainterface.h:37
void AddDomainData(std::string name, const VECTOR *new_data)
Definition: integratormixeddims.h:710
Definition: dopeexception.h:35
SCALAR ComputeBoundaryScalar(PROBLEM &pde)
Definition: integratormixeddims.h:464
SCALAR ComputePointScalar(PROBLEM &pde)
Definition: integratormixeddims.h:444
void AddParamData(std::string name, const dealii::Vector< SCALAR > *new_data)
Definition: integratormixeddims.h:752