* 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;
* 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;
* 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;
* 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;
* 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;
/**
* The FEFaceValues object for the current cell.
*/
- FEFaceValues<dim> internal_fe_face_values;
+ FEFaceValues<dim, spacedim> internal_fe_face_values;
/**
* The FEFaceValues object for the current cell if the cell is refined.
*/
- FESubfaceValues<dim> internal_fe_subface_values;
+ FESubfaceValues<dim, spacedim> internal_fe_subface_values;
/**
* The FEFaceValues object for the neighboring cell.
*/
- FEFaceValues<dim> internal_fe_face_values_neighbor;
+ FEFaceValues<dim, spacedim> internal_fe_face_values_neighbor;
/**
* The FEFaceValues object for the neighboring cell if the cell is refined.
*/
- FESubfaceValues<dim> internal_fe_subface_values_neighbor;
+ FESubfaceValues<dim, spacedim> internal_fe_subface_values_neighbor;
/**
* Pointer to internal_fe_face_values or internal_fe_subface_values,
* respectively as determined in reinit().
*/
- FEFaceValuesBase<dim> *fe_face_values;
+ FEFaceValuesBase<dim, spacedim> *fe_face_values;
/**
* Pointer to internal_fe_face_values_neighbor,
* internal_fe_subface_values_neighbor, or nullptr, respectively
* as determined in reinit().
*/
- FEFaceValuesBase<dim> *fe_face_values_neighbor;
+ FEFaceValuesBase<dim, spacedim> *fe_face_values_neighbor;
};
template <int dim, int spacedim>
-Tensor<1, dim>
+Tensor<1, spacedim>
FEInterfaceValues<dim, spacedim>::average_gradient(
const unsigned int interface_dof_index,
const unsigned int q_point,
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],
template <int dim, int spacedim>
-Tensor<2, dim>
+Tensor<2, spacedim>
FEInterfaceValues<dim, spacedim>::average_hessian(
const unsigned int interface_dof_index,
const unsigned int q_point,
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],
template <int dim, int spacedim>
-Tensor<1, dim>
+Tensor<1, spacedim>
FEInterfaceValues<dim, spacedim>::jump_gradient(
const unsigned int interface_dof_index,
const unsigned int q_point,
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],
template <int dim, int spacedim>
-Tensor<2, dim>
+Tensor<2, spacedim>
FEInterfaceValues<dim, spacedim>::jump_hessian(
const unsigned int interface_dof_index,
const unsigned int q_point,
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],
template <int dim, int spacedim>
-Tensor<3, dim>
+Tensor<3, spacedim>
FEInterfaceValues<dim, spacedim>::jump_3rd_derivative(
const unsigned int interface_dof_index,
const unsigned int q_point,
q_point,
component);
- Tensor<3, dim> value;
+ Tensor<3, spacedim> value;
if (dof_pair[0] != numbers::invalid_unsigned_int)
value +=
#include <deal.II/differentiation/ad.h>
+#include <deal.II/fe/fe_interface_values.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
* 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),
* 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
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<dim, spacedim> &
+ reinit(const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_face_no,
+ const typename DoFHandler<dim, spacedim>::active_cell_iterator
+ & cell_neighbor,
+ const unsigned int face_no_neighbor,
+ const unsigned int sub_face_no_neighbor);
+
/**
* Get the currently initialized FEValues.
*
* 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<dim, spacedim> &
reinit_neighbor(
* 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<dim, spacedim> &
reinit_neighbor(
* 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.
*/
std::unique_ptr<FESubfaceValues<dim, spacedim>> neighbor_fe_subface_values;
+ /**
+ * Interface values on facets.
+ */
+ std::unique_ptr<FEInterfaceValues<dim, spacedim>> interface_fe_values;
+
/**
* Dof indices on the current cell.
*/
* A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues
* object on this cell.
*/
- SmartPointer<FEValuesBase<dim, spacedim>> current_fe_values;
+ SmartPointer<const FEValuesBase<dim, spacedim>> current_fe_values;
/**
* A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues
* object on the neighbor cell.
*/
- SmartPointer<FEValuesBase<dim, spacedim>> current_neighbor_fe_values;
+ SmartPointer<const FEValuesBase<dim, spacedim>> current_neighbor_fe_values;
};
#ifndef DOXYGEN
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);
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;
*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;
fe_subface_values = std::make_unique<FESubfaceValues<dim, spacedim>>(
*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();
+ template <int dim, int spacedim>
+ const FEInterfaceValues<dim, spacedim> &
+ ScratchData<dim, spacedim>::reinit(
+ const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_face_no,
+ const typename DoFHandler<dim, spacedim>::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<FEInterfaceValues<dim, spacedim>>(
+ *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 <int dim, int spacedim>
const FEValues<dim, spacedim> &
ScratchData<dim, spacedim>::reinit_neighbor(
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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 <deal.II/base/function_lib.h>
+#include <deal.II/base/patterns.h>
+#include <deal.II/base/thread_management.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <unordered_map>
+
+#include "../tests.h"
+
+using namespace MeshWorker;
+
+template <int dim, int spacedim>
+void
+test()
+{
+ Triangulation<dim, spacedim> tria;
+ FE_DGQ<dim, spacedim> fe(1);
+ DoFHandler<dim, spacedim> dh(tria);
+
+ Functions::ConstantFunction<spacedim> rhs_function(1);
+ Functions::ConstantFunction<spacedim> boundary_function(0);
+
+ AffineConstraints<double> 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<double> matrix;
+ matrix.reinit(sparsity);
+
+ Vector<double> solution(dh.n_dofs());
+ Vector<double> rhs(dh.n_dofs());
+
+ QGauss<dim> quad(3);
+ QGauss<dim - 1> 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<dim, spacedim>;
+
+ auto cell = dh.begin_active();
+ auto endc = dh.end();
+
+ typedef decltype(cell) Iterator;
+
+ struct CopyDataFace
+ {
+ FullMatrix<double> cell_matrix;
+ std::vector<types::global_dof_index> joint_dof_indices;
+ };
+
+ struct CopyData
+ {
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+ std::vector<types::global_dof_index> local_dof_indices;
+ std::vector<CopyDataFace> face_data;
+
+ 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>();
+}
--- /dev/null
+
+DEAL::Linfty norm of solution 0.0736360
projects / dealii.git / commitdiff