DOpE
integrator_multimesh.h
Go to the documentation of this file.
1 
24 #ifndef IntegratorMultiMesh_H_
25 #define IntegratorMultiMesh_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/grid/grid_tools.h>
33 #include <deal.II/base/function.h>
34 
35 #include <vector>
36 
37 #include "multimesh_elementdatacontainer.h"
38 #include "multimesh_facedatacontainer.h"
39 
40 namespace DOpE
41 {
59  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
60  int dim>
61  class IntegratorMultiMesh
62  {
63  public:
64  IntegratorMultiMesh(INTEGRATORDATACONT& idc);
65 
66  ~IntegratorMultiMesh();
67 
68  void
69  ReInit();
70 
71  template<typename PROBLEM>
72  void
73  ComputeNonlinearResidual(PROBLEM& pde, VECTOR &residual,
74  bool apply_boundary_values = true);
75  template<typename PROBLEM>
76  void
77  ComputeNonlinearLhs(PROBLEM& pde, VECTOR &residual,
78  bool apply_boundary_values = true);
79  template<typename PROBLEM>
80  void
81  ComputeNonlinearRhs(PROBLEM& pde, VECTOR &residual,
82  bool apply_boundary_values = true);
83  template<typename PROBLEM, typename MATRIX>
84  void
85  ComputeMatrix(PROBLEM& pde, MATRIX &matrix);
86 
87  template<typename PROBLEM>
88  void ComputeNonlinearAlgebraicResidual (PROBLEM& pde, VECTOR &residual);
89  template<typename PROBLEM>
90  void ComputeLocalControlConstraints (PROBLEM& pde, VECTOR &constraints);
91 
92  template<typename PROBLEM>
93  SCALAR
94  ComputeDomainScalar(PROBLEM& pde);
95  template<typename PROBLEM>
96  SCALAR
97  ComputePointScalar(PROBLEM& pde);
98  template<typename PROBLEM>
99  SCALAR
100  ComputeBoundaryScalar(PROBLEM& pde);
101  template<typename PROBLEM>
102  SCALAR
103  ComputeFaceScalar(PROBLEM& pde);
104  template<typename PROBLEM>
105  SCALAR
106  ComputeAlgebraicScalar(PROBLEM& pde);
107 
108  template<typename PROBLEM>
109  void
110  ApplyInitialBoundaryValues(PROBLEM& pde, VECTOR &u);
111  template<typename PROBLEM>
112  void
113  ApplyNewtonBoundaryValues(PROBLEM& pde, VECTOR &u);
114  template<typename PROBLEM, typename MATRIX>
115  void
116  ApplyNewtonBoundaryValues(PROBLEM& pde, MATRIX &matrix, VECTOR &rhs,
117  VECTOR &sol);
118 
119  inline void
120  AddDomainData(std::string name, const VECTOR* new_data);
121  inline void
122  DeleteDomainData(std::string name);
123 
124  inline void
125  AddParamData(std::string name, const dealii::Vector<SCALAR>* new_data);
126  inline void
127  DeleteParamData(std::string name);
128 
129  protected:
130  inline const std::map<std::string, const dealii::Vector<SCALAR>*>&
131  GetParamData() const;
132 
133  inline INTEGRATORDATACONT&
134  GetIntegratorDataContainer() const;
135  inline const std::map<std::string, const VECTOR*>&
136  GetDomainData() const;
137  inline void AddPresetRightHandSide(double s, VECTOR &residual) const;
138 
139  private:
140  template<template<int, int> class DH>
141  void
142  InterpolateBoundaryValues(
143  const DOpEWrapper::DoFHandler<dim, DH>* dof_handler,
144  const unsigned int color, const dealii::Function<dim>& function,
145  std::map<unsigned int, SCALAR>& boundary_values,
146  const std::vector<bool>& comp_mask) const;
147 
152  template<typename PROBLEM, template<int, int> class DH>
153  inline void ComputeNonlinearResidual_Recursive(
154  PROBLEM& pde,
155  VECTOR& residual,
156  typename std::vector<typename DH<dim, dim>::cell_iterator>& element_iter,
157  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& tria_element_iter,
158  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
159  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc,
160  Multimesh_FaceDataContainer<DH, VECTOR, dim>& fdc);
161 
166  template<typename PROBLEM, template<int, int> class DH>
167  inline void ComputeNonlinearRhs_Recursive(
168  PROBLEM& pde,
169  VECTOR& residual,
170  typename std::vector<typename DH<dim, dim>::cell_iterator>& element_iter,
171  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& tria_element_iter,
172  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
173  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc,
174  Multimesh_FaceDataContainer<DH, VECTOR, dim>& fdc);
175 
180  template<typename PROBLEM, typename MATRIX, template<int, int> class DH>
181  inline void ComputeMatrix_Recursive(
182  PROBLEM& pde,
183  MATRIX& matrix,
184  typename std::vector<typename DH<dim, dim>::cell_iterator>& element_iter,
185  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& tria_element_iter,
186  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
187  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc,
188  Multimesh_FaceDataContainer<DH, VECTOR, dim>& fdc);
189 
194  template<typename PROBLEM, template<int, int> class DH>
195  inline SCALAR ComputeDomainScalar_Recursive(
196  PROBLEM& pde,
197  typename std::vector<typename DH<dim, dim>::cell_iterator>& element_iter,
198  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& tria_element_iter,
199  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
200  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc);
201 
206  template<typename PROBLEM, template<int, int> class DH>
207  inline SCALAR ComputeBoundaryScalar_Recursive(
208  PROBLEM& pde,
209  typename std::vector<typename DH<dim, dim>::cell_iterator>& element_iter,
210  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& tria_element_iter,
211  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
212  Multimesh_FaceDataContainer<DH, VECTOR, dim>& edc);
213 
214  INTEGRATORDATACONT & idc_;
215 
216  std::map<std::string, const VECTOR*> domain_data_;
217  std::map<std::string, const dealii::Vector<SCALAR>*> param_data_;
218  };
219 
220  /**********************************Implementation*******************************************/
221 
222  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
223  int dim>
224  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::IntegratorMultiMesh(
225  INTEGRATORDATACONT& idc) :
226  idc_(idc)
227  {
228  }
229 
230  /**********************************Implementation*******************************************/
231 
232  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
233  int dim>
234  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::~IntegratorMultiMesh()
235  {
236 
237  }
238 
239  /*******************************************************************************************/
240 
241  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
242  int dim>
243  void
244  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ReInit()
245  {
246 
247  }
248 
249  /*******************************************************************************************/
250 
251  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
252  int dim>
253  template<typename PROBLEM>
254  void
255  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearResidual(
256  PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
257  {
258  residual = 0.;
259 
260  const auto& dof_handler =
261  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
262 
263  assert(dof_handler.size() == 2);
264 
265  const auto tria_element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler().get_tria(),
266  dof_handler[1]->GetDEALDoFHandler().get_tria());
267  const auto element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler(),
268  dof_handler[1]->GetDEALDoFHandler());
269 
270  auto element_iter = element_list.begin();
271  auto tria_element_iter = tria_element_list.begin();
272 
273  std::vector<decltype(tria_element_iter->first)> tria_element(2);
274  tria_element[0] = tria_element_iter->first;
275  tria_element[1] = tria_element_iter->second;
276 
277  std::vector<decltype(element_iter->first)> element(2);
278  element[0] = element_iter->first;
279  element[1] = element_iter->second;
280  int coarse_index = 0; //element[coarse_index] is the coarser of the two.
281  int fine_index = 0; //element[fine_index] is the finer of the two (both indices are 2 = element.size() if they are equally refined.
282 
283  // Generate the data containers.
284  idc_.InitializeMMEDC(pde.GetUpdateFlags(),
285  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element, tria_element,
286  this->GetParamData(), this->GetDomainData());
287  auto& edc = idc_.GetMultimeshElementDataContainer();
288 
289  bool need_interfaces = pde.HasInterfaces();
290  idc_.InitializeMMFDC(pde.GetFaceUpdateFlags(),
291  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
292  element,
293  tria_element,
294  this->GetParamData(),
295  this->GetDomainData(),
296  need_interfaces);
297  auto& fdc = idc_.GetMultimeshFaceDataContainer();
298 
299  for(; element_iter != element_list.end(); element_iter++)
300  {
301  element[0] = element_iter->first;
302  element[1] = element_iter->second;
303  tria_element[0] = tria_element_iter->first;
304  tria_element[1] = tria_element_iter->second;
305  FullMatrix<SCALAR> prolong_matrix;
306 
307  if(element[0]->has_children())
308  {
309  prolong_matrix = IdentityMatrix(element[1]->get_fe().dofs_per_cell);
310  coarse_index =1;
311  fine_index = 0;
312  }
313  else
314  {
315  if(element[1]->has_children())
316  {
317  prolong_matrix = IdentityMatrix(element[0]->get_fe().dofs_per_cell);
318  coarse_index =0;
319  fine_index = 1;
320  }
321  else
322  {
323  assert(element.size() ==2);
324  coarse_index = fine_index = 2;
325  }
326  }
327  ComputeNonlinearResidual_Recursive(pde,residual,element,tria_element,prolong_matrix,coarse_index,fine_index,edc,fdc);
328  tria_element_iter++;
329  }
330  //Check if some preset righthandside exists.
331  AddPresetRightHandSide(-1.,residual);
332 
333  if (apply_boundary_values)
334  {
335  ApplyNewtonBoundaryValues(pde, residual);
336  }
337  }
338 
339  /*******************************************************************************************/
340 
341  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
342  int dim>
343  template<typename PROBLEM>
344  void
345  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearLhs(
346  PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
347  {
348  {
349  throw DOpEException("This function needs to be implemented!", "IntegratorMultiMesh::ComputeNonlinearLhs");
350  }
351  }
352 
353  /*******************************************************************************************/
354 
355  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
356  int dim>
357  template<typename PROBLEM>
358  void
359  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearRhs(
360  PROBLEM& pde, VECTOR &residual, bool apply_boundary_values)
361  {
362  residual = 0.;
363 
364  const auto& dof_handler =
365  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
366 
367  assert(dof_handler.size() == 2);
368 
369  const auto tria_element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler().get_tria(),
370  dof_handler[1]->GetDEALDoFHandler().get_tria());
371  const auto element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler(),
372  dof_handler[1]->GetDEALDoFHandler());
373 
374  auto element_iter = element_list.begin();
375  auto tria_element_iter = tria_element_list.begin();
376 
377  std::vector<decltype(tria_element_iter->first)> tria_element(2);
378  tria_element[0] = tria_element_iter->first;
379  tria_element[1] = tria_element_iter->second;
380 
381  std::vector<decltype(element_iter->first)> element(2);
382  element[0] = element_iter->first;
383  element[1] = element_iter->second;
384  int coarse_index = 0; //element[coarse_index] is the coarser of the two.
385  int fine_index = 0; //element[fine_index] is the finer of the two (both indices are 2 = element.size() if they are equally refined.
386 
387  // Generate the data containers.
388  idc_.InitializeMMEDC(pde.GetUpdateFlags(),
389  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element, tria_element,
390  this->GetParamData(), this->GetDomainData());
391  auto& edc = idc_.GetMultimeshElementDataContainer();
392 
393  bool need_interfaces = pde.HasInterfaces();
394  idc_.InitializeMMFDC(pde.GetFaceUpdateFlags(),
395  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
396  element,
397  tria_element,
398  this->GetParamData(),
399  this->GetDomainData(),
400  need_interfaces);
401  auto& fdc = idc_.GetMultimeshFaceDataContainer();
402 
403  for(; element_iter != element_list.end(); element_iter++)
404  {
405  element[0] = element_iter->first;
406  element[1] = element_iter->second;
407  tria_element[0] = tria_element_iter->first;
408  tria_element[1] = tria_element_iter->second;
409  FullMatrix<SCALAR> prolong_matrix;
410 
411  if(element[0]->has_children())
412  {
413  prolong_matrix = IdentityMatrix(element[1]->get_fe().dofs_per_cell);
414  coarse_index =1;
415  fine_index = 0;
416  }
417  else
418  {
419  if(element[1]->has_children())
420  {
421  prolong_matrix = IdentityMatrix(element[0]->get_fe().dofs_per_cell);
422  coarse_index =0;
423  fine_index = 1;
424  }
425  else
426  {
427  assert(element.size() ==2);
428  coarse_index = fine_index = 2;
429  }
430  }
431  ComputeNonlinearRhs_Recursive(pde,residual,element,tria_element,prolong_matrix,coarse_index,fine_index,edc,fdc);
432  tria_element_iter++;
433  }
434  //Check if some preset righthandside exists.
435  AddPresetRightHandSide(1.,residual);
436 
437  if (apply_boundary_values)
438  {
439  ApplyNewtonBoundaryValues(pde, residual);
440  }
441  }
442 
443  /*******************************************************************************************/
444 
445  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
446  int dim>
447  template<typename PROBLEM, typename MATRIX>
448  void
449  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeMatrix(
450  PROBLEM& pde, MATRIX &matrix)
451  {
452  matrix = 0.;
453 
454  const auto& dof_handler =
455  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
456 
457  assert(dof_handler.size() == 2);
458 
459  const auto tria_element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler().get_tria(),
460  dof_handler[1]->GetDEALDoFHandler().get_tria());
461  const auto element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler(),
462  dof_handler[1]->GetDEALDoFHandler());
463 
464  auto element_iter = element_list.begin();
465  auto tria_element_iter = tria_element_list.begin();
466 
467  std::vector<decltype(tria_element_iter->first)> tria_element(2);
468  tria_element[0] = tria_element_iter->first;
469  tria_element[1] = tria_element_iter->second;
470 
471  std::vector<decltype(element_iter->first)> element(2);
472  element[0] = element_iter->first;
473  element[1] = element_iter->second;
474  int coarse_index = 0; //element[coarse_index] is the coarser of the two.
475  int fine_index = 0; //element[fine_index] is the finer of the two (both indices are 2 = element.size() if they are equally refined.
476 
477  // Generate the data containers.
478  idc_.InitializeMMEDC(pde.GetUpdateFlags(),
479  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element, tria_element,
480  this->GetParamData(), this->GetDomainData());
481  auto& edc = idc_.GetMultimeshElementDataContainer();
482 
483  bool need_interfaces = pde.HasInterfaces();
484  idc_.InitializeMMFDC(pde.GetFaceUpdateFlags(),
485  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
486  element,
487  tria_element,
488  this->GetParamData(),
489  this->GetDomainData(),
490  need_interfaces);
491  auto & fdc = idc_.GetMultimeshFaceDataContainer();
492 
493  for(; element_iter != element_list.end(); element_iter++)
494  {
495  element[0] = element_iter->first;
496  element[1] = element_iter->second;
497  tria_element[0] = tria_element_iter->first;
498  tria_element[1] = tria_element_iter->second;
499  FullMatrix<SCALAR> prolong_matrix;
500 
501  if(element[0]->has_children())
502  {
503  prolong_matrix = IdentityMatrix(element[1]->get_fe().dofs_per_cell);
504  coarse_index =1;
505  fine_index = 0;
506  }
507  else
508  {
509  if(element[1]->has_children())
510  {
511  prolong_matrix = IdentityMatrix(element[0]->get_fe().dofs_per_cell);
512  coarse_index =0;
513  fine_index = 1;
514  }
515  else
516  {
517  assert(element.size() ==2);
518  coarse_index = fine_index = 2;
519  }
520  }
521  ComputeMatrix_Recursive(pde,matrix,element,tria_element,prolong_matrix,coarse_index,fine_index,edc,fdc);
522  tria_element_iter++;
523  }
524  }
525 
526  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
527  int dim>
528  template<typename PROBLEM>
529  SCALAR
530  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeDomainScalar(
531  PROBLEM& pde)
532  {
533  SCALAR ret = 0.;
534 
535  const auto& dof_handler =
536  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
537 
538  assert(dof_handler.size() == 2);
539 
540  const auto tria_element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler().get_tria(),
541  dof_handler[1]->GetDEALDoFHandler().get_tria());
542  const auto element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler(),
543  dof_handler[1]->GetDEALDoFHandler());
544 
545  auto element_iter = element_list.begin();
546  auto tria_element_iter = tria_element_list.begin();
547 
548  std::vector<decltype(tria_element_iter->first)> tria_element(2);
549  tria_element[0] = tria_element_iter->first;
550  tria_element[1] = tria_element_iter->second;
551 
552  std::vector<decltype(element_iter->first)> element(2);
553  element[0] = element_iter->first;
554  element[1] = element_iter->second;
555  int coarse_index = 0; //element[coarse_index] is the coarser of the two.
556  int fine_index = 0; //element[fine_index] is the finer of the two (both indices are 2 = element.size() if they are equally refined.
557 
558  if (pde.HasFaces())
559  {
560  throw DOpEException("This function should not be called when faces are needed!",
561  "IntegratorMultiMesh::ComputeDomainScalar");
562  }
563 
564  // Generate the data containers.
565  idc_.InitializeMMEDC(pde.GetUpdateFlags(),
566  *(pde.GetBaseProblem().GetSpaceTimeHandler()), element, tria_element,
567  this->GetParamData(), this->GetDomainData());
568  auto& edc = idc_.GetMultimeshElementDataContainer();
569 
570  for(; element_iter != element_list.end(); element_iter++)
571  {
572  element[0] = element_iter->first;
573  element[1] = element_iter->second;
574  tria_element[0] = tria_element_iter->first;
575  tria_element[1] = tria_element_iter->second;
576  FullMatrix<SCALAR> prolong_matrix;
577 
578  if(element[0]->has_children())
579  {
580  prolong_matrix = IdentityMatrix(element[1]->get_fe().dofs_per_cell);
581  coarse_index =1;
582  fine_index = 0;
583  }
584  else
585  {
586  if(element[1]->has_children())
587  {
588  prolong_matrix = IdentityMatrix(element[0]->get_fe().dofs_per_cell);
589  coarse_index =0;
590  fine_index = 1;
591  }
592  else
593  {
594  assert(element.size() ==2);
595  coarse_index = fine_index = 2;
596  }
597  }
598  ret += ComputeDomainScalar_Recursive(pde,element,tria_element,prolong_matrix,coarse_index,fine_index,edc);
599  tria_element_iter++;
600  }
601  return ret;
602  }
603 
604 
605  /*******************************************************************************************/
606 
607  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
608  int dim>
609  template<typename PROBLEM>
610  SCALAR
611  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputePointScalar(
612  PROBLEM& pde)
613  {
614 
615  {
616  SCALAR ret = 0.;
617  ret += pde.PointFunctional(this->GetParamData(),
618  this->GetDomainData());
619 
620  return ret;
621  }
622  }
623  /*******************************************************************************************/
624 
625  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
626  int dim>
627  template<typename PROBLEM>
628  SCALAR
629  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeBoundaryScalar(
630  PROBLEM& pde
631  )
632  {
633  SCALAR ret = 0.;
634  // Begin integration
635  const auto& dof_handler =
636  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
637 
638  const auto tria_element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler().get_tria(),
639  dof_handler[1]->GetDEALDoFHandler().get_tria());
640  const auto element_list = GridTools::get_finest_common_cells (dof_handler[0]->GetDEALDoFHandler(),
641  dof_handler[1]->GetDEALDoFHandler());
642 
643  auto element_iter = element_list.begin();
644  auto tria_element_iter = tria_element_list.begin();
645 
646  std::vector<decltype(tria_element_iter->first)> tria_element(2);
647  tria_element[0] = tria_element_iter->first;
648  tria_element[1] = tria_element_iter->second;
649 
650  std::vector<decltype(element_iter->first)> element(2);
651  element[0] = element_iter->first;
652  element[1] = element_iter->second;
653  int coarse_index = 0; //element[coarse_index] is the coarser of the two.
654  int fine_index = 0; //element[fine_index] is the finer of the two (both indices are 2 = element.size() if they are equally refined.
655 
656 
657  idc_.InitializeMMFDC(pde.GetFaceUpdateFlags(),
658  *(pde.GetBaseProblem().GetSpaceTimeHandler()),
659  element, tria_element,
660  this->GetParamData(),
661  this->GetDomainData());
662  auto & fdc = idc_.GetMultimeshFaceDataContainer();
663 
664  std::vector<unsigned int> boundary_functional_colors = pde.GetBoundaryFunctionalColors();
665  bool need_boundary_integrals = (boundary_functional_colors.size() > 0);
666  if(!need_boundary_integrals)
667  {
668  throw DOpEException("No boundary colors given!","IntegratorMultiMesh::ComputeBoundaryScalar");
669  }
670 
671  for(; element_iter != element_list.end(); element_iter++)
672  {
673  element[0] = element_iter->first;
674  element[1] = element_iter->second;
675  tria_element[0] = tria_element_iter->first;
676  tria_element[1] = tria_element_iter->second;
677  FullMatrix<SCALAR> prolong_matrix;
678 
679  if(element[0]->has_children())
680  {
681  prolong_matrix = IdentityMatrix(element[1]->get_fe().dofs_per_cell);
682  coarse_index =1;
683  fine_index = 0;
684  }
685  else
686  {
687  if(element[1]->has_children())
688  {
689  prolong_matrix = IdentityMatrix(element[0]->get_fe().dofs_per_cell);
690  coarse_index =0;
691  fine_index = 1;
692  }
693  else
694  {
695  assert(element.size() ==2);
696  coarse_index = fine_index = 2;
697  }
698  }
699  ret += ComputeBoundaryScalar_Recursive(pde,element,tria_element,prolong_matrix,coarse_index,fine_index,fdc);
700  tria_element_iter++;
701  }
702 
703  return ret;
704  }
705  /*******************************************************************************************/
706 
707  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
708  int dim>
709  template<typename PROBLEM>
710  SCALAR
711  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeFaceScalar(
712  PROBLEM& /*pde*/
713  )
714  {
715  throw DOpEException("This function needs to be implemented!", "IntegratorMultiMesh::ComputeFaceScalar");
716 // {
717 // SCALAR ret = 0.;
718 //#if deal_II_dimension == 2 || deal_II_dimension == 3
719 // // Begin integration
720 // const auto& dof_handler =
721 // pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler();
722 // auto element = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerBeginActive();
723 // auto endc = pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandlerEnd();
724 //
725 // idc_.InitializeFDC(pde.GetFaceUpdateFlags(),
726 // *(pde.GetBaseProblem().GetSpaceTimeHandler()),
727 // element,
728 // this->GetParamData(),
729 // this->GetDomainData());
730 // auto & fdc = idc_.GetFaceDataContainer();
731 //
732 // bool need_faces = pde.HasFaces();
733 // if(!need_faces)
734 // {
735 // throw DOpEException("No faces required!","IntegratorMultiMesh::ComputeFaceScalar");
736 // }
737 //
738 // for (;element[0]!=endc[0]; element[0]++)
739 // {
740 // for(unsigned int dh=1; dh<dof_handler.size(); dh++)
741 // {
742 // if( element[dh] == endc[dh])
743 // {
744 // throw DOpEException("Elementnumbers in DoFHandlers are not matching!","IntegratorMultiMesh::ComputeFaceScalar");
745 // }
746 // }
747 //
748 // if(need_faces)
749 // {
750 // for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
751 // {
752 // if (element[0]->neighbor_index(face) != -1)
753 // {
754 // fdc.ReInit(face);
755 // ret +=pde.FaceFunctional(fdc);
756 // }
757 // }
758 // }
759 // for(unsigned int dh=1; dh<dof_handler.size(); dh++)
760 // {
761 // element[dh]++;
762 // }
763 // }
764 //#else
765 // throw DOpEException("Not implemented in this dimension!",
766 // "IntegratorMultiMesh::ComputeFaceScalar");
767 //#endif
768 // return ret;
769 // }
770  }
771  /*******************************************************************************************/
772 
773  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
774  int dim>
775  template<typename PROBLEM>
776  SCALAR
777  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeAlgebraicScalar(
778  PROBLEM& pde)
779  {
780 
781  {
782  SCALAR ret = 0.;
783  ret = pde.AlgebraicFunctional(this->GetParamData(),
784  this->GetDomainData());
785  return ret;
786  }
787  }
788 
789  /*******************************************************************************************/
790 
791  template <typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,int dim>
792  template<typename PROBLEM>
793  void IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>
794  ::ComputeNonlinearAlgebraicResidual (PROBLEM& pde, VECTOR &residual)
795  {
796  residual = 0.;
797  pde.AlgebraicResidual(residual,this->GetParamData(),this->GetDomainData());
798  //Check if some preset righthandside exists.
799  AddPresetRightHandSide(-1.,residual);
800  }
801 
802  /*******************************************************************************************/
803 
804  template <typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,int dim>
805  template<typename PROBLEM>
806  void IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>
807  ::ComputeLocalControlConstraints (PROBLEM& pde, VECTOR &constraints)
808  {
809  constraints = 0.;
810  pde.ComputeLocalControlConstraints(constraints,this->GetParamData(),this->GetDomainData());
811  }
812 
813  /*******************************************************************************************/
814 
815  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
816  int dim>
817  template<typename PROBLEM>
818  void
819  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyInitialBoundaryValues(
820  PROBLEM& pde, VECTOR &u)
821  {
822  //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
823  // e.g. sometimes we need zero sometimes we need other values.
824 
825  //Never Condense Nodes Here ! Or All will fail if the state is not initialized with zero!
826  //pde.GetDoFConstraints().condense(u);
827  std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
828  for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
829  {
830  unsigned int color = dirichlet_colors[i];
831  std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
832  std::map<unsigned int, SCALAR> boundary_values;
833 
834  InterpolateBoundaryValues(pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0], color, pde.GetDirichletValues(color, this->GetParamData(),
835  this->GetDomainData()), boundary_values, comp_mask );
836 
837  for (typename std::map<unsigned int, SCALAR>::const_iterator p =
838  boundary_values.begin(); p != boundary_values.end(); p++)
839  {
840  u(p->first) = p->second;
841  }
842  }
843  }
844 
845  /*******************************************************************************************/
846 
847  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
848  int dim>
849  template<typename PROBLEM>
850  void
851  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyNewtonBoundaryValues(
852  PROBLEM& pde, VECTOR &u)
853  {
854  //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
855  // e.g. sometimes we need zero sometimes we need other values.
856 
857  pde.GetDoFConstraints().condense(u);
858  std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
859  for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
860  {
861  unsigned int color = dirichlet_colors[i];
862  std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
863  std::map<unsigned int, SCALAR> boundary_values;
864 
865 
866  InterpolateBoundaryValues(
867  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
868  color, dealii::ZeroFunction<dim>(comp_mask.size()),
869  boundary_values, comp_mask);
870 
871  for (typename std::map<unsigned int, SCALAR>::const_iterator p =
872  boundary_values.begin(); p != boundary_values.end(); p++)
873  {
874  u(p->first) = p->second;
875  }
876  }
877  }
878  /*******************************************************************************************/
879 
880  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
881  int dim>
882  template<typename PROBLEM, typename MATRIX>
883  void
884  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ApplyNewtonBoundaryValues(
885  PROBLEM& pde, MATRIX& matrix, VECTOR &rhs, VECTOR &sol)
886  {
887  //TODO Apply constraints locally, see, e.g., dealii step-27 ? But howto do this in the newton iter
888  // e.g. sometimes we need zero sometimes we need other values.
889  pde.GetDoFConstraints().condense(rhs);
890  pde.GetDoFConstraints().condense(matrix);
891  std::vector<unsigned int> dirichlet_colors = pde.GetDirichletColors();
892  for (unsigned int i = 0; i < dirichlet_colors.size(); i++)
893  {
894  unsigned int color = dirichlet_colors[i];
895  std::vector<bool> comp_mask = pde.GetDirichletCompMask(color);
896  std::map<unsigned int, SCALAR> boundary_values;
897 
898  InterpolateBoundaryValues(
899  pde.GetBaseProblem().GetSpaceTimeHandler()->GetDoFHandler()[0],
900  color, dealii::ZeroFunction<dim>(comp_mask.size()),
901  boundary_values, comp_mask);
902 
903  dealii::MatrixTools::apply_boundary_values(boundary_values, matrix,
904  sol, rhs);
905  }
906  }
907 
908  /*******************************************************************************************/
909 
910  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
911  int dim>
912  void
913  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddDomainData(
914  std::string name, const VECTOR* new_data)
915  {
916  if (domain_data_.find(name) != domain_data_.end())
917  {
918  throw DOpEException(
919  "Adding multiple Data with name " + name + " is prohibited!",
920  "IntegratorMultiMesh::AddDomainData");
921  }
922  domain_data_.insert(std::pair<std::string, const VECTOR*>(name, new_data));
923  }
924 
925  /*******************************************************************************************/
926 
927  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
928  int dim>
929  void
930  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::DeleteDomainData(
931  std::string name)
932  {
933  typename std::map<std::string, const VECTOR *>::iterator it =
934  domain_data_.find(name);
935  if (it == domain_data_.end())
936  {
937  throw DOpEException(
938  "Deleting Data " + name + " is impossible! Data not found",
939  "IntegratorMultiMesh::DeleteDomainData");
940  }
941  domain_data_.erase(it);
942  }
943  /*******************************************************************************************/
944 
945  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
946  int dim>
947  void
948  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddPresetRightHandSide(double s,
949  VECTOR &residual) const
950  {
951  typename std::map<std::string, const VECTOR *>::const_iterator it =
952  domain_data_.find("fixed_rhs");
953  if (it != domain_data_.end())
954  {
955  assert(residual.size() == it->second->size());
956  residual.add(s,*(it->second));
957  }
958  }
959 
960  /*******************************************************************************************/
961 
962  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
963  int dim>
964  const std::map<std::string, const VECTOR*>&
965  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetDomainData() const
966  {
967  return domain_data_;
968  }
969 
970  /*******************************************************************************************/
971 
972  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
973  int dim>
974  void
975  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::AddParamData(
976  std::string name, const dealii::Vector<SCALAR>* new_data)
977  {
978  if (param_data_.find(name) != param_data_.end())
979  {
980  throw DOpEException(
981  "Adding multiple Data with name " + name + " is prohibited!",
982  "IntegratorMultiMesh::AddParamData");
983  }
984  param_data_.insert(
985  std::pair<std::string, const dealii::Vector<SCALAR>*>(name, new_data));
986  }
987 
988  /*******************************************************************************************/
989 
990  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
991  int dim>
992  void
993  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::DeleteParamData(
994  std::string name)
995  {
996  typename std::map<std::string, const dealii::Vector<SCALAR>*>::iterator
997  it = param_data_.find(name);
998  if (it == param_data_.end())
999  {
1000  throw DOpEException(
1001  "Deleting Data " + name + " is impossible! Data not found",
1002  "IntegratorMultiMesh::DeleteParamData");
1003  }
1004  param_data_.erase(it);
1005  }
1006 
1007  /*******************************************************************************************/
1008 
1009  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1010  int dim>
1011  const std::map<std::string, const dealii::Vector<SCALAR>*>&
1012  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetParamData() const
1013  {
1014  return param_data_;
1015  }
1016 
1017  /*******************************************************************************************/
1018 
1019  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1020  int dim>
1021  INTEGRATORDATACONT&
1022  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::GetIntegratorDataContainer() const
1023  {
1024  return idc_;
1025  }
1026 
1027 
1028  /*******************************************************************************************/
1029 
1030  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1031  int dim>
1032  template<template<int, int> class DH>
1033  void
1034  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::InterpolateBoundaryValues(
1035  const DOpEWrapper::DoFHandler<dim, DH>* dof_handler,
1036  const unsigned int color, const dealii::Function<dim>& function,
1037  std::map<unsigned int, SCALAR>& boundary_values,
1038  const std::vector<bool>& comp_mask) const
1039  {
1040  dealii::VectorTools::interpolate_boundary_values(
1041  dof_handler->GetDEALDoFHandler(), color, function,
1042  boundary_values, comp_mask);
1043  }
1044 
1045  /*******************************************************************************************/
1046 
1047  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,int dim>
1048  template<typename PROBLEM, template<int, int> class DH>
1049  void IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearResidual_Recursive(
1050  PROBLEM& pde,
1051  VECTOR& residual,
1052  typename std::vector<typename DH<dim, dim>::cell_iterator> &element,
1053  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator> &tria_element,
1054  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
1055  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc,Multimesh_FaceDataContainer<DH, VECTOR, dim>& fdc)
1056  {
1057  if(!element[0]->has_children() && ! element[1]->has_children())
1058  {
1059  unsigned int dofs_per_element;
1060  dealii::Vector<SCALAR> local_vector;
1061  std::vector<unsigned int> local_dof_indices;
1062 
1063  bool need_faces = pde.HasFaces();
1064  std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
1065  bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
1066  bool need_interfaces = pde.HasInterfaces();
1067 
1068  edc.ReInit(coarse_index,fine_index,prolong_matrix);
1069 
1070  dofs_per_element = element[0]->get_fe().dofs_per_cell;
1071 
1072  local_vector.reinit(dofs_per_element);
1073  local_vector = 0;
1074 
1075  local_dof_indices.resize(0);
1076  local_dof_indices.resize(dofs_per_element, 0);
1077 
1078  //the second '1' plays only a role in the stationary case. In the non-stationary
1079  //case, scale_ico is set by the time-stepping-scheme
1080  pde.ElementEquation(edc, local_vector, 1., 1.);
1081  pde.ElementRhs(edc, local_vector, -1.);
1082 
1083  //FIXME Integrate on Faces of Element[0] that contain a fine-element face.
1084  if(need_faces || need_interfaces)
1085  {
1086  throw DOpEException(" Faces on multiple meshes not implemented yet!",
1087  "IntegratorMultiMesh::ComputeNonlinearResidual_Recursive");
1088  }
1089  unsigned int b_index = fine_index%2; //This takes care that if fine_index ==2 then we select the
1090  //zeros entry in the element vector
1091 
1092  if(need_boundary_integrals && element[b_index]->at_boundary())
1093  {
1094  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1095  {
1096 #if DEAL_II_VERSION_GTE(8,3,0)
1097  if (element[b_index]->face(face)->at_boundary()
1098  &&
1099  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1100  element[b_index]->face(face)->boundary_id()) != boundary_equation_colors.end()))
1101 #else
1102  if (element[b_index]->face(face)->at_boundary()
1103  &&
1104  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1105  element[b_index]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
1106 #endif
1107  {
1108  fdc.ReInit(coarse_index,fine_index,prolong_matrix,face);
1109  pde.BoundaryEquation(fdc,local_vector, 1., 1.);
1110  pde.BoundaryRhs(fdc,local_vector,-1.);
1111  }
1112  }
1113  }
1114  // if(need_faces)
1115  // {
1116  // for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1117  // {
1118  // if (element[0]->neighbor_index(face) != -1)
1119  // {
1120  // fdc.ReInit(face);
1121  // pde.FaceEquation(fdc, local_vector);
1122  // pde.FaceRhs(fdc, local_vector,-1.);
1123  // }
1124  // }
1125  // }
1126  // if( need_interfaces)
1127  // {
1128  // for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1129  // {
1130  // fdc.ReInit(face);
1131  // if (element[0]->neighbor_index(face) != -1
1132  // &&
1133  // fdc.GetMaterialId()!= fdc.GetNbrMaterialId())
1134  // {
1135  // fdc.ReInitNbr();
1136  // pde.InterfaceEquation(fdc, local_vector);
1137  // }
1138  // }
1139  // }
1140  if(coarse_index == 0) //Need to transfer the computed residual to the one on the coarse element
1141  {
1142  dealii::Vector<SCALAR> tmp(dofs_per_element);
1143  prolong_matrix.Tvmult(tmp,local_vector);
1144 
1145  //LocalToGlobal
1146  element[0]->get_dof_indices(local_dof_indices);
1147  for (unsigned int i = 0; i < dofs_per_element; ++i)
1148  {
1149  residual(local_dof_indices[i]) += tmp(i);
1150  }
1151  }
1152  else //Testfunctions are already the right ones...
1153  {
1154  //LocalToGlobal
1155  element[0]->get_dof_indices(local_dof_indices);
1156  for (unsigned int i = 0; i < dofs_per_element; ++i)
1157  {
1158  residual(local_dof_indices[i]) += local_vector(i);
1159  }
1160  }
1161 
1162  }//Endof the case on the finest level
1163  else
1164  {
1165  assert(fine_index != coarse_index);
1166  assert(element[fine_index]->has_children());
1167  assert(!element[coarse_index]->has_children());
1168  assert(tria_element[fine_index]->has_children());
1169  assert(!tria_element[coarse_index]->has_children());
1170 
1171  unsigned int local_n_dofs = element[coarse_index]->get_fe().dofs_per_cell;
1172 
1173  typename DH<dim, dim>::cell_iterator dofh_fine = element[fine_index];
1174  typename dealii::Triangulation<dim>::cell_iterator tria_fine = tria_element[fine_index];
1175 
1176  for (unsigned int child=0; child<GeometryInfo<dim>::max_children_per_cell;++child)
1177  {
1178  FullMatrix<SCALAR> new_matrix(local_n_dofs);
1179  element[coarse_index]->get_fe().get_prolongation_matrix(child).mmult (new_matrix,
1180  prolong_matrix);
1181  element[fine_index] = dofh_fine->child(child);
1182  tria_element[fine_index] = tria_fine->child(child);
1183 
1184  ComputeNonlinearResidual_Recursive(pde,residual,element,tria_element,new_matrix, coarse_index, fine_index, edc, fdc);
1185  }
1186  }
1187  }
1188  /*******************************************************************************************/
1189 
1190  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,int dim>
1191  template<typename PROBLEM, template<int, int> class DH>
1192  void IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeNonlinearRhs_Recursive(
1193  PROBLEM& pde,
1194  VECTOR& residual,
1195  typename std::vector<typename DH<dim, dim>::cell_iterator> &element,
1196  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator> &tria_element,
1197  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
1198  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc,Multimesh_FaceDataContainer<DH, VECTOR, dim>& fdc)
1199  {
1200  if(!element[0]->has_children() && ! element[1]->has_children())
1201  {
1202  unsigned int dofs_per_element;
1203  dealii::Vector<SCALAR> local_vector;
1204  std::vector<unsigned int> local_dof_indices;
1205 
1206  bool need_faces = pde.HasFaces();
1207  std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
1208  bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
1209  bool need_interfaces = pde.HasInterfaces();
1210 
1211  edc.ReInit(coarse_index,fine_index,prolong_matrix);
1212 
1213  dofs_per_element = element[0]->get_fe().dofs_per_cell;
1214 
1215  local_vector.reinit(dofs_per_element);
1216  local_vector = 0;
1217 
1218  local_dof_indices.resize(0);
1219  local_dof_indices.resize(dofs_per_element, 0);
1220 
1221  //the second '1' plays only a role in the stationary case. In the non-stationary
1222  //case, scale_ico is set by the time-stepping-scheme
1223  pde.ElementRhs(edc, local_vector, 1.);
1224 
1225  //FIXME Integrate on Faces of Element[0] that contain a fine-element face.
1226  if(need_faces )
1227  {
1228  throw DOpEException(" Faces on multiple meshes not implemented yet!",
1229  "IntegratorMultiMesh::ComputeNonlinearRhs_Recursive");
1230  }
1231  if(need_interfaces )
1232  {
1233  throw DOpEException(" Faces on multiple meshes not implemented yet!",
1234  "IntegratorMultiMesh::ComputeNonlinearRhs_Recursive");
1235  }
1236  unsigned int b_index = fine_index%2; //This takes care that if fine_index ==2 then we select the
1237  //zeros entry in the element vector
1238 
1239  if(need_boundary_integrals && element[b_index]->at_boundary())
1240  {
1241  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1242  {
1243 #if DEAL_II_VERSION_GTE(8,3,0)
1244  if (element[b_index]->face(face)->at_boundary()
1245  &&
1246  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1247  element[b_index]->face(face)->boundary_id()) != boundary_equation_colors.end()))
1248 #else
1249  if (element[b_index]->face(face)->at_boundary()
1250  &&
1251  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1252  element[b_index]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
1253 #endif
1254  {
1255  fdc.ReInit(coarse_index,fine_index,prolong_matrix,face);
1256  pde.BoundaryRhs(fdc,local_vector,1.);
1257  }
1258  }
1259  }
1260  // if(need_faces)
1261  // {
1262  // for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1263  // {
1264  // if (element[0]->neighbor_index(face) != -1)
1265  // {
1266  // fdc.ReInit(face);
1267  // pde.FaceRhs(fdc, local_vector,1.);
1268  // }
1269  // }
1270  // }
1271  if(coarse_index == 0) //Need to transfer the computed residual to the one on the coarse element
1272  {
1273  dealii::Vector<SCALAR> tmp(dofs_per_element);
1274  prolong_matrix.Tvmult(tmp,local_vector);
1275 
1276  //LocalToGlobal
1277  element[0]->get_dof_indices(local_dof_indices);
1278  for (unsigned int i = 0; i < dofs_per_element; ++i)
1279  {
1280  residual(local_dof_indices[i]) += tmp(i);
1281  }
1282  }
1283  else //Testfunctions are already the right ones...
1284  {
1285  //LocalToGlobal
1286  element[0]->get_dof_indices(local_dof_indices);
1287  for (unsigned int i = 0; i < dofs_per_element; ++i)
1288  {
1289  residual(local_dof_indices[i]) += local_vector(i);
1290  }
1291  }
1292 
1293  }//Endof the case on the finest level
1294  else
1295  {
1296  assert(fine_index != coarse_index);
1297  assert(element[fine_index]->has_children());
1298  assert(!element[coarse_index]->has_children());
1299  assert(tria_element[fine_index]->has_children());
1300  assert(!tria_element[coarse_index]->has_children());
1301 
1302  unsigned int local_n_dofs = element[coarse_index]->get_fe().dofs_per_cell;
1303 
1304  typename DH<dim, dim>::cell_iterator dofh_fine = element[fine_index];
1305  typename dealii::Triangulation<dim>::cell_iterator tria_fine = tria_element[fine_index];
1306 
1307  for (unsigned int child=0; child<GeometryInfo<dim>::max_children_per_cell;++child)
1308  {
1309  FullMatrix<SCALAR> new_matrix(local_n_dofs);
1310  element[coarse_index]->get_fe().get_prolongation_matrix(child).mmult (new_matrix,
1311  prolong_matrix);
1312  element[fine_index] = dofh_fine->child(child);
1313  tria_element[fine_index] = tria_fine->child(child);
1314 
1315  ComputeNonlinearRhs_Recursive(pde,residual,element,tria_element,new_matrix, coarse_index, fine_index, edc, fdc);
1316  }
1317  }
1318  }
1319 
1320  /*******************************************************************************************/
1321 
1322  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,
1323  int dim>
1324  template<typename PROBLEM, typename MATRIX, template<int, int> class DH>
1325  void
1326  IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeMatrix_Recursive(
1327  PROBLEM& pde, MATRIX &matrix, typename std::vector<typename DH<dim, dim>::cell_iterator>& element,
1328  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator>& tria_element,
1329  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
1330  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc,
1331  Multimesh_FaceDataContainer<DH, VECTOR, dim>& fdc)
1332  {
1333 
1334  if(!element[0]->has_children() && ! element[1]->has_children())
1335  {
1336  unsigned int dofs_per_element;
1337  std::vector<unsigned int> local_dof_indices;
1338 
1339  bool need_faces = pde.HasFaces();
1340  bool need_interfaces = pde.HasInterfaces();
1341  std::vector<unsigned int> boundary_equation_colors = pde.GetBoundaryEquationColors();
1342  bool need_boundary_integrals = (boundary_equation_colors.size() > 0);
1343 
1344  edc.ReInit(coarse_index,fine_index,prolong_matrix);
1345  dofs_per_element = element[0]->get_fe().dofs_per_cell;
1346 
1347  dealii::FullMatrix<SCALAR> local_matrix(dofs_per_element,
1348  dofs_per_element);
1349  local_matrix = 0;
1350 
1351  local_dof_indices.resize(0);
1352  local_dof_indices.resize(dofs_per_element, 0);
1353  pde.ElementMatrix(edc, local_matrix);
1354 
1355  //FIXME Integrate on Faces of Element[0] that contain a fine-element face.
1356  if(need_faces || need_interfaces)
1357  {
1358  throw DOpEException(" Faces on multiple meshes not implemented yet!",
1359  "IntegratorMultiMesh::ComputeMatrix_Recursive");
1360  }
1361  unsigned int b_index = fine_index%2; //This takes care that if fine_index ==2 then we select the
1362  //zeros entry in the element vector
1363  if(need_boundary_integrals && element[b_index]->at_boundary())
1364  {
1365 
1366  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1367  {
1368 #if DEAL_II_VERSION_GTE(8,3,0)
1369  if (element[b_index]->face(face)->at_boundary()
1370  &&
1371  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1372  element[b_index]->face(face)->boundary_id()) != boundary_equation_colors.end()))
1373 #else
1374  if (element[b_index]->face(face)->at_boundary()
1375  &&
1376  (find(boundary_equation_colors.begin(),boundary_equation_colors.end(),
1377  element[b_index]->face(face)->boundary_indicator()) != boundary_equation_colors.end()))
1378 #endif
1379  {
1380  fdc.ReInit(coarse_index,fine_index,prolong_matrix,face);
1381  pde.BoundaryMatrix(fdc, local_matrix);
1382  }
1383  }
1384  }
1385  // if(need_faces)
1386  // {
1387  // for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1388  // {
1389  // if (element[0]->neighbor_index(face) != -1)
1390  // {
1391  // fdc.ReInit(face);
1392  // pde.FaceMatrix(fdc, local_matrix);
1393  // }
1394  // }
1395  // }
1396  // if( need_interfaces)
1397  // {
1398  // for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1399  // {
1400  // fdc.ReInit(face);
1401  // if (element[0]->neighbor_index(face) != -1
1402  // &&
1403  // fdc.GetMaterialId()!= fdc.GetNbrMaterialId())
1404  // {
1405  // fdc.ReInitNbr();
1406  //
1407  // nbr_dofs_per_element = fdc.GetNbrNDoFsPerElement();
1408  // nbr_local_dof_indices.resize(0);
1409  // nbr_local_dof_indices.resize(nbr_dofs_per_element, 0);
1410  // dealii::FullMatrix<SCALAR> local_interface_matrix(dofs_per_element,nbr_dofs_per_element );
1411  // local_interface_matrix = 0;
1412  //
1413  // pde.InterfaceMatrix(fdc, local_interface_matrix);
1414  //
1415  // element[0]->get_dof_indices(local_dof_indices);
1416  // element[0]->neighbor(face)->get_dof_indices(nbr_local_dof_indices);
1417  //
1418  // for (unsigned int i = 0; i < dofs_per_element; ++i)
1419  // {
1420  // for (unsigned int j = 0; j < nbr_dofs_per_element; ++j)
1421  // {
1422  // matrix.add(local_dof_indices[i], nbr_local_dof_indices[j],
1423  // local_interface_matrix(i, j));
1424  // } //endfor j
1425  // } //endfor i
1426  // }
1427  // }
1428  // }
1429  if(coarse_index == 0) //Need to transfer the computed residual to the one on the coarse element
1430  {
1431  dealii::FullMatrix<SCALAR> tmp(dofs_per_element);
1432  tmp = 0.;
1433  prolong_matrix.Tmmult(tmp,local_matrix);
1434  local_matrix = 0.;
1435  tmp.mmult(local_matrix,prolong_matrix);
1436 
1437  //LocalToGlobal
1438  element[0]->get_dof_indices(local_dof_indices);
1439  for (unsigned int i = 0; i < dofs_per_element; ++i)
1440  {
1441  for (unsigned int j = 0; j < dofs_per_element; ++j)
1442  {
1443  matrix.add(local_dof_indices[i], local_dof_indices[j],
1444  local_matrix(i, j));
1445  }
1446  }
1447  }
1448  else
1449  {
1450  //LocalToGlobal
1451  element[0]->get_dof_indices(local_dof_indices);
1452  for (unsigned int i = 0; i < dofs_per_element; ++i)
1453  {
1454  for (unsigned int j = 0; j < dofs_per_element; ++j)
1455  {
1456  matrix.add(local_dof_indices[i], local_dof_indices[j],
1457  local_matrix(i, j));
1458  }
1459  }
1460  }
1461 
1462  }//Endof the case on the finest level
1463  else
1464  {
1465  assert(fine_index != coarse_index);
1466  assert(element[fine_index]->has_children());
1467  assert(!element[coarse_index]->has_children());
1468  assert(tria_element[fine_index]->has_children());
1469  assert(!tria_element[coarse_index]->has_children());
1470 
1471  unsigned int local_n_dofs = element[coarse_index]->get_fe().dofs_per_cell;
1472 
1473  typename DH<dim, dim>::cell_iterator dofh_fine = element[fine_index];
1474  typename dealii::Triangulation<dim>::cell_iterator tria_fine = tria_element[fine_index];
1475 
1476  for (unsigned int child=0; child<GeometryInfo<dim>::max_children_per_cell;++child)
1477  {
1478  FullMatrix<SCALAR> new_matrix(local_n_dofs);
1479  element[coarse_index]->get_fe().get_prolongation_matrix(child).mmult (new_matrix,
1480  prolong_matrix);
1481  element[fine_index] = dofh_fine->child(child);
1482  tria_element[fine_index] = tria_fine->child(child);
1483 
1484  ComputeMatrix_Recursive(pde,matrix,element,tria_element,new_matrix, coarse_index, fine_index, edc, fdc);
1485  }
1486  }
1487  }
1488 
1489  /*******************************************************************************************/
1490 
1491  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,int dim>
1492  template<typename PROBLEM, template<int, int> class DH>
1493  SCALAR IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeDomainScalar_Recursive(
1494  PROBLEM& pde,
1495  typename std::vector<typename DH<dim, dim>::cell_iterator> &element,
1496  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator> &tria_element,
1497  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
1498  Multimesh_ElementDataContainer<DH, VECTOR, dim>& edc)
1499  {
1500  if(!element[0]->has_children() && ! element[1]->has_children())
1501  {
1502  SCALAR ret = 0.;
1503  edc.ReInit(coarse_index,fine_index,prolong_matrix);
1504  ret += pde.ElementFunctional(edc);
1505  return ret;
1506  } //Endof the case on the finest level
1507  else
1508  {
1509  assert(fine_index != coarse_index);
1510  assert(element[fine_index]->has_children());
1511  assert(!element[coarse_index]->has_children());
1512  assert(tria_element[fine_index]->has_children());
1513  assert(!tria_element[coarse_index]->has_children());
1514 
1515  unsigned int local_n_dofs = element[coarse_index]->get_fe().dofs_per_cell;
1516 
1517  typename DH<dim, dim>::cell_iterator dofh_fine = element[fine_index];
1518  typename dealii::Triangulation<dim>::cell_iterator tria_fine = tria_element[fine_index];
1519  SCALAR ret = 0.;
1520  for (unsigned int child=0; child<GeometryInfo<dim>::max_children_per_cell;++child)
1521  {
1522  FullMatrix<SCALAR> new_matrix(local_n_dofs);
1523  element[coarse_index]->get_fe().get_prolongation_matrix(child).mmult (new_matrix,
1524  prolong_matrix);
1525  element[fine_index] = dofh_fine->child(child);
1526  tria_element[fine_index] = tria_fine->child(child);
1527 
1528  ret += ComputeDomainScalar_Recursive(pde,element,tria_element,new_matrix, coarse_index, fine_index, edc);
1529  }
1530  return ret;
1531  }
1532  }
1533  /*******************************************************************************************/
1534 
1535  template<typename INTEGRATORDATACONT, typename VECTOR, typename SCALAR,int dim>
1536  template<typename PROBLEM, template<int, int> class DH>
1537  SCALAR IntegratorMultiMesh<INTEGRATORDATACONT, VECTOR, SCALAR, dim>::ComputeBoundaryScalar_Recursive(
1538  PROBLEM& pde,
1539  typename std::vector<typename DH<dim, dim>::cell_iterator> &element,
1540  typename std::vector<typename dealii::Triangulation<dim>::cell_iterator> &tria_element,
1541  const FullMatrix<SCALAR>& prolong_matrix,unsigned int coarse_index,unsigned int fine_index,
1542  Multimesh_FaceDataContainer<DH, VECTOR, dim>& fdc)
1543  {
1544  if(!element[0]->has_children() && ! element[1]->has_children())
1545  {
1546  SCALAR ret = 0.;
1547  std::vector<unsigned int> boundary_functional_colors = pde.GetBoundaryFunctionalColors();
1548  unsigned int b_index = fine_index%2; //This takes care that if fine_index ==2 then we select the
1549  //zeros entry in the element vector
1550  for (unsigned int face=0; face < dealii::GeometryInfo<dim>::faces_per_cell; ++face)
1551  {
1552 #if DEAL_II_VERSION_GTE(8,3,0)
1553  if (element[b_index]->face(face)->at_boundary()
1554  &&
1555  (find(boundary_functional_colors.begin(),boundary_functional_colors.end(),
1556  element[b_index]->face(face)->boundary_id()) != boundary_functional_colors.end()))
1557 #else
1558  if (element[b_index]->face(face)->at_boundary()
1559  &&
1560  (find(boundary_functional_colors.begin(),boundary_functional_colors.end(),
1561  element[b_index]->face(face)->boundary_indicator()) != boundary_functional_colors.end()))
1562 #endif
1563  {
1564  fdc.ReInit(coarse_index,fine_index,prolong_matrix,face);
1565  ret += pde.BoundaryFunctional(fdc);
1566  }
1567  }
1568  return ret;
1569  } //Endof the case on the finest level
1570  else
1571  {
1572  assert(fine_index != coarse_index);
1573  assert(element[fine_index]->has_children());
1574  assert(!element[coarse_index]->has_children());
1575  assert(tria_element[fine_index]->has_children());
1576  assert(!tria_element[coarse_index]->has_children());
1577 
1578  unsigned int local_n_dofs = element[coarse_index]->get_fe().dofs_per_cell;
1579 
1580  typename DH<dim, dim>::cell_iterator dofh_fine = element[fine_index];
1581  typename dealii::Triangulation<dim>::cell_iterator tria_fine = tria_element[fine_index];
1582  SCALAR ret = 0.;
1583  for (unsigned int child=0; child<GeometryInfo<dim>::max_children_per_cell;++child)
1584  {
1585  FullMatrix<SCALAR> new_matrix(local_n_dofs);
1586  element[coarse_index]->get_fe().get_prolongation_matrix(child).mmult (new_matrix,
1587  prolong_matrix);
1588  element[fine_index] = dofh_fine->child(child);
1589  tria_element[fine_index] = tria_fine->child(child);
1590 
1591  ret += ComputeBoundaryScalar_Recursive(pde,element,tria_element,new_matrix, coarse_index, fine_index, fdc);
1592  }
1593  return ret;
1594  }
1595  }
1596 //ENDOF NAMESPACE DOpE
1597 }
1598 #endif
1599 
DOpE::IntegratorMultiMesh::GetDomainData
const std::map< std::string, const VECTOR * > & GetDomainData() const
Definition: integrator_multimesh.h:965
DOpE::IntegratorMultiMesh::ComputeLocalControlConstraints
void ComputeLocalControlConstraints(PROBLEM &pde, VECTOR &constraints)
Definition: integrator_multimesh.h:807
multimesh_facedatacontainer.h
DOpE::IntegratorMultiMesh::ComputeBoundaryScalar
SCALAR ComputeBoundaryScalar(PROBLEM &pde)
Definition: integrator_multimesh.h:629
DOpE::IntegratorMultiMesh
Definition: integrator_multimesh.h:61
DOpE::IntegratorMultiMesh::ComputeFaceScalar
SCALAR ComputeFaceScalar(PROBLEM &pde)
Definition: integrator_multimesh.h:711
DOpE::IntegratorMultiMesh::ComputeDomainScalar
SCALAR ComputeDomainScalar(PROBLEM &pde)
Definition: integrator_multimesh.h:530
DOpE::IntegratorMultiMesh::ComputeNonlinearLhs
void ComputeNonlinearLhs(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integrator_multimesh.h:345
DOpE::Multimesh_ElementDataContainer
Definition: multimesh_elementdatacontainer.h:57
DOpE::IntegratorMultiMesh::ComputeNonlinearRhs
void ComputeNonlinearRhs(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integrator_multimesh.h:359
DOpE::IntegratorMultiMesh::ComputeNonlinearResidual
void ComputeNonlinearResidual(PROBLEM &pde, VECTOR &residual, bool apply_boundary_values=true)
Definition: integrator_multimesh.h:255
DOpE::IntegratorMultiMesh::ComputeAlgebraicScalar
SCALAR ComputeAlgebraicScalar(PROBLEM &pde)
Definition: integrator_multimesh.h:777
DOpE::IntegratorMultiMesh::AddParamData
void AddParamData(std::string name, const dealii::Vector< SCALAR > *new_data)
Definition: integrator_multimesh.h:975
DOpEWrapper::DoFHandler::GetDEALDoFHandler
const DOFHANDLER< dim, dim > & GetDEALDoFHandler() const
Definition: dofhandler_wrapper.h:69
DOpE::IntegratorMultiMesh::ComputePointScalar
SCALAR ComputePointScalar(PROBLEM &pde)
Definition: integrator_multimesh.h:611
DOpE::IntegratorMultiMesh::ApplyNewtonBoundaryValues
void ApplyNewtonBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integrator_multimesh.h:851
DOpE::IntegratorMultiMesh::ComputeMatrix
void ComputeMatrix(PROBLEM &pde, MATRIX &matrix)
Definition: integrator_multimesh.h:449
DOpE::IntegratorMultiMesh::~IntegratorMultiMesh
~IntegratorMultiMesh()
Definition: integrator_multimesh.h:234
DOpE::IntegratorMultiMesh::DeleteParamData
void DeleteParamData(std::string name)
Definition: integrator_multimesh.h:993
DOpE::IntegratorMultiMesh::ReInit
void ReInit()
Definition: integrator_multimesh.h:244
DOpE::IntegratorMultiMesh::DeleteDomainData
void DeleteDomainData(std::string name)
Definition: integrator_multimesh.h:930
multimesh_elementdatacontainer.h
DOpE::IntegratorMultiMesh::ApplyInitialBoundaryValues
void ApplyInitialBoundaryValues(PROBLEM &pde, VECTOR &u)
Definition: integrator_multimesh.h:819
DOpE::IntegratorMultiMesh::IntegratorMultiMesh
IntegratorMultiMesh(INTEGRATORDATACONT &idc)
Definition: integrator_multimesh.h:224
DOpE::IntegratorMultiMesh::AddDomainData
void AddDomainData(std::string name, const VECTOR *new_data)
Definition: integrator_multimesh.h:913
DOpE::IntegratorMultiMesh::ComputeNonlinearAlgebraicResidual
void ComputeNonlinearAlgebraicResidual(PROBLEM &pde, VECTOR &residual)
Definition: integrator_multimesh.h:794
DOpE::Multimesh_FaceDataContainer
Definition: multimesh_facedatacontainer.h:53
DOpE::IntegratorMultiMesh::AddPresetRightHandSide
void AddPresetRightHandSide(double s, VECTOR &residual) const
Definition: integrator_multimesh.h:948
DOpE::IntegratorMultiMesh::GetIntegratorDataContainer
INTEGRATORDATACONT & GetIntegratorDataContainer() const
Definition: integrator_multimesh.h:1022
DOpE::DOpEException
Definition: dopeexception.h:35
DOpE::IntegratorMultiMesh::GetParamData
const std::map< std::string, const dealii::Vector< SCALAR > * > & GetParamData() const
Definition: integrator_multimesh.h:1012
Generated on Mon Jan 25 2016 11:32:19 for DOpE by   doxygen 1.8.6