API tools faq
paste
Advertisement
wq6

Untitled

wq6
Jun 19th, 2019
97
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 23.94 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.  * O
  3.  *
  4.  * Author: Gabriel Brandao de Miranda
  5.  * Date: 20/04/19.
  6. */
  7. #include <deal.II/base/quadrature_lib.h>
  8. #include <deal.II/base/logstream.h>
  9. #include <deal.II/base/function.h>
  10. #include <deal.II/base/tensor_function.h>
  11. #include <deal.II/base/convergence_table.h>
  12.  
  13. #include <deal.II/lac/full_matrix.h>
  14. #include <deal.II/lac/sparse_direct.h>
  15. #include <deal.II/lac/precondition.h>
  16. #include <deal.II/lac/solver_bicgstab.h>
  17. #include <deal.II/lac/solver_gmres.h>
  18. #include <deal.II/lac/precondition.h>
  19.  
  20. #include <deal.II/grid/tria.h>
  21. #include <deal.II/grid/grid_generator.h>
  22. #include <deal.II/grid/tria_accessor.h>
  23. #include <deal.II/grid/tria_iterator.h>
  24. #include <deal.II/grid/grid_tools.h>
  25.  
  26. #include <deal.II/dofs/dof_handler.h>
  27. #include <deal.II/dofs/dof_renumbering.h>
  28. #include <deal.II/dofs/dof_accessor.h>
  29. #include <deal.II/dofs/dof_tools.h>
  30.  
  31. #include <deal.II/fe/fe_dgq.h>
  32. #include <deal.II/fe/fe_system.h>
  33. #include <deal.II/fe/fe_values.h>
  34. #include <deal.II/fe/fe_raviart_thomas.h>
  35. #include <deal.II/fe/fe_dg_vector.h>
  36. #include <deal.II/fe/fe_face.h>   // Para problema com multiplicador descontinuo
  37. #include <deal.II/fe/fe_trace.h>  // Para problema com multiplicador continuo
  38.  
  39. #include <deal.II/numerics/vector_tools.h>
  40. #include <deal.II/numerics/matrix_tools.h>
  41. #include <deal.II/numerics/data_out.h>
  42. #include <deal.II/numerics/data_out_faces.h>
  43.  
  44. #include <fstream>
  45. #include <iostream>
  46. #include <cmath>
  47. #include <math.h>
  48.  
  49. namespace PoissonProblem {
  50.     using namespace dealii;
  51.     const double epsilon = 1.0;
  52.     const double coef_reat = 1.0;
  53.     const double coef_conv = 0.0;
  54.  
  55.     template <int dim>
  56.     class SourceTerm : public Function<dim> {
  57.     public:
  58.  
  59.         SourceTerm() : Function<dim>() { }
  60.  
  61.         virtual double value(const Point<dim> &p, const unsigned int component = 0) const {
  62.             double return_value = epsilon * dim * pow(M_PI,2) * exp( -dim * pow(M_PI,2) * this->get_time() );
  63.             return_value -=     coef_reat * dim * pow(M_PI,2) * exp( -dim * pow(M_PI,2) * this->get_time() );
  64.  
  65.             for (int i = 0; i < dim; ++i) {
  66.                 return_value *= sin(M_PI * p[i]);
  67.             }
  68.  
  69.             for (unsigned int i = 0; i < dim; ++i) {
  70.                 double gradient = coef_conv * exp( -dim * this->get_time() * pow(M_PI,2) );
  71.                 for (unsigned int j = 0; j < dim; ++j) {
  72.                     if (i == j) {
  73.                         gradient *= cos(M_PI * p[j]);
  74.                     } else {
  75.                         gradient *= sin(M_PI * p[j]);
  76.                     }
  77.                 }
  78.                 return_value += gradient;
  79.             }
  80.             return return_value;
  81.         }
  82.     };
  83.  
  84.     template <int dim>
  85.     class InitialCondition : public Function<dim> {
  86.     public:
  87.         InitialCondition() : Function<dim>() { }
  88.  
  89.         virtual double value(const Point<dim> &p, const unsigned int component = 0) const {
  90.             double return_value = 1.0;
  91.  
  92.             for (int i = 0; i < dim; ++i)       return_value *= sin(M_PI * p[i]);
  93.  
  94.             return return_value;
  95.         }
  96.     };
  97.  
  98.  
  99.     template <int dim>
  100.     class ExactSolution : public Function<dim> {
  101.     public:
  102.  
  103.         ExactSolution() : Function<dim>() {
  104.         }
  105.  
  106.         virtual double value(const Point<dim> &p, const unsigned int component = 0) const {
  107.             double return_value = exp( -dim * pow(M_PI,2) * this->get_time() );
  108.             for (int i = 0; i < dim; ++i)       return_value *= sin(M_PI * p[i]);
  109.             return return_value;
  110.         }
  111.  
  112.         virtual Tensor<1, dim> gradient(const Point<dim> &p, const unsigned int component = 0) const {
  113.             Tensor<1, dim> return_value;
  114.  
  115.             for (unsigned int i = 0; i < dim; ++i) {
  116.                 return_value[i] = M_PI * exp( -dim * this->get_time() *pow(M_PI,2) );
  117.  
  118.                 for (unsigned int j = 0; j < dim; ++j) {
  119.                     if (i == j) {
  120.                         return_value[i] *= cos(M_PI * p[j]);
  121.                     } else {
  122.                         return_value[i] *= sin(M_PI * p[j]);
  123.                     }
  124.                 }
  125.             }
  126.             return return_value;
  127.         }
  128.     };
  129.  
  130.     template <int dim>
  131.     class HybridPoissonProblem {
  132.     public:
  133.  
  134.         HybridPoissonProblem(const unsigned int degree) :
  135.                 degree(degree), fe_local(degree), dof_handler_local(triangulation), fe(degree), dof_handler(triangulation) {
  136.         }
  137.  
  138.         void run(const int refinements, ConvergenceTable &convergence_table, const double final_time ) {
  139.             double current_time = 0.0;
  140.             InitialCondition<dim> initial_condition;
  141.             make_grid_and_dofs( refinements, current_time );
  142.             VectorTools::interpolate( dof_handler_local, initial_condition, solution_local );
  143.  
  144. //            DataOut<dim> data_out;
  145. //            data_out.attach_dof_handler( dof_handler_local );
  146. //            data_out.add_data_vector( solution_local_last, "inicial" );
  147. //            data_out.build_patches();
  148. //            std::string file_name = "inicial.vtk";
  149. //            std::ofstream output(file_name.c_str());
  150. //            data_out.write_vtk(output);
  151.  
  152.             while( current_time < final_time ) {
  153.                 current_time += dt;
  154.                 assemble_system( true, current_time);
  155.                 change_boundary( current_time );
  156.                 solve();
  157.                 assemble_system( false, current_time );
  158. //                solution_local_last = solution_local;      /// Altera solucao por solucao anterior
  159.                 std::cout << "   Current time:  " << current_time << std::endl;
  160.  
  161.                 system_rhs.reinit(dof_handler.n_dofs());
  162.                 system_matrix.reinit(sparsity_pattern);
  163.             }
  164.  
  165.             ExactSolution<dim> solution_function;
  166.             solution_function.set_time( current_time );
  167.             output_results(refinements);
  168.             compute_errors(convergence_table, current_time);
  169.         }
  170.  
  171.     private:
  172.         const unsigned int degree;
  173.         const double epsilon = 1.0;
  174.         double h_mesh;
  175.         double dt;
  176.         Triangulation<dim> triangulation;
  177.  
  178.         FE_DGQ<dim> fe_local;
  179.         DoFHandler<dim> dof_handler_local;
  180.         Vector<double> solution_local;
  181.         Vector<double> solution_local_last;
  182.  
  183.         FE_FaceQ<dim>   fe;
  184.         DoFHandler<dim> dof_handler;
  185.         Vector<double>  solution;
  186.         Vector<double>  system_rhs;
  187.  
  188.         SparsityPattern sparsity_pattern;
  189.         SparseMatrix<double> system_matrix;
  190.  
  191.         void change_boundary( const double current_time ){
  192.  
  193.             std::map<types::global_dof_index,double> dir_boundary_values;
  194.             ExactSolution<dim> solution_function;
  195.  
  196.             solution_function.set_time( current_time );
  197.             VectorTools::interpolate_boundary_values( dof_handler, 0, solution_function, dir_boundary_values);
  198.  
  199.             MatrixTools::apply_boundary_values(dir_boundary_values, system_matrix, solution, system_rhs);
  200.  
  201. //            std::map<types::boundary_id, const Function<dim> *> boundary_functions;
  202. //            ExactSolution<dim> solution_function;
  203. //            solution_function.set_time(current_time);
  204. //            boundary_functions[0] = &solution_function;
  205. //            VectorTools::project_boundary_values(dof_handler, boundary_functions, QGauss < dim - 1 > (fe.degree + 1), constraints);
  206.         }
  207.  
  208.         void make_grid_and_dofs(int i, const double current_time ) {
  209.  
  210.             GridGenerator::hyper_cube(triangulation, 0, 1.0);
  211.             triangulation.refine_global(i);
  212.  
  213.             h_mesh = GridTools::maximal_cell_diameter(triangulation);
  214.             dt = std::pow( h_mesh, degree + 1 );
  215.  
  216.             dof_handler.distribute_dofs(fe);
  217.             dof_handler_local.distribute_dofs(fe_local);
  218.  
  219.             std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl
  220.                       << "Total number of cells: " << triangulation.n_cells() << std::endl
  221.                       << "Number of degrees of freedom for the multiplier: " << dof_handler.n_dofs() << std::endl << std::endl;
  222.  
  223.             solution.reinit(dof_handler.n_dofs());
  224.             solution_local.reinit(dof_handler_local.n_dofs());
  225.             solution_local_last.reinit(dof_handler_local.n_dofs());
  226.  
  227.             DynamicSparsityPattern dsp(dof_handler.n_dofs());
  228.             DoFTools::make_sparsity_pattern( dof_handler, dsp );
  229.             sparsity_pattern.copy_from(dsp);
  230.  
  231.             system_rhs.reinit(dof_handler.n_dofs());
  232.             system_matrix.reinit(sparsity_pattern);
  233.  
  234. //            DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  235.         }
  236.  
  237.         void assemble_system(bool globalProblem , const double current_time) {
  238.  
  239.             QGauss < dim >          quadrature_formula(degree + 2);
  240.             QGauss < dim - 1 > face_quadrature_formula(degree + 2);
  241.  
  242.             FEValues<dim>     fe_values_local(fe_local, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values);
  243.             FEFaceValues<dim> fe_face_values(fe, face_quadrature_formula, update_values | update_normal_vectors | update_quadrature_points | update_JxW_values);
  244.             FEFaceValues<dim> fe_face_values_local(fe_local, face_quadrature_formula, update_values | update_gradients | update_normal_vectors | update_quadrature_points | update_JxW_values);
  245.  
  246.             const unsigned int dofs_local_per_cell = fe_local.dofs_per_cell;
  247.             const unsigned int n_q_points = quadrature_formula.size();
  248.             const unsigned int n_face_q_points = face_quadrature_formula.size();
  249.  
  250.             std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
  251.  
  252.             FullMatrix<double> A_matrix(dofs_local_per_cell, dofs_local_per_cell);
  253.             FullMatrix<double> B_matrix(dofs_local_per_cell, fe.dofs_per_cell);
  254.             FullMatrix<double> BT_matrix(fe.dofs_per_cell, dofs_local_per_cell);
  255.             FullMatrix<double> C_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
  256.             Vector<double> F_vector(dofs_local_per_cell);
  257.             Vector<double> U_vector(dofs_local_per_cell);
  258.             Vector<double> cell_vector(fe.dofs_per_cell);
  259.  
  260.             FullMatrix<double> aux_matrix(dofs_local_per_cell, fe.dofs_per_cell);
  261.             Vector<double> aux_vector(dofs_local_per_cell);
  262.  
  263.             std::vector<double> mult_values(face_quadrature_formula.size());
  264.  
  265.             std::vector<double> st_values(n_q_points);
  266.             std::vector<double> prev_solution( n_q_points );
  267.             std::vector<double> boundary_values(n_face_q_points);
  268.  
  269.             for (const auto &cell : dof_handler.active_cell_iterators()) {
  270.  
  271.                 Tensor<1, dim> b_velocity;
  272.                 for( int i = 0; i < dim; i++)
  273.                     b_velocity[i] = coef_conv;
  274.  
  275.                 const double beta_stab = degree*degree*100/h_mesh ;
  276.  
  277.                 typename DoFHandler<dim>::active_cell_iterator loc_cell(&triangulation, cell->level(), cell->index(), &dof_handler_local);
  278.                 fe_values_local.reinit(loc_cell);
  279.  
  280.                 A_matrix  = 0.;
  281.                 F_vector  = 0.;
  282.                 B_matrix  = 0.;
  283.                 BT_matrix = 0.;
  284.                 C_matrix  = 0.;
  285.                 U_vector  = 0.;
  286.  
  287.                 SourceTerm<dim> source_term;
  288.                 source_term.set_time( current_time );
  289.                 source_term.value_list(fe_values_local.get_quadrature_points(), st_values);
  290.  
  291.                 fe_values_local.get_function_values( solution_local , prev_solution );
  292.  
  293.                 for (unsigned int q = 0; q < n_q_points; ++q) {
  294.                     const double JxW = fe_values_local.JxW(q);
  295.                     for (unsigned int i = 0; i < dofs_local_per_cell; ++i) {
  296.                         const double phi_i_u = fe_values_local.shape_value(i, q);
  297.                         const Tensor<1, dim> grad_phi_i_u = fe_values_local.shape_grad(i,q);
  298.  
  299.                         for (unsigned int j = 0; j < dofs_local_per_cell; ++j) {
  300.                             const double phi_j_u = fe_values_local.shape_value(j, q);
  301.                             const Tensor<1, dim> grad_phi_j_u = fe_values_local.shape_grad(j, q);
  302.  
  303.                             A_matrix(i, j)  +=  (   dt*epsilon * grad_phi_j_u * grad_phi_i_u   // < grad(u) , grad(v) >
  304.                                                 +   dt*b_velocity*grad_phi_j_u * phi_i_u       // < b*grad(u) , v >
  305.                                                 +   coef_reat * phi_j_u * phi_i_u      // < u , v > / dt
  306.                                                 ) * JxW ;
  307.                         }
  308.                         F_vector(i) +=  (dt*phi_i_u * st_values[q] + phi_i_u * prev_solution[q]/dt ) * JxW;
  309.                     }
  310.                 }
  311.                 /// Loop nas faces dos elementos
  312.                 for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell; ++face_n) {
  313.                     fe_face_values.reinit(cell, face_n);
  314.                     fe_face_values_local.reinit(loc_cell, face_n);
  315.                     for (unsigned int q = 0; q < n_face_q_points; ++q) {
  316.  
  317.                         const double JxW = fe_face_values.JxW(q);
  318.                         const Tensor<1, dim> normal = fe_face_values.normal_vector(q);
  319.  
  320.                         for (unsigned int i = 0; i < dofs_local_per_cell; ++i){
  321.                             const Tensor<1, dim> grad_phi_i_u = fe_face_values_local.shape_grad(i, q);
  322.                             const double phi_i_u = fe_face_values_local.shape_value(i,q);
  323.  
  324.                             for (unsigned int j = 0; j < dofs_local_per_cell; ++j) {
  325.                                 const Tensor<1, dim> grad_phi_j_u = fe_face_values_local.shape_grad(j, q);
  326.                                 const double phi_j_u = fe_face_values_local.shape_value(j,q);
  327.  
  328.                                 A_matrix(i, j)  +=  (-   phi_j_u * grad_phi_i_u*normal
  329.                                                      -   grad_phi_j_u*normal * phi_i_u
  330.                                                      +   beta_stab * phi_j_u * phi_i_u
  331.                                                      +   phi_j_u * fmax(-b_velocity * normal,0) * phi_i_u    // <u,[b.n]-v>
  332.                                                     ) * dt * JxW ;
  333.                             }
  334.                         }
  335.                     }
  336.                 }
  337.                 if (globalProblem) {
  338.  
  339.                     /// Loop nas faces dos elementos
  340.                     for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell; ++face_n) {
  341.                         fe_face_values.reinit(cell, face_n);
  342.                         fe_face_values_local.reinit(loc_cell, face_n);
  343.                         for (unsigned int q = 0; q < n_face_q_points; ++q) {
  344.  
  345.                             const double JxW = fe_face_values.JxW(q);
  346.                             const Tensor<1, dim> normal = fe_face_values.normal_vector(q);
  347.  
  348.                             for (unsigned int i = 0; i < dofs_local_per_cell; ++i) {
  349.                                 const Tensor<1, dim> grad_phi_i_u = fe_face_values_local.shape_grad(i,q);
  350.                                 const double phi_i_u = fe_face_values_local.shape_value(i,q);
  351.  
  352.                                 for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) {
  353.                                     const double phi_j_m = fe_face_values.shape_value(j, q);
  354.  
  355.                                     B_matrix(i, j) +=   (   phi_j_m * (grad_phi_i_u * normal)
  356.                                                         -   beta_stab * phi_j_m * phi_i_u
  357.                                                         -   phi_j_m * fmax(-b_velocity * normal,0) * phi_i_u    // <lmbda,[b.n]-v>
  358.                                                         ) * dt * JxW;
  359.                                 }
  360.                             }
  361.                             for (unsigned int j = 0; j < dofs_local_per_cell; ++j) {
  362.                                 const Tensor<1, dim> grad_phi_j_u = fe_face_values_local.shape_grad(j,q);
  363.                                 const double phi_j_u = fe_face_values_local.shape_value(j,q);
  364.  
  365.                                 for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) {
  366.                                     const double phi_i_m = fe_face_values.shape_value(i, q);
  367.  
  368.  
  369.                                     BT_matrix(i, j) +=   (  phi_i_m * (grad_phi_j_u * normal)
  370.                                                         -   beta_stab * phi_i_m * phi_j_u
  371.                                                         -   phi_i_m * fmax(b_velocity * normal,0) * phi_j_u    // <lmbda,[b.n]-v>
  372.                                                          ) * JxW;
  373.                                 }
  374.                             }
  375.                             for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) {
  376.                                 const double phi_i_m = fe_face_values.shape_value(i, q);
  377.  
  378.                                 for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) {
  379.                                     const double phi_j_m = fe_face_values.shape_value(j, q);
  380.                                     /// Negativo para facilitar a condensacao estatica
  381.                                     C_matrix(i, j) -=   (   beta_stab * phi_j_m * phi_i_m
  382.                                                         +   phi_j_m * fmax(b_velocity * normal,0) * phi_i_m
  383.                                                         ) * JxW;
  384.                                 }
  385.                             }
  386.                         }
  387.                     }
  388.                 } else {
  389.  
  390.                     for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell; ++face_n) {
  391.                         fe_face_values.reinit(cell, face_n);
  392.                         fe_face_values_local.reinit(loc_cell, face_n);
  393.  
  394.                         // solution_function.value_list(fe_face_values.get_quadrature_points(), mult_values);  // Rotina0
  395.                         fe_face_values.get_function_values( solution, mult_values);
  396.  
  397.                         for (unsigned int q = 0; q < n_face_q_points; ++q) {
  398.                             const double JxW = fe_face_values_local.JxW(q);
  399.                             const Tensor<1, dim> normal = fe_face_values_local.normal_vector(q);
  400.  
  401.                             for (unsigned int i = 0; i < dofs_local_per_cell; ++i) {
  402.                                 const Tensor<1, dim> grad_phi_i_u = fe_face_values_local.shape_grad(i,q);
  403.                                 const double phi_i_u = fe_face_values_local.shape_value(i, q);
  404.  
  405.                                 F_vector(i) += (-   grad_phi_i_u * normal
  406.                                                 +   beta_stab * phi_i_u
  407.                                                 +   fmax(-b_velocity * normal,0) * phi_i_u
  408.                                                ) * dt * mult_values[q] * JxW;
  409.                             }
  410.                         }
  411.                     }
  412.                 }
  413.  
  414.  
  415.                 A_matrix.gauss_jordan(); // A = A^{-1}
  416.  
  417.  
  418.                 if (globalProblem) {
  419.  
  420.                     A_matrix.vmult(aux_vector, F_vector, false);            //  A^{-1} * F
  421.                     BT_matrix.vmult(cell_vector, aux_vector, false);        //  B.T * A^{-1} * F
  422.                     A_matrix.mmult(aux_matrix, B_matrix, false);            //  A^{-1} * B
  423.                     BT_matrix.mmult(C_matrix, aux_matrix, true);            // -C + B.T * A^{-1} * B
  424.                     cell->get_dof_indices(dof_indices);
  425.  
  426.                     std::vector<types::global_dof_index> local_dof_indices(fe.dofs_per_cell);
  427.                     cell->get_dof_indices(local_dof_indices);
  428.                     for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) {
  429.                         for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) {
  430.                             system_matrix.add(local_dof_indices[i], local_dof_indices[j], C_matrix(i, j));
  431.                         }
  432.                         system_rhs(local_dof_indices[i]) += cell_vector(i);
  433.                     }
  434.                 } else {
  435.                     A_matrix.vmult(U_vector, F_vector, false);              // A^{-1} * F
  436.                     loc_cell->set_dof_values(U_vector, solution_local);
  437.                 }
  438.             }
  439.         }
  440.  
  441.         void solve() {
  442.             PreconditionJacobi<SparseMatrix<double> > precondition;
  443.             precondition.initialize(system_matrix, .6);
  444.  
  445.             SolverControl solver_control( system_matrix.m()*10, 1e-11 * system_rhs.l2_norm() );
  446.             SolverBicgstab solver(solver_control);
  447.             solver.solve(system_matrix, solution, system_rhs, precondition);
  448.  
  449.             std::cout << std::endl << "   Number of BiCGStab iterations: " << solver_control.last_step() << std::endl;
  450.  
  451.         }
  452.  
  453. //        void solve() {
  454. //            SolverControl solver_control(1000, 1e-12);
  455. //            SolverCG<> solver(solver_control);
  456. //            solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
  457. //
  458. //            std::cout << "   " << solver_control.last_step()
  459. //                      << " CG iterations needed to obtain convergence."
  460. //                      << std::endl;
  461. //        }
  462.  
  463.         void output_results(const unsigned int refinements) const {
  464.             DataOut<dim> data_out;
  465.  
  466.             data_out.attach_dof_handler(dof_handler_local);
  467.             data_out.add_data_vector(solution_local, "solution");
  468.  
  469.             data_out.build_patches();
  470.  
  471.             std::string file_name = "solution_" + dealii::Utilities::to_string(refinements, 0) + ".vtk";
  472.             std::ofstream output(file_name.c_str());
  473.             data_out.write_vtk(output);
  474.         }
  475.  
  476.         void compute_errors(ConvergenceTable &convergence_table, const double current_time) const {
  477.             ExactSolution<dim> exact_solution;
  478.             exact_solution.set_time( current_time );
  479.  
  480.             Vector<double> cellwise_errors(triangulation.n_active_cells());
  481.  
  482.             QGauss<dim> quadrature_formula(degree + 2);
  483.  
  484.             Vector<double> elementwise_l2_errors;
  485.             VectorTools::integrate_difference(dof_handler_local, solution_local, exact_solution, cellwise_errors, quadrature_formula, VectorTools::L2_norm );
  486.  
  487.             double L2_error = VectorTools::compute_global_error(triangulation, cellwise_errors, VectorTools::L2_norm);
  488.  
  489.             Vector<double> elementwise_h1_errors;
  490.             VectorTools::integrate_difference(dof_handler_local, solution_local, exact_solution, cellwise_errors, quadrature_formula, VectorTools::H1_norm );
  491.             double H1_error = VectorTools::compute_global_error(triangulation, cellwise_errors, VectorTools::H1_norm);
  492.  
  493.             convergence_table.add_value("cells", triangulation.n_active_cells());
  494.             convergence_table.add_value("L2", L2_error);
  495.             convergence_table.add_value("H1", H1_error);
  496.         }
  497.     };
  498. }
  499.  
  500. int main() {
  501.     using namespace dealii;
  502.     using namespace PoissonProblem;
  503.  
  504.     const int dim           = 2;
  505.     const int degree        = 1;
  506.     const double final_time = 0.5;
  507.  
  508.     ConvergenceTable convergence_table;
  509.  
  510.     for (int refinement = 1; refinement < 5; ++refinement) {
  511.         std::cout << " --- Refinement #" << refinement << " --- " << std::endl;
  512.         HybridPoissonProblem<dim> hyb_poisson(degree);
  513.         hyb_poisson.run( refinement, convergence_table, final_time );
  514.     }
  515.  
  516.     convergence_table.set_precision("L2", 3);
  517.     convergence_table.set_scientific("L2", true);
  518.     convergence_table.evaluate_convergence_rates("L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
  519.  
  520.     convergence_table.set_precision("H1", 3);
  521.     convergence_table.set_scientific("H1", true);
  522.     convergence_table.evaluate_convergence_rates("H1", "cells", ConvergenceTable::reduction_rate_log2, dim);
  523.  
  524.     convergence_table.write_text(std::cout);
  525.  
  526.     std::ofstream tex_output("taxas.tex");
  527.     convergence_table.write_tex(tex_output);
  528.  
  529.     std::ofstream data_output("taxas.dat");
  530.     convergence_table.write_text(data_output);
  531.  
  532.     return 0;
  533. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!