#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/q_collection.h>
DEAL_II_NAMESPACE_OPEN
const Quadrature<dim - 1> & quadrature,
const UpdateFlags update_flags);
+ /**
+ * Construct the FEInterfaceValues object with different FiniteElements
+ * assigned to different faces.
+ */
+ FEInterfaceValues(
+ const hp::MappingCollection<dim, spacedim> &mapping_collection,
+ const hp::FECollection<dim, spacedim> & fe_collection,
+ const hp::QCollection<dim - 1> & quadrature_collection,
+ const UpdateFlags update_flags);
+
+ /**
+ * Same as above, but using a Q1 mapping.
+ */
+ FEInterfaceValues(const hp::FECollection<dim, spacedim> &fe_collection,
+ const hp::QCollection<dim - 1> &quadrature_collection,
+ const UpdateFlags update_flags);
+
/**
* Re-initialize this object to be used on a new interface given by two faces
* of two neighboring cells. The `cell` and `cell_neighbor` cells will be
/**
* The FEFaceValues object for the current cell.
*/
- FEFaceValues<dim, spacedim> internal_fe_face_values;
+ std::unique_ptr<FEFaceValues<dim, spacedim>> internal_fe_face_values;
/**
* The FEFaceValues object for the current cell if the cell is refined.
*/
- FESubfaceValues<dim, spacedim> internal_fe_subface_values;
+ std::unique_ptr<FESubfaceValues<dim, spacedim>> internal_fe_subface_values;
/**
* The FEFaceValues object for the neighboring cell.
*/
- FEFaceValues<dim, spacedim> internal_fe_face_values_neighbor;
+ std::unique_ptr<FEFaceValues<dim, spacedim>> internal_fe_face_values_neighbor;
/**
* The FEFaceValues object for the neighboring cell if the cell is refined.
*/
- FESubfaceValues<dim, spacedim> internal_fe_subface_values_neighbor;
+ std::unique_ptr<FESubfaceValues<dim, spacedim>>
+ internal_fe_subface_values_neighbor;
/**
* Pointer to internal_fe_face_values or internal_fe_subface_values,
*/
FEFaceValuesBase<dim, spacedim> *fe_face_values_neighbor;
- /* Make the view classes friends of this class, since they access internal
- * data.
+ /**
+ * An hp::FEValues object for the FEFaceValues on the
+ * present cell.
+ */
+ std::unique_ptr<hp::FEFaceValues<dim, spacedim>> internal_hp_fe_face_values;
+
+ /**
+ * An hp::FEValues object for the FESubfaceValues on the
+ * present cell.
+ */
+ std::unique_ptr<hp::FESubfaceValues<dim, spacedim>>
+ internal_hp_fe_subface_values;
+
+ /**
+ * An hp::FEValues object for the FEFaceValues on the
+ * neighbor of the present cell.
+ */
+ std::unique_ptr<hp::FEFaceValues<dim, spacedim>>
+ internal_hp_fe_face_values_neighbor;
+
+ /**
+ * An hp::FEValues object for the FESubfaceValues on the
+ * neighboring cell.
+ */
+ std::unique_ptr<hp::FESubfaceValues<dim, spacedim>>
+ internal_hp_fe_subface_values_neighbor;
+
+
+ /*
+ * Make the view classes friends of this class, since they
+ * access internal data.
*/
template <int, int>
friend class FEInterfaceViews::Scalar;
const Quadrature<dim - 1> & quadrature,
const UpdateFlags update_flags)
: n_quadrature_points(quadrature.size())
- , internal_fe_face_values(mapping, fe, quadrature, update_flags)
- , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
- , internal_fe_face_values_neighbor(mapping, fe, quadrature, update_flags)
- , internal_fe_subface_values_neighbor(mapping, fe, quadrature, update_flags)
+ , internal_fe_face_values(
+ std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature,
+ update_flags))
+ , internal_fe_subface_values(
+ std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature,
+ update_flags))
+ , internal_fe_face_values_neighbor(
+ std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature,
+ update_flags))
+ , internal_fe_subface_values_neighbor(
+ std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature,
+ update_flags))
, fe_face_values(nullptr)
, fe_face_values_neighbor(nullptr)
{}
const hp::QCollection<dim - 1> & quadrature,
const UpdateFlags update_flags)
: n_quadrature_points(quadrature.max_n_quadrature_points())
- , internal_fe_face_values(mapping, fe, quadrature, update_flags)
- , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
- , internal_fe_face_values_neighbor(mapping, fe, quadrature[0], update_flags)
- , internal_fe_subface_values_neighbor(mapping,
- fe,
- quadrature[0],
- update_flags)
+
+ , internal_fe_face_values(
+ std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature,
+ update_flags))
+ , internal_fe_subface_values(
+ std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature,
+ update_flags))
+ , internal_fe_face_values_neighbor(
+ std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature[0],
+ update_flags))
+ , internal_fe_subface_values_neighbor(
+ std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+ fe,
+ quadrature[0],
+ update_flags))
, fe_face_values(nullptr)
, fe_face_values_neighbor(nullptr)
{}
+template <int dim, int spacedim>
+FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
+ const hp::MappingCollection<dim, spacedim> &mapping_collection,
+ const hp::FECollection<dim, spacedim> & fe_collection,
+ const hp::QCollection<dim - 1> & quadrature_collection,
+ const UpdateFlags update_flags)
+ : n_quadrature_points(quadrature_collection.max_n_quadrature_points())
+ , fe_face_values(nullptr)
+ , fe_face_values_neighbor(nullptr)
+ , internal_hp_fe_face_values(
+ std::make_unique<hp::FEFaceValues<dim, spacedim>>(mapping_collection,
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+ , internal_hp_fe_subface_values(
+ std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+ mapping_collection,
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+ , internal_hp_fe_face_values_neighbor(
+ std::make_unique<hp::FEFaceValues<dim, spacedim>>(mapping_collection,
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+ , internal_hp_fe_subface_values_neighbor(
+ std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+ mapping_collection,
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+{
+ AssertDimension(dim, spacedim);
+}
+
+
+
+template <int dim, int spacedim>
+FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
+ const hp::FECollection<dim, spacedim> &fe_collection,
+ const hp::QCollection<dim - 1> & quadrature_collection,
+ const UpdateFlags update_flags)
+ : n_quadrature_points(quadrature_collection.max_n_quadrature_points())
+ , fe_face_values(nullptr)
+ , fe_face_values_neighbor(nullptr)
+ , internal_hp_fe_face_values(
+ std::make_unique<hp::FEFaceValues<dim, spacedim>>(
+ fe_collection.get_reference_cell_default_linear_mapping(),
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+ , internal_hp_fe_subface_values(
+ std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+ fe_collection.get_reference_cell_default_linear_mapping(),
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+ , internal_hp_fe_face_values_neighbor(
+ std::make_unique<hp::FEFaceValues<dim, spacedim>>(
+ fe_collection.get_reference_cell_default_linear_mapping(),
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+ , internal_hp_fe_subface_values_neighbor(
+ std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+ fe_collection.get_reference_cell_default_linear_mapping(),
+ fe_collection,
+ quadrature_collection,
+ update_flags))
+{
+ AssertDimension(dim, spacedim);
+}
+
+
+
template <int dim, int spacedim>
FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim - 1> & quadrature,
const UpdateFlags update_flags)
: n_quadrature_points(quadrature.size())
- , internal_fe_face_values(
+ , internal_fe_face_values(std::make_unique<FEFaceValues<dim, spacedim>>(
fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
fe,
quadrature,
- update_flags)
- , internal_fe_subface_values(
+ update_flags))
+ , internal_fe_subface_values(std::make_unique<FESubfaceValues<dim, spacedim>>(
fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
fe,
quadrature,
- update_flags)
+ update_flags))
, internal_fe_face_values_neighbor(
- fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
- fe,
- quadrature,
- update_flags)
+ std::make_unique<FEFaceValues<dim, spacedim>>(
+ fe.reference_cell()
+ .template get_default_linear_mapping<dim, spacedim>(),
+ fe,
+ quadrature,
+ update_flags))
, internal_fe_subface_values_neighbor(
- fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
- fe,
- quadrature,
- update_flags)
+ std::make_unique<FESubfaceValues<dim, spacedim>>(
+ fe.reference_cell()
+ .template get_default_linear_mapping<dim, spacedim>(),
+ fe,
+ quadrature,
+ update_flags))
, fe_face_values(nullptr)
, fe_face_values_neighbor(nullptr)
{}
const unsigned int face_no_neighbor,
const unsigned int sub_face_no_neighbor)
{
- if (sub_face_no == numbers::invalid_unsigned_int)
- {
- internal_fe_face_values.reinit(cell, face_no);
- fe_face_values = &internal_fe_face_values;
- }
- else
- {
- internal_fe_subface_values.reinit(cell, face_no, sub_face_no);
- fe_face_values = &internal_fe_subface_values;
- }
- if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
+ Assert(internal_fe_face_values || internal_hp_fe_face_values,
+ ExcNotInitialized());
+
+ if (internal_fe_face_values)
{
- internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor);
- fe_face_values_neighbor = &internal_fe_face_values_neighbor;
+ if (sub_face_no == numbers::invalid_unsigned_int)
+ {
+ internal_fe_face_values->reinit(cell, face_no);
+ fe_face_values = internal_fe_face_values.get();
+ }
+ else
+ {
+ internal_fe_subface_values->reinit(cell, face_no, sub_face_no);
+ fe_face_values = internal_fe_subface_values.get();
+ }
+ if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
+ {
+ internal_fe_face_values_neighbor->reinit(cell_neighbor,
+ face_no_neighbor);
+ fe_face_values_neighbor = internal_fe_face_values_neighbor.get();
+ }
+ else
+ {
+ internal_fe_subface_values_neighbor->reinit(cell_neighbor,
+ face_no_neighbor,
+ sub_face_no_neighbor);
+ fe_face_values_neighbor = internal_fe_subface_values_neighbor.get();
+ }
+
+ AssertDimension(fe_face_values->n_quadrature_points,
+ fe_face_values_neighbor->n_quadrature_points);
+
+ const_cast<unsigned int &>(this->n_quadrature_points) =
+ fe_face_values->n_quadrature_points;
}
- else
+ else if (internal_hp_fe_face_values)
{
- internal_fe_subface_values_neighbor.reinit(cell_neighbor,
- face_no_neighbor,
- sub_face_no_neighbor);
- fe_face_values_neighbor = &internal_fe_subface_values_neighbor;
- }
+ const unsigned int dominated_fe_index =
+ internal_hp_fe_face_values->get_fe_collection().find_dominated_fe(
+ {cell->active_fe_index(), cell_neighbor->active_fe_index()});
+
+ // Same as if above, but when hp is enabled.
+ if (sub_face_no == numbers::invalid_unsigned_int)
+ {
+ internal_hp_fe_face_values->reinit(cell,
+ face_no,
+ dominated_fe_index,
+ dominated_fe_index);
+ fe_face_values = &const_cast<FEFaceValues<dim, spacedim> &>(
+ internal_hp_fe_face_values->get_present_fe_values());
+ }
+ else
+ {
+ internal_hp_fe_subface_values->reinit(
+ cell, face_no, sub_face_no, dominated_fe_index, dominated_fe_index);
+
+ fe_face_values = &const_cast<FESubfaceValues<dim, spacedim> &>(
+ internal_hp_fe_subface_values->get_present_fe_values());
+ }
+ if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
+ {
+ internal_hp_fe_face_values_neighbor->reinit(cell_neighbor,
+ face_no_neighbor,
+ dominated_fe_index,
+ dominated_fe_index);
+
+ fe_face_values_neighbor = &const_cast<FEFaceValues<dim, spacedim> &>(
+ internal_hp_fe_face_values_neighbor->get_present_fe_values());
+ }
+ else
+ {
+ internal_hp_fe_subface_values_neighbor->reinit(cell_neighbor,
+ face_no_neighbor,
+ sub_face_no_neighbor,
+ dominated_fe_index,
+ dominated_fe_index);
+
+ fe_face_values_neighbor =
+ &const_cast<FESubfaceValues<dim, spacedim> &>(
+ internal_hp_fe_subface_values_neighbor->get_present_fe_values());
+ }
- AssertDimension(fe_face_values->n_quadrature_points,
- fe_face_values_neighbor->n_quadrature_points);
+ AssertDimension(fe_face_values->n_quadrature_points,
+ fe_face_values_neighbor->n_quadrature_points);
- const_cast<unsigned int &>(this->n_quadrature_points) =
- fe_face_values->n_quadrature_points;
+ const_cast<unsigned int &>(this->n_quadrature_points) =
+ fe_face_values->n_quadrature_points;
+ }
// Set up dof mapping and remove duplicates (for continuous elements).
{
FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
const unsigned int face_no)
{
- internal_fe_face_values.reinit(cell, face_no);
- fe_face_values = &internal_fe_face_values;
- fe_face_values_neighbor = nullptr;
+ Assert(internal_fe_face_values || internal_hp_fe_face_values,
+ ExcNotInitialized());
- interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
- cell->get_active_or_mg_dof_indices(interface_dof_indices);
+ if (internal_fe_face_values)
+ {
+ internal_fe_face_values->reinit(cell, face_no);
+ fe_face_values = internal_fe_face_values.get();
+ fe_face_values_neighbor = nullptr;
- dofmap.resize(interface_dof_indices.size());
+ interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
+ cell->get_active_or_mg_dof_indices(interface_dof_indices);
+ }
+ else if (internal_hp_fe_face_values)
+ {
+ internal_hp_fe_face_values->reinit(cell, face_no);
+ fe_face_values = &const_cast<FEFaceValues<dim> &>(
+ internal_hp_fe_face_values->get_present_fe_values());
+ fe_face_values_neighbor = nullptr;
+ interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
+ cell->get_active_or_mg_dof_indices(interface_dof_indices);
+ }
+
+ dofmap.resize(interface_dof_indices.size());
for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
{
dofmap[i] = {{i, numbers::invalid_unsigned_int}};
const Mapping<dim, spacedim> &
FEInterfaceValues<dim, spacedim>::get_mapping() const
{
- return internal_fe_face_values.get_mapping();
+ return internal_fe_face_values->get_mapping();
}
const FiniteElement<dim, spacedim> &
FEInterfaceValues<dim, spacedim>::get_fe() const
{
- return internal_fe_face_values.get_fe();
+ return internal_fe_face_values->get_fe();
}
const Quadrature<dim - 1> &
FEInterfaceValues<dim, spacedim>::get_quadrature() const
{
- return internal_fe_face_values.get_quadrature();
+ return internal_fe_face_values->get_quadrature();
}
UpdateFlags
FEInterfaceValues<dim, spacedim>::get_update_flags() const
{
- return internal_fe_face_values.get_update_flags();
+ return internal_fe_face_values->get_update_flags();
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test basic properties of FEInterfaceValues in the hp case.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+inspect_fiv(FEInterfaceValues<dim> &fiv)
+{
+ deallog << "at_boundary(): " << fiv.at_boundary() << "\n"
+ << "n_current_interface_dofs(): " << fiv.n_current_interface_dofs()
+ << "\n";
+
+ std::vector<types::global_dof_index> indices =
+ fiv.get_interface_dof_indices();
+ Assert(indices.size() == fiv.n_current_interface_dofs(), ExcInternalError());
+
+ deallog << "interface_dof_indices: ";
+ for (auto i : indices)
+ deallog << i << ' ';
+ deallog << "\n";
+
+
+ unsigned int idx = 0;
+ for (auto v : indices)
+ {
+ deallog << " index " << idx << " global_dof_index:" << v << ":\n";
+
+ const auto pair = fiv.interface_dof_to_dof_indices(idx);
+ deallog << " dof indices: " << static_cast<int>(pair[0]) << " | "
+ << static_cast<int>(pair[1]) << "\n";
+
+ ++idx;
+ }
+
+ deallog << std::endl;
+}
+
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void
+make_2_cells<2>(Triangulation<2> &tria)
+{
+ const unsigned int dim = 2;
+ std::vector<unsigned int> repetitions = {2, 1};
+ Point<dim> p1;
+ Point<dim> p2(2.0, 1.0);
+
+ GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void
+make_2_cells<3>(Triangulation<3> &tria)
+{
+ const unsigned int dim = 3;
+ std::vector<unsigned int> repetitions = {2, 1, 1};
+ Point<dim> p1;
+ Point<dim> p2(2.0, 1.0, 1.0);
+
+ GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(const unsigned int p)
+{
+ Triangulation<dim> tria;
+ make_2_cells(tria);
+
+ DoFHandler<dim> dofh(tria);
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim - 1> q_collection;
+
+ fe_collection.push_back(FE_DGQ<dim>(p));
+ fe_collection.push_back(FE_DGQ<dim>(p + 1));
+
+ q_collection.push_back(QGauss<dim - 1>(p));
+ q_collection.push_back(QGauss<dim - 1>(p + 1));
+
+ // Set different finite elements spaces on the two cells.
+ unsigned int fe_index = 0;
+ for (const auto &cell : dofh.active_cell_iterators())
+ {
+ cell->set_active_fe_index(fe_index);
+ ++fe_index;
+ }
+
+ dofh.distribute_dofs(fe_collection);
+
+ UpdateFlags update_flags = update_values | update_gradients |
+ update_quadrature_points | update_JxW_values;
+
+ FEInterfaceValues<dim> fiv(fe_collection, q_collection, update_flags);
+
+
+ auto cell = dofh.begin();
+
+ deallog << "** interface between cell 0 and 1 **\n";
+
+ for (const unsigned int f : GeometryInfo<dim>::face_indices())
+ if (!cell->at_boundary(f))
+ {
+ fiv.reinit(cell,
+ f,
+ numbers::invalid_unsigned_int,
+ cell->neighbor(f),
+ cell->neighbor_of_neighbor(f),
+ numbers::invalid_unsigned_int);
+
+ Assert(fiv.get_fe_face_values(0).get_cell() == cell,
+ ExcInternalError());
+ Assert(fiv.get_fe_face_values(1).get_cell() == cell->neighbor(f),
+ ExcInternalError());
+ Assert(fiv.n_current_interface_dofs() ==
+ fe_collection[0].n_dofs_per_cell() +
+ fe_collection[1].n_dofs_per_cell(),
+ ExcInternalError());
+ Assert(!fiv.at_boundary(), ExcInternalError());
+
+ auto mycell = cell;
+ for (unsigned int c = 0; c < 2; ++c)
+ {
+ std::vector<types::global_dof_index> indices(
+ fe_collection[c].n_dofs_per_cell());
+ mycell->get_dof_indices(indices);
+ deallog << "cell " << c << ": ";
+ for (auto i : indices)
+ deallog << i << ' ';
+ deallog << "\n";
+ ++mycell;
+ }
+
+ inspect_fiv(fiv);
+ }
+
+ deallog << "** boundary interface on cell 0 **\n" << std::endl;
+
+ {
+ fiv.reinit(cell, 0);
+ Assert(fiv.get_fe_face_values(0).get_cell() == cell, ExcInternalError());
+ Assert(fiv.n_current_interface_dofs() == fe_collection[0].n_dofs_per_cell(),
+ ExcInternalError());
+ Assert(fiv.at_boundary(), ExcInternalError());
+ inspect_fiv(fiv);
+ }
+}
+
+
+
+int
+main()
+{
+ initlog();
+ for (unsigned int p : {1, 2})
+ test<2>(p);
+ for (unsigned int p : {1})
+ test<3>(p);
+}
--- /dev/null
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3
+cell 1: 4 5 6 7 8 9 10 11 12
+at_boundary(): 0
+n_current_interface_dofs(): 13
+interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12
+ index 0 global_dof_index:0:
+ dof indices: 0 | -1
+ index 1 global_dof_index:1:
+ dof indices: 1 | -1
+ index 2 global_dof_index:2:
+ dof indices: 2 | -1
+ index 3 global_dof_index:3:
+ dof indices: 3 | -1
+ index 4 global_dof_index:4:
+ dof indices: -1 | 0
+ index 5 global_dof_index:5:
+ dof indices: -1 | 1
+ index 6 global_dof_index:6:
+ dof indices: -1 | 2
+ index 7 global_dof_index:7:
+ dof indices: -1 | 3
+ index 8 global_dof_index:8:
+ dof indices: -1 | 4
+ index 9 global_dof_index:9:
+ dof indices: -1 | 5
+ index 10 global_dof_index:10:
+ dof indices: -1 | 6
+ index 11 global_dof_index:11:
+ dof indices: -1 | 7
+ index 12 global_dof_index:12:
+ dof indices: -1 | 8
+
+DEAL::** boundary interface on cell 0 **
+
+DEAL::at_boundary(): 1
+n_current_interface_dofs(): 4
+interface_dof_indices: 0 1 2 3
+ index 0 global_dof_index:0:
+ dof indices: 0 | -1
+ index 1 global_dof_index:1:
+ dof indices: 1 | -1
+ index 2 global_dof_index:2:
+ dof indices: 2 | -1
+ index 3 global_dof_index:3:
+ dof indices: 3 | -1
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3 4 5 6 7 8
+cell 1: 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
+at_boundary(): 0
+n_current_interface_dofs(): 25
+interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
+ index 0 global_dof_index:0:
+ dof indices: 0 | -1
+ index 1 global_dof_index:1:
+ dof indices: 1 | -1
+ index 2 global_dof_index:2:
+ dof indices: 2 | -1
+ index 3 global_dof_index:3:
+ dof indices: 3 | -1
+ index 4 global_dof_index:4:
+ dof indices: 4 | -1
+ index 5 global_dof_index:5:
+ dof indices: 5 | -1
+ index 6 global_dof_index:6:
+ dof indices: 6 | -1
+ index 7 global_dof_index:7:
+ dof indices: 7 | -1
+ index 8 global_dof_index:8:
+ dof indices: 8 | -1
+ index 9 global_dof_index:9:
+ dof indices: -1 | 0
+ index 10 global_dof_index:10:
+ dof indices: -1 | 1
+ index 11 global_dof_index:11:
+ dof indices: -1 | 2
+ index 12 global_dof_index:12:
+ dof indices: -1 | 3
+ index 13 global_dof_index:13:
+ dof indices: -1 | 4
+ index 14 global_dof_index:14:
+ dof indices: -1 | 5
+ index 15 global_dof_index:15:
+ dof indices: -1 | 6
+ index 16 global_dof_index:16:
+ dof indices: -1 | 7
+ index 17 global_dof_index:17:
+ dof indices: -1 | 8
+ index 18 global_dof_index:18:
+ dof indices: -1 | 9
+ index 19 global_dof_index:19:
+ dof indices: -1 | 10
+ index 20 global_dof_index:20:
+ dof indices: -1 | 11
+ index 21 global_dof_index:21:
+ dof indices: -1 | 12
+ index 22 global_dof_index:22:
+ dof indices: -1 | 13
+ index 23 global_dof_index:23:
+ dof indices: -1 | 14
+ index 24 global_dof_index:24:
+ dof indices: -1 | 15
+
+DEAL::** boundary interface on cell 0 **
+
+DEAL::at_boundary(): 1
+n_current_interface_dofs(): 9
+interface_dof_indices: 0 1 2 3 4 5 6 7 8
+ index 0 global_dof_index:0:
+ dof indices: 0 | -1
+ index 1 global_dof_index:1:
+ dof indices: 1 | -1
+ index 2 global_dof_index:2:
+ dof indices: 2 | -1
+ index 3 global_dof_index:3:
+ dof indices: 3 | -1
+ index 4 global_dof_index:4:
+ dof indices: 4 | -1
+ index 5 global_dof_index:5:
+ dof indices: 5 | -1
+ index 6 global_dof_index:6:
+ dof indices: 6 | -1
+ index 7 global_dof_index:7:
+ dof indices: 7 | -1
+ index 8 global_dof_index:8:
+ dof indices: 8 | -1
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3 4 5 6 7
+cell 1: 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
+at_boundary(): 0
+n_current_interface_dofs(): 35
+interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
+ index 0 global_dof_index:0:
+ dof indices: 0 | -1
+ index 1 global_dof_index:1:
+ dof indices: 1 | -1
+ index 2 global_dof_index:2:
+ dof indices: 2 | -1
+ index 3 global_dof_index:3:
+ dof indices: 3 | -1
+ index 4 global_dof_index:4:
+ dof indices: 4 | -1
+ index 5 global_dof_index:5:
+ dof indices: 5 | -1
+ index 6 global_dof_index:6:
+ dof indices: 6 | -1
+ index 7 global_dof_index:7:
+ dof indices: 7 | -1
+ index 8 global_dof_index:8:
+ dof indices: -1 | 0
+ index 9 global_dof_index:9:
+ dof indices: -1 | 1
+ index 10 global_dof_index:10:
+ dof indices: -1 | 2
+ index 11 global_dof_index:11:
+ dof indices: -1 | 3
+ index 12 global_dof_index:12:
+ dof indices: -1 | 4
+ index 13 global_dof_index:13:
+ dof indices: -1 | 5
+ index 14 global_dof_index:14:
+ dof indices: -1 | 6
+ index 15 global_dof_index:15:
+ dof indices: -1 | 7
+ index 16 global_dof_index:16:
+ dof indices: -1 | 8
+ index 17 global_dof_index:17:
+ dof indices: -1 | 9
+ index 18 global_dof_index:18:
+ dof indices: -1 | 10
+ index 19 global_dof_index:19:
+ dof indices: -1 | 11
+ index 20 global_dof_index:20:
+ dof indices: -1 | 12
+ index 21 global_dof_index:21:
+ dof indices: -1 | 13
+ index 22 global_dof_index:22:
+ dof indices: -1 | 14
+ index 23 global_dof_index:23:
+ dof indices: -1 | 15
+ index 24 global_dof_index:24:
+ dof indices: -1 | 16
+ index 25 global_dof_index:25:
+ dof indices: -1 | 17
+ index 26 global_dof_index:26:
+ dof indices: -1 | 18
+ index 27 global_dof_index:27:
+ dof indices: -1 | 19
+ index 28 global_dof_index:28:
+ dof indices: -1 | 20
+ index 29 global_dof_index:29:
+ dof indices: -1 | 21
+ index 30 global_dof_index:30:
+ dof indices: -1 | 22
+ index 31 global_dof_index:31:
+ dof indices: -1 | 23
+ index 32 global_dof_index:32:
+ dof indices: -1 | 24
+ index 33 global_dof_index:33:
+ dof indices: -1 | 25
+ index 34 global_dof_index:34:
+ dof indices: -1 | 26
+
+DEAL::** boundary interface on cell 0 **
+
+DEAL::at_boundary(): 1
+n_current_interface_dofs(): 8
+interface_dof_indices: 0 1 2 3 4 5 6 7
+ index 0 global_dof_index:0:
+ dof indices: 0 | -1
+ index 1 global_dof_index:1:
+ dof indices: 1 | -1
+ index 2 global_dof_index:2:
+ dof indices: 2 | -1
+ index 3 global_dof_index:3:
+ dof indices: 3 | -1
+ index 4 global_dof_index:4:
+ dof indices: 4 | -1
+ index 5 global_dof_index:5:
+ dof indices: 5 | -1
+ index 6 global_dof_index:6:
+ dof indices: 6 | -1
+ index 7 global_dof_index:7:
+ dof indices: 7 | -1
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate jump_in_shape_values(), average_of_shape_values(), shape_value() of
+// FEInterfaceValues on an adaptive mesh in the hp scenario.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void
+make_2_cells<2>(Triangulation<2> &tria)
+{
+ const unsigned int dim = 2;
+ std::vector<unsigned int> repetitions = {2, 1};
+ Point<dim> p1;
+ Point<dim> p2(2.0, 1.0);
+
+ GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void
+make_2_cells<3>(Triangulation<3> &tria)
+{
+ const unsigned int dim = 3;
+ std::vector<unsigned int> repetitions = {2, 1, 1};
+ Point<dim> p1;
+ Point<dim> p2(2.0, 1.0, 1.0);
+
+ GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(const unsigned int fe_degree0, const unsigned int fe_degree1 = 0)
+{
+ Triangulation<dim> tria;
+ make_2_cells(tria);
+
+ tria.begin()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+
+
+ UpdateFlags update_flags = update_values | update_gradients |
+ update_quadrature_points | update_JxW_values;
+ DoFHandler<dim> dofh(tria);
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim - 1> q_collection;
+ fe_collection.push_back(FE_DGQ<dim>(fe_degree0));
+ fe_collection.push_back(FE_DGQ<dim>(fe_degree1));
+
+ q_collection.push_back(QGauss<dim - 1>(fe_degree0 + 1));
+ // Set different finite elements spaces on the two cells.
+ unsigned int fe_index = 0;
+ for (const auto &cell : dofh.active_cell_iterators_on_level(0))
+ {
+ cell->set_active_fe_index(fe_index);
+ ++fe_index;
+ }
+
+ deallog << fe_collection[0].get_name() << "-" << fe_collection[1].get_name()
+ << std::endl;
+ dofh.distribute_dofs(fe_collection);
+
+ FEInterfaceValues<dim> fiv(fe_collection, q_collection, update_flags);
+
+ auto cell = dofh.begin(1);
+ ++cell;
+
+ for (const unsigned int f : GeometryInfo<dim>::face_indices())
+ if (!cell->at_boundary(f))
+ {
+ if (!cell->neighbor_is_coarser(f))
+ continue;
+
+ auto nn = cell->neighbor_of_coarser_neighbor(f);
+ fiv.reinit(cell,
+ f,
+ numbers::invalid_unsigned_int,
+ cell->neighbor(f),
+ nn.first,
+ nn.second);
+
+ const unsigned int n_dofs = fiv.n_current_interface_dofs();
+ Vector<double> cell_vector(n_dofs);
+
+ const auto &q_points = fiv.get_quadrature_points();
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ deallog << "qpoint " << qpoint << ": " << q_points[qpoint]
+ << std::endl;
+
+ for (unsigned int idx = 0; idx < n_dofs; ++idx)
+ {
+ const auto pair = fiv.interface_dof_to_dof_indices(idx);
+ deallog << " idx: " << idx
+ << " global: " << fiv.get_interface_dof_indices()[idx]
+ << " dof indices: " << static_cast<int>(pair[0]) << " | "
+ << static_cast<int>(pair[1]) << std::endl;
+ }
+
+ cell_vector = 0.0;
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ cell_vector(i) +=
+ fiv.shape_value(true, i, qpoint) * fiv.get_JxW_values()[qpoint];
+ deallog << "shape_value(true): " << cell_vector << std::endl;
+
+ cell_vector = 0.0;
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ cell_vector(i) +=
+ fiv.shape_value(false, i, qpoint) * fiv.get_JxW_values()[qpoint];
+ deallog << "shape_value(false): " << cell_vector << std::endl;
+
+ cell_vector = 0.0;
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ cell_vector(i) += fiv.jump_in_shape_values(i, qpoint) *
+ fiv.get_JxW_values()[qpoint];
+ deallog << "jump_in_shape_values(): " << cell_vector << std::endl;
+
+ cell_vector = 0.0;
+ for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ cell_vector(i) += fiv.average_of_shape_values(i, qpoint) *
+ fiv.get_JxW_values()[qpoint];
+ deallog << "average_of_shape_values(): " << cell_vector << std::endl;
+ }
+}
+
+
+
+int
+main()
+{
+ initlog();
+ test<2>(0);
+ test<2>(1, 2);
+ test<3>(0);
+ test<3>(1, 2);
+}
--- /dev/null
+
+DEAL::FE_DGQ<2>(0)-FE_DGQ<2>(0)
+DEAL::qpoint 0: 1.00000 0.250000
+DEAL:: idx: 0 global: 0 dof indices: -1 | 0
+DEAL:: idx: 1 global: 2 dof indices: 0 | -1
+DEAL::shape_value(true): 0.00000 0.500000
+DEAL::shape_value(false): 0.500000 0.00000
+DEAL::jump_in_shape_values(): -0.500000 0.500000
+DEAL::average_of_shape_values(): 0.250000 0.250000
+DEAL::FE_DGQ<2>(1)-FE_DGQ<2>(2)
+DEAL::qpoint 0: 1.00000 0.105662
+DEAL::qpoint 1: 1.00000 0.394338
+DEAL:: idx: 0 global: 0 dof indices: -1 | 0
+DEAL:: idx: 1 global: 1 dof indices: -1 | 1
+DEAL:: idx: 2 global: 2 dof indices: -1 | 2
+DEAL:: idx: 3 global: 3 dof indices: -1 | 3
+DEAL:: idx: 4 global: 8 dof indices: 0 | -1
+DEAL:: idx: 5 global: 9 dof indices: 1 | -1
+DEAL:: idx: 6 global: 10 dof indices: 2 | -1
+DEAL:: idx: 7 global: 11 dof indices: 3 | -1
+DEAL::shape_value(true): 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000
+DEAL::shape_value(false): 0.375000 0.00000 0.125000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::jump_in_shape_values(): -0.375000 0.00000 -0.125000 0.00000 0.00000 0.250000 0.00000 0.250000
+DEAL::average_of_shape_values(): 0.187500 0.00000 0.0625000 0.00000 0.00000 0.125000 0.00000 0.125000
+DEAL::FE_DGQ<3>(0)-FE_DGQ<3>(0)
+DEAL::qpoint 0: 1.00000 0.250000 0.250000
+DEAL:: idx: 0 global: 0 dof indices: -1 | 0
+DEAL:: idx: 1 global: 2 dof indices: 0 | -1
+DEAL::shape_value(true): 0.00000 0.250000
+DEAL::shape_value(false): 0.250000 0.00000
+DEAL::jump_in_shape_values(): -0.250000 0.250000
+DEAL::average_of_shape_values(): 0.125000 0.125000
+DEAL::FE_DGQ<3>(1)-FE_DGQ<3>(2)
+DEAL::qpoint 0: 1.00000 0.105662 0.105662
+DEAL::qpoint 1: 1.00000 0.394338 0.105662
+DEAL::qpoint 2: 1.00000 0.105662 0.394338
+DEAL::qpoint 3: 1.00000 0.394338 0.394338
+DEAL:: idx: 0 global: 0 dof indices: -1 | 0
+DEAL:: idx: 1 global: 1 dof indices: -1 | 1
+DEAL:: idx: 2 global: 2 dof indices: -1 | 2
+DEAL:: idx: 3 global: 3 dof indices: -1 | 3
+DEAL:: idx: 4 global: 4 dof indices: -1 | 4
+DEAL:: idx: 5 global: 5 dof indices: -1 | 5
+DEAL:: idx: 6 global: 6 dof indices: -1 | 6
+DEAL:: idx: 7 global: 7 dof indices: -1 | 7
+DEAL:: idx: 8 global: 16 dof indices: 0 | -1
+DEAL:: idx: 9 global: 17 dof indices: 1 | -1
+DEAL:: idx: 10 global: 18 dof indices: 2 | -1
+DEAL:: idx: 11 global: 19 dof indices: 3 | -1
+DEAL:: idx: 12 global: 20 dof indices: 4 | -1
+DEAL:: idx: 13 global: 21 dof indices: 5 | -1
+DEAL:: idx: 14 global: 22 dof indices: 6 | -1
+DEAL:: idx: 15 global: 23 dof indices: 7 | -1
+DEAL::shape_value(true): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000
+DEAL::shape_value(false): 0.140625 0.00000 0.0468750 0.00000 0.0468750 0.00000 0.0156250 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::jump_in_shape_values(): -0.140625 0.00000 -0.0468750 0.00000 -0.0468750 0.00000 -0.0156250 0.00000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000
+DEAL::average_of_shape_values(): 0.0703125 0.00000 0.0234375 0.00000 0.0234375 0.00000 0.00781250 0.00000 0.00000 0.0312500 0.00000 0.0312500 0.00000 0.0312500 0.00000 0.0312500
projects / dealii.git / commitdiff