#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/parallel.h>
-#include <deal.II/base/vectorization.h>
#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
* Sets the vector size to the given size and initializes all elements with
* T().
*/
- AlignedVector (const size_type size);
+ AlignedVector (const size_type size,
+ const T &init = T());
/**
* Destructor.
static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
public:
/**
- * Constructor. Issues a parallel call if
- * there are sufficiently many elements,
- * otherwise work in serial. Copies the data
- * from source to destination and then calls
- * destructor on the source. If the optional
- * argument is set to true, the source is left
- * untouched instead.
+ * Constructor. Issues a parallel call if there are sufficiently many
+ * elements, otherwise work in serial. Copies the data from source to
+ * destination and then calls destructor on the source. If the optional
+ * argument is set to true, the source is left untouched instead.
*/
AlignedVectorMove (T *source_begin,
T *source_end,
}
/**
- * This method moves elements from the source
- * to the destination given in the constructor
- * on a subrange given by two integers.
+ * This method moves elements from the source to the destination given in
+ * the constructor on a subrange given by two integers.
*/
virtual void apply_to_subrange (const std::size_t begin,
const std::size_t end) const
{
- // for classes trivial assignment can use
- // memcpy
+ // for classes trivial assignment can use memcpy
if (std_cxx1x::is_trivial<T>::value == true)
std::memcpy (destination_+begin, source_+begin, (end-begin)*sizeof(T));
else if (copy_only_ == false)
for (std::size_t i=begin; i<end; ++i)
{
- // initialize memory, copy, and destruct
- new (&destination_[i]) T;
- destination_[i] = source_[i];
+ // initialize memory (copy construct), and destruct
+ new (&destination_[i]) T(source_[i]);
source_[i].~T();
}
else
for (std::size_t i=begin; i<end; ++i)
{
- new (&destination_[i]) T;
- destination_[i] = source_[i];
+ new (&destination_[i]) T(source_[i]);
}
}
static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1;
public:
/**
- * Constructor. Issues a parallel call if
- * there are sufficiently many elements,
- * otherwise work in serial.
+ * Constructor. Issues a parallel call if there are sufficiently many
+ * elements, otherwise work in serial.
*/
AlignedVectorSet (const std::size_t size,
const T &element,
private:
/**
- * This sets elements on a subrange given by
- * two integers.
+ * This sets elements on a subrange given by two integers.
*/
virtual void apply_to_subrange (const std::size_t begin,
const std::size_t end) const
{
- // for classes with trivial assignment of zero
- // can use memset
+ // for classes with trivial assignment of zero can use memset
if (std_cxx1x::is_trivial<T>::value == true && trivial_element)
std::memset (destination_+begin, 0, (end-begin)*sizeof(T));
else
+ // initialize memory and set
for (std::size_t i=begin; i<end; ++i)
- {
- // initialize memory and set
- new (&destination_[i]) T;
- destination_[i] = element_;
- }
+ new (&destination_[i]) T(element_);
}
const T &element_;
template < class T >
inline
-AlignedVector<T>::AlignedVector (const size_type size)
+AlignedVector<T>::AlignedVector (const size_type size,
+ const T &init)
:
_data (0),
_end_data (0),
_end_allocated (0)
{
if (size > 0)
- resize (size);
+ resize (size, init);
}
#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
- // allocate and align along boundaries of the size of
- // VectorizedArray<double>, which is 16 bytes for SSE and 32 bytes for
- // AVX
- T *new_data = static_cast<T *>(_mm_malloc (size_actual_allocate,
- sizeof(VectorizedArray<double>)));
+ // allocate and align along 64-byte boundaries (this is enough for all
+ // levels of vectorization currently supported by deal.II)
+ T *new_data = static_cast<T *>(_mm_malloc (size_actual_allocate, 64));
#else
T *new_data = static_cast<T *>(malloc (size_actual_allocate));
#endif
* @name 1: General operations
*/
//@{
- /**
- * Destructor.
- */
- ~FEEvaluationBase();
-
/**
* Initializes the operation pointer to the current cell. Unlike the
* FEValues::reinit function, where the information related to a particular
const DoFHandler<dim> &dof_handler,
const unsigned int base_element = 0);
+ /**
+ * Copy constructor
+ */
+ FEEvaluationBase (const FEEvaluationBase &other);
+
/**
* A unified function to read from and write into vectors based on the given
* template operation. It can perform the operation for @p read_dof_values,
* Stores a reference to the unit cell data, i.e., values, gradients and
* Hessians in 1D at the quadrature points that constitute the tensor
* product. Also contained in matrix_info, but it simplifies code if we
- * store a reference to it. If initialized without MatrixFree object, this
- * call actually initializes the evaluation.
+ * store a reference to it. If the object is initialized without MatrixFree
+ * object, the constructor creates this data structure.
*/
- const internal::MatrixFreeFunctions::ShapeInfo<Number> *data;
+ std_cxx1x::shared_ptr<const internal::MatrixFreeFunctions::ShapeInfo<Number> > data;
/**
* A pointer to the Cartesian Jacobian information of the present cell. Only
const unsigned int dofs_per_cell,
const unsigned int n_q_points);
-
/**
* Constructor that comes with reduced functionality and works similar as
* FEValues. The user has to provide a structure of type EvaluatedGeometry
FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
const unsigned int base_element = 0);
+
+ /**
+ * Copy constructor
+ */
+ FEEvaluationAccess (const FEEvaluationAccess &other);
};
FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
const unsigned int base_element = 0);
+
+ /**
+ * Copy constructor
+ */
+ FEEvaluationAccess (const FEEvaluationAccess &other);
};
FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
const DoFHandler<dim> &dof_handler,
const unsigned int base_element = 0);
+
+ /**
+ * Copy constructor
+ */
+ FEEvaluationAccess (const FEEvaluationAccess &other);
};
const DoFHandler<dim> &dof_handler,
const unsigned int base_element = 0);
+ /**
+ * Copy constructor
+ */
+ FEEvaluationGeneral (const FEEvaluationGeneral &other);
+
/**
* Evaluates the function values, the gradients, and the Laplacians of the
* FE function given at the DoF values in the input vector at the quadrature
* Internally stored variables for the different data fields.
*/
VectorizedArray<Number> my_data_array[n_components*(dofs_per_cell+(dim*dim+2*dim+1)*n_q_points)];
+
+private:
+ /**
+ * Sets the pointers from the data array to values_dof, etc.
+ */
+ void set_data_pointers();
};
const unsigned int fe_no = 0,
const unsigned int quad_no = 0);
+ /**
+ * Copy constructor
+ */
+ FEEvaluation (const FEEvaluation &other);
+
/**
* Constructor that comes with reduced functionality and works similar as
* FEValues. The user has to provide a structure of type EvaluatedGeometry
const DoFHandler<dim> &dof_handler,
const unsigned int base_element = 0);
+ /**
+ * Copy constructor
+ */
+ FEEvaluationGL (const FEEvaluationGL &other);
+
/**
* Evaluates the function values, the gradients, and the Hessians of the FE
* function given at the DoF values in the input vector at the quadrature
const DoFHandler<dim> &dof_handler,
const unsigned int base_element = 0);
+ /**
+ * Copy constructor
+ */
+ FEEvaluationDGP (const FEEvaluationDGP &other);
+
/**
* Evaluates the function values, the gradients, and the Hessians of the FE
* function given at the DoF values in the input vector at the quadrature
#ifndef DOXYGEN
+namespace internal
+{
+ namespace MatrixFreeFunctions
+ {
+ // a small class that gives control over the delete behavior of
+ // std::shared_ptr: we need to disable it when we initialize a pointer
+ // from another structure.
+ template <typename CLASS>
+ struct DummyDeleter
+ {
+ DummyDeleter (const bool do_delete = false)
+ :
+ do_delete(do_delete)
+ {}
+
+ void operator () (CLASS *pointer)
+ {
+ if (do_delete)
+ delete pointer;
+ }
+
+ const bool do_delete;
+ };
+ }
+}
+
+
+
/*----------------------- FEEvaluationBase ----------------------------------*/
template <int dim, int n_components_, typename Number>
mapping_info (&data_in.get_mapping_info()),
data (&data_in.get_shape_info
(fe_no_in, quad_no_in, active_fe_index,
- active_quad_index)),
+ active_quad_index),
+ internal::MatrixFreeFunctions::DummyDeleter
+ <const internal::MatrixFreeFunctions::ShapeInfo<Number> >(false)),
cartesian_data (0),
jacobian (0),
J_value (0),
template <int dim, int n_components_, typename Number>
inline
-FEEvaluationBase<dim,n_components_,Number>::~FEEvaluationBase()
+FEEvaluationBase<dim,n_components_,Number>
+::FEEvaluationBase (const FEEvaluationBase<dim,n_components_,Number> &other)
+ :
+ quad_no (other.quad_no),
+ n_fe_components (other.n_fe_components),
+ active_fe_index (other.active_fe_index),
+ active_quad_index (other.active_quad_index),
+ matrix_info (other.matrix_info),
+ dof_info (other.dof_info),
+ mapping_info (other.mapping_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),
+ evaluated_geometry (other.evaluated_geometry),
+ dof_handler (other.dof_handler)
{
- // delete memory held by data in case this structure was initialized without
- // a MatrixFree object which could hold the data structure.
- if (matrix_info == 0)
- delete data;
+ for (unsigned int c=0; c<n_components_; ++c)
+ {
+ values_dofs[c] = 0;
+ values_quad[c] = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ gradients_quad[c][d] = 0;
+ for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
+ hessians_quad[c][d] = 0;
+ }
}
+
template <int dim, int n_components_, typename Number>
inline
void
+template <int dim, int n_components_, typename Number>
+inline
+FEEvaluationAccess<dim,n_components_,Number>
+::FEEvaluationAccess (const FEEvaluationAccess<dim,n_components_,Number> &other)
+ :
+ FEEvaluationBase <dim,n_components_,Number>(other)
+{}
+
+
+
/*-------------------- FEEvaluationAccess scalar ----------------------------*/
+template <int dim, typename Number>
+inline
+FEEvaluationAccess<dim,1,Number>
+::FEEvaluationAccess (const FEEvaluationAccess<dim,1,Number>&other)
+ :
+ FEEvaluationBase <dim,1,Number>(other)
+{}
+
+
+
template <int dim, typename Number>
inline
VectorizedArray<Number>
+template <int dim, typename Number>
+inline
+FEEvaluationAccess<dim,dim,Number>
+::FEEvaluationAccess (const FEEvaluationAccess<dim,dim,Number>&other)
+ :
+ FEEvaluationBase <dim,dim,Number>(other)
+{}
+
+
+
template <int dim, typename Number>
inline
Tensor<2,dim,VectorizedArray<Number> >
:
BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points)
{
- // set the pointers to the correct position in the data array
- for (unsigned int c=0; c<n_components_; ++c)
- {
- this->values_dofs[c] = &my_data_array[c*dofs_per_cell];
- this->values_quad[c] = &my_data_array[n_components*dofs_per_cell+c*n_q_points];
- for (unsigned int d=0; d<dim; ++d)
- this->gradients_quad[c][d] = &my_data_array[n_components*(dofs_per_cell+n_q_points)
- +
- (c*dim+d)*n_q_points];
- for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
- this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+dofs_per_cell)
- +
- (c*(dim*dim+dim)+d)*n_q_points];
- }
+ set_data_pointers();
+
#ifdef DEBUG
// print error message when the dimensions do not match. Propose a possible
// fix
const unsigned int base_element)
:
BaseClass (geometry, dof_handler, base_element)
+{
+ set_data_pointers();
+}
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationGeneral (const FEEvaluationGeneral &other)
+ :
+ BaseClass (other)
+{
+ set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::set_data_pointers()
{
// set the pointers to the correct position in the data array
for (unsigned int c=0; c<n_components_; ++c)
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const FEEvaluation &other)
+ :
+ BaseClass (other)
+{
+ compute_even_odd_factors();
+}
+
+
+
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
if (fe_degree > 1 || n_q_points_1d > 3)
internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
VectorizedArray<Number>, direction, dof_to_quad, add, 0>
- (shape_val_evenodd, in, out);
+ (&shape_val_evenodd[0], in, out);
else
internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
VectorizedArray<Number>, direction, dof_to_quad, add>
if (fe_degree > 1 || n_q_points_1d > 3)
internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
VectorizedArray<Number>, direction, dof_to_quad, add, 1>
- (shape_gra_evenodd, in, out);
+ (&shape_gra_evenodd[0], in, out);
else
internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
VectorizedArray<Number>, direction, dof_to_quad, add>
if (fe_degree > 1 || n_q_points_1d > 3)
internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
VectorizedArray<Number>, direction, dof_to_quad, add, 2>
- (shape_hes_evenodd, in, out);
+ (&shape_hes_evenodd[0], in, out);
else
internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
VectorizedArray<Number>, direction, dof_to_quad, add>
+template <int dim, int fe_degree, int n_components_, typename Number>
+inline
+FEEvaluationGL<dim,fe_degree,n_components_,Number>
+::FEEvaluationGL (const FEEvaluationGL &other)
+ :
+ BaseClass (other)
+{}
+
+
+
template <int dim, int fe_degree, int n_components_, typename Number>
inline
void
{}
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationDGP (const FEEvaluationDGP &other)
+ :
+ BaseClass (other)
+{}
+
+
+
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
typename Number>
inline
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// same test as matrix_vector_stokes_noflux, but allocating FEEvaluation on
+// the heap (using AlignedVector) instead of allocating it on the stack. Tests
+// also copy constructors of FEEvaluation.
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/aligned_vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+
+
+template <int dim, int degree_p, typename VectorType>
+class MatrixFreeTest
+{
+public:
+ typedef typename DoFHandler<dim>::active_cell_iterator CellIterator;
+ typedef double Number;
+
+ MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+ data (data_in)
+ {};
+
+ void
+ local_apply (const MatrixFree<dim,Number> &data,
+ VectorType &dst,
+ const VectorType &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ typedef VectorizedArray<Number> vector_t;
+ // allocate FEEvaluation. This test will test proper alignment
+ AlignedVector<FEEvaluation<dim,degree_p+1,degree_p+2,dim,Number> > velocity
+ (1, FEEvaluation<dim,degree_p+1,degree_p+2,dim,Number>(data, 0));
+ AlignedVector<FEEvaluation<dim,degree_p, degree_p+2,1, Number> > pressure
+ (1, FEEvaluation<dim,degree_p, degree_p+2,1, Number>(data, 1));
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ velocity[0].reinit (cell);
+ velocity[0].read_dof_values (src.block(0));
+ velocity[0].evaluate (false,true,false);
+ pressure[0].reinit (cell);
+ pressure[0].read_dof_values (src.block(1));
+ pressure[0].evaluate (true,false,false);
+
+ for (unsigned int q=0; q<velocity[0].n_q_points; ++q)
+ {
+ SymmetricTensor<2,dim,vector_t> sym_grad_u =
+ velocity[0].get_symmetric_gradient (q);
+ vector_t pres = pressure[0].get_value(q);
+ vector_t div = -velocity[0].get_divergence(q);
+ pressure[0].submit_value (div, q);
+
+ // subtract p * I
+ for (unsigned int d=0; d<dim; ++d)
+ sym_grad_u[d][d] -= pres;
+
+ velocity[0].submit_symmetric_gradient(sym_grad_u, q);
+ }
+
+ velocity[0].integrate (false,true);
+ velocity[0].distribute_local_to_global (dst.block(0));
+ pressure[0].integrate (true,false);
+ pressure[0].distribute_local_to_global (dst.block(1));
+ }
+ }
+
+
+ void vmult (VectorType &dst,
+ const VectorType &src) const
+ {
+ dst = 0;
+ data.cell_loop (&MatrixFreeTest<dim,degree_p,VectorType>::local_apply,
+ this, dst, src);
+ };
+
+private:
+ const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_shell (triangulation, Point<dim>(),
+ 0.5, 1., 96, true);
+ static HyperShellBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+ triangulation.set_boundary (1, boundary);
+ triangulation.begin_active()->set_refine_flag();
+ triangulation.last()->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+ triangulation.refine_global (3-dim);
+ triangulation.last()->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+
+ MappingQ<dim> mapping (3);
+ FE_Q<dim> fe_u_scal (fe_degree+1);
+ FESystem<dim> fe_u (fe_u_scal,dim);
+ FE_Q<dim> fe_p (fe_degree);
+ FESystem<dim> fe (fe_u_scal, dim, fe_p, 1);
+ DoFHandler<dim> dof_handler_u (triangulation);
+ DoFHandler<dim> dof_handler_p (triangulation);
+ DoFHandler<dim> dof_handler (triangulation);
+
+ MatrixFree<dim,double> mf_data;
+
+ ConstraintMatrix constraints, constraints_u, constraints_p;
+
+ BlockSparsityPattern sparsity_pattern;
+ BlockSparseMatrix<double> system_matrix;
+
+ BlockVector<double> solution;
+ BlockVector<double> system_rhs;
+ BlockVector<double> mf_solution;
+
+ dof_handler.distribute_dofs (fe);
+ dof_handler_u.distribute_dofs (fe_u);
+ dof_handler_p.distribute_dofs (fe_p);
+ std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
+ stokes_sub_blocks[dim] = 1;
+ DoFRenumbering::component_wise (dof_handler, stokes_sub_blocks);
+
+ std::set<unsigned char> no_normal_flux_boundaries;
+ no_normal_flux_boundaries.insert (0);
+ no_normal_flux_boundaries.insert (1);
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+ VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
+ no_normal_flux_boundaries,
+ constraints, mapping);
+ constraints.close ();
+ DoFTools::make_hanging_node_constraints (dof_handler_u,
+ constraints_u);
+ VectorTools::compute_no_normal_flux_constraints (dof_handler_u, 0,
+ no_normal_flux_boundaries,
+ constraints_u, mapping);
+ constraints_u.close ();
+ DoFTools::make_hanging_node_constraints (dof_handler_p,
+ constraints_p);
+ constraints_p.close ();
+
+ std::vector<types::global_dof_index> dofs_per_block (2);
+ DoFTools::count_dofs_per_block (dof_handler, dofs_per_block,
+ stokes_sub_blocks);
+
+ //std::cout << "Number of active cells: "
+ // << triangulation.n_active_cells()
+ // << std::endl
+ // << "Number of degrees of freedom: "
+ // << dof_handler.n_dofs()
+ // << " (" << n_u << '+' << n_p << ')'
+ // << std::endl;
+
+ {
+ BlockCompressedSimpleSparsityPattern csp (2,2);
+
+ for (unsigned int d=0; d<2; ++d)
+ for (unsigned int e=0; e<2; ++e)
+ csp.block(d,e).reinit (dofs_per_block[d], dofs_per_block[e]);
+
+ csp.collect_sizes();
+
+ DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+ sparsity_pattern.copy_from (csp);
+ }
+
+ system_matrix.reinit (sparsity_pattern);
+
+ // this is from step-22
+ {
+ QGauss<dim> quadrature_formula(fe_degree+2);
+
+ FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+ update_values |
+ update_JxW_values |
+ update_gradients);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ const FEValuesExtractors::Vector velocities (0);
+ const FEValuesExtractors::Scalar pressure (dim);
+
+ std::vector<SymmetricTensor<2,dim> > phi_grads_u (dofs_per_cell);
+ std::vector<double> div_phi_u (dofs_per_cell);
+ std::vector<double> phi_p (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit (cell);
+ local_matrix = 0;
+
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+ div_phi_u[k] = fe_values[velocities].divergence (k, q);
+ phi_p[k] = fe_values[pressure].value (k, q);
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<=i; ++j)
+ {
+ local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
+ - div_phi_u[i] * phi_p[j]
+ - phi_p[i] * div_phi_u[j])
+ * fe_values.JxW(q);
+ }
+ }
+ }
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+ local_matrix(i,j) = local_matrix(j,i);
+
+ cell->get_dof_indices (local_dof_indices);
+ constraints.distribute_local_to_global (local_matrix,
+ local_dof_indices,
+ system_matrix);
+ }
+ }
+
+
+ solution.reinit (2);
+ for (unsigned int d=0; d<2; ++d)
+ solution.block(d).reinit (dofs_per_block[d]);
+ solution.collect_sizes ();
+
+ system_rhs.reinit (solution);
+ mf_solution.reinit (solution);
+
+ // fill system_rhs with random numbers
+ for (unsigned int j=0; j<system_rhs.block(0).size(); ++j)
+ if (constraints_u.is_constrained(j) == false)
+ {
+ const double val = -1 + 2.*(double)Testing::rand()/double(RAND_MAX);
+ system_rhs.block(0)(j) = val;
+ }
+ for (unsigned int j=0; j<system_rhs.block(1).size(); ++j)
+ if (constraints_p.is_constrained(j) == false)
+ {
+ const double val = -1 + 2.*(double)Testing::rand()/double(RAND_MAX);
+ system_rhs.block(1)(j) = val;
+ }
+
+ // setup matrix-free structure
+ {
+ std::vector<const DoFHandler<dim>*> dofs;
+ dofs.push_back(&dof_handler_u);
+ dofs.push_back(&dof_handler_p);
+ std::vector<const ConstraintMatrix *> constraints;
+ constraints.push_back (&constraints_u);
+ constraints.push_back (&constraints_p);
+ QGauss<1> quad(fe_degree+2);
+ // no parallelism
+ mf_data.reinit (mapping, dofs, constraints, quad,
+ typename MatrixFree<dim>::AdditionalData
+ (MPI_COMM_WORLD,
+ MatrixFree<dim>::AdditionalData::none));
+ }
+
+ system_matrix.vmult (solution, system_rhs);
+
+ MatrixFreeTest<dim,fe_degree,BlockVector<double> > mf (mf_data);
+ mf.vmult (mf_solution, system_rhs);
+
+ // Verification
+ mf_solution -= solution;
+ const double error = mf_solution.linfty_norm();
+ const double relative = solution.linfty_norm();
+ deallog << "Verification fe degree " << fe_degree << ": "
+ << error/relative << std::endl << std::endl;
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+
+ {
+ deallog << std::endl << "Test with doubles" << std::endl << std::endl;
+ deallog.threshold_double(1.e-12);
+ deallog.push("2d");
+ test<2,1>();
+ test<2,2>();
+ test<2,3>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3,1>();
+ deallog.pop();
+ }
+}
+
--- /dev/null
+
+DEAL::
+DEAL::Test with doubles
+DEAL::
+DEAL:2d::Verification fe degree 1: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 2: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 3: 0
+DEAL:2d::
+DEAL:3d::Verification fe degree 1: 0
+DEAL:3d::
projects / dealii-svn.git / commitdiff