DUBeat 1.0.1
High-order discontinuous Galerkin methods and applications to cardiac electrophysiology
Loading...
Searching...
No Matches
model_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 MODEL_DG_HPP_
28#define MODEL_DG_HPP_
29
30#include <deal.II/base/exceptions.h>
31#include <deal.II/base/parameter_handler.h>
32
33#include <deal.II/fe/mapping_q1_eulerian.h>
34
35#include <deal.II/lac/full_matrix.h>
36
37#include <memory>
38#include <string>
39#include <vector>
40
41#include "DUB_FEM_handler.hpp"
42#include "assemble_DG.hpp"
43#include "compute_errors_DG.hpp"
44#include "dof_handler_DG.hpp"
45#include "face_handler_DG.hpp"
46#include "source/core_model.hpp"
47#include "source/geometry/mesh_handler.hpp"
48#include "source/init.hpp"
49#include "source/io/data_writer.hpp"
50#include "source/numerics/bc_handler.hpp"
51#include "source/numerics/linear_solver_handler.hpp"
52#include "source/numerics/preconditioner_handler.hpp"
53#include "source/numerics/tools.hpp"
54#include "volume_handler_DG.hpp"
55
61template <class basis>
62class ModelDG : public lifex::CoreModel
63{
64public:
66 ModelDG(std::string model_name)
67 : CoreModel(model_name)
70 std::make_shared<lifex::utils::MeshHandler>(prm_subsection_path,
71 mpi_comm))
72 , linear_solver(prm_subsection_path + " / Linear solver",
73 {"CG", "GMRES", "BiCGStab"},
74 "GMRES")
75 , preconditioner(prm_subsection_path + " / Preconditioner", true)
76 {}
77
80
82 ModelDG<basis>(const ModelDG<basis> &ModelDG) = default;
83
86
88 virtual void
89 declare_parameters(lifex::ParamHandler &params) const override;
90
92 virtual void
93 parse_parameters(lifex::ParamHandler &params) override;
94
96 unsigned int
97 get_dofs_per_cell() const;
98
100 virtual void
101 run() override;
102
104 virtual ~ModelDG() = default;
105
106protected:
108 virtual void
109 setup_system();
110
114 void
116 dealii::DynamicSparsityPattern &sparsity,
117 const dealii::AffineConstraints<double> &constraints =
118 dealii::AffineConstraints<double>(),
119 const bool keep_constrained_dofs = true,
120 const dealii::types::subdomain_id subdomain_id =
121 dealii::numbers::invalid_subdomain_id);
122
124 void
125 initialize_solution(lifex::LinAlg::MPI::Vector &solution_owned,
126 lifex::LinAlg::MPI::Vector &solution);
127
130 virtual void
132
134 void
135 create_mesh();
136
138 void
139 create_mesh(std::string mesh_path);
140
142 void
143 solve_system();
144
148 void
149 compute_errors(const lifex::LinAlg::MPI::Vector &solution_owned,
150 const lifex::LinAlg::MPI::Vector &solution_ex_owned,
151 const std::shared_ptr<dealii::Function<lifex::dim>> &u_ex,
152 const std::shared_ptr<dealii::Function<lifex::dim>> &grad_u_ex,
153 const char *solution_name = (char *)"u") const;
154
160 virtual void
161 output_results(std::string output_name = "solution") const;
162
167 virtual void
168 conversion_to_fem(lifex::LinAlg::MPI::Vector &sol_owned);
169
172 virtual lifex::LinAlg::MPI::Vector
173 conversion_to_fem(const lifex::LinAlg::MPI::Vector &sol_owned) const;
174
175 // As conversion_to_fem but used when converting to more refined FEM meshes.
176 // See the description of the equivalent method in DUBFEMHandler. Notice that
177 // it works only when using Dubiner basis and it requires that the Dubiner
178 // mesh is nested in the FEM mesh.
179 virtual void
180 conversion_to_fem(lifex::LinAlg::MPI::Vector &sol_owned,
181 const std::string fem_file_path,
182 const unsigned int degree_fem = 1,
183 const double scaling_factor = 1) const;
184
187 virtual void
188 conversion_to_dub(lifex::LinAlg::MPI::Vector &sol_owned);
189
192 virtual lifex::LinAlg::MPI::Vector
193 conversion_to_dub(const lifex::LinAlg::MPI::Vector &sol_owned) const;
194
195 // As conversion_to_dub but used when converting from more refined FEM meshes.
196 // See the description of the equivalent method in DUBFEMHandler. Notice that
197 // it works only when using Dubiner basis and it requires that the Dubiner
198 // mesh is nested in the FEM mesh.
199 virtual void
200 conversion_to_dub(lifex::LinAlg::MPI::Vector &sol_owned,
201 const std::string fem_file_path,
202 const unsigned int degree_fem = 1,
203 const double scaling_factor = 1) const;
204
206 void
208 const std::shared_ptr<dealii::Function<lifex::dim>> &u_analytical,
209 lifex::LinAlg::MPI::Vector &sol_owned);
210
212 const std::string model_name;
214 unsigned int prm_fe_degree;
216 unsigned int prm_n_refinements;
222 double scaling_factor = 1;
224 std::shared_ptr<lifex::utils::MeshHandler> triangulation;
226 unsigned int dofs_per_cell;
231 std::shared_ptr<DUBFEMHandler<basis>> dub_fem_values;
233 std::unique_ptr<AssembleDG<basis>> assemble;
235 lifex::utils::LinearSolverHandler<lifex::LinAlg::MPI::Vector> linear_solver;
237 lifex::utils::PreconditionerHandler preconditioner;
239 lifex::LinAlg::MPI::SparseMatrix matrix;
241 lifex::LinAlg::MPI::Vector rhs;
243 lifex::LinAlg::MPI::Vector solution_owned;
245 lifex::LinAlg::MPI::Vector solution;
247 lifex::LinAlg::MPI::Vector solution_ex_owned;
249 lifex::LinAlg::MPI::Vector solution_ex;
251 std::shared_ptr<lifex::utils::FunctionDirichlet> u_ex;
253 std::shared_ptr<dealii::Function<lifex::dim>> grad_u_ex;
255 std::shared_ptr<dealii::Function<lifex::dim>> f_ex;
257 std::shared_ptr<dealii::Function<lifex::dim>> g_n;
258};
259
260template <class basis>
261void
262ModelDG<basis>::declare_parameters(lifex::ParamHandler &params) const
263{
264 // Default parameters.
265 linear_solver.declare_parameters(params);
266 preconditioner.declare_parameters(params);
267
268 // Extra parameters.
269 params.enter_subsection("Mesh and space discretization");
270 {
271 params.declare_entry(
272 "Number of refinements",
273 "2",
274 dealii::Patterns::Integer(0),
275 "Number of global mesh refinement steps applied to initial grid.");
276 params.declare_entry("FE space degree",
277 "1",
278 dealii::Patterns::Integer(1),
279 "Degree of the FE space.");
280 }
281 params.leave_subsection();
282
283 params.enter_subsection("Discontinuous Galerkin");
284 {
285 params.declare_entry(
286 "Penalty coefficient",
287 "1",
288 dealii::Patterns::Double(-1, 1),
289 "Penalty coefficient in the Discontinuous Galerkin formulation.");
290 params.declare_entry(
291 "Stability coefficient",
292 "10",
293 dealii::Patterns::Double(0),
294 "Stabilization term in the Discontinuous Galerkin formulation.");
295 }
296 params.leave_subsection();
297}
298
299template <class basis>
300void
301ModelDG<basis>::parse_parameters(lifex::ParamHandler &params)
302{
303 // Parse input file.
304 params.parse();
305 // Read input parameters.
306 linear_solver.parse_parameters(params);
307 preconditioner.parse_parameters(params);
308
309 // Extra parameters.
310 params.enter_subsection("Mesh and space discretization");
311 prm_n_refinements = params.get_integer("Number of refinements");
312
313 prm_fe_degree = params.get_integer("FE space degree");
314 params.leave_subsection();
315
316 params.enter_subsection("Discontinuous Galerkin");
317 prm_penalty_coeff = params.get_double("Penalty coefficient");
318 AssertThrow(prm_penalty_coeff == 1. || prm_penalty_coeff == 0. ||
319 prm_penalty_coeff == -1.,
320 dealii::StandardExceptions::ExcMessage(
321 "Penalty coefficient must be 1 (SIP method) or 0 (IIP method) "
322 "or -1 (NIP method)."));
323
324 prm_stability_coeff = params.get_double("Stability coefficient");
325 params.leave_subsection();
326}
327
328template <class basis>
329unsigned int
331{
332 // The analytical formula is:
333 // n_dof_per_cell = (p+1)*(p+2)*...(p+d) / d!,
334 // where p is the space order and d the space dimension..
335
336 unsigned int denominator = 1;
337 unsigned int nominator = 1;
338
339 for (unsigned int i = 1; i <= lifex::dim; i++)
340 {
341 denominator *= i;
342 nominator *= prm_fe_degree + i;
343 }
344
345 return (int)(nominator / denominator);
346}
347
348
349template <class basis>
350void
352{
353 // Initialization
354 create_mesh();
355 setup_system();
356 initialize_solution(solution_owned, solution);
357 initialize_solution(solution_ex_owned, solution_ex);
358 discretize_analytical_solution(u_ex, solution_ex_owned);
359
360
361 // Initial guess.
362 solution_ex = solution_ex_owned;
363 solution = solution_owned = 0;
364
365 // Computation of the numerical solution.
366 assemble_system();
367 solve_system();
368
369
370 // Computation and output of the errors.
371 compute_errors(solution_owned, solution_ex_owned, u_ex, grad_u_ex, "u");
372
373 // Generation of the graphical output.
374 conversion_to_fem(solution_ex);
375 solution = solution_owned;
376 conversion_to_fem(solution);
377 output_results();
378}
379
380template <class basis>
381void
383{
384 assemble = std::make_unique<AssembleDG<basis>>(prm_fe_degree);
385
386 dof_handler.reinit(triangulation->get());
387 dof_handler.distribute_dofs(prm_fe_degree);
388 dub_fem_values =
389 std::make_shared<DUBFEMHandler<basis>>(prm_fe_degree, dof_handler);
390
391 triangulation->get_info().print(prm_subsection_path,
392 dof_handler.n_dofs(),
393 true);
394
395 dealii::IndexSet owned_dofs = dof_handler.locally_owned_dofs();
396 dealii::IndexSet relevant_dofs = owned_dofs;
397 dealii::DynamicSparsityPattern dsp(relevant_dofs);
398
399
400 // Add (dof, dof_neigh) to dsp, so to the matrix
401 dofs_per_cell = get_dofs_per_cell();
402 std::vector<lifex::types::global_dof_index> dof_indices(dofs_per_cell);
403 std::vector<lifex::types::global_dof_index> dof_indices_neigh(dofs_per_cell);
404
405
406 for (const auto &cell : dof_handler.active_cell_iterators())
407 {
408 dof_indices = dof_handler.get_dof_indices(cell);
409
410 for (const auto &edge : cell->face_indices())
411 {
412 if (!cell->at_boundary(edge))
413 {
414 const auto neighcell = cell->neighbor(edge);
415 dof_indices_neigh = dof_handler.get_dof_indices(neighcell);
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 dsp.add(dof_indices[i], dof_indices_neigh[j]);
421 }
422 }
423 }
424 }
425 }
426
427 make_sparsity_pattern(dof_handler, dsp);
428 lifex::SparsityTools::distribute_sparsity_pattern(dsp,
429 owned_dofs,
430 mpi_comm,
431 relevant_dofs);
432 lifex::utils::initialize_matrix(matrix, owned_dofs, dsp);
433
434 rhs.reinit(owned_dofs, mpi_comm);
435}
436
437
438template <class basis>
439void
441 const DoFHandlerDG<basis> &dof,
442 dealii::DynamicSparsityPattern &sparsity,
443 const dealii::AffineConstraints<double> &constraints,
444 const bool keep_constrained_dofs,
445 const dealii::types::subdomain_id subdomain_id)
446{
447 Assert(sparsity.n_rows() == n_dofs,
448 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
449 Assert(sparsity.n_cols() == n_dofs,
450 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
451
452 // If we have a distributed Triangulation only allow locally_owned
453 // subdomain. Not setting a subdomain is also okay, because we skip
454 // ghost cells in the loop below.
455 if (const auto *triangulation = dynamic_cast<
456 const dealii::parallel::DistributedTriangulationBase<lifex::dim> *>(
457 &dof.get_triangulation()))
458 {
459 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
460 (subdomain_id == triangulation->locally_owned_subdomain()),
461 ExcMessage(
462 "For distributed Triangulation objects and associated "
463 "DoFHandler objects, asking for any subdomain other than the "
464 "locally owned one does not make sense."));
465 }
466 std::vector<dealii::types::global_dof_index> dofs_on_this_cell;
467 dofs_on_this_cell.reserve(dof.n_dofs_per_cell());
468
469 // In case we work with a distributed sparsity pattern of Trilinos
470 // type, we only have to do the work if the current cell is owned by
471 // the calling processor. Otherwise, just continue.
472 for (const auto &cell : dof.active_cell_iterators())
473 if (((subdomain_id == dealii::numbers::invalid_subdomain_id) ||
474 (subdomain_id == cell->subdomain_id())) &&
475 cell->is_locally_owned())
476 {
477 dofs_on_this_cell.resize(dof.n_dofs_per_cell());
478 dofs_on_this_cell = dof.get_dof_indices(cell);
479
480 // make sparsity pattern for this cell. if no constraints pattern
481 // was given, then the following call acts as if simply no
482 // constraints existed
483 constraints.add_entries_local_to_global(dofs_on_this_cell,
484 sparsity,
485 keep_constrained_dofs);
486 }
487}
488
489
490template <class basis>
491void
492ModelDG<basis>::initialize_solution(lifex::LinAlg::MPI::Vector &solution_owned,
493 lifex::LinAlg::MPI::Vector &solution)
494{
495 dealii::IndexSet owned_dofs = dof_handler.locally_owned_dofs();
496 dealii::IndexSet relevant_dofs = owned_dofs;
497
498 solution_owned.reinit(owned_dofs, mpi_comm);
499 solution.reinit(owned_dofs, relevant_dofs, mpi_comm);
500}
501
502template <class basis>
503void
505{
506 std::string mesh_path = "../meshes/" + std::to_string(lifex::dim) + "D_" +
507 std::to_string(prm_n_refinements) + ".msh";
508 AssertThrow(std::filesystem::exists(mesh_path),
509 dealii::StandardExceptions::ExcMessage(
510 "This mesh file/directory does not exist."));
511
512 // deal.II does not provide runtime generation of tetrahedral meshes, hence
513 // they can currently be imported only from file. This version of create_mesh
514 // picks the mesh file from the default path.
515 triangulation->initialize_from_file(mesh_path, this->scaling_factor);
516 triangulation->set_element_type(lifex::utils::MeshHandler::ElementType::Tet);
517 triangulation->create_mesh();
518}
519
520template <class basis>
521void
522ModelDG<basis>::create_mesh(std::string mesh_path)
523{
524 AssertThrow(std::filesystem::exists(mesh_path),
525 dealii::StandardExceptions::ExcMessage(
526 "This mesh file/directory does not exist."));
527
528 // deal.II does not provide runtime generation of tetrahedral meshes, hence
529 // they can currently be imported only from file. This version of create_mesh
530 // picks the mesh file from a user-defined path.
531 triangulation->initialize_from_file(mesh_path, this->scaling_factor);
532 triangulation->set_element_type(lifex::utils::MeshHandler::ElementType::Tet);
533 triangulation->create_mesh();
534}
535
536template <class basis>
537void
539{
540 preconditioner.initialize(matrix);
541 linear_solver.solve(matrix, solution_owned, rhs, preconditioner);
542 solution = solution_owned;
543}
544
545template <class basis>
546void
548 const lifex::LinAlg::MPI::Vector &solution_owned,
549 const lifex::LinAlg::MPI::Vector &solution_ex_owned,
550 const std::shared_ptr<dealii::Function<lifex::dim>> &u_ex,
551 const std::shared_ptr<dealii::Function<lifex::dim>> &grad_u_ex,
552 const char *solution_name) const
553{
554 ComputeErrorsDG<basis> error_calculator(prm_fe_degree,
555 prm_stability_coeff,
556 dofs_per_cell,
557 dof_handler,
558 prm_n_refinements,
559 model_name);
560 error_calculator.reinit(
561 solution_owned, solution_ex_owned, u_ex, grad_u_ex, solution_name);
562 error_calculator.compute_errors();
563 std::vector<double> errors = error_calculator.output_errors();
564
565 pcout << std::endl
566 << solution_name << ":" << std::endl
567 << "Error L^inf: " << std::setw(6) << std::fixed << std::setprecision(6)
568 << errors[0] << std::endl
569 << "Error L^2: " << std::setw(6) << std::fixed << std::setprecision(6)
570 << errors[1] << std::endl
571 << "Error H^1: " << std::setw(6) << std::fixed << std::setprecision(6)
572 << errors[2] << std::endl
573 << "Error DG: " << std::setw(6) << std::fixed << std::setprecision(6)
574 << errors[3] << std::endl;
575
576 error_calculator.update_datafile();
577}
578
579template <class basis>
580void
581ModelDG<basis>::output_results(std::string output_name) const
582{
583 lifex::DataOut<lifex::dim> data_out;
584
585 // To output results, we need the deal.II DoFHandler instead of the DUBeat
586 // DGDoFHandler since we are required here to use deal.II functions. On the
587 // other hand, the deal.II DofHandler is limited to the 2nd order polynomials.
588 const unsigned int output_fe_degree = (prm_fe_degree < 3) ? prm_fe_degree : 2;
589 dealii::FE_SimplexDGP<lifex::dim> fe(output_fe_degree);
590 dealii::DoFHandler<lifex::dim> dof_handler_fem;
591 dof_handler_fem.reinit(triangulation->get());
592 dof_handler_fem.distribute_dofs(fe);
593
594 data_out.add_data_vector(dof_handler_fem, solution, "u");
595 data_out.build_patches();
596 data_out.add_data_vector(dof_handler_fem, solution_ex, "u_ex");
597 data_out.build_patches();
598 lifex::utils::dataout_write_hdf5(data_out, output_name, false);
599
600 data_out.clear();
601}
602
605template <class basis>
606void
607ModelDG<basis>::conversion_to_fem(lifex::LinAlg::MPI::Vector &sol_owned)
608{
609 return;
610}
611
614template <>
615void
617 lifex::LinAlg::MPI::Vector &sol_owned)
618{
619 sol_owned = dub_fem_values->dubiner_to_fem(sol_owned);
620}
621
623template <class basis>
624lifex::LinAlg::MPI::Vector
626 const lifex::LinAlg::MPI::Vector &sol_owned) const
627{
628 return sol_owned;
629}
630
632template <>
633lifex::LinAlg::MPI::Vector
635 const lifex::LinAlg::MPI::Vector &sol_owned) const
636{
637 lifex::LinAlg::MPI::Vector sol_fem =
638 dub_fem_values->dubiner_to_fem(sol_owned);
639 return sol_fem;
640}
641
642template <class basis>
643void
644ModelDG<basis>::conversion_to_fem(lifex::LinAlg::MPI::Vector &sol_owned,
645 const std::string fem_file_path,
646 const unsigned int degree_fem,
647 const double scaling_factor) const
648{
649 return;
650}
651
652template <>
653void
655 lifex::LinAlg::MPI::Vector &sol_owned,
656 const std::string fem_file_path,
657 const unsigned int degree_fem,
658 const double scaling_factor) const
659{
660 sol_owned = dub_fem_values->dubiner_to_fem(sol_owned,
661 fem_file_path,
662 prm_subsection_path,
663 mpi_comm,
664 degree_fem,
665 scaling_factor);
666}
667
670template <class basis>
671void
672ModelDG<basis>::conversion_to_dub(lifex::LinAlg::MPI::Vector &sol_owned)
673{
674 return;
675}
676
679template <>
680void
682 lifex::LinAlg::MPI::Vector &sol_owned)
683{
684 sol_owned = dub_fem_values->fem_to_dubiner(sol_owned);
685}
686
688template <class basis>
689lifex::LinAlg::MPI::Vector
691 const lifex::LinAlg::MPI::Vector &sol_owned) const
692{
693 return sol_owned;
694}
695
697template <>
698lifex::LinAlg::MPI::Vector
700 const lifex::LinAlg::MPI::Vector &sol_owned) const
701{
702 lifex::LinAlg::MPI::Vector sol_dub =
703 dub_fem_values->fem_to_dubiner(sol_owned);
704 return sol_dub;
705}
706
707template <class basis>
708void
709ModelDG<basis>::conversion_to_dub(lifex::LinAlg::MPI::Vector &sol_owned,
710 const std::string fem_file_path,
711 const unsigned int degree_fem,
712 const double scaling_factor) const
713{
714 return;
715}
716
717template <>
718void
720 lifex::LinAlg::MPI::Vector &sol_owned,
721 const std::string fem_file_path,
722 const unsigned int degree_fem,
723 const double scaling_factor) const
724{
725 sol_owned = dub_fem_values->fem_to_dubiner(sol_owned,
726 fem_file_path,
727 prm_subsection_path,
728 mpi_comm,
729 degree_fem,
730 scaling_factor);
731}
732
735template <>
736void
737ModelDG<dealii::FE_SimplexDGP<lifex::dim>>::discretize_analytical_solution(
738 const std::shared_ptr<dealii::Function<lifex::dim>> &u_analytical,
739 lifex::LinAlg::MPI::Vector &sol_owned)
740{
741 dealii::VectorTools::interpolate(dof_handler, *u_analytical, sol_owned);
742}
743
746template <>
747void
748ModelDG<DUBValues<lifex::dim>>::discretize_analytical_solution(
749 const std::shared_ptr<dealii::Function<lifex::dim>> &u_analytical,
750 lifex::LinAlg::MPI::Vector &sol_owned)
751{
752 sol_owned = dub_fem_values->analytical_to_dubiner(sol_owned, u_analytical);
753}
754
755#endif /* MODEL_DG_HPP_*/
Class to compute the errors between the numerical solution (solution_owned) and the exact solution (s...
void update_datafile() const
Update the datafile with the new errors.
void reinit(const lifex::LinAlg::MPI::Vector &sol_owned, const lifex::LinAlg::MPI::Vector &sol_ex_owned, const std::shared_ptr< dealii::Function< lifex::dim > > &u_ex_input, const std::shared_ptr< dealii::Function< lifex::dim > > &grad_u_ex_input, const char *solution_name_input)
Reinitialization with new computed and exact solutions.
void compute_errors(std::list< const char * > errors_defs={"inf", "L2", "H1", "DG"})
Compute errors following the preferences in the list.
std::vector< double > output_errors(std::list< const char * > errors_defs={"inf", "L2", "H1", "DG"}) const
Output of the errors following the preferences in the list.
Class to work with global and local degrees of freedom and their mapping.
unsigned int n_dofs_per_cell() const
Return copy of n_dofs_per_cell.
std::vector< lifex::types::global_dof_index > get_dof_indices(active_cell_iterator cell) const
Returns the global dofs referred to the input cell.
Class representing the resolution of problems using discontinuous Galerkin methods.
Definition model_DG.hpp:63
const std::string model_name
Name of the class/problem.
Definition model_DG.hpp:212
lifex::LinAlg::MPI::SparseMatrix matrix
Distributed matrix of the linear system.
Definition model_DG.hpp:239
std::shared_ptr< dealii::Function< lifex::dim > > g_n
Neumann boundary conditions.
Definition model_DG.hpp:257
unsigned int prm_n_refinements
Mesh refinement level (>=1).
Definition model_DG.hpp:216
virtual void conversion_to_dub(lifex::LinAlg::MPI::Vector &sol_owned)
To convert a discretized solution in Dubiner basis (only for problems using Dubiner basis).
Definition model_DG.hpp:672
void initialize_solution(lifex::LinAlg::MPI::Vector &solution_owned, lifex::LinAlg::MPI::Vector &solution)
To inizialize the solutions using the deal.II reinit.
Definition model_DG.hpp:492
virtual void assemble_system()=0
Assembly of the linear system, pure virtual.
unsigned int prm_fe_degree
Polynomials degree.
Definition model_DG.hpp:214
double prm_stability_coeff
DG stabilty coefficient.
Definition model_DG.hpp:220
void make_sparsity_pattern(const DoFHandlerDG< basis > &dof, dealii::DynamicSparsityPattern &sparsity, const dealii::AffineConstraints< double > &constraints=dealii::AffineConstraints< double >(), const bool keep_constrained_dofs=true, const dealii::types::subdomain_id subdomain_id=dealii::numbers::invalid_subdomain_id)
Creation of the sparsity pattern to assign to the system matrix before assembling.
Definition model_DG.hpp:440
std::shared_ptr< DUBFEMHandler< basis > > dub_fem_values
Member used for conversions between analytical, nodal and modal representations of the solutions.
Definition model_DG.hpp:231
lifex::LinAlg::MPI::Vector solution
Distributed solution vector, with ghost entries.
Definition model_DG.hpp:245
std::unique_ptr< AssembleDG< basis > > assemble
Matrix assembler.
Definition model_DG.hpp:233
lifex::LinAlg::MPI::Vector solution_ex
Distributed exact solution vector, without ghost entries.
Definition model_DG.hpp:249
virtual void declare_parameters(lifex::ParamHandler &params) const override
Declare main parameters.
Definition model_DG.hpp:262
DoFHandlerDG< basis > dof_handler
DoFHandler (internal use for useful already implemented methods).
Definition model_DG.hpp:228
virtual void conversion_to_fem(lifex::LinAlg::MPI::Vector &sol_owned)
To convert a discretized solution from modal to nodal basis (does nothing if problem is already in no...
Definition model_DG.hpp:607
void create_mesh()
Load the mesh from the default path.
Definition model_DG.hpp:504
lifex::LinAlg::MPI::Vector solution_owned
Distributed solution vector, without ghost entries.
Definition model_DG.hpp:243
ModelDG(std::string model_name)
Constructor.
Definition model_DG.hpp:66
std::shared_ptr< dealii::Function< lifex::dim > > f_ex
Known forcing term.
Definition model_DG.hpp:255
virtual ~ModelDG()=default
Destructor.
lifex::utils::PreconditionerHandler preconditioner
Preconditioner handler.
Definition model_DG.hpp:237
double prm_penalty_coeff
DG Penalty coefficient.
Definition model_DG.hpp:218
std::shared_ptr< dealii::Function< lifex::dim > > grad_u_ex
Pointer to exact gradient solution Function.
Definition model_DG.hpp:253
virtual void output_results(std::string output_name="solution") const
Output of results.
Definition model_DG.hpp:581
void compute_errors(const lifex::LinAlg::MPI::Vector &solution_owned, const lifex::LinAlg::MPI::Vector &solution_ex_owned, const std::shared_ptr< dealii::Function< lifex::dim > > &u_ex, const std::shared_ptr< dealii::Function< lifex::dim > > &grad_u_ex, const char *solution_name=(char *)"u") const
To compute the error, the error, the error and the error at the end of system solving,...
Definition model_DG.hpp:547
unsigned int dofs_per_cell
Number of degrees of freedom per cell.
Definition model_DG.hpp:226
void solve_system()
System solving.
Definition model_DG.hpp:538
std::shared_ptr< lifex::utils::MeshHandler > triangulation
Triangulation (internal use for useful already implemented methods).
Definition model_DG.hpp:224
double scaling_factor
Scaling factor.
Definition model_DG.hpp:222
std::shared_ptr< lifex::utils::FunctionDirichlet > u_ex
Pointer to exact solution function.
Definition model_DG.hpp:251
lifex::LinAlg::MPI::Vector solution_ex_owned
Distributed exact solution vector, without ghost entries.
Definition model_DG.hpp:247
virtual void run() override
Run the simulation.
Definition model_DG.hpp:351
lifex::LinAlg::MPI::Vector rhs
Distributed right hand side vector of the linear system.
Definition model_DG.hpp:241
unsigned int get_dofs_per_cell() const
Return the number of degrees of freedom per element.
Definition model_DG.hpp:330
void discretize_analytical_solution(const std::shared_ptr< dealii::Function< lifex::dim > > &u_analytical, lifex::LinAlg::MPI::Vector &sol_owned)
Conversion of an analytical solution from FEM to basis coefficients.
lifex::utils::LinearSolverHandler< lifex::LinAlg::MPI::Vector > linear_solver
Linear solver handler.
Definition model_DG.hpp:235
virtual void setup_system()
Setup of the problem before the resolution.
Definition model_DG.hpp:382
virtual void parse_parameters(lifex::ParamHandler &params) override
Parse parameters from .prm file.
Definition model_DG.hpp:301