--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// Check that ScratchData returns the correct jumps in, and averages of,
+// solution values, gradients, etc. at interfaces.
+// - Scalar valued FE
+// - hp variant
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+template <int dim>
+class DiscontinuousFunction : public Function<dim>
+{
+public:
+ DiscontinuousFunction(const unsigned int n_components = 1)
+ : Function<dim>(n_components)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override
+ {
+ const double value = (p[0] * p[1] > 0 ? +3 : -1); // Non-trivial average
+ return value * (1.0 + std::sin(p[0]) * std::cos(p[1]));
+ }
+};
+
+
+template <int dim,
+ int spacedim = dim,
+ typename NumberType = double,
+ typename ExtractorType>
+void
+run(const ExtractorType &extractor)
+{
+ LogStream::Prefix prefix("Dim " + Utilities::to_string(dim));
+ std::cout << "Dim: " << dim << std::endl;
+
+ hp::FECollection<dim, spacedim> fe(FE_DGQ<dim, spacedim>(3),
+ FE_DGQ<dim, spacedim>(4));
+ hp::QCollection<dim> qf_cell(QGauss<dim>(fe.max_degree()),
+ QGauss<dim>(fe.max_degree() + 1));
+ hp::QCollection<dim - 1> qf_face(QGauss<dim - 1>(fe.max_degree()),
+ QGauss<dim - 1>(fe.max_degree() + 1));
+
+ Triangulation<dim, spacedim> triangulation;
+ GridGenerator::hyper_cube(triangulation, -1, 1);
+ triangulation.refine_global(1);
+
+ DoFHandler<dim, spacedim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> solution(dof_handler.n_dofs());
+ // Don't interpolate, as points precisely on the interfaces get evaluated and
+ // we don't end up with a jump in the values across it.
+ AffineConstraints<double> constraints;
+ constraints.close();
+ VectorTools::project(
+ dof_handler,
+ constraints,
+ hp::QCollection<dim>(QGauss<spacedim>(fe.max_degree() + 2),
+ QGauss<spacedim>(fe.max_degree() + 3)),
+ DiscontinuousFunction<spacedim>(fe.n_components()),
+ solution);
+
+ const UpdateFlags update_flags =
+ update_values | update_gradients | update_hessians | update_3rd_derivatives;
+ MeshWorker::ScratchData<dim, spacedim> scratch_data(
+ fe, qf_cell, update_flags, qf_face, update_flags);
+
+ const auto cell = dof_handler.begin_active();
+ unsigned int face = 0;
+ while (cell->face(face)->at_boundary())
+ ++face;
+
+ scratch_data.reinit(cell,
+ face,
+ cell->neighbor(face),
+ cell->neighbor_face_no(face));
+ scratch_data.extract_local_dof_values("solution", solution);
+
+ deallog << "Jumps in values: "
+ << scratch_data.get_jumps_in_values("solution", extractor)[0]
+ << std::endl;
+ deallog << "Jumps in gradients: "
+ << scratch_data.get_jumps_in_gradients("solution", extractor)[0]
+ << std::endl;
+ deallog << "Jumps in Hessians: "
+ << scratch_data.get_jumps_in_hessians("solution", extractor)[0]
+ << std::endl;
+ deallog << "Jumps in third derivatives: "
+ << scratch_data.get_jumps_in_third_derivatives("solution",
+ extractor)[0]
+ << std::endl;
+
+ deallog << "Averages of values: "
+ << scratch_data.get_averages_of_values("solution", extractor)[0]
+ << std::endl;
+ deallog << "Averages of gradients: "
+ << scratch_data.get_averages_of_gradients("solution", extractor)[0]
+ << std::endl;
+ deallog << "Averages of Hessians: "
+ << scratch_data.get_averages_of_hessians("solution", extractor)[0]
+ << std::endl;
+
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main(int argc, char *argv[])
+{
+ initlog();
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, testing_max_num_threads());
+
+ FEValuesExtractors::Scalar extractor(0);
+
+ run<2>(extractor);
+ run<3>(extractor);
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL:Dim 2::Jumps in values: 4.00030
+DEAL:Dim 2::Jumps in gradients: 2.40085 0.000407768
+DEAL:Dim 2::Jumps in Hessians: 0.0456300 3.24213 3.24213 -0.000332829
+DEAL:Dim 2::Jumps in third derivatives: -2.06802 0.0616191 0.0616191 -2.64630 0.0616191 -2.64630 -2.64630 -0.000238992
+DEAL:Dim 2::Averages of values: 1.00030
+DEAL:Dim 2::Averages of gradients: 0.600214 0.000407768
+DEAL:Dim 2::Averages of Hessians: 0.0456300 0.810533 0.810533 -0.000332829
+DEAL:Dim 2::OK
+DEAL:Dim 3::Jumps in values: 4.00030
+DEAL:Dim 3::Jumps in gradients: 2.40085 0.000407768 -1.31631e-14
+DEAL:Dim 3::Jumps in Hessians: 0.0456300 3.24213 2.93071e-13 3.24213 -0.000332829 -6.39488e-13 2.93071e-13 -6.39488e-13 -2.08722e-13
+DEAL:Dim 3::Jumps in third derivatives: -2.06802 0.0616191 -4.40314e-13 0.0616191 -2.64630 -7.67264e-12 -4.40314e-13 -7.67264e-12 -3.20455e-12 0.0616191 -2.64630 -7.67264e-12 -2.64630 -0.000238992 4.35918e-12 -7.67264e-12 4.35918e-12 4.07763e-12 -4.40314e-13 -7.67264e-12 -3.20455e-12 -7.67264e-12 4.35918e-12 4.07763e-12 -3.20455e-12 4.07763e-12 7.54952e-13
+DEAL:Dim 3::Averages of values: 1.00030
+DEAL:Dim 3::Averages of gradients: 0.600214 0.000407768 5.69336e-15
+DEAL:Dim 3::Averages of Hessians: 0.0456300 0.810533 3.11001e-14 0.810533 -0.000332829 -4.40536e-13 3.11001e-14 -4.40536e-13 -1.72751e-13
+DEAL:Dim 3::OK
+DEAL::OK
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// Check that ScratchData returns the correct jumps in, and averages of,
+// solution values, gradients, etc. at interfaces.
+// - Vector valued FE
+// - hp variant
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+template <int dim>
+class DiscontinuousFunction : public Function<dim>
+{
+public:
+ DiscontinuousFunction(const unsigned int n_components = 1)
+ : Function<dim>(n_components)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override
+ {
+ const double value = (p[0] * p[1] > 0 ? +3 : -1); // Non-trivial average
+ return value * (1.0 + std::sin(p[0]) * std::cos(p[1]));
+ }
+};
+
+
+template <int dim,
+ int spacedim = dim,
+ typename NumberType = double,
+ typename ExtractorType>
+void
+run(const ExtractorType &extractor)
+{
+ LogStream::Prefix prefix("Dim " + Utilities::to_string(dim));
+ std::cout << "Dim: " << dim << std::endl;
+
+ hp::FECollection<dim, spacedim> fe(
+ FESystem<dim, spacedim>(FE_DGQ<dim, spacedim>(3), dim),
+ FESystem<dim, spacedim>(FE_DGQ<dim, spacedim>(4), dim));
+ hp::QCollection<dim> qf_cell(QGauss<dim>(fe.max_degree()),
+ QGauss<dim>(fe.max_degree() + 1));
+ hp::QCollection<dim - 1> qf_face(QGauss<dim - 1>(fe.max_degree()),
+ QGauss<dim - 1>(fe.max_degree() + 1));
+
+ Triangulation<dim, spacedim> triangulation;
+ GridGenerator::hyper_cube(triangulation, -1, 1);
+ triangulation.refine_global(1);
+
+ DoFHandler<dim, spacedim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> solution(dof_handler.n_dofs());
+ // Don't interpolate, as points precisely on the interfaces get evaluated and
+ // we don't end up with a jump in the values across it.
+ AffineConstraints<double> constraints;
+ constraints.close();
+ VectorTools::project(
+ dof_handler,
+ constraints,
+ hp::QCollection<dim>(QGauss<spacedim>(fe.max_degree() + 2),
+ QGauss<spacedim>(fe.max_degree() + 3)),
+ DiscontinuousFunction<spacedim>(fe.n_components()),
+ solution);
+
+ const UpdateFlags update_flags =
+ update_values | update_gradients | update_hessians | update_3rd_derivatives;
+ MeshWorker::ScratchData<dim, spacedim> scratch_data(
+ fe, qf_cell, update_flags, qf_face, update_flags);
+
+ const auto cell = dof_handler.begin_active();
+ unsigned int face = 0;
+ while (cell->face(face)->at_boundary())
+ ++face;
+
+ scratch_data.reinit(cell,
+ face,
+ cell->neighbor(face),
+ cell->neighbor_face_no(face));
+ scratch_data.extract_local_dof_values("solution", solution);
+
+ deallog << "Jumps in values: "
+ << scratch_data.get_jumps_in_values("solution", extractor)[0]
+ << std::endl;
+ deallog << "Jumps in gradients: "
+ << scratch_data.get_jumps_in_gradients("solution", extractor)[0]
+ << std::endl;
+ deallog << "Jumps in Hessians: "
+ << scratch_data.get_jumps_in_hessians("solution", extractor)[0]
+ << std::endl;
+ deallog << "Jumps in third derivatives: "
+ << scratch_data.get_jumps_in_third_derivatives("solution",
+ extractor)[0]
+ << std::endl;
+
+ deallog << "Averages of values: "
+ << scratch_data.get_averages_of_values("solution", extractor)[0]
+ << std::endl;
+ deallog << "Averages of gradients: "
+ << scratch_data.get_averages_of_gradients("solution", extractor)[0]
+ << std::endl;
+ deallog << "Averages of Hessians: "
+ << scratch_data.get_averages_of_hessians("solution", extractor)[0]
+ << std::endl;
+
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main(int argc, char *argv[])
+{
+ initlog();
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, testing_max_num_threads());
+
+ FEValuesExtractors::Vector extractor(0);
+
+ run<2>(extractor);
+ run<3>(extractor);
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL:Dim 2::Jumps in values: 4.00030 4.00030
+DEAL:Dim 2::Jumps in gradients: 2.40085 0.000407768 2.40085 0.000407768
+DEAL:Dim 2::Jumps in Hessians: 0.0456300 3.24213 3.24213 -0.000332829 0.0456300 3.24213 3.24213 -0.000332829
+DEAL:Dim 2::Jumps in third derivatives: -2.06802 0.0616191 0.0616191 -2.64630 0.0616191 -2.64630 -2.64630 -0.000238992 -2.06802 0.0616191 0.0616191 -2.64630 0.0616191 -2.64630 -2.64630 -0.000238992
+DEAL:Dim 2::Averages of values: 1.00030 1.00030
+DEAL:Dim 2::Averages of gradients: 0.600214 0.000407768 0.600214 0.000407768
+DEAL:Dim 2::Averages of Hessians: 0.0456300 0.810533 0.810533 -0.000332829 0.0456300 0.810533 0.810533 -0.000332829
+DEAL:Dim 2::OK
+DEAL:Dim 3::Jumps in values: 4.00030 4.00030 4.00030
+DEAL:Dim 3::Jumps in gradients: 2.40085 0.000407768 -9.48894e-15 2.40085 0.000407768 -9.48894e-15 2.40085 0.000407768 -9.48894e-15
+DEAL:Dim 3::Jumps in Hessians: 0.0456300 3.24213 2.80921e-13 3.24213 -0.000332829 -6.32244e-13 2.80921e-13 -6.32244e-13 -2.22489e-13 0.0456300 3.24213 2.80921e-13 3.24213 -0.000332829 -6.32244e-13 2.80921e-13 -6.32244e-13 -2.22489e-13 0.0456300 3.24213 2.80921e-13 3.24213 -0.000332829 -6.32244e-13 2.80921e-13 -6.32244e-13 -2.22489e-13
+DEAL:Dim 3::Jumps in third derivatives: -2.06802 0.0616191 -3.77309e-13 0.0616191 -2.64630 -7.53642e-12 -3.77309e-13 -7.53642e-12 -3.14782e-12 0.0616191 -2.64630 -7.53642e-12 -2.64630 -0.000238992 4.33609e-12 -7.53642e-12 4.33609e-12 4.07230e-12 -3.77309e-13 -7.53642e-12 -3.14782e-12 -7.53642e-12 4.33609e-12 4.07230e-12 -3.14782e-12 4.07230e-12 7.48956e-13 -2.06802 0.0616191 -3.77309e-13 0.0616191 -2.64630 -7.53642e-12 -3.77309e-13 -7.53642e-12 -3.14782e-12 0.0616191 -2.64630 -7.53642e-12 -2.64630 -0.000238992 4.33609e-12 -7.53642e-12 4.33609e-12 4.07230e-12 -3.77309e-13 -7.53642e-12 -3.14782e-12 -7.53642e-12 4.33609e-12 4.07230e-12 -3.14782e-12 4.07230e-12 7.48956e-13 -2.06802 0.0616191 -3.77309e-13 0.0616191 -2.64630 -7.53642e-12 -3.77309e-13 -7.53642e-12 -3.14782e-12 0.0616191 -2.64630 -7.53642e-12 -2.64630 -0.000238992 4.33609e-12 -7.53642e-12 4.33609e-12 4.07230e-12 -3.77309e-13 -7.53642e-12 -3.14782e-12 -7.53642e-12 4.33609e-12 4.07230e-12 -3.14782e-12 4.07230e-12 7.48956e-13
+DEAL:Dim 3::Averages of values: 1.00030 1.00030 1.00030
+DEAL:Dim 3::Averages of gradients: 0.600214 0.000407768 4.74447e-15 0.600214 0.000407768 4.74447e-15 0.600214 0.000407768 4.74447e-15
+DEAL:Dim 3::Averages of Hessians: 0.0456300 0.810533 3.71751e-14 0.810533 -0.000332829 -4.44159e-13 3.71751e-14 -4.44159e-13 -1.74749e-13 0.0456300 0.810533 3.71751e-14 0.810533 -0.000332829 -4.44159e-13 3.71751e-14 -4.44159e-13 -1.74749e-13 0.0456300 0.810533 3.71751e-14 0.810533 -0.000332829 -4.44159e-13 3.71751e-14 -4.44159e-13 -1.74749e-13
+DEAL:Dim 3::OK
+DEAL::OK
projects / dealii.git / commitdiff