--- /dev/null
+Improved: The `Meshworker::mesh_loop()` function is now capable of working with an
+IteratorRange, and also supports iterator ranges constructed from
+FilteredIterators.
+<br>
+(Jean-Paul Pelteret, 2019/05/02)
#include <deal.II/meshworker/loop.h>
#include <functional>
+#include <type_traits>
DEAL_II_NAMESPACE_OPEN
namespace MeshWorker
{
+ namespace internal
+ {
+ /**
+ * A helper class to provide a type definition for the underlying cell
+ * iterator type.
+ */
+ template <class CellIteratorType>
+ struct CellIteratorBaseType
+ {
+ /**
+ * Type definition for the cell iterator type.
+ */
+ using type = CellIteratorType;
+ };
+
+ /**
+ * A helper class to provide a type definition for the underlying cell
+ * iterator type.
+ *
+ * This specialization is for IteratorRange, which may have either a
+ * TriaActiveIterator or a FilteredIterator as its base type.
+ */
+ template <class CellIteratorType>
+ struct CellIteratorBaseType<
+ typename IteratorRange<CellIteratorType>::IteratorOverIterators>
+ {
+ /**
+ * Type definition for the cell iterator type.
+ */
+ using type = typename IteratorRange<
+ CellIteratorType>::IteratorOverIterators::BaseIterator;
+
+ static_assert(std::is_same<type, CellIteratorType>::value,
+ "Expected another type as the base iterator.");
+ };
+
+ /**
+ * A helper class to provide a type definition for the underlying cell
+ * iterator type.
+ *
+ * This specialization is for FilteredIterator, which may have either a
+ * TriaActiveIterator as its base type, or may be nested with another
+ * FilteredIterator as the type to iterate over.
+ */
+ template <class CellIteratorType>
+ struct CellIteratorBaseType<FilteredIterator<CellIteratorType>>
+ {
+ /**
+ * Type definition for the cell iterator type.
+ */
+ // Since we can have nested filtered iterators, we recursively
+ // remove the template layers to retrieve the underlying iterator type.
+ using type = typename CellIteratorBaseType<CellIteratorType>::type;
+ };
+
+ /**
+ * A helper class to provide a type definition for the underlying cell
+ * iterator type.
+ *
+ * This specialization is for an IteratorRange that iterates over
+ * FilteredIterators.
+ */
+ template <class CellIteratorType>
+ struct CellIteratorBaseType<typename IteratorRange<
+ FilteredIterator<CellIteratorType>>::IteratorOverIterators>
+ {
+ /**
+ * Type definition for the cell iterator type.
+ */
+ // Since we can have nested filtered iterators, we recursively
+ // remove the template layers to retrieve the underlying iterator type.
+ using type =
+ typename CellIteratorBaseType<FilteredIterator<CellIteratorType>>::type;
+ };
+ } // namespace internal
+
/**
* This function extends the WorkStream concept to work on meshes
* (cells and/or faces) and handles the complicated logic for
* @ingroup MeshWorker
* @author Luca Heltai and Timo Heister, 2017
*/
- template <class CellIteratorType, class ScratchData, class CopyData>
+ template <class CellIteratorType,
+ class ScratchData,
+ class CopyData,
+ class CellIteratorBaseType =
+ typename internal::CellIteratorBaseType<CellIteratorType>::type>
void
- mesh_loop(const CellIteratorType & begin,
- const typename identity<CellIteratorType>::type &end,
+ mesh_loop(
+ const CellIteratorType & begin,
+ const typename identity<CellIteratorType>::type &end,
- const typename identity<std::function<
- void(const CellIteratorType &, ScratchData &, CopyData &)>>::type
- &cell_worker,
- const typename identity<std::function<void(const CopyData &)>>::type
- &copier,
-
- const ScratchData &sample_scratch_data,
- const CopyData & sample_copy_data,
-
- const AssembleFlags flags = assemble_own_cells,
-
- const typename identity<std::function<void(const CellIteratorType &,
- const unsigned int,
- ScratchData &,
- CopyData &)>>::type
- &boundary_worker = std::function<void(const CellIteratorType &,
- const unsigned int,
- ScratchData &,
- CopyData &)>(),
-
- const typename identity<std::function<void(const CellIteratorType &,
- const unsigned int,
- const unsigned int,
- const CellIteratorType &,
- const unsigned int,
- const unsigned int,
- ScratchData &,
- CopyData &)>>::type
- &face_worker = std::function<void(const CellIteratorType &,
- const unsigned int,
- const unsigned int,
- const CellIteratorType &,
- const unsigned int,
- const unsigned int,
- ScratchData &,
- CopyData &)>(),
+ const typename identity<std::function<
+ void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
+ &cell_worker,
+ const typename identity<std::function<void(const CopyData &)>>::type
+ &copier,
- const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
- const unsigned int chunk_size = 8)
+ const ScratchData &sample_scratch_data,
+ const CopyData & sample_copy_data,
+
+ const AssembleFlags flags = assemble_own_cells,
+
+ const typename identity<std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>>::type
+ &boundary_worker = std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>(),
+
+ const typename identity<std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>>::type
+ &face_worker = std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>(),
+
+ const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
+ const unsigned int chunk_size = 8)
{
Assert(
(!cell_worker) == !(flags & work_on_cells),
ExcMessage(
"If you specify a boundary_worker, assemble_boundary_faces needs to be set."));
- auto cell_action = [&](const CellIteratorType &cell,
- ScratchData & scratch,
- CopyData & copy) {
+ auto cell_action = [&](const CellIteratorBaseType &cell,
+ ScratchData & scratch,
+ CopyData & copy) {
// First reset the CopyData class to the empty copy_data given by the
// user.
copy = sample_copy_data;
if (flags & (work_on_faces | work_on_boundary))
for (unsigned int face_no = 0;
- face_no < GeometryInfo<CellIteratorType::AccessorType::Container::
- dimension>::faces_per_cell;
+ face_no < GeometryInfo<CellIteratorBaseType::AccessorType::
+ Container::dimension>::faces_per_cell;
++face_no)
{
if (cell->at_boundary(face_no) &&
else
{
// interior face, potentially assemble
- TriaIterator<typename CellIteratorType::AccessorType> neighbor =
- cell->neighbor_or_periodic_neighbor(face_no);
+ TriaIterator<typename CellIteratorBaseType::AccessorType>
+ neighbor = cell->neighbor_or_periodic_neighbor(face_no);
types::subdomain_id neighbor_subdomain_id =
numbers::artificial_subdomain_id;
{
// If iterator is active and neighbor is refined, skip
// internal face.
- if (internal::is_active_iterator(cell) &&
+ if (dealii::internal::is_active_iterator(cell) &&
neighbor->has_children())
continue;
chunk_size);
}
+ /**
+ * Same as the function above, but for iterator ranges (and, therefore,
+ * filtered iterators).
+ */
+ template <class CellIteratorType,
+ class ScratchData,
+ class CopyData,
+ class CellIteratorBaseType =
+ typename internal::CellIteratorBaseType<CellIteratorType>::type>
+ void
+ mesh_loop(
+ IteratorRange<CellIteratorType> iterator_range,
+ const typename identity<std::function<
+ void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
+ &cell_worker,
+ const typename identity<std::function<void(const CopyData &)>>::type
+ &copier,
+
+ const ScratchData &sample_scratch_data,
+ const CopyData & sample_copy_data,
+
+ const AssembleFlags flags = assemble_own_cells,
+
+ const typename identity<std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>>::type
+ &boundary_worker = std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>(),
+
+ const typename identity<std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>>::type
+ &face_worker = std::function<void(const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ ScratchData &,
+ CopyData &)>(),
+
+ const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
+ const unsigned int chunk_size = 8)
+ {
+ // Call the function above
+ mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
+ ScratchData,
+ CopyData,
+ CellIteratorBaseType>(iterator_range.begin(),
+ iterator_range.end(),
+ cell_worker,
+ copier,
+ sample_scratch_data,
+ sample_copy_data,
+ flags,
+ boundary_worker,
+ face_worker,
+ queue_length,
+ chunk_size);
+ }
+
/**
* This is a variant of the mesh_loop() function, that can be used for worker
* and copier functions that are member functions of a class.
queue_length,
chunk_size);
}
+
+ /**
+ * Same as the function above, but for iterator ranges (and, therefore,
+ * filtered iterators).
+ */
+ template <class CellIteratorType,
+ class ScratchData,
+ class CopyData,
+ class MainClass,
+ class CellIteratorBaseType =
+ typename internal::CellIteratorBaseType<CellIteratorType>::type>
+ void
+ mesh_loop(IteratorRange<CellIteratorType> iterator_range,
+ MainClass & main_class,
+ void (MainClass::*cell_worker)(const CellIteratorBaseType &,
+ ScratchData &,
+ CopyData &),
+ void (MainClass::*copier)(const CopyData &),
+ const ScratchData & sample_scratch_data,
+ const CopyData & sample_copy_data,
+ const AssembleFlags flags = assemble_own_cells,
+ void (MainClass::*boundary_worker)(const CellIteratorBaseType &,
+ const unsigned int,
+ ScratchData &,
+ CopyData &) = nullptr,
+ void (MainClass::*face_worker)(const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ const CellIteratorBaseType &,
+ const unsigned int,
+ const unsigned int,
+ ScratchData &,
+ CopyData &) = nullptr,
+ const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
+ const unsigned int chunk_size = 8)
+ {
+ // Call the function above
+ mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
+ ScratchData,
+ CopyData,
+ MainClass,
+ CellIteratorBaseType>(iterator_range.begin(),
+ iterator_range.end(),
+ main_class,
+ cell_worker,
+ copier,
+ sample_scratch_data,
+ sample_copy_data,
+ flags,
+ boundary_worker,
+ face_worker,
+ queue_length,
+ chunk_size);
+ }
} // namespace MeshWorker
DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 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 mesh_loop for iterator ranges
+// This test is based off meshworker/scratch_data_01.cc, but for iterator ranges
+
+#include <deal.II/base/iterator_range.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.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_Q<dim, spacedim> fe(1);
+ DoFHandler<dim, spacedim> dh(tria);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(1);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ dh.distribute_dofs(fe);
+
+ QGauss<dim> quad(3);
+ QGauss<dim - 1> face_quad(3);
+
+ UpdateFlags cell_flags = update_JxW_values;
+ UpdateFlags face_flags = update_JxW_values;
+
+ using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
+ using CopyData = MeshWorker::CopyData<0, 1, 0>;
+
+ ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags);
+ CopyData copy(3); // store measures
+
+ using Iterator = decltype(dh.begin_active());
+
+ auto measures = std::make_tuple(0.0, 0.0, 0.0);
+
+
+ auto cell_worker = [](const Iterator &cell, ScratchData &s, CopyData &c) {
+ const auto &fev = s.reinit(cell);
+ const auto &JxW = s.get_JxW_values();
+ c = 0;
+ for (auto w : JxW)
+ c.vectors[0][0] += w;
+ };
+
+ auto boundary_worker = [](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();
+ for (auto w : JxW)
+ c.vectors[0][1] += w;
+ };
+
+ auto face_worker = [](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);
+ const auto &nfev = s.reinit_neighbor(ncell, nf, nsf);
+
+ const auto &JxW = s.get_JxW_values();
+ const auto &nJxW = s.get_neighbor_JxW_values();
+ for (auto w : JxW)
+ c.vectors[0][2] += w;
+ };
+
+ auto copier = [&measures](const CopyData &c) {
+ std::get<0>(measures) += c.vectors[0][0];
+ std::get<1>(measures) += c.vectors[0][1];
+ std::get<2>(measures) += c.vectors[0][2];
+ };
+
+ mesh_loop(dh.active_cell_iterators(),
+ cell_worker,
+ copier,
+ scratch,
+ copy,
+ assemble_own_cells | assemble_boundary_faces |
+ assemble_own_interior_faces_once,
+ boundary_worker,
+ face_worker);
+
+ deallog << "Testing <" << dim << "," << spacedim << ">" << std::endl;
+
+ deallog << "Volume: " << std::get<0>(measures) << std::endl;
+
+ deallog << "Boundary surface: " << std::get<1>(measures) << std::endl;
+
+ deallog << "Interior surface: " << std::get<2>(measures) << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ test<2, 2>();
+ test<2, 3>();
+ test<3, 3>();
+}
--- /dev/null
+
+DEAL::Testing <2,2>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 3.00000
+DEAL::Testing <2,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 3.00000
+DEAL::Testing <3,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 6.00000
+DEAL::Interior surface: 3.75000
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 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 mesh_loop for filtered iterator ranges
+// This test is based off base/work_stream_06.cc, but for parallel distributed
+// Triangulations.
+
+#include <deal.II/base/iterator_range.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <unordered_map>
+
+#include "../tests.h"
+
+using namespace MeshWorker;
+
+template <int dim, int spacedim>
+void
+test()
+{
+ parallel::distributed::Triangulation<dim, spacedim> tria(
+ MPI_COMM_WORLD,
+ Triangulation<dim, spacedim>::limit_level_difference_at_vertices);
+ FE_Q<dim, spacedim> fe(1);
+ DoFHandler<dim, spacedim> dh(tria);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(1);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ dh.distribute_dofs(fe);
+
+ QGauss<dim> quad(3);
+ QGauss<dim - 1> face_quad(3);
+
+ UpdateFlags cell_flags = update_JxW_values;
+ UpdateFlags face_flags = update_JxW_values;
+
+ using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
+ using CopyData = MeshWorker::CopyData<0, 1, 0>;
+
+ ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags);
+ CopyData copy(3); // store measures
+
+ using Iterator = decltype(dh.begin_active());
+
+ auto measures = std::make_tuple(0.0, 0.0, 0.0);
+
+
+ auto cell_worker = [](const Iterator &cell, ScratchData &s, CopyData &c) {
+ const auto &fev = s.reinit(cell);
+ const auto &JxW = s.get_JxW_values();
+ c = 0;
+ for (auto w : JxW)
+ c.vectors[0][0] += w;
+ };
+
+ auto boundary_worker = [](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();
+ for (auto w : JxW)
+ c.vectors[0][1] += w;
+ };
+
+ auto face_worker = [](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);
+ const auto &nfev = s.reinit_neighbor(ncell, nf, nsf);
+
+ const auto &JxW = s.get_JxW_values();
+ const auto &nJxW = s.get_neighbor_JxW_values();
+ for (auto w : JxW)
+ c.vectors[0][2] += w;
+ };
+
+ auto copier = [&measures](const CopyData &c) {
+ std::get<0>(measures) += c.vectors[0][0];
+ std::get<1>(measures) += c.vectors[0][1];
+ std::get<2>(measures) += c.vectors[0][2];
+ };
+
+ const auto filtered_iterator_range =
+ filter_iterators(dh.active_cell_iterators(),
+ IteratorFilters::LocallyOwnedCell());
+
+ mesh_loop(filtered_iterator_range,
+ cell_worker,
+ copier,
+ scratch,
+ copy,
+ assemble_own_cells | assemble_boundary_faces |
+ assemble_own_interior_faces_once,
+ boundary_worker,
+ face_worker);
+
+ std::get<0>(measures) =
+ Utilities::MPI::sum(std::get<0>(measures), MPI_COMM_WORLD);
+ std::get<1>(measures) =
+ Utilities::MPI::sum(std::get<1>(measures), MPI_COMM_WORLD);
+ std::get<2>(measures) =
+ Utilities::MPI::sum(std::get<2>(measures), MPI_COMM_WORLD);
+
+ deallog << "Testing <" << dim << "," << spacedim << ">" << std::endl;
+
+ deallog << "Volume: " << std::get<0>(measures) << std::endl;
+
+ deallog << "Boundary surface: " << std::get<1>(measures) << std::endl;
+
+ deallog << "Interior surface: " << std::get<2>(measures) << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+ mpi_initlog();
+
+ test<2, 2>();
+ test<2, 3>();
+ test<3, 3>();
+}
--- /dev/null
+
+DEAL::Testing <2,2>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 2.00000
+DEAL::Testing <2,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 2.00000
+DEAL::Testing <3,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 6.00000
+DEAL::Interior surface: 3.00000
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 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 mesh_loop for filtered iterator ranges
+// This test is based off base/work_stream_06.cc, but for filtered iterator
+// ranges. In particular, we test the material ID predicate.
+
+#include <deal.II/base/iterator_range.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <unordered_map>
+
+#include "../tests.h"
+
+using namespace MeshWorker;
+
+template <int dim, int spacedim>
+void
+test(const types::material_id test_id)
+{
+ Assert(test_id == 0 || test_id == 1,
+ ExcMessage("Only material IDs 0 and 1 are valid."));
+
+ Triangulation<dim, spacedim> tria;
+ FE_Q<dim, spacedim> fe(1);
+ DoFHandler<dim, spacedim> dh(tria);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(1);
+ tria.begin_active()->set_material_id(1);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ dh.distribute_dofs(fe);
+
+ QGauss<dim> quad(3);
+ QGauss<dim - 1> face_quad(3);
+
+ UpdateFlags cell_flags = update_JxW_values;
+ UpdateFlags face_flags = update_JxW_values;
+
+ using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
+ using CopyData = MeshWorker::CopyData<0, 1, 0>;
+
+ ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags);
+ CopyData copy(3); // store measures
+
+ using Iterator = decltype(dh.begin_active());
+
+ auto measures = std::make_tuple(0.0, 0.0, 0.0);
+
+
+ auto cell_worker = [](const Iterator &cell, ScratchData &s, CopyData &c) {
+ const auto &fev = s.reinit(cell);
+ const auto &JxW = s.get_JxW_values();
+ c = 0;
+ for (auto w : JxW)
+ c.vectors[0][0] += w;
+ };
+
+ auto boundary_worker = [](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();
+ for (auto w : JxW)
+ c.vectors[0][1] += w;
+ };
+
+ auto face_worker = [](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);
+ const auto &nfev = s.reinit_neighbor(ncell, nf, nsf);
+
+ const auto &JxW = s.get_JxW_values();
+ const auto &nJxW = s.get_neighbor_JxW_values();
+ for (auto w : JxW)
+ c.vectors[0][2] += w;
+ };
+
+ auto copier = [&measures](const CopyData &c) {
+ std::get<0>(measures) += c.vectors[0][0];
+ std::get<1>(measures) += c.vectors[0][1];
+ std::get<2>(measures) += c.vectors[0][2];
+ };
+
+ const auto filtered_iterator_range =
+ filter_iterators(dh.active_cell_iterators(),
+ IteratorFilters::MaterialIdEqualTo(test_id));
+
+ mesh_loop(filtered_iterator_range,
+ cell_worker,
+ copier,
+ scratch,
+ copy,
+ assemble_own_cells | assemble_boundary_faces |
+ assemble_own_interior_faces_once,
+ boundary_worker,
+ face_worker);
+
+ deallog << "Testing <" << dim << "," << spacedim << ">" << std::endl;
+
+ deallog << "Volume: " << std::get<0>(measures) << std::endl;
+
+ deallog << "Boundary surface: " << std::get<1>(measures) << std::endl;
+
+ deallog << "Interior surface: " << std::get<2>(measures) << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ for (auto mat_id : {0, 1})
+ {
+ deallog.push("Material ID " + Utilities::to_string(mat_id));
+ {
+ test<2, 2>(mat_id);
+ test<2, 3>(mat_id);
+ test<3, 3>(mat_id);
+ }
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:Material ID 0::Testing <2,2>
+DEAL:Material ID 0::Volume: 0.750000
+DEAL:Material ID 0::Boundary surface: 3.00000
+DEAL:Material ID 0::Interior surface: 1.00000
+DEAL:Material ID 0::Testing <2,3>
+DEAL:Material ID 0::Volume: 0.750000
+DEAL:Material ID 0::Boundary surface: 3.00000
+DEAL:Material ID 0::Interior surface: 1.00000
+DEAL:Material ID 0::Testing <3,3>
+DEAL:Material ID 0::Volume: 0.875000
+DEAL:Material ID 0::Boundary surface: 5.25000
+DEAL:Material ID 0::Interior surface: 2.25000
+DEAL:Material ID 1::Testing <2,2>
+DEAL:Material ID 1::Volume: 0.250000
+DEAL:Material ID 1::Boundary surface: 1.00000
+DEAL:Material ID 1::Interior surface: 2.00000
+DEAL:Material ID 1::Testing <2,3>
+DEAL:Material ID 1::Volume: 0.250000
+DEAL:Material ID 1::Boundary surface: 1.00000
+DEAL:Material ID 1::Interior surface: 2.00000
+DEAL:Material ID 1::Testing <3,3>
+DEAL:Material ID 1::Volume: 0.125000
+DEAL:Material ID 1::Boundary surface: 0.750000
+DEAL:Material ID 1::Interior surface: 1.50000
projects / dealii.git / commitdiff