DOpE
integrator.h
Go to the documentation of this file.
1 
24 #ifndef Integrator_H_
25 #define Integrator_H_
26 
27 #include <deal.II/lac/vector.h>
28 #include <deal.II/lac/block_sparsity_pattern.h>
29 #include <deal.II/lac/block_sparse_matrix.h>
30 #include <deal.II/numerics/matrix_tools.h>
31 #include <deal.II/numerics/vector_tools.h>
32 #include <deal.II/base/function.h>
33 
34 #include <vector>
35 
36 #include "elementdatacontainer.h"
37 #include "facedatacontainer.h"
38 #include "dwrdatacontainer.h"
39 #include "residualestimator.h"
40 #include "dopetypes.h"
41 
42 namespace DOpE
43 {
55  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
56  int dim>
57  class Integrator
58  {
59  public:
60  Integrator(INTEGRATORDATACONT& idc1);
66  Integrator(INTEGRATORDATACONT& idc1, INTEGRATORDATACONT& idc2);
67 
68  ~Integrator();
69 
74  void
75  ReInit();
76 
93  template<typename PROBLEM>
94  void
95  ComputeNonlinearResidual(PROBLEM& pde, VECTOR &residual,
96  bool apply_boundary_values = true);
115  template<typename PROBLEM>
116  void
117  ComputeNonlinearLhs(PROBLEM& pde, VECTOR &residual,
118  bool apply_boundary_values = true);
137  template<typename PROBLEM>
138  void
139  ComputeNonlinearRhs(PROBLEM& pde, VECTOR &residual,
140  bool apply_boundary_values = true);
154  template<typename PROBLEM, typename MATRIX>
155  void
156  ComputeMatrix(PROBLEM& pde, MATRIX &matrix);
157 
170  template<typename PROBLEM>
171  void
172  ComputeNonlinearAlgebraicResidual(PROBLEM& pde, VECTOR &residual);
187  template<typename PROBLEM>
188  void
189  ComputeLocalControlConstraints(PROBLEM& pde, VECTOR &constraints);
190 
202  template<typename PROBLEM>
203  SCALAR
204  ComputeDomainScalar(PROBLEM& pde);
217  template<typename PROBLEM>
218  SCALAR
219  ComputePointScalar(PROBLEM& pde);
234  template<typename PROBLEM>
235  SCALAR
236  ComputeBoundaryScalar(PROBLEM& pde);
251  template<typename PROBLEM>
252  SCALAR
253  ComputeFaceScalar(PROBLEM& pde);
268  template<typename PROBLEM>
269  SCALAR
270  ComputeAlgebraicScalar(PROBLEM& pde);
271 
284  template<typename PROBLEM>
285  void
286  ApplyInitialBoundaryValues(PROBLEM& pde, VECTOR &u);
287 
303  template<typename PROBLEM>
304  void
305  ApplyTransposedInitialBoundaryValues(PROBLEM& pde, VECTOR &u);
319  template<typename PROBLEM>
320  void
321  ApplyNewtonBoundaryValues(PROBLEM& pde, VECTOR &u);
338  template<typename PROBLEM, typename MATRIX>
339  void
340  ApplyNewtonBoundaryValues(PROBLEM& pde, MATRIX &matrix, VECTOR &rhs,
341  VECTOR &sol);
342 
354  inline void
355  AddDomainData(std::string name, const VECTOR* new_data);
366  inline void
367  DeleteDomainData(std::string name);
378  inline void
379  AddParamData(std::string name, const dealii::Vector<SCALAR>* new_data);inline void
390  DeleteParamData(std::string name);
391 
408  template<typename PROBLEM, class STH, class EDC, class FDC>
409  void
410  ComputeRefinementIndicators(PROBLEM& pde,
411  DWRDataContainer<STH, INTEGRATORDATACONT, EDC, FDC, VECTOR>& dwrc);
423  template<typename PROBLEM>
424  void
425  ComputeRefinementIndicators(PROBLEM& pde,
426  ResidualErrorContainer<VECTOR>& dwrc);
427 
428  inline INTEGRATORDATACONT&
429  GetIntegratorDataContainer() const;
430 
431  inline INTEGRATORDATACONT&
432  GetIntegratorDataContainerFunc() const;
433 
434  protected:
440  inline const std::map<std::string, const VECTOR*>&
441  GetDomainData() const;
447  inline const std::map<std::string, const dealii::Vector<SCALAR>*>&
448  GetParamData() const;
449 
463  inline void AddPresetRightHandSide(double s, VECTOR &residual) const;
464 
465  private:
466  template<template<int, int> class DH>
467  void
468  InterpolateBoundaryValues(
469  const DOpEWrapper::Mapping<dim, DH>& mapping,
470  const DOpEWrapper::DoFHandler<dim, DH>* dof_handler,
471  const unsigned int color, const dealii::Function<dim>& function,
472  std::map<unsigned int, SCALAR>& boundary_values,
473  const std::vector<bool>& comp_mask) const;
474 
475 // /**
476 // * Given a vector of active element iterators and a facenumber, checks if the face
477 // * belongs to an 'interface' (i.e. the adjoining elements have different material ids).
478 // *
479 // * @template ELEMENTITERATOR Class of the elementiterator.
480 // *
481 // * @param element The element in question.
482 // * @param face Local number of the face for which we ask if it is
483 // * at the interface.
484 // */
485 // template<typename ELEMENTITERATOR>
486 // bool
487 // AtInterface(ELEMENTITERATOR& element, unsigned int face)
488 // {
489 // if (element[0]->neighbor_index(face) != -1)
490 // if (element[0]->material_id()
491 // != element[0]->neighbor(face)->material_id())
492 // return true;
493 // return false;
494 // }
495 
496  INTEGRATORDATACONT & idc1_;
497  INTEGRATORDATACONT & idc2_;
498 
499  std::map<std::string, const VECTOR*> domain_data_;
500  std::map<std::string, const dealii::Vector<SCALAR>*> param_data_;
501  };
502 
503  /**********************************Implementation*******************************************/
504 
505  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
506  int dim>
507  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::Integrator(
508  INTEGRATORDATACONT& idc)
509  : idc1_(idc), idc2_(idc)
510  {
511  }
512 
513  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
514  int dim>
515  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::Integrator(
516  INTEGRATORDATACONT& idc1, INTEGRATORDATACONT& idc2)
517  : idc1_(idc1), idc2_(idc2)
518  {
519  }
520 
521  /**********************************Implementation*******************************************/
522 
523  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
524  int dim>
525  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::~Integrator()
526  {
527 
528  }
529 
530  /*******************************************************************************************/
531 
532  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
533  int dim>
534  void
535  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ReInit()
536  {
537 
538  }
539 
540  /*******************************************************************************************/
541 
542  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
543  int dim>
544  template<typename PROBLEM>
545  void
546  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearResidual(
547  PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
548  {
549  residual = 0.;
550  // Begin integration
551  unsigned int dofs_per_element;
552  dealii::Vector<SCALAR> local_vector;
553  std::vector<unsigned int> local_dof_indices;
554 
555  const bool need_point_rhs = pde.HasPoints();
556 
557  const auto& dof_handler =
558  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
559  auto element =
560  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
561  auto endc =
562  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
563 
564  // Generate the data containers.
565  GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
566  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
567  this->GetParamData(), this->GetDomainData());
568  auto& edc = GetIntegratorDataContainer().GetElementDataContainer();
569 
570  bool need_faces = pde.HasFaces();
571  bool need_interfaces = pde.HasInterfaces();
572  std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
573  bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
574  GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
575  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
576  element,
577  this->GetParamData(),
578  this->GetDomainData(),
579  need_interfaces);
580  auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();
581 
582  for (; element[0] != endc[0]; element[0]++)
583  {
584  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
585  {
586  if (element[dh] == endc[dh])
587  {
588  throw DOpEException(
589  "Elementnumbers in DoFHandlers are not matching!",
590  "Integrator::ComputeNonlinearResidual");
591  }
592  }
593 
594  edc.ReInit();
595 
596  dofs_per_element = element[0]->get_fe().dofs_per_cell;
597 
598  local_vector.reinit(dofs_per_element);
599  local_vector = 0;
600 
601  local_dof_indices.resize(0);
602  local_dof_indices.resize(dofs_per_element, 0);
603 
604  //the second '1' plays only a role in the stationary case. In the non-stationary
605  //case, scale_ico is set by the time-stepping-scheme
606  pde.ElementEquation(edc, local_vector, 1., 1.);
607  pde.ElementRhs(edc, local_vector, -1.);
608 
609  if(need_boundary_integrals && element[0]->at_boundary())
610  {
611  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
612  {
613 #if DEAL_II_VERSION_GTE(8,3,0)
614  if (element[0]->face(face)->at_boundary()
615  &&
616  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
617  element[0]->face(face)->boundary_id()) != boundary_equation_colors.end()))
618 #else
619  if (element[0]->face(face)->at_boundary()
620  &&
621  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
622  element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
623 #endif
624  {
625  fdc.ReInit(face);
626  pde.BoundaryEquation(fdc,local_vector, 1., 1.);
627  pde.BoundaryRhs(fdc,local_vector,-1.);
628  }
629  }
630  }
631  if(need_faces)
632  {
633  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
634  {
635  if (element[0]->neighbor_index(face) != -1)
636  {
637  fdc.ReInit(face);
638  pde.FaceEquation(fdc, local_vector, 1., 1.);
639  pde.FaceRhs(fdc, local_vector,-1.);
640  }
641  }
642  }
643  if( need_interfaces)
644  {
645 
646  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
647  {
648  // auto face_it = element[0]->face(face);
649  // first, check if we are at an interface, i.e. not the neighbour exists and
650  // it has a different material_id than the actual element
651  if(pde.AtInterface(element, face))
652  {
653  //There exist now 3 different scenarios, given the actual element and face:
654  // The neighbour behind this face is [ more | as much | less] refined
655  // than/as the actual element. We have to distinguish here only between the case 1
656  // and the other two, because these will be distinguished in in the FaceDataContainer.
657 
658  if (element[0]->neighbor(face)->has_children())
659  {
660  //first: neighbour is finer
661 
662  for (unsigned int subface_no=0;
663  subface_no < element[0]->face(face)->n_children();
664  ++subface_no)
665  {
666  //TODO Now here we have to initialise the subface_values on the
667  // actual element and then the facevalues of the neighbours
668  fdc.ReInit(face, subface_no);
669  fdc.ReInitNbr();
670 
671  pde.InterfaceEquation(fdc, local_vector, 1., 1.);
672 
673  }
674  }
675  else
676  {
677  // either neighbor is as fine as this element or
678  // it is coarser
679 
680  fdc.ReInit(face);
681  fdc.ReInitNbr();
682  pde.InterfaceEquation(fdc, local_vector, 1., 1.);
683  }
684 
685  } //endif atinterface
686  } //endfor faces
687  } //endif need_interfaces
688  //LocalToGlobal
689  element[0]->get_dof_indices(local_dof_indices);
690  for (unsigned int i = 0; i < dofs_per_element; ++i)
691  {
692  residual(local_dof_indices[i]) += local_vector(i);
693  }
694 
695  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
696  {
697  element[dh]++;
698  }
699  }
700 
701  //check if we need the evaluation of PointRhs
702  if (need_point_rhs)
703  {
704  VECTOR point_rhs;
705  pde.PointRhs(this->GetParamData(), this->GetDomainData(), point_rhs,
706  -1.);
707  residual += point_rhs;
708  }
709  //Check if some preset righthandside exists.
710  AddPresetRightHandSide(-1.,residual);
711 
712  if (apply_boundary_values)
713  {
714  ApplyNewtonBoundaryValues(pde, residual);
715  }
716  }
717 
718  /*******************************************************************************************/
719 
720  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
721  int dim>
722  template<typename PROBLEM>
723  void
724  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearLhs(
725  PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
726  {
727  {
728 
729  residual = 0.;
730  // Begin integration
731  unsigned int dofs_per_element;
732 
733  dealii::Vector<SCALAR> local_vector;
734 
735  std::vector<unsigned int> local_dof_indices;
736 
737  const auto& dof_handler =
738  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
739  auto element =
740  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
741  auto endc =
742  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
743 
744  // Generate the data containers.
745  GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
746  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
747  this->GetParamData(), this->GetDomainData());
748  auto& edc = GetIntegratorDataContainer().GetElementDataContainer();
749 
750  bool need_faces = pde.HasFaces();
751  bool need_interfaces = pde.HasInterfaces();
752  std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
753  bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
754 
755  GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
756  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
757  element,
758  this->GetParamData(),
759  this->GetDomainData(),
760  need_interfaces);
761  auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();
762 
763  for (; element[0] != endc[0]; element[0]++)
764  {
765  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
766  {
767  if (element[dh] == endc[dh])
768  {
769  throw DOpEException(
770  "Elementnumbers in DoFHandlers are not matching!",
771  "mIntegrator::ComputeNonlinearLhs");
772  }
773  }
774 
775  edc.ReInit();
776  dofs_per_element = element[0]->get_fe().dofs_per_cell;
777 
778  local_vector.reinit(dofs_per_element);
779  local_vector = 0;
780 
781  local_dof_indices.resize(0);
782  local_dof_indices.resize(dofs_per_element, 0);
783 
784  //the second '1' plays only a role in the stationary case. In the non-stationary
785  //case, scale_ico is set by the time-stepping-scheme
786  pde.ElementEquation(edc, local_vector, 1., 1.);
787 
788  if(need_boundary_integrals)
789  {
790  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
791  {
792 #if DEAL_II_VERSION_GTE(8,3,0)
793  if (element[0]->face(face)->at_boundary()
794  &&
795  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
796  element[0]->face(face)->boundary_id()) != boundary_equation_colors.end()))
797 #else
798  if (element[0]->face(face)->at_boundary()
799  &&
800  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
801  element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
802 #endif
803  {
804  fdc.ReInit(face);
805  pde.BoundaryEquation(fdc,local_vector, 1., 1.);
806  }
807  }
808  }
809  if(need_faces)
810  {
811  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
812  {
813  if (element[0]->neighbor_index(face) != -1)
814  {
815  fdc.ReInit(face);
816  pde.FaceEquation(fdc, local_vector, 1., 1.);
817  }
818  }
819  }
820 
821  if( need_interfaces)
822  {
823 
824  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
825  {
826  //auto face_it = element[0]->face(face);
827  // first, check if we are at an interface, i.e. not the neighbour exists and
828  // it has a different material_id than the actual element
829  if(pde.AtInterface(element, face))
830  {
831  //There exist now 3 different scenarios, given the actual element and face:
832  // The neighbour behind this face is [ more | as much | less] refined
833  // than/as the actual element. We have to distinguish here only between the case 1
834  // and the other two, because these will be distinguished in in the FaceDataContainer.
835 
836  if (element[0]->neighbor(face)->has_children())
837  {
838  //first: neighbour is finer
839 
840  for (unsigned int subface_no=0;
841  subface_no < element[0]->face(face)->n_children();
842  ++subface_no)
843  {
844  //TODO Now here we have to initialise the subface_values on the
845  // actual element and then the facevalues of the neighbours
846  fdc.ReInit(face, subface_no);
847  fdc.ReInitNbr();
848 
849  pde.InterfaceEquation(fdc, local_vector, 1., 1.);
850 
851  }
852  }
853  else
854  {
855  // either neighbor is as fine as this element or
856  // it is coarser
857 
858  fdc.ReInit(face);
859  fdc.ReInitNbr();
860  pde.InterfaceEquation(fdc, local_vector, 1., 1.);
861  }
862  } //endif atinterface
863  } //endfor face
864  } //endif need_interfaces
865  //LocalToGlobal
866  element[0]->get_dof_indices(local_dof_indices);
867  for (unsigned int i = 0; i < dofs_per_element; ++i)
868  {
869  residual(local_dof_indices[i]) += local_vector(i);
870  }
871 
872  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
873  {
874  element[dh]++;
875  }
876  }
877 
878  if (apply_boundary_values)
879  {
880  ApplyNewtonBoundaryValues(pde, residual);
881  }
882  }
883  }
884 
885  /*******************************************************************************************/
886 
887  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
888  int dim>
889  template<typename PROBLEM>
890  void
891  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearRhs(
892  PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
893  {
894  residual = 0.;
895  // Begin integration
896  unsigned int dofs_per_element;
897  dealii::Vector<SCALAR> local_vector;
898  std::vector<unsigned int> local_dof_indices;
899 
900  const bool need_point_rhs = pde.HasPoints();
901 
902  const auto& dof_handler =
903  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
904  auto element =
905  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
906  auto endc =
907  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
908 
909  // Initialize the data containers.
910  GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
911  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
912  this->GetParamData(), this->GetDomainData());
913  auto& edc = GetIntegratorDataContainer().GetElementDataContainer();
914 
915 // We don't have interface terms in the Rhs! They are all to be included in the Equation!
916 // bool need_interfaces = pde.HasInterfaces();
917 // if(need_interfaces )
918 // {
919 // throw DOpEException(" Interfaces not implemented yet!",
920 // "Integrator::ComputeNonlinearRhs");
921 // }
922  bool need_faces = pde.HasFaces();
923  std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
924  bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
925 
926  GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
927  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
928  element,
929  this->GetParamData(),
930  this->GetDomainData());
931  auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();
932 
933  for (; element[0] != endc[0]; element[0]++)
934  {
935  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
936  {
937  if (element[dh] == endc[dh])
938  {
939  throw DOpEException(
940  "Elementnumbers in DoFHandlers are not matching!",
941  "Integrator::ComputeNonlinearRhs");
942  }
943  }
944 
945  edc.ReInit();
946  dofs_per_element = element[0]->get_fe().dofs_per_cell;
947 
948  local_vector.reinit(dofs_per_element);
949  local_vector = 0;
950 
951  local_dof_indices.resize(0);
952  local_dof_indices.resize(dofs_per_element, 0);
953  pde.ElementRhs(edc, local_vector, 1.);
954 
955  if(need_boundary_integrals)
956  {
957  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
958  {
959 #if DEAL_II_VERSION_GTE(8,3,0)
960  if (element[0]->face(face)->at_boundary()
961  &&
962  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
963  element[0]->face(face)->boundary_id()) != boundary_equation_colors.end()))
964 #else
965  if (element[0]->face(face)->at_boundary()
966  &&
967  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
968  element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
969 #endif
970  {
971  fdc.ReInit(face);
972  pde.BoundaryRhs(fdc,local_vector,1.);
973  }
974  }
975  }
976  if(need_faces)
977  {
978  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
979  {
980  if (element[0]->neighbor_index(face) != -1)
981  {
982  fdc.ReInit(face);
983  pde.FaceRhs(fdc, local_vector);
984  }
985  }
986  }
987  //LocalToGlobal
988  element[0]->get_dof_indices(local_dof_indices);
989  for (unsigned int i = 0; i < dofs_per_element; ++i)
990  {
991  residual(local_dof_indices[i]) += local_vector(i);
992  }
993 
994  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
995  {
996  element[dh]++;
997  }
998  }
999 
1000  //check if we need the evaluation of PointRhs
1001  if (need_point_rhs)
1002  {
1003  VECTOR point_rhs;
1004  pde.PointRhs(this->GetParamData(), this->GetDomainData(), point_rhs,
1005  1.);
1006  residual += point_rhs;
1007  }
1008  //Check if some preset righthandside exists.
1009  AddPresetRightHandSide(1.,residual);
1010 
1011  if (apply_boundary_values)
1012  {
1013  ApplyNewtonBoundaryValues(pde, residual);
1014  }
1015  }
1016 
1017  /*******************************************************************************************/
1018 
1019  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1020  int dim>
1021  template<typename PROBLEM, typename MATRIX>
1022  void
1023  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeMatrix(
1024  PROBLEM& pde, MATRIX &matrix)
1025  {
1026  matrix = 0.;
1027  // Begin integration
1028  unsigned int dofs_per_element;
1029  std::vector<unsigned int> local_dof_indices;
1030 
1031  const auto& dof_handler =
1032  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
1033  auto element =
1034  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
1035  auto endc =
1036  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
1037 
1038  GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
1039  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
1040  this->GetParamData(), this->GetDomainData());
1041  auto& edc = GetIntegratorDataContainer().GetElementDataContainer();
1042 
1043  //for the interface-case
1044  unsigned int nbr_dofs_per_element;
1045  std::vector<unsigned int> nbr_local_dof_indices;
1046 
1047  bool need_faces = pde.HasFaces();
1048  bool need_interfaces = pde.HasInterfaces();
1049  std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
1050  bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
1051 
1052  GetIntegratorDataContainer().InitializeFDC(pde.GetFaceUpdateFlags(),
1053  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
1054  element,
1055  this->GetParamData(),
1056  this->GetDomainData(),
1057  need_interfaces);
1058  auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();
1059 
1060  for (; element[0] != endc[0]; element[0]++)
1061  {
1062  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
1063  {
1064  if (element[dh] == endc[dh])
1065  {
1066  throw DOpEException(
1067  "Elementnumbers in DoFHandlers are not matching!",
1068  "Integrator::ComputeMatrix");
1069  }
1070  }
1071  edc.ReInit();
1072  dofs_per_element = element[0]->get_fe().dofs_per_cell;
1073 
1074  dealii::FullMatrix<SCALAR> local_matrix(dofs_per_element,
1075  dofs_per_element);
1076  local_matrix = 0;
1077 
1078  local_dof_indices.resize(0);
1079  local_dof_indices.resize(dofs_per_element, 0);
1080  pde.ElementMatrix(edc, local_matrix);
1081 
1082  if(need_boundary_integrals)
1083  {
1084  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1085  {
1086 #if DEAL_II_VERSION_GTE(8,3,0)
1087  if (element[0]->face(face)->at_boundary()
1088  &&
1089  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1090  element[0]->face(face)->boundary_id()) != boundary_equation_colors.end()))
1091 #else
1092  if (element[0]->face(face)->at_boundary()
1093  &&
1094  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1095  element[0]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
1096 #endif
1097  {
1098  fdc.ReInit(face);
1099  pde.BoundaryMatrix(fdc, local_matrix);
1100  }
1101  }
1102  }
1103  if(need_faces)
1104  {
1105  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1106  {
1107  if (element[0]->neighbor_index(face) != -1)
1108  {
1109  fdc.ReInit(face);
1110  pde.FaceMatrix(fdc, local_matrix);
1111  }
1112  }
1113  }
1114 
1115  if( need_interfaces)
1116  {
1117  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1118  {
1119  //auto face_it = element[0]->face(face);
1120  // first, check if we are at an interface, i.e. not the neighbour exists and
1121  // it has a different material_id than the actual element
1122  if(pde.AtInterface(element, face))
1123  {
1124  //There exist now 3 different scenarios, given the actual element and face:
1125  // The neighbour behind this face is [ more | as much | less] refined
1126  // than/as the actual element. We have to distinguish here only between the case 1
1127  // and the other two, because these will be distinguished in in the FaceDataContainer.
1128 
1129  if (element[0]->neighbor(face)->has_children())
1130  {
1131  //first: neighbour is finer
1132 
1133  for (unsigned int subface_no=0;
1134  subface_no < element[0]->face(face)->n_children();
1135  ++subface_no)
1136  {
1137  //TODO Now here we have to initialise the subface_values on the
1138  // actual element and then the facevalues of the neighbours
1139  fdc.ReInit(face, subface_no);
1140  fdc.ReInitNbr();
1141 
1142  //TODO to be swapped out?
1143  nbr_dofs_per_element = fdc.GetNbrNDoFsPerElement();
1144  nbr_local_dof_indices.resize(0);
1145  nbr_local_dof_indices.resize(nbr_dofs_per_element, 0);
1146  dealii::FullMatrix<SCALAR> local_interface_matrix(dofs_per_element,nbr_dofs_per_element );
1147  local_interface_matrix = 0;
1148 
1149  pde.InterfaceMatrix(fdc, local_interface_matrix);
1150 
1151  element[0]->get_dof_indices(local_dof_indices);
1152  element[0]->neighbor(face)->get_dof_indices(nbr_local_dof_indices);
1153 
1154  for (unsigned int i = 0; i < dofs_per_element; ++i)
1155  {
1156  for (unsigned int j = 0; j < nbr_dofs_per_element; ++j)
1157  {
1158  matrix.add(local_dof_indices[i], nbr_local_dof_indices[j],
1159  local_interface_matrix(i, j));
1160  } //endfor j
1161  } //endfor i
1162 
1163  }
1164  }
1165  else
1166  {
1167  // either neighbor is as fine as this element or it is coarser
1168  fdc.ReInit(face);
1169  fdc.ReInitNbr();
1170 
1171  //TODO to be swapped out?
1172  nbr_dofs_per_element = fdc.GetNbrNDoFsPerElement();
1173  nbr_local_dof_indices.resize(0);
1174  nbr_local_dof_indices.resize(nbr_dofs_per_element, 0);
1175  dealii::FullMatrix<SCALAR> local_interface_matrix(dofs_per_element,nbr_dofs_per_element );
1176  local_interface_matrix = 0;
1177 
1178  pde.InterfaceMatrix(fdc, local_interface_matrix);
1179 
1180  element[0]->get_dof_indices(local_dof_indices);
1181  element[0]->neighbor(face)->get_dof_indices(nbr_local_dof_indices);
1182 
1183  for (unsigned int i = 0; i < dofs_per_element; ++i)
1184  {
1185  for (unsigned int j = 0; j < nbr_dofs_per_element; ++j)
1186  {
1187  matrix.add(local_dof_indices[i], nbr_local_dof_indices[j],
1188  local_interface_matrix(i, j));
1189  } //endfor j
1190  } //endfor i
1191  }
1192  } //endif atinterface
1193 
1194  }
1195  } //endif need_interfaces
1196 
1197  //LocalToGlobal
1198  element[0]->get_dof_indices(local_dof_indices);
1199  for (unsigned int i = 0; i < dofs_per_element; ++i)
1200  {
1201  for (unsigned int j = 0; j < dofs_per_element; ++j)
1202  {
1203  matrix.add(local_dof_indices[i], local_dof_indices[j],
1204  local_matrix(i, j));
1205  }
1206  }
1207 
1208  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
1209  {
1210  element[dh]++;
1211  }
1212  }
1213  }
1214 
1215  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1216  int dim>
1217  template<typename PROBLEM>
1218  SCALAR
1219  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeDomainScalar(
1220  PROBLEM& pde)
1221  {
1222  {
1223  SCALAR ret = 0.;
1224 
1225  const auto& dof_handler =
1226  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
1227  auto element =
1228  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
1229  auto endc =
1230  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
1231  GetIntegratorDataContainerFunc().InitializeEDC(pde.GetUpdateFlags(),
1232  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
1233  this->GetParamData(), this->GetDomainData());
1234  auto& edc = GetIntegratorDataContainerFunc().GetElementDataContainer();
1235 
1236  for (; element[0] != endc[0]; element[0]++)
1237  {
1238  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
1239  {
1240  if (element[dh] == endc[dh])
1241  {
1242  throw DOpEException(
1243  "Elementnumbers in DoFHandlers are not matching!",
1244  "Integrator::ComputeDomainScalar");
1245  }
1246  }
1247 
1248  edc.ReInit();
1249  ret += pde.ElementFunctional(edc);
1250 
1251  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
1252  {
1253  element[dh]++;
1254  }
1255  }
1256  return ret;
1257  }
1258  }
1259  /*******************************************************************************************/
1260 
1261  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1262  int dim>
1263  template<typename PROBLEM>
1264  SCALAR
1265  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputePointScalar(
1266  PROBLEM& pde)
1267  {
1268 
1269  {
1270  SCALAR ret = 0.;
1271  ret += pde.PointFunctional(this->GetParamData(),
1272  this->GetDomainData());
1273 
1274  return ret;
1275  }
1276  }
1277  /*******************************************************************************************/
1278 
1279  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1280  int dim>
1281  template<typename PROBLEM>
1282  SCALAR
1283  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeBoundaryScalar(
1284  PROBLEM& pde
1285  )
1286  {
1287  {
1288  SCALAR ret = 0.;
1289  // Begin integration
1290  const auto& dof_handler =
1291  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
1292  auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
1293  auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
1294 
1295  GetIntegratorDataContainerFunc().InitializeFDC(pde.GetFaceUpdateFlags(),
1296  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
1297  element,
1298  this->GetParamData(),
1299  this->GetDomainData());
1300  auto & fdc = GetIntegratorDataContainerFunc().GetFaceDataContainer();
1301 
1302  std::vector<unsigned int> boundary_functional_colors = pde.GetBoundaryFunctionalColors();
1303  bool need_boundary_integrals = (boundary_functional_colors.size() > 0);
1304  if(!need_boundary_integrals)
1305  {
1306  throw DOpEException("No boundary colors given!","Integrator::ComputeBoundaryScalar");
1307  }
1308 
1309  for (;element[0]!=endc[0]; element[0]++)
1310  {
1311  for(unsigned int dh=1; dh<dof_handler.size(); dh++)
1312  {
1313  if( element[dh] == endc[dh])
1314  {
1315  throw DOpEException("Elementnumbers in DoFHandlers are not matching!","Integrator::ComputeBoundaryScalar");
1316  }
1317  }
1318 
1319  if(need_boundary_integrals)
1320  {
1321  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1322  {
1323 #if DEAL_II_VERSION_GTE(8,3,0)
1324  if (element[0]->face(face)->at_boundary()
1325  &&
1326  (find(boundary_functional_colors.begin(),boundary_functional_colors.end(),
1327  element[0]->face(face)->boundary_id()) != boundary_functional_colors.end()))
1328 #else
1329  if (element[0]->face(face)->at_boundary()
1330  &&
1331  (find(boundary_functional_colors.begin(),boundary_functional_colors.end(),
1332  element[0]->face(face)->boundary_indicator()) != boundary_functional_colors.end()))
1333 #endif
1334  {
1335  fdc.ReInit(face);
1336  ret += pde.BoundaryFunctional(fdc);
1337  }
1338  }
1339  }
1340  for(unsigned int dh=1; dh<dof_handler.size(); dh++)
1341  {
1342  element[dh]++;
1343  }
1344  }
1345 
1346  return ret;
1347 
1348  }
1349  }
1350  /*******************************************************************************************/
1351 
1352  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1353  int dim>
1354  template<typename PROBLEM>
1355  SCALAR
1356  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeFaceScalar(
1357  PROBLEM& pde
1358  )
1359  {
1360 
1361  {
1362  SCALAR ret = 0.;
1363  // Begin integration
1364  const auto& dof_handler =
1365  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
1366  auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
1367  auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
1368 
1369  GetIntegratorDataContainerFunc().InitializeFDC(pde.GetFaceUpdateFlags(),
1370  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
1371  element,
1372  this->GetParamData(),
1373  this->GetDomainData());
1374  auto & fdc = GetIntegratorDataContainerFunc().GetFaceDataContainer();
1375 
1376  bool need_faces = pde.HasFaces();
1377  if(!need_faces)
1378  {
1379  throw DOpEException("No faces required!","Integrator::ComputeFaceScalar");
1380  }
1381 
1382  for (;element[0]!=endc[0]; element[0]++)
1383  {
1384  for(unsigned int dh=1; dh<dof_handler.size(); dh++)
1385  {
1386  if( element[dh] == endc[dh])
1387  {
1388  throw DOpEException("Elementnumbers in DoFHandlers are not matching!","Integrator::ComputeFaceScalar");
1389  }
1390  }
1391 
1392  if(need_faces)
1393  {
1394  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1395  {
1396  if (element[0]->neighbor_index(face) != -1)
1397  {
1398  fdc.ReInit(face);
1399  ret +=pde.FaceFunctional(fdc);
1400  }
1401  }
1402  }
1403  for(unsigned int dh=1; dh<dof_handler.size(); dh++)
1404  {
1405  element[dh]++;
1406  }
1407  }
1408  return ret;
1409  }
1410  }
1411  /*******************************************************************************************/
1412 
1413  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1414  int dim>
1415  template<typename PROBLEM>
1416  SCALAR
1417  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeAlgebraicScalar(
1418  PROBLEM& pde)
1419  {
1420 
1421  {
1422  SCALAR ret = 0.;
1423  ret = pde.AlgebraicFunctional(this->GetParamData(),
1424  this->GetDomainData());
1425  return ret;
1426  }
1427  }
1428 
1429  /*******************************************************************************************/
1430 
1431  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1432  int dim>
1433  template<typename PROBLEM>
1434  void
1435  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearAlgebraicResidual(
1436  PROBLEM& pde, VECTOR &residual)
1437  {
1438  residual = 0.;
1439  pde.AlgebraicResidual(residual, this->GetParamData(),
1440  this->GetDomainData());
1441  //Check if some preset righthandside exists.
1442  AddPresetRightHandSide(-1.,residual);
1443  }
1444 
1445  /*******************************************************************************************/
1446 
1447  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1448  int dim>
1449  template<typename PROBLEM>
1450  void
1451  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeLocalControlConstraints(
1452  PROBLEM& pde, VECTOR &constraints)
1453  {
1454  constraints = 0.;
1455  pde.ComputeLocalControlConstraints(constraints, this->GetParamData(),
1456  this->GetDomainData());
1457  }
1458 
1459  /*******************************************************************************************/
1460 
1461  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1462  int dim>
1463  template<typename PROBLEM>
1464  void
1465  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyTransposedInitialBoundaryValues(
1466  PROBLEM& /*pde*/, VECTOR &/*u*/)
1467  {
1468  // Is not required here ...
1469  }
1470 
1471  /*******************************************************************************************/
1472 
1473  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1474  int dim>
1475  template<typename PROBLEM>
1476  void
1477  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyInitialBoundaryValues(
1478  PROBLEM& pde, VECTOR &u)
1479  {
1480  //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
1481  // e.g. sometimes we need zero sometimes we need other values.
1482 
1483  //Never Condense Nodes Here ! Or All will fail if the state is not initialized with zero!
1484  //pde.GetDoFConstraints().condense(u);
1485  std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
1486  for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
1487  {
1488  unsigned int color = dirichlet_colors[i];
1489  std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
1490  std::map<unsigned int, SCALAR> boundary_values;
1491 
1492  InterpolateBoundaryValues( pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapping(),
1493  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
1494  color,
1495  pde.GetDirichletValues(color, this->GetParamData(),
1496  this->GetDomainData()), boundary_values, comp_mask);
1497 
1498  for (typename std::map<unsigned int, SCALAR>::const_iterator p =
1499  boundary_values.begin(); p != boundary_values.end(); p++)
1500  {
1501  u(p->first) = p->second;
1502  }
1503  }
1504  }
1505 
1506  /*******************************************************************************************/
1507 
1508  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1509  int dim>
1510  template<typename PROBLEM>
1511  void
1512  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyNewtonBoundaryValues(
1513  PROBLEM& pde, VECTOR &u)
1514  {
1515  //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
1516  // e.g. sometimes we need zero sometimes we need other values.
1517 
1518  pde.GetDoFConstraints().condense(u);
1519  std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
1520  for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
1521  {
1522  unsigned int color = dirichlet_colors[i];
1523  std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
1524  std::map<unsigned int, SCALAR> boundary_values;
1525 
1526  InterpolateBoundaryValues( pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapping(),
1527  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
1528  color, dealii::ZeroFunction<dim>(comp_mask.size()),
1529  boundary_values, comp_mask);
1530 
1531  for (typename std::map<unsigned int, SCALAR>::const_iterator p =
1532  boundary_values.begin(); p != boundary_values.end(); p++)
1533  {
1534  u(p->first) = p->second;
1535  }
1536  }
1537  }
1538  /*******************************************************************************************/
1539 
1540  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1541  int dim>
1542  template<typename PROBLEM, typename MATRIX>
1543  void
1544  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyNewtonBoundaryValues(
1545  PROBLEM& pde, MATRIX& matrix, VECTOR &rhs, VECTOR &sol)
1546  {
1547  //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
1548  // e.g. sometimes we need zero sometimes we need other values.
1549  pde.GetDoFConstraints().condense(rhs);
1550  pde.GetDoFConstraints().condense(matrix);
1551  std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
1552  for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
1553  {
1554  unsigned int color = dirichlet_colors[i];
1555  std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
1556  std::map<unsigned int, SCALAR> boundary_values;
1557 
1558  InterpolateBoundaryValues( pde.GetBaseProblem().GetSpaceTimeHandler()->GetMapping(),
1559  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
1560  color, dealii::ZeroFunction<dim>(comp_mask.size()),
1561  boundary_values, comp_mask);
1562 
1563  dealii::MatrixTools::apply_boundary_values(boundary_values, matrix,
1564  sol, rhs);
1565  }
1566  }
1567 
1568  /*******************************************************************************************/
1569 
1570  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1571  int dim>
1572  void
1573  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddDomainData(
1574  std::string name, const VECTOR* new_data)
1575  {
1576  if (domain_data_.find(name) != domain_data_.end())
1577  {
1578  throw DOpEException(
1579  "Adding multiple Data with name " + name + " is prohibited!",
1580  "Integrator::AddDomainData");
1581  }
1582  domain_data_.insert(
1583  std::pair<std::string, const VECTOR*>(name, new_data));
1584  }
1585 
1586  /*******************************************************************************************/
1587 
1588  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1589  int dim>
1590  void
1591  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::DeleteDomainData(
1592  std::string name)
1593  {
1594  typename std::map<std::string, const VECTOR *>::iterator it =
1595  domain_data_.find(name);
1596  if (it == domain_data_.end())
1597  {
1598  throw DOpEException(
1599  "Deleting Data " + name + " is impossible! Data not found",
1600  "Integrator::DeleteDomainData");
1601  }
1602  domain_data_.erase(it);
1603  }
1604 
1605  /*******************************************************************************************/
1606 
1607  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1608  int dim>
1609  const std::map<std::string, const VECTOR*>&
1610  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetDomainData() const
1611  {
1612  return domain_data_;
1613  }
1614 
1615  /*******************************************************************************************/
1616 
1617  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1618  int dim>
1619  void
1620  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddParamData(
1621  std::string name, const dealii::Vector<SCALAR>* new_data)
1622  {
1623  if (param_data_.find(name) != param_data_.end())
1624  {
1625  throw DOpEException(
1626  "Adding multiple Data with name " + name + " is prohibited!",
1627  "Integrator::AddParamData");
1628  }
1629  param_data_.insert(
1630  std::pair<std::string, const dealii::Vector<SCALAR>*>(name,
1631  new_data));
1632  }
1633 
1634  /*******************************************************************************************/
1635 
1636  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1637  int dim>
1638  void
1639  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::DeleteParamData(
1640  std::string name)
1641  {
1642  typename std::map<std::string, const dealii::Vector<SCALAR>*>::iterator it =
1643  param_data_.find(name);
1644  if (it == param_data_.end())
1645  {
1646  throw DOpEException(
1647  "Deleting Data " + name + " is impossible! Data not found",
1648  "Integrator::DeleteParamData");
1649  }
1650  param_data_.erase(it);
1651  }
1652 
1653  /*******************************************************************************************/
1654 
1655  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1656  int dim>
1657  const std::map<std::string, const dealii::Vector<SCALAR>*>&
1658  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetParamData() const
1659  {
1660  return param_data_;
1661  }
1662 
1663  /*******************************************************************************************/
1664 
1665  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1666  int dim>
1667  template<typename PROBLEM, class STH, class EDC, class FDC>
1668  void
1669  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeRefinementIndicators(
1670  PROBLEM& pde,
1671  DWRDataContainer<STH, INTEGRATORDATACONT, EDC, FDC, VECTOR>& dwrc)
1672  {
1673  unsigned int n_error_comps = dwrc.GetNErrorComps();
1674  //for primal and dual part of the error
1675  std::vector<double> element_sum(n_error_comps, 0);
1676  element_sum.resize(n_error_comps, 0);
1677 
1678  // Begin integration
1679  const auto& dof_handler =
1680  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
1681  auto element =
1682  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
1683  auto endc =
1684  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
1685 
1686  const auto& dof_handler_weight = dwrc.GetWeightSTH().GetDoFHandler();
1687  auto element_weight = dwrc.GetWeightSTH().GetDoFHandlerBeginActive();
1688  auto endc_high = dwrc.GetWeightSTH().GetDoFHandlerEnd();
1689 
1690  // Generate the data containers. Notice that we use the quadrature
1691  //formula from the higher order idc!.
1692  GetIntegratorDataContainer().InitializeEDC(
1693  dwrc.GetWeightIDC().GetQuad(), pde.GetUpdateFlags(),
1694  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
1695  this->GetParamData(), this->GetDomainData());
1696  auto& edc = GetIntegratorDataContainer().GetElementDataContainer();
1697 
1698  dwrc.GetWeightIDC().InitializeEDC(pde.GetUpdateFlags(),
1699  dwrc.GetWeightSTH(), element_weight, this->GetParamData(),
1700  dwrc.GetWeightData());
1701  auto& edc_weight = dwrc.GetElementWeight();
1702 
1703  // we want to integrate the face-terms only once, so
1704  // we store the values on each face in this map
1705  // and distribute it at the end to the adjacent elements.
1706  typename std::map<typename dealii::Triangulation<dim>::face_iterator,std::vector<double> >
1707  face_integrals;
1708  // initialize the map with a big value to make sure
1709  // that we take notice if we forget to add a face
1710  // during the error estimation process
1711  auto element_it = element[0];
1712  std::vector<double> face_init(n_error_comps,-1e20);
1713  for (; element_it != endc[0]; element_it++)
1714  {
1715  for (unsigned int face_no=0;face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
1716  {
1717  face_integrals[element_it->face(face_no)] = face_init;
1718  }
1719  }
1720 
1721 // bool need_faces = pde.HasFaces();
1722 // bool need_interfaces = pde.HasInterfaces();
1723 // std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
1724 // bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
1725 
1726  GetIntegratorDataContainer().InitializeFDC(dwrc.GetWeightIDC().GetFaceQuad(),
1727  pde.GetFaceUpdateFlags(),
1728  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
1729  element,
1730  this->GetParamData(),
1731  this->GetDomainData(),
1732  true);
1733  auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();
1734 
1735  dwrc.GetWeightIDC().InitializeFDC(pde.GetFaceUpdateFlags(),
1736  dwrc.GetWeightSTH(),
1737  element_weight,
1738  this->GetParamData(),
1739  dwrc.GetWeightData(),
1740  true);
1741 
1742  for (unsigned int element_index = 0; element[0] != endc[0];
1743  element[0]++, element_index++)
1744  {
1745  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
1746  {
1747  if (element[dh] == endc[dh])
1748  {
1749  throw DOpEException(
1750  "Elementnumbers in DoFHandlers are not matching!",
1751  "Integrator::ComputeRefinementIndicators");
1752  }
1753  }
1754  for (unsigned int dh = 0; dh < dof_handler_weight.size(); dh++)
1755  {
1756  if (element_weight[dh] == endc_high[dh])
1757  {
1758  throw DOpEException(
1759  "Elementnumbers in DoFHandlers are not matching!",
1760  "Integrator::ComputeRefinementIndicators");
1761  }
1762  }
1763  element_sum.clear();
1764  element_sum.resize(n_error_comps, 0);
1765 
1766  edc.ReInit();
1767  edc_weight.ReInit();
1768 
1769  //first the element-residual
1770  pde.ElementErrorContribution(edc, dwrc, element_sum, 1.);
1771  for(unsigned int l =0; l < n_error_comps; l ++)
1772  {
1773  dwrc.GetErrorIndicators(l)(element_index) = element_sum[l];
1774  }
1775  //dwrc.GetPrimalErrorIndicators()(element_index) = element_sum[0];
1776  //dwrc.GetDualErrorIndicators()(element_index) = element_sum[1];
1777  element_sum.clear();
1778  element_sum.resize(n_error_comps, 0);
1779 
1780 
1781  //Now to the face terms. We compute them only once for each face and distribute the
1782  //afterwards. We choose always to work from the coarser element, if both neigbors of the
1783  //face are on the same level, we pick the one with the lower index
1784  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1785  {
1786  auto face_it = element[0]->face(face);
1787 
1788  //check if the face lies at a boundary
1789  if(face_it->at_boundary())
1790  {
1791  fdc.ReInit(face);
1792  dwrc.GetFaceWeight().ReInit(face);
1793  pde.BoundaryErrorContribution(fdc, dwrc, element_sum, 1.);
1794 
1795  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
1796  ExcInternalError());
1797  Assert (face_integrals[element[0]->face(face)] == face_init,
1798  ExcInternalError());
1799 
1800  face_integrals[element[0]->face(face)] = element_sum;
1801  element_sum.clear();
1802  element_sum.resize(n_error_comps,0.);
1803  }
1804  else
1805  {
1806  //There exist now 3 different scenarios, given the actual element and face:
1807  // The neighbour behind this face is [ more | as much | less] refined
1808  // than/as the actual element. We have to distinguish here only between the case 1
1809  // and the other two, because these will be distinguished in in the FaceDataContainer.
1810  if (element[0]->neighbor(face)->has_children())
1811  {
1812  //first: neighbour is finer
1813  std::vector<double> sum(n_error_comps,0.);
1814  for (unsigned int subface_no=0;
1815  subface_no < element[0]->face(face)->n_children();
1816  ++subface_no)
1817  {
1818  //TODO Now here we have to initialise the subface_values on the
1819  // actual element and then the facevalues of the neighbours
1820  fdc.ReInit(face, subface_no);
1821  fdc.ReInitNbr();
1822  dwrc.GetFaceWeight().ReInit(face, subface_no);
1823 
1824  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
1825  for(unsigned int l =0; l < n_error_comps; l ++)
1826  {
1827  sum[l] += element_sum[l];
1828  }
1829  face_integrals[element[0]->neighbor_child_on_subface(face, subface_no)
1830  ->face(element[0]->neighbor_of_neighbor(face))] = element_sum;
1831  element_sum.clear();
1832  element_sum.resize(n_error_comps,0);
1833  }
1834 
1835  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
1836  ExcInternalError());
1837  Assert (face_integrals[element[0]->face(face)] == face_init,
1838  ExcInternalError());
1839 
1840  face_integrals[element[0]->face(face)] = sum;
1841  element_sum.clear();
1842  element_sum.resize(n_error_comps,0);
1843  }
1844  else
1845  {
1846  // either neighbor is as fine as this element or
1847  // it is coarser
1848  Assert(element[0]->neighbor(face)->level() <= element[0]->level(),ExcInternalError());
1849  //now we work always from the coarser element. if both elements
1850  //are on the same level, we pick the one with the lower index
1851  if(element[0]->level() == element[0]->neighbor(face)->level()
1852  && element[0]->index() < element[0]->neighbor(face)->index())
1853  {
1854  fdc.ReInit(face);
1855  fdc.ReInitNbr();
1856  dwrc.GetFaceWeight().ReInit(face);
1857 
1858  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
1859  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
1860  ExcInternalError());
1861  Assert (face_integrals[element[0]->face(face)] == face_init,
1862  ExcInternalError());
1863 
1864  face_integrals[element[0]->face(face)] = element_sum;
1865  element_sum.clear();
1866  element_sum.resize(n_error_comps,0);
1867  }
1868  }
1869  }
1870  } //endfor faces
1871  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
1872  {
1873  element[dh]++;
1874  }
1875  for (unsigned int dh = 0; dh < dof_handler_weight.size(); dh++)
1876  {
1877  element_weight[dh]++;
1878  }
1879  } //endfor element
1880  //now we have to incorporate the face and boundary_values
1881  //into
1882  unsigned int present_element = 0;
1883  element =
1884  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
1885  for (;
1886  element[0] !=endc[0]; element[0]++, ++present_element)
1887  {
1888  for (unsigned int face_no = 0;
1889  face_no < GeometryInfo<dim>::faces_per_cell; ++face_no)
1890  {
1891  Assert(
1892  face_integrals.find(element[0]->face(face_no)) != face_integrals.end(),
1893  ExcInternalError());
1894  if(element[0]->face(face_no)->at_boundary())
1895  {
1896  for(unsigned int l =0; l < n_error_comps; l ++)
1897  {
1898  dwrc.GetErrorIndicators(l)(present_element) +=
1899  face_integrals[element[0]->face(face_no)][l];
1900  }
1901 // dwrc.GetPrimalErrorIndicators()(present_element) +=
1902 // face_integrals[element[0]->face(face_no)][0];
1903 // dwrc.GetDualErrorIndicators()(present_element) +=
1904 // face_integrals[element[0]->face(face_no)][1];
1905  }
1906  else
1907  {
1908  for(unsigned int l =0; l < n_error_comps; l ++)
1909  {
1910  dwrc.GetErrorIndicators(l)(present_element) +=
1911  0.5*face_integrals[element[0]->face(face_no)][l];
1912  }
1913 // dwrc.GetPrimalErrorIndicators()(present_element) +=
1914 // 0.5 * face_integrals[element[0]->face(face_no)][0];
1915 // dwrc.GetDualErrorIndicators()(present_element) +=
1916 // 0.5 * face_integrals[element[0]->face(face_no)][1];
1917  }
1918 
1919  }
1920  }
1921  }
1922  /*******************************************************************************************/
1923 
1924  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1925  int dim>
1926  template<typename PROBLEM>
1927  void
1928  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeRefinementIndicators(
1929  PROBLEM& pde,
1930  ResidualErrorContainer<VECTOR>& dwrc)
1931  {
1932  //for primal and dual part of the error
1933  std::vector<double> element_sum(2, 0);
1934  element_sum.resize(2, 0);
1935  // Begin integration
1936  const auto& dof_handler =
1937  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
1938  auto element =
1939  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
1940  auto endc =
1941  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
1942 
1943  {
1944  //Add Weights
1945  auto wd = dwrc.GetWeightData().begin();
1946  auto wend = dwrc.GetWeightData().end();
1947  for (; wd != wend; wd++)
1948  {
1949  AddDomainData(wd->first, wd->second);
1950  }
1951  }
1952 
1953  // Generate the data containers. Notice that we use the quadrature
1954  // formula from the higher order idc!.
1955  GetIntegratorDataContainer().InitializeEDC(pde.GetUpdateFlags(),
1956  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element,
1957  this->GetParamData(), this->GetDomainData());
1958  auto& edc = GetIntegratorDataContainer().GetElementDataContainer();
1959 
1960  //we want to integrate the face-terms only once
1961  typename std::map<typename dealii::Triangulation<dim>::face_iterator,std::vector<double> >
1962  face_integrals;
1963  //initialize the map
1964  auto element_it = element[0];
1965  std::vector<double> face_init(2,-1e20);
1966  for (; element_it != endc[0]; element_it++)
1967  {
1968  for (unsigned int face_no=0;face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
1969  {
1970  face_integrals[element_it->face(face_no)] = face_init;
1971  }
1972  }
1973 
1974 // bool need_faces = pde.HasFaces();
1975 // bool need_interfaces = pde.HasInterfaces();
1976 // std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
1977 // bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
1978 
1979  GetIntegratorDataContainer().InitializeFDC(
1980  pde.GetFaceUpdateFlags(),
1981  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
1982  element,
1983  this->GetParamData(),
1984  this->GetDomainData(),
1985  true);
1986  auto & fdc = GetIntegratorDataContainer().GetFaceDataContainer();
1987 
1988  for (unsigned int element_index = 0; element[0] != endc[0];
1989  element[0]++, element_index++)
1990  {
1991  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
1992  {
1993  if (element[dh] == endc[dh])
1994  {
1995  throw DOpEException(
1996  "Elementnumbers in DoFHandlers are not matching!",
1997  "Integrator::ComputeRefinementIndicators");
1998  }
1999  }
2000 
2001  element_sum.clear();
2002  element_sum.resize(2, 0);
2003 
2004  edc.ReInit();
2005  dwrc.InitElement(element[0]->diameter());
2006  //first the element-residual
2007  pde.ElementErrorContribution(edc, dwrc, element_sum, 1.);
2008  dwrc.GetPrimalErrorIndicators()(element_index) = element_sum[0];
2009  dwrc.GetDualErrorIndicators()(element_index) = element_sum[1];
2010  element_sum.clear();
2011  element_sum.resize(2, 0);
2012  //Now to the face terms. We compute them only once for each face and distribute the
2013  //afterwards. We choose always to work from the coarser element, if both neigbors of the
2014  //face are on the same level, we pick the one with the lower index
2015 
2016  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
2017  {
2018  auto face_it = element[0]->face(face);
2019 
2020  //check if the face lies at a boundary
2021  if(face_it->at_boundary())
2022  {
2023  fdc.ReInit(face);
2024 #if deal_II_dimension > 1
2025  dwrc.InitFace(element[0]->face(face)->diameter());
2026 #else
2027  dwrc.InitFace(element[0]->diameter());
2028 #endif
2029  pde.BoundaryErrorContribution(fdc, dwrc, element_sum, 1.);
2030 
2031  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
2032  ExcInternalError());
2033  Assert (face_integrals[element[0]->face(face)] == face_init,
2034  ExcInternalError());
2035  face_integrals[element[0]->face(face)] = element_sum;
2036  element_sum.clear();
2037  element_sum.resize(2,0.);
2038  }
2039  else
2040  {
2041  //There exist now 3 different scenarios, given the actual element and face:
2042  // The neighbour behind this face is [ more | as much | less] refined
2043  // than/as the actual element. We have to distinguish here only between the case 1
2044  // and the other two, because these will be distinguished in in the FaceDataContainer.
2045  if (element[0]->neighbor(face)->has_children())
2046  {
2047  //first: neighbour is finer
2048  std::vector<double> sum(2,0.);
2049  for (unsigned int subface_no=0;
2050  subface_no < element[0]->face(face)->n_children();
2051  ++subface_no)
2052  {
2053  //TODO Now here we have to initialise the subface_values on the
2054  // actual element and then the facevalues of the neighbours
2055  fdc.ReInit(face, subface_no);
2056  fdc.ReInitNbr();
2057 #if deal_II_dimension > 1
2058  dwrc.InitFace(element[0]->face(face)->diameter());
2059 #else
2060  dwrc.InitFace(element[0]->diameter());
2061 #endif
2062 
2063  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
2064  sum[0]= element_sum[0];
2065  sum[1]= element_sum[1];
2066  element_sum.clear();
2067  element_sum.resize(2,0);
2068  face_integrals[element[0]->neighbor_child_on_subface(face, subface_no)
2069  ->face(element[0]->neighbor_of_neighbor(face))] = element_sum;
2070  element_sum.clear();
2071  element_sum.resize(2,0.);
2072  }
2073 
2074  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
2075  ExcInternalError());
2076  Assert (face_integrals[element[0]->face(face)] == face_init,
2077  ExcInternalError());
2078 
2079  face_integrals[element[0]->face(face)] = sum;
2080 
2081  }
2082  else
2083  {
2084  // either neighbor is as fine as this element or
2085  // it is coarser
2086  Assert(element[0]->neighbor(face)->level() <= element[0]->level(),ExcInternalError());
2087  //now we work always from the coarser element. if both elements
2088  //are on the same level, we pick the one with the lower index
2089  if(element[0]->level() == element[0]->neighbor(face)->level()
2090  && element[0]->index() < element[0]->neighbor(face)->index())
2091  {
2092  fdc.ReInit(face);
2093  fdc.ReInitNbr();
2094 #if deal_II_dimension > 1
2095  dwrc.InitFace(element[0]->face(face)->diameter());
2096 #else
2097  dwrc.InitFace(element[0]->diameter());
2098 #endif
2099 
2100  pde.FaceErrorContribution(fdc, dwrc, element_sum, 1.);
2101  Assert (face_integrals.find (element[0]->face(face)) != face_integrals.end(),
2102  ExcInternalError());
2103  Assert (face_integrals[element[0]->face(face)] == face_init,
2104  ExcInternalError());
2105 
2106  face_integrals[element[0]->face(face)] = element_sum;
2107  element_sum.clear();
2108  element_sum.resize(2,0);
2109  }
2110  }
2111  }
2112  } //endfor faces
2113 // }//end else
2114 
2115  for (unsigned int dh = 1; dh < dof_handler.size(); dh++)
2116  {
2117  element[dh]++;
2118  }
2119  } //endfor element
2120  //now we have to incorporate the face and boundary_values
2121  //into
2122  unsigned int present_element = 0;
2123  element =
2124  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
2125  for (;
2126  element[0] !=endc[0]; element[0]++, ++present_element)
2127  for (unsigned int face_no = 0;
2128  face_no < GeometryInfo<dim>::faces_per_cell; ++face_no)
2129  {
2130  Assert(
2131  face_integrals.find(element[0]->face(face_no)) != face_integrals.end(),
2132  ExcInternalError());
2133  dwrc.GetPrimalErrorIndicators()(present_element) +=
2134  0.5 * face_integrals[element[0]->face(face_no)][0];
2135  dwrc.GetDualErrorIndicators()(present_element) +=
2136  0.5 * face_integrals[element[0]->face(face_no)][1];
2137  }
2138  {
2139  //Remove Weights
2140  auto wd = dwrc.GetWeightData().begin();
2141  auto wend = dwrc.GetWeightData().end();
2142  for (; wd != wend; wd++)
2143  {
2144  DeleteDomainData(wd->first);
2145  }
2146  }
2147  }
2148 
2149  /*******************************************************************************************/
2150 
2151  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
2152  int dim>
2153  INTEGRATORDATACONT&
2154  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetIntegratorDataContainer() const
2155  {
2156  return idc1_;
2157  }
2158 
2159  /*******************************************************************************************/
2160 
2161  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
2162  int dim>
2163  INTEGRATORDATACONT&
2164  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetIntegratorDataContainerFunc() const
2165  {
2166  return idc2_;
2167  }
2168 
2169  /*******************************************************************************************/
2170 
2171  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
2172  int dim>
2173  template<template<int, int> class DH>
2174  void
2175  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::InterpolateBoundaryValues(
2176  const DOpEWrapper::Mapping<dim, DH>& mapping,
2177  const DOpEWrapper::DoFHandler<dim, DH>* dof_handler,
2178  const unsigned int color, const dealii::Function<dim>& function,
2179  std::map<unsigned int, SCALAR>& boundary_values,
2180  const std::vector<bool>& comp_mask) const
2181  {
2182  //TODO: mapping[0] is a workaround, as deal does not support interpolate
2183  // boundary_values with a mapping collection at this point.
2184  dealii::VectorTools::interpolate_boundary_values(mapping[0],
2185  dof_handler->GetDEALDoFHandler(), color, function,
2186  boundary_values, comp_mask);
2187  }
2188  /*******************************************************************************************/
2189 
2190  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
2191  int dim>
2192  void
2193  Integrator<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddPresetRightHandSide(double s,
2194  VECTOR &residual) const
2195  {
2196  typename std::map<std::string, const VECTOR *>::const_iterator it =
2197  domain_data_.find("fixed_rhs");
2198  if (it != domain_data_.end())
2199  {
2200  assert(residual.size() == it->second->size());
2201  residual.add(s,*(it->second));
2202  }
2203  }
2204 
2205 
2206 }
2207 #endif
2208 
DOpE::Integrator::ComputeDomainScalar
SCALAR ComputeDomainScalar(PROBLEM &pde)
Definition: integrator.h:1219
DOpEWrapper::Mapping
Definition: mapping_wrapper.h:48
DOpE::DWRDataContainerBase::GetWeightData
const std::map< std::string, const VECTOR * > & GetWeightData() const
Definition: dwrdatacontainer.h:398
DOpE::Integrator::ComputeBoundaryScalar
SCALAR ComputeBoundaryScalar(PROBLEM &pde)
Definition: integrator.h:1283
DOpE::Integrator::AddParamData
void AddParamData(std::string name, const dealii::Vector< SCALAR > *new_data)
Definition: integrator.h:1620
DOpE::Integrator::ComputeNonlinearLhs
void ComputeNonlinearLhs(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integrator.h:724
DOpE::DWRDataContainer::GetElementWeight
virtual EDC & GetElementWeight() const =0
DOpE::Integrator::~Integrator
~Integrator()
Definition: integrator.h:525
DOpE::Integrator::ComputePointScalar
SCALAR ComputePointScalar(PROBLEM &pde)
Definition: integrator.h:1265
DOpE::Integrator
Definition: integrator.h:57
DOpE::Integrator::ApplyInitialBoundaryValues
void ApplyInitialBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integrator.h:1477
DOpE::DWRDataContainerBase::GetPrimalErrorIndicators
Vector< double > & GetPrimalErrorIndicators()
Definition: dwrdatacontainer.h:233
DOpE::Integrator::ComputeLocalControlConstraints
void ComputeLocalControlConstraints(PROBLEM &pde, VECTOR &constraints)
Definition: integrator.h:1451
DOpE::Integrator::ComputeMatrix
void ComputeMatrix(PROBLEM &pde, MATRIX &matrix)
Definition: integrator.h:1023
DOpEWrapper::DoFHandler::GetDEALDoFHandler
const DOFHANDLER< dim, dim > & GetDEALDoFHandler() const
Definition: dofhandler_wrapper.h:69
DOpE::Integrator::ComputeRefinementIndicators
void ComputeRefinementIndicators(PROBLEM &pde, DWRDataContainer< STH, INTEGRATORDATACONT, EDC, FDC, VECTOR > &dwrc)
Definition: integrator.h:1669
DOpE::DWRDataContainerBase::GetErrorIndicators
const Vector< double > & GetErrorIndicators() const
Definition: dwrdatacontainer.h:213
DOpE::Integrator::GetParamData
const std::map< std::string, const dealii::Vector< SCALAR > * > & GetParamData() const
Definition: integrator.h:1658
DOpE::Integrator::GetIntegratorDataContainerFunc
INTEGRATORDATACONT & GetIntegratorDataContainerFunc() const
Definition: integrator.h:2164
DOpE::DWRDataContainerBase::GetDualErrorIndicators
Vector< double > & GetDualErrorIndicators()
Definition: dwrdatacontainer.h:251
DOpE::DWRDataContainer::GetWeightSTH
virtual STH & GetWeightSTH()=0
DOpE::Integrator::ComputeFaceScalar
SCALAR ComputeFaceScalar(PROBLEM &pde)
Definition: integrator.h:1356
DOpE::Integrator::AddPresetRightHandSide
void AddPresetRightHandSide(double s, VECTOR &residual) const
Definition: integrator.h:2193
DOpE::Integrator::DeleteParamData
void DeleteParamData(std::string name)
Definition: integrator.h:1639
DOpE::Integrator::DeleteDomainData
void DeleteDomainData(std::string name)
Definition: integrator.h:1591
DOpE::Integrator::Integrator
Integrator(INTEGRATORDATACONT &idc1)
Definition: integrator.h:507
DOpE::ResidualErrorContainer
Definition: residualestimator.h:38
DOpE::Integrator::ComputeAlgebraicScalar
SCALAR ComputeAlgebraicScalar(PROBLEM &pde)
Definition: integrator.h:1417
DOpE::Integrator::AddDomainData
void AddDomainData(std::string name, const VECTOR *new_data)
Definition: integrator.h:1573
DOpE::Integrator::ApplyTransposedInitialBoundaryValues
void ApplyTransposedInitialBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integrator.h:1465
DOpE::DWRDataContainer::GetFaceWeight
virtual FDC & GetFaceWeight() const =0
DOpE::Integrator::GetIntegratorDataContainer
INTEGRATORDATACONT & GetIntegratorDataContainer() const
Definition: integrator.h:2154
DOpE::DWRDataContainerBase::GetNErrorComps
unsigned int GetNErrorComps() const
Definition: dwrdatacontainer.h:317
DOpE::ResidualErrorContainer::InitFace
virtual void InitFace(double)=0
DOpE::Integrator::ReInit
void ReInit()
Definition: integrator.h:535
DOpE::DWRDataContainer::GetWeightIDC
virtual IDC & GetWeightIDC()=0
DOpE::DOpEException
Definition: dopeexception.h:35
DOpE::Integrator::ApplyNewtonBoundaryValues
void ApplyNewtonBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integrator.h:1512
DOpE::Integrator::ComputeNonlinearAlgebraicResidual
void ComputeNonlinearAlgebraicResidual(PROBLEM &pde, VECTOR &residual)
Definition: integrator.h:1435
DOpE::DWRDataContainer
Definition: dwrdatacontainer.h:545
DOpE::Integrator::ComputeNonlinearRhs
void ComputeNonlinearRhs(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integrator.h:891
DOpE::ResidualErrorContainer::InitElement
virtual void InitElement(double)=0
DOpE::Integrator::ComputeNonlinearResidual
void ComputeNonlinearResidual(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integrator.h:546
DOpE::Integrator::GetDomainData
const std::map< std::string, const VECTOR * > & GetDomainData() const
Definition: integrator.h:1610
Generated on Mon Jan 25 2016 11:32:19 for DOpE by   doxygen 1.8.6