Skip to content

Commit

Permalink
Merge pull request dealii#17244 from dominiktassilostill/Face_to_cell…
Browse files Browse the repository at this point in the history 
…_index_nodal_simplex

Add `face_to_cell_index_nodal` for simplex elements
  • Loading branch information
@kronbichler
kronbichler authored Jul 17, 2024
2 parents 3b075af + a682183 commit 008872e
Show file tree
Hide file tree
Showing 3 changed files with 235 additions and 0 deletions.
78 changes: 78 additions & 0 deletions include/deal.II/matrix_free/shape_info.templates.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -496,6 +496,84 @@ namespace internal

univariate_shape_data.nodal_at_cell_boundaries = true;

const ReferenceCell reference_cell = fe.reference_cell();
if (reference_cell.is_simplex())
{
if (dim == 2)
dofs_per_component_on_face = fe.degree + 1;
else
dofs_per_component_on_face =
(fe.degree + 1) * (fe.degree + 2) / 2;

face_to_cell_index_nodal.reinit(reference_cell.n_faces(),
dofs_per_component_on_face);


for (unsigned int face = 0; face < reference_cell.n_faces();
++face)
{
// first get info from reference cell, i.e. the linear case
unsigned int d = 0;
for (; d < dim; ++d)
face_to_cell_index_nodal[face][d] =
reference_cell.face_to_cell_vertices(face, d, 1);

// now fill the rest of the indices, start with the lines
if (fe.degree == 2)
for (; d < dofs_per_component_on_face; ++d)
face_to_cell_index_nodal[face][d] =
reference_cell.n_vertices() +
reference_cell.face_to_cell_lines(face, d - dim, 1);

// in the cubic case it is more complicated as more DoFs are
// on the lines
else if (fe.degree == 3)
{
for (unsigned int line = 0;
d < dofs_per_component_on_face - 1;
++line, d += 2)
{
const unsigned int face_to_cell_lines =
reference_cell.face_to_cell_lines(face, line, 1);
// check the direction of the line
// is it 0 -> 1 or 1 -> 0
// as DoFs on the line are ordered differently
if (reference_cell.line_to_cell_vertices(
face_to_cell_lines, 0) ==
reference_cell.face_to_cell_vertices(face,
line,
1))
{
face_to_cell_index_nodal[face][d] =
reference_cell.n_vertices() +
2 * face_to_cell_lines;
face_to_cell_index_nodal[face][d + 1] =
reference_cell.n_vertices() +
2 * face_to_cell_lines + 1;
}
else
{
face_to_cell_index_nodal[face][d + 1] =
reference_cell.n_vertices() +
2 * face_to_cell_lines;
face_to_cell_index_nodal[face][d] =
reference_cell.n_vertices() +
2 * face_to_cell_lines + 1;
}
}
// in 3D we also need the DoFs on the quads
if (dim == 3)
{
face_to_cell_index_nodal
[face][dofs_per_component_on_face - 1] =
reference_cell.n_vertices() +
2 * reference_cell.n_lines() + face;
}
}
else if (fe.degree > 3)
DEAL_II_NOT_IMPLEMENTED();
}
}
// TODO: set up face_to_cell_index_nodal, face_to_cell_index_hermite,
// face_orientations

Expand Down
102 changes: 102 additions & 0 deletions tests/matrix_free/face_to_cell_index_nodal_simplex.cc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2020 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------



// test the correctness of face_to_cell_index_nodal
// in internal::MatrixFreeFunctions::ShapeInfo
// for simplex elements

#include <deal.II/base/quadrature_lib.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_system.h>

#include <deal.II/grid/grid_generator.h>

#include <deal.II/matrix_free/shape_info.h>

#include <iostream>

#include "../tests.h"



template <int dim>
void
test(const FiniteElement<dim> &fe, const Quadrature<dim> &quad)
{
for (unsigned int i = 0; i < fe.n_base_elements(); ++i)
{
internal::MatrixFreeFunctions::ShapeInfo<double> shape_info;
shape_info.reinit(quad, fe, i);
deallog << "Testing: " << fe.get_name() << std::endl;

Triangulation<dim> tria;
GridGenerator::reference_cell(tria, fe.reference_cell());

FE_SimplexP<dim> fe_continous(fe.degree);

DoFHandler<dim> dof_handler(tria);
dof_handler.distribute_dofs(fe_continous);

std::vector<types::global_dof_index> dof_indices(
fe_continous.n_dofs_per_face());
for (const auto &cell : dof_handler.active_cell_iterators())
{
unsigned int face_counter = 0;
for (const auto &face : cell->face_iterators())
{
face->get_dof_indices(dof_indices);

unsigned int dof_counter = 0;
for (const auto i : dof_indices)
{
if (i == shape_info.face_to_cell_index_nodal[face_counter]
[dof_counter])
deallog << "Ok ";
else
deallog << "Not ok";
++dof_counter;
}
deallog << std::endl;
++face_counter;
}
}
}
}


int
main()
{
initlog();
for (unsigned int degree = 1; degree < 4; ++degree)
{
FE_SimplexP<2> fe(degree);
FE_SimplexDGP<2> fe_dg(degree);
QGaussSimplex<2> quad(1);
test(fe, quad);
test(fe_dg, quad);


FE_SimplexP<3> fe_3(degree);
FE_SimplexDGP<3> fe_dg_3(degree);
QGaussSimplex<3> quad_3(1);
test(fe_3, quad_3);
test(fe_dg_3, quad_3);
}
}
55 changes: 55 additions & 0 deletions tests/matrix_free/face_to_cell_index_nodal_simplex.output
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@

DEAL::Testing: FE_SimplexP<2>(1)
DEAL::Ok Ok
DEAL::Ok Ok
DEAL::Ok Ok
DEAL::Testing: FE_SimplexDGP<2>(1)
DEAL::Ok Ok
DEAL::Ok Ok
DEAL::Ok Ok
DEAL::Testing: FE_SimplexP<3>(1)
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Testing: FE_SimplexDGP<3>(1)
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Testing: FE_SimplexP<2>(2)
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Testing: FE_SimplexDGP<2>(2)
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Ok Ok Ok
DEAL::Testing: FE_SimplexP<3>(2)
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Testing: FE_SimplexDGP<3>(2)
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok
DEAL::Testing: FE_SimplexP<2>(3)
DEAL::Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok
DEAL::Testing: FE_SimplexDGP<2>(3)
DEAL::Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok
DEAL::Testing: FE_SimplexP<3>(3)
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok
DEAL::Testing: FE_SimplexDGP<3>(3)
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok
DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok

0 comments on commit 008872e

Please sign in to comment.