const unsigned int fe_index,
const unsigned int dof_handler_index = 0) const;
+ /**
+ * In the hp adaptive case, a subrange of internal faces as computed during
+ * loop() might contain internal faces with elements of different active
+ * fe indices. Use this function to compute what the subrange for a given pair
+ * of active fe indices is.
+ */
+ std::pair<unsigned int, unsigned int>
+ create_inner_face_subrange_hp_by_index(
+ const std::pair<unsigned int, unsigned int> &range,
+ const unsigned int fe_index_interior,
+ const unsigned int fe_index_exterior,
+ const unsigned int dof_handler_index = 0) const;
+
+ /**
+ * In the hp adaptive case, a subrange of boundary faces as computed during
+ * loop() might contain boundary faces with elements of different active
+ * fe indices. Use this function to compute what the subrange for a given
+ * active fe indices is.
+ */
+ std::pair<unsigned int, unsigned int>
+ create_boundary_face_subrange_hp_by_index(
+ const std::pair<unsigned int, unsigned int> &range,
+ const unsigned int fe_index,
+ const unsigned int dof_handler_index = 0) const;
+
+ /**
+ * In the hp adaptive case, return number of active_fe_indices.
+ */
+ unsigned int
+ n_active_fe_indices() const;
+
+ /**
+ * In the hp adaptive case, return the active_fe_index of a cell range.
+ */
+ unsigned int
+ get_cell_active_fe_index(
+ const std::pair<unsigned int, unsigned int> range) const;
+
+ /**
+ * In the hp adaptive case, return the active_fe_index of a face range.
+ */
+ unsigned int
+ get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
+ const bool is_interior_face = true) const;
+
//@}
/**
+template <int dim, typename Number, typename VectorizedArrayType>
+unsigned int
+MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
+{
+ return shape_info.size(2);
+}
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+unsigned int
+MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
+ const std::pair<unsigned int, unsigned int> range) const
+{
+ const auto &fe_indices = dof_info[0].cell_active_fe_index;
+
+ if (fe_indices.empty() == true)
+ return 0;
+
+ const auto index = fe_indices[range.first];
+
+ for (unsigned int i = range.first; i < range.second; ++i)
+ AssertDimension(index, fe_indices[i]);
+
+ return index;
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+unsigned int
+MatrixFree<dim, Number, VectorizedArrayType>::get_face_active_fe_index(
+ const std::pair<unsigned int, unsigned int> range,
+ const bool is_interior_face) const
+{
+ const auto &fe_indices = dof_info[0].cell_active_fe_index;
+
+ if (fe_indices.empty() == true)
+ return 0;
+
+ if (is_interior_face)
+ {
+ const unsigned int index =
+ fe_indices[face_info.faces[range.first].cells_interior[0] /
+ VectorizedArrayType::size()];
+
+ for (unsigned int i = range.first; i < range.second; ++i)
+ AssertDimension(index,
+ fe_indices[face_info.faces[i].cells_interior[0] /
+ VectorizedArrayType::size()]);
+
+ return index;
+ }
+ else
+ {
+ const unsigned int index =
+ fe_indices[face_info.faces[range.first].cells_exterior[0] /
+ VectorizedArrayType::size()];
+
+ for (unsigned int i = range.first; i < range.second; ++i)
+ AssertDimension(index,
+ fe_indices[face_info.faces[i].cells_exterior[0] /
+ VectorizedArrayType::size()]);
+
+ return index;
+ }
+}
+
+
+
template <int dim, typename Number, typename VectorizedArrayType>
inline unsigned int
MatrixFree<dim, Number, VectorizedArrayType>::n_components_filled(
cell(const std::pair<unsigned int, unsigned int> &cell_range) override
{
if (cell_function != nullptr && cell_range.second > cell_range.first)
- (container.*
- cell_function)(matrix_free, this->dst, this->src, cell_range);
+ for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
+ {
+ const auto cell_subrange =
+ matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
+
+ if (cell_subrange.second <= cell_subrange.first)
+ continue;
+
+ (container.*
+ cell_function)(matrix_free, this->dst, this->src, cell_subrange);
+ }
}
// Runs the assembler on interior faces. If no function is given, nothing
face(const std::pair<unsigned int, unsigned int> &face_range) override
{
if (face_function != nullptr && face_range.second > face_range.first)
- (container.*
- face_function)(matrix_free, this->dst, this->src, face_range);
+ for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
+ for (unsigned int j = 0; j < matrix_free.n_active_fe_indices(); ++j)
+ {
+ const auto face_subrange =
+ matrix_free.create_inner_face_subrange_hp_by_index(face_range,
+ i,
+ j);
+
+ if (face_subrange.second <= face_subrange.first)
+ continue;
+ (container.*
+ face_function)(matrix_free, this->dst, this->src, face_subrange);
+ }
}
// Runs the assembler on boundary faces. If no function is given, nothing
boundary(const std::pair<unsigned int, unsigned int> &face_range) override
{
if (boundary_function != nullptr && face_range.second > face_range.first)
- (container.*
- boundary_function)(matrix_free, this->dst, this->src, face_range);
+ for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
+ {
+ const auto face_subrange =
+ matrix_free.create_boundary_face_subrange_hp_by_index(face_range,
+ i);
+
+ if (face_subrange.second <= face_subrange.first)
+ continue;
+
+ (container.*boundary_function)(matrix_free,
+ this->dst,
+ this->src,
+ face_subrange);
+ }
}
// Starts the communication for the update ghost values operation. We
MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp_by_index(
const std::pair<unsigned int, unsigned int> &range,
const unsigned int fe_index,
- const unsigned int vector_component) const
+ const unsigned int dof_handler_index) const
{
- AssertIndexRange(fe_index, dof_info[vector_component].max_fe_index);
+ if (dof_info[dof_handler_index].max_fe_index == 0)
+ return range;
+
+ AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index);
const std::vector<unsigned int> &fe_indices =
- dof_info[vector_component].cell_active_fe_index;
- if (fe_indices.empty() == true)
+ dof_info[dof_handler_index].cell_active_fe_index;
+
+ if (fe_indices.empty() == true ||
+ dof_handlers[dof_handler_index]->get_fe_collection().size() == 1)
return range;
else
{
+namespace
+{
+ class FaceRangeCompartor
+ {
+ public:
+ FaceRangeCompartor(const std::vector<unsigned int> &fe_indices,
+ const bool include)
+ : fe_indices(fe_indices)
+ , include(include)
+ {}
+
+ template <int vectorization_width>
+ bool
+ operator()(const internal::MatrixFreeFunctions::FaceToCellTopology<
+ vectorization_width> &face,
+ const unsigned int & fe_index)
+ {
+ const unsigned int fe_index_face =
+ fe_indices[face.cells_interior[0] / vectorization_width];
+
+ return fe_index_face < fe_index ||
+ (include ? (fe_index_face == fe_index) : false);
+ }
+
+ template <int vectorization_width>
+ bool
+ operator()(const internal::MatrixFreeFunctions::FaceToCellTopology<
+ vectorization_width> & face,
+ const std::pair<unsigned int, unsigned int> &fe_index)
+ {
+ const std::pair<unsigned int, unsigned int> fe_index_face = {
+ fe_indices[face.cells_interior[0] / vectorization_width],
+ fe_indices[face.cells_exterior[0] / vectorization_width]};
+
+ return include ? (fe_index_face <= fe_index) : (fe_index_face < fe_index);
+ }
+
+ private:
+ const std::vector<unsigned int> &fe_indices;
+ const bool include;
+ };
+} // namespace
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+std::pair<unsigned int, unsigned int>
+MatrixFree<dim, Number, VectorizedArrayType>::
+ create_inner_face_subrange_hp_by_index(
+ const std::pair<unsigned int, unsigned int> &range,
+ const unsigned int fe_index_interior,
+ const unsigned int fe_index_exterior,
+ const unsigned int dof_handler_index) const
+{
+ if (dof_info[dof_handler_index].max_fe_index == 0 ||
+ dof_handlers[dof_handler_index]->get_fe_collection().size() == 1)
+ return range;
+
+ AssertIndexRange(fe_index_interior, dof_info[dof_handler_index].max_fe_index);
+ AssertIndexRange(fe_index_exterior, dof_info[dof_handler_index].max_fe_index);
+ const std::vector<unsigned int> &fe_indices =
+ dof_info[dof_handler_index].cell_active_fe_index;
+ if (fe_indices.empty() == true)
+ return range;
+ else
+ {
+ std::pair<unsigned int, unsigned int> return_range;
+ return_range.first =
+ std::lower_bound(face_info.faces.begin() + range.first,
+ face_info.faces.begin() + range.second,
+ std::pair<unsigned int, unsigned int>{
+ fe_index_interior, fe_index_exterior},
+ FaceRangeCompartor(fe_indices, false)) -
+ face_info.faces.begin();
+ return_range.second =
+ std::lower_bound(face_info.faces.begin() + return_range.first,
+ face_info.faces.begin() + range.second,
+ std::pair<unsigned int, unsigned int>{
+ fe_index_interior, fe_index_exterior},
+ FaceRangeCompartor(fe_indices, true)) -
+ face_info.faces.begin();
+ Assert(return_range.first >= range.first &&
+ return_range.second <= range.second,
+ ExcInternalError());
+ return return_range;
+ }
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+std::pair<unsigned int, unsigned int>
+MatrixFree<dim, Number, VectorizedArrayType>::
+ create_boundary_face_subrange_hp_by_index(
+ const std::pair<unsigned int, unsigned int> &range,
+ const unsigned int fe_index,
+ const unsigned int dof_handler_index) const
+{
+ if (dof_info[dof_handler_index].max_fe_index == 0 ||
+ dof_handlers[dof_handler_index]->get_fe_collection().size() == 1)
+ return range;
+
+ AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index);
+ const std::vector<unsigned int> &fe_indices =
+ dof_info[dof_handler_index].cell_active_fe_index;
+ if (fe_indices.empty() == true)
+ return range;
+ else
+ {
+ std::pair<unsigned int, unsigned int> return_range;
+ return_range.first =
+ std::lower_bound(face_info.faces.begin() + range.first,
+ face_info.faces.begin() + range.second,
+ fe_index,
+ FaceRangeCompartor(fe_indices, false)) -
+ face_info.faces.begin();
+ return_range.second =
+ std::lower_bound(face_info.faces.begin() + return_range.first,
+ face_info.faces.begin() + range.second,
+ fe_index,
+ FaceRangeCompartor(fe_indices, true)) -
+ face_info.faces.begin();
+ Assert(return_range.first >= range.first &&
+ return_range.second <= range.second,
+ ExcInternalError());
+ return return_range;
+ }
+}
+
+
+
template <int dim, typename Number, typename VectorizedArrayType>
void
MatrixFree<dim, Number, VectorizedArrayType>::renumber_dofs(
std::vector<types::global_dof_index> &renumbering,
- const unsigned int vector_component)
+ const unsigned int dof_handler_index)
{
- AssertIndexRange(vector_component, dof_info.size());
- dof_info[vector_component].compute_dof_renumbering(renumbering);
+ AssertIndexRange(dof_handler_index, dof_info.size());
+ dof_info[dof_handler_index].compute_dof_renumbering(renumbering);
}
#ifdef DEAL_II_WITH_MPI
{
// non-buffering mode is only supported if the indices of all cells are
- // contiguous for all dof_info objects and hp is not enabled.
+ // contiguous for all dof_info objects.
bool is_non_buffering_sm_supported = true;
for (const auto &di : dof_info)
{
#include "../tests.h"
-#include "./simplex_grids.h"
-
using namespace dealii;
+namespace dealii
+{
+ namespace GridGenerator
+ {
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_simplices_mix(
+ Triangulation<dim, spacedim> & tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<dim> & p1,
+ const Point<dim> & p2,
+ const bool colorize = false)
+ {
+ AssertDimension(dim, spacedim);
+
+ AssertThrow(colorize == false, ExcNotImplemented());
+
+ std::vector<Point<spacedim>> vertices;
+ std::vector<CellData<dim>> cells;
+
+ if (dim == 2)
+ {
+ // determine cell sizes
+ const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+ (p2[1] - p1[1]) / repetitions[1]);
+
+ // create vertices
+ for (unsigned int j = 0; j <= repetitions[1]; ++j)
+ for (unsigned int i = 0; i <= repetitions[0]; ++i)
+ vertices.push_back(
+ Point<spacedim>(p1[0] + dx[0] * i, p1[1] + dx[1] * j));
+
+ // create cells
+ for (unsigned int j = 0; j < repetitions[1]; ++j)
+ for (unsigned int i = 0; i < repetitions[0]; ++i)
+ {
+ // create reference QUAD cell
+ std::array<unsigned int, 4> quad{{
+ (j + 0) * (repetitions[0] + 1) + i + 0, //
+ (j + 0) * (repetitions[0] + 1) + i + 1, //
+ (j + 1) * (repetitions[0] + 1) + i + 0, //
+ (j + 1) * (repetitions[0] + 1) + i + 1 //
+ }}; //
+
+ if (j < repetitions[1] / 2 && i < repetitions[0] / 2)
+ {
+ CellData<dim> quad_;
+ quad_.vertices = {quad[0], quad[1], quad[2], quad[3]};
+ cells.push_back(quad_);
+
+ continue;
+ }
+
+ // TRI cell 0
+ {
+ CellData<dim> tri;
+ tri.vertices = {quad[0], quad[1], quad[2]};
+ cells.push_back(tri);
+ }
+
+ // TRI cell 1
+ {
+ CellData<dim> tri;
+ tri.vertices = {quad[3], quad[2], quad[1]};
+ cells.push_back(tri);
+ }
+ }
+ }
+ else
+ {
+ AssertThrow(colorize == false, ExcNotImplemented());
+ }
+
+ // actually create triangulation
+ tria.create_triangulation(vertices, cells, SubCellData());
+ }
+
+
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_cube_with_simplices_mix(Triangulation<dim, spacedim> &tria,
+ const unsigned int repetitions,
+ const double p1 = 0.0,
+ const double p2 = 1.0,
+ const bool colorize = false)
+ {
+ if (dim == 2)
+ {
+ subdivided_hyper_rectangle_with_simplices_mix(
+ tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize);
+ }
+ else if (dim == 3)
+ {
+ subdivided_hyper_rectangle_with_simplices_mix(
+ tria,
+ {{repetitions, repetitions, repetitions}},
+ {p1, p1, p1},
+ {p2, p2, p2},
+ colorize);
+ }
+ else
+ {
+ AssertThrow(false, ExcNotImplemented())
+ }
+ }
+ } // namespace GridGenerator
+} // namespace dealii
+
template <int dim>
class PoissonOperator
{
public:
- using VectorType = LinearAlgebra::distributed::Vector<double>;
+ using VectorType = LinearAlgebra::distributed::Vector<double>;
+ using FECellIntegrator = FEEvaluation<dim, -1, 0, 1, double>;
PoissonOperator(const MatrixFree<dim, double> &matrix_free,
const bool do_helmholtz)
const int dummy = 0;
matrix_free.template cell_loop<VectorType, int>(
- [&](const auto &data, auto &dst, const auto &, const auto cell_range) {
- for (unsigned int i = 0; i < 2; ++i)
- {
- const auto cell_subrange =
- data.create_cell_subrange_hp_by_index(cell_range, i);
+ [&](const auto &data, auto &dst, const auto &, const auto range) {
+ const auto i = data.get_cell_active_fe_index(range);
+ FECellIntegrator phi(matrix_free, 0, 0, 0, i, i);
- FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 0, 0, 0, i, i);
-
- for (unsigned int cell = cell_subrange.first;
- cell < cell_subrange.second;
- ++cell)
- {
- phi.reinit(cell);
- for (unsigned int q = 0; q < phi.n_q_points; ++q)
- phi.submit_value(1.0, q);
+ for (unsigned int cell = range.first; cell < range.second; ++cell)
+ {
+ phi.reinit(cell);
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ phi.submit_value(1.0, q);
- phi.integrate_scatter(true, false, dst);
- }
+ phi.integrate_scatter(true, false, dst);
}
},
vec,
vmult(VectorType &dst, const VectorType &src) const
{
matrix_free.template cell_loop<VectorType, VectorType>(
- [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
- for (unsigned int i = 0; i < 2; ++i)
+ [&](const auto &data, auto &dst, const auto &src, const auto range) {
+ const auto i = data.get_cell_active_fe_index(range);
+ FECellIntegrator phi(matrix_free, 0, 0, 0, i, i);
+
+ for (unsigned int cell = range.first; cell < range.second; ++cell)
{
- const auto cell_subrange =
- data.create_cell_subrange_hp_by_index(cell_range, i);
+ phi.reinit(cell);
+ phi.gather_evaluate(src, do_helmholtz, true);
- FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 0, 0, 0, i, i);
- for (unsigned int cell = cell_subrange.first;
- cell < cell_subrange.second;
- ++cell)
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
- phi.reinit(cell);
- phi.gather_evaluate(src, do_helmholtz, true);
+ if (do_helmholtz)
+ phi.submit_value(phi.get_value(q), q);
- for (unsigned int q = 0; q < phi.n_q_points; ++q)
- {
- if (do_helmholtz)
- phi.submit_value(phi.get_value(q), q);
-
- phi.submit_gradient(phi.get_gradient(q), q);
- }
-
- phi.integrate_scatter(do_helmholtz, true, dst);
+ phi.submit_gradient(phi.get_gradient(q), q);
}
+
+ phi.integrate_scatter(do_helmholtz, true, dst);
}
},
dst,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test ShapeData for Simplex::FE_P and Simplex::QGauss
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <ios>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test()
+{
+ const unsigned int degree = 1;
+
+ Triangulation<dim> tria;
+ GridGenerator::subdivided_hyper_cube(tria, 4);
+
+ hp::FECollection<dim> fe{FE_Q<2>(degree), FE_Q<2>(degree)};
+ MappingQGeneric<dim> mapping(1);
+ QGauss<dim> quadrature(degree + 1);
+
+ DoFHandler<dim> dof_handler(tria);
+
+ unsigned int counter = 0;
+
+ for (auto &cell : dof_handler)
+ if (counter++ < tria.n_cells() / 2)
+ cell.set_active_fe_index(1);
+
+ dof_handler.distribute_dofs(fe);
+
+
+ typename MatrixFree<dim, double>::AdditionalData additional_data;
+ additional_data.mapping_update_flags_boundary_faces =
+ update_gradients | update_values;
+ additional_data.mapping_update_flags_inner_faces =
+ update_gradients | update_values;
+ additional_data.initialize_mapping = false;
+ additional_data.tasks_parallel_scheme =
+ MatrixFree<dim, double>::AdditionalData::none;
+
+ AffineConstraints<double> constraints;
+
+ MatrixFree<dim, double> matrix_free;
+ matrix_free.reinit(
+ mapping, dof_handler, constraints, quadrature, additional_data);
+
+ using VectorType = Vector<double>;
+
+ VectorType src, dst;
+ matrix_free.initialize_dof_vector(src);
+ matrix_free.initialize_dof_vector(dst);
+
+ std::vector<std::vector<CellId>> cells(2);
+ std::vector<std::vector<std::pair<CellId, unsigned int>>> boundary_faces(2);
+
+ std::vector<std::vector<std::vector<std::pair<CellId, CellId>>>> inner_faces(
+ 2);
+ for (unsigned int i = 0; i < 2; ++i)
+ inner_faces[i].resize(2);
+
+ matrix_free.template loop<VectorType, VectorType>(
+ [&](const auto &data, auto &dst, const auto &src, const auto range) {
+ const auto i = data.get_cell_active_fe_index(range);
+
+ for (unsigned int cell = range.first; cell < range.second; ++cell)
+ for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
+ ++v)
+ cells[i].push_back(data.get_cell_iterator(cell, v)->id());
+ },
+ [&](const auto &data, auto &dst, const auto &src, const auto range) {
+ const auto i = data.get_face_active_fe_index(range, true);
+ const auto j = data.get_face_active_fe_index(range, false);
+
+ for (unsigned int face = range.first; face < range.second; ++face)
+ for (unsigned int v = 0; v < data.n_active_entries_per_face_batch(face);
+ ++v)
+ inner_faces[i][j].emplace_back(
+ data.get_face_iterator(face, v, true).first->id(),
+ data.get_face_iterator(face, v, false).first->id());
+ },
+ [&](const auto &data, auto &dst, const auto &src, const auto range) {
+ const auto i = data.get_face_active_fe_index(range);
+
+ for (unsigned int face = range.first; face < range.second; ++face)
+ for (unsigned int v = 0; v < data.n_active_entries_per_face_batch(face);
+ ++v)
+ boundary_faces[i].emplace_back(
+ data.get_face_iterator(face, v).first->id(),
+ data.get_face_iterator(face, v).second);
+ },
+ dst,
+ src);
+
+ deallog << "CELL categories:" << std::endl;
+ for (auto &i : cells)
+ {
+ std::sort(i.begin(), i.end());
+
+ for (const auto &j : i)
+ deallog << j << " ";
+ deallog << std::endl;
+ }
+
+ deallog << "BOUNDARY_FACE categories:" << std::endl;
+ for (auto &i : boundary_faces)
+ {
+ std::sort(i.begin(), i.end());
+ for (const auto &j : i)
+ deallog << j.first << "@" << j.second << " ";
+ deallog << std::endl;
+ }
+
+ deallog << "INNER_FACE categories:" << std::endl;
+ for (auto &k : inner_faces)
+ {
+ for (auto &i : k)
+ {
+ std::sort(i.begin(), i.end());
+ for (const auto &j : i)
+ deallog << j.first << "@" << j.second << " ";
+ deallog << std::endl;
+ }
+ }
+}
+
+int
+main()
+{
+ initlog();
+
+ test<2>();
+}
--- /dev/null
+
+DEAL::CELL categories:
+DEAL::8_0: 9_0: 10_0: 11_0: 12_0: 13_0: 14_0: 15_0:
+DEAL::0_0: 1_0: 2_0: 3_0: 4_0: 5_0: 6_0: 7_0:
+DEAL::BOUNDARY_FACE categories:
+DEAL::8_0:@0 11_0:@1 12_0:@0 12_0:@3 13_0:@3 14_0:@3 15_0:@1 15_0:@3
+DEAL::0_0:@0 0_0:@2 1_0:@2 2_0:@2 3_0:@1 3_0:@2 4_0:@0 7_0:@1
+DEAL::INNER_FACE categories:
+DEAL::9_0:@8_0: 10_0:@9_0: 11_0:@10_0: 12_0:@8_0: 13_0:@9_0: 13_0:@12_0: 14_0:@10_0: 14_0:@13_0: 15_0:@11_0: 15_0:@14_0:
+DEAL::
+DEAL::4_0:@8_0: 5_0:@9_0: 6_0:@10_0: 7_0:@11_0:
+DEAL::1_0:@0_0: 2_0:@1_0: 3_0:@2_0: 4_0:@0_0: 5_0:@1_0: 5_0:@4_0: 6_0:@2_0: 6_0:@5_0: 7_0:@3_0: 7_0:@6_0: