From: Luca Heltai Date: Tue, 9 Apr 2019 21:05:08 +0000 (+0200) Subject: Test ScratchData for Laplacian (SIPDG) X-Git-Tag: v9.1.0-rc1~210^2~7 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7d6357304090cccb8407443d95cf90be94012eed;p=dealii.git Test ScratchData for Laplacian (SIPDG) --- diff --git a/tests/meshworker/scratch_data_03.cc b/tests/meshworker/scratch_data_03.cc new file mode 100644 index 0000000000..c8dd6ebf9a --- /dev/null +++ b/tests/meshworker/scratch_data_03.cc @@ -0,0 +1,240 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2018 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 Laplacian using SIPG + mesh_loop + ScratchData + +#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 MeshWorker; + +template +void +test() +{ + Triangulation tria; + FE_DGQ fe(1); + DoFHandler dh(tria); + + FunctionParser rhs_function("1"); + FunctionParser boundary_function("0"); + + AffineConstraints constraints; + constraints.close(); + + GridGenerator::hyper_cube(tria); + tria.refine_global(5); + tria.execute_coarsening_and_refinement(); + dh.distribute_dofs(fe); + + + SparsityPattern sparsity; + + { + DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dh, dsp); + sparsity.copy_from(dsp); + } + + SparseMatrix matrix; + matrix.reinit(sparsity); + + Vector solution(dh.n_dofs()); + Vector rhs(dh.n_dofs()); + + QGauss quad(3); + QGauss face_quad(3); + + UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + UpdateFlags face_flags = update_values | update_gradients | + update_quadrature_points | + update_face_normal_vectors | update_JxW_values; + + // Stabilization for SIPG + double gamma = 100; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<1 + GeometryInfo::faces_per_cell, + 1, + 1 + GeometryInfo::faces_per_cell>; + + ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags); + CopyData copy(fe.dofs_per_cell); + + auto cell = dh.begin_active(); + auto endc = dh.end(); + + typedef decltype(cell) Iterator; + + auto cell_worker = + [&rhs_function](const Iterator &cell, ScratchData &s, CopyData &c) { + const auto &fev = s.reinit(cell); + const auto &JxW = s.get_JxW_values(); + const auto &p = s.get_quadrature_points(); + + c = 0; + + c.local_dof_indices[0] = s.get_local_dof_indices(); + + for (unsigned int q = 0; q < p.size(); ++q) + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + { + c.matrices[0](i, j) += + fev.shape_grad(i, q) * fev.shape_grad(j, q) * JxW[q]; + } + c.vectors[0](i) += + fev.shape_value(i, q) * rhs_function.value(p[q]) * JxW[q]; + } + }; + + auto boundary_worker = [gamma, &boundary_function](const Iterator & cell, + const unsigned int &f, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f); + const auto &JxW = s.get_JxW_values(); + const auto &p = s.get_quadrature_points(); + const auto &n = s.get_normal_vectors(); + + for (unsigned int q = 0; q < p.size(); ++q) + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + { + c.matrices[0](i, j) += + (-fev.shape_grad(i, q) * n[q] * fev.shape_value(j, q) + + -fev.shape_grad(j, q) * n[q] * fev.shape_value(i, q) + + gamma / cell->face(f)->diameter() * fev.shape_value(i, q) * + fev.shape_value(j, q)) * + JxW[q]; + } + c.vectors[0](i) += + ((gamma / cell->face(f)->diameter() * fev.shape_value(i, q) - + fev.shape_grad(i, q) * n[q]) * + boundary_function.value(p[q])) * + JxW[q]; + } + }; + + auto face_worker = [gamma](const Iterator & cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator & ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f, sf); + const auto &JxW = s.get_JxW_values(); + const auto &nfev = s.reinit_neighbor(ncell, nf, nsf); + + c.local_dof_indices[f + 1] = s.get_neighbor_dof_indices(); + + const auto &p = s.get_quadrature_points(); + const auto &n = s.get_normal_vectors(); + const auto &nn = s.get_neighbor_normal_vectors(); + + const double gh = gamma / cell->face(f)->diameter(); + + for (unsigned int q = 0; q < p.size(); ++q) + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + { + c.matrices[0](i, j) += + (-.5 * fev.shape_grad(i, q) * n[q] * fev.shape_value(j, q) + + -.5 * fev.shape_value(i, q) * n[q] * fev.shape_grad(j, q) + + gh * fev.shape_value(i, q) * fev.shape_value(j, q)) * + JxW[q]; + + c.matrices[f + 1](i, j) += + (-.5 * fev.shape_grad(i, q) * nn[q] * nfev.shape_value(j, q) + + -.5 * fev.shape_value(i, q) * n[q] * nfev.shape_grad(j, q) - + gh * fev.shape_value(i, q) * nfev.shape_value(j, q)) * + JxW[q]; + } + }; + + auto copier = [&constraints, &matrix, &rhs](const CopyData &c) { + constraints.distribute_local_to_global( + c.matrices[0], c.vectors[0], c.local_dof_indices[0], matrix, rhs); + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + constraints.distribute_local_to_global(c.matrices[1 + f], + c.local_dof_indices[0], + c.local_dof_indices[1 + f], + matrix); + }; + + mesh_loop(cell, + endc, + cell_worker, + copier, + scratch, + copy, + assemble_own_cells | assemble_boundary_faces | + assemble_own_interior_faces_both, + boundary_worker, + face_worker); + + SparseDirectUMFPACK inv; + inv.initialize(matrix); + + inv.vmult(solution, rhs); + constraints.distribute(solution); + + deallog << "Linfty norm of solution " << solution.linfty_norm() << std::endl; +} + + +int +main() +{ + initlog(); + test<2, 2>(); +} diff --git a/tests/meshworker/scratch_data_03.with_muparser=on.with_umfpack=on.output b/tests/meshworker/scratch_data_03.with_muparser=on.with_umfpack=on.output new file mode 100644 index 0000000000..1b2777af16 --- /dev/null +++ b/tests/meshworker/scratch_data_03.with_muparser=on.with_umfpack=on.output @@ -0,0 +1,2 @@ + +DEAL::Linfty norm of solution 0.0737281