From de321ee3897c798121a06f67f9fc043dfb0cf248 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 14 Oct 2016 20:38:16 +0200 Subject: [PATCH] Add JxW() method to FEEvaluation. --- include/deal.II/matrix_free/fe_evaluation.h | 56 +++++++--- tests/matrix_free/jxw_values.cc | 108 ++++++++++++++++++++ tests/matrix_free/jxw_values.output | 5 + 3 files changed, 156 insertions(+), 13 deletions(-) create mode 100644 tests/matrix_free/jxw_values.cc create mode 100644 tests/matrix_free/jxw_values.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index f7f5302966..c72adc88cb 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -464,6 +464,12 @@ public: */ value_type integrate_value () const; + /** + * Return the determinant of the Jacobian from the unit to the real cell + * times the quadrature weight. + */ + VectorizedArray JxW(const unsigned int q_point) const; + //@} /** @@ -2022,7 +2028,7 @@ FEEvaluationBase jacobian_grad_upper(0), cell (numbers::invalid_unsigned_int), cell_type (internal::MatrixFreeFunctions::undefined), - cell_data_number (0), + cell_data_number (numbers::invalid_unsigned_int), first_selected_component (0) { for (unsigned int c=0; c jacobian_grad_upper(0), cell (0), cell_type (internal::MatrixFreeFunctions::general), - cell_data_number (0), + cell_data_number (numbers::invalid_unsigned_int), // keep the number of the selected component within the current base element // for reading dof values first_selected_component (fe.component_to_base_index(first_selected_component).second) @@ -2137,16 +2143,20 @@ FEEvaluationBase mapping_info (other.mapping_info), stored_shape_info (other.stored_shape_info), data (other.data), - cartesian_data (other.cartesian_data), - jacobian (other.jacobian), - J_value (other.J_value), - quadrature_weights (other.quadrature_weights), - quadrature_points (other.quadrature_points), - jacobian_grad (other.jacobian_grad), - jacobian_grad_upper(other.jacobian_grad_upper), - cell (other.cell), - cell_type (other.cell_type), - cell_data_number (other.cell_data_number), + cartesian_data (0), + jacobian (0), + J_value (0), + quadrature_weights (mapping_info != 0 ? + mapping_info->mapping_data_gen[quad_no]. + quadrature_weights[active_quad_index].begin() + : + 0), + quadrature_points (0), + jacobian_grad (0), + jacobian_grad_upper(0), + cell (numbers::invalid_unsigned_int), + cell_type (internal::MatrixFreeFunctions::general), + cell_data_number (numbers::invalid_unsigned_int), first_selected_component (other.first_selected_component) { for (unsigned int c=0; c jacobian = mapped_geometry->get_inverse_jacobians().begin(); J_value = mapped_geometry->get_JxW_values().begin(); quadrature_points = mapped_geometry->get_quadrature_points().begin(); + cell = 0; } } @@ -2336,7 +2347,7 @@ FEEvaluationBase ::fill_JxW_values(AlignedVector > &JxW_values) const { AssertDimension(JxW_values.size(), data->n_q_points); - Assert (this->J_value != 0, ExcNotImplemented()); + Assert (this->J_value != 0, ExcNotInitialized()); if (this->cell_type == internal::MatrixFreeFunctions::cartesian || this->cell_type == internal::MatrixFreeFunctions::affine) { @@ -2352,6 +2363,24 @@ FEEvaluationBase +template +inline +VectorizedArray +FEEvaluationBase::JxW(const unsigned int q_point) const +{ + Assert (this->J_value != 0, ExcNotInitialized()); + if (this->cell_type == internal::MatrixFreeFunctions::cartesian || + this->cell_type == internal::MatrixFreeFunctions::affine) + { + Assert (this->mapping_info != 0, ExcInternalError()); + return this->J_value[0] * this->quadrature_weights[q_point]; + } + else + return this->J_value[q_point]; +} + + + namespace internal { // write access to generic vectors that have operator (). @@ -6928,6 +6957,7 @@ FEEvaluation { Assert (this->mapping_info->quadrature_points_initialized == true, ExcNotInitialized()); + Assert (this->quadrature_points != 0, ExcNotInitialized()); AssertIndexRange (q, n_q_points); // Cartesian mesh: not all quadrature points are stored, only the diff --git a/tests/matrix_free/jxw_values.cc b/tests/matrix_free/jxw_values.cc new file mode 100644 index 0000000000..b4cec1ff97 --- /dev/null +++ b/tests/matrix_free/jxw_values.cc @@ -0,0 +1,108 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// this function tests the correctness of JxW values returned by FEEvaluation +// when compared to FEValues + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + +#include "create_mesh.h" + +std::ofstream logfile("output"); + + +template +void test () +{ + Triangulation tria; + create_mesh (tria); + tria.refine_global(4-dim); + + // refine a few cells + for (unsigned int i=0; i<10-3*dim; ++i) + { + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end(); + unsigned int counter = 0; + for (; cell!=endc; ++cell, ++counter) + if (counter % (7-i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_Q fe (1); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + constraints.close(); + + MatrixFree mf_data; + { + const QGauss<1> quad (2); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_JxW_values; + mf_data.reinit (dof, constraints, quad, data); + } + + double error = 0, error2 = 0, abs = 0; + + QGauss quad(2); + FEValues fe_values(fe, quad, update_JxW_values); + FEEvaluation fe_eval(mf_data); + AlignedVector > jxw_values_manual(fe_eval.n_q_points); + for (unsigned int cell=0; cell(); + test<3>(); +} diff --git a/tests/matrix_free/jxw_values.output b/tests/matrix_free/jxw_values.output new file mode 100644 index 0000000000..eb1bf849db --- /dev/null +++ b/tests/matrix_free/jxw_values.output @@ -0,0 +1,5 @@ + +DEAL::Norm of difference: 0 0 +DEAL:: +DEAL::Norm of difference: 0 0 +DEAL:: -- 2.39.5