#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q_generic.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
+template <int dim>
+class DeformedManifold : public ChartManifold<dim, dim, dim>
+{
+public:
+ DeformedManifold() = default;
+
+ virtual std::unique_ptr<Manifold<dim, dim>>
+ clone() const
+ {
+ return std_cxx14::make_unique<DeformedManifold<dim>>();
+ }
+
+ virtual Point<dim>
+ push_forward(const Point<dim> &chart_point) const
+ {
+ Point<dim> result = chart_point;
+ result[0] = std::tanh(2.0 * result[0]) / std::tanh(2.0);
+ return result;
+ }
+
+ virtual Point<dim>
+ pull_back(const Point<dim> &space_point) const
+ {
+ Point<dim> result = space_point;
+ result[0] = 0.5 * std::atanh(result[0] * std::tanh(2.0));
+ return result;
+ }
+};
+
+
+
+template <int dim>
+void
+test_deformed_cube()
+{
+ deallog << "Deformed hyper cube" << std::endl;
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria, -1., 1.);
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, DeformedManifold<dim>());
+ tria.refine_global(6 - dim);
+
+ FE_Q<dim> fe(2);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+ constraints.close();
+
+ const QGauss<1> quad(3);
+ MatrixFree<dim> mf;
+ typename MatrixFree<dim>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim>::AdditionalData::none;
+ data.mapping_update_flags_inner_faces =
+ update_gradients | update_normal_vectors;
+ data.mapping_update_flags_boundary_faces =
+ update_gradients | update_normal_vectors;
+
+ mf.reinit(dof, constraints, quad, data);
+ const unsigned int n_macro_cells = mf.n_macro_cells();
+ Assert(n_macro_cells > 1, ExcInternalError());
+
+ {
+ std::vector<unsigned int> n_cell_types(4, 0);
+ for (unsigned int i = 0; i < n_macro_cells; ++i)
+ n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
+
+ // should have all Cartesian type and no other cell type
+ AssertDimension(n_cell_types[0], n_macro_cells);
+
+ // should have as many different Jacobians as we have cell batches in x
+ // direction; as the mesh is a cube, we can easily calculate it
+ deallog << "Number of different Jacobians: "
+ << mf.get_mapping_info().cell_data[0].jacobians[0].size() / 2
+ << std::endl;
+ }
+
+ // check again, now using a mapping that displaces points
+ {
+ MappingQGeneric<dim> mapping(3);
+ mf.reinit(mapping, dof, constraints, quad, data);
+
+ std::vector<unsigned int> n_cell_types(4, 0);
+ for (unsigned int i = 0; i < n_macro_cells; ++i)
+ n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
+
+ // should have all general type and no other cell type
+ AssertDimension(n_cell_types[3], n_macro_cells);
+
+ // should have as many different Jacobians as we have cell batches in x
+ // direction times the number of quadrature points; as the mesh is a cube,
+ // we can easily calculate it
+ deallog << "Number of different Jacobians: "
+ << mf.get_mapping_info().cell_data[0].jacobians[0].size()
+ << std::endl;
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
int
main()
{
test<2>();
test_cube<2>();
test_parallelogram<2>();
+ test_deformed_cube<2>();
deallog.pop();
deallog.push("3d");
test<3>();
test_cube<3>();
test_parallelogram<3>();
+ test_deformed_cube<3>();
deallog.pop();
}
}
DEAL:2d::OK
DEAL:2d::Parallelogram
DEAL:2d::OK
+DEAL:2d::Deformed hyper cube
+DEAL:2d::Number of different Jacobians: 8
+DEAL:2d::Number of different Jacobians: 72
+DEAL:2d::OK
DEAL:3d::General mesh
DEAL:3d::OK
DEAL:3d::Hyper cube
DEAL:3d::OK
DEAL:3d::Parallelogram
DEAL:3d::OK
+DEAL:3d::Deformed hyper cube
+DEAL:3d::Number of different Jacobians: 4
+DEAL:3d::Number of different Jacobians: 108
+DEAL:3d::OK
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Let MatrixFree initialize the geometry for ridiculously high numbers of
+// quadrature points and compute the volume and surface of a hyper ball. In an
+// initial implementation we used to consume a lot of memory and take a long
+// time, whereas the new state is relatively cheap.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball(tria, Point<dim>(), 1.);
+
+ FE_DGQ<dim> fe(1);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ MappingQGeneric<dim> mapping(degree);
+
+ MatrixFree<dim, double> mf_data;
+ const QGauss<1> quad(degree);
+ typename MatrixFree<dim, double>::AdditionalData data;
+ data.mapping_update_flags_boundary_faces =
+ (update_gradients | update_JxW_values);
+
+ mf_data.reinit(
+ mapping,
+ std::vector<const DoFHandler<dim> *>({{&dof}}),
+ std::vector<const AffineConstraints<double> *>({{&constraints}}),
+ std::vector<Quadrature<1>>({{QGauss<1>(std::max(degree / 2, 1U)),
+ QGauss<1>(degree + 1),
+ QGauss<1>(3 * degree / 2)}}),
+ data);
+
+ deallog << std::setw(5) << degree;
+ for (unsigned int index = 0; index < 3; ++index)
+ {
+ double volume = 0;
+ FEEvaluation<dim, -1> fe_eval(mf_data, 0, index);
+ for (unsigned int cell = 0; cell < mf_data.n_cell_batches(); ++cell)
+ {
+ fe_eval.reinit(cell);
+ VectorizedArray<double> local = 0;
+ for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+ local += fe_eval.JxW(q);
+ for (unsigned int v = 0;
+ v < mf_data.n_active_entries_per_cell_batch(cell);
+ ++v)
+ volume += local[v];
+ }
+ deallog << std::setw(11)
+ << std::abs((dim == 2 ? numbers::PI : 4. / 3. * numbers::PI) -
+ volume);
+ }
+ deallog << " ";
+ for (unsigned int index = 0; index < 3; ++index)
+ {
+ double area = 0;
+ FEFaceEvaluation<dim, -1> fe_eval(mf_data, true, 0, index);
+ for (unsigned int face = mf_data.n_inner_face_batches();
+ face <
+ mf_data.n_boundary_face_batches() + mf_data.n_inner_face_batches();
+ ++face)
+ {
+ fe_eval.reinit(face);
+ VectorizedArray<double> local = 0;
+ for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+ local += fe_eval.JxW(q);
+ for (unsigned int v = 0;
+ v < mf_data.n_active_entries_per_face_batch(face);
+ ++v)
+ area += local[v];
+ }
+ deallog << std::setw(11) << std::abs(2. * (dim - 1) * numbers::PI - area);
+ }
+ deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ deallog << std::setprecision(4);
+ deallog << "2D volume area" << std::endl;
+ deallog << "deg err deg/2 err deg+1 err 3deg/2 "
+ << "err deg/2 err deg+1 err 3deg/2" << std::endl;
+ for (unsigned int q = 1; q < 20; ++q)
+ test<2>(q);
+ deallog << std::endl;
+
+ deallog << "3D volume area" << std::endl;
+ deallog << "deg err deg/2 err deg+1 err 3deg/2 "
+ << "err deg/2 err deg+1 err 3deg/2" << std::endl;
+ for (unsigned int q = 1; q < 17; ++q)
+ test<3>(q);
+ deallog << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::2D volume area
+DEAL::deg err deg/2 err deg+1 err 3deg/2 err deg/2 err deg+1 err 3deg/2
+DEAL::1 1.142 1.142 1.142 0.6263 0.6263 0.6263
+DEAL::2 0.02998 0.03702 0.03702 0.6263 0.03493 0.03493
+DEAL::3 0.2367 0.0005893 0.0005893 0.003915 0.0001972 0.0001972
+DEAL::4 0.004112 5.461e-06 5.461e-06 0.008015 5.006e-06 4.602e-06
+DEAL::5 0.006630 3.288e-08 3.288e-08 8.248e-05 8.603e-09 1.122e-08
+DEAL::6 3.905e-05 1.386e-10 1.386e-10 4.224e-05 1.254e-10 1.136e-10
+DEAL::7 5.266e-05 4.334e-13 4.370e-13 3.522e-07 9.948e-14 1.368e-13
+DEAL::8 1.580e-07 8.882e-16 3.553e-15 1.179e-07 1.776e-15 2.665e-15
+DEAL::9 1.964e-07 6.661e-15 1.776e-15 8.182e-10 3.553e-15 2.665e-15
+DEAL::10 3.578e-10 3.553e-15 2.665e-15 2.036e-10 0.000 8.882e-16
+DEAL::11 4.246e-10 2.665e-15 3.553e-15 1.210e-12 1.776e-15 1.776e-15
+DEAL::12 5.302e-13 8.438e-15 1.332e-15 2.327e-13 8.882e-16 1.776e-15
+DEAL::13 5.844e-13 4.441e-15 1.776e-15 7.105e-15 1.776e-15 8.882e-16
+DEAL::14 8.882e-16 1.776e-15 1.155e-14 8.882e-16 1.776e-15 8.882e-16
+DEAL::15 7.105e-15 5.329e-15 7.994e-15 6.217e-15 8.882e-15 4.441e-15
+DEAL::16 1.865e-14 1.243e-14 4.885e-15 8.882e-16 0.000 1.776e-15
+DEAL::17 1.243e-14 1.332e-15 1.332e-14 8.882e-16 8.882e-16 1.776e-15
+DEAL::18 4.441e-16 4.885e-15 9.326e-15 2.665e-15 1.776e-15 8.882e-16
+DEAL::19 3.997e-15 5.773e-15 7.994e-15 4.441e-15 0.000 2.665e-15
+DEAL::
+DEAL::3D volume area
+DEAL::deg err deg/2 err deg+1 err 3deg/2 err deg/2 err deg+1 err 3deg/2
+DEAL::1 2.747 2.649 2.747 4.566 4.566 4.566
+DEAL::2 0.1217 0.1120 0.1120 0.5664 0.2217 0.2217
+DEAL::3 0.8702 0.01703 0.01703 3.949 0.02024 0.02024
+DEAL::4 0.03972 0.0005453 0.0005447 0.03028 0.001079 0.001038
+DEAL::5 0.1274 0.0001358 0.0001358 0.3096 5.472e-05 6.758e-05
+DEAL::6 0.002403 5.127e-06 5.121e-06 0.006959 1.007e-05 9.541e-06
+DEAL::7 0.01271 1.081e-06 1.080e-06 0.03144 1.443e-07 2.279e-08
+DEAL::8 0.0001101 4.117e-08 4.111e-08 0.001177 8.077e-08 7.626e-08
+DEAL::9 0.001234 9.165e-09 9.161e-09 0.002962 7.696e-09 5.659e-09
+DEAL::10 7.936e-07 3.567e-10 3.561e-10 0.0001491 6.945e-10 6.492e-10
+DEAL::11 0.0001194 7.771e-11 7.768e-11 0.0002851 1.123e-10 9.018e-11
+DEAL::12 7.006e-07 3.033e-12 3.001e-12 1.740e-05 5.846e-12 5.437e-12
+DEAL::13 1.151e-05 6.777e-13 6.688e-13 2.715e-05 1.450e-12 1.204e-12
+DEAL::14 1.345e-07 3.553e-14 9.770e-15 1.934e-06 4.619e-14 3.908e-14
+DEAL::15 1.112e-06 1.688e-14 4.885e-14 2.602e-06 2.842e-14 3.197e-14
+DEAL::16 1.865e-08 3.464e-14 1.332e-14 2.088e-07 1.421e-14 1.066e-14
+DEAL::
projects / dealii.git / commitdiff