data.mapping_support_points = this->compute_mapping_support_points(cell);
data.cell_of_current_support_points = cell;
+ // if the order of the mapping is greater than 1, then do not reuse any cell
+ // similarity information. This is necessary because the cell similarity
+ // value is computed with just cell vertices and does not take into account
+ // cell curvature.
+ const CellSimilarity::Similarity computed_cell_similarity =
+ (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
+
internal::maybe_compute_q_points<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.quadrature_points);
- internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
+ internal::maybe_update_Jacobians<dim,spacedim> (computed_cell_similarity,
QProjector<dim>::DataSetDescriptor::cell (),
data);
ExcDimensionMismatch(output_data.normal_vectors.size(), n_q_points));
- if (cell_similarity != CellSimilarity::translation)
+ if (computed_cell_similarity != CellSimilarity::translation)
for (unsigned int point=0; point<n_q_points; ++point)
{
output_data.JxW_values[point]
= sqrt(determinant(G)) * weights[point];
- if (cell_similarity == CellSimilarity::inverted_translation)
+ if (computed_cell_similarity == CellSimilarity::inverted_translation)
{
// we only need to flip the normal
if (update_flags & update_normal_vectors)
if (update_flags & update_jacobians)
{
AssertDimension (output_data.jacobians.size(), n_q_points);
- if (cell_similarity != CellSimilarity::translation)
+ if (computed_cell_similarity != CellSimilarity::translation)
for (unsigned int point=0; point<n_q_points; ++point)
output_data.jacobians[point] = data.contravariant[point];
}
if (update_flags & update_inverse_jacobians)
{
AssertDimension (output_data.inverse_jacobians.size(), n_q_points);
- if (cell_similarity != CellSimilarity::translation)
+ if (computed_cell_similarity != CellSimilarity::translation)
for (unsigned int point=0; point<n_q_points; ++point)
output_data.inverse_jacobians[point] = data.covariant[point].transpose();
}
- internal::maybe_update_jacobian_grads<dim,spacedim> (cell_similarity,
+ internal::maybe_update_jacobian_grads<dim,spacedim> (computed_cell_similarity,
QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.jacobian_grads);
- internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (cell_similarity,
+ internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (computed_cell_similarity,
QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.jacobian_pushed_forward_grads);
- internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (cell_similarity,
+ internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (computed_cell_similarity,
QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.jacobian_2nd_derivatives);
- internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (cell_similarity,
+ internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (computed_cell_similarity,
QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.jacobian_pushed_forward_2nd_derivatives);
- internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (cell_similarity,
+ internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (computed_cell_similarity,
QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.jacobian_3rd_derivatives);
- internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (cell_similarity,
+ internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (computed_cell_similarity,
QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.jacobian_pushed_forward_3rd_derivatives);
- return cell_similarity;
+ return computed_cell_similarity;
}
--- /dev/null
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx11/shared_ptr.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "../tests.h"
+
+// Test that, for a single thread, CellSimilarity calculations return 'none'
+// when we are using curved cells. This is a regression test for a bug where,
+// for this particular manifold, some cells would be marked as translations
+// even though they were not (they varied in curvature of faces). This (small)
+// error causes the Laplace solver to converge at order 2 instead of order 5,
+// resulting in much higher L2 errors when the grid is refined.
+
+using namespace dealii;
+
+static const unsigned int fe_order = 4;
+static const dealii::types::boundary_id boundary_id = 0;
+static const dealii::types::manifold_id cubic_manifold_id = 1;
+static const double pi = M_PI;
+
+// ----------------------------------------------------------------------------
+// Manufactured solution and manufactured forcing
+// ----------------------------------------------------------------------------
+
+template <int dim>
+class HardManufacturedSolution : public Function<dim>
+{
+public:
+ virtual double value(const Point<dim> &point,
+ const unsigned int) const
+ {
+ const double &x = point[0];
+ const double &y = point[1];
+
+ return (Utilities::fixed_power<3>(y) + std::exp(-Utilities::fixed_power<2>(y))
+ + std::sin(4.5*Utilities::fixed_power<2>(y))
+ + std::sin(20*y))*(20*std::cos(4*pi*x) + 0.1*std::sin(20*pi*x)
+ - 80*std::sin(6*pi*x));
+ }
+};
+
+
+
+template <int dim>
+class HardManufacturedForcing : public Function<dim>
+{
+public:
+ virtual double value(const Point<dim> &point,
+ const unsigned int) const
+ {
+ const double &x = point[0];
+ const double &y = point[1];
+ return -40.0*(Utilities::fixed_power<3>(y) + std::exp(-Utilities::fixed_power<2>(y)) + std::sin(4.5*Utilities::fixed_power<2>(y)) + std::sin(20*y))*(-8.0*pi*pi*std::cos(4*pi*x) - 1.0*pi*pi*std::sin(20*pi*x) + 72.0*pi*pi*std::sin(6*pi*x)) - (4*Utilities::fixed_power<2>(y)*std::exp(-Utilities::fixed_power<2>(y)) - 81.0*Utilities::fixed_power<2>(y)*std::sin(4.5*Utilities::fixed_power<2>(y)) + 6*y + 9.0*std::cos(4.5*Utilities::fixed_power<2>(y)) - 2*std::exp(-Utilities::fixed_power<2>(y)) - 400*std::sin(20*y))*(20*std::cos(4*pi*x) + 0.1*std::sin(20*pi*x) - 80*std::sin(6*pi*x));
+ }
+};
+
+
+// ----------------------------------------------------------------------------
+// Description of the curved geometry
+// ----------------------------------------------------------------------------
+
+template <int dim>
+class PushForward : public Function<dim>
+{
+public:
+ PushForward() : Function<dim>(dim)
+ {}
+
+ double value(const Point<dim> &point,
+ const unsigned int component = 0) const
+ {
+ switch (component)
+ {
+ case 0:
+ return point[0];
+ case 1:
+ {
+ const double &x = point[0];
+ return point[1] + 0.25*(2*x - 1.0)*(x - 1.0)*x;
+ }
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return std::numeric_limits<double>::quiet_NaN();
+ }
+};
+
+template <int dim>
+class PullBack : public Function<dim>
+{
+public:
+ PullBack() : Function<dim>(dim)
+ {}
+
+ double value(const Point<dim> &point,
+ const unsigned int component = 0) const
+ {
+ switch (component)
+ {
+ case 0:
+ return point[0];
+ case 1:
+ {
+ const double &x = point[0];
+ return point[1] - 0.25*(2*x - 1.0)*(x - 1.0)*x;
+ }
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return std::numeric_limits<double>::quiet_NaN();
+ }
+};
+
+/**
+ * A struct containing the push forward and pull back functions. This is
+ * only used with CubicRoofManifold: the only reason this struct exists is
+ * so that CubicRoofManifold holds instances of these functions (which are
+ * fed into FunctionManifold, so they must exist at that point in the
+ * constructor call) and destroys them after destroying its FunctionManifold
+ * part.
+ */
+template <int dim>
+struct CubicRoofFunctions
+{
+ PushForward<dim> forward;
+ PullBack<dim> backward;
+};
+
+/**
+ * A manifold describing the cubic roof.
+ */
+template <int dim>
+class CubicRoofManifold : private CubicRoofFunctions<dim>,
+ public FunctionManifold<dim>
+{
+public:
+ CubicRoofManifold() : FunctionManifold<dim>(this->forward, this->backward)
+ {}
+};
+
+template <int dim>
+std_cxx11::shared_ptr<Manifold<dim> >
+cubic_roof(Triangulation<dim> &triangulation)
+{
+ std_cxx11::shared_ptr<Manifold<dim> > boundary(new CubicRoofManifold<dim>());
+ GridGenerator::hyper_cube(triangulation);
+
+ triangulation.set_all_manifold_ids(cubic_manifold_id);
+ triangulation.set_manifold(cubic_manifold_id, *boundary);
+ return boundary;
+}
+
+// ----------------------------------------------------------------------------
+// The Laplace solver that has trouble with 1 thread, but works well with 2
+// ----------------------------------------------------------------------------
+template <int dim>
+class JxWError
+{
+public:
+ JxWError(const unsigned int n_global_refines);
+
+ double run();
+
+protected:
+ std_cxx11::shared_ptr<Function<dim> > manufactured_solution;
+ std_cxx11::shared_ptr<Function<dim> > manufactured_forcing;
+
+ std_cxx11::shared_ptr<Manifold<dim> > boundary_manifold;
+ Triangulation<dim> triangulation;
+ FE_Q<dim> finite_element;
+ DoFHandler<dim> dof_handler;
+ QGauss<dim> cell_quadrature;
+ MappingQGeneric<dim> cell_mapping;
+
+ ConstraintMatrix all_constraints;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+ Vector<double> system_rhs;
+ Vector<double> solution;
+
+ void setup_dofs();
+ void setup_matrices();
+ double solve();
+};
+
+
+
+template <int dim>
+JxWError<dim>::JxWError(const unsigned int n_global_refines) :
+ manufactured_solution(new HardManufacturedSolution<dim>()),
+ manufactured_forcing(new HardManufacturedForcing<dim>()),
+ finite_element(fe_order),
+ dof_handler(triangulation),
+ cell_quadrature(fe_order + 1),
+ cell_mapping(fe_order)
+
+{
+ boundary_manifold = cubic_roof(triangulation);
+ triangulation.refine_global(n_global_refines);
+}
+
+
+
+template <int dim>
+void JxWError<dim>::setup_dofs()
+{
+ dof_handler.distribute_dofs(finite_element);
+ VectorTools::interpolate_boundary_values(cell_mapping,
+ dof_handler,
+ boundary_id,
+ *manufactured_solution,
+ all_constraints);
+ all_constraints.close();
+
+ DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
+ dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler,
+ dynamic_sparsity_pattern,
+ all_constraints,
+ /*keep_constrained_dofs=*/false);
+ sparsity_pattern.copy_from(dynamic_sparsity_pattern);
+}
+
+
+template <int dim>
+void JxWError<dim>::setup_matrices()
+{
+ system_matrix.reinit(sparsity_pattern);
+ system_rhs.reinit(dof_handler.n_dofs());
+
+ const UpdateFlags flags = update_values | update_gradients | update_JxW_values
+ | update_quadrature_points;
+ FEValues<dim> fe_values(cell_mapping, finite_element, cell_quadrature, flags);
+
+ const unsigned int dofs_per_cell = finite_element.dofs_per_cell;
+ FullMatrix<double> cell_system(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell != endc; ++cell)
+ {
+ cell->get_dof_indices(local_dof_indices);
+ cell_system = 0.0;
+ cell_rhs = 0.0;
+ fe_values.reinit(cell);
+
+ for (unsigned int q_point_n = 0;
+ q_point_n < fe_values.n_quadrature_points; ++q_point_n)
+ {
+ const double point_forcing =
+ manufactured_forcing->value(fe_values.quadrature_point(q_point_n));
+
+ for (unsigned int test_n = 0; test_n < dofs_per_cell; ++test_n)
+ {
+ for (unsigned int trial_n = 0; trial_n < dofs_per_cell; ++trial_n)
+ {
+ cell_system(test_n, trial_n) +=
+ fe_values.JxW(q_point_n)*
+ (fe_values.shape_grad(test_n, q_point_n)
+ *fe_values.shape_grad(trial_n, q_point_n));
+ }
+
+ cell_rhs[test_n] += fe_values.JxW(q_point_n)
+ *fe_values.shape_value(test_n, q_point_n)
+ *point_forcing;
+ }
+ }
+
+ all_constraints.distribute_local_to_global(cell_system, cell_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+ }
+}
+
+
+
+template <int dim>
+double JxWError<dim>::solve()
+{
+ {
+ SolverControl solver_control(std::max(types::global_dof_index(100),
+ static_cast<types::global_dof_index>(system_rhs.size())),
+ 1e-14*system_rhs.l2_norm(),
+ false,
+ false);
+ SolverCG<> solver(solver_control);
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.2);
+
+ solution.reinit(system_rhs);
+ solver.solve(system_matrix, solution, system_rhs, preconditioner);
+ all_constraints.distribute(solution);
+ }
+
+ Vector<double> cell_l2_error(triangulation.n_cells());
+ VectorTools::integrate_difference(cell_mapping,
+ dof_handler,
+ solution,
+ *manufactured_solution,
+ cell_l2_error,
+ // use QIterated to avoid spurious superconvergence
+ QIterated<dim>(QGauss<1>(finite_element.degree), 2),
+ VectorTools::L2_norm);
+
+ const double l2_error = VectorTools::compute_global_error
+ (triangulation, cell_l2_error, VectorTools::L2_norm);
+ return l2_error;
+}
+
+
+
+template <int dim>
+double JxWError<dim>::run()
+{
+ setup_dofs();
+ setup_matrices();
+ const double error = solve();
+ return error;
+}
+
+
+
+int main(int argc, char **argv)
+{
+ // Use exactly one thread so that the CellSimilarity checks are not disabled.
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ static const int dim = 2;
+
+ std::ofstream logfile ("output");
+ deallog << std::setprecision(10);
+ deallog.attach(logfile);
+ for (unsigned int n_global_refines = 3; n_global_refines < 6; ++n_global_refines)
+ {
+ JxWError<dim> solver(n_global_refines);
+ deallog << "L2 error: "
+ << solver.run()
+ << std::endl;
+ }
+}
projects / dealii.git / commitdiff