From 5cdac2e064ba4b488bfae6fda25201bf0384b4fc Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 28 Jul 2020 11:59:26 +0200 Subject: [PATCH] Add mini simplex application tests --- source/base/qprojector.cc | 10 +- tests/simplex/data_out_write_vtk_01.cc | 2 +- tests/simplex/poisson_01.cc | 455 +++++++ ...linos=false.with_simplex_support=on.output | 25 + ...ilinos=true.with_simplex_support=on.output | 25 + ...ilinos=true.with_simplex_support=on.output | 25 + tests/simplex/poisson_02.cc | 670 +++++++++++ .../poisson_02.with_simplex_support=on.output | 37 + tests/simplex/step-12.cc | 563 +++++++++ .../step-12.with_simplex_support=on.output | 15 + tests/simplex/step-18.cc | 1053 +++++++++++++++++ ...petsc=false.with_simplex_support=on.output | 94 ++ ..._petsc=true.with_simplex_support=on.output | 94 ++ 13 files changed, 3065 insertions(+), 3 deletions(-) create mode 100644 tests/simplex/poisson_01.cc create mode 100644 tests/simplex/poisson_01.mpirun=1.with_trilinos=false.with_simplex_support=on.output create mode 100644 tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output create mode 100644 tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output create mode 100644 tests/simplex/poisson_02.cc create mode 100644 tests/simplex/poisson_02.with_simplex_support=on.output create mode 100644 tests/simplex/step-12.cc create mode 100644 tests/simplex/step-12.with_simplex_support=on.output create mode 100644 tests/simplex/step-18.cc create mode 100644 tests/simplex/step-18.mpirun=1.with_petsc=false.with_simplex_support=on.output create mode 100644 tests/simplex/step-18.mpirun=1.with_petsc=true.with_simplex_support=on.output diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index d9c87099c0..e3e7b4e133 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -906,8 +906,11 @@ QProjector<2>::project_to_all_subfaces( 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; @@ -964,8 +967,11 @@ QProjector<3>::project_to_all_subfaces( 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); diff --git a/tests/simplex/data_out_write_vtk_01.cc b/tests/simplex/data_out_write_vtk_01.cc index 027654e519..f830e658a2 100644 --- a/tests/simplex/data_out_write_vtk_01.cc +++ b/tests/simplex/data_out_write_vtk_01.cc @@ -56,7 +56,7 @@ void test(const FiniteElement &fe, const unsigned int n_components) { Triangulation tria; - Simplex::GridGenerator::subdivided_hyper_cube(tria, dim == 2 ? 4 : 2); + GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 4 : 2); DoFHandler dof_handler(tria); diff --git a/tests/simplex/poisson_01.cc b/tests/simplex/poisson_01.cc new file mode 100644 index 0000000000..84410336fb --- /dev/null +++ b/tests/simplex/poisson_01.cc @@ -0,0 +1,455 @@ +// --------------------------------------------------------------------- +// +// 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 +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +struct Parameters +{ + unsigned int degree = 2; + + + // GridGenerator + bool use_grid_generator = true; + std::vector repetitions; + Point p1; + Point p2; + + // GridIn + std::string file_name_in = ""; + + // GridOut + std::string file_name_out = ""; +}; + +template +MPI_Comm +get_communicator(const Triangulation &tria) +{ + if (auto tria_ = + dynamic_cast *>(&tria)) + return tria_->get_communicator(); + + return MPI_COMM_SELF; +} + +template +void +test(const Triangulation &tria, + const FiniteElement &fe, + const Quadrature & quad, + const Quadrature & face_quad, + const Mapping & mapping, + const double r_boundary) +{ + std::string label = + (dynamic_cast *>( + &tria) ? + "parallel::shared::Triangulation" : + (dynamic_cast< + const parallel::fullydistributed::Triangulation *>( + &tria) ? + "parallel::fullydistributed::Triangulation" : + (dynamic_cast< + const parallel::distributed::Triangulation *>( + &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 dof_handler(tria); + dof_handler.distribute_dofs(fe); + + AffineConstraints 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; + TrilinosWrappers::SparseMatrix system_matrix; + VectorType solution; + VectorType system_rhs; +#else + using VectorType = Vector; + SparsityPattern sparsity_pattern; + SparseMatrix 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 fe_values(mapping, fe, quad, flag); + + std::shared_ptr> fe_face_values; + + fe_face_values.reset( + new FEFaceValues(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 dof_indices(dofs_per_cell); + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector 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 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::vertices_per_cell); + + deallog << std::endl; +} + +template +void +test_tet(const MPI_Comm &comm, const Parameters ¶ms) +{ + const unsigned int tria_type = 2; + + // 1) Create triangulation... + Triangulation *tria; + + // a) serial triangulation + Triangulation tr_1; + + // b) shared triangulation (with artificial cells) + parallel::shared::Triangulation tr_2( + MPI_COMM_WORLD, + ::Triangulation::none, + true, + parallel::shared::Triangulation::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 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 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 fe(params.degree); + + Simplex::PGauss quad(dim == 2 ? (params.degree == 1 ? 3 : 7) : + (params.degree == 1 ? 4 : 10)); + + Simplex::PGauss face_quad(dim == 2 ? (params.degree == 1 ? 2 : 3) : + (params.degree == 1 ? 3 : 7)); + + Simplex::FE_P fe_mapping(1); + MappingFE mapping(fe_mapping); + + // 4) Perform test (independent of mesh type) + test(*tria, fe, quad, face_quad, mapping, params.p2[0]); +} + +template +void +test_hex(const MPI_Comm &comm, const Parameters ¶ms) +{ + // 1) Create triangulation... + parallel::distributed::Triangulation 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 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 fe(params.degree); + + QGauss quad(params.degree + 1); + + QGauss quad_face(params.degree + 1); + + MappingQ 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{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{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); + } + } +} diff --git a/tests/simplex/poisson_01.mpirun=1.with_trilinos=false.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=1.with_trilinos=false.with_simplex_support=on.output new file mode 100644 index 0000000000..02c1cc8169 --- /dev/null +++ b/tests/simplex/poisson_01.mpirun=1.with_trilinos=false.with_simplex_support=on.output @@ -0,0 +1,25 @@ + +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:: diff --git a/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output new file mode 100644 index 0000000000..a45012d350 --- /dev/null +++ b/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output @@ -0,0 +1,25 @@ + +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:: diff --git a/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output new file mode 100644 index 0000000000..0e41dea463 --- /dev/null +++ b/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output @@ -0,0 +1,25 @@ + +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:: diff --git a/tests/simplex/poisson_02.cc b/tests/simplex/poisson_02.cc new file mode 100644 index 0000000000..03110d8288 --- /dev/null +++ b/tests/simplex/poisson_02.cc @@ -0,0 +1,670 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +//#define HEX + +using namespace dealii; + +template +struct ScratchData +{ + ScratchData(const Mapping & mapping, + const FiniteElement & fe, + const Quadrature & quad, + const Quadrature &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 &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 fe_values; + FEInterfaceValues fe_interface_values; +}; + + + +struct CopyDataFace +{ + FullMatrix cell_matrix; + std::vector joint_dof_indices; +}; + + + +struct CopyData +{ + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + + template + 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 +class RightHandSideFunction : public Function +{ +public: + RightHandSideFunction() + {} + + virtual double + value(const Point &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 +class DGHeat +{ +public: + DGHeat(const bool hex, + FiniteElement * fe, + Mapping * mapping, + Quadrature * quad, + Quadrature *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> + HEX(unsigned int degree, + unsigned int initial_refinement, + unsigned int number_refinement) + { + return std::make_unique>(true, + new FE_DGQ(degree), + new MappingQ(1), + new QGauss(degree + 1), + new QGauss(degree + 1), + initial_refinement, + number_refinement); + } + + static std::unique_ptr> + TET(unsigned int degree, + unsigned int initial_refinement, + unsigned int number_refinement) + { + return std::make_unique>( + false, + new Simplex::FE_DGP(degree), + new MappingFE(Simplex::FE_P(1)), + new Simplex::PGauss(dim == 2 ? (degree == 1 ? 3 : 7) : + (degree == 1 ? 4 : 10)), + new Simplex::PGauss(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 triangulation; + + bool hex; + + std::unique_ptr> fe; + const std::unique_ptr> mapping; + std::unique_ptr> quad; + std::unique_ptr> face_quad; + + DoFHandler dof_handler; + + + RightHandSideFunction right_hand_side; + + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + Point center; + + ConvergenceTable error_table; + + unsigned int initial_refinement_level; + unsigned int number_refinement; +}; + + +template +void +DGHeat::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 +void +DGHeat::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 +void +DGHeat::assemble_system() +{ + using Iterator = typename DoFHandler::active_cell_iterator; + + auto cell_worker = [&](const Iterator & cell, + ScratchData &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 & fe_v = scratch_data.fe_values; + const std::vector &JxW = fe_v.get_JxW_values(); + + std::vector 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 & scratch_data, + CopyData & copy_data) { + scratch_data.fe_interface_values.reinit(cell, face_no); + + const FEFaceValuesBase &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 &JxW = fe_face.get_JxW_values(); + + const std::vector> &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 & scratch_data, + CopyData & copy_data) { + FEInterfaceValues &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 & JxW = fe_iv.get_JxW_values(); + const std::vector> &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 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 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 +void +DGHeat::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 +void +DGHeat::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 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 +void +DGHeat::calculateL2Error() +{ + FEValues 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 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 +void +DGHeat::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; +} diff --git a/tests/simplex/poisson_02.with_simplex_support=on.output b/tests/simplex/poisson_02.with_simplex_support=on.output new file mode 100644 index 0000000000..825396daad --- /dev/null +++ b/tests/simplex/poisson_02.with_simplex_support=on.output @@ -0,0 +1,37 @@ + + 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:: diff --git a/tests/simplex/step-12.cc b/tests/simplex/step-12.cc new file mode 100644 index 0000000000..e1f12330b6 --- /dev/null +++ b/tests/simplex/step-12.cc @@ -0,0 +1,563 @@ +// --------------------------------------------------------------------- +// +// 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 +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +// Finally, the new include file for using the mesh_loop from the MeshWorker +// framework +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +//#define HEX + +namespace Step12 +{ + using namespace dealii; + + template + class BoundaryValues : public Function + { + public: + BoundaryValues() = default; + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const override; + }; + + template + void + BoundaryValues::value_list(const std::vector> &points, + std::vector & 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 + Tensor<1, dim> + beta(const Point &p) + { + Assert(dim >= 2, ExcNotImplemented()); + + Point 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 + struct ScratchData + { + ScratchData(const Mapping & mapping, + const FiniteElement & fe, + const Quadrature & quad, + const Quadrature &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 &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 fe_values; + FEInterfaceValues fe_interface_values; + }; + + + + struct CopyDataFace + { + FullMatrix cell_matrix; + std::vector joint_dof_indices; + }; + + + + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + + template + 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 + 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 triangulation; +#ifdef HEX + const MappingQ1 mapping; +#else + Simplex::FE_P fe_mapping; + const MappingFE mapping; +#endif + + // Furthermore we want to use DG elements. +#ifdef HEX + FE_DGQ fe; +#else + Simplex::FE_DGP fe; +#endif + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector right_hand_side; + }; + + + template + AdvectionProblem::AdvectionProblem() +#ifdef HEX + : mapping() + , fe(2) +#else + : fe_mapping(1) + , mapping(fe_mapping) + , fe(2) +#endif + , dof_handler(triangulation) + {} + + + template + void + AdvectionProblem::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 + void + AdvectionProblem::assemble_system() + { + using Iterator = typename DoFHandler::active_cell_iterator; + const BoundaryValues boundary_function; + + auto cell_worker = [&](const Iterator & cell, + ScratchData &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 & fe_v = scratch_data.fe_values; + const std::vector &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 & scratch_data, + CopyData & copy_data) { + scratch_data.fe_interface_values.reinit(cell, face_no); + const FEFaceValuesBase &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 & JxW = fe_face.get_JxW_values(); + const std::vector> &normals = fe_face.get_normal_vectors(); + + std::vector 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 & scratch_data, + CopyData & copy_data) { + FEInterfaceValues &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 & JxW = fe_iv.get_JxW_values(); + const std::vector> &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 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 quad(degree + 1); + + QGauss face_quad(degree + 1); +#else + Simplex::PGauss quad(dim == 2 ? (degree == 1 ? 3 : 7) : + (degree == 1 ? 4 : 10)); + + Simplex::PGauss face_quad(dim == 2 ? (degree == 1 ? 2 : 3) : + (degree == 1 ? 3 : 7)); +#endif + + ScratchData 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 + void + AdvectionProblem::solve() + { + SolverControl solver_control(1000, 1e-12); + SolverRichardson> solver(solver_control); + + PreconditionBlockSSOR> 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 + void + AdvectionProblem::refine_grid() + { + Vector 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 + void + AdvectionProblem::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 data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "u", DataOut::type_dof_data); + + data_out.build_patches(mapping, 2); + + data_out.write_vtk(output); +#endif + + { + Vector values(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Functions::ZeroFunction(), + values, + QGauss(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 + void + AdvectionProblem::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; +} diff --git a/tests/simplex/step-12.with_simplex_support=on.output b/tests/simplex/step-12.with_simplex_support=on.output new file mode 100644 index 0000000000..7fbe5eb0f7 --- /dev/null +++ b/tests/simplex/step-12.with_simplex_support=on.output @@ -0,0 +1,15 @@ + +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 diff --git a/tests/simplex/step-18.cc b/tests/simplex/step-18.cc new file mode 100644 index 0000000000..a146186b15 --- /dev/null +++ b/tests/simplex/step-18.cc @@ -0,0 +1,1053 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +// simplex +#include + +#include +#include +#include + +//#define HEX + +const unsigned int degree = 1; + +namespace Step18 +{ + using namespace dealii; + + template + struct PointHistory + { + SymmetricTensor<2, dim> old_stress; + }; + + template + 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 + inline SymmetricTensor<2, dim> + get_strain(const FEValues &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 + inline SymmetricTensor<2, dim> + get_strain(const std::vector> &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> &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> &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 + 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 triangulation; + + FESystem fe; + + DoFHandler dof_handler; + + AffineConstraints hanging_node_constraints; + + const Quadrature quadrature_formula; + +#ifdef HEX + MappingQGeneric mapping; +#else + MappingFE mapping; +#endif + + std::vector> quadrature_point_history; + +#ifdef DEAL_II_WITH_PETSC + PETScWrappers::MPI::SparseMatrix system_matrix; + PETScWrappers::MPI::Vector system_rhs; +#else + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector system_rhs; +#endif + + Vector 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 + class BodyForce : public Function + { + public: + BodyForce(); + + virtual void + vector_value(const Point &p, Vector &values) const override; + + virtual void + vector_value_list(const std::vector> &points, + std::vector> & value_list) const override; + }; + + + template + BodyForce::BodyForce() + : Function(dim) + {} + + + template + inline void + BodyForce::vector_value(const Point & /*p*/, + Vector &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 + void + BodyForce::vector_value_list( + const std::vector> &points, + std::vector> & 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::vector_value(points[p], value_list[p]); + } + + + + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues(const double present_time, + const double present_timestep); + + virtual void + vector_value(const Point &p, Vector &values) const override; + + virtual void + vector_value_list(const std::vector> &points, + std::vector> & value_list) const override; + + private: + const double velocity; + const double present_time; + const double present_timestep; + }; + + + template + IncrementalBoundaryValues::IncrementalBoundaryValues( + const double present_time, + const double present_timestep) + : Function(dim) + , velocity(.08) + , present_time(present_time) + , present_timestep(present_timestep) + {} + + + template + void + IncrementalBoundaryValues::vector_value(const Point & /*p*/, + Vector &values) const + { + Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim)); + + values = 0; + values(2) = -present_timestep * velocity; + } + + + + template + void + IncrementalBoundaryValues::vector_value_list( + const std::vector> &points, + std::vector> & 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::vector_value(points[p], value_list[p]); + } + + + template + const SymmetricTensor<4, dim> TopLevel::stress_strain_tensor = + get_stress_strain_tensor(/*lambda = */ 9.695e10, + /*mu = */ 7.617e10); + + +#ifdef HEX + template + TopLevel::TopLevel() + : triangulation() + , fe(FE_Q(degree), dim) + , dof_handler(triangulation) + , quadrature_formula(QGauss(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 + TopLevel::TopLevel() + : triangulation() + , fe(Simplex::FE_P(degree), dim) + , dof_handler(triangulation) + , quadrature_formula(Simplex::PGauss(fe.degree == 1 ? 4 : 10)) + , mapping(Simplex::FE_P(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 + TopLevel::~TopLevel() + { + dof_handler.clear(); + } + + + template + void + TopLevel::run() + { + do_initial_timestep(false); + + while (present_time < end_time) + do_timestep(std::abs(end_time - present_time - present_timestep) < 10e-5); + } + + + template + void + TopLevel::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 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 + void + TopLevel::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 + void + TopLevel::assemble_system() + { + system_rhs = 0; + system_matrix = 0; + + FEValues 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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + BodyForce body_force; + std::vector> body_force_values(n_q_points, + Vector(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 *local_quadrature_points_data = + reinterpret_cast *>(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 boundary_values; + VectorTools::interpolate_boundary_values(mapping, + dof_handler, + 0, + Functions::ZeroFunction(dim), + boundary_values); + VectorTools::interpolate_boundary_values( + mapping, + dof_handler, + 1, + IncrementalBoundaryValues(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 tmp(dof_handler.n_dofs()); +#endif + MatrixTools::apply_boundary_values( + boundary_values, system_matrix, tmp, system_rhs, false); + incremental_displacement = tmp; + } + + + template + void + TopLevel::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 + unsigned int + TopLevel::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 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> 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 + void + TopLevel::output_results() const + { + return; + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + + std::vector solution_names; + std::vector + 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, dim, dim>:: + DataVectorType::type_automatic, + solution_interpretation); + + + Vector 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 *>(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 partition_int( + triangulation.n_active_cells()); + GridTools::get_subdomain_association(triangulation, partition_int); + const Vector 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 + void + TopLevel::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 + void + TopLevel::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 + void + TopLevel::refine_initial_grid() + { + // First, let each process compute error indicators for the cells it owns: + Vector error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate( + dof_handler, + QGauss(fe.degree + 1), + std::map *>(), + 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 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 + void + TopLevel::move_mesh() + { + deallog << " Moving mesh..." << std::endl; + + std::vector 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 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 + void + TopLevel::setup_quadrature_point_history() + { + triangulation.clear_user_data(); + + { + std::vector> 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 + void + TopLevel::update_quadrature_point_history() + { + FEValues fe_values(mapping, + fe, + quadrature_formula, + update_values | update_gradients); + + std::vector>> displacement_increment_grads( + quadrature_formula.size(), std::vector>(dim)); + + for (auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history = + reinterpret_cast *>(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>(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; +} diff --git a/tests/simplex/step-18.mpirun=1.with_petsc=false.with_simplex_support=on.output b/tests/simplex/step-18.mpirun=1.with_petsc=false.with_simplex_support=on.output new file mode 100644 index 0000000000..0672e697e3 --- /dev/null +++ b/tests/simplex/step-18.mpirun=1.with_petsc=false.with_simplex_support=on.output @@ -0,0 +1,94 @@ + +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:: diff --git a/tests/simplex/step-18.mpirun=1.with_petsc=true.with_simplex_support=on.output b/tests/simplex/step-18.mpirun=1.with_petsc=true.with_simplex_support=on.output new file mode 100644 index 0000000000..8fb41aab10 --- /dev/null +++ b/tests/simplex/step-18.mpirun=1.with_petsc=true.with_simplex_support=on.output @@ -0,0 +1,94 @@ + +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:: -- 2.39.5