*/
value_type integrate_value () const;
+ /**
+ * Return the determinant of the Jacobian from the unit to the real cell
+ * times the quadrature weight.
+ */
+ VectorizedArray<Number> JxW(const unsigned int q_point) const;
+
//@}
/**
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<n_components_; ++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)
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<n_components_; ++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;
}
}
::fill_JxW_values(AlignedVector<VectorizedArray<Number> > &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)
{
+template <int dim, int n_components_, typename Number>
+inline
+VectorizedArray<Number>
+FEEvaluationBase<dim,n_components_,Number>::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 ().
{
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include "create_mesh.h"
+
+std::ofstream logfile("output");
+
+
+template <int dim>
+void test ()
+{
+ Triangulation<dim> 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<dim>::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<dim> fe (1);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ constraints.close();
+
+ MatrixFree<dim> mf_data;
+ {
+ const QGauss<1> quad (2);
+ typename MatrixFree<dim>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim>::AdditionalData::none;
+ data.mapping_update_flags = update_JxW_values;
+ mf_data.reinit (dof, constraints, quad, data);
+ }
+
+ double error = 0, error2 = 0, abs = 0;
+
+ QGauss<dim> quad(2);
+ FEValues<dim> fe_values(fe, quad, update_JxW_values);
+ FEEvaluation<dim,1> fe_eval(mf_data);
+ AlignedVector<VectorizedArray<double> > jxw_values_manual(fe_eval.n_q_points);
+ for (unsigned int cell=0; cell<mf_data.n_macro_cells(); ++cell)
+ {
+ fe_eval.reinit(cell);
+ fe_eval.fill_JxW_values(jxw_values_manual);
+ for (unsigned int v=0; v<mf_data.n_components_filled(cell); ++v)
+ {
+ fe_values.reinit(mf_data.get_cell_iterator(cell,v));
+ for (unsigned int q=0; q<quad.size(); ++q)
+ {
+ abs += fe_values.JxW(q);
+ error += std::abs(fe_values.JxW(q) - fe_eval.JxW(q)[v]);
+ error2 += std::abs(fe_values.JxW(q) - jxw_values_manual[q][v]);
+ }
+ }
+ }
+
+ deallog << "Norm of difference: " << error/abs << " " << error2/abs
+ << std::endl << std::endl;
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog << std::setprecision (3);
+
+ deallog.threshold_double(1.e-12);
+ test<2>();
+ test<3>();
+}
projects / dealii.git / commitdiff