DUBeat  1.0.1
High-order discontinuous Galerkin methods and applications to cardiac electrophysiology
assemble_DG.hpp
Go to the documentation of this file.
1 /********************************************************************************
2  Copyright (C) 2024 by the DUBeat authors.
3 
4  This file is part of DUBeat.
5 
6  DUBeat is free software; you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  DUBeat is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with DUBeat. If not, see <http://www.gnu.org/licenses/>.
18 ********************************************************************************/
19 
27 #ifndef ASSEMBLE_DG_HPP_
28 #define ASSEMBLE_DG_HPP_
29 
30 #include <deal.II/base/parameter_handler.h>
31 
32 #include <deal.II/fe/mapping_q1_eulerian.h>
33 
34 #include <deal.II/lac/full_matrix.h>
35 
36 #include <memory>
37 #include <utility>
38 #include <vector>
39 
40 #include "DUB_FEM_handler.hpp"
41 #include "face_handler_DG.hpp"
42 #include "source/core_model.hpp"
43 #include "source/geometry/mesh_handler.hpp"
44 #include "source/init.hpp"
45 #include "source/io/data_writer.hpp"
46 #include "source/numerics/bc_handler.hpp"
47 #include "source/numerics/linear_solver_handler.hpp"
48 #include "source/numerics/preconditioner_handler.hpp"
49 #include "source/numerics/tools.hpp"
50 #include "volume_handler_DG.hpp"
51 
67 template <class basis>
68 class AssembleDG
69 {
70 protected:
72  const std::unique_ptr<basis> basis_ptr;
73 
75  VolumeHandlerDG<lifex::dim> vol_handler;
76 
78  FaceHandlerDG<lifex::dim> face_handler;
79 
81  FaceHandlerDG<lifex::dim> face_handler_neigh;
82 
84  unsigned int dofs_per_cell;
85 
87  const unsigned int n_quad_points_1D;
88 
90  const unsigned int n_quad_points;
91 
93  const unsigned int n_quad_points_face;
94 
96  const unsigned int poly_degree;
97 
99  typename dealii::DoFHandler<lifex::dim>::active_cell_iterator cell;
100 
101 public:
103  AssembleDG<basis>(const unsigned int degree)
104  : basis_ptr(std::make_unique<basis>(degree))
105  , vol_handler(VolumeHandlerDG<lifex::dim>(degree))
106  , face_handler(FaceHandlerDG<lifex::dim>(degree))
107  , face_handler_neigh(FaceHandlerDG<lifex::dim>(degree))
108  , n_quad_points_1D(degree + 2)
109  , n_quad_points(static_cast<int>(std::pow(n_quad_points_1D, lifex::dim)))
110  , n_quad_points_face(
111  static_cast<int>(std::pow(n_quad_points_1D, lifex::dim - 1)))
112  , poly_degree(degree)
113  {
114  dofs_per_cell = this->get_dofs_per_cell();
115  }
116 
118  AssembleDG<basis>(AssembleDG<basis> &) = default;
119 
121  AssembleDG<basis>(const AssembleDG<basis> &) = default;
122 
124  AssembleDG<basis>(AssembleDG<basis> &&) = default;
125 
127  void
128  reinit(const typename dealii::DoFHandler<lifex::dim>::active_cell_iterator
129  &new_cell,
130  const unsigned int new_edge);
131 
133  void
134  reinit(const typename dealii::DoFHandler<lifex::dim>::active_cell_iterator
135  &new_cell);
136 
138  unsigned int
139  get_dofs_per_cell() const;
140 
143  dealii::FullMatrix<double>
144  local_V() const;
145 
149  dealii::FullMatrix<double>
150  local_V(const dealii::Tensor<2, lifex::dim> &Sigma) const;
151 
155  dealii::FullMatrix<double>
156  local_M() const;
157 
163  dealii::Vector<double>
164  local_rhs(const std::shared_ptr<dealii::Function<lifex::dim>> &f_ex) const;
165 
171  dealii::Vector<double>
172  local_rhs_edge_dirichlet(
173  const double stability_coefficient,
174  const double theta,
175  const std::shared_ptr<lifex::utils::FunctionDirichlet> &u_ex) const;
176 
180  dealii::Vector<double>
181  local_rhs_edge_neumann(
182  const std::shared_ptr<dealii::Function<lifex::dim>> &g) const;
183 
188  dealii::FullMatrix<double>
189  local_SC(const double stability_coefficient) const;
190 
196  dealii::FullMatrix<double>
197  local_SC(const double stability_coefficient,
198  const dealii::Tensor<2, lifex::dim> &Sigma) const;
199 
204  dealii::FullMatrix<double>
205  local_SN(const double stability_coefficient) const;
206 
212  dealii::FullMatrix<double>
213  local_SN(const double stability_coefficient,
214  const dealii::Tensor<2, lifex::dim> &Sigma) const;
215 
225  std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
226  local_IC(const double theta) const;
227 
240  std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
241  local_IC(const double theta,
242  const dealii::Tensor<2, lifex::dim> &Sigma) const;
243 
251  std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
252  local_IB(const double theta) const;
253 
264  std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
265  local_IB(const double theta,
266  const dealii::Tensor<2, lifex::dim> &Sigma) const;
267 
275  std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
276  local_IN(const double theta) const;
277 
287  std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
288  local_IN(const double theta,
289  const dealii::Tensor<2, lifex::dim> &Sigma) const;
290 
295  dealii::Vector<double>
296  local_u0_M_rhs(
297  const lifex::LinAlg::MPI::Vector &u0,
298  const std::vector<dealii::types::global_dof_index> &dof_indices) const;
299 
306  dealii::Vector<double>
307  local_w0_M_rhs(
308  const lifex::LinAlg::MPI::Vector &u0,
309  const std::vector<dealii::types::global_dof_index> &dof_indices) const;
310 
316  dealii::FullMatrix<double>
317  local_non_linear_fitzhugh(
318  const lifex::LinAlg::MPI::Vector &u0,
319  const double a,
320  const std::vector<dealii::types::global_dof_index> &dof_indices) const;
321 
323  virtual ~AssembleDG() = default;
324 };
325 
326 template <class basis>
327 void
328 AssembleDG<basis>::reinit(
329  const typename dealii::DoFHandler<lifex::dim>::active_cell_iterator &new_cell,
330  const unsigned int new_edge)
331 {
332  cell = new_cell;
333  face_handler.reinit(new_cell, new_edge);
334 
335  if (!cell->at_boundary(new_edge))
336  {
337  const auto neighcell = cell->neighbor(new_edge);
338  const auto neighedge = cell->neighbor_face_no(new_edge);
339 
340  face_handler_neigh.reinit(neighcell, neighedge);
341  }
342 }
343 
344 template <class basis>
345 void
346 AssembleDG<basis>::reinit(
347  const typename dealii::DoFHandler<lifex::dim>::active_cell_iterator &new_cell)
348 {
349  cell = new_cell;
350  vol_handler.reinit(new_cell);
351 }
352 
353 template <class basis>
354 unsigned int
355 AssembleDG<basis>::get_dofs_per_cell() const
356 {
357  // The analytical formula is:
358  // n_dof_per_cell = (p+1)*(p+2)*...(p+d) / d!,
359  // where p is the space order and d the space dimension..
360 
361  unsigned int denominator = 1;
362  unsigned int nominator = 1;
363 
364  for (unsigned int i = 1; i <= lifex::dim; i++)
365  {
366  denominator *= i;
367  nominator *= poly_degree + i;
368  }
369 
370  return (int)(nominator / denominator);
371 }
372 
373 
374 template <class basis>
375 dealii::FullMatrix<double>
376 AssembleDG<basis>::local_V() const
377 {
378  dealii::FullMatrix<double> V(dofs_per_cell, dofs_per_cell);
379  const dealii::Tensor<2, lifex::dim> BJinv =
380  vol_handler.get_jacobian_inverse();
381  const double det = 1 / determinant(BJinv);
382 
383  for (unsigned int q = 0; q < n_quad_points; ++q)
384  {
385  for (unsigned int i = 0; i < dofs_per_cell; ++i)
386  {
387  for (unsigned int j = 0; j < dofs_per_cell; ++j)
388  {
389  V(i, j) +=
390  (basis_ptr->shape_grad(i, vol_handler.quadrature_ref(q)) *
391  BJinv) *
392  (basis_ptr->shape_grad(j, vol_handler.quadrature_ref(q)) *
393  BJinv) *
394  vol_handler.quadrature_weight(q) * det;
395  }
396  }
397  }
398 
399  return V;
400 }
401 
402 
403 template <class basis>
404 dealii::FullMatrix<double>
405 AssembleDG<basis>::local_V(const dealii::Tensor<2, lifex::dim> &Sigma) const
406 {
407  dealii::FullMatrix<double> V(dofs_per_cell, dofs_per_cell);
408  const dealii::Tensor<2, lifex::dim> BJinv =
409  vol_handler.get_jacobian_inverse();
410  const double det = 1 / determinant(BJinv);
411 
412  const dealii::Tensor<2, lifex::dim> Sigma_tr = transpose(Sigma);
413 
414  for (unsigned int q = 0; q < n_quad_points; ++q)
415  {
416  for (unsigned int i = 0; i < dofs_per_cell; ++i)
417  {
418  for (unsigned int j = 0; j < dofs_per_cell; ++j)
419  {
420  V(i, j) +=
421  (basis_ptr->shape_grad(i, vol_handler.quadrature_ref(q)) *
422  BJinv) *
423  (basis_ptr->shape_grad(j, vol_handler.quadrature_ref(q)) *
424  BJinv * Sigma_tr) *
425  vol_handler.quadrature_weight(q) * det;
426  }
427  }
428  }
429 
430  return V;
431 }
432 
433 template <class basis>
434 dealii::FullMatrix<double>
435 AssembleDG<basis>::local_M() const
436 {
437  dealii::FullMatrix<double> M(dofs_per_cell, dofs_per_cell);
438  const dealii::Tensor<2, lifex::dim> BJinv =
439  vol_handler.get_jacobian_inverse();
440  const double det = 1 / determinant(BJinv);
441 
442  for (unsigned int q = 0; q < n_quad_points; ++q)
443  {
444  for (unsigned int i = 0; i < dofs_per_cell; ++i)
445  {
446  for (unsigned int j = 0; j < dofs_per_cell; ++j)
447  {
448  M(i, j) +=
449  (basis_ptr->shape_value(i, vol_handler.quadrature_ref(q))) *
450  (basis_ptr->shape_value(j, vol_handler.quadrature_ref(q))) *
451  vol_handler.quadrature_weight(q) * det;
452  }
453  }
454  }
455 
456  return M;
457 }
458 
459 template <class basis>
460 dealii::Vector<double>
461 AssembleDG<basis>::local_rhs(
462  const std::shared_ptr<dealii::Function<lifex::dim>> &f_ex) const
463 {
464  AssertThrow(f_ex != nullptr,
465  dealii::StandardExceptions::ExcMessage(
466  "Not valid pointer to the source term."));
467 
468  dealii::Vector<double> cell_rhs(dofs_per_cell);
469  const dealii::Tensor<2, lifex::dim> BJinv =
470  vol_handler.get_jacobian_inverse();
471  const double det = 1 / determinant(BJinv);
472 
473  for (unsigned int q = 0; q < n_quad_points; ++q)
474  {
475  for (unsigned int i = 0; i < dofs_per_cell; ++i)
476  {
477  cell_rhs(i) +=
478  f_ex->value(vol_handler.quadrature_real(q)) *
479  basis_ptr->shape_value(i, vol_handler.quadrature_ref(q)) *
480  vol_handler.quadrature_weight(q) * det;
481  }
482  }
483 
484  return cell_rhs;
485 }
486 
487 template <class basis>
488 dealii::Vector<double>
489 AssembleDG<basis>::local_rhs_edge_dirichlet(
490  const double stability_coefficient,
491  const double theta,
492  const std::shared_ptr<lifex::utils::FunctionDirichlet> &u_ex) const
493 {
494  AssertThrow(theta == 1. || theta == 0. || theta == -1.,
495  dealii::StandardExceptions::ExcMessage(
496  "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
497  "or -1 (NIP method)."));
498 
499  AssertThrow(u_ex != nullptr,
500  dealii::StandardExceptions::ExcMessage(
501  "Not valid pointer to the exact solution."));
502 
503  dealii::Vector<double> cell_rhs(dofs_per_cell);
504  const dealii::Tensor<2, lifex::dim> BJinv =
505  face_handler.get_jacobian_inverse();
506 
507  const double face_measure = face_handler.get_measure();
508  const double unit_face_measure = (4.0 - lifex::dim) / 2;
509  const double measure_ratio = face_measure / unit_face_measure;
510  const double h_local = (cell->measure()) / face_measure;
511 
512  const double local_pen_coeff =
513  (stability_coefficient * pow(poly_degree, 2)) / h_local;
514 
515  for (unsigned int q = 0; q < n_quad_points_face; ++q)
516  {
517  for (unsigned int i = 0; i < dofs_per_cell; ++i)
518  {
519  cell_rhs(i) +=
520  local_pen_coeff *
521  basis_ptr->shape_value(i, face_handler.quadrature_ref(q)) *
522  u_ex->value(face_handler.quadrature_real(q)) *
523  face_handler.quadrature_weight(q) * measure_ratio;
524 
525  cell_rhs(i) -=
526  theta *
527  ((basis_ptr->shape_grad(i, face_handler.quadrature_ref(q)) *
528  BJinv) *
529  face_handler.get_normal()) *
530  u_ex->value(face_handler.quadrature_real(q)) *
531  face_handler.quadrature_weight(q) * measure_ratio;
532  }
533  }
534 
535  return cell_rhs;
536 }
537 
538 template <class basis>
539 dealii::Vector<double>
540 AssembleDG<basis>::local_rhs_edge_neumann(
541  const std::shared_ptr<dealii::Function<lifex::dim>> &g) const
542 {
543  AssertThrow(g != nullptr,
544  dealii::StandardExceptions::ExcMessage(
545  "Not valid pointer to the Neumann function."));
546 
547  dealii::Vector<double> cell_rhs(dofs_per_cell);
548 
549  const double face_measure = face_handler.get_measure();
550  const double unit_face_measure = (4.0 - lifex::dim) / 2;
551  const double measure_ratio = face_measure / unit_face_measure;
552 
553  for (unsigned int q = 0; q < n_quad_points_face; ++q)
554  {
555  for (unsigned int i = 0; i < dofs_per_cell; ++i)
556  {
557  cell_rhs(i) +=
558  basis_ptr->shape_value(i, face_handler.quadrature_ref(q)) *
559  g->value(face_handler.quadrature_real(q)) *
560  face_handler.quadrature_weight(q) * measure_ratio;
561  }
562  }
563 
564  return cell_rhs;
565 }
566 
567 template <class basis>
568 dealii::FullMatrix<double>
569 AssembleDG<basis>::local_SC(const double stability_coefficient) const
570 {
571  const double face_measure = face_handler.get_measure();
572  const double unit_face_measure = (4.0 - lifex::dim) / 2;
573  const double measure_ratio = face_measure / unit_face_measure;
574  const double h_local = (cell->measure()) / face_measure;
575 
576  const double local_pen_coeff =
577  (stability_coefficient * pow(poly_degree, 2)) / h_local;
578 
579  dealii::FullMatrix<double> SC(dofs_per_cell, dofs_per_cell);
580 
581  for (unsigned int q = 0; q < n_quad_points_face; ++q)
582  {
583  for (unsigned int i = 0; i < dofs_per_cell; ++i)
584  {
585  for (unsigned int j = 0; j < dofs_per_cell; ++j)
586  {
587  SC(i, j) +=
588  local_pen_coeff *
589  basis_ptr->shape_value(i, face_handler.quadrature_ref(q)) *
590  basis_ptr->shape_value(j, face_handler.quadrature_ref(q)) *
591  face_handler.quadrature_weight(q) * measure_ratio;
592  }
593  }
594  }
595 
596  return SC;
597 }
598 
599 template <class basis>
600 dealii::FullMatrix<double>
601 AssembleDG<basis>::local_SC(const double stability_coefficient,
602  const dealii::Tensor<2, lifex::dim> &Sigma) const
603 {
604  const double face_measure = face_handler.get_measure();
605  const double unit_face_measure = (4.0 - lifex::dim) / 2;
606  const double measure_ratio = face_measure / unit_face_measure;
607  const double h_local = (cell->measure()) / face_measure;
608 
609  const double local_pen_coeff =
610  (stability_coefficient * pow(poly_degree, 2)) / h_local;
611 
612  const double penalty_sigma =
613  std::abs(face_handler.get_normal() * Sigma * face_handler.get_normal());
614 
615  dealii::FullMatrix<double> SC(dofs_per_cell, dofs_per_cell);
616 
617  for (unsigned int q = 0; q < n_quad_points_face; ++q)
618  {
619  for (unsigned int i = 0; i < dofs_per_cell; ++i)
620  {
621  for (unsigned int j = 0; j < dofs_per_cell; ++j)
622  {
623  SC(i, j) +=
624  local_pen_coeff * penalty_sigma *
625  basis_ptr->shape_value(i, face_handler.quadrature_ref(q)) *
626  basis_ptr->shape_value(j, face_handler.quadrature_ref(q)) *
627  face_handler.quadrature_weight(q) * measure_ratio;
628  }
629  }
630  }
631 
632  return SC;
633 }
634 
635 
636 template <class basis>
637 dealii::FullMatrix<double>
638 AssembleDG<basis>::local_SN(const double stability_coefficient) const
639 {
640  const double face_measure = face_handler.get_measure();
641  const double unit_face_measure = (4.0 - lifex::dim) / 2;
642  const double measure_ratio = face_measure / unit_face_measure;
643  const double h_local = (cell->measure()) / face_measure;
644 
645  const double local_pen_coeff =
646  (stability_coefficient * pow(poly_degree, 2)) / h_local;
647 
648  dealii::FullMatrix<double> SN(dofs_per_cell, dofs_per_cell);
649 
650  for (unsigned int q = 0; q < n_quad_points_face; ++q)
651  {
652  for (unsigned int i = 0; i < dofs_per_cell; ++i)
653  {
654  for (unsigned int j = 0; j < dofs_per_cell; ++j)
655  {
656  const unsigned int nq =
657  face_handler.corresponding_neigh_index(q, face_handler_neigh);
658 
659  SN(i, j) -=
660  local_pen_coeff *
661  basis_ptr->shape_value(i, face_handler.quadrature_ref(q)) *
662  basis_ptr->shape_value(j,
663  face_handler_neigh.quadrature_ref(nq)) *
664  face_handler.quadrature_weight(q) * measure_ratio;
665  }
666  }
667  }
668 
669  return SN;
670 }
671 
672 template <class basis>
673 dealii::FullMatrix<double>
674 AssembleDG<basis>::local_SN(const double stability_coefficient,
675  const dealii::Tensor<2, lifex::dim> &Sigma) const
676 {
677  const double face_measure = face_handler.get_measure();
678  const double unit_face_measure = (4.0 - lifex::dim) / 2;
679  const double measure_ratio = face_measure / unit_face_measure;
680  const double h_local = (cell->measure()) / face_measure;
681 
682  const double local_pen_coeff =
683  (stability_coefficient * pow(poly_degree, 2)) / h_local;
684 
685  const double penalty_sigma =
686  std::abs(face_handler.get_normal() * Sigma * face_handler.get_normal());
687 
688  dealii::FullMatrix<double> SN(dofs_per_cell, dofs_per_cell);
689 
690  for (unsigned int q = 0; q < n_quad_points_face; ++q)
691  {
692  for (unsigned int i = 0; i < dofs_per_cell; ++i)
693  {
694  for (unsigned int j = 0; j < dofs_per_cell; ++j)
695  {
696  const unsigned int nq =
697  face_handler.corresponding_neigh_index(q, face_handler_neigh);
698 
699  SN(i, j) -=
700  local_pen_coeff * penalty_sigma *
701  basis_ptr->shape_value(i, face_handler.quadrature_ref(q)) *
702  basis_ptr->shape_value(j,
703  face_handler_neigh.quadrature_ref(nq)) *
704  face_handler.quadrature_weight(q) * measure_ratio;
705  }
706  }
707  }
708 
709  return SN;
710 }
711 
712 template <class basis>
713 std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
714 AssembleDG<basis>::local_IC(const double theta) const
715 {
716  AssertThrow(theta == 1. || theta == 0. || theta == -1.,
717  dealii::StandardExceptions::ExcMessage(
718  "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
719  "or -1 (NIP method)."));
720 
721  const dealii::Tensor<2, lifex::dim> BJinv =
722  face_handler.get_jacobian_inverse();
723 
724  const double face_measure = face_handler.get_measure();
725  const double unit_face_measure = (4.0 - lifex::dim) / 2;
726  const double measure_ratio = face_measure / unit_face_measure;
727 
728  dealii::FullMatrix<double> IC(dofs_per_cell, dofs_per_cell);
729  dealii::FullMatrix<double> IC_t(dofs_per_cell, dofs_per_cell);
730 
731  for (unsigned int q = 0; q < n_quad_points_face; ++q)
732  {
733  for (unsigned int i = 0; i < dofs_per_cell; ++i)
734  {
735  for (unsigned int j = 0; j < dofs_per_cell; ++j)
736  {
737  IC(i, j) +=
738  0.5 *
739  ((basis_ptr->shape_grad(i, face_handler.quadrature_ref(q)) *
740  BJinv) *
741  face_handler.get_normal()) *
742  basis_ptr->shape_value(j, face_handler.quadrature_ref(q)) *
743  face_handler.quadrature_weight(q) * measure_ratio;
744  }
745  }
746  }
747 
748  IC *= (-1);
749  IC_t.copy_transposed(IC);
750  IC *= theta;
751 
752  return {IC, IC_t};
753 }
754 
755 template <class basis>
756 std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
757 AssembleDG<basis>::local_IC(const double theta,
758  const dealii::Tensor<2, lifex::dim> &Sigma) const
759 {
760  AssertThrow(theta == 1. || theta == 0. || theta == -1.,
761  dealii::StandardExceptions::ExcMessage(
762  "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
763  "or -1 (NIP method)."));
764 
765  const dealii::Tensor<2, lifex::dim> BJinv =
766  face_handler.get_jacobian_inverse();
767 
768  const double face_measure = face_handler.get_measure();
769  const double unit_face_measure = (4.0 - lifex::dim) / 2;
770  const double measure_ratio = face_measure / unit_face_measure;
771 
772  const dealii::Tensor<2, lifex::dim> Sigma_tr = transpose(Sigma);
773 
774  dealii::FullMatrix<double> IC(dofs_per_cell, dofs_per_cell);
775  dealii::FullMatrix<double> IC_t(dofs_per_cell, dofs_per_cell);
776 
777  for (unsigned int q = 0; q < n_quad_points_face; ++q)
778  {
779  for (unsigned int i = 0; i < dofs_per_cell; ++i)
780  {
781  for (unsigned int j = 0; j < dofs_per_cell; ++j)
782  {
783  IC(i, j) +=
784  0.5 *
785  ((basis_ptr->shape_grad(i, face_handler.quadrature_ref(q)) *
786  BJinv * Sigma_tr) *
787  face_handler.get_normal()) *
788  basis_ptr->shape_value(j, face_handler.quadrature_ref(q)) *
789  face_handler.quadrature_weight(q) * measure_ratio;
790  }
791  }
792  }
793 
794  IC *= (-1);
795  IC_t.copy_transposed(IC);
796  IC *= theta;
797 
798  return {IC, IC_t};
799 }
800 
801 template <class basis>
802 std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
803 AssembleDG<basis>::local_IB(const double theta) const
804 {
805  AssertThrow(theta == 1. || theta == 0. || theta == -1.,
806  dealii::StandardExceptions::ExcMessage(
807  "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
808  "or -1 (NIP method)."));
809 
810  const dealii::Tensor<2, lifex::dim> BJinv =
811  face_handler.get_jacobian_inverse();
812 
813  const double face_measure = face_handler.get_measure();
814  const double unit_face_measure = (4.0 - lifex::dim) / 2;
815  const double measure_ratio = face_measure / unit_face_measure;
816 
817  dealii::FullMatrix<double> IB(dofs_per_cell, dofs_per_cell);
818  dealii::FullMatrix<double> IB_t(dofs_per_cell, dofs_per_cell);
819 
820  for (unsigned int q = 0; q < n_quad_points_face; ++q)
821  {
822  for (unsigned int i = 0; i < dofs_per_cell; ++i)
823  {
824  for (unsigned int j = 0; j < dofs_per_cell; ++j)
825  {
826  IB(i, j) +=
827  ((basis_ptr->shape_grad(i, face_handler.quadrature_ref(q)) *
828  BJinv) *
829  face_handler.get_normal()) *
830  basis_ptr->shape_value(j, face_handler.quadrature_ref(q)) *
831  face_handler.quadrature_weight(q) * measure_ratio;
832  }
833  }
834  }
835 
836  IB *= (-1);
837  IB_t.copy_transposed(IB);
838  IB *= theta;
839 
840  return {IB, IB_t};
841 }
842 
843 template <class basis>
844 std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
845 AssembleDG<basis>::local_IB(const double theta,
846  const dealii::Tensor<2, lifex::dim> &Sigma) const
847 {
848  AssertThrow(theta == 1. || theta == 0. || theta == -1.,
849  dealii::StandardExceptions::ExcMessage(
850  "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
851  "or -1 (NIP method)."));
852 
853  const dealii::Tensor<2, lifex::dim> BJinv =
854  face_handler.get_jacobian_inverse();
855 
856  const double face_measure = face_handler.get_measure();
857  const double unit_face_measure = (4.0 - lifex::dim) / 2;
858  const double measure_ratio = face_measure / unit_face_measure;
859 
860  const dealii::Tensor<2, lifex::dim> Sigma_tr = transpose(Sigma);
861 
862  dealii::FullMatrix<double> IB(dofs_per_cell, dofs_per_cell);
863  dealii::FullMatrix<double> IB_t(dofs_per_cell, dofs_per_cell);
864 
865  for (unsigned int q = 0; q < n_quad_points_face; ++q)
866  {
867  for (unsigned int i = 0; i < dofs_per_cell; ++i)
868  {
869  for (unsigned int j = 0; j < dofs_per_cell; ++j)
870  {
871  IB(i, j) +=
872  ((basis_ptr->shape_grad(i, face_handler.quadrature_ref(q)) *
873  BJinv * Sigma_tr) *
874  face_handler.get_normal()) *
875  basis_ptr->shape_value(j, face_handler.quadrature_ref(q)) *
876  face_handler.quadrature_weight(q) * measure_ratio;
877  }
878  }
879  }
880 
881  IB *= (-1);
882  IB_t.copy_transposed(IB);
883  IB *= theta;
884 
885  return {IB, IB_t};
886 }
887 
888 template <class basis>
889 std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
890 AssembleDG<basis>::local_IN(const double theta) const
891 {
892  AssertThrow(theta == 1. || theta == 0. || theta == -1.,
893  dealii::StandardExceptions::ExcMessage(
894  "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
895  "or -1 (NIP method)."));
896 
897  const dealii::Tensor<2, lifex::dim> BJinv =
898  face_handler.get_jacobian_inverse();
899 
900  const double face_measure = face_handler.get_measure();
901  const double unit_face_measure = (4.0 - lifex::dim) / 2;
902  const double measure_ratio = face_measure / unit_face_measure;
903 
904  dealii::FullMatrix<double> IN(dofs_per_cell, dofs_per_cell);
905  dealii::FullMatrix<double> IN_t(dofs_per_cell, dofs_per_cell);
906 
907  for (unsigned int q = 0; q < n_quad_points_face; ++q)
908  {
909  for (unsigned int i = 0; i < dofs_per_cell; ++i)
910  {
911  for (unsigned int j = 0; j < dofs_per_cell; ++j)
912  {
913  const unsigned int nq =
914  face_handler.corresponding_neigh_index(q, face_handler_neigh);
915 
916  IN(i, j) +=
917  0.5 *
918  ((basis_ptr->shape_grad(i, face_handler.quadrature_ref(q)) *
919  BJinv) *
920  face_handler.get_normal()) *
921  basis_ptr->shape_value(j,
922  face_handler_neigh.quadrature_ref(nq)) *
923  face_handler.quadrature_weight(q) * measure_ratio;
924  }
925  }
926  }
927 
928  IN_t.copy_transposed(IN);
929  IN *= theta;
930 
931  return {IN, IN_t};
932 }
933 
934 template <class basis>
935 std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
936 AssembleDG<basis>::local_IN(const double theta,
937  const dealii::Tensor<2, lifex::dim> &Sigma) const
938 {
939  AssertThrow(theta == 1. || theta == 0. || theta == -1.,
940  dealii::StandardExceptions::ExcMessage(
941  "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
942  "or -1 (NIP method)."));
943 
944  const dealii::Tensor<2, lifex::dim> BJinv =
945  face_handler.get_jacobian_inverse();
946 
947  const double face_measure = face_handler.get_measure();
948  const double unit_face_measure = (4.0 - lifex::dim) / 2;
949  const double measure_ratio = face_measure / unit_face_measure;
950 
951  const dealii::Tensor<2, lifex::dim> Sigma_tr = transpose(Sigma);
952 
953  dealii::FullMatrix<double> IN(dofs_per_cell, dofs_per_cell);
954  dealii::FullMatrix<double> IN_t(dofs_per_cell, dofs_per_cell);
955 
956  for (unsigned int q = 0; q < n_quad_points_face; ++q)
957  {
958  for (unsigned int i = 0; i < dofs_per_cell; ++i)
959  {
960  for (unsigned int j = 0; j < dofs_per_cell; ++j)
961  {
962  const unsigned int nq =
963  face_handler.corresponding_neigh_index(q, face_handler_neigh);
964 
965  IN(i, j) +=
966  0.5 *
967  ((basis_ptr->shape_grad(i, face_handler.quadrature_ref(q)) *
968  BJinv * Sigma_tr) *
969  face_handler.get_normal()) *
970  basis_ptr->shape_value(j,
971  face_handler_neigh.quadrature_ref(nq)) *
972  face_handler.quadrature_weight(q) * measure_ratio;
973  }
974  }
975  }
976 
977  IN_t.copy_transposed(IN);
978  IN *= theta;
979 
980  return {IN, IN_t};
981 }
982 
983 template <class basis>
984 dealii::Vector<double>
985 AssembleDG<basis>::local_u0_M_rhs(
986  const lifex::LinAlg::MPI::Vector &u0,
987  const std::vector<dealii::types::global_dof_index> &dof_indices) const
988 {
989  std::vector<double> u_bdf_loc(n_quad_points);
990  std::fill(u_bdf_loc.begin(), u_bdf_loc.end(), 0);
991 
992  for (unsigned int q = 0; q < n_quad_points; ++q)
993  {
994  for (unsigned int i = 0; i < dofs_per_cell; ++i)
995  {
996  u_bdf_loc[q] +=
997  u0[dof_indices[i]] *
998  basis_ptr->shape_value(i, vol_handler.quadrature_ref(q));
999  }
1000  }
1001 
1002  dealii::Vector<double> cell_rhs(dofs_per_cell);
1003  const dealii::Tensor<2, lifex::dim> BJinv =
1004  vol_handler.get_jacobian_inverse();
1005  const double det = 1 / determinant(BJinv);
1006 
1007  for (unsigned int q = 0; q < n_quad_points; ++q)
1008  {
1009  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1010  {
1011  cell_rhs(i) +=
1012  u_bdf_loc[q] *
1013  basis_ptr->shape_value(i, vol_handler.quadrature_ref(q)) *
1014  vol_handler.quadrature_weight(q) * det;
1015  }
1016  }
1017 
1018  return cell_rhs;
1019 }
1020 
1021 template <class basis>
1022 dealii::Vector<double>
1023 AssembleDG<basis>::local_w0_M_rhs(
1024  const lifex::LinAlg::MPI::Vector &u0,
1025  const std::vector<dealii::types::global_dof_index> &dof_indices) const
1026 {
1027  dealii::Vector<double> cell_rhs(dofs_per_cell);
1028  const dealii::Tensor<2, lifex::dim> BJinv =
1029  vol_handler.get_jacobian_inverse();
1030  const double det = 1 / determinant(BJinv);
1031 
1032  for (unsigned int q = 0; q < n_quad_points; ++q)
1033  {
1034  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1035  {
1036  for (unsigned int j = 0; j < dofs_per_cell; ++j)
1037  {
1038  cell_rhs(i) +=
1039  (basis_ptr->shape_value(i, vol_handler.quadrature_ref(q))) *
1040  (basis_ptr->shape_value(j, vol_handler.quadrature_ref(q))) *
1041  vol_handler.quadrature_weight(q) * det * u0[dof_indices[j]];
1042  }
1043  }
1044  }
1045 
1046  return cell_rhs;
1047 }
1048 
1049 template <class basis>
1050 dealii::FullMatrix<double>
1051 AssembleDG<basis>::local_non_linear_fitzhugh(
1052  const lifex::LinAlg::MPI::Vector &u0,
1053  const double a,
1054  const std::vector<dealii::types::global_dof_index> &dof_indices) const
1055 {
1056  dealii::FullMatrix<double> C(dofs_per_cell, dofs_per_cell);
1057  const dealii::Tensor<2, lifex::dim> BJinv =
1058  vol_handler.get_jacobian_inverse();
1059  const double det = 1 / determinant(BJinv);
1060 
1061 
1062  double non_lu;
1063  double non_lin;
1064 
1065  for (unsigned int q = 0; q < n_quad_points; ++q)
1066  {
1067  non_lu = 0.0;
1068 
1069  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1070  {
1071  non_lu += u0[dof_indices[i]] *
1072  basis_ptr->shape_value(i, vol_handler.quadrature_ref(q));
1073  }
1074 
1075  non_lin = (non_lu - 1) * (non_lu - a);
1076 
1077  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1078  {
1079  for (unsigned int j = 0; j < dofs_per_cell; ++j)
1080  {
1081  C(i, j) +=
1082  non_lin *
1083  basis_ptr->shape_value(i, vol_handler.quadrature_ref(q)) *
1084  basis_ptr->shape_value(j, vol_handler.quadrature_ref(q)) *
1085  vol_handler.quadrature_weight(q) * det;
1086  }
1087  }
1088  }
1089 
1090  return C;
1091 }
1092 
1093 #endif /* ASSEMBLE_DG_HPP_*/
DUB_FEM_handler.hpp
AssembleDG
Class for the assembly of the main local matrices for discontinuous Galerkin methods.
Definition: assemble_DG.hpp:69
AssembleDG::face_handler
FaceHandlerDG< lifex::dim > face_handler
Face element handler.
Definition: assemble_DG.hpp:78
AssembleDG::local_M
dealii::FullMatrix< double > local_M() const
Assembly of local mass matrix:
Definition: assemble_DG.hpp:435
AssembleDG::local_SN
dealii::FullMatrix< double > local_SN(const double stability_coefficient) const
Assembly of the component of the local matrix S that is evaluated on the two sides of the edge:
Definition: assemble_DG.hpp:638
AssembleDG::local_non_linear_fitzhugh
dealii::FullMatrix< double > local_non_linear_fitzhugh(const lifex::LinAlg::MPI::Vector &u0, const double a, const std::vector< dealii::types::global_dof_index > &dof_indices) const
Assembly of the non-linear local matrix of the Fitzhugh-Nagumo model:
Definition: assemble_DG.hpp:1051
AssembleDG::local_V
dealii::FullMatrix< double > local_V() const
Assembly of stiffness matrix:
Definition: assemble_DG.hpp:376
AssembleDG::local_SC
dealii::FullMatrix< double > local_SC(const double stability_coefficient) const
Assembly of the component of the local matrix S that is evaluated on the same side of the edge:
Definition: assemble_DG.hpp:569
AssembleDG::local_IN
std::pair< dealii::FullMatrix< double >, dealii::FullMatrix< double > > local_IN(const double theta) const
Assembly of the component of the local matrix I that is evaluated on the two sides of the edge.
Definition: assemble_DG.hpp:890
AssembleDG::local_rhs
dealii::Vector< double > local_rhs(const std::shared_ptr< dealii::Function< lifex::dim >> &f_ex) const
Assembly of the local forcing term:
Definition: assemble_DG.hpp:461
AssembleDG::n_quad_points_1D
const unsigned int n_quad_points_1D
Number of quadrature points in 1D.
Definition: assemble_DG.hpp:87
AssembleDG::local_rhs_edge_dirichlet
dealii::Vector< double > local_rhs_edge_dirichlet(const double stability_coefficient, const double theta, const std::shared_ptr< lifex::utils::FunctionDirichlet > &u_ex) const
Assembly of the local matrix associated to non-homogeneous Dirichlet boundary conditions:
Definition: assemble_DG.hpp:489
AssembleDG::face_handler_neigh
FaceHandlerDG< lifex::dim > face_handler_neigh
Face element handler of the neighbor face.
Definition: assemble_DG.hpp:81
AssembleDG::~AssembleDG
virtual ~AssembleDG()=default
Destructor.
AssembleDG::dofs_per_cell
unsigned int dofs_per_cell
Degree of freedom for each cell.
Definition: assemble_DG.hpp:84
AssembleDG::local_IC
std::pair< dealii::FullMatrix< double >, dealii::FullMatrix< double > > local_IC(const double theta) const
Assembly of the component of the local matrix I that is evaluated on the same side of the edge.
Definition: assemble_DG.hpp:714
AssembleDG::basis_ptr
const std::unique_ptr< basis > basis_ptr
Basis function class.
Definition: assemble_DG.hpp:72
AssembleDG::poly_degree
const unsigned int poly_degree
Polynomial degree.
Definition: assemble_DG.hpp:96
AssembleDG::local_u0_M_rhs
dealii::Vector< double > local_u0_M_rhs(const lifex::LinAlg::MPI::Vector &u0, const std::vector< dealii::types::global_dof_index > &dof_indices) const
Assembly of the right hand side term associated to the previous time-step solution in time-dependent ...
Definition: assemble_DG.hpp:985
AssembleDG::local_rhs_edge_neumann
dealii::Vector< double > local_rhs_edge_neumann(const std::shared_ptr< dealii::Function< lifex::dim >> &g) const
Assembly of the right hand side term associated to non-homogeneous Neumann boundary conditions:
Definition: assemble_DG.hpp:540
AssembleDG::cell
dealii::DoFHandler< lifex::dim >::active_cell_iterator cell
Cell object.
Definition: assemble_DG.hpp:99
AssembleDG::local_w0_M_rhs
dealii::Vector< double > local_w0_M_rhs(const lifex::LinAlg::MPI::Vector &u0, const std::vector< dealii::types::global_dof_index > &dof_indices) const
Assembly of the right hand side term associated to the previous time-step gating variable solution (f...
Definition: assemble_DG.hpp:1023
AssembleDG::n_quad_points_face
const unsigned int n_quad_points_face
Number of quadrature points in the face: (n_quad_points_1D)^(dim-1).
Definition: assemble_DG.hpp:93
AssembleDG::vol_handler
VolumeHandlerDG< lifex::dim > vol_handler
Volume element handler.
Definition: assemble_DG.hpp:75
AssembleDG::get_dofs_per_cell
unsigned int get_dofs_per_cell() const
Return the number of degrees of freedom per element.
Definition: assemble_DG.hpp:355
AssembleDG::reinit
void reinit(const typename dealii::DoFHandler< lifex::dim >::active_cell_iterator &new_cell, const unsigned int new_edge)
Reinitialize object on the current new_edge of the new_cell.
Definition: assemble_DG.hpp:328
AssembleDG::local_IB
std::pair< dealii::FullMatrix< double >, dealii::FullMatrix< double > > local_IB(const double theta) const
Assembly of the component of the local matrix I that is evaluated on the boundary edges.
Definition: assemble_DG.hpp:803
AssembleDG::n_quad_points
const unsigned int n_quad_points
Number of quadrature points in the volume: (n_quad_points_1D)^(dim).
Definition: assemble_DG.hpp:90
active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler_DG.hpp:49
face_handler_DG.hpp
volume_handler_DG.hpp