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