]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add JxW() method to FEEvaluation.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 14 Oct 2016 18:38:16 +0000 (20:38 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 14 Oct 2016 19:17:40 +0000 (21:17 +0200)
include/deal.II/matrix_free/fe_evaluation.h
tests/matrix_free/jxw_values.cc [new file with mode: 0644]
tests/matrix_free/jxw_values.output [new file with mode: 0644]

index f7f5302966728d8cbb284649aa27ca248dcd0921..c72adc88cbad3b0e88b59c658ea2113979a5e796 100644 (file)
@@ -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<Number> JxW(const unsigned int q_point) const;
+
   //@}
 
   /**
@@ -2022,7 +2028,7 @@ FEEvaluationBase<dim,n_components_,Number>
   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)
@@ -2086,7 +2092,7 @@ FEEvaluationBase<dim,n_components_,Number>
   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<dim,n_components_,Number>
   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)
@@ -2170,6 +2180,7 @@ FEEvaluationBase<dim,n_components_,Number>
       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<dim,n_components_,Number>
 ::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)
     {
@@ -2352,6 +2363,24 @@ FEEvaluationBase<dim,n_components_,Number>
 
 
 
+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 ().
@@ -6928,6 +6957,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 {
   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 (file)
index 0000000..b4cec1f
--- /dev/null
@@ -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 <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>();
+}
diff --git a/tests/matrix_free/jxw_values.output b/tests/matrix_free/jxw_values.output
new file mode 100644 (file)
index 0000000..eb1bf84
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Norm of difference: 0 0
+DEAL::
+DEAL::Norm of difference: 0 0
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.