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;
124 INTEGRATORDATACONT & idc_;
126 std::map<std::string, const VECTOR*> domain_data_;
127 std::map<std::string, const dealii::Vector<SCALAR>*> param_data_;
132 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
133 int dimlow,
int dimhigh>
142 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
151 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
160 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
162 template<
typename PROBLEM>
165 bool apply_boundary_values)
170 const unsigned int dofs =
171 pde.GetBaseProblem().GetSpaceTimeHandler()->GetControlNDoFs();
173 dealii::Vector<SCALAR> local_vector(dofs);
174 const auto& dof_handler =
175 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
176 auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
177 auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
180 idc_.InitializeEDC(pde.GetUpdateFlags(),
181 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
182 this->GetParamData(), this->GetDomainData());
183 auto& edc = idc_.GetElementDataContainer();
185 bool need_faces = pde.HasFaces();
186 std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
187 bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
188 idc_.InitializeFDC(pde.GetFaceUpdateFlags(),
189 *(pde.GetBaseProblem().GetSpaceTimeHandler()),
191 this->GetParamData(),
192 this->GetDomainData());
193 auto & fdc = idc_.GetFaceDataContainer();
195 for (; element[0] != endc[0]; element[0]++)
197 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
199 if (element[dh] == endc[dh])
201 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
202 "IntegratorMixedDimensions::ComputeNonlinearResidual");
209 pde.ElementRhs(edc,local_vector, -1.);
211 if(need_boundary_integrals)
213 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
215 if (element[0]->face(face)->at_boundary()
217 (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
220 pde.BoundaryRhs(fdc,local_vector,-1.);
226 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
228 if (element[0]->neighbor_index(face) != -1)
231 pde.FaceRhs(fdc,local_vector,-1.);
236 for (
unsigned int i = 0; i < dofs; ++i)
238 residual(i) += local_vector(i);
241 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
248 pde.ElementEquation(edc, local_vector, 1., 1.);
250 for (
unsigned int i = 0; i < dofs; ++i)
252 residual(i) += local_vector(i);
255 if (pde.HasControlInDirichletData())
257 ApplyTransposedInitialBoundaryValues(pde,residual, -1.);
261 AddPresetRightHandSide(-1.,local_vector);
262 for (
unsigned int i = 0; i < dofs; ++i)
264 residual(i) += local_vector(i);
267 if (apply_boundary_values)
269 ApplyNewtonBoundaryValues(pde,residual);
275 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
277 template<
typename PROBLEM>
280 bool apply_boundary_values)
285 const unsigned int dofs =
286 pde.GetBaseProblem().GetSpaceTimeHandler()->GetControlNDoFs();
288 dealii::Vector<SCALAR> local_vector(dofs);
289 const auto& dof_handler =
290 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
291 auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
292 auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
295 idc_.InitializeEDC(pde.GetUpdateFlags(),
296 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
297 this->GetParamData(), this->GetDomainData());
298 auto& edc = idc_.GetElementDataContainer();
300 bool need_faces = pde.HasFaces();
301 std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
302 bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
303 idc_.InitializeFDC(pde.GetFaceUpdateFlags(),
304 *(pde.GetBaseProblem().GetSpaceTimeHandler()),
306 this->GetParamData(),
307 this->GetDomainData());
308 auto & fdc = idc_.GetFaceDataContainer();
310 for (; element[0] != endc[0]; element[0]++)
312 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
314 if (element[dh] == endc[dh])
316 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
317 "IntegratorMixedDimensions::ComputeNonlinearResidual");
324 pde.ElementRhs(edc,local_vector, 1.);
326 if(need_boundary_integrals)
328 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
330 if (element[0]->face(face)->at_boundary()
332 (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
335 pde.BoundaryRhs(fdc,local_vector,1.);
341 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
343 if (element[0]->neighbor_index(face) != -1)
346 pde.FaceRhs(fdc,local_vector,1.);
351 for (
unsigned int i = 0; i < dofs; ++i)
353 residual(i) += local_vector(i);
356 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
362 if (pde.HasControlInDirichletData())
364 ApplyTransposedInitialBoundaryValues(pde,residual, -1.);
368 AddPresetRightHandSide(-1.,local_vector);
369 for (
unsigned int i = 0; i < dofs; ++i)
371 residual(i) += local_vector(i);
374 if (apply_boundary_values)
376 ApplyNewtonBoundaryValues(pde,residual);
383 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
385 template<
typename PROBLEM,
typename MATRIX>
388 throw DOpEException(
"You should not use this function, try VoidLinearSolver instead.",
389 "IntegratorMixedDimensions::ComputeMatrix");
394 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
396 template<
typename PROBLEM>
400 pde.ComputeLocalControlConstraints(constraints,this->GetParamData(),this->GetDomainData());
404 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
405 int dimlow,
int dimhigh>
406 template<
typename PROBLEM>
409 dimhigh>::ComputeDomainScalar(PROBLEM& pde)
413 const unsigned int n_q_points = this->GetQuadratureFormula()->size();
415 const auto& dof_handler =
416 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
417 auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
418 auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
420 idc_.InitializeEDC(pde.GetUpdateFlags(),
421 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
422 this->GetParamData(), this->GetDomainData());
423 auto& edc = idc_.GetElementDataContainer();
425 for (; element[0] != endc[0]; element[0]++)
427 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
429 if (element[dh] == endc[dh])
431 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
432 "IntegratorMixedDimensions::ComputeDomainScalar");
436 ret += pde.ElementFunctional(edc);
438 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
448 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
450 template<
typename PROBLEM>
453 if (pde.GetFEValuesNeededToBeInitialized())
455 this->InitializeFEValues();
461 ret += pde.PointFunctional(this->GetParamData(), this->GetDomainData());
468 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
470 template<
typename PROBLEM>
478 const std::vector<const DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler >*>& dof_handler =
479 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
480 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler >::active_element_iterator>
481 element(dof_handler.size());
482 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler>::active_element_iterator>
483 endc(dof_handler.size());
485 for (
unsigned int dh = 0; dh < dof_handler.size(); dh++)
487 element[dh] = dof_handler[dh]->begin_active();
488 endc[dh] = dof_handler[dh]->end();
493 pde.GetFaceUpdateFlags(),
494 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
495 this->GetParamData(), this->GetDomainData());
497 std::vector<unsigned int> boundary_functional_colors =
498 pde.GetBoundaryFunctionalColors();
499 bool need_boundary_integrals = (boundary_functional_colors.size() > 0);
500 if (!need_boundary_integrals)
503 "IntegratorMixedDimensions::ComputeBoundaryScalar");
506 for (; element[0] != endc[0]; element[0]++)
508 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
510 if (element[dh] == endc[dh])
512 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
513 "IntegratorMixedDimensions::ComputeBoundaryScalar");
517 if(need_boundary_integrals)
519 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
521 if (element[0]->face(face)->at_boundary()
523 (find(boundary_functional_colors.begin(),boundary_functional_colors.end(),
524 element[0]->face(face)->boundary_indicator()) != boundary_functional_colors.end()))
528 ret += pde.BoundaryFunctional(fdc);
532 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
542 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
544 template<
typename PROBLEM>
551 const std::vector<const DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler >*>& dof_handler =
552 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
553 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler>::active_element_iterator>
554 element(dof_handler.size());
555 std::vector<typename DOpEWrapper::DoFHandler<dimhigh, dealii::DoFHandler>::active_element_iterator>
556 endc(dof_handler.size());
558 for (
unsigned int dh = 0; dh < dof_handler.size(); dh++)
560 element[dh] = dof_handler[dh]->begin_active();
561 endc[dh] = dof_handler[dh]->end();
565 pde.GetFaceUpdateFlags(),
566 *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
567 this->GetParamData(), this->GetDomainData());
569 bool need_faces = pde.HasFaces();
572 throw DOpEException(
"No faces required!",
"IntegratorMixedDimensions::ComputeFaceScalar");
575 for (; element[0] != endc[0]; element[0]++)
577 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
579 if (element[dh] == endc[dh])
581 throw DOpEException(
"Elementnumbers in DoFHandlers are not matching!",
582 "IntegratorMixedDimensions::ComputeFaceScalar");
588 for (
unsigned int face=0; face < dealii::GeometryInfo<dimhigh>::faces_per_cell; ++face)
591 ret +=pde.FaceFunctional(fdc);
594 for (
unsigned int dh = 1; dh < dof_handler.size(); dh++)
605 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
607 template<
typename PROBLEM>
611 ret = pde.AlgebraicFunctional(this->GetParamData(), this->GetDomainData());
615 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
617 template<
typename PROBLEM>
624 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
626 template<
typename PROBLEM>
634 pde.GetBaseProblem().GetSpaceTimeHandler()->GetControlNDoFs();
635 dealii::Vector<SCALAR> local_vector(dofs);
637 std::vector<unsigned int> dirichlet_colors = pde.GetTransposedDirichletColors();
638 std::vector<bool> selected_components;
639 if (dirichlet_colors.size() > 0)
641 selected_components.resize(
642 pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0]->n_dofs());
643 const std::vector<Point<dimhigh> >& support_points =
644 pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapDoFToSupportPoints();
646 for (
unsigned int i = 0; i < dirichlet_colors.size(); i++)
648 unsigned int color = dirichlet_colors[i];
649 std::vector<bool> comp_mask = pde.GetTransposedDirichletCompMask(color);
650 std::vector<bool> current_comp(comp_mask.size(),
false);
651 #if DEAL_II_VERSION_GTE(7,3)
652 std::set<types::boundary_id> boundary_indicators;
654 std::set<unsigned char> boundary_indicators;
656 boundary_indicators.insert(color);
657 for (
unsigned int j = 0; j < comp_mask.size(); j++)
660 current_comp[j - 1] =
false;
663 current_comp[j] =
true;
665 DoFTools::extract_boundary_dofs(pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0]->GetDEALDoFHandler(),
666 current_comp, selected_components, boundary_indicators);
669 pde.GetTransposedDirichletValues(color,
670 this->GetParamData(),
671 this->GetDomainData());
672 for (
unsigned int k = 0; k < selected_components.size(); k++)
674 if (selected_components[k])
677 DD.
value(support_points[k], j, k, local_vector);
678 for (
unsigned int l = 0; l < dofs; ++l)
680 u(l) += scale * local_vector(l);
692 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
694 template<
typename PROBLEM>
702 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
704 template<
typename PROBLEM,
typename MATRIX>
715 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
719 const VECTOR* new_data)
721 if (domain_data_.find(name) != domain_data_.end())
723 throw DOpEException(
"Adding multiple Data with name " + name +
" is prohibited!",
724 "IntegratorMixedDimensions::AddDomainData");
726 domain_data_.insert(std::pair<std::string, const VECTOR*>(name, new_data));
731 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
736 typename std::map<std::string, const VECTOR *>::iterator it = domain_data_.find(name);
737 if (it == domain_data_.end())
739 throw DOpEException(
"Deleting Data " + name +
" is impossible! Data not found",
740 "IntegratorMixedDimensions::DeleteDomainData");
742 domain_data_.erase(it);
747 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
750 SCALAR, dimlow, dimhigh>::GetDomainData()
const
757 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
761 const dealii::Vector<
764 if (param_data_.find(name) != param_data_.end())
766 throw DOpEException(
"Adding multiple Data with name " + name +
" is prohibited!",
767 "IntegratorMixedDimensions::AddParamData");
769 param_data_.insert(std::pair<std::string,
const dealii::Vector<SCALAR>*>(name, new_data));
774 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
779 typename std::map<std::string, const dealii::Vector<SCALAR>*>::iterator it =
780 param_data_.find(name);
781 if (it == param_data_.end())
783 throw DOpEException(
"Deleting Data " + name +
" is impossible! Data not found",
784 "IntegratorMixedDimensions::DeleteParamData");
786 param_data_.erase(it);
791 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
int dimlow,
794 INTEGRATORDATACONT, VECTOR, SCALAR, dimlow, dimhigh>::GetParamData()
const
800 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
801 int dimlow,
int dimhigh>
802 const INTEGRATORDATACONT&
810 template<
typename INTEGRATORDATACONT,
typename VECTOR,
typename SCALAR,
811 int dimlow,
int dimhigh>
814 dealii::Vector<SCALAR> &residual)
const
816 typename std::map<std::string, const dealii::Vector<SCALAR>*>::const_iterator it =
817 param_data_.find(
"fixed_rhs");
818 if (it != param_data_.end())
820 assert(residual.size() == it->second->size());
821 residual.add(s,*(it->second));
void ComputeMatrix(PROBLEM &pde, MATRIX &matrix)
Definition: integratormixeddims.h:386
const std::map< std::string, const VECTOR * > & GetDomainData() const
Definition: integratormixeddims.h:750
IntegratorMixedDimensions(INTEGRATORDATACONT &idc)
Definition: integratormixeddims.h:135
const INTEGRATORDATACONT & GetIntegratorDataContainer() const
Definition: integratormixeddims.h:803
~IntegratorMixedDimensions()
Definition: integratormixeddims.h:144
const std::map< std::string, const dealii::Vector< SCALAR > * > & GetParamData() const
Definition: integratormixeddims.h:794
void ApplyNewtonBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integratormixeddims.h:695
void ComputeNonlinearRhs(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integratormixeddims.h:278
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:153
void DeleteDomainData(std::string name)
Definition: integratormixeddims.h:733
void ComputeLocalControlConstraints(PROBLEM &pde, VECTOR &constraints)
Definition: integratormixeddims.h:397
Definition: facedatacontainer.h:48
SCALAR ComputeAlgebraicScalar(PROBLEM &pde)
Definition: integratormixeddims.h:608
void AddPresetRightHandSide(double s, dealii::Vector< SCALAR > &residual) const
Definition: integratormixeddims.h:813
SCALAR ComputeDomainScalar(PROBLEM &pde)
Definition: integratormixeddims.h:409
void DeleteParamData(std::string name)
Definition: integratormixeddims.h:776
void ApplyInitialBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integratormixeddims.h:618
Definition: integratormixeddims.h:70
void ApplyTransposedInitialBoundaryValues(PROBLEM &pde, VECTOR &u, SCALAR scale)
Definition: integratormixeddims.h:628
void ComputeNonlinearResidual(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integratormixeddims.h:163
SCALAR ComputeFaceScalar(PROBLEM &pde)
Definition: integratormixeddims.h:545
Definition: transposeddirichletdatainterface.h:37
void AddDomainData(std::string name, const VECTOR *new_data)
Definition: integratormixeddims.h:717
Definition: dopeexception.h:35
SCALAR ComputeBoundaryScalar(PROBLEM &pde)
Definition: integratormixeddims.h:471
SCALAR ComputePointScalar(PROBLEM &pde)
Definition: integratormixeddims.h:451
void AddParamData(std::string name, const dealii::Vector< SCALAR > *new_data)
Definition: integratormixeddims.h:759