528 *
double return_value = 0.0;
532 *
return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) +
533 *
+24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0])) +
534 *
2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
535 *
(2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]);
539 *
return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) *
540 *
p[2] * (1.0 - p[2])) +
541 *
24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) *
542 *
p[2] * (1.0 - p[2])) +
543 *
24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) *
544 *
p[1] * (1.0 - p[1])) +
545 *
2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
546 *
(2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
547 *
Utilities::fixed_power<2>(p[2] * (1.0 - p[2])) +
548 *
2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
549 *
(2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
550 *
Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) +
551 *
2.0 * (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
552 *
(2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
553 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
558 *
return return_value;
578 *
const unsigned int component = 0)
const override;
582 *
const unsigned int component = 0)
const override;
586 *
const unsigned int component = 0)
const override;
593 *
const unsigned int )
const
595 *
double return_value = 0.0;
600 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
604 *
return_value = Utilities::fixed_power<2>(
605 *
p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
610 *
return return_value;
618 *
const unsigned int )
const
625 *
(2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
626 *
4.0 * Utilities::fixed_power<3>(p[0])) *
627 *
Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
629 *
(2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
630 *
4.0 * Utilities::fixed_power<3>(p[1])) *
631 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
636 *
(2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
637 *
4.0 * Utilities::fixed_power<3>(p[0])) *
638 *
Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
640 *
(2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
641 *
4.0 * Utilities::fixed_power<3>(p[1])) *
642 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2]));
644 *
(2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
645 *
4.0 * Utilities::fixed_power<3>(p[2])) *
646 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
659 *
const unsigned int )
const
665 *
return_hessian[0][0] = (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
666 *
Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
668 *
(2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
669 *
4.0 * Utilities::fixed_power<3>(p[0])) *
670 *
(2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
671 *
4.0 * Utilities::fixed_power<3>(p[1]));
672 *
return_hessian[1][1] = (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
673 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
678 *
(2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) *
679 *
Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2]));
681 *
(2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
682 *
4.0 * Utilities::fixed_power<3>(p[0])) *
683 *
(2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
684 *
4.0 * Utilities::fixed_power<3>(p[1])) *
685 *
Utilities::fixed_power<2>(p[2] * (1.0 - p[2]));
687 *
(2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) +
688 *
4.0 * Utilities::fixed_power<3>(p[0])) *
689 *
(2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
690 *
4.0 * Utilities::fixed_power<3>(p[2])) *
691 *
Utilities::fixed_power<2>(p[1] * (1.0 - p[1]));
693 *
(2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) *
694 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2]));
696 *
(2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) +
697 *
4.0 * Utilities::fixed_power<3>(p[1])) *
698 *
(2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) +
699 *
4.0 * Utilities::fixed_power<3>(p[2])) *
700 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]));
702 *
(2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) *
703 *
Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]));
716 * <a name=
"step_82-ImplementationofthecodeBiLaplacianLDGLiftcodeclass"></a>
722 * <a name=
"step_82-BiLaplacianLDGLiftBiLaplacianLDGLift"></a>
734 *
const unsigned int fe_degree,
750 * <a name=
"step_82-BiLaplacianLDGLiftmake_grid"></a>
756 * <
code>GridGenerator::hyper_cube</code>
and then refine
it using
764 *
std::cout <<
"Building the mesh............." << std::endl;
770 *
std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
779 * <a name=
"step_82-BiLaplacianLDGLiftsetup_system"></a>
786 * right-
hand side vectors
801 *
dof_handler.distribute_dofs(fe);
803 *
std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
808 *
const auto dofs_per_cell = fe.dofs_per_cell;
810 *
for (
const auto &cell : dof_handler.active_cell_iterators())
813 *
cell->get_dof_indices(dofs);
815 *
for (
unsigned int f = 0; f < cell->n_faces(); ++f)
816 *
if (!cell->face(f)->at_boundary())
820 *
std::vector<types::global_dof_index> tmp(dofs_per_cell);
823 *
dofs.insert(std::end(dofs), std::begin(tmp), std::end(tmp));
826 *
for (
const auto i : dofs)
834 *
sparsity_pattern.copy_from(
dsp);
837 *
matrix.reinit(sparsity_pattern);
838 *
rhs.reinit(dof_handler.n_dofs());
840 *
solution.reinit(dof_handler.n_dofs());
849 *
std::ofstream out(
"sparsity-pattern.svg");
850 *
sparsity_pattern.print_svg(out);
858 * <a name=
"step_82-BiLaplacianLDGLiftassemble_system"></a>
870 *
std::cout <<
"Assembling the system............." << std::endl;
875 *
std::cout <<
"Done. " << std::endl;
883 * <a name=
"step_82-BiLaplacianLDGLiftassemble_matrix"></a>
901 *
const unsigned int n_q_points = quad.size();
902 *
const unsigned int n_q_points_face =
quad_face.size();
912 *
const unsigned int n_dofs = fe_values.dofs_per_cell;
914 *
std::vector<types::global_dof_index> local_dof_indices(n_dofs);
953 *
std::vector<std::vector<std::vector<Tensor<2, dim>>>>
957 *
for (
const auto &cell : dof_handler.active_cell_iterators())
960 *
cell->get_dof_indices(local_dof_indices);
981 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
983 *
const double dx = fe_values.JxW(
q);
985 *
for (
unsigned int i = 0; i < n_dofs; ++i)
986 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
995 *
for (
unsigned int i = 0; i < n_dofs; ++i)
996 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
998 *
matrix(local_dof_indices[i], local_dof_indices[
j]) +=
1015 *
const bool at_boundary = face->at_boundary();
1036 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1038 *
const double dx = fe_values.JxW(
q);
1040 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1042 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1062 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1064 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1066 *
matrix(local_dof_indices[i],
1070 *
local_dof_indices[
j]) +=
1093 *
const bool at_boundary = face->at_boundary();
1131 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1133 *
const double dx = fe_values.JxW(
q);
1135 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1136 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1155 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1156 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1181 *
const double mesh_inv = 1.0 / face->diameter();
1183 *
1.0 / Utilities::fixed_power<3>(face->diameter());
1189 *
const bool at_boundary = face->at_boundary();
1192 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
1196 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1197 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1214 *
cell->neighbor_of_neighbor(
face_no);
1233 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
1237 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1239 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1283 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1285 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1287 *
matrix(local_dof_indices[i], local_dof_indices[
j]) +=
1294 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1296 *
for (
unsigned int j = 0;
j < n_dofs; ++
j)
1298 *
matrix(local_dof_indices[i],
1319 * <a name=
"step_82-BiLaplacianLDGLiftassemble_rhs"></a>
1330 *
template <
int dim>
1339 *
const unsigned int n_dofs = fe_values.dofs_per_cell;
1340 *
const unsigned int n_quad_pts = quad.size();
1345 *
std::vector<types::global_dof_index> local_dof_indices(n_dofs);
1347 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1349 *
fe_values.
reinit(cell);
1350 *
cell->get_dof_indices(local_dof_indices);
1355 *
const double dx = fe_values.JxW(
q);
1357 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1361 *
fe_values.shape_value(i,
q) *
dx;
1365 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1375 * <a name=
"step_82-BiLaplacianLDGLiftsolve"></a>
1376 * <
h4>BiLaplacianLDGLift::solve</
h4>
1386 *
template <
int dim>
1399 * <a name=
"step_82-BiLaplacianLDGLiftcompute_errors"></a>
1410 *
template <
int dim>
1434 *
const unsigned int n_q_points = quad.size();
1435 *
const unsigned int n_q_points_face =
quad_face.size();
1453 *
std::vector<double> solution_values(n_q_points_face);
1455 *
std::vector<Tensor<1, dim>> solution_gradients(n_q_points_face);
1458 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1460 *
fe_values.
reinit(cell);
1471 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1473 *
const double dx = fe_values.JxW(
q);
1483 *
error_L2 += Utilities::fixed_power<2>(
1484 *
u_exact.value(fe_values.quadrature_point(
q)) -
1499 *
const double mesh_inv = 1.0 / face->diameter();
1501 *
1.0 / Utilities::fixed_power<3>(face->diameter());
1505 *
fe_face.get_function_values(solution, solution_values);
1506 *
fe_face.get_function_gradients(solution, solution_gradients);
1508 *
const bool at_boundary = face->at_boundary();
1511 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
1525 *
solution_values[
q]) *
1529 *
solution_values[
q]) *
1539 *
cell->neighbor_of_neighbor(
face_no);
1553 *
fe_face.get_function_values(solution, solution_values);
1556 *
fe_face.get_function_gradients(solution,
1557 *
solution_gradients);
1561 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
1581 *
solution_values[
q]) *
1586 *
solution_values[
q]) *
1601 *
std::cout <<
"DG H2 norm of the error: " <<
error_H2 << std::endl;
1602 *
std::cout <<
"DG H1 norm of the error: " <<
error_H1 << std::endl;
1603 *
std::cout <<
" L2 norm of the error: " <<
error_L2 << std::endl;
1611 * <a name=
"step_82-BiLaplacianLDGLiftoutput_results"></a>
1620 *
template <
int dim>
1624 *
data_out.attach_dof_handler(dof_handler);
1625 *
data_out.add_data_vector(solution,
"solution");
1626 *
data_out.build_patches();
1628 *
std::ofstream output(
"solution.vtk");
1629 *
data_out.write_vtk(output);
1637 * <a name=
"step_82-BiLaplacianLDGLiftassemble_local_matrix"></a>
1648 *
template <
int dim>
1651 *
const unsigned int n_q_points,
1659 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1663 *
for (
unsigned int m = 0; m < n_dofs; ++m)
1664 *
for (
unsigned int n = 0; n < n_dofs; ++n)
1679 * <a name=
"step_82-BiLaplacianLDGLiftcompute_discrete_hessians"></a>
1730 *
const unsigned int n_q_points = quad.size();
1731 *
const unsigned int n_q_points_face =
quad_face.size();
1751 *
const unsigned int n_dofs = fe_values.dofs_per_cell;
1783 *
fe_values.reinit(cell);
1794 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1795 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1799 *
for (
unsigned int face_no = 0;
1816 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1826 *
const bool at_boundary = face->at_boundary();
1848 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
1872 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
1902 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1937 *
const bool at_boundary = face->at_boundary();
1954 *
cell->neighbor_of_neighbor(
face_no);
1957 *
for (
unsigned int i = 0; i < n_dofs; ++i)
1965 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
1988 *
for (
unsigned int q = 0;
q < n_q_points_face; ++
q)
2011 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
2036 * <a name=
"step_82-BiLaplacianLDGLiftrun"></a>
2037 * <
h4>BiLaplacianLDGLift::run</
h4>
2040 *
template <
int dim>
2061 * <a name=
"step_82-Thecodemaincodefunction"></a>
2076 *
const unsigned int n_ref = 3;
2078 *
const unsigned int degree =
2093 *
catch (std::exception &exc)
2095 *
std::cerr << std::endl
2097 *
<<
"----------------------------------------------------"
2099 *
std::cerr <<
"Exception on processing: " << std::endl
2100 *
<< exc.what() << std::endl
2101 *
<<
"Aborting!" << std::endl
2102 *
<<
"----------------------------------------------------"
2108 *
std::cerr << std::endl
2110 *
<<
"----------------------------------------------------"
2112 *
std::cerr <<
"Unknown exception!" << std::endl
2113 *
<<
"Aborting!" << std::endl
2114 *
<<
"----------------------------------------------------"
2142<table
align=
"center" class=
"doxtable">
2231<table
align=
"center" class=
"doxtable">
2280<a name=
"step_82-PlainProg"></a>
friend class FEValuesViews::Tensor
const ObserverPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
const Quadrature< dim > quadrature
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
@ 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()
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)