From: David Wells Date: Sun, 8 Jan 2017 01:26:59 +0000 (-0500) Subject: Do not reuse CellSimilarity for higher order MappingQGeneric maps. X-Git-Tag: v8.5.0-rc1~285^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F3738%2Fhead;p=dealii.git Do not reuse CellSimilarity for higher order MappingQGeneric maps. This commit fixes bug where two cells, both alike in vertices (but varying in curvature) would sometimes be treated as translations of each-other by MappingQGeneric. This lead to an artificial limit of second order convergence since Jacobian weights (and other values) were not correctly updated. --- diff --git a/doc/news/changes/minor/20170105DavidWells b/doc/news/changes/minor/20170105DavidWells new file mode 100644 index 0000000000..297cc39073 --- /dev/null +++ b/doc/news/changes/minor/20170105DavidWells @@ -0,0 +1,3 @@ +Fixed: MappingQGeneric, when computing values on a given cell, sometimes reused information computed on a different cell whose faces were curved in a different way. MappingQGeneric now recomputes all values when its order is greater than one (that is, when it can represent curved cells). +
+(David Wells, 2017/01/05) diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index ae48043d15..333dae05e7 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -2589,10 +2589,17 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, 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 (QProjector::DataSetDescriptor::cell (), data, output_data.quadrature_points); - internal::maybe_update_Jacobians (cell_similarity, + internal::maybe_update_Jacobians (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), data); @@ -2612,7 +2619,7 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, 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::cell_iterator &cell, 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) @@ -2693,7 +2700,7 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, 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::cell_iterator &cell, 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 (cell_similarity, + internal::maybe_update_jacobian_grads (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), data, output_data.jacobian_grads); - internal::maybe_update_jacobian_pushed_forward_grads (cell_similarity, + internal::maybe_update_jacobian_pushed_forward_grads (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), data, output_data.jacobian_pushed_forward_grads); - internal::maybe_update_jacobian_2nd_derivatives (cell_similarity, + internal::maybe_update_jacobian_2nd_derivatives (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), data, output_data.jacobian_2nd_derivatives); - internal::maybe_update_jacobian_pushed_forward_2nd_derivatives (cell_similarity, + internal::maybe_update_jacobian_pushed_forward_2nd_derivatives (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), data, output_data.jacobian_pushed_forward_2nd_derivatives); - internal::maybe_update_jacobian_3rd_derivatives (cell_similarity, + internal::maybe_update_jacobian_3rd_derivatives (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), data, output_data.jacobian_3rd_derivatives); - internal::maybe_update_jacobian_pushed_forward_3rd_derivatives (cell_similarity, + internal::maybe_update_jacobian_pushed_forward_3rd_derivatives (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), data, output_data.jacobian_pushed_forward_3rd_derivatives); - return cell_similarity; + return computed_cell_similarity; } diff --git a/tests/fe/fe_values_function_manifold.cc b/tests/fe/fe_values_function_manifold.cc new file mode 100644 index 0000000000..f277963813 --- /dev/null +++ b/tests/fe/fe_values_function_manifold.cc @@ -0,0 +1,376 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#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 +class HardManufacturedSolution : public Function +{ +public: + virtual double value(const Point &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 +class HardManufacturedForcing : public Function +{ +public: + virtual double value(const Point &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 +class PushForward : public Function +{ +public: + PushForward() : Function(dim) + {} + + double value(const Point &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::quiet_NaN(); + } +}; + +template +class PullBack : public Function +{ +public: + PullBack() : Function(dim) + {} + + double value(const Point &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::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 +struct CubicRoofFunctions +{ + PushForward forward; + PullBack backward; +}; + +/** + * A manifold describing the cubic roof. + */ +template +class CubicRoofManifold : private CubicRoofFunctions, + public FunctionManifold +{ +public: + CubicRoofManifold() : FunctionManifold(this->forward, this->backward) + {} +}; + +template +std_cxx11::shared_ptr > +cubic_roof(Triangulation &triangulation) +{ + std_cxx11::shared_ptr > boundary(new CubicRoofManifold()); + 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 +class JxWError +{ +public: + JxWError(const unsigned int n_global_refines); + + double run(); + +protected: + std_cxx11::shared_ptr > manufactured_solution; + std_cxx11::shared_ptr > manufactured_forcing; + + std_cxx11::shared_ptr > boundary_manifold; + Triangulation triangulation; + FE_Q finite_element; + DoFHandler dof_handler; + QGauss cell_quadrature; + MappingQGeneric cell_mapping; + + ConstraintMatrix all_constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector system_rhs; + Vector solution; + + void setup_dofs(); + void setup_matrices(); + double solve(); +}; + + + +template +JxWError::JxWError(const unsigned int n_global_refines) : + manufactured_solution(new HardManufacturedSolution()), + manufactured_forcing(new HardManufacturedForcing()), + 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 +void JxWError::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 +void JxWError::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 fe_values(cell_mapping, finite_element, cell_quadrature, flags); + + const unsigned int dofs_per_cell = finite_element.dofs_per_cell; + FullMatrix cell_system(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + + typename DoFHandler::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 +double JxWError::solve() +{ + { + SolverControl solver_control(std::max(types::global_dof_index(100), + static_cast(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 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(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 +double JxWError::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 solver(n_global_refines); + deallog << "L2 error: " + << solver.run() + << std::endl; + } +} diff --git a/tests/fe/fe_values_function_manifold.output b/tests/fe/fe_values_function_manifold.output new file mode 100644 index 0000000000..7e92c1be5e --- /dev/null +++ b/tests/fe/fe_values_function_manifold.output @@ -0,0 +1,4 @@ + +DEAL::L2 error: 0.1233048027 +DEAL::L2 error: 0.004171790017 +DEAL::L2 error: 0.0001334762670