From 07af92a315cbe64494ecc6e38912d74c846fde10 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 9 Apr 2019 23:05:35 +0200 Subject: [PATCH] ScratchData + Sacado applied to energies --- tests/meshworker/scratch_data_06.cc | 177 ++++++++++++++++++ ...on.with_muparser=on.with_umfpack=on.output | 2 + 2 files changed, 179 insertions(+) create mode 100644 tests/meshworker/scratch_data_06.cc create mode 100644 tests/meshworker/scratch_data_06.with_trilinos=on.with_muparser=on.with_umfpack=on.output diff --git a/tests/meshworker/scratch_data_06.cc b/tests/meshworker/scratch_data_06.cc new file mode 100644 index 0000000000..9dd081d082 --- /dev/null +++ b/tests/meshworker/scratch_data_06.cc @@ -0,0 +1,177 @@ +/* --------------------------------------------------------------------- + * + * 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. + * + * --------------------------------------------------------------------- + */ + +// Compute Laplacian using mesh_loop, sacado, ScratchData, with +// second derivatives. + +#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 MeshWorker; + +typedef Sacado::Fad::DFad Sdouble; +typedef Sacado::Fad::DFad SSdouble; + +template +void +test() +{ + Triangulation tria; + FE_Q fe(1); + DoFHandler dh(tria); + + FunctionParser rhs_function("1"); + FunctionParser boundary_function("0"); + + AffineConstraints constraints; + + GridGenerator::hyper_cube(tria); + tria.refine_global(5); + tria.execute_coarsening_and_refinement(); + dh.distribute_dofs(fe); + + + DoFTools::make_hanging_node_constraints(dh, constraints); + VectorTools::interpolate_boundary_values(dh, + 0, + boundary_function, + constraints); + + constraints.close(); + + 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); + + UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<1, 1, 1>; + + ScratchData scratch(fe, quad, cell_flags); + CopyData copy(fe.dofs_per_cell); + + auto cell = dh.begin_active(); + auto endc = dh.end(); + + typedef decltype(cell) Iterator; + + FEValuesExtractors::Scalar scalar(0); + + auto cell_worker = [&rhs_function, &solution, &scalar](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(); + + SSdouble energy(0); + + s.extract_local_dof_values("solution", solution, energy); + const auto &u = s.get_values("solution", scalar, energy); + const auto &grad_u = s.get_gradients("solution", scalar, energy); + + c = 0; + + c.local_dof_indices[0] = s.get_local_dof_indices(); + + for (unsigned int q = 0; q < p.size(); ++q) + { + energy += + (.5 * grad_u[q] * grad_u[q] - rhs_function.value(p[q]) * u[q]) * + JxW[q]; + } + + // Copy from Sacado to normal vectors + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + c.vectors[0][i] = -energy.dx(i).val(); + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + c.matrices[0](i, j) = energy.dx(i).dx(j); + } + }; + + 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); + }; + + mesh_loop(cell, endc, cell_worker, copier, scratch, copy, assemble_own_cells); + + 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(1); + MultithreadInfo::set_thread_limit(1); // to make output deterministic + + test<2, 2>(); +} diff --git a/tests/meshworker/scratch_data_06.with_trilinos=on.with_muparser=on.with_umfpack=on.output b/tests/meshworker/scratch_data_06.with_trilinos=on.with_muparser=on.with_umfpack=on.output new file mode 100644 index 0000000000..b5ad838ee8 --- /dev/null +++ b/tests/meshworker/scratch_data_06.with_trilinos=on.with_muparser=on.with_umfpack=on.output @@ -0,0 +1,2 @@ + +DEAL::Linfty norm of solution: 0.0737281 -- 2.39.5