const ReferenceCell::Type reference_cell_type,
const SubQuadrature & quadrature)
{
+ if (reference_cell_type == ReferenceCell::Type::Tri ||
+ reference_cell_type == ReferenceCell::Type::Tet)
+ return Quadrature<2>(); // nothing to do
+
Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
- (void)reference_cell_type;
const unsigned int dim = 2;
const ReferenceCell::Type reference_cell_type,
const SubQuadrature & quadrature)
{
+ if (reference_cell_type == ReferenceCell::Type::Tri ||
+ reference_cell_type == ReferenceCell::Type::Tet)
+ return Quadrature<3>(); // nothing to do
+
Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
- (void)reference_cell_type;
const unsigned int dim = 3;
SubQuadrature q_reflected = reflect(quadrature);
test(const FiniteElement<dim, spacedim> &fe, const unsigned int n_components)
{
Triangulation<dim, spacedim> tria;
- Simplex::GridGenerator::subdivided_hyper_cube(tria, dim == 2 ? 4 : 2);
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 4 : 2);
DoFHandler<dim> dof_handler(tria);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Solve Poisson problem on a tet mesh and on a quad mesh with the same number
+// of subdivisions.
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/distributed/fully_distributed_tria.h>
+#include <deal.II/distributed/shared_tria.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+struct Parameters
+{
+ unsigned int degree = 2;
+
+
+ // GridGenerator
+ bool use_grid_generator = true;
+ std::vector<unsigned int> repetitions;
+ Point<dim> p1;
+ Point<dim> p2;
+
+ // GridIn
+ std::string file_name_in = "";
+
+ // GridOut
+ std::string file_name_out = "";
+};
+
+template <int dim, int spacedim>
+MPI_Comm
+get_communicator(const Triangulation<dim, spacedim> &tria)
+{
+ if (auto tria_ =
+ dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
+ return tria_->get_communicator();
+
+ return MPI_COMM_SELF;
+}
+
+template <int dim, int spacedim = dim>
+void
+test(const Triangulation<dim, spacedim> &tria,
+ const FiniteElement<dim, spacedim> &fe,
+ const Quadrature<dim> & quad,
+ const Quadrature<dim - 1> & face_quad,
+ const Mapping<dim, spacedim> & mapping,
+ const double r_boundary)
+{
+ std::string label =
+ (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
+ &tria) ?
+ "parallel::shared::Triangulation" :
+ (dynamic_cast<
+ const parallel::fullydistributed::Triangulation<dim, spacedim> *>(
+ &tria) ?
+ "parallel::fullydistributed::Triangulation" :
+ (dynamic_cast<
+ const parallel::distributed::Triangulation<dim, spacedim> *>(
+ &tria) ?
+ "parallel::distributed::Triangulation" :
+ "Triangulation")));
+
+ deallog << " on " << label << std::endl;
+
+
+ for (const auto &cell : tria.active_cell_iterators())
+ for (const auto &face : cell->face_iterators())
+ if (face->at_boundary() &&
+ (std::abs(face->center()[0] - r_boundary) < 1e-6))
+ face->set_boundary_id(1);
+ else if (face->at_boundary() && face->center()[1] == 0.0)
+ face->set_boundary_id(2);
+ else if (face->at_boundary() && face->center()[1] == 1.0)
+ face->set_boundary_id(2);
+ else if (dim == 3 && face->at_boundary() && face->center()[2] == 0.0)
+ face->set_boundary_id(2);
+ else if (dim == 3 && face->at_boundary() && face->center()[2] == 1.0)
+ face->set_boundary_id(2);
+ else if (face->at_boundary())
+ face->set_boundary_id(0);
+
+
+ DoFHandler<dim, spacedim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ AffineConstraints<double> constraint_matrix;
+ DoFTools::make_zero_boundary_constraints(dof_handler, 0, constraint_matrix);
+ constraint_matrix.close();
+
+ // constraint_matrix.print(std::cout);
+
+ const MPI_Comm comm = get_communicator(dof_handler.get_triangulation());
+
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+ using VectorType = LinearAlgebra::distributed::Vector<double>;
+ TrilinosWrappers::SparseMatrix system_matrix;
+ VectorType solution;
+ VectorType system_rhs;
+#else
+ using VectorType = Vector<double>;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+ VectorType solution(dof_handler.n_dofs());
+ VectorType system_rhs(dof_handler.n_dofs());
+#endif
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+ TrilinosWrappers::SparsityPattern dsp(dof_handler.locally_owned_dofs(), comm);
+#else
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+#endif
+ DoFTools::make_sparsity_pattern(dof_handler, dsp, constraint_matrix);
+#ifdef DEAL_II_WITH_TRILINOS
+ dsp.compress();
+ system_matrix.reinit(dsp);
+
+
+ solution.reinit(dof_handler.locally_owned_dofs(),
+ locally_relevant_dofs,
+ comm);
+ system_rhs.reinit(dof_handler.locally_owned_dofs(),
+ locally_relevant_dofs,
+ comm);
+#else
+ sparsity_pattern.copy_from(dsp);
+ system_matrix.reinit(sparsity_pattern);
+#endif
+
+ const UpdateFlags flag = update_JxW_values | update_values |
+ update_gradients | update_quadrature_points;
+ FEValues<dim, spacedim> fe_values(mapping, fe, quad, flag);
+
+ std::shared_ptr<FEFaceValues<dim, spacedim>> fe_face_values;
+
+ fe_face_values.reset(
+ new FEFaceValues<dim, spacedim>(mapping, fe, face_quad, flag));
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quad.size();
+
+ std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
+
+ for (const auto &cell : dof_handler.cell_iterators())
+ {
+ if (!cell->is_locally_owned())
+ continue;
+
+ fe_values.reinit(cell);
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ cell_matrix(i, j) +=
+ (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
+ fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
+ fe_values.JxW(q_index)); // dx
+ cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
+ 1.0 * // 1.0
+ fe_values.JxW(q_index)); // dx
+ }
+
+ if (fe_face_values)
+ for (const auto &face : cell->face_iterators())
+ if (face->at_boundary() && (face->boundary_id() == 1))
+ {
+ fe_face_values->reinit(cell, face);
+ for (unsigned int q = 0; q < face_quad.size(); ++q)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ cell_rhs(i) +=
+ (1.0 * // 1.0
+ fe_face_values->shape_value(i, q) * // phi_i(x_q)
+ fe_face_values->JxW(q)); // dx
+ }
+
+ cell->get_dof_indices(dof_indices);
+
+ constraint_matrix.distribute_local_to_global(
+ cell_matrix, cell_rhs, dof_indices, system_matrix, system_rhs);
+ }
+
+ system_matrix.compress(VectorOperation::add);
+ system_rhs.compress(VectorOperation::add);
+
+ SolverControl solver_control(1000, 1e-12);
+ SolverCG<VectorType> solver(solver_control);
+ solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+ deallog << " with " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence" << std::endl;
+
+ // system_rhs.print(std::cout);
+ // solution.print(std::cout);
+
+ bool hex_mesh = true;
+
+ for (const auto &cell : tria.active_cell_iterators())
+ hex_mesh &= (cell->n_vertices() == GeometryInfo<dim>::vertices_per_cell);
+
+ deallog << std::endl;
+}
+
+template <int dim, int spacedim = dim>
+void
+test_tet(const MPI_Comm &comm, const Parameters<dim> ¶ms)
+{
+ const unsigned int tria_type = 2;
+
+ // 1) Create triangulation...
+ Triangulation<dim, spacedim> *tria;
+
+ // a) serial triangulation
+ Triangulation<dim, spacedim> tr_1;
+
+ // b) shared triangulation (with artificial cells)
+ parallel::shared::Triangulation<dim> tr_2(
+ MPI_COMM_WORLD,
+ ::Triangulation<dim>::none,
+ true,
+ parallel::shared::Triangulation<dim>::partition_custom_signal);
+
+ tr_2.signals.create.connect([&]() {
+ GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm),
+ tr_2);
+ });
+
+ // c) distributed triangulation
+ parallel::fullydistributed::Triangulation<dim> tr_3(comm);
+
+
+ // ... choose the right triangulation
+ if (tria_type == 0 || tria_type == 2)
+ tria = &tr_1;
+ else if (tria_type == 1)
+ tria = &tr_2;
+
+ // ... create triangulation
+ if (params.use_grid_generator)
+ {
+ // ...via Simplex::GridGenerator
+ GridGenerator::subdivided_hyper_rectangle_with_simplices(
+ *tria, params.repetitions, params.p1, params.p2, false);
+ }
+ else
+ {
+ // ...via GridIn
+ GridIn<dim, spacedim> grid_in;
+ grid_in.attach_triangulation(*tria);
+ std::ifstream input_file(params.file_name_in);
+ grid_in.read_ucd(input_file);
+ // std::ifstream input_file("test_tet_geometry.unv");
+ // grid_in.read_unv(input_file);
+ }
+
+ // ... partition serial triangulation and create distributed triangulation
+ if (tria_type == 0 || tria_type == 2)
+ {
+ GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm),
+ tr_1);
+
+ auto construction_data = TriangulationDescription::Utilities::
+ create_description_from_triangulation(tr_1, comm);
+
+ tr_3.create_triangulation(construction_data);
+
+ tria = &tr_3;
+ }
+
+ // 2) Output generated triangulation via GridOut
+ GridOut grid_out;
+ std::ofstream out(params.file_name_out + "." +
+ std::to_string(Utilities::MPI::this_mpi_process(comm)) +
+ ".vtk");
+ grid_out.write_vtk(*tria, out);
+
+ // 3) Select components
+ Simplex::FE_P<dim> fe(params.degree);
+
+ Simplex::PGauss<dim> quad(dim == 2 ? (params.degree == 1 ? 3 : 7) :
+ (params.degree == 1 ? 4 : 10));
+
+ Simplex::PGauss<dim - 1> face_quad(dim == 2 ? (params.degree == 1 ? 2 : 3) :
+ (params.degree == 1 ? 3 : 7));
+
+ Simplex::FE_P<dim> fe_mapping(1);
+ MappingFE<dim> mapping(fe_mapping);
+
+ // 4) Perform test (independent of mesh type)
+ test(*tria, fe, quad, face_quad, mapping, params.p2[0]);
+}
+
+template <int dim, int spacedim = dim>
+void
+test_hex(const MPI_Comm &comm, const Parameters<dim> ¶ms)
+{
+ // 1) Create triangulation...
+ parallel::distributed::Triangulation<dim, spacedim> tria(comm);
+
+ if (params.use_grid_generator)
+ {
+ // ...via GridGenerator
+ GridGenerator::subdivided_hyper_rectangle(
+ tria, params.repetitions, params.p1, params.p2, false);
+ }
+ else
+ {
+ // ...via GridIn
+ GridIn<dim, spacedim> grid_in;
+ grid_in.attach_triangulation(tria);
+ std::ifstream input_file(params.file_name_in);
+ grid_in.read_ucd(input_file);
+ }
+
+ // 2) Output generated triangulation via GridOut
+ GridOut grid_out;
+ std::ofstream out(params.file_name_out + "." +
+ std::to_string(Utilities::MPI::this_mpi_process(comm)) +
+ ".vtk");
+ grid_out.write_vtk(tria, out);
+
+ // 3) Select components
+ FE_Q<dim> fe(params.degree);
+
+ QGauss<dim> quad(params.degree + 1);
+
+ QGauss<dim - 1> quad_face(params.degree + 1);
+
+ MappingQ<dim, spacedim> mapping(1);
+
+ // 4) Perform test (independent of mesh type)
+ test(tria, fe, quad, quad_face, mapping, params.p2[0]);
+}
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ initlog();
+
+ const MPI_Comm comm = MPI_COMM_WORLD;
+
+ // 2D
+ {
+ Parameters<2> params;
+ params.use_grid_generator = true;
+ params.repetitions = std::vector<unsigned int>{10, 10};
+
+ // test TRI
+ {
+ deallog << "Solve problem on TRI mesh:" << std::endl;
+
+ params.file_name_out = "mesh-tri";
+ params.p1 = Point<2>(0, 0);
+ params.p2 = Point<2>(1, 1);
+ test_tet(comm, params);
+ }
+
+ // test QUAD
+ {
+ deallog << "Solve problem on QUAD mesh:" << std::endl;
+
+ params.file_name_out = "mesh-quad";
+ params.p1 = Point<2>(1.1, 0); // shift to the right for
+ params.p2 = Point<2>(2.1, 1); // visualization purposes
+ test_hex(comm, params);
+ }
+ }
+
+ // 3D
+ {
+ Parameters<3> params;
+ params.use_grid_generator = true;
+ params.repetitions = std::vector<unsigned int>{10, 10, 10};
+
+ // test TET
+ {
+ deallog << "Solve problem on TET mesh:" << std::endl;
+
+ params.file_name_out = "mesh-tet";
+ params.p1 = Point<3>(0, 0, 0);
+ params.p2 = Point<3>(1, 1, 1);
+ test_tet(comm, params);
+ }
+
+ // test HEX
+ {
+ deallog << "Solve problem on HEX mesh:" << std::endl;
+
+ params.file_name_out = "mesh-hex";
+ params.p1 = Point<3>(1.1, 0, 0);
+ params.p2 = Point<3>(2.1, 1, 1);
+ test_hex(comm, params);
+ }
+ }
+}
--- /dev/null
+
+DEAL:0::Solve problem on TRI mesh:
+DEAL:0:: on parallel::fullydistributed::Triangulation
+DEAL:0:cg::Starting value 0.245798
+DEAL:0:cg::Convergence step 114 value 6.62438e-13
+DEAL:0:: with 114 CG iterations needed to obtain convergence
+DEAL:0::
+DEAL:0::Solve problem on QUAD mesh:
+DEAL:0:: on parallel::distributed::Triangulation
+DEAL:0:cg::Starting value 0.244628
+DEAL:0:cg::Convergence step 98 value 8.55703e-13
+DEAL:0:: with 98 CG iterations needed to obtain convergence
+DEAL:0::
+DEAL:0::Solve problem on TET mesh:
+DEAL:0:: on parallel::fullydistributed::Triangulation
+DEAL:0:cg::Starting value 0.0607616
+DEAL:0:cg::Convergence step 156 value 9.07624e-13
+DEAL:0:: with 156 CG iterations needed to obtain convergence
+DEAL:0::
+DEAL:0::Solve problem on HEX mesh:
+DEAL:0:: on parallel::distributed::Triangulation
+DEAL:0:cg::Starting value 0.0573704
+DEAL:0:cg::Convergence step 134 value 8.23917e-13
+DEAL:0:: with 134 CG iterations needed to obtain convergence
+DEAL:0::
--- /dev/null
+
+DEAL::Solve problem on TRI mesh:
+DEAL:: on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.245798
+DEAL:cg::Convergence step 114 value 6.95996e-13
+DEAL:: with 114 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on QUAD mesh:
+DEAL:: on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.244628
+DEAL:cg::Convergence step 98 value 8.55703e-13
+DEAL:: with 98 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on TET mesh:
+DEAL:: on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.0607616
+DEAL:cg::Convergence step 156 value 9.07635e-13
+DEAL:: with 156 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on HEX mesh:
+DEAL:: on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.0573704
+DEAL:cg::Convergence step 134 value 8.23917e-13
+DEAL:: with 134 CG iterations needed to obtain convergence
+DEAL::
--- /dev/null
+
+DEAL::Solve problem on TRI mesh:
+DEAL:: on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.245798
+DEAL:cg::Convergence step 114 value 6.45793e-13
+DEAL:: with 114 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on QUAD mesh:
+DEAL:: on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.244628
+DEAL:cg::Convergence step 98 value 8.55703e-13
+DEAL:: with 98 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on TET mesh:
+DEAL:: on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.0607616
+DEAL:cg::Convergence step 156 value 9.07640e-13
+DEAL:: with 156 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on HEX mesh:
+DEAL:: on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.0573704
+DEAL:cg::Convergence step 134 value 8.23917e-13
+DEAL:: with 134 CG iterations needed to obtain convergence
+DEAL::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Solve Poisson problem on a tet mesh with DG.
+
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+//#define HEX
+
+using namespace dealii;
+
+template <int dim>
+struct ScratchData
+{
+ ScratchData(const Mapping<dim> & mapping,
+ const FiniteElement<dim> & fe,
+ const Quadrature<dim> & quad,
+ const Quadrature<dim - 1> &quad_face,
+ const UpdateFlags update_flags = update_values |
+ update_gradients |
+ update_quadrature_points |
+ update_JxW_values,
+ const UpdateFlags interface_update_flags =
+ update_values | update_gradients | update_quadrature_points |
+ update_JxW_values | update_normal_vectors)
+ : fe_values(mapping, fe, quad, update_flags)
+ , fe_interface_values(mapping, fe, quad_face, interface_update_flags)
+ {}
+
+
+ ScratchData(const ScratchData<dim> &scratch_data)
+ : fe_values(scratch_data.fe_values.get_mapping(),
+ scratch_data.fe_values.get_fe(),
+ scratch_data.fe_values.get_quadrature(),
+ scratch_data.fe_values.get_update_flags())
+ , fe_interface_values(scratch_data.fe_values.get_mapping(),
+ scratch_data.fe_values.get_fe(),
+ scratch_data.fe_interface_values.get_quadrature(),
+ scratch_data.fe_interface_values.get_update_flags())
+ {}
+
+ FEValues<dim> fe_values;
+ FEInterfaceValues<dim> fe_interface_values;
+};
+
+
+
+struct CopyDataFace
+{
+ FullMatrix<double> cell_matrix;
+ std::vector<types::global_dof_index> joint_dof_indices;
+};
+
+
+
+struct CopyData
+{
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+ std::vector<types::global_dof_index> local_dof_indices;
+ std::vector<CopyDataFace> face_data;
+
+ template <class Iterator>
+ void
+ reinit(const Iterator &cell, unsigned int dofs_per_cell)
+ {
+ cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+ cell_rhs.reinit(dofs_per_cell);
+
+ local_dof_indices.resize(dofs_per_cell);
+ cell->get_dof_indices(local_dof_indices);
+ }
+};
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+ RightHandSideFunction()
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
+ {
+ if (dim == 2)
+ return -2. * M_PI * M_PI * std::sin(M_PI * p(0)) * std::sin(M_PI * p(1));
+ else /* if(dim == 3)*/
+ return -3. * M_PI * M_PI * std::sin(M_PI * p(0)) * std::sin(M_PI * p(1)) *
+ std::sin(M_PI * p(2));
+ }
+};
+
+template <int dim>
+class DGHeat
+{
+public:
+ DGHeat(const bool hex,
+ FiniteElement<dim> * fe,
+ Mapping<dim> * mapping,
+ Quadrature<dim> * quad,
+ Quadrature<dim - 1> *face_quad,
+ unsigned int initial_refinement,
+ unsigned int number_refinement)
+ : hex(hex)
+ , fe(fe)
+ , mapping(mapping)
+ , quad(quad)
+ , face_quad(face_quad)
+ , dof_handler(triangulation)
+ , initial_refinement_level(initial_refinement)
+ , number_refinement(number_refinement)
+ {}
+
+ static std::unique_ptr<DGHeat<dim>>
+ HEX(unsigned int degree,
+ unsigned int initial_refinement,
+ unsigned int number_refinement)
+ {
+ return std::make_unique<DGHeat<dim>>(true,
+ new FE_DGQ<dim>(degree),
+ new MappingQ<dim>(1),
+ new QGauss<dim>(degree + 1),
+ new QGauss<dim - 1>(degree + 1),
+ initial_refinement,
+ number_refinement);
+ }
+
+ static std::unique_ptr<DGHeat<dim>>
+ TET(unsigned int degree,
+ unsigned int initial_refinement,
+ unsigned int number_refinement)
+ {
+ return std::make_unique<DGHeat<dim>>(
+ false,
+ new Simplex::FE_DGP<dim>(degree),
+ new MappingFE<dim>(Simplex::FE_P<dim>(1)),
+ new Simplex::PGauss<dim>(dim == 2 ? (degree == 1 ? 3 : 7) :
+ (degree == 1 ? 4 : 10)),
+ new Simplex::PGauss<dim - 1>(dim == 2 ? (degree == 1 ? 2 : 3) :
+ (degree == 1 ? 3 : 7)),
+ initial_refinement,
+ number_refinement);
+ }
+
+
+
+ void
+ run();
+
+
+
+private:
+ void
+ make_grid(int refinements = -1);
+ void
+ setup_system();
+ void
+ assemble_system();
+ void
+ solve();
+ void
+ output_results(unsigned int it) const;
+ void
+ calculateL2Error();
+
+ Triangulation<dim> triangulation;
+
+ bool hex;
+
+ std::unique_ptr<FiniteElement<dim>> fe;
+ const std::unique_ptr<Mapping<dim>> mapping;
+ std::unique_ptr<Quadrature<dim>> quad;
+ std::unique_ptr<Quadrature<dim - 1>> face_quad;
+
+ DoFHandler<dim> dof_handler;
+
+
+ RightHandSideFunction<dim> right_hand_side;
+
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+ Point<dim> center;
+
+ ConvergenceTable error_table;
+
+ unsigned int initial_refinement_level;
+ unsigned int number_refinement;
+};
+
+
+template <int dim>
+void
+DGHeat<dim>::make_grid(int refinements)
+{
+ triangulation.clear();
+
+ const unsigned int ref =
+ refinements == -1 ? initial_refinement_level : refinements;
+
+ if (hex)
+ GridGenerator::subdivided_hyper_cube(triangulation,
+ Utilities::pow(2, ref),
+ -1.0,
+ +1.0);
+ else
+ GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+ Utilities::pow(2, ref),
+ -1.0,
+ +1.0);
+
+ // deallog << " Number of active cells: " <<
+ // triangulation.n_active_cells()
+ // << std::endl
+ // << " Total number of cells: " << triangulation.n_cells()
+ // << std::endl;
+}
+
+
+template <int dim>
+void
+DGHeat<dim>::setup_system()
+{
+ dof_handler.distribute_dofs(*fe);
+
+ // deallog << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ // << std::endl;
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+ sparsity_pattern.copy_from(dsp);
+
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(dof_handler.n_dofs());
+ system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void
+DGHeat<dim>::assemble_system()
+{
+ using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+
+ auto cell_worker = [&](const Iterator & cell,
+ ScratchData<dim> &scratch_data,
+ CopyData & copy_data) {
+ const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell;
+ copy_data.reinit(cell, n_dofs);
+ scratch_data.fe_values.reinit(cell);
+
+ const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+ const FEValues<dim> & fe_v = scratch_data.fe_values;
+ const std::vector<double> &JxW = fe_v.get_JxW_values();
+
+ std::vector<double> f(q_points.size());
+ right_hand_side.value_list(q_points, f);
+
+ for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+ {
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ {
+ for (unsigned int j = 0; j < n_dofs; ++j)
+ {
+ copy_data.cell_matrix(i, j) +=
+ fe_v.shape_grad(i, point) // \nabla \phi_i
+ * fe_v.shape_grad(j, point) // \nabla \phi_j
+ * JxW[point]; // dx
+ }
+
+ // Right Hand Side
+ copy_data.cell_rhs(i) +=
+ (fe_v.shape_value(i, point) * f[point] * JxW[point]);
+ }
+ }
+ };
+
+ auto boundary_worker = [&](const Iterator & cell,
+ const unsigned int &face_no,
+ ScratchData<dim> & scratch_data,
+ CopyData & copy_data) {
+ scratch_data.fe_interface_values.reinit(cell, face_no);
+
+ const FEFaceValuesBase<dim> &fe_face =
+ scratch_data.fe_interface_values.get_fe_face_values(0);
+
+ const auto & q_points = fe_face.get_quadrature_points();
+ const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
+ const std::vector<double> &JxW = fe_face.get_JxW_values();
+
+ const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
+
+ double h;
+ if (dim == 2)
+ {
+ if (hex)
+ h = std::sqrt(4. * cell->measure() / M_PI);
+ else
+ h = std::sqrt(4. * (4.0 / triangulation.n_cells()) / M_PI);
+ }
+ else if (dim == 3)
+ {
+ if (hex)
+ h = pow(6 * cell->measure() / M_PI, 1. / 3.);
+ else
+ h = pow(6 * (8.0 / triangulation.n_cells()) / M_PI, 1. / 3.);
+ }
+
+
+
+ const double beta = 10.;
+
+ for (unsigned int point = 0; point < q_points.size(); ++point)
+ for (unsigned int i = 0; i < n_facet_dofs; ++i)
+ for (unsigned int j = 0; j < n_facet_dofs; ++j)
+ {
+ copy_data.cell_matrix(i, j) +=
+ -normals[point] * fe_face.shape_grad(i, point) // n*\nabla \phi_i
+ * fe_face.shape_value(j, point) // \phi_j
+ * JxW[point]; // dx
+
+ copy_data.cell_matrix(i, j) +=
+ -fe_face.shape_value(i, point) // \phi_i
+ * fe_face.shape_grad(j, point) * normals[point] // n*\nabla \phi_j
+ * JxW[point]; // dx
+
+ copy_data.cell_matrix(i, j) +=
+ beta * 1. / h * fe_face.shape_value(i, point) // \phi_i
+ * fe_face.shape_value(j, point) * JxW[point]; // dx
+ }
+ };
+
+ auto face_worker = [&](const Iterator & cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const Iterator & ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ ScratchData<dim> & scratch_data,
+ CopyData & copy_data) {
+ FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
+
+ fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
+
+ const auto &q_points = fe_iv.get_quadrature_points();
+
+ copy_data.face_data.emplace_back();
+ CopyDataFace ©_data_face = copy_data.face_data.back();
+
+ const unsigned int n_dofs = fe_iv.n_current_interface_dofs();
+ copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
+
+ copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+ const std::vector<double> & JxW = fe_iv.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+
+ double h;
+ if (dim == 2)
+ {
+ if (hex)
+ h = std::sqrt(4. * cell->measure() / M_PI);
+ else
+ h = std::sqrt(4. * (4.0 / triangulation.n_cells()) / M_PI);
+ }
+ else if (dim == 3)
+ {
+ if (hex)
+ h = pow(6 * cell->measure() / M_PI, 1. / 3.);
+ else
+ h = pow(6 * (8.0 / triangulation.n_cells()) / M_PI, 1. / 3.);
+ }
+
+ const double beta = 10.;
+
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ {
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ {
+ for (unsigned int j = 0; j < n_dofs; ++j)
+ {
+ copy_data_face.cell_matrix(i, j) +=
+ -normals[qpoint] * fe_iv.average_gradient(i, qpoint) *
+ fe_iv.jump(j, qpoint) * JxW[qpoint];
+
+ copy_data_face.cell_matrix(i, j) +=
+ -fe_iv.jump(i, qpoint) // \phi_i
+ * fe_iv.average_gradient(j, qpoint) *
+ normals[qpoint] // n*\nabla \phi_j
+ * JxW[qpoint]; // dx
+
+ copy_data_face.cell_matrix(i, j) +=
+ beta * 1. / h * fe_iv.jump(i, qpoint) *
+ fe_iv.jump(j, qpoint) * JxW[qpoint];
+ }
+ }
+ }
+ };
+
+ AffineConstraints<double> constraints;
+
+ auto copier = [&](const CopyData &c) {
+ constraints.distribute_local_to_global(c.cell_matrix,
+ c.cell_rhs,
+ c.local_dof_indices,
+ system_matrix,
+ system_rhs);
+
+ for (auto &cdf : c.face_data)
+ {
+ constraints.distribute_local_to_global(cdf.cell_matrix,
+ cdf.joint_dof_indices,
+ system_matrix);
+ }
+ };
+
+
+ ScratchData<dim> scratch_data(*mapping, *fe, *quad, *face_quad);
+ CopyData copy_data;
+
+ MeshWorker::mesh_loop(dof_handler.begin_active(),
+ dof_handler.end(),
+ cell_worker,
+ copier,
+ scratch_data,
+ copy_data,
+ MeshWorker::assemble_own_cells |
+ MeshWorker::assemble_boundary_faces |
+ MeshWorker::assemble_own_interior_faces_once,
+ boundary_worker,
+ face_worker);
+}
+
+template <int dim>
+void
+DGHeat<dim>::solve()
+{
+ SolverControl solver_control(10000, 1e-8);
+ SolverCG<> solver(solver_control);
+ solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+ // We have made one addition, though: since we suppress output from the
+ // linear solvers, we have to print the number of iterations by hand.
+ // deallog << " " << solver_control.last_step()
+ // << " CG iterations needed to obtain convergence." << std::endl;
+
+ // error_table.add_value("iterations", solver_control.last_step());
+}
+
+template <int dim>
+void
+DGHeat<dim>::output_results(unsigned int it) const
+{
+ return;
+
+ std::string type = hex ? "hex" : "tet";
+
+ std::string dimension(dim == 2 ? "solution-2d-" + type + "-case-" :
+ "solution-3d-" + type + "-case-");
+
+ std::string fname = dimension + Utilities::int_to_string(it) + ".vtk";
+
+ deallog << " Writing solution to <" << fname << ">" << std::endl;
+
+ std::ofstream output(fname.c_str());
+
+ if (false)
+ {
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "solution");
+
+ data_out.build_patches(*mapping);
+ data_out.write_vtk(output);
+ }
+}
+
+// Find the l2 norm of the error between the finite element sol'n and the exact
+// sol'n
+template <int dim>
+void
+DGHeat<dim>::calculateL2Error()
+{
+ FEValues<dim> fe_values(*mapping,
+ *fe,
+ *quad,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell =
+ fe->dofs_per_cell; // This gives you dofs per cell
+ std::vector<types::global_dof_index> local_dof_indices(
+ dofs_per_cell); // Local connectivity
+
+ const unsigned int n_q_points = quad->size();
+
+ double l2error = 0.;
+
+ // loop over elements
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ fe_values.reinit(cell);
+
+ cell->get_dof_indices(local_dof_indices);
+
+ for (unsigned int q = 0; q < n_q_points; q++)
+ {
+ const double u_exact =
+ dim == 2 ? -std::sin(M_PI * fe_values.quadrature_point(q)[0]) *
+ std::sin(M_PI * fe_values.quadrature_point(q)[1]) :
+ -std::sin(M_PI * fe_values.quadrature_point(q)[0]) *
+ std::sin(M_PI * fe_values.quadrature_point(q)[1]) *
+ std::sin(M_PI * fe_values.quadrature_point(q)[2]);
+
+ double u_sim = 0;
+
+ // Find the values of x and u_h (the finite element solution) at the
+ // quadrature points
+ for (unsigned int i = 0; i < dofs_per_cell; i++)
+ {
+ u_sim +=
+ fe_values.shape_value(i, q) * solution[local_dof_indices[i]];
+ }
+ l2error += (u_sim - u_exact) * (u_sim - u_exact) * fe_values.JxW(q);
+ // deallog << " x = " << x << " y = " << y << " r = " << r <<
+ // " u_exact = " << u_exact << " u_sim=" << u_sim <<
+ // std::endl;
+ }
+ }
+
+
+ // deallog << "L2Error is : " << std::sqrt(l2error) << std::endl;
+ error_table.add_value("error", std::sqrt(l2error));
+ error_table.add_value("cells", triangulation.n_global_active_cells());
+ error_table.add_value("dofs", dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void
+DGHeat<dim>::run()
+{
+ for (unsigned int it = 0; it < number_refinement; ++it)
+ {
+ make_grid(initial_refinement_level + it);
+ setup_system();
+ assemble_system();
+ solve();
+ output_results(it);
+ calculateL2Error();
+ }
+
+ // error_table.omit_column_from_convergence_rate_evaluation("iterations");
+ error_table.omit_column_from_convergence_rate_evaluation("cells");
+ error_table.evaluate_all_convergence_rates(
+ ConvergenceTable::reduction_rate_log2);
+
+ error_table.set_scientific("error", true);
+
+ error_table.write_text(deallog.get_file_stream());
+ deallog << std::endl;
+}
+
+int
+main()
+{
+ initlog();
+
+ deallog.depth_file(1);
+
+
+ {
+ auto problem = DGHeat<2>::TET(1 /*=degree*/, 2, 3);
+ problem->run();
+ }
+ {
+ auto problem = DGHeat<2>::TET(2 /*=degree*/, 2, 3);
+ problem->run();
+ }
+ {
+ auto problem = DGHeat<3>::TET(1 /*=degree*/, 2, 2);
+ problem->run();
+ }
+ {
+ auto problem = DGHeat<3>::TET(2 /*=degree*/, 2, 2);
+ problem->run();
+ }
+
+ {
+ auto problem = DGHeat<2>::HEX(1 /*=degree*/, 2, 3);
+ problem->run();
+ }
+ {
+ auto problem = DGHeat<2>::HEX(2 /*=degree*/, 2, 3);
+ problem->run();
+ }
+ {
+ auto problem = DGHeat<3>::HEX(1 /*=degree*/, 2, 2);
+ problem->run();
+ }
+ {
+ auto problem = DGHeat<3>::HEX(2 /*=degree*/, 2, 2);
+ problem->run();
+ }
+
+ return 0;
+}
--- /dev/null
+
+ error cells dofs
+3.0924e-01 - 32 96 -
+1.1363e-01 1.44 128 384 -2.00
+3.2531e-02 1.80 512 1536 -2.00
+DEAL::
+ error cells dofs
+3.2745e-02 - 32 192 -
+3.6158e-03 3.18 128 768 -2.00
+4.1566e-04 3.12 512 3072 -2.00
+DEAL::
+ error cells dofs
+1.6346e-01 - 320 1280 -
+1.0816e-01 0.60 2560 10240 -3.00
+DEAL::
+ error cells dofs
+8.4685e-02 - 320 3200 -
+9.4765e-03 3.16 2560 25600 -3.00
+DEAL::
+ error cells dofs
+1.6226e-01 - 16 64 -
+4.8634e-02 1.74 64 256 -2.00
+1.2684e-02 1.94 256 1024 -2.00
+DEAL::
+ error cells dofs
+1.6600e-02 - 16 144 -
+1.6552e-03 3.33 64 576 -2.00
+1.7799e-04 3.22 256 2304 -2.00
+DEAL::
+ error cells dofs
+1.6067e-01 - 64 512 -
+4.8569e-02 1.73 512 4096 -3.00
+DEAL::
+ error cells dofs
+2.2671e-02 - 64 1728 -
+2.2931e-03 3.31 512 13824 -3.00
+DEAL::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Step-12 with tetrahedron mesh.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/derivative_approximation.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+// Finally, the new include file for using the mesh_loop from the MeshWorker
+// framework
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include <fstream>
+#include <iostream>
+
+//#define HEX
+
+namespace Step12
+{
+ using namespace dealii;
+
+ template <int dim>
+ class BoundaryValues : public Function<dim>
+ {
+ public:
+ BoundaryValues() = default;
+ virtual void
+ value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & values,
+ const unsigned int component = 0) const override;
+ };
+
+ template <int dim>
+ void
+ BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & values,
+ const unsigned int component) const
+ {
+ (void)component;
+ AssertIndexRange(component, 1);
+ Assert(values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i = 0; i < values.size(); ++i)
+ {
+ if (points[i](0) < 0.5)
+ values[i] = 1.;
+ else
+ values[i] = 0.;
+ }
+ }
+
+
+ template <int dim>
+ Tensor<1, dim>
+ beta(const Point<dim> &p)
+ {
+ Assert(dim >= 2, ExcNotImplemented());
+
+ Point<dim> wind_field;
+ wind_field(0) = -p(1);
+ wind_field(1) = p(0);
+
+ if (wind_field.norm() > 1e-6)
+ wind_field /= wind_field.norm();
+
+ return wind_field;
+ }
+
+
+ template <int dim>
+ struct ScratchData
+ {
+ ScratchData(const Mapping<dim> & mapping,
+ const FiniteElement<dim> & fe,
+ const Quadrature<dim> & quad,
+ const Quadrature<dim - 1> &quad_face,
+ const UpdateFlags update_flags = update_values |
+ update_gradients |
+ update_quadrature_points |
+ update_JxW_values,
+ const UpdateFlags interface_update_flags =
+ update_values | update_gradients | update_quadrature_points |
+ update_JxW_values | update_normal_vectors)
+ : fe_values(mapping, fe, quad, update_flags)
+ , fe_interface_values(mapping, fe, quad_face, interface_update_flags)
+ {}
+
+
+ ScratchData(const ScratchData<dim> &scratch_data)
+ : fe_values(scratch_data.fe_values.get_mapping(),
+ scratch_data.fe_values.get_fe(),
+ scratch_data.fe_values.get_quadrature(),
+ scratch_data.fe_values.get_update_flags())
+ , fe_interface_values(
+ scratch_data.fe_values
+ .get_mapping(), // TODO: implement for fe_interface_values
+ scratch_data.fe_values.get_fe(),
+ scratch_data.fe_interface_values.get_quadrature(),
+ scratch_data.fe_interface_values.get_update_flags())
+ {}
+
+ FEValues<dim> fe_values;
+ FEInterfaceValues<dim> fe_interface_values;
+ };
+
+
+
+ struct CopyDataFace
+ {
+ FullMatrix<double> cell_matrix;
+ std::vector<types::global_dof_index> joint_dof_indices;
+ };
+
+
+
+ struct CopyData
+ {
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+ std::vector<types::global_dof_index> local_dof_indices;
+ std::vector<CopyDataFace> face_data;
+
+ template <class Iterator>
+ void
+ reinit(const Iterator &cell, unsigned int dofs_per_cell)
+ {
+ cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+ cell_rhs.reinit(dofs_per_cell);
+
+ local_dof_indices.resize(dofs_per_cell);
+ cell->get_dof_indices(local_dof_indices);
+ }
+ };
+
+
+ template <int dim>
+ class AdvectionProblem
+ {
+ public:
+ AdvectionProblem();
+ void
+ run();
+
+ private:
+ void
+ setup_system();
+ void
+ assemble_system();
+ void
+ solve();
+ void
+ refine_grid();
+ void
+ output_results(const unsigned int cycle) const;
+
+ Triangulation<dim> triangulation;
+#ifdef HEX
+ const MappingQ1<dim> mapping;
+#else
+ Simplex::FE_P<dim> fe_mapping;
+ const MappingFE<dim> mapping;
+#endif
+
+ // Furthermore we want to use DG elements.
+#ifdef HEX
+ FE_DGQ<dim> fe;
+#else
+ Simplex::FE_DGP<dim> fe;
+#endif
+ DoFHandler<dim> dof_handler;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> right_hand_side;
+ };
+
+
+ template <int dim>
+ AdvectionProblem<dim>::AdvectionProblem()
+#ifdef HEX
+ : mapping()
+ , fe(2)
+#else
+ : fe_mapping(1)
+ , mapping(fe_mapping)
+ , fe(2)
+#endif
+ , dof_handler(triangulation)
+ {}
+
+
+ template <int dim>
+ void
+ AdvectionProblem<dim>::setup_system()
+ {
+ dof_handler.distribute_dofs(fe);
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+ sparsity_pattern.copy_from(dsp);
+
+ system_matrix.reinit(sparsity_pattern);
+ solution.reinit(dof_handler.n_dofs());
+ right_hand_side.reinit(dof_handler.n_dofs());
+ }
+
+ template <int dim>
+ void
+ AdvectionProblem<dim>::assemble_system()
+ {
+ using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+ const BoundaryValues<dim> boundary_function;
+
+ auto cell_worker = [&](const Iterator & cell,
+ ScratchData<dim> &scratch_data,
+ CopyData & copy_data) {
+ const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell;
+ copy_data.reinit(cell, n_dofs);
+ scratch_data.fe_values.reinit(cell);
+
+ const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+ const FEValues<dim> & fe_v = scratch_data.fe_values;
+ const std::vector<double> &JxW = fe_v.get_JxW_values();
+
+ for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+ {
+ auto beta_q = beta(q_points[point]);
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ for (unsigned int j = 0; j < n_dofs; ++j)
+ {
+ copy_data.cell_matrix(i, j) +=
+ -beta_q // -\beta
+ * fe_v.shape_grad(i, point) // \nabla \phi_i
+ * fe_v.shape_value(j, point) // \phi_j
+ * JxW[point]; // dx
+ }
+ }
+ };
+
+ auto boundary_worker = [&](const Iterator & cell,
+ const unsigned int &face_no,
+ ScratchData<dim> & scratch_data,
+ CopyData & copy_data) {
+ scratch_data.fe_interface_values.reinit(cell, face_no);
+ const FEFaceValuesBase<dim> &fe_face =
+ scratch_data.fe_interface_values.get_fe_face_values(0);
+
+ const auto &q_points = fe_face.get_quadrature_points();
+
+ const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
+ const std::vector<double> & JxW = fe_face.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
+
+ std::vector<double> g(q_points.size());
+ boundary_function.value_list(q_points, g);
+
+ for (unsigned int point = 0; point < q_points.size(); ++point)
+ {
+ const double beta_dot_n = beta(q_points[point]) * normals[point];
+
+ if (beta_dot_n > 0)
+ {
+ for (unsigned int i = 0; i < n_facet_dofs; ++i)
+ for (unsigned int j = 0; j < n_facet_dofs; ++j)
+ copy_data.cell_matrix(i, j) +=
+ fe_face.shape_value(i, point) // \phi_i
+ * fe_face.shape_value(j, point) // \phi_j
+ * beta_dot_n // \beta . n
+ * JxW[point]; // dx
+ }
+ else
+ for (unsigned int i = 0; i < n_facet_dofs; ++i)
+ copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
+ * g[point] // g
+ * beta_dot_n // \beta . n
+ * JxW[point]; // dx
+ }
+ };
+
+ auto face_worker = [&](const Iterator & cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const Iterator & ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ ScratchData<dim> & scratch_data,
+ CopyData & copy_data) {
+ FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
+ fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
+ const auto &q_points = fe_iv.get_quadrature_points();
+
+ copy_data.face_data.emplace_back();
+ CopyDataFace ©_data_face = copy_data.face_data.back();
+
+ const unsigned int n_dofs = fe_iv.n_current_interface_dofs();
+ copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
+
+ copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+ const std::vector<double> & JxW = fe_iv.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ {
+ const double beta_dot_n = beta(q_points[qpoint]) * normals[qpoint];
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ for (unsigned int j = 0; j < n_dofs; ++j)
+ copy_data_face.cell_matrix(i, j) +=
+ fe_iv.jump(i, qpoint) // [\phi_i]
+ *
+ fe_iv.shape_value((beta_dot_n > 0), j, qpoint) // phi_j^{upwind}
+ * beta_dot_n // (\beta . n)
+ * JxW[qpoint]; // dx
+ }
+ };
+
+ AffineConstraints<double> constraints;
+
+ auto copier = [&](const CopyData &c) {
+ constraints.distribute_local_to_global(c.cell_matrix,
+ c.cell_rhs,
+ c.local_dof_indices,
+ system_matrix,
+ right_hand_side);
+
+ for (auto &cdf : c.face_data)
+ {
+ constraints.distribute_local_to_global(cdf.cell_matrix,
+ cdf.joint_dof_indices,
+ system_matrix);
+ }
+ };
+
+ const unsigned int degree = dof_handler.get_fe().degree;
+
+#ifdef HEX
+ QGauss<dim> quad(degree + 1);
+
+ QGauss<dim - 1> face_quad(degree + 1);
+#else
+ Simplex::PGauss<dim> quad(dim == 2 ? (degree == 1 ? 3 : 7) :
+ (degree == 1 ? 4 : 10));
+
+ Simplex::PGauss<dim - 1> face_quad(dim == 2 ? (degree == 1 ? 2 : 3) :
+ (degree == 1 ? 3 : 7));
+#endif
+
+ ScratchData<dim> scratch_data(mapping, fe, quad, face_quad);
+ CopyData copy_data;
+
+ MeshWorker::mesh_loop(dof_handler.begin_active(),
+ dof_handler.end(),
+ cell_worker,
+ copier,
+ scratch_data,
+ copy_data,
+ MeshWorker::assemble_own_cells |
+ MeshWorker::assemble_boundary_faces |
+ MeshWorker::assemble_own_interior_faces_once,
+ boundary_worker,
+ face_worker);
+ }
+
+ template <int dim>
+ void
+ AdvectionProblem<dim>::solve()
+ {
+ SolverControl solver_control(1000, 1e-12);
+ SolverRichardson<Vector<double>> solver(solver_control);
+
+ PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
+
+ preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+
+ solver.solve(system_matrix, solution, right_hand_side, preconditioner);
+
+ deallog << " Solver converged in " << solver_control.last_step()
+ << " iterations." << std::endl;
+ }
+
+
+ template <int dim>
+ void
+ AdvectionProblem<dim>::refine_grid()
+ {
+ Vector<float> gradient_indicator(triangulation.n_active_cells());
+
+ DerivativeApproximation::approximate_gradient(mapping,
+ dof_handler,
+ solution,
+ gradient_indicator);
+
+ unsigned int cell_no = 0;
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ gradient_indicator(cell_no++) *=
+ std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
+
+ GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+ gradient_indicator,
+ 0.3,
+ 0.1);
+
+ triangulation.execute_coarsening_and_refinement();
+ }
+
+
+ template <int dim>
+ void
+ AdvectionProblem<dim>::output_results(const unsigned int cycle) const
+ {
+#if false
+# ifdef HEX
+ const std::string filename =
+ dim == 2 ? ("step12-quad-" + std::to_string(cycle) + ".vtk") :
+ ("step12-hex-" + std::to_string(cycle) + ".vtk");
+# else
+ const std::string filename =
+ dim == 2 ? ("step12-tri-" + std::to_string(cycle) + ".vtk") :
+ ("step12-tet-" + std::to_string(cycle) + ".vtk");
+# endif
+ deallog << " Writing solution to <" << filename << ">" << std::endl;
+ std::ofstream output(filename);
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
+
+ data_out.build_patches(mapping, 2);
+
+ data_out.write_vtk(output);
+#endif
+
+ {
+ Vector<float> values(triangulation.n_active_cells());
+ VectorTools::integrate_difference(mapping,
+ dof_handler,
+ solution,
+ Functions::ZeroFunction<dim>(),
+ values,
+ QGauss<dim>(fe.degree + 1),
+ VectorTools::Linfty_norm);
+ const double l_infty =
+ VectorTools::compute_global_error(triangulation,
+ values,
+ VectorTools::Linfty_norm);
+ deallog << " L-infinity norm: " << l_infty << std::endl;
+ }
+ }
+
+
+ template <int dim>
+ void
+ AdvectionProblem<dim>::run()
+ {
+ //#ifdef HEX
+ // for (unsigned int cycle = 0; cycle < 6; ++cycle)
+ //#else
+ for (unsigned int cycle = 0; cycle < 1; ++cycle)
+ //#endif
+ {
+ deallog << "Cycle " << cycle << std::endl;
+
+ if (cycle == 0)
+ {
+#ifdef HEX
+ // GridGenerator::hyper_cube(triangulation);
+ // triangulation.refine_global(3);
+ GridGenerator::subdivided_hyper_cube(triangulation, 16);
+#else
+ GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+ dim == 2 ? 32 :
+ 8);
+#endif
+ }
+ else
+ refine_grid();
+
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells() << std::endl;
+
+ setup_system();
+
+ deallog << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+
+ assemble_system();
+ solve();
+
+ output_results(cycle);
+ }
+ }
+} // namespace Step12
+
+
+int
+main()
+{
+ initlog();
+
+ {
+ Step12::AdvectionProblem<2> dgmethod;
+ dgmethod.run();
+ }
+ {
+ Step12::AdvectionProblem<3> dgmethod;
+ dgmethod.run();
+ }
+
+ return 0;
+}
--- /dev/null
+
+DEAL::Cycle 0
+DEAL:: Number of active cells: 2048
+DEAL:: Number of degrees of freedom: 12288
+DEAL:Richardson::Starting value 0.0883883
+DEAL:Richardson::Convergence step 16 value 9.71637e-17
+DEAL:: Solver converged in 16 iterations.
+DEAL:: L-infinity norm: 2.33280
+DEAL::Cycle 0
+DEAL:: Number of active cells: 2560
+DEAL:: Number of degrees of freedom: 25600
+DEAL:Richardson::Starting value 0.0360844
+DEAL:Richardson::Convergence step 13 value 2.52800e-14
+DEAL:: Solver converged in 13 iterations.
+DEAL:: L-infinity norm: 13.9260
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Step-18 with tetrahedron mesh.
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/physics/transformations.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+// simplex
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+//#define HEX
+
+const unsigned int degree = 1;
+
+namespace Step18
+{
+ using namespace dealii;
+
+ template <int dim>
+ struct PointHistory
+ {
+ SymmetricTensor<2, dim> old_stress;
+ };
+
+ template <int dim>
+ SymmetricTensor<4, dim>
+ get_stress_strain_tensor(const double lambda, const double mu)
+ {
+ SymmetricTensor<4, dim> tmp;
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ for (unsigned int k = 0; k < dim; ++k)
+ for (unsigned int l = 0; l < dim; ++l)
+ tmp[i][j][k][l] = (((i == k) && (j == l) ? mu : 0.0) +
+ ((i == l) && (j == k) ? mu : 0.0) +
+ ((i == j) && (k == l) ? lambda : 0.0));
+ return tmp;
+ }
+
+ template <int dim>
+ inline SymmetricTensor<2, dim>
+ get_strain(const FEValues<dim> &fe_values,
+ const unsigned int shape_func,
+ const unsigned int q_point)
+ {
+ SymmetricTensor<2, dim> tmp;
+
+ for (unsigned int i = 0; i < dim; ++i)
+ tmp[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i];
+
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = i + 1; j < dim; ++j)
+ tmp[i][j] =
+ (fe_values.shape_grad_component(shape_func, q_point, i)[j] +
+ fe_values.shape_grad_component(shape_func, q_point, j)[i]) /
+ 2;
+
+ return tmp;
+ }
+
+
+ template <int dim>
+ inline SymmetricTensor<2, dim>
+ get_strain(const std::vector<Tensor<1, dim>> &grad)
+ {
+ Assert(grad.size() == dim, ExcInternalError());
+
+ SymmetricTensor<2, dim> strain;
+ for (unsigned int i = 0; i < dim; ++i)
+ strain[i][i] = grad[i][i];
+
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = i + 1; j < dim; ++j)
+ strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
+
+ return strain;
+ }
+
+
+ Tensor<2, 2>
+ get_rotation_matrix(const std::vector<Tensor<1, 2>> &grad_u)
+ {
+ const double curl = (grad_u[1][0] - grad_u[0][1]);
+
+ const double angle = std::atan(curl);
+
+ return Physics::Transformations::Rotations::rotation_matrix_2d(-angle);
+ }
+
+
+ Tensor<2, 3>
+ get_rotation_matrix(const std::vector<Tensor<1, 3>> &grad_u)
+ {
+ const Point<3> curl(grad_u[2][1] - grad_u[1][2],
+ grad_u[0][2] - grad_u[2][0],
+ grad_u[1][0] - grad_u[0][1]);
+
+ const double tan_angle = std::sqrt(curl * curl);
+ const double angle = std::atan(tan_angle);
+
+ if (std::abs(angle) < 1e-9)
+ {
+ static const double rotation[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+ static const Tensor<2, 3> rot(rotation);
+ return rot;
+ }
+
+ const Point<3> axis = curl / tan_angle;
+ return Physics::Transformations::Rotations::rotation_matrix_3d(axis,
+ -angle);
+ }
+
+
+
+ template <int dim>
+ class TopLevel
+ {
+ public:
+ TopLevel();
+ ~TopLevel();
+ void
+ run();
+
+ private:
+ void
+ create_coarse_grid();
+
+ void
+ setup_system();
+
+ void
+ assemble_system();
+
+ void
+ solve_timestep();
+
+ unsigned int
+ solve_linear_problem();
+
+ void
+ output_results() const;
+
+ void
+ do_initial_timestep(const bool do_output = false);
+
+ void
+ do_timestep(const bool do_output = false);
+
+ void
+ refine_initial_grid();
+
+ void
+ move_mesh();
+
+ void
+ setup_quadrature_point_history();
+
+ void
+ update_quadrature_point_history();
+
+ Triangulation<dim> triangulation;
+
+ FESystem<dim> fe;
+
+ DoFHandler<dim> dof_handler;
+
+ AffineConstraints<double> hanging_node_constraints;
+
+ const Quadrature<dim> quadrature_formula;
+
+#ifdef HEX
+ MappingQGeneric<dim, dim> mapping;
+#else
+ MappingFE<dim, dim> mapping;
+#endif
+
+ std::vector<PointHistory<dim>> quadrature_point_history;
+
+#ifdef DEAL_II_WITH_PETSC
+ PETScWrappers::MPI::SparseMatrix system_matrix;
+ PETScWrappers::MPI::Vector system_rhs;
+#else
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+ Vector<double> system_rhs;
+#endif
+
+ Vector<double> incremental_displacement;
+
+ double present_time;
+ double present_timestep;
+ double end_time;
+ unsigned int timestep_no;
+
+ MPI_Comm mpi_communicator;
+
+ const unsigned int n_mpi_processes;
+
+ const unsigned int this_mpi_process;
+
+ ConditionalOStream pcout;
+
+ IndexSet locally_owned_dofs;
+ IndexSet locally_relevant_dofs;
+
+ static const SymmetricTensor<4, dim> stress_strain_tensor;
+ };
+
+
+ template <int dim>
+ class BodyForce : public Function<dim>
+ {
+ public:
+ BodyForce();
+
+ virtual void
+ vector_value(const Point<dim> &p, Vector<double> &values) const override;
+
+ virtual void
+ vector_value_list(const std::vector<Point<dim>> &points,
+ std::vector<Vector<double>> & value_list) const override;
+ };
+
+
+ template <int dim>
+ BodyForce<dim>::BodyForce()
+ : Function<dim>(dim)
+ {}
+
+
+ template <int dim>
+ inline void
+ BodyForce<dim>::vector_value(const Point<dim> & /*p*/,
+ Vector<double> &values) const
+ {
+ Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
+
+ const double g = 9.81;
+ const double rho = 7700;
+
+ values = 0;
+ values(dim - 1) = -rho * g;
+ }
+
+
+
+ template <int dim>
+ void
+ BodyForce<dim>::vector_value_list(
+ const std::vector<Point<dim>> &points,
+ std::vector<Vector<double>> & value_list) const
+ {
+ const unsigned int n_points = points.size();
+
+ Assert(value_list.size() == n_points,
+ ExcDimensionMismatch(value_list.size(), n_points));
+
+ for (unsigned int p = 0; p < n_points; ++p)
+ BodyForce<dim>::vector_value(points[p], value_list[p]);
+ }
+
+
+
+ template <int dim>
+ class IncrementalBoundaryValues : public Function<dim>
+ {
+ public:
+ IncrementalBoundaryValues(const double present_time,
+ const double present_timestep);
+
+ virtual void
+ vector_value(const Point<dim> &p, Vector<double> &values) const override;
+
+ virtual void
+ vector_value_list(const std::vector<Point<dim>> &points,
+ std::vector<Vector<double>> & value_list) const override;
+
+ private:
+ const double velocity;
+ const double present_time;
+ const double present_timestep;
+ };
+
+
+ template <int dim>
+ IncrementalBoundaryValues<dim>::IncrementalBoundaryValues(
+ const double present_time,
+ const double present_timestep)
+ : Function<dim>(dim)
+ , velocity(.08)
+ , present_time(present_time)
+ , present_timestep(present_timestep)
+ {}
+
+
+ template <int dim>
+ void
+ IncrementalBoundaryValues<dim>::vector_value(const Point<dim> & /*p*/,
+ Vector<double> &values) const
+ {
+ Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
+
+ values = 0;
+ values(2) = -present_timestep * velocity;
+ }
+
+
+
+ template <int dim>
+ void
+ IncrementalBoundaryValues<dim>::vector_value_list(
+ const std::vector<Point<dim>> &points,
+ std::vector<Vector<double>> & value_list) const
+ {
+ const unsigned int n_points = points.size();
+
+ Assert(value_list.size() == n_points,
+ ExcDimensionMismatch(value_list.size(), n_points));
+
+ for (unsigned int p = 0; p < n_points; ++p)
+ IncrementalBoundaryValues<dim>::vector_value(points[p], value_list[p]);
+ }
+
+
+ template <int dim>
+ const SymmetricTensor<4, dim> TopLevel<dim>::stress_strain_tensor =
+ get_stress_strain_tensor<dim>(/*lambda = */ 9.695e10,
+ /*mu = */ 7.617e10);
+
+
+#ifdef HEX
+ template <int dim>
+ TopLevel<dim>::TopLevel()
+ : triangulation()
+ , fe(FE_Q<dim>(degree), dim)
+ , dof_handler(triangulation)
+ , quadrature_formula(QGauss<dim>(fe.degree + 1))
+ , mapping(1)
+ , present_time(0.0)
+ , present_timestep(1.0)
+ , end_time(10.0)
+ , timestep_no(0)
+ , mpi_communicator(MPI_COMM_WORLD)
+ , n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
+ , this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
+ , pcout(std::cout, this_mpi_process == 0)
+ {}
+#else
+ template <int dim>
+ TopLevel<dim>::TopLevel()
+ : triangulation()
+ , fe(Simplex::FE_P<dim>(degree), dim)
+ , dof_handler(triangulation)
+ , quadrature_formula(Simplex::PGauss<dim>(fe.degree == 1 ? 4 : 10))
+ , mapping(Simplex::FE_P<dim>(1))
+ , present_time(0.0)
+ , present_timestep(1.0)
+ , end_time(10.0)
+ , timestep_no(0)
+ , mpi_communicator(MPI_COMM_WORLD)
+ , n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
+ , this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
+ , pcout(std::cout, this_mpi_process == 0)
+ {}
+#endif
+
+
+
+ template <int dim>
+ TopLevel<dim>::~TopLevel()
+ {
+ dof_handler.clear();
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::run()
+ {
+ do_initial_timestep(false);
+
+ while (present_time < end_time)
+ do_timestep(std::abs(end_time - present_time - present_timestep) < 10e-5);
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::create_coarse_grid()
+ {
+ const unsigned int n = 5;
+
+#ifdef HEX
+ GridGenerator::subdivided_hyper_rectangle(triangulation,
+ {1 * n, 1 * n, 3 * n},
+ {-0.5, -0.5, 0},
+ {+0.5, +0.5, +3});
+#else
+ GridGenerator::subdivided_hyper_rectangle_with_simplices(
+ triangulation, {1 * n, 1 * n, 3 * n}, {-0.5, -0.5, 0}, {+0.5, +0.5, +3});
+#endif
+
+ for (const auto &cell : triangulation.active_cell_iterators())
+ for (const auto &face : cell->face_iterators())
+ if (face->at_boundary())
+ {
+ const Point<dim> face_center = face->center();
+
+ if (face_center[2] == 0)
+ face->set_boundary_id(0);
+ else if (face_center[2] == 3)
+ face->set_boundary_id(1);
+ else
+ face->set_boundary_id(2);
+ }
+
+ setup_quadrature_point_history();
+ }
+
+ template <int dim>
+ void
+ TopLevel<dim>::setup_system()
+ {
+ dof_handler.distribute_dofs(fe);
+ locally_owned_dofs = dof_handler.locally_owned_dofs();
+ DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+ // The next step is to set up constraints due to hanging nodes. This has
+ // been handled many times before:
+ hanging_node_constraints.clear();
+ DoFTools::make_hanging_node_constraints(dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close();
+
+#ifdef DEAL_II_WITH_PETSC
+ DynamicSparsityPattern sparsity_pattern(locally_relevant_dofs);
+ DoFTools::make_sparsity_pattern(dof_handler,
+ sparsity_pattern,
+ hanging_node_constraints,
+ /*keep constrained dofs*/ false);
+ SparsityTools::distribute_sparsity_pattern(sparsity_pattern,
+ locally_owned_dofs,
+ mpi_communicator,
+ locally_relevant_dofs);
+
+ system_matrix.reinit(locally_owned_dofs,
+ locally_owned_dofs,
+ sparsity_pattern,
+ mpi_communicator);
+
+ system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+ incremental_displacement.reinit(dof_handler.n_dofs());
+
+#else
+ // DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ // DoFTools::make_sparsity_pattern(dof_handler, dsp,
+ // hanging_node_constraints, false);
+ // sparsity_pattern.copy_from(dsp);
+ // system_matrix.reinit(sparsity_pattern);
+ // sparsity_pattern.copy_from(dsp);
+ // system_matrix.reinit(sparsity_pattern);
+
+ DynamicSparsityPattern dsp(locally_relevant_dofs);
+ DoFTools::make_sparsity_pattern(dof_handler,
+ dsp,
+ hanging_node_constraints,
+ /*keep constrained dofs*/ false);
+ SparsityTools::distribute_sparsity_pattern(dsp,
+ locally_owned_dofs,
+ mpi_communicator,
+ locally_relevant_dofs);
+ sparsity_pattern.copy_from(dsp);
+ system_matrix.reinit(sparsity_pattern);
+
+ system_rhs.reinit(dof_handler.n_dofs());
+ incremental_displacement.reinit(dof_handler.n_dofs());
+#endif
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::assemble_system()
+ {
+ system_rhs = 0;
+ system_matrix = 0;
+
+ FEValues<dim> fe_values(mapping,
+ fe,
+ quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+ BodyForce<dim> body_force;
+ std::vector<Vector<double>> body_force_values(n_q_points,
+ Vector<double>(dim));
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ fe_values.reinit(cell);
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ const SymmetricTensor<2, dim>
+ eps_phi_i = get_strain(fe_values, i, q_point),
+ eps_phi_j = get_strain(fe_values, j, q_point);
+
+ cell_matrix(i, j) += (eps_phi_i * //
+ stress_strain_tensor * //
+ eps_phi_j //
+ ) * //
+ fe_values.JxW(q_point); //
+ }
+
+
+ const PointHistory<dim> *local_quadrature_points_data =
+ reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+
+ body_force.vector_value_list(fe_values.get_quadrature_points(),
+ body_force_values);
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ const SymmetricTensor<2, dim> &old_stress =
+ local_quadrature_points_data[q_point].old_stress;
+
+ cell_rhs(i) +=
+ (body_force_values[q_point](component_i) *
+ fe_values.shape_value(i, q_point) -
+ old_stress * get_strain(fe_values, i, q_point)) *
+ fe_values.JxW(q_point);
+ }
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+
+ hanging_node_constraints.distribute_local_to_global(cell_matrix,
+ cell_rhs,
+ local_dof_indices,
+ system_matrix,
+ system_rhs);
+ }
+
+ // Now compress the vector and the system matrix:
+#ifdef DEAL_II_WITH_PETSC
+ system_matrix.compress(VectorOperation::add);
+ system_rhs.compress(VectorOperation::add);
+#endif
+
+
+ FEValuesExtractors::Scalar z_component(dim - 1);
+ std::map<types::global_dof_index, double> boundary_values;
+ VectorTools::interpolate_boundary_values(mapping,
+ dof_handler,
+ 0,
+ Functions::ZeroFunction<dim>(dim),
+ boundary_values);
+ VectorTools::interpolate_boundary_values(
+ mapping,
+ dof_handler,
+ 1,
+ IncrementalBoundaryValues<dim>(present_time, present_timestep),
+ boundary_values,
+ fe.component_mask(z_component));
+
+#ifdef DEAL_II_WITH_PETSC
+ PETScWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);
+#else
+ Vector<double> tmp(dof_handler.n_dofs());
+#endif
+ MatrixTools::apply_boundary_values(
+ boundary_values, system_matrix, tmp, system_rhs, false);
+ incremental_displacement = tmp;
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::solve_timestep()
+ {
+ deallog << " Assembling system..." << std::flush;
+ assemble_system();
+ deallog << " norm of rhs is " << system_rhs.l2_norm() << std::endl;
+
+ const unsigned int n_iterations = solve_linear_problem();
+
+ deallog << " Solver converged in " << n_iterations << " iterations."
+ << std::endl;
+
+ deallog << " Updating quadrature point data..." << std::flush;
+ update_quadrature_point_history();
+ deallog << std::endl;
+ }
+
+
+ template <int dim>
+ unsigned int
+ TopLevel<dim>::solve_linear_problem()
+ {
+#ifdef DEAL_II_WITH_PETSC
+ PETScWrappers::MPI::Vector distributed_incremental_displacement(
+ locally_owned_dofs, mpi_communicator);
+ distributed_incremental_displacement = incremental_displacement;
+#else
+ Vector<double> distributed_incremental_displacement(dof_handler.n_dofs());
+ distributed_incremental_displacement = incremental_displacement;
+#endif
+
+ SolverControl solver_control(dof_handler.n_dofs(),
+ 1e-16 * system_rhs.l2_norm());
+
+#ifdef DEAL_II_WITH_PETSC
+ PETScWrappers::SolverCG cg(solver_control, mpi_communicator);
+
+ PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
+
+ cg.solve(system_matrix,
+ distributed_incremental_displacement,
+ system_rhs,
+ preconditioner);
+#else
+ SolverCG<Vector<double>> solver(solver_control);
+ solver.solve(system_matrix,
+ distributed_incremental_displacement,
+ system_rhs,
+ PreconditionIdentity());
+#endif
+
+ deallog << "norm: " << distributed_incremental_displacement.linfty_norm()
+ << " " << distributed_incremental_displacement.l1_norm() << " "
+ << distributed_incremental_displacement.l2_norm() << std::endl;
+
+ incremental_displacement = distributed_incremental_displacement;
+
+ hanging_node_constraints.distribute(incremental_displacement);
+
+ return solver_control.last_step();
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::output_results() const
+ {
+ return;
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+
+ std::vector<std::string> solution_names;
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ solution_interpretation;
+
+ solution_names.assign(dim, "delta");
+ solution_interpretation.assign(
+ dim,
+ DataComponentInterpretation::DataComponentInterpretation::
+ component_is_part_of_vector);
+
+ data_out.add_data_vector(incremental_displacement,
+ solution_names,
+ DataOut_DoFData<DoFHandler<dim, dim>, dim, dim>::
+ DataVectorType::type_automatic,
+ solution_interpretation);
+
+
+ Vector<double> norm_of_stress(triangulation.n_active_cells());
+ {
+ // Loop over all the cells...
+ for (auto &cell : triangulation.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ // On these cells, add up the stresses over all quadrature
+ // points...
+ SymmetricTensor<2, dim> accumulated_stress;
+ for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
+ accumulated_stress +=
+ reinterpret_cast<PointHistory<dim> *>(cell->user_pointer())[q]
+ .old_stress;
+
+ // ...then write the norm of the average to their destination:
+ norm_of_stress(cell->active_cell_index()) =
+ (accumulated_stress / quadrature_formula.size()).norm();
+ }
+ else
+ norm_of_stress(cell->active_cell_index()) = -1e+20;
+ }
+
+ data_out.add_data_vector(norm_of_stress, "norm_of_stress");
+
+ std::vector<types::subdomain_id> partition_int(
+ triangulation.n_active_cells());
+ GridTools::get_subdomain_association(triangulation, partition_int);
+ const Vector<double> partitioning(partition_int.begin(),
+ partition_int.end());
+ data_out.add_data_vector(partitioning, "partitioning");
+
+ data_out.build_patches(mapping, 2);
+
+#if false
+ std::ofstream output("step18." + std::to_string(this_mpi_process) + "." +
+ std::to_string(timestep_no) + ".vtk");
+ data_out.write_vtk(output);
+#endif
+ }
+
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::do_initial_timestep(const bool do_output)
+ {
+ present_time += present_timestep;
+ ++timestep_no;
+ deallog << "Timestep " << timestep_no << " at time " << present_time
+ << std::endl;
+
+ for (unsigned int cycle = 0; cycle < 1; ++cycle)
+ {
+ deallog << " Cycle " << cycle << ':' << std::endl;
+
+ if (cycle == 0)
+ create_coarse_grid();
+ // else
+ // refine_initial_grid();
+
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells() << " (by partition:";
+ for (unsigned int p = 0; p < n_mpi_processes; ++p)
+ deallog << (p == 0 ? ' ' : '+')
+ << (GridTools::count_cells_with_subdomain_association(
+ triangulation, p));
+ deallog << ")" << std::endl;
+
+ setup_system();
+
+ deallog << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ << " (by partition:";
+ for (unsigned int p = 0; p < n_mpi_processes; ++p)
+ deallog << (p == 0 ? ' ' : '+')
+ << (DoFTools::count_dofs_with_subdomain_association(
+ dof_handler, p));
+ deallog << ")" << std::endl;
+
+ solve_timestep();
+ }
+
+ move_mesh();
+
+ if (do_output)
+ output_results();
+
+ deallog << std::endl;
+ }
+
+ template <int dim>
+ void
+ TopLevel<dim>::do_timestep(const bool do_output)
+ {
+ present_time += present_timestep;
+ ++timestep_no;
+ deallog << "Timestep " << timestep_no << " at time " << present_time
+ << std::endl;
+ if (present_time > end_time)
+ {
+ present_timestep -= (present_time - end_time);
+ present_time = end_time;
+ }
+
+
+ solve_timestep();
+
+ move_mesh();
+
+ if (do_output)
+ output_results();
+
+ deallog << std::endl;
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::refine_initial_grid()
+ {
+ // First, let each process compute error indicators for the cells it owns:
+ Vector<float> error_per_cell(triangulation.n_active_cells());
+ KellyErrorEstimator<dim>::estimate(
+ dof_handler,
+ QGauss<dim - 1>(fe.degree + 1),
+ std::map<types::boundary_id, const Function<dim> *>(),
+ incremental_displacement,
+ error_per_cell,
+ ComponentMask(),
+ nullptr,
+ MultithreadInfo::n_threads(),
+ this_mpi_process);
+
+ const unsigned int n_local_cells = triangulation.n_active_cells();
+
+#ifdef DEAL_II_WITH_PETSC
+ PETScWrappers::MPI::Vector distributed_error_per_cell(
+ mpi_communicator, triangulation.n_active_cells(), n_local_cells);
+#else
+ Vector<double> distributed_error_per_cell(n_local_cells);
+#endif
+
+ for (unsigned int i = 0; i < error_per_cell.size(); ++i)
+ if (error_per_cell(i) != 0)
+ distributed_error_per_cell(i) = error_per_cell(i);
+ distributed_error_per_cell.compress(VectorOperation::insert);
+
+ error_per_cell = distributed_error_per_cell;
+ GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+ error_per_cell,
+ 0.35,
+ 0.03);
+ triangulation.execute_coarsening_and_refinement();
+
+ setup_quadrature_point_history();
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::move_mesh()
+ {
+ deallog << " Moving mesh..." << std::endl;
+
+ std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
+ for (auto &cell : dof_handler.active_cell_iterators())
+ for (unsigned int v = 0; v < cell->n_vertices(); ++v)
+ if (vertex_touched[cell->vertex_index(v)] == false)
+ {
+ vertex_touched[cell->vertex_index(v)] = true;
+
+ Point<dim> vertex_displacement;
+ for (unsigned int d = 0; d < dim; ++d)
+ vertex_displacement[d] =
+ incremental_displacement(cell->vertex_dof_index(v, d));
+
+ cell->vertex(v) += vertex_displacement;
+ }
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::setup_quadrature_point_history()
+ {
+ triangulation.clear_user_data();
+
+ {
+ std::vector<PointHistory<dim>> tmp;
+ quadrature_point_history.swap(tmp);
+ }
+ quadrature_point_history.resize(triangulation.n_active_cells() *
+ quadrature_formula.size());
+
+ unsigned int history_index = 0;
+ for (auto &cell : triangulation.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ cell->set_user_pointer(&quadrature_point_history[history_index]);
+ history_index += quadrature_formula.size();
+ }
+
+ Assert(history_index == quadrature_point_history.size(),
+ ExcInternalError());
+ }
+
+
+ template <int dim>
+ void
+ TopLevel<dim>::update_quadrature_point_history()
+ {
+ FEValues<dim> fe_values(mapping,
+ fe,
+ quadrature_formula,
+ update_values | update_gradients);
+
+ std::vector<std::vector<Tensor<1, dim>>> displacement_increment_grads(
+ quadrature_formula.size(), std::vector<Tensor<1, dim>>(dim));
+
+ for (auto &cell : dof_handler.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ PointHistory<dim> *local_quadrature_points_history =
+ reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+ Assert(local_quadrature_points_history >=
+ &quadrature_point_history.front(),
+ ExcInternalError());
+ Assert(local_quadrature_points_history <=
+ &quadrature_point_history.back(),
+ ExcInternalError());
+
+ fe_values.reinit(cell);
+ fe_values.get_function_gradients(incremental_displacement,
+ displacement_increment_grads);
+
+ for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
+ {
+ const SymmetricTensor<2, dim> new_stress =
+ (local_quadrature_points_history[q].old_stress +
+ (stress_strain_tensor *
+ get_strain(displacement_increment_grads[q])));
+
+ const Tensor<2, dim> rotation =
+ get_rotation_matrix(displacement_increment_grads[q]);
+
+ const SymmetricTensor<2, dim> rotated_new_stress =
+ symmetrize(transpose(rotation) *
+ static_cast<Tensor<2, dim>>(new_stress) * rotation);
+
+ local_quadrature_points_history[q].old_stress =
+ rotated_new_stress;
+ }
+ }
+ }
+} // namespace Step18
+
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ try
+ {
+ using namespace dealii;
+ using namespace Step18;
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ TopLevel<3> elastic_problem;
+ elastic_problem.run();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
--- /dev/null
+
+DEAL::Timestep 1 at time 1.00000
+DEAL:: Cycle 0:
+DEAL:: Number of active cells: 1875 (by partition: 1875)
+DEAL:: Number of degrees of freedom: 1728 (by partition: 1728)
+DEAL:: Assembling system... norm of rhs is 2.45750e+10
+DEAL:cg::Starting value 1.95453e+10
+DEAL:cg::Convergence step 126 value 2.13799e-06
+DEAL::norm: 0.0800000 25.2034 1.12390
+DEAL:: Solver converged in 126 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 2 at time 2.00000
+DEAL:: Assembling system... norm of rhs is 2.49844e+10
+DEAL:cg::Starting value 2.03490e+10
+DEAL:cg::Convergence step 127 value 1.74037e-06
+DEAL::norm: 0.0800000 25.3187 1.12498
+DEAL:: Solver converged in 127 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 3 at time 3.00000
+DEAL:: Assembling system... norm of rhs is 2.54414e+10
+DEAL:cg::Starting value 2.12110e+10
+DEAL:cg::Convergence step 127 value 2.31040e-06
+DEAL::norm: 0.0800000 25.4405 1.12611
+DEAL:: Solver converged in 127 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 4 at time 4.00000
+DEAL:: Assembling system... norm of rhs is 2.59506e+10
+DEAL:cg::Starting value 2.21374e+10
+DEAL:cg::Convergence step 127 value 2.50495e-06
+DEAL::norm: 0.0800000 25.5693 1.12729
+DEAL:: Solver converged in 127 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 5 at time 5.00000
+DEAL:: Assembling system... norm of rhs is 2.65175e+10
+DEAL:cg::Starting value 2.31346e+10
+DEAL:cg::Convergence step 128 value 2.04181e-06
+DEAL::norm: 0.0800000 25.7059 1.12851
+DEAL:: Solver converged in 128 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 6 at time 6.00000
+DEAL:: Assembling system... norm of rhs is 2.71482e+10
+DEAL:cg::Starting value 2.42104e+10
+DEAL:cg::Convergence step 129 value 2.29400e-06
+DEAL::norm: 0.0800000 25.8512 1.12979
+DEAL:: Solver converged in 129 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 7 at time 7.00000
+DEAL:: Assembling system... norm of rhs is 2.78496e+10
+DEAL:cg::Starting value 2.53733e+10
+DEAL:cg::Convergence step 131 value 2.19345e-06
+DEAL::norm: 0.0800000 26.0062 1.13114
+DEAL:: Solver converged in 131 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 8 at time 8.00000
+DEAL:: Assembling system... norm of rhs is 2.86298e+10
+DEAL:cg::Starting value 2.66332e+10
+DEAL:cg::Convergence step 132 value 2.63526e-06
+DEAL::norm: 0.0800000 26.1724 1.13256
+DEAL:: Solver converged in 132 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 9 at time 9.00000
+DEAL:: Assembling system... norm of rhs is 2.94979e+10
+DEAL:cg::Starting value 2.80016e+10
+DEAL:cg::Convergence step 135 value 2.78060e-06
+DEAL::norm: 0.0800000 26.3512 1.13406
+DEAL:: Solver converged in 135 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 10 at time 10.0000
+DEAL:: Assembling system... norm of rhs is 3.04646e+10
+DEAL:cg::Starting value 2.94916e+10
+DEAL:cg::Convergence step 137 value 2.61758e-06
+DEAL::norm: 0.0800000 26.5447 1.13566
+DEAL:: Solver converged in 137 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
--- /dev/null
+
+DEAL::Timestep 1 at time 1.00000
+DEAL:: Cycle 0:
+DEAL:: Number of active cells: 1875 (by partition: 1875)
+DEAL:: Number of degrees of freedom: 1728 (by partition: 1728)
+DEAL:: Assembling system... norm of rhs is 6.42608e+09
+DEAL::Starting value 0.406456
+DEAL::Convergence step 49 value 2.89749e-07
+DEAL::norm: 0.0800000 25.2034 1.12390
+DEAL:: Solver converged in 49 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 2 at time 2.00000
+DEAL:: Assembling system... norm of rhs is 6.38314e+09
+DEAL::Starting value 0.413459
+DEAL::Convergence step 48 value 2.33590e-07
+DEAL::norm: 0.0800000 25.3187 1.12498
+DEAL:: Solver converged in 48 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 3 at time 3.00000
+DEAL:: Assembling system... norm of rhs is 6.34811e+09
+DEAL::Starting value 0.420747
+DEAL::Convergence step 47 value 2.56013e-07
+DEAL::norm: 0.0800000 25.4405 1.12611
+DEAL:: Solver converged in 47 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 4 at time 4.00000
+DEAL:: Assembling system... norm of rhs is 6.32248e+09
+DEAL::Starting value 0.428346
+DEAL::Convergence step 45 value 6.05128e-07
+DEAL::norm: 0.0800000 25.5693 1.12729
+DEAL:: Solver converged in 45 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 5 at time 5.00000
+DEAL:: Assembling system... norm of rhs is 6.30808e+09
+DEAL::Starting value 0.436284
+DEAL::Convergence step 45 value 1.99224e-07
+DEAL::norm: 0.0800000 25.7059 1.12851
+DEAL:: Solver converged in 45 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 6 at time 6.00000
+DEAL:: Assembling system... norm of rhs is 6.30720e+09
+DEAL::Starting value 0.444594
+DEAL::Convergence step 43 value 5.57039e-07
+DEAL::norm: 0.0800000 25.8512 1.12979
+DEAL:: Solver converged in 43 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 7 at time 7.00000
+DEAL:: Assembling system... norm of rhs is 6.32271e+09
+DEAL::Starting value 0.453310
+DEAL::Convergence step 42 value 3.35520e-07
+DEAL::norm: 0.0800000 26.0062 1.13114
+DEAL:: Solver converged in 42 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 8 at time 8.00000
+DEAL:: Assembling system... norm of rhs is 6.35827e+09
+DEAL::Starting value 0.462469
+DEAL::Convergence step 41 value 2.38274e-07
+DEAL::norm: 0.0800000 26.1724 1.13256
+DEAL:: Solver converged in 41 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 9 at time 9.00000
+DEAL:: Assembling system... norm of rhs is 6.41862e+09
+DEAL::Starting value 0.472114
+DEAL::Convergence step 40 value 2.20795e-07
+DEAL::norm: 0.0800000 26.3512 1.13406
+DEAL:: Solver converged in 40 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
+DEAL::Timestep 10 at time 10.0000
+DEAL:: Assembling system... norm of rhs is 6.50998e+09
+DEAL::Starting value 0.482291
+DEAL::Convergence step 38 value 4.09099e-07
+DEAL::norm: 0.0800000 26.5447 1.13566
+DEAL:: Solver converged in 38 iterations.
+DEAL:: Updating quadrature point data...
+DEAL:: Moving mesh...
+DEAL::
projects / dealii.git / commitdiff