DUBeat 1.0.1
High-order discontinuous Galerkin methods and applications to cardiac electrophysiology
Loading...
Searching...
No Matches
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
67template <class basis>
68class AssembleDG
69{
70protected:
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
101public:
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
326template <class basis>
327void
328AssembleDG<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
344template <class basis>
345void
346AssembleDG<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
353template <class basis>
354unsigned int
355AssembleDG<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
374template <class basis>
375dealii::FullMatrix<double>
376AssembleDG<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
403template <class basis>
404dealii::FullMatrix<double>
405AssembleDG<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
433template <class basis>
434dealii::FullMatrix<double>
435AssembleDG<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
459template <class basis>
460dealii::Vector<double>
461AssembleDG<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
487template <class basis>
488dealii::Vector<double>
489AssembleDG<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
538template <class basis>
539dealii::Vector<double>
540AssembleDG<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
567template <class basis>
568dealii::FullMatrix<double>
569AssembleDG<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
599template <class basis>
600dealii::FullMatrix<double>
601AssembleDG<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
636template <class basis>
637dealii::FullMatrix<double>
638AssembleDG<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
672template <class basis>
673dealii::FullMatrix<double>
674AssembleDG<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
712template <class basis>
713std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
714AssembleDG<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
755template <class basis>
756std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
757AssembleDG<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
801template <class basis>
802std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
803AssembleDG<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
843template <class basis>
844std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
845AssembleDG<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
888template <class basis>
889std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
890AssembleDG<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
934template <class basis>
935std::pair<dealii::FullMatrix<double>, dealii::FullMatrix<double>>
936AssembleDG<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
983template <class basis>
984dealii::Vector<double>
985AssembleDG<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
1021template <class basis>
1022dealii::Vector<double>
1023AssembleDG<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
1049template <class basis>
1050dealii::FullMatrix<double>
1051AssembleDG<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_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::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_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::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::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::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
FaceHandlerDG
Class for the main operations on the faces of a discontinuous Galerkin element.
Definition face_handler_DG.hpp:49
VolumeHandlerDG
Class for the main operations on a discontinuous Galerkin volume element.
Definition volume_handler_DG.hpp:51
face_handler_DG.hpp
volume_handler_DG.hpp