Logbook  (07-04-2025)
Static problems
static_vector_solver_ii.hpp
1 /******************************************************************************
2  * Copyright (C) Siarhei Uzunbajakau, 2023.
3  *
4  * This program is free software. You can use, modify, and redistribute it under
5  * the terms of the GNU Lesser General Public License as published by the Free
6  * Software Foundation, either version 3 or (at your option) any later version.
7  * This program is distributed without any warranty.
8  *
9  * Refer to COPYING.LESSER for more details.
10  ******************************************************************************/
11 
12 #ifndef StaticVectorSolverII_H__
13 #define StaticVectorSolverII_H__
14 
15 #include <deal.II/base/exceptions.h>
16 #include <deal.II/base/tensor.h>
17 #include <deal.II/base/thread_management.h>
18 #include <deal.II/base/timer.h>
19 #include <deal.II/base/types.h>
20 #include <deal.II/base/work_stream.h>
21 
22 #include <deal.II/grid/tria.h>
23 
24 #include <deal.II/dofs/dof_handler.h>
25 #include <deal.II/dofs/dof_renumbering.h>
26 #include <deal.II/dofs/dof_tools.h>
27 
28 #include <deal.II/fe/fe_nedelec.h>
29 #include <deal.II/fe/fe_values.h>
30 
31 #include <deal.II/lac/affine_constraints.h>
32 #include <deal.II/lac/dynamic_sparsity_pattern.h>
33 #include <deal.II/lac/precondition.h>
34 #include <deal.II/lac/solver_cg.h>
35 #include <deal.II/lac/solver_control.h>
36 #include <deal.II/lac/sparse_matrix.h>
37 #include <deal.II/lac/vector.h>
38 
39 #include <deal.II/numerics/data_out.h>
40 #include <deal.II/numerics/vector_tools.h>
41 #include <deal.II/numerics/vector_tools_common.h>
42 #include <deal.II/numerics/vector_tools_project.h>
43 
44 #include <fstream>
45 #include <iomanip>
46 #include <ios>
47 #include <iostream>
48 #include <map>
49 
50 #include "constants.hpp"
51 #include "settings.hpp"
52 #include "static_vector_input.hpp"
53 
54 #define VE scratch_data.ve
55 #define SE scratch_data.se
56 
57 #define TMR(__name) TimerOutput::Scope timer_section(timer, __name)
58 
59 using namespace dealii;
60 
61 namespace StaticVectorSolver {
222 template<int dim, int stage = 1>
223 class Solver2
224 {
225 public:
226  Solver2() = delete;
227 
265  Solver2(unsigned int p,
266  unsigned int mapping_degree,
267  const Triangulation<dim>& triangulation_T,
268  const DoFHandler<dim>& dof_handler_T,
269  const Vector<double>& solution_T,
270  unsigned int type_of_pde_rhs = 0,
271  double eta_squared = 0.0,
272  std::string fname = "data",
273  const Function<dim>* exact_solution = nullptr,
274  bool print_time_tables = false,
275  bool project_exact_solution = false,
276  bool write_higher_order_cells = false)
277  : triangulation_T(triangulation_T)
278  , dof_handler_T(dof_handler_T)
279  , solution_T(solution_T)
280  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281  // The public attribute fe.degree is the maximal polynomial degree of a
282  // shape function in a single coordinate direction, not the degree of
283  // the finite element. For FE_Nedelec and FE_RaviartThomas degree of
284  // the finite element is: degree_of_element = fe.degree - 1.
285  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
286  // fe(dof_handler_T.get_fe().degree-1),
287  , fe(p)
288  , mapping_degree(mapping_degree)
289  , type_of_pde_rhs(type_of_pde_rhs)
290  , eta_squared(eta_squared)
291  , fname(fname)
292  , exact_solution(exact_solution)
293  , print_time_tables(print_time_tables)
294  , project_exact_solution(project_exact_solution)
295  , write_higher_order_cells(write_higher_order_cells)
296  {
297  }
298 
328  virtual void fill_dirichlet_stack() = 0;
329 
333  virtual void solve() = 0;
334 
344  void setup();
345 
349  void assemble();
350 
354  void compute_error_norms();
355 
366  void project_exact_solution_fcn();
367 
400  void save() const;
401 
417  void save_matrix_and_rhs_to_csv(std::string fname) const;
418 
423  void clear()
424  {
425  system_matrix.clear();
426  system_rhs.reinit(0);
427  }
428 
432  const DoFHandler<dim>& get_dof_handler() const { return dof_handler; }
433 
437  const Vector<double>& get_solution() const { return solution; }
438 
442  unsigned int get_n_cells() const
443  {
444  return static_cast<unsigned int>(triangulation_T.n_active_cells());
445  }
446 
450  unsigned int get_n_dofs() const
451  {
452  return static_cast<unsigned int>(dof_handler.n_dofs());
453  }
454 
458  double get_L2_norm() const { return L2_norm; }
459 
463  double get_Linfty_norm() const { return L2_norm; }
464 
469  unsigned int get_mapping_degree() const { return mapping_degree; }
470 
478  void run()
479  {
480  TimerOutput::OutputFrequency tf =
481  (print_time_tables) ? TimerOutput::summary : TimerOutput::never;
482 
483  TimerOutput timer(std::cout, tf, TimerOutput::cpu_and_wall_times_grouped);
484 
485  {
486  TMR("Fill Dirichlet stack");
487  fill_dirichlet_stack();
488  }
489  {
490  TMR("Setup");
491  setup();
492  }
493  {
494  TMR("Assemble");
495  assemble();
496  }
497  {
498  TMR("Solve");
499  solve();
500  }
501 
502  if (exact_solution) {
503  if (project_exact_solution) {
504  TMR("Project exact solution");
505  project_exact_solution_fcn();
506  }
507 
508  {
509  TMR("Compute error norms");
510  compute_error_norms();
511  }
512  }
513 
514  {
515  TMR("Save");
516  save();
517  }
518  };
519 
520  virtual ~Solver2() = default;
521 
522 protected:
526  const Triangulation<dim>& triangulation_T;
527 
532  const DoFHandler<dim>& dof_handler_T;
533 
538  const Vector<double>& solution_T;
539 
544  std::map<types::boundary_id, const Function<dim>*> dirichlet_stack;
545 
549  const FE_Nedelec<dim> fe;
550 
554  DoFHandler<dim> dof_handler;
555 
560  Vector<double> solution;
561 
565  Vector<double> projected_exact_solution;
566 
570  AffineConstraints<double> constraints;
571 
575  SparsityPattern sparsity_pattern;
576 
580  SparseMatrix<double> system_matrix;
581 
585  Vector<double> system_rhs;
586 
590  double L2_norm;
591 
595  double Linfty_norm;
596 
597 private:
598  const unsigned int mapping_degree;
599  const unsigned int type_of_pde_rhs;
600  const double eta_squared;
601  const std::string fname;
602  const Function<dim>* exact_solution;
603  const bool print_time_tables;
604  const bool project_exact_solution;
605  const bool write_higher_order_cells;
606 
607  Vector<float> L2_per_cell;
608  Vector<float> Linfty_per_cell;
609 
610  // ----------------------------------------------------------------------------
611  // These structures and functions are related to the Work Stream algorithm.
612  // See article "WorkStream – A Design Pattern for Multicore-Enabled Finite
613  // Element Computations." by BRUNO TURCKSIN, MARTIN KRONBICHLER,
614  // WOLFGANG BANGERTH for more details.
615  // ----------------------------------------------------------------------------
616 
617  using IteratorTuple =
618  std::tuple<typename DoFHandler<dim>::active_cell_iterator,
619  typename DoFHandler<dim>::active_cell_iterator>;
620 
621  using IteratorPair = SynchronousIterators<IteratorTuple>;
622 
623  struct AssemblyScratchData
624  {
625  AssemblyScratchData(const FiniteElement<dim>& fe,
626  const DoFHandler<dim>& dof_hand_T,
627  const Vector<double>& dofs_T,
628  unsigned int mapping_degree,
629  unsigned int type_of_pde_rhs,
630  double eta_squared);
631 
632  AssemblyScratchData(const AssemblyScratchData& scratch_data);
633 
634  TheCoefficient<dim, stage> the_coefficient;
635  Gamma<dim, stage> gamma;
636  RobinRhs<dim, stage> robin_rhs;
637  FreeSurfaceCurrent<dim, stage> free_surface_current;
638 
639  MappingQ<dim> mapping;
640  Constants::QuadratureTableVector<dim> qt;
641  FEValues<dim> fe_values;
642  FEFaceValues<dim> fe_face_values;
643 
644  FEValues<dim> fe_values_T;
645  FEFaceValues<dim> fe_face_values_T;
646 
647  const unsigned int dofs_per_cell;
648  const unsigned int n_q_points;
649  const unsigned int n_q_points_face;
650 
651  std::vector<double> the_coefficient_list;
652  std::vector<Tensor<1, dim>> vector_values;
653  std::vector<Tensor<1, dim>> vector_values_face;
654  std::vector<double> values;
655  std::vector<double> values_face;
656  std::vector<double> gamma_list;
657  std::vector<Tensor<1, dim>> robin_rhs_list;
658  std::vector<Tensor<1, dim>> free_surface_current_list;
659 
660  const DoFHandler<dim>& dof_hand_T;
661  const Vector<double>& dofs_T;
662 
663  const unsigned int type_of_pde_rhs;
664  const double eta_squared;
665  const FEValuesExtractors::Vector ve;
666  const FEValuesExtractors::Scalar se;
667  bool do_robin;
668  bool do_K;
669  bool do_T_on_boundary;
670  };
671 
672  struct AssemblyCopyData
673  {
674  FullMatrix<double> cell_matrix;
675  Vector<double> cell_rhs;
676  std::vector<types::global_dof_index> local_dof_indices;
677  };
678 
679  void system_matrix_local(const IteratorPair& IP,
680  AssemblyScratchData& scratch_data,
681  AssemblyCopyData& copy_data);
682 
683  void copy_local_to_global(const AssemblyCopyData& copy_data);
684  //-----------------------------------------------------------------------------
685  //-----------------------------------------------------------------------------
686  //-----------------------------------------------------------------------------
687 };
688 
689 template<int dim, int stage>
690 void
691 Solver2<dim, stage>::setup()
692 {
693  dof_handler.reinit(triangulation_T);
694  dof_handler.distribute_dofs(fe);
695 
696  constraints.clear();
697  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
698 
699 #pragma GCC diagnostic push
700 #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
701  for (auto item : dirichlet_stack) {
702  Assert(item.first % 2 == 1, ExcInternalError());
703  }
704 #pragma GCC diagnostic pop
705 
706  for (auto item : dirichlet_stack)
707  VectorTools::project_boundary_values_curl_conforming_l2(
708  dof_handler,
709  0, // first vector component
710  *item.second, // boundary function
711  item.first, // boundary id
712  constraints, // constraints
713  MappingQ<dim>(mapping_degree));
714 
715  constraints.close();
716 
717  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
718  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
719 
720  sparsity_pattern.copy_from(dsp);
721  system_matrix.reinit(sparsity_pattern);
722  solution.reinit(dof_handler.n_dofs());
723  system_rhs.reinit(dof_handler.n_dofs());
724 
725  if (project_exact_solution)
726  projected_exact_solution.reinit(dof_handler.n_dofs());
727 
728  if (exact_solution) {
729  L2_per_cell.reinit(triangulation_T.n_active_cells());
730  Linfty_per_cell.reinit(triangulation_T.n_active_cells());
731  }
732 }
733 
734 template<int dim, int stage>
735 void
736 Solver2<dim, stage>::assemble()
737 {
738  WorkStream::run(
739  IteratorPair(
740  IteratorTuple(dof_handler.begin_active(), dof_handler_T.begin_active())),
741  IteratorPair(IteratorTuple(dof_handler.end(), dof_handler_T.end())),
742  *this,
743  &Solver2::system_matrix_local,
744  &Solver2::copy_local_to_global,
745  AssemblyScratchData(fe,
746  dof_handler_T,
747  solution_T,
748  mapping_degree,
749  type_of_pde_rhs,
750  eta_squared),
751  AssemblyCopyData());
752 }
753 
754 template<int dim, int stage>
755 Solver2<dim, stage>::AssemblyScratchData::AssemblyScratchData(
756  const FiniteElement<dim>& fe,
757  const DoFHandler<dim>& dof_hand_T,
758  const Vector<double>& dofs_T,
759  unsigned int mapping_degree,
760  unsigned int type_of_pde_rhs,
761  double eta_squared)
762  : mapping(mapping_degree)
763  , qt(fe.degree - 1)
764  , fe_values(mapping,
765  fe,
766  QGauss<dim>(qt.sim()),
767  update_gradients | update_values | update_quadrature_points |
768  update_JxW_values)
769  , fe_face_values(mapping,
770  fe,
771  QGauss<dim - 1>(qt.sim()),
772  update_values | update_normal_vectors |
773  update_quadrature_points | update_JxW_values)
774  , fe_values_T(mapping,
775  dof_hand_T.get_fe(),
776  QGauss<dim>(qt.sim()),
777  update_values | update_gradients)
778  , fe_face_values_T(mapping,
779  dof_hand_T.get_fe(),
780  QGauss<dim - 1>(qt.sim()),
781  update_values | update_gradients)
782  , dofs_per_cell(fe_values.dofs_per_cell)
783  , n_q_points(fe_values.get_quadrature().size())
784  , n_q_points_face(fe_face_values.get_quadrature().size())
785  , the_coefficient_list(n_q_points)
786  , vector_values(n_q_points)
787  , vector_values_face(n_q_points_face)
788  , values(n_q_points)
789  , values_face(n_q_points_face)
790  , gamma_list(n_q_points_face)
791  , robin_rhs_list(n_q_points_face)
792  , free_surface_current_list(n_q_points_face)
793  , dof_hand_T(dof_hand_T)
794  , dofs_T(dofs_T)
795  , type_of_pde_rhs(type_of_pde_rhs)
796  , eta_squared(eta_squared)
797  , ve(0)
798  , se(0)
799 {
800 }
801 
802 template<int dim, int stage>
803 Solver2<dim, stage>::AssemblyScratchData::AssemblyScratchData(
804  const AssemblyScratchData& scratch_data)
805  : mapping(scratch_data.mapping.get_degree())
806  , qt(scratch_data.qt)
807  , fe_values(mapping,
808  scratch_data.fe_values.get_fe(),
809  scratch_data.fe_values.get_quadrature(),
810  update_gradients | update_values | update_quadrature_points |
811  update_JxW_values)
812  , fe_face_values(mapping,
813  scratch_data.fe_face_values.get_fe(),
814  scratch_data.fe_face_values.get_quadrature(),
815  update_values | update_normal_vectors |
816  update_quadrature_points | update_JxW_values)
817  , fe_values_T(mapping,
818  scratch_data.fe_values_T.get_fe(),
819  scratch_data.fe_values_T.get_quadrature(),
820  update_values | update_gradients)
821  , fe_face_values_T(mapping,
822  scratch_data.fe_face_values_T.get_fe(),
823  scratch_data.fe_face_values_T.get_quadrature(),
824  update_values | update_gradients)
825  , dofs_per_cell(fe_values.dofs_per_cell)
826  , n_q_points(fe_values.get_quadrature().size())
827  , n_q_points_face(fe_face_values.get_quadrature().size())
828  , the_coefficient_list(n_q_points)
829  , vector_values(n_q_points)
830  , vector_values_face(n_q_points_face)
831  , values(n_q_points)
832  , values_face(n_q_points_face)
833  , gamma_list(n_q_points_face)
834  , robin_rhs_list(n_q_points_face)
835  , free_surface_current_list(n_q_points_face)
836  , dof_hand_T(scratch_data.dof_hand_T)
837  , dofs_T(scratch_data.dofs_T)
838  , type_of_pde_rhs(scratch_data.type_of_pde_rhs)
839  , eta_squared(scratch_data.eta_squared)
840  , ve(0)
841  , se(0)
842 {
843 }
844 
845 template<int dim, int stage>
846 void
847 Solver2<dim, stage>::system_matrix_local(const IteratorPair& IP,
848  AssemblyScratchData& scratch_data,
849  AssemblyCopyData& copy_data)
850 {
851  // See the following boxes:
852  // (1) Recipe for static vector solver in 3D
853  // (2) Recipe for static vector solver in 2D
854 
855  // The comments below refer to these recipes by number, i.e., recipe (1)
856  // and recipe (2).
857 
858  copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
859  scratch_data.dofs_per_cell);
860 
861  copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
862 
863  copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
864 
865  auto cell = std::get<0>(*IP);
866  auto cell_T = std::get<1>(*IP);
867 
868  scratch_data.fe_values.reinit(cell);
869  scratch_data.fe_values_T.reinit(cell_T);
870 
871  scratch_data.the_coefficient.value_list(
872  scratch_data.fe_values.get_quadrature_points(),
873  cell->material_id(),
874  cell->user_index(),
875  scratch_data.the_coefficient_list);
876 
877  if (dim == 2) {
878  scratch_data.fe_values_T[SE].get_function_values(scratch_data.dofs_T,
879  scratch_data.values);
880  } else if (dim == 3) {
881  scratch_data.fe_values_T[VE].get_function_values(
882  scratch_data.dofs_T, scratch_data.vector_values);
883  }
884 
885  for (unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index) {
886  for (unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i) {
887  for (unsigned int j = 0; j < scratch_data.dofs_per_cell;
888  ++j) { // Integral I_a1+I_a3 in recipes (1) and (2).
889  copy_data.cell_matrix(i, j) +=
890  ((1.0 / scratch_data.the_coefficient_list[q_index]) * // 1 / mu
891  scratch_data.fe_values[VE].curl(i, q_index) * // curl N_i
892  scratch_data.fe_values[VE].curl(j, q_index) // curl N_j
893  + scratch_data.eta_squared * // eta^2
894  scratch_data.fe_values[VE].value(i, q_index) * // N_i
895  scratch_data.fe_values[VE].value(j, q_index) // N_j
896  ) *
897  scratch_data.fe_values.JxW(q_index); // dV (dS in 2D)
898  }
899 
900  if (dim == 2) { // Integral I_b3-1 in recipe (2).
901  copy_data.cell_rhs(i) +=
902  scratch_data.values[q_index] * // T
903  scratch_data.fe_values[VE].curl(i, q_index)[0] * // curl_s N_i
904  scratch_data.fe_values.JxW(q_index); // dS
905  } else if (dim == 3) { // Integral I_b3-1 in recipe (1).
906  copy_data.cell_rhs(i) +=
907  (scratch_data.vector_values[q_index][0] *
908  scratch_data.fe_values[VE].curl(i, q_index)[0] +
909  scratch_data.vector_values[q_index][1] *
910  scratch_data.fe_values[VE].curl(i, q_index)[1] +
911  scratch_data.vector_values[q_index][2] *
912  scratch_data.fe_values[VE].curl(i, q_index)[2] // T
913  ) * // curl N_i
914  scratch_data.fe_values.JxW(q_index); // dV
915  }
916  }
917  }
918 
919  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) {
920  scratch_data.do_robin = ((cell->face(f)->at_boundary()) &&
921  (cell->face(f)->boundary_id() % 2 == 0) &&
922  (cell->face(f)->boundary_id() != 0));
923 
924  scratch_data.do_K =
925  ((cell->user_index() > 0) && (cell->face(f)->user_index() > 0));
926 
927  scratch_data.do_T_on_boundary =
928  ((cell->face(f)->at_boundary()) && (type_of_pde_rhs != 2));
929 
930  Assert(!(scratch_data.do_robin && scratch_data.do_K), ExcInternalError());
931 
932  if (scratch_data.do_robin || scratch_data.do_K ||
933  scratch_data.do_T_on_boundary) {
934 
935  scratch_data.fe_face_values.reinit(cell, f);
936  scratch_data.fe_face_values_T.reinit(cell_T, f);
937 
938  if (scratch_data.do_robin) {
939  scratch_data.gamma.value_list(
940  scratch_data.fe_face_values.get_quadrature_points(),
941  scratch_data.fe_face_values.get_normal_vectors(),
942  cell->face(f)->boundary_id(),
943  cell->material_id(),
944  cell->user_index(),
945  cell->face(f)->user_index(),
946  scratch_data.gamma_list);
947 
948  scratch_data.robin_rhs.value_list(
949  scratch_data.fe_face_values.get_quadrature_points(),
950  scratch_data.fe_face_values.get_normal_vectors(),
951  cell->face(f)->boundary_id(),
952  cell->material_id(),
953  cell->user_index(),
954  cell->face(f)->user_index(),
955  scratch_data.robin_rhs_list);
956  }
957 
958  if (scratch_data.do_K) {
959  scratch_data.free_surface_current.value_list(
960  scratch_data.fe_face_values.get_quadrature_points(),
961  scratch_data.fe_face_values.get_normal_vectors(),
962  cell->material_id(),
963  cell->user_index(),
964  cell->face(f)->user_index(),
965  scratch_data.free_surface_current_list);
966  }
967 
968  if (scratch_data.do_T_on_boundary) {
969  if (dim == 2) {
970  scratch_data.fe_face_values_T[SE].get_function_values(
971  scratch_data.dofs_T, scratch_data.values_face);
972  } else if (dim == 3) {
973  scratch_data.fe_face_values_T[VE].get_function_values(
974  scratch_data.dofs_T, scratch_data.vector_values_face);
975  }
976  }
977 
978  for (unsigned int q_index_face = 0;
979  q_index_face < scratch_data.n_q_points_face;
980  ++q_index_face) {
981  for (unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i) {
982  if (scratch_data.do_robin) {
983  for (unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j) {
984  if (dim == 2) { // Integral I_a2 in recipe (2)
985  copy_data.cell_matrix(i, j) +=
986  scratch_data.gamma_list[q_index_face] * // gamma
987  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
988  scratch_data.fe_face_values[VE].value(i, q_index_face)[1] -
989  scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
990  scratch_data.fe_face_values[VE].value(i,
991  q_index_face)[0]) *
992  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
993  scratch_data.fe_face_values[VE].value(j, q_index_face)[1] -
994  scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
995  scratch_data.fe_face_values[VE].value(j,
996  q_index_face)[0]) *
997  scratch_data.fe_face_values.JxW(
998  q_index_face); // (n xs N_i)(n xs N_j) dl
999  } else if (dim == 3) { // Integral I_a2 in recipe (1).
1000  copy_data.cell_matrix(i, j) +=
1001  scratch_data.gamma_list[q_index_face] * // gamma
1002  ((scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
1003  scratch_data.fe_face_values[VE].value(i,
1004  q_index_face)[2] -
1005  scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
1006  scratch_data.fe_face_values[VE].value(i,
1007  q_index_face)[1]) *
1008  (scratch_data.fe_face_values.normal_vector(
1009  q_index_face)[1] *
1010  scratch_data.fe_face_values[VE].value(j,
1011  q_index_face)[2] -
1012  scratch_data.fe_face_values.normal_vector(
1013  q_index_face)[2] *
1014  scratch_data.fe_face_values[VE].value(
1015  j, q_index_face)[1]) +
1016  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
1017  scratch_data.fe_face_values[VE].value(i,
1018  q_index_face)[2] -
1019  scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
1020  scratch_data.fe_face_values[VE].value(i,
1021  q_index_face)[0]) *
1022  (scratch_data.fe_face_values.normal_vector(
1023  q_index_face)[0] *
1024  scratch_data.fe_face_values[VE].value(j,
1025  q_index_face)[2] -
1026  scratch_data.fe_face_values.normal_vector(
1027  q_index_face)[2] *
1028  scratch_data.fe_face_values[VE].value(
1029  j, q_index_face)[0]) +
1030  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
1031  scratch_data.fe_face_values[VE].value(i,
1032  q_index_face)[1] -
1033  scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
1034  scratch_data.fe_face_values[VE].value(i,
1035  q_index_face)[0]) *
1036  (scratch_data.fe_face_values.normal_vector(
1037  q_index_face)[0] *
1038  scratch_data.fe_face_values[VE].value(j,
1039  q_index_face)[1] -
1040  scratch_data.fe_face_values.normal_vector(
1041  q_index_face)[1] *
1042  scratch_data.fe_face_values[VE].value(
1043  j, q_index_face)[0])) // (n x N_i) . (n x N_j)
1044  * scratch_data.fe_face_values.JxW(q_index_face); // dS
1045  } else {
1046  Assert(false, ExcInternalError());
1047  }
1048  }
1049 
1050  // Integral I_b1 in recipes (1) and (2).
1051  copy_data.cell_rhs(i) +=
1052  -scratch_data.robin_rhs_list[q_index_face] * // Q
1053  scratch_data.fe_face_values[VE].value(i, q_index_face) * // N_i
1054  scratch_data.fe_face_values.JxW(q_index_face); // dS (dl in 2D)
1055  } // if (scratch_data.do_robin)
1056 
1057  if (scratch_data.do_K) { // Integral I_b2 in recipes (1) and (2).
1058  copy_data.cell_rhs(i) +=
1059  scratch_data.free_surface_current_list[q_index_face] * // K_f
1060  scratch_data.fe_face_values[VE].value(i, q_index_face) * // N_i
1061  scratch_data.fe_face_values.JxW(q_index_face); // dS (dl in 2D)
1062  }
1063 
1064  if (scratch_data.do_T_on_boundary) {
1065  if (dim == 2) { // Integral I_b3-2 in recipe (2).
1066  copy_data.cell_rhs(i) -=
1067  scratch_data.values_face[q_index_face] * // T
1068  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
1069  scratch_data.fe_face_values[VE].value(i, q_index_face)[1] -
1070  scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
1071  scratch_data.fe_face_values[VE].value(
1072  i, q_index_face)[0]) * // s
1073  scratch_data.fe_face_values.JxW(q_index_face); // (n x N_i) dl
1074  } else if (dim == 3) { // Integral I_b3-2 in recipe (1).
1075  copy_data.cell_rhs(i) -=
1076  (scratch_data.vector_values_face[q_index_face][0] *
1077  (scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
1078  scratch_data.fe_face_values[VE].value(i,
1079  q_index_face)[2] -
1080  scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
1081  scratch_data.fe_face_values[VE].value(i,
1082  q_index_face)[1]) -
1083  scratch_data.vector_values_face[q_index_face][1] *
1084  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
1085  scratch_data.fe_face_values[VE].value(i,
1086  q_index_face)[2] -
1087  scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
1088  scratch_data.fe_face_values[VE].value(i,
1089  q_index_face)[0]) +
1090  scratch_data.vector_values_face[q_index_face][2] *
1091  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
1092  scratch_data.fe_face_values[VE].value(i,
1093  q_index_face)[1] -
1094  scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
1095  scratch_data.fe_face_values[VE].value(
1096  i, q_index_face)[0]) // T . (n x N_i )
1097  ) *
1098  scratch_data.fe_face_values.JxW(q_index_face); // dS
1099  } else {
1100  Assert(false, ExcInternalError());
1101  }
1102  }
1103  }
1104  }
1105  }
1106  }
1107  cell->get_dof_indices(copy_data.local_dof_indices);
1108 }
1109 
1110 template<int dim, int stage>
1111 void
1112 Solver2<dim, stage>::copy_local_to_global(const AssemblyCopyData& copy_data)
1113 {
1114  constraints.distribute_local_to_global(copy_data.cell_matrix,
1115  copy_data.cell_rhs,
1116  copy_data.local_dof_indices,
1117  system_matrix,
1118  system_rhs);
1119 }
1120 
1121 template<int dim, int stage>
1122 void
1123 Solver2<dim, stage>::compute_error_norms()
1124 {
1125  Weight<dim, stage> weight;
1126  const Function<dim, double>* mask = &weight;
1127 
1128  Constants::QuadratureTableVector<dim> qt(dof_handler.get_fe().degree - 1);
1129  QGauss<dim> quadrature(qt.enorm());
1130 
1131  VectorTools::integrate_difference(MappingQ<dim>(mapping_degree),
1132  dof_handler,
1133  solution,
1134  *exact_solution,
1135  L2_per_cell,
1136  quadrature,
1137  VectorTools::L2_norm,
1138  mask);
1139 
1140  L2_norm = VectorTools::compute_global_error(
1141  triangulation_T, L2_per_cell, VectorTools::L2_norm);
1142 
1143  VectorTools::integrate_difference(MappingQ<dim>(mapping_degree),
1144  dof_handler,
1145  solution,
1146  *exact_solution,
1147  Linfty_per_cell,
1148  QGauss<dim>(1),
1149  VectorTools::Linfty_norm,
1150  mask);
1151 
1152  Linfty_norm = Linfty_per_cell.linfty_norm();
1153 }
1154 
1155 template<int dim, int stage>
1156 void
1157 Solver2<dim, stage>::project_exact_solution_fcn()
1158 {
1159  Constants::QuadratureTableVector<dim> qt(fe.degree - 1);
1160 
1161  AffineConstraints<double> constraints_empty;
1162  constraints_empty.close();
1163 
1164  VectorTools::project(MappingQ<dim>(mapping_degree),
1165  dof_handler,
1166  constraints_empty,
1167  QGauss<dim>(qt.sim()),
1168  *exact_solution,
1169  projected_exact_solution);
1170 }
1171 
1172 template<int dim, int stage>
1173 void
1174 Solver2<dim, stage>::save() const
1175 {
1176  std::vector<std::string> solution_names(dim, "VectorField");
1177  std::vector<DataComponentInterpretation::DataComponentInterpretation>
1178  interpretation(dim,
1179  DataComponentInterpretation::component_is_part_of_vector);
1180 
1181  DataOut<dim> data_out;
1182 
1183  data_out.add_data_vector(
1184  dof_handler, solution, solution_names, interpretation);
1185 
1186  if (project_exact_solution && exact_solution) {
1187  std::vector<std::string> solution_names_ex(dim, "VectorFieldExact");
1188 
1189  data_out.add_data_vector(
1190  dof_handler, projected_exact_solution, solution_names_ex, interpretation);
1191  }
1192 
1193  if (exact_solution) {
1194  data_out.add_data_vector(L2_per_cell, "L2norm");
1195  data_out.add_data_vector(Linfty_per_cell, "LinftyNorm");
1196  }
1197 
1198  std::ofstream ofs;
1199 
1200  if (write_higher_order_cells) {
1201  DataOutBase::VtkFlags flags;
1202  flags.write_higher_order_cells = true;
1203  data_out.set_flags(flags);
1204 
1205  const MappingQ<dim> mapping(mapping_degree);
1206 
1207  data_out.build_patches(mapping,
1208  fe.degree + 2,
1209  DataOut<dim>::CurvedCellRegion::curved_inner_cells);
1210 
1211  ofs.open(fname + ".vtu");
1212  data_out.write_vtu(ofs);
1213 
1214  } else {
1215 
1216  data_out.build_patches();
1217 
1218  ofs.open(fname + ".vtk");
1219  data_out.write_vtk(ofs);
1220  }
1221 
1222  ofs.close();
1223 }
1224 
1225 template<int dim, int stage>
1226 void
1227 Solver2<dim, stage>::save_matrix_and_rhs_to_csv(std::string fname) const
1228 {
1229  std::ofstream ofs_matrix(fname + "_matrix.csv");
1230  std::ofstream ofs_rhs(fname + "_rhs.csv");
1231 
1232  for (unsigned int i = 0; i < system_matrix.m(); ++i) {
1233  ofs_rhs << system_rhs(i);
1234  if (i < (system_matrix.m() - 1))
1235  ofs_rhs << "\n";
1236 
1237  for (unsigned int j = 0; j < system_matrix.n(); ++j) {
1238  ofs_matrix << std::scientific << std::setprecision(16)
1239  << system_matrix.el(i, j);
1240 
1241  if (j < (system_matrix.m() - 1))
1242  ofs_matrix << ", ";
1243  }
1244  if (i < (system_matrix.m() - 1))
1245  ofs_matrix << "\n";
1246  }
1247 
1248  ofs_rhs.close();
1249  ofs_matrix.close();
1250 }
1251 
1252 } // namespace StaticVectorSolver
1253 
1254 #endif
Constants::QuadratureTableVector
The tables that contain the amount of quadrature points used in vector problems.
Definition: constants.hpp:101
Constants::QuadratureTableVector::enorm
unsigned int enorm() const
Returns the amount of quadrature points used when calculating the error norms.
Constants::QuadratureTableVector::sim
unsigned int sim() const
Returns the amount of quadrature points used when assembling system if linear equations.
StaticVectorSolver::Gamma
Implements the coefficient on the left-hand side of the Robin boundary condition (iii) of the static...
Definition: static_vector_input.hpp:221
StaticVectorSolver::RobinRhs
Implements the vector field, , on the right-hand side of the Robin boundary condition (iii) of the st...
Definition: static_vector_input.hpp:312
StaticVectorSolver::Solver2
Solves static vector boundary value problem.
Definition: static_vector_solver_ii.hpp:224
StaticVectorSolver::Solver2::save_matrix_and_rhs_to_csv
void save_matrix_and_rhs_to_csv(std::string fname) const
Saves the system matrix and the right-hand side into a csv file.
Definition: static_vector_solver_ii.hpp:1227
StaticVectorSolver::Solver2::assemble
void assemble()
Assembles the system matrix and the right-hand side vector.
Definition: static_vector_solver_ii.hpp:736
StaticVectorSolver::Solver2::dirichlet_stack
std::map< types::boundary_id, const Function< dim > * > dirichlet_stack
A map that contains pairs of boundary IDs and the corresponding Dirichlet boundary conditions....
Definition: static_vector_solver_ii.hpp:544
StaticVectorSolver::Solver2::clear
void clear()
Releases computer memory associated with system matrix and right-hand side.
Definition: static_vector_solver_ii.hpp:423
StaticVectorSolver::Solver2::get_mapping_degree
unsigned int get_mapping_degree() const
Returns degree of the interpolating Lagrange polynomials used for mapping from the reference cell to ...
Definition: static_vector_solver_ii.hpp:469
StaticVectorSolver::Solver2::fe
const FE_Nedelec< dim > fe
The finite elements.
Definition: static_vector_solver_ii.hpp:549
StaticVectorSolver::Solver2::system_matrix
SparseMatrix< double > system_matrix
The system matrix.
Definition: static_vector_solver_ii.hpp:580
StaticVectorSolver::Solver2::get_solution
const Vector< double > & get_solution() const
Returns a reference to solution.
Definition: static_vector_solver_ii.hpp:437
StaticVectorSolver::Solver2::solve
virtual void solve()=0
Solves the system of linear equations.
StaticVectorSolver::Solver2::dof_handler
DoFHandler< dim > dof_handler
The finite elements.
Definition: static_vector_solver_ii.hpp:554
StaticVectorSolver::Solver2::solution
Vector< double > solution
The solution vector, that is, degrees of freedom yielded by the simulation.
Definition: static_vector_solver_ii.hpp:560
StaticVectorSolver::Solver2::Solver2
Solver2(unsigned int p, unsigned int mapping_degree, const Triangulation< dim > &triangulation_T, const DoFHandler< dim > &dof_handler_T, const Vector< double > &solution_T, unsigned int type_of_pde_rhs=0, double eta_squared=0.0, std::string fname="data", const Function< dim > *exact_solution=nullptr, bool print_time_tables=false, bool project_exact_solution=false, bool write_higher_order_cells=false)
The only constructor.
Definition: static_vector_solver_ii.hpp:265
StaticVectorSolver::Solver2::system_rhs
Vector< double > system_rhs
The system right-hand side vector.
Definition: static_vector_solver_ii.hpp:585
StaticVectorSolver::Solver2::triangulation_T
const Triangulation< dim > & triangulation_T
Reference to the mesh.
Definition: static_vector_solver_ii.hpp:526
StaticVectorSolver::Solver2::compute_error_norms
void compute_error_norms()
Computes error norms.
Definition: static_vector_solver_ii.hpp:1123
StaticVectorSolver::Solver2::project_exact_solution_fcn
void project_exact_solution_fcn()
Projects exact solution.
Definition: static_vector_solver_ii.hpp:1157
StaticVectorSolver::Solver2::setup
void setup()
Initializes system matrix and the right-hand side vector.
Definition: static_vector_solver_ii.hpp:691
StaticVectorSolver::Solver2::get_n_dofs
unsigned int get_n_dofs() const
Returns the total amount of the degrees of freedom.
Definition: static_vector_solver_ii.hpp:450
StaticVectorSolver::Solver2::get_dof_handler
const DoFHandler< dim > & get_dof_handler() const
Returns a reference to a dof handler.
Definition: static_vector_solver_ii.hpp:432
StaticVectorSolver::Solver2::get_n_cells
unsigned int get_n_cells() const
Returns the number of active cells in the mesh.
Definition: static_vector_solver_ii.hpp:442
StaticVectorSolver::Solver2::sparsity_pattern
SparsityPattern sparsity_pattern
The sparsity pattern of the system matrix.
Definition: static_vector_solver_ii.hpp:575
StaticVectorSolver::Solver2::solution_T
const Vector< double > & solution_T
Reference to the degrees of freedom that describe the current vector potential, .
Definition: static_vector_solver_ii.hpp:538
StaticVectorSolver::Solver2::get_Linfty_norm
double get_Linfty_norm() const
Returns error norm.
Definition: static_vector_solver_ii.hpp:463
StaticVectorSolver::Solver2::run
void run()
Runs the simulation.
Definition: static_vector_solver_ii.hpp:478
StaticVectorSolver::Solver2::save
void save() const
Saves simulation results into a vtk file.
Definition: static_vector_solver_ii.hpp:1174
StaticVectorSolver::Solver2::dof_handler_T
const DoFHandler< dim > & dof_handler_T
Reference to the dof handler that describes the current vector potential, .
Definition: static_vector_solver_ii.hpp:532
StaticVectorSolver::Solver2::fill_dirichlet_stack
virtual void fill_dirichlet_stack()=0
Initializes the data member StaticVectorSolver::Solver2::dirichlet_stack.
StaticVectorSolver::Solver2::projected_exact_solution
Vector< double > projected_exact_solution
The projected exact solution vector.
Definition: static_vector_solver_ii.hpp:565
StaticVectorSolver::Solver2::constraints
AffineConstraints< double > constraints
The constraints associated with the Dirichlet boundary conditions.
Definition: static_vector_solver_ii.hpp:570
StaticVectorSolver::Solver2::get_L2_norm
double get_L2_norm() const
Returns error norm.
Definition: static_vector_solver_ii.hpp:458
StaticVectorSolver::TheCoefficient
Implements the permeability, , in the partial differential equation (i) of the vector boundary value ...
Definition: static_vector_input.hpp:53
StaticVectorSolver::Weight
Implements the weight function for calculating the and error norms.
Definition: static_vector_input.hpp:475