deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-74.h
Go to the documentation of this file.
1) const
275 *   {
276 *   using numbers::PI;
277 *   for (unsigned int i = 0; i < values.size(); ++i)
278 *   values[i] =
279 *   std::sin(2. * PI * points[i][0]) * std::sin(2. * PI * points[i][1]);
280 *   }
281 *  
282 *  
283 *  
284 *   template <int dim>
287 *   const unsigned int /*component*/) const
288 *   {
289 *   Tensor<1, dim> return_value;
290 *   using numbers::PI;
291 *   return_value[0] =
292 *   2. * PI * std::cos(2. * PI * point[0]) * std::sin(2. * PI * point[1]);
293 *   return_value[1] =
294 *   2. * PI * std::sin(2. * PI * point[0]) * std::cos(2. * PI * point[1]);
295 *   return return_value;
296 *   }
297 *  
298 *  
299 *  
300 * @endcode
301 *
302 * The corresponding right-hand side of the smooth function:
303 *
304 * @code
305 *   template <int dim>
306 *   class SmoothRightHandSide : public Function<dim>
307 *   {
308 *   public:
310 *   : Function<dim>()
311 *   {}
312 *  
313 *   virtual void value_list(const std::vector<Point<dim>> &points,
314 *   std::vector<double> &values,
315 *   const unsigned int /*component*/) const override;
316 *   };
317 *  
318 *  
319 *  
320 *   template <int dim>
321 *   void
322 *   SmoothRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
323 *   std::vector<double> &values,
324 *   const unsigned int /*component*/) const
325 *   {
326 *   using numbers::PI;
327 *   for (unsigned int i = 0; i < values.size(); ++i)
328 *   values[i] = 8. * PI * PI * std::sin(2. * PI * points[i][0]) *
329 *   std::sin(2. * PI * points[i][1]);
330 *   }
331 *  
332 *  
333 *  
334 * @endcode
335 *
336 * The right-hand side that corresponds to the function
339 *
340 * @code
341 *   template <int dim>
342 *   class SingularRightHandSide : public Function<dim>
343 *   {
344 *   public:
346 *   : Function<dim>()
347 *   {}
348 *  
349 *   virtual void value_list(const std::vector<Point<dim>> &points,
350 *   std::vector<double> &values,
351 *   const unsigned int /*component*/) const override;
352 *  
353 *   private:
355 *   };
356 *  
357 *  
358 *  
359 *   template <int dim>
360 *   void
361 *   SingularRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
362 *   std::vector<double> &values,
363 *   const unsigned int /*component*/) const
364 *   {
365 *   for (unsigned int i = 0; i < values.size(); ++i)
366 *   values[i] = -ref.laplacian(points[i]);
367 *   }
368 *  
369 *  
370 *  
371 * @endcode
372 *
373 *
374 * <a name="step_74-Auxiliaryfunctions"></a>
376 * This function computes the penalty @f$\sigma@f$.
377 *
378 * @code
379 *   double get_penalty_factor(const unsigned int fe_degree,
380 *   const double cell_extent_left,
381 *   const double cell_extent_right)
382 *   {
383 *   const unsigned int degree = std::max(1U, fe_degree);
384 *   return degree * (degree + 1.) * 0.5 *
385 *   (1. / cell_extent_left + 1. / cell_extent_right);
386 *   }
387 *  
388 *  
389 * @endcode
390 *
391 *
392 * <a name="step_74-TheCopyData"></a>
393 * <h3>The CopyData</h3>
394 * In the following, we define "Copy" objects for the MeshWorker::mesh_loop(),
395 * which is essentially the same as @ref step_12 "step-12". Note that the
396 * "Scratch" object is not defined here because we use
397 * MeshWorker::ScratchData<dim> instead. (The use of "Copy" and "Scratch"
398 * objects is extensively explained in the WorkStream namespace documentation.
399 *
400 * @code
401 *   struct CopyDataFace
402 *   {
404 *   std::vector<types::global_dof_index> joint_dof_indices;
405 *   std::array<double, 2> values;
406 *   std::array<unsigned int, 2> cell_indices;
407 *   };
408 *  
409 *  
410 *  
411 *   struct CopyData
412 *   {
415 *   std::vector<types::global_dof_index> local_dof_indices;
416 *   std::vector<CopyDataFace> face_data;
417 *   double value;
418 *   unsigned int cell_index;
419 *  
420 *  
421 *   template <class Iterator>
422 *   void reinit(const Iterator &cell, const unsigned int dofs_per_cell)
423 *   {
424 *   cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
425 *   cell_rhs.reinit(dofs_per_cell);
426 *   local_dof_indices.resize(dofs_per_cell);
427 *   cell->get_dof_indices(local_dof_indices);
428 *   }
429 *   };
430 *  
431 *  
432 *  
433 * @endcode
434 *
435 *
436 * <a name="step_74-TheSIPGLaplaceclass"></a>
437 * <h3>The SIPGLaplace class</h3>
442 * assemble face terms.
443 *
444 * @code
445 *   template <int dim>
446 *   class SIPGLaplace
447 *   {
448 *   public:
450 *   void run();
451 *  
452 *   private:
453 *   void setup_system();
454 *   void assemble_system();
455 *   void solve();
456 *   void refine_grid();
457 *   void output_results(const unsigned int cycle) const;
458 *  
459 *   void compute_errors();
460 *   void compute_error_estimate();
461 *   double compute_energy_norm_error();
462 *  
464 *   const unsigned int degree;
465 *   const QGauss<dim> quadrature;
466 *   const QGauss<dim - 1> face_quadrature;
469 *   const MappingQ1<dim> mapping;
470 *  
471 *   using ScratchData = MeshWorker::ScratchData<dim>;
472 *  
473 *   const FE_DGQ<dim> fe;
474 *   DoFHandler<dim> dof_handler;
475 *  
476 *   SparsityPattern sparsity_pattern;
477 *   SparseMatrix<double> system_matrix;
478 *   Vector<double> solution;
480 *  
481 * @endcode
482 *
483 * The remainder of the class's members are used for the following:
484 * - Vectors to store error estimator square and energy norm square per
485 * cell.
486 * - Print convergence rate and errors on the screen.
487 * - The fiffusion coefficient @f$\nu@f$ is set to 1.
488 * - Members that store information about the test case to be computed.
489 *
490 * @code
491 *   Vector<double> estimated_error_square_per_cell;
492 *   Vector<double> energy_norm_square_per_cell;
493 *  
494 *   ConvergenceTable convergence_table;
495 *  
496 *   const double diffusion_coefficient = 1.;
497 *  
498 *   const TestCase test_case;
499 *   std::unique_ptr<const Function<dim>> exact_solution;
500 *   std::unique_ptr<const Function<dim>> rhs_function;
501 *   };
502 *  
503 * @endcode
504 *
505 * The constructor here takes the test case as input and then
506 * determines the correct solution and right-hand side classes. The
507 * remaining member variables are initialized in the obvious way.
508 *
509 * @code
510 *   template <int dim>
511 *   SIPGLaplace<dim>::SIPGLaplace(const TestCase &test_case)
512 *   : degree(3)
513 *   , quadrature(degree + 1)
514 *   , face_quadrature(degree + 1)
515 *   , quadrature_overintegration(degree + 2)
516 *   , face_quadrature_overintegration(degree + 2)
517 *   , mapping()
518 *   , fe(degree)
519 *   , dof_handler(triangulation)
520 *   , test_case(test_case)
521 *   {
522 *   if (test_case == TestCase::convergence_rate)
523 *   {
524 *   exact_solution = std::make_unique<const SmoothSolution<dim>>();
525 *   rhs_function = std::make_unique<const SmoothRightHandSide<dim>>();
526 *   }
527 *  
528 *   else if (test_case == TestCase::l_singularity)
529 *   {
530 *   exact_solution =
531 *   std::make_unique<const Functions::LSingularityFunction>();
532 *   rhs_function = std::make_unique<const SingularRightHandSide<dim>>();
533 *   }
534 *   else
535 *   AssertThrow(false, ExcNotImplemented());
536 *   }
537 *  
538 *  
539 *  
540 *   template <int dim>
541 *   void SIPGLaplace<dim>::setup_system()
542 *   {
543 *   dof_handler.distribute_dofs(fe);
544 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
545 *   DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
546 *   sparsity_pattern.copy_from(dsp);
547 *  
548 *   system_matrix.reinit(sparsity_pattern);
549 *   solution.reinit(dof_handler.n_dofs());
550 *   system_rhs.reinit(dof_handler.n_dofs());
551 *   }
552 *  
553 *  
554 *  
555 * @endcode
556 *
557 *
558 * <a name="step_74-Theassemble_systemfunction"></a>
559 * <h3>The assemble_system function</h3>
560 * The assemble function here is similar to that in @ref step_12 "step-12" and @ref step_47 "step-47".
561 * Different from assembling by hand, we just need to focus
562 * on assembling on each cell, each boundary face, and each
563 * interior face. The loops over cells and faces are handled
564 * automatically by MeshWorker::mesh_loop().
565 *
566
567 *
568 * The function starts by defining a local (lambda) function that is
569 * used to integrate the cell terms:
570 *
571 * @code
572 *   template <int dim>
573 *   void SIPGLaplace<dim>::assemble_system()
574 *   {
575 *   const auto cell_worker =
576 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
577 *   ScratchData &scratch_data,
578 *   CopyData &copy_data) {
579 *   const FEValues<dim> &fe_v = scratch_data.reinit(cell);
580 *   const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
581 *   copy_data.reinit(cell, dofs_per_cell);
582 *  
583 *   const std::vector<Point<dim>> &q_points =
584 *   scratch_data.get_quadrature_points();
585 *   const unsigned int n_q_points = q_points.size();
586 *   const std::vector<double> &JxW = scratch_data.get_JxW_values();
587 *  
588 *   std::vector<double> rhs(n_q_points);
589 *   rhs_function->value_list(q_points, rhs);
590 *  
591 *   for (unsigned int point = 0; point < n_q_points; ++point)
592 *   for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
593 *   {
594 *   for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
595 *   copy_data.cell_matrix(i, j) +=
596 *   diffusion_coefficient * // nu
597 *   fe_v.shape_grad(i, point) * // grad v_h
598 *   fe_v.shape_grad(j, point) * // grad u_h
599 *   JxW[point]; // dx
600 *  
601 *   copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
602 *   rhs[point] * // f
603 *   JxW[point]; // dx
604 *   }
605 *   };
606 *  
607 * @endcode
608 *
609 * Next, we need a function that assembles face integrals on the boundary:
610 *
611 * @code
612 *   const auto boundary_worker =
613 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
614 *   const unsigned int &face_no,
615 *   ScratchData &scratch_data,
616 *   CopyData &copy_data) {
617 *   const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
618 *  
619 *   const std::vector<Point<dim>> &q_points =
620 *   scratch_data.get_quadrature_points();
621 *   const unsigned int n_q_points = q_points.size();
622 *   const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
623 *  
624 *   const std::vector<double> &JxW = scratch_data.get_JxW_values();
625 *   const std::vector<Tensor<1, dim>> &normals =
626 *   scratch_data.get_normal_vectors();
627 *  
628 *   std::vector<double> g(n_q_points);
629 *   exact_solution->value_list(q_points, g);
630 *  
631 *   const double extent1 = cell->measure() / cell->face(face_no)->measure();
632 *   const double penalty = get_penalty_factor(degree, extent1, extent1);
633 *  
634 *   for (unsigned int point = 0; point < n_q_points; ++point)
635 *   {
636 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
637 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
638 *   copy_data.cell_matrix(i, j) +=
639 *   (-diffusion_coefficient * // - nu
640 *   fe_fv.shape_value(i, point) * // v_h
641 *   (fe_fv.shape_grad(j, point) * // (grad u_h .
642 *   normals[point]) // n)
643 *  
644 *   - diffusion_coefficient * // - nu
645 *   (fe_fv.shape_grad(i, point) * // (grad v_h .
646 *   normals[point]) * // n)
647 *   fe_fv.shape_value(j, point) // u_h
648 *  
649 *   + diffusion_coefficient * penalty * // + nu sigma
650 *   fe_fv.shape_value(i, point) * // v_h
651 *   fe_fv.shape_value(j, point) // u_h
652 *  
653 *   ) *
654 *   JxW[point]; // dx
655 *  
656 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
657 *   copy_data.cell_rhs(i) +=
658 *   (-diffusion_coefficient * // - nu
659 *   (fe_fv.shape_grad(i, point) * // (grad v_h .
660 *   normals[point]) * // n)
661 *   g[point] // g
662 *  
663 *  
664 *   + diffusion_coefficient * penalty * // + nu sigma
665 *   fe_fv.shape_value(i, point) * g[point] // v_h g
666 *  
667 *   ) *
668 *   JxW[point]; // dx
669 *   }
670 *   };
671 *  
672 * @endcode
673 *
674 * Finally, a function that assembles face integrals on interior
675 * faces. To reinitialize FEInterfaceValues, we need to pass
676 * cells, face and subface indices (for adaptive refinement) to
677 * the reinit() function of FEInterfaceValues:
678 *
679 * @code
680 *   const auto face_worker =
681 *   [&](const typename DoFHandler<dim>::cell_iterator &cell,
682 *   const unsigned int &f,
683 *   const unsigned int &sf,
684 *   const typename DoFHandler<dim>::cell_iterator &ncell,
685 *   const unsigned int &nf,
686 *   const unsigned int &nsf,
687 *   ScratchData &scratch_data,
688 *   CopyData &copy_data) {
689 *   const FEInterfaceValues<dim> &fe_iv =
690 *   scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
691 *  
692 *   copy_data.face_data.emplace_back();
693 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
694 *   const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs();
695 *   copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
696 *   copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
697 *  
698 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
699 *   const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
700 *  
701 *   const double extent1 = cell->measure() / cell->face(f)->measure();
702 *   const double extent2 = ncell->measure() / ncell->face(nf)->measure();
703 *   const double penalty = get_penalty_factor(degree, extent1, extent2);
704 *  
705 *   for (const unsigned int point : fe_iv.quadrature_point_indices())
706 *   {
707 *   for (const unsigned int i : fe_iv.dof_indices())
708 *   for (const unsigned int j : fe_iv.dof_indices())
709 *   copy_data_face.cell_matrix(i, j) +=
710 *   (-diffusion_coefficient * // - nu
711 *   fe_iv.jump_in_shape_values(i, point) * // [v_h]
712 *   (fe_iv.average_of_shape_gradients(j,
713 *   point) * // ({grad u_h} .
714 *   normals[point]) // n)
715 *  
716 *   - diffusion_coefficient * // - nu
717 *   (fe_iv.average_of_shape_gradients(i,
718 *   point) * // (grad v_h .
719 *   normals[point]) * // n)
720 *   fe_iv.jump_in_shape_values(j, point) // [u_h]
721 *  
722 *   + diffusion_coefficient * penalty * // + nu sigma
723 *   fe_iv.jump_in_shape_values(i, point) * // [v_h]
724 *   fe_iv.jump_in_shape_values(j, point) // [u_h]
725 *  
726 *   ) *
727 *   JxW[point]; // dx
728 *   }
729 *   };
730 *  
731 * @endcode
732 *
733 * The following lambda function will then copy data into the
734 * global matrix and right-hand side. Though there are no hanging
735 * node constraints in DG discretization, we define an empty
736 * AffineConstraints object that allows us to use the
737 * AffineConstraints::distribute_local_to_global() functionality.
738 *
739 * @code
740 *   AffineConstraints<double> constraints;
741 *   constraints.close();
742 *   const auto copier = [&](const CopyData &c) {
743 *   constraints.distribute_local_to_global(c.cell_matrix,
744 *   c.cell_rhs,
745 *   c.local_dof_indices,
746 *   system_matrix,
747 *   system_rhs);
748 *  
749 * @endcode
750 *
751 * Copy data from interior face assembly to the global matrix.
752 *
753 * @code
754 *   for (const CopyDataFace &cdf : c.face_data)
755 *   {
756 *   constraints.distribute_local_to_global(cdf.cell_matrix,
757 *   cdf.joint_dof_indices,
758 *   system_matrix);
759 *   }
760 *   };
761 *  
762 *  
763 * @endcode
764 *
765 * With the assembly functions defined, we can now create
766 * ScratchData and CopyData objects, and pass them together with
767 * the lambda functions above to MeshWorker::mesh_loop(). In
768 * addition, we need to specify that we want to assemble on
769 * interior faces exactly once.
770 *
771 * @code
772 *   const UpdateFlags cell_flags = update_values | update_gradients |
773 *   update_quadrature_points | update_JxW_values;
774 *   const UpdateFlags face_flags = update_values | update_gradients |
775 *   update_quadrature_points |
776 *   update_normal_vectors | update_JxW_values;
777 *  
778 *   ScratchData scratch_data(
779 *   mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
780 *   CopyData copy_data;
781 *  
782 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
783 *   dof_handler.end(),
784 *   cell_worker,
785 *   copier,
786 *   scratch_data,
787 *   copy_data,
788 *   MeshWorker::assemble_own_cells |
789 *   MeshWorker::assemble_boundary_faces |
790 *   MeshWorker::assemble_own_interior_faces_once,
791 *   boundary_worker,
792 *   face_worker);
793 *   }
794 *  
795 *  
796 *  
797 * @endcode
798 *
799 *
800 * <a name="step_74-Thesolveandoutput_resultsfunction"></a>
801 * <h3>The solve() and output_results() function</h3>
802 * The following two functions are entirely standard and without difficulty.
803 *
804 * @code
805 *   template <int dim>
806 *   void SIPGLaplace<dim>::solve()
807 *   {
808 *   SparseDirectUMFPACK A_direct;
809 *   A_direct.initialize(system_matrix);
810 *   A_direct.vmult(solution, system_rhs);
811 *   }
812 *  
813 *  
814 *  
815 *   template <int dim>
816 *   void SIPGLaplace<dim>::output_results(const unsigned int cycle) const
817 *   {
818 *   const std::string filename = "sol_Q" + Utilities::int_to_string(degree, 1) +
819 *   "-" + Utilities::int_to_string(cycle, 2) +
820 *   ".vtu";
821 *   std::ofstream output(filename);
822 *  
823 *   DataOut<dim> data_out;
824 *   data_out.attach_dof_handler(dof_handler);
825 *   data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
826 *   data_out.build_patches(mapping);
827 *   data_out.write_vtu(output);
828 *   }
829 *  
830 *  
831 * @endcode
832 *
833 *
834 * <a name="step_74-Thecompute_error_estimatefunction"></a>
835 * <h3>The compute_error_estimate() function</h3>
836 * The assembly of the error estimator here is quite similar to
837 * that of the global matrix and right-had side and can be handled
838 * by the MeshWorker::mesh_loop() framework. To understand what
839 * each of the local (lambda) functions is doing, recall first that
840 * the local cell residual is defined as
841 * @f$h_K^2 \left\| f + \nu \Delta u_h \right\|_K^2@f$:
842 *
843 * @code
844 *   template <int dim>
845 *   void SIPGLaplace<dim>::compute_error_estimate()
846 *   {
847 *   const auto cell_worker =
848 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
849 *   ScratchData &scratch_data,
850 *   CopyData &copy_data) {
851 *   const FEValues<dim> &fe_v = scratch_data.reinit(cell);
852 *  
853 *   copy_data.cell_index = cell->active_cell_index();
854 *  
855 *   const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
856 *   const unsigned int n_q_points = q_points.size();
857 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
858 *  
859 *   std::vector<Tensor<2, dim>> hessians(n_q_points);
860 *   fe_v.get_function_hessians(solution, hessians);
861 *  
862 *   std::vector<double> rhs(n_q_points);
863 *   rhs_function->value_list(q_points, rhs);
864 *  
865 *   const double hk = cell->diameter();
866 *   double residual_norm_square = 0;
867 *  
868 *   for (unsigned int point = 0; point < n_q_points; ++point)
869 *   {
870 *   const double residual =
871 *   rhs[point] + diffusion_coefficient * trace(hessians[point]);
872 *   residual_norm_square += residual * residual * JxW[point];
873 *   }
874 *   copy_data.value = hk * hk * residual_norm_square;
875 *   };
876 *  
877 * @endcode
878 *
879 * Next compute boundary terms @f$\sum_{f\in \partial K \cap \partial \Omega}
880 * \sigma \left\| [ u_h-g_D ] \right\|_f^2 @f$:
881 *
882 * @code
883 *   const auto boundary_worker =
884 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
885 *   const unsigned int &face_no,
886 *   ScratchData &scratch_data,
887 *   CopyData &copy_data) {
888 *   const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
889 *  
890 *   const std::vector<Point<dim>> &q_points = fe_fv.get_quadrature_points();
891 *   const unsigned n_q_points = q_points.size();
892 *  
893 *   const std::vector<double> &JxW = fe_fv.get_JxW_values();
894 *  
895 *   std::vector<double> g(n_q_points);
896 *   exact_solution->value_list(q_points, g);
897 *  
898 *   std::vector<double> sol_u(n_q_points);
899 *   fe_fv.get_function_values(solution, sol_u);
900 *  
901 *   const double extent1 = cell->measure() / cell->face(face_no)->measure();
902 *   const double penalty = get_penalty_factor(degree, extent1, extent1);
903 *  
904 *   double difference_norm_square = 0.;
905 *   for (unsigned int point = 0; point < q_points.size(); ++point)
906 *   {
907 *   const double diff = (g[point] - sol_u[point]);
908 *   difference_norm_square += diff * diff * JxW[point];
909 *   }
910 *   copy_data.value += penalty * difference_norm_square;
911 *   };
912 *  
913 * @endcode
914 *
915 * And finally interior face terms @f$\sum_{f\in \partial K}\lbrace \sigma
916 * \left\| [u_h] \right\|_f^2 + h_f \left\| [\nu \nabla u_h \cdot
917 * \mathbf n ] \right\|_f^2 \rbrace@f$:
918 *
919 * @code
920 *   const auto face_worker =
921 *   [&](const typename DoFHandler<dim>::cell_iterator &cell,
922 *   const unsigned int &f,
923 *   const unsigned int &sf,
924 *   const typename DoFHandler<dim>::cell_iterator &ncell,
925 *   const unsigned int &nf,
926 *   const unsigned int &nsf,
927 *   ScratchData &scratch_data,
928 *   CopyData &copy_data) {
929 *   const FEInterfaceValues<dim> &fe_iv =
930 *   scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
931 *  
932 *   copy_data.face_data.emplace_back();
933 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
934 *  
935 *   copy_data_face.cell_indices[0] = cell->active_cell_index();
936 *   copy_data_face.cell_indices[1] = ncell->active_cell_index();
937 *  
938 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
939 *   const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
940 *  
941 *   const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
942 *   const unsigned int n_q_points = q_points.size();
943 *  
944 *   std::vector<double> jump(n_q_points);
945 *   fe_iv.get_jump_in_function_values(solution, jump);
946 *  
947 *   std::vector<Tensor<1, dim>> grad_jump(n_q_points);
948 *   fe_iv.get_jump_in_function_gradients(solution, grad_jump);
949 *  
950 *   const double h = cell->face(f)->diameter();
951 *  
952 *   const double extent1 = cell->measure() / cell->face(f)->measure();
953 *   const double extent2 = ncell->measure() / ncell->face(nf)->measure();
954 *   const double penalty = get_penalty_factor(degree, extent1, extent2);
955 *  
956 *   double flux_jump_square = 0;
957 *   double u_jump_square = 0;
958 *   for (unsigned int point = 0; point < n_q_points; ++point)
959 *   {
960 *   u_jump_square += jump[point] * jump[point] * JxW[point];
961 *   const double flux_jump = grad_jump[point] * normals[point];
962 *   flux_jump_square +=
963 *   diffusion_coefficient * flux_jump * flux_jump * JxW[point];
964 *   }
965 *   copy_data_face.values[0] =
966 *   0.5 * h * (flux_jump_square + penalty * u_jump_square);
967 *   copy_data_face.values[1] = copy_data_face.values[0];
968 *   };
969 *  
970 * @endcode
971 *
972 * Having computed local contributions for each cell, we still
973 * need a way to copy these into the global vector that will hold
974 * the error estimators for all cells:
975 *
976 * @code
977 *   const auto copier = [&](const CopyData &copy_data) {
978 *   if (copy_data.cell_index != numbers::invalid_unsigned_int)
979 *   estimated_error_square_per_cell[copy_data.cell_index] +=
980 *   copy_data.value;
981 *   for (const CopyDataFace &cdf : copy_data.face_data)
982 *   for (unsigned int j = 0; j < 2; ++j)
983 *   estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
984 *   };
985 *  
986 * @endcode
987 *
988 * After all of this set-up, let's do the actual work: We resize
989 * the vector into which the results will be written, and then
990 * drive the whole process using the MeshWorker::mesh_loop()
991 * function.
992 *
993 * @code
994 *   estimated_error_square_per_cell.reinit(triangulation.n_active_cells());
995 *  
996 *   const UpdateFlags cell_flags =
998 *   const UpdateFlags face_flags = update_values | update_gradients |
1001 *  
1002 *   ScratchData scratch_data(
1003 *   mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
1004 *  
1005 *   CopyData copy_data;
1006 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1007 *   dof_handler.end(),
1008 *   cell_worker,
1009 *   copier,
1010 *   scratch_data,
1011 *   copy_data,
1016 *   face_worker);
1017 *   }
1018 *  
1019 * @endcode
1020 *
1021 *
1022 * <a name="step_74-Thecompute_energy_norm_errorfunction"></a>
1023 * <h3>The compute_energy_norm_error() function</h3>
1024 * Next, we evaluate the accuracy in terms of the energy norm.
1026 * Here we compute the square of the energy norm defined by
1027 * @f[
1028 * \|u \|_{1,h}^2 = \sum_{K \in \Gamma_h} \nu\|\nabla u \|_K^2 +
1029 * \sum_{f \in F_i} \sigma \| [ u ] \|_f^2 +
1030 * \sum_{f \in F_b} \sigma \|u\|_f^2.
1031 * @f]
1033 * @f[
1034 * \|u -u_h \|_{1,h}^2 = \sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2
1035 * + \sum_{f \in F_i} \sigma \|[ u_h ] \|_f^2 + \sum_{f \in F_b}\sigma
1036 * \|u_h-g_D\|_f^2.
1037 * @f]
1038 *
1039 * @code
1040 *   template <int dim>
1042 *   {
1043 *   energy_norm_square_per_cell.reinit(triangulation.n_active_cells());
1044 *  
1045 * @endcode
1046 *
1047 * Assemble @f$\sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2 @f$.
1048 *
1049 * @code
1050 *   const auto cell_worker =
1051 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
1052 *   ScratchData &scratch_data,
1053 *   CopyData &copy_data) {
1054 *   const FEValues<dim> &fe_v = scratch_data.reinit(cell);
1055 *  
1056 *   copy_data.cell_index = cell->active_cell_index();
1057 *  
1058 *   const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
1059 *   const unsigned int n_q_points = q_points.size();
1060 *   const std::vector<double> &JxW = fe_v.get_JxW_values();
1061 *  
1062 *   std::vector<Tensor<1, dim>> grad_u(n_q_points);
1063 *   fe_v.get_function_gradients(solution, grad_u);
1064 *  
1065 *   std::vector<Tensor<1, dim>> grad_exact(n_q_points);
1066 *   exact_solution->gradient_list(q_points, grad_exact);
1067 *  
1068 *   double norm_square = 0;
1069 *   for (unsigned int point = 0; point < n_q_points; ++point)
1070 *   {
1071 *   norm_square +=
1072 *   (grad_u[point] - grad_exact[point]).norm_square() * JxW[point];
1073 *   }
1074 *   copy_data.value = diffusion_coefficient * norm_square;
1075 *   };
1076 *  
1077 * @endcode
1078 *
1080 *
1081 * @code
1082 *   const auto boundary_worker =
1083 *   [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
1084 *   const unsigned int &face_no,
1085 *   ScratchData &scratch_data,
1086 *   CopyData &copy_data) {
1087 *   const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
1088 *  
1089 *   const std::vector<Point<dim>> &q_points = fe_fv.get_quadrature_points();
1090 *   const unsigned n_q_points = q_points.size();
1091 *  
1092 *   const std::vector<double> &JxW = fe_fv.get_JxW_values();
1093 *  
1094 *   std::vector<double> g(n_q_points);
1095 *   exact_solution->value_list(q_points, g);
1096 *  
1097 *   std::vector<double> sol_u(n_q_points);
1098 *   fe_fv.get_function_values(solution, sol_u);
1099 *  
1100 *   const double extent1 = cell->measure() / cell->face(face_no)->measure();
1101 *   const double penalty = get_penalty_factor(degree, extent1, extent1);
1102 *  
1103 *   double difference_norm_square = 0.;
1104 *   for (unsigned int point = 0; point < q_points.size(); ++point)
1105 *   {
1106 *   const double diff = (g[point] - sol_u[point]);
1107 *   difference_norm_square += diff * diff * JxW[point];
1108 *   }
1109 *   copy_data.value += penalty * difference_norm_square;
1110 *   };
1111 *  
1112 * @endcode
1113 *
1114 * Assemble @f$\sum_{f \in F_i} \sigma \| [ u_h ] \|_f^2@f$.
1115 *
1116 * @code
1117 *   const auto face_worker =
1118 *   [&](const typename DoFHandler<dim>::cell_iterator &cell,
1119 *   const unsigned int &f,
1120 *   const unsigned int &sf,
1121 *   const typename DoFHandler<dim>::cell_iterator &ncell,
1122 *   const unsigned int &nf,
1123 *   const unsigned int &nsf,
1124 *   ScratchData &scratch_data,
1125 *   CopyData &copy_data) {
1127 *   scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
1128 *  
1129 *   copy_data.face_data.emplace_back();
1130 *   CopyDataFace &copy_data_face = copy_data.face_data.back();
1131 *  
1132 *   copy_data_face.cell_indices[0] = cell->active_cell_index();
1133 *   copy_data_face.cell_indices[1] = ncell->active_cell_index();
1134 *  
1135 *   const std::vector<double> &JxW = fe_iv.get_JxW_values();
1136 *  
1137 *   const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
1138 *   const unsigned int n_q_points = q_points.size();
1139 *  
1140 *   std::vector<double> jump(n_q_points);
1141 *   fe_iv.get_jump_in_function_values(solution, jump);
1142 *  
1143 *   const double extent1 = cell->measure() / cell->face(f)->measure();
1144 *   const double extent2 = ncell->measure() / ncell->face(nf)->measure();
1145 *   const double penalty = get_penalty_factor(degree, extent1, extent2);
1146 *  
1147 *   double u_jump_square = 0;
1148 *   for (unsigned int point = 0; point < n_q_points; ++point)
1149 *   {
1150 *   u_jump_square += jump[point] * jump[point] * JxW[point];
1151 *   }
1154 *   };
1155 *  
1156 *   const auto copier = [&](const CopyData &copy_data) {
1157 *   if (copy_data.cell_index != numbers::invalid_unsigned_int)
1158 *   energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
1159 *   for (const CopyDataFace &cdf : copy_data.face_data)
1160 *   for (unsigned int j = 0; j < 2; ++j)
1161 *   energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
1162 *   };
1163 *  
1164 *   const UpdateFlags cell_flags =
1166 *   UpdateFlags face_flags =
1168 *  
1169 *   const ScratchData scratch_data(mapping,
1170 *   fe,
1172 *   cell_flags,
1174 *   face_flags);
1175 *  
1176 *   CopyData copy_data;
1177 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1178 *   dof_handler.end(),
1179 *   cell_worker,
1180 *   copier,
1181 *   scratch_data,
1182 *   copy_data,
1187 *   face_worker);
1188 *   const double energy_error =
1190 *   return energy_error;
1191 *   }
1192 *  
1193 *  
1194 *  
1195 * @endcode
1196 *
1197 *
1198 * <a name="step_74-Therefine_gridfunction"></a>
1199 * <h3>The refine_grid() function</h3>
1200 *
1201 * @code
1202 *   template <int dim>
1203 *   void SIPGLaplace<dim>::refine_grid()
1204 *   {
1205 *   const double refinement_fraction = 0.1;
1206 *  
1209 *  
1210 *   triangulation.execute_coarsening_and_refinement();
1211 *   }
1212 *  
1213 *  
1214 *  
1215 * @endcode
1216 *
1217 *
1218 * <a name="step_74-Thecompute_errorsfunction"></a>
1219 * <h3>The compute_errors() function</h3>
1220 * We compute three errors in the @f$L_2@f$ norm, @f$H_1@f$ seminorm, and
1223 * with mesh refinement and which can be output in one step at the
1224 * end of the program.
1225 *
1226 * @code
1227 *   template <int dim>
1228 *   void SIPGLaplace<dim>::compute_errors()
1229 *   {
1230 *   double L2_error, H1_error, energy_error;
1231 *  
1232 *   {
1235 *   dof_handler,
1236 *   solution,
1237 *   *(exact_solution.get()),
1241 *  
1245 *   convergence_table.add_value("L2", L2_error);
1246 *   }
1247 *  
1248 *   {
1251 *   dof_handler,
1252 *   solution,
1253 *   *(exact_solution.get()),
1257 *  
1261 *   convergence_table.add_value("H1", H1_error);
1262 *   }
1263 *  
1264 *   {
1266 *   convergence_table.add_value("Energy", energy_error);
1267 *   }
1268 *  
1269 *   std::cout << " Error in the L2 norm : " << L2_error << std::endl
1270 *   << " Error in the H1 seminorm : " << H1_error << std::endl
1271 *   << " Error in the energy norm : " << energy_error
1272 *   << std::endl;
1273 *   }
1274 *  
1275 *  
1276 *  
1277 * @endcode
1278 *
1279 *
1280 * <a name="step_74-Therunfunction"></a>
1281 * <h3>The run() function</h3>
1282 *
1283 * @code
1284 *   template <int dim>
1285 *   void SIPGLaplace<dim>::run()
1286 *   {
1287 *   const unsigned int max_cycle =
1288 *   (test_case == TestCase::convergence_rate ? 6 : 20);
1289 *   for (unsigned int cycle = 0; cycle < max_cycle; ++cycle)
1290 *   {
1291 *   std::cout << "Cycle " << cycle << std::endl;
1292 *  
1293 *   switch (test_case)
1294 *   {
1296 *   {
1297 *   if (cycle == 0)
1298 *   {
1300 *  
1301 *   triangulation.refine_global(2);
1302 *   }
1303 *   else
1304 *   {
1305 *   triangulation.refine_global(1);
1306 *   }
1307 *   break;
1308 *   }
1309 *  
1311 *   {
1312 *   if (cycle == 0)
1313 *   {
1315 *   triangulation.refine_global(3);
1316 *   }
1317 *   else
1318 *   {
1319 *   refine_grid();
1320 *   }
1321 *   break;
1322 *   }
1323 *  
1324 *   default:
1325 *   {
1327 *   }
1328 *   }
1329 *  
1330 *   std::cout << " Number of active cells : "
1331 *   << triangulation.n_active_cells() << std::endl;
1332 *   setup_system();
1333 *  
1334 *   std::cout << " Number of degrees of freedom : " << dof_handler.n_dofs()
1335 *   << std::endl;
1336 *  
1337 *   assemble_system();
1338 *   solve();
1339 *   output_results(cycle);
1340 *   {
1341 *   convergence_table.add_value("cycle", cycle);
1342 *   convergence_table.add_value("cells", triangulation.n_active_cells());
1343 *   convergence_table.add_value("dofs", dof_handler.n_dofs());
1344 *   }
1345 *   compute_errors();
1346 *  
1348 *   {
1350 *   std::cout << " Estimated error : "
1352 *   << std::endl;
1353 *  
1354 *   convergence_table.add_value(
1355 *   "Estimator",
1357 *   }
1358 *   std::cout << std::endl;
1359 *   }
1360 *  
1361 * @endcode
1362 *
1364 * table how to format its data and output it to screen:
1365 *
1366 * @code
1367 *   convergence_table.set_precision("L2", 3);
1368 *   convergence_table.set_precision("H1", 3);
1369 *   convergence_table.set_precision("Energy", 3);
1370 *  
1371 *   convergence_table.set_scientific("L2", true);
1372 *   convergence_table.set_scientific("H1", true);
1373 *   convergence_table.set_scientific("Energy", true);
1374 *  
1376 *   {
1377 *   convergence_table.evaluate_convergence_rates(
1379 *   convergence_table.evaluate_convergence_rates(
1381 *   }
1383 *   {
1384 *   convergence_table.set_precision("Estimator", 3);
1385 *   convergence_table.set_scientific("Estimator", true);
1386 *   }
1387 *  
1388 *   std::cout << "degree = " << degree << std::endl;
1389 *   convergence_table.write_text(
1391 *   }
1392 *   } // namespace Step74
1393 *  
1394 *  
1395 *  
1396 * @endcode
1397 *
1398 *
1399 * <a name="step_74-Themainfunction"></a>
1400 * <h3>The main() function</h3>
1403 *
1404 * @code
1405 *   int main()
1406 *   {
1407 *   try
1408 *   {
1409 *   using namespace dealii;
1410 *   using namespace Step74;
1411 *  
1413 *  
1415 *   problem.run();
1416 *   }
1417 *   catch (std::exception &exc)
1418 *   {
1419 *   std::cerr << std::endl
1420 *   << std::endl
1421 *   << "----------------------------------------------------"
1422 *   << std::endl;
1423 *   std::cerr << "Exception on processing: " << std::endl
1424 *   << exc.what() << std::endl
1425 *   << "Aborting!" << std::endl
1426 *   << "----------------------------------------------------"
1427 *   << std::endl;
1428 *   return 1;
1429 *   }
1430 *   catch (...)
1431 *   {
1432 *   std::cerr << std::endl
1433 *   << std::endl
1434 *   << "----------------------------------------------------"
1435 *   << std::endl;
1436 *   std::cerr << "Unknown exception!" << std::endl
1437 *   << "Aborting!" << std::endl
1438 *   << "----------------------------------------------------"
1439 *   << std::endl;
1440 *   return 1;
1441 *   };
1442 *  
1443 *   return 0;
1444 *   }
1445 * @endcode
1446<a name="step_74-Results"></a><h1>Results</h1>
1447
1448
1449The output of this program consist of the console output and
1451
1453@code
1454Cycle 0
1455 Number of active cells : 16
1456 Number of degrees of freedom : 256
1457 Error in the L2 norm : 0.00193285
1458 Error in the H1 seminorm : 0.106087
1459 Error in the energy norm : 0.150625
1460
1461Cycle 1
1462 Number of active cells : 64
1463 Number of degrees of freedom : 1024
1464 Error in the L2 norm : 9.60497e-05
1465 Error in the H1 seminorm : 0.0089954
1466 Error in the energy norm : 0.0113265
1467
1468Cycle 2
1469.
1470.
1471.
1472@endcode
1473
1474When using the smooth case with polynomial degree 3, the convergence
1475table will look like this:
1476<table align="center" class="doxtable">
1477 <tr>
1478 <th>cycle</th>
1479 <th>n_cells</th>
1480 <th>n_dofs</th>
1481 <th>L2 </th>
1482 <th>rate</th>
1483 <th>H1</th>
1484 <th>rate</th>
1485 <th>Energy</th>
1486 </tr>
1487 <tr>
1488 <td align="center">0</td>
1489 <td align="right">16</td>
1490 <td align="right">256</td>
1491 <td align="center">1.933e-03</td>
1492 <td>&nbsp;</td>
1493 <td align="center">1.061e-01</td>
1494 <td>&nbsp;</td>
1495 <td align="center">1.506e-01</td>
1496 </tr>
1497 <tr>
1498 <td align="center">1</td>
1499 <td align="right">64</td>
1500 <td align="right">1024</td>
1501 <td align="center">9.605e-05</td>
1502 <td align="center">4.33</td>
1503 <td align="center">8.995e-03</td>
1504 <td align="center">3.56</td>
1505 <td align="center">1.133e-02</td>
1506 </tr>
1507 <tr>
1508 <td align="center">2</td>
1509 <td align="right">256</td>
1510 <td align="right">4096</td>
1511 <td align="center">5.606e-06</td>
1512 <td align="center">4.10</td>
1513 <td align="center">9.018e-04</td>
1514 <td align="center">3.32</td>
1515 <td align="center">9.736e-04</td>
1516 </tr>
1517 <tr>
1518 <td align="center">3</td>
1519 <td align="right">1024</td>
1520 <td align="right">16384</td>
1521 <td align="center">3.484e-07</td>
1522 <td align="center">4.01</td>
1523 <td align="center">1.071e-04</td>
1524 <td align="center">3.07</td>
1525 <td align="center">1.088e-04</td>
1526 </tr>
1527 <tr>
1528 <td align="center">4</td>
1529 <td align="right">4096</td>
1530 <td align="right">65536</td>
1531 <td align="center">2.179e-08</td>
1532 <td align="center">4.00</td>
1533 <td align="center">1.327e-05</td>
1534 <td align="center">3.01</td>
1535 <td align="center">1.331e-05</td>
1536 </tr>
1537 <tr>
1538 <td align="center">5</td>
1539 <td align="right">16384</td>
1540 <td align="right">262144</td>
1541 <td align="center">1.363e-09</td>
1542 <td align="center">4.00</td>
1543 <td align="center">1.656e-06</td>
1544 <td align="center">3.00</td>
1545 <td align="center">1.657e-06</td>
1546 </tr>
1547</table>
1548
1549Theoretically, for polynomial degree @f$p@f$, the order of convergence in @f$L_2@f$
1551results are in good agreement with theory.
1552
1554@code
1555Cycle 0
1556 Number of active cells : 192
1557 Number of degrees of freedom : 3072
1558 Error in the L2 norm : 0.000323585
1559 Error in the H1 seminorm : 0.0296202
1560 Error in the energy norm : 0.0420478
1561 Estimated error : 0.136067
1562
1563Cycle 1
1564 Number of active cells : 249
1565 Number of degrees of freedom : 3984
1566 Error in the L2 norm : 0.000114739
1567 Error in the H1 seminorm : 0.0186571
1568 Error in the energy norm : 0.0264879
1569 Estimated error : 0.0857186
1570
1571Cycle 2
1572.
1573.
1574.
1575@endcode
1576
1578the number of degrees of freedom for this test case on the L-shaped
1579domain. In order to interpret it, let @f$n@f$ be the number of degrees of
1581@f$1/\sqrt{n}@f$ in 2D. Combining the theoretical results in the previous case,
1582we see that if the solution is sufficiently smooth,
1583we can expect the error in the @f$L_2@f$ norm to be of order @f$O(n^{-\frac{p+1}{2}})@f$
1584and in @f$H^1@f$ seminorm to be @f$O(n^{-\frac{p}{2}})@f$. It is not a priori
1585clear that one would get the same kind of behavior as a function of
1587test case, but one can certainly hope. Indeed, from the figure, we see
1589the kinds of hoped-for results:
1590
1591<img width="600px" src="https://www.dealii.org/images/steps/developer/step-74.log-log-plot.png" alt="">
1592
1594at almost the same rate as the errors in the energy norm and @f$H^1@f$ seminorm,
1597
1599large-scale solver in terms of computing time with matrix-free solution techniques.
1601because the multigrid interface matrices are not as easily determined,
1603 *
1604 *
1605<a name="step_74-PlainProg"></a>
1606<h1> The plain program</h1>
1607@include "step-74.cc"
1608*/
numbers::NumberTraits< Number >::real_type norm() const
constexpr void clear()
friend class Tensor
Definition tensor.h:865
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3020
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:851
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
unsigned int cell_index
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:281
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
Definition mpi.cc:740
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
@ matrix
Contents is actually a matrix.
constexpr char L
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition l2.h:159
@ assemble_boundary_faces
@ assemble_own_interior_faces_once
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
double compute_global_error(const Triangulation< dim, spacedim > &tria, const InVector &cellwise_error, const NormType &norm, const double exponent=2.)
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const ReadVector< Number > &fe_function, const Function< spacedim, Number > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14889
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr double PI
Definition numbers.h:241
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation