From: Luca Heltai Date: Tue, 23 Jun 2020 20:24:58 +0000 (+0200) Subject: Added FEInterfaceValues to ScratchData. X-Git-Tag: v9.3.0-rc1~1382^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7471a504c573652e24bdb86cf34e5f8576d370dd;p=dealii.git Added FEInterfaceValues to ScratchData. --- diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index c2c20978fb..7412b0a608 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -340,7 +340,7 @@ public: * If this is a boundary face (at_boundary() returns true), then * $\average{\nabla u}=\nabla u_{\text{cell0}}$. */ - Tensor<1, dim> + Tensor<1, spacedim> average_gradient(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component = 0) const; @@ -355,7 +355,7 @@ public: * If this is a boundary face (at_boundary() returns true), then * $\average{\nabla^2 u}=\nabla^2 u_{\text{cell0}}$. */ - Tensor<2, dim> + Tensor<2, spacedim> average_hessian(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component = 0) const; @@ -369,7 +369,7 @@ public: * If this is a boundary face (at_boundary() returns true), then * $\jump{\nabla u}=\nabla u_{\text{cell0}}$. */ - Tensor<1, dim> + Tensor<1, spacedim> jump_gradient(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component = 0) const; @@ -384,7 +384,7 @@ public: * If this is a boundary face (at_boundary() returns true), then * $\jump{\nabla^2 u} = \nabla^2 u_{\text{cell0}}$. */ - Tensor<2, dim> + Tensor<2, spacedim> jump_hessian(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component = 0) const; @@ -398,7 +398,7 @@ public: * If this is a boundary face (at_boundary() returns true), then * $\jump{\nabla^3 u} = \nabla^3 u_{\text{cell0}}$. */ - Tensor<3, dim> + Tensor<3, spacedim> jump_3rd_derivative(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component = 0) const; @@ -423,35 +423,35 @@ private: /** * The FEFaceValues object for the current cell. */ - FEFaceValues internal_fe_face_values; + FEFaceValues internal_fe_face_values; /** * The FEFaceValues object for the current cell if the cell is refined. */ - FESubfaceValues internal_fe_subface_values; + FESubfaceValues internal_fe_subface_values; /** * The FEFaceValues object for the neighboring cell. */ - FEFaceValues internal_fe_face_values_neighbor; + FEFaceValues internal_fe_face_values_neighbor; /** * The FEFaceValues object for the neighboring cell if the cell is refined. */ - FESubfaceValues internal_fe_subface_values_neighbor; + FESubfaceValues internal_fe_subface_values_neighbor; /** * Pointer to internal_fe_face_values or internal_fe_subface_values, * respectively as determined in reinit(). */ - FEFaceValuesBase *fe_face_values; + FEFaceValuesBase *fe_face_values; /** * Pointer to internal_fe_face_values_neighbor, * internal_fe_subface_values_neighbor, or nullptr, respectively * as determined in reinit(). */ - FEFaceValuesBase *fe_face_values_neighbor; + FEFaceValuesBase *fe_face_values_neighbor; }; @@ -818,7 +818,7 @@ FEInterfaceValues::average( template -Tensor<1, dim> +Tensor<1, spacedim> FEInterfaceValues::average_gradient( const unsigned int interface_dof_index, const unsigned int q_point, @@ -831,7 +831,7 @@ FEInterfaceValues::average_gradient( q_point, component); - Tensor<1, dim> value; + Tensor<1, spacedim> value; if (dof_pair[0] != numbers::invalid_unsigned_int) value += 0.5 * get_fe_face_values(0).shape_grad_component(dof_pair[0], @@ -848,7 +848,7 @@ FEInterfaceValues::average_gradient( template -Tensor<2, dim> +Tensor<2, spacedim> FEInterfaceValues::average_hessian( const unsigned int interface_dof_index, const unsigned int q_point, @@ -861,7 +861,7 @@ FEInterfaceValues::average_hessian( q_point, component); - Tensor<2, dim> value; + Tensor<2, spacedim> value; if (dof_pair[0] != numbers::invalid_unsigned_int) value += 0.5 * get_fe_face_values(0).shape_hessian_component(dof_pair[0], @@ -878,7 +878,7 @@ FEInterfaceValues::average_hessian( template -Tensor<1, dim> +Tensor<1, spacedim> FEInterfaceValues::jump_gradient( const unsigned int interface_dof_index, const unsigned int q_point, @@ -891,7 +891,7 @@ FEInterfaceValues::jump_gradient( q_point, component); - Tensor<1, dim> value; + Tensor<1, spacedim> value; if (dof_pair[0] != numbers::invalid_unsigned_int) value += 1.0 * get_fe_face_values(0).shape_grad_component(dof_pair[0], @@ -908,7 +908,7 @@ FEInterfaceValues::jump_gradient( template -Tensor<2, dim> +Tensor<2, spacedim> FEInterfaceValues::jump_hessian( const unsigned int interface_dof_index, const unsigned int q_point, @@ -921,7 +921,7 @@ FEInterfaceValues::jump_hessian( q_point, component); - Tensor<2, dim> value; + Tensor<2, spacedim> value; if (dof_pair[0] != numbers::invalid_unsigned_int) value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0], @@ -937,7 +937,7 @@ FEInterfaceValues::jump_hessian( template -Tensor<3, dim> +Tensor<3, spacedim> FEInterfaceValues::jump_3rd_derivative( const unsigned int interface_dof_index, const unsigned int q_point, @@ -950,7 +950,7 @@ FEInterfaceValues::jump_3rd_derivative( q_point, component); - Tensor<3, dim> value; + Tensor<3, spacedim> value; if (dof_pair[0] != numbers::invalid_unsigned_int) value += diff --git a/include/deal.II/meshworker/scratch_data.h b/include/deal.II/meshworker/scratch_data.h index e83f97cb6f..3a18fee1b0 100644 --- a/include/deal.II/meshworker/scratch_data.h +++ b/include/deal.II/meshworker/scratch_data.h @@ -22,6 +22,7 @@ #include +#include #include #include @@ -54,9 +55,9 @@ namespace MeshWorker * function and with the MeshWorker::mesh_loop function(). * * The ScratchData class has three main goals: - * - to create FEValues, FEFaceValues, and FESubfaceValues for the current - * cell and for its neighbor cell *on demand* (only if they are - * necessary for the algorithm provided by the user), and to provide a + * - to create FEValues, FEFaceValues, FESubfaceValues, and FEInterfaceValues + * for the current cell and for its neighbor cell *on demand* (only if they + * are necessary for the algorithm provided by the user), and to provide a * uniform interface to access the FEValues objects when assembling cell, * face, or subface contributions * - to store arbitrary data types (or references to arbitrary data types), @@ -70,10 +71,10 @@ namespace MeshWorker * values, gradients, etc., of already computed solution vectors. * * The methods in the section "Methods to work on current cell" - * initialize on demand internal FEValues, - * FEFaceValues, and FESubfaceValues objects on the current cell, allowing the - * use of this class as a single substitute for three different objects used - * to integrate and query finite element values on cells, faces, and subfaces. + * initialize on demand internal FEValues, FEFaceValues, FESubfaceValues, and + * FEInterfaceValues objects on the current cell, allowing the use of this + * class as a single substitute for four different objects used to integrate + * and query finite element values on cells, faces, and subfaces. * * Similarly, the methods in the section "Methods to work on neighbor cell" * initialize on demand @@ -372,6 +373,29 @@ namespace MeshWorker const unsigned int face_no, const unsigned int subface_no); + /** + * Initialize the internal FEInterfaceValues with the given arguments, and + * return a reference to it. + * + * After calling this function, get_local_dof_indices(), + * get_quadrature_points(), get_normal_vectors(), and get_JxW_values() will + * be forwarded to the local FEInterfaceValues object. The methods + * get_current_fe_values() will return the FEValuesBase associated to the + * current cell, while get_neighbor_fe_values() will be associated with the + * neighbor cell. The method get_local_dof_indices() will return the + * same result of FEInterfaceValues::get_interface_dof_to_dof_indices(), + * while the get_neighbor_dof_indices() will return the local dof indices + * of the neighbor cell. + */ + const FEInterfaceValues & + reinit(const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_face_no, + const typename DoFHandler::active_cell_iterator + & cell_neighbor, + const unsigned int face_no_neighbor, + const unsigned int sub_face_no_neighbor); + /** * Get the currently initialized FEValues. * @@ -421,8 +445,8 @@ namespace MeshWorker * Initialize the internal neighbor FEValues to use the given @p cell, and * return a reference to it. * - * After calling this function, get_current_fe_values() will return the - * same object of this method, as an FEValuesBase reference. + * After calling this function, get_current_neighbor_fe_values() will return + * the same object of this method, as an FEValuesBase reference. */ const FEValues & reinit_neighbor( @@ -432,8 +456,8 @@ namespace MeshWorker * Initialize the internal FEFaceValues to use the given @p face_no on the * given @p cell, and return a reference to it. * - * After calling this function, get_current_fe_values() will return the - * same object of this method, as an FEValuesBase reference. + * After calling this function, get_current_neighbor_fe_values() will return + * the same object of this method, as an FEValuesBase reference. */ const FEFaceValues & reinit_neighbor( @@ -444,8 +468,8 @@ namespace MeshWorker * Initialize the internal FESubfaceValues to use the given @p subface_no, * on @p face_no, on the given @p cell, and return a reference to it. * - * After calling this function, get_current_fe_values() will return the - * same object of this method, as an FEValuesBase reference. + * After calling this function, get_current_neighbor_fe_values() will return + * the same object of this method, as an FEValuesBase reference. * * If @p subface_no is numbers::invalid_unsigned_int, the reinit() function * that takes only the @p cell and the @p face_no is called. @@ -862,6 +886,11 @@ namespace MeshWorker */ std::unique_ptr> neighbor_fe_subface_values; + /** + * Interface values on facets. + */ + std::unique_ptr> interface_fe_values; + /** * Dof indices on the current cell. */ @@ -886,13 +915,13 @@ namespace MeshWorker * A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues * object on this cell. */ - SmartPointer> current_fe_values; + SmartPointer> current_fe_values; /** * A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues * object on the neighbor cell. */ - SmartPointer> current_neighbor_fe_values; + SmartPointer> current_neighbor_fe_values; }; #ifndef DOXYGEN @@ -936,7 +965,7 @@ namespace MeshWorker const VectorType & input_vector, const Number dummy) { - const unsigned int n_dofs = get_current_fe_values().get_fe().dofs_per_cell; + const unsigned int n_dofs = local_dof_indices.size(); const std::string name = get_unique_dofs_name(global_vector_name, n_dofs, dummy); diff --git a/source/meshworker/scratch_data.cc b/source/meshworker/scratch_data.cc index f26a4a095c..c2229cce32 100644 --- a/source/meshworker/scratch_data.cc +++ b/source/meshworker/scratch_data.cc @@ -136,6 +136,7 @@ namespace MeshWorker cell_update_flags); fe_values->reinit(cell); + local_dof_indices.resize(fe_values->dofs_per_cell); cell->get_dof_indices(local_dof_indices); current_fe_values = fe_values.get(); return *fe_values; @@ -154,6 +155,7 @@ namespace MeshWorker *mapping, *fe, face_quadrature, face_update_flags); fe_face_values->reinit(cell, face_no); + local_dof_indices.resize(fe->dofs_per_cell); cell->get_dof_indices(local_dof_indices); current_fe_values = fe_face_values.get(); return *fe_face_values; @@ -174,6 +176,7 @@ namespace MeshWorker fe_subface_values = std::make_unique>( *mapping, *fe, face_quadrature, face_update_flags); fe_subface_values->reinit(cell, face_no, subface_no); + local_dof_indices.resize(fe->dofs_per_cell); cell->get_dof_indices(local_dof_indices); current_fe_values = fe_subface_values.get(); @@ -185,6 +188,37 @@ namespace MeshWorker + template + const FEInterfaceValues & + ScratchData::reinit( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_face_no, + const typename DoFHandler::active_cell_iterator + & cell_neighbor, + const unsigned int face_no_neighbor, + const unsigned int sub_face_no_neighbor) + { + if (!interface_fe_values) + interface_fe_values = std::make_unique>( + *mapping, *fe, face_quadrature, face_update_flags); + interface_fe_values->reinit(cell, + face_no, + sub_face_no, + cell_neighbor, + face_no_neighbor, + sub_face_no_neighbor); + + current_fe_values = &interface_fe_values->get_fe_face_values(0); + current_neighbor_fe_values = &interface_fe_values->get_fe_face_values(1); + + cell_neighbor->get_dof_indices(neighbor_dof_indices); + local_dof_indices = interface_fe_values->get_interface_dof_indices(); + return *interface_fe_values; + } + + + template const FEValues & ScratchData::reinit_neighbor( diff --git a/tests/meshworker/scratch_data_08.cc b/tests/meshworker/scratch_data_08.cc new file mode 100644 index 0000000000..6c2c5260a7 --- /dev/null +++ b/tests/meshworker/scratch_data_08.cc @@ -0,0 +1,257 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2018 - 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 Laplacian using SIPG + mesh_loop + ScratchData + FEInterfaceData + +#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); + + Functions::ConstantFunction rhs_function(1); + Functions::ConstantFunction 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; + + auto cell = dh.begin_active(); + auto endc = dh.end(); + + typedef decltype(cell) Iterator; + + 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; + + 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); + } + }; + + ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags); + CopyData copy; + + 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.reinit(cell, s.get_local_dof_indices().size()); + + 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.cell_matrix(i, j) += + fev.shape_grad(i, q) * fev.shape_grad(j, q) * JxW[q]; + } + c.cell_rhs(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.cell_matrix(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.cell_rhs(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, ncell, nf, nsf); + const auto &JxW = s.get_JxW_values(); + + const auto &p = s.get_quadrature_points(); + const auto &n = s.get_normal_vectors(); + + + c.face_data.emplace_back(); + auto ©_data_face = c.face_data.back(); + auto &face_matrix = copy_data_face.cell_matrix; + copy_data_face.joint_dof_indices = fev.get_interface_dof_indices(); + const auto n_dofs = fev.n_current_interface_dofs(); + face_matrix.reinit(n_dofs, n_dofs); + + const double gh = gamma / cell->face(f)->diameter(); + + for (unsigned int q = 0; q < p.size(); ++q) + for (unsigned int i = 0; i < n_dofs; ++i) + for (unsigned int j = 0; j < n_dofs; ++j) + { + face_matrix(i, j) += + (-fev.jump_gradient(i, q) * n[q] * fev.average(j, q) - + fev.average(i, q) * fev.jump_gradient(j, q) * n[q] + + gh * fev.jump(i, q) * fev.jump(j, q)) * + JxW[q]; + } + }; + + auto copier = [&constraints, &matrix, &rhs](const CopyData &c) { + constraints.distribute_local_to_global( + c.cell_matrix, c.cell_rhs, c.local_dof_indices, matrix, rhs); + + for (auto &cdf : c.face_data) + constraints.distribute_local_to_global(cdf.cell_matrix, + cdf.joint_dof_indices, + matrix); + }; + + + + mesh_loop(cell, + endc, + cell_worker, + copier, + scratch, + copy, + assemble_own_cells | assemble_boundary_faces | + assemble_own_interior_faces_once, + boundary_worker, + face_worker); + + SparseDirectUMFPACK inv; + inv.initialize(matrix); + + inv.vmult(solution, rhs); + + deallog << "Linfty norm of solution " << solution.linfty_norm() << std::endl; +} + + +int +main() +{ + initlog(); + test<2, 2>(); +} diff --git a/tests/meshworker/scratch_data_08.output b/tests/meshworker/scratch_data_08.output new file mode 100644 index 0000000000..200ebcc8cd --- /dev/null +++ b/tests/meshworker/scratch_data_08.output @@ -0,0 +1,2 @@ + +DEAL::Linfty norm of solution 0.0736360