<h3>Specific improvements</h3>
<ol>
+ <li> New: There is now a class FEEvaluationQ_DG0 that does
+ optimized matrix-free evaluation for FE_Q_DG0 elements.
+ <br>
+ (Martin Kronbichler, 2014/06/13)
+ </li>
+
<li> Bugfix: Filter libclang_rt* from the PETSc link line.
<br>
(Matthias Maier, 2014/06/04)
}
+
template < class T >
inline
void
#include <deal.II/base/qprojector.h>
#include <deal.II/base/polynomial_space.h>
#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/tensor_product_polynomials_const.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_poly.h>
{
return poly.get_numbering_inverse();
}
+
+ template <int dim>
+ inline
+ std::vector<unsigned int>
+ get_poly_space_numbering (const TensorProductPolynomialsConst<dim> &poly)
+ {
+ return poly.get_numbering();
+ }
+
+ template <int dim>
+ inline
+ std::vector<unsigned int>
+ get_poly_space_numbering_inverse (const TensorProductPolynomialsConst<dim> &poly)
+ {
+ return poly.get_numbering_inverse();
+ }
}
/**
* 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)];
+ VectorizedArray<Number> my_data_array[n_components*(dofs_per_cell+1+(dim*dim+2*dim+1)*n_q_points)];
+
+ /**
+ * Checks the number of cell dofs with the value that we expect
+ */
+ void check_dofs_per_cell (const unsigned int dofs_per_cell) const;
private:
/**
* Sets the pointers from the data array to values_dof, etc.
*/
void set_data_pointers();
+
+ /**
+ * Checks that the element number of degrees of freedom given by the
+ * template arguments (via fe_degree) coincides with the number of degrees
+ * of freedom in the stored shape values.
+ */
+ void check_template_arguments(const unsigned int fe_no);
};
private:
/**
- * Fills the fields shapve_???_evenodd, called in the constructor.
+ * Fills the fields shape_???_evenodd, called in the constructor.
*/
void compute_even_odd_factors();
};
};
+
+/**
+ * The class that provides all functions necessary to evaluate functions at
+ * quadrature points and cell integrations. In functionality, this class is
+ * similar to FEValues<dim>, however, it includes a lot of specialized
+ * functions that make it much faster (between 5 and 500 times as fast,
+ * depending on the polynomial order). Access to the data fields is provided
+ * through functionality in the class FEEvaluationAccess.
+ *
+ * This class is an extension of FEEvaluation to work with continuous elements
+ * supplemented with a single discontinuous degree of freedom, i.e., FE_Q_DG0.
+ *
+ * @author Martin Kronbichler, 2014
+ */
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+ int n_components_ = 1, typename Number = double >
+class FEEvaluationQ_DG0 :
+ public FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+{
+public:
+ typedef FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number> BaseClass;
+ typedef Number number_type;
+ typedef typename BaseClass::value_type value_type;
+ typedef typename BaseClass::gradient_type gradient_type;
+ static const unsigned int dimension = dim;
+ static const unsigned int n_components = n_components_;
+ static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell + 1;
+ static const unsigned int n_q_points = BaseClass::n_q_points;
+
+ /**
+ * Constructor. Takes all data stored in MatrixFree. If applied to problems
+ * with more than one finite element or more than one quadrature formula
+ * selected during construction of @p matrix_free, @p fe_no and @p quad_no
+ * allow to select the appropriate components.
+ */
+ FEEvaluationQ_DG0 (const MatrixFree<dim,Number> &matrix_free,
+ const unsigned int fe_no = 0,
+ const unsigned int quad_no = 0);
+
+ /**
+ * Constructor that comes with reduced functionality and works similar as
+ * FEValues. The user has to provide a structure of type MappingFEEvaluation
+ * and a DoFHandler in order to allow for reading out the finite element
+ * data. It uses the data provided by dof_handler.get_fe(). If the element
+ * is vector-valued, the optional argument allows to specify the index of
+ * the base element (as long as the element is primitive, non-primitive are
+ * not supported currently).
+ *
+ * With initialization from a FEValues object, no call to a reinit method of
+ * this class is necessary. Instead, it is enough if the geometry is
+ * initialized to a given cell iterator. It can also read from or write to
+ * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
+ * types (which is less efficient with MPI since index translation has to be
+ * done), but of course only for one cell at a time. Hence, a kernel using
+ * this method does not vectorize over several elements (which is most
+ * efficient for vector operations), but only possibly within the element if
+ * the evaluate/integrate routines are combined (e.g. for computing cell
+ * matrices).
+ */
+ FEEvaluationQ_DG0 (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component = 0);
+
+ /**
+ * Copy constructor
+ */
+ FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &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
+ * points of the unit cell. The function arguments specify which parts shall
+ * actually be computed. Needs to be called before the functions @p
+ * get_value(), @p get_gradient() or @p get_laplacian give useful
+ * information (unless these values have been set manually).
+ */
+ void evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl = false);
+
+ /**
+ * This function takes the values and/or gradients that are stored on
+ * quadrature points, tests them by all the basis functions/gradients on the
+ * cell and performs the cell integration. The two function arguments @p
+ * integrate_val and @p integrate_grad are used to enable/disable some of
+ * values or gradients.
+ */
+ void integrate (const bool integrate_val,
+ const bool integrate_grad);
+};
+
+
/*----------------------- Inline functions ----------------------------------*/
#ifndef DOXYGEN
const unsigned int quad_no)
:
BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points)
+{
+ check_template_arguments(fe_no);
+ 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 MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component)
+ :
+ BaseClass (geometry, dof_handler, first_selected_component)
+{
+ check_template_arguments(numbers::invalid_unsigned_int);
+ 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>
+::check_template_arguments(const unsigned int fe_no)
+{
#ifdef DEBUG
// print error message when the dimensions do not match. Propose a possible
// fix
- if ((dofs_per_cell != this->data->dofs_per_cell &&
- internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value !=
- this->data->dofs_per_cell)
+ if (fe_degree != this->data->fe_degree
||
n_q_points != this->data->n_q_points)
{
message += "," + Utilities::int_to_string(n_components);
message += ",Number>(data, ";
message += Utilities::int_to_string(fe_no) + ", ";
- message += Utilities::int_to_string(quad_no) + ")\n";
+ message += Utilities::int_to_string(this->quad_no) + ")\n";
// check whether some other vector component has the correct number of
// points
unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
proposed_quad_comp = numbers::invalid_unsigned_int;
- if (dofs_per_cell == this->data->dofs_per_cell)
- proposed_dof_comp = fe_no;
- else
- for (unsigned int no=0; no<this->matrix_info->n_components(); ++no)
- if (this->matrix_info->get_dof_info(no).dofs_per_cell[this->active_fe_index]
- == dofs_per_cell)
- {
- proposed_dof_comp = no;
- break;
- }
- if (n_q_points ==
- this->mapping_info->mapping_data_gen[quad_no].n_q_points[this->active_quad_index])
- proposed_quad_comp = quad_no;
- else
- for (unsigned int no=0; no<this->mapping_info->mapping_data_gen.size(); ++no)
- if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index]
- == n_q_points)
- {
- proposed_quad_comp = no;
- break;
- }
+ if (fe_no != numbers::invalid_unsigned_int)
+ {
+ if (fe_degree == this->data->fe_degree)
+ proposed_dof_comp = fe_no;
+ else
+ for (unsigned int no=0; no<this->matrix_info->n_components(); ++no)
+ if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree
+ == fe_degree)
+ {
+ proposed_dof_comp = no;
+ break;
+ }
+ if (n_q_points ==
+ this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index])
+ proposed_quad_comp = this->quad_no;
+ else
+ for (unsigned int no=0; no<this->mapping_info->mapping_data_gen.size(); ++no)
+ if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index]
+ == n_q_points)
+ {
+ proposed_quad_comp = no;
+ break;
+ }
+ }
if (proposed_dof_comp != numbers::invalid_unsigned_int &&
proposed_quad_comp != numbers::invalid_unsigned_int)
{
correct_pos = " ^ ";
else
correct_pos = " ";
- if (proposed_quad_comp != quad_no)
+ if (proposed_quad_comp != this->quad_no)
correct_pos += " ^\n";
else
correct_pos += " \n";
}
// ok, did not find the numbers specified by the template arguments in
// the given list. Suggest correct template arguments
- const unsigned int proposed_fe_degree = static_cast<unsigned int>(std::pow(1.001*this->data->dofs_per_cell,1./dim))-1;
const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(std::pow(1.001*this->data->n_q_points,1./dim));
message += "Wrong template arguments:\n";
message += " Did you mean FEEvaluation<dim,";
- message += Utilities::int_to_string(proposed_fe_degree) + ",";
+ message += Utilities::int_to_string(this->data->fe_degree) + ",";
message += Utilities::int_to_string(proposed_n_q_points_1d);
message += "," + Utilities::int_to_string(n_components);
message += ",Number>(data, ";
message += Utilities::int_to_string(fe_no) + ", ";
- message += Utilities::int_to_string(quad_no) + ")?\n";
+ message += Utilities::int_to_string(this->quad_no) + ")?\n";
std::string correct_pos;
- if (proposed_fe_degree != fe_degree)
+ if (this->data->fe_degree != fe_degree)
correct_pos = " ^";
else
correct_pos = " ";
correct_pos += " \n";
message += " " + correct_pos;
- Assert (dofs_per_cell == this->data->dofs_per_cell &&
+ Assert (fe_degree == this->data->fe_degree &&
n_q_points == this->data->n_q_points,
ExcMessage(message));
}
- AssertDimension (n_q_points,
- this->mapping_info->mapping_data_gen[this->quad_no].
- n_q_points[this->active_quad_index]);
- AssertDimension (this->data->dofs_per_cell * this->n_fe_components,
- this->dof_info->dofs_per_cell[this->active_fe_index]);
+ if (fe_no != numbers::invalid_unsigned_int)
+ {
+ AssertDimension (n_q_points,
+ this->mapping_info->mapping_data_gen[this->quad_no].
+ n_q_points[this->active_quad_index]);
+ AssertDimension (this->data->dofs_per_cell * this->n_fe_components,
+ this->dof_info->dofs_per_cell[this->active_fe_index]);
+ }
#endif
}
-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 MappingFEEvaluation<dim,Number> &geometry,
- const DoFHandler<dim> &dof_handler,
- const unsigned int first_selected_component)
- :
- BaseClass (geometry, dof_handler, first_selected_component)
-{
- set_data_pointers();
-
- Assert ((dofs_per_cell == this->data->dofs_per_cell ||
- internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value ==
- this->data->dofs_per_cell
- )
- &&
- n_q_points == this->data->n_q_points,
- ExcMessage("Underlying element and template arguments do not match"));
-}
-
-
-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
FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
::set_data_pointers()
{
+ AssertIndexRange(this->data->dofs_per_cell, dofs_per_cell+2);
+ const unsigned int desired_dofs_per_cell = this->data->dofs_per_cell;
+
// 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];
+ this->values_dofs[c] = &my_data_array[c*desired_dofs_per_cell];
+ this->values_quad[c] = &my_data_array[n_components*desired_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)
+ this->gradients_quad[c][d] = &my_data_array[n_components*(desired_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)
+ this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+
+ desired_dofs_per_cell)
+
(c*(dim*dim+dim)+d)*n_q_points];
}
}
+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>
+::check_dofs_per_cell (const unsigned int given_dofs_per_cell) const
+{
+#ifdef DEBUG
+ if (given_dofs_per_cell != this->data->dofs_per_cell)
+ {
+ std::ostringstream str;
+ str << "Dofs per cell in FEEvaluation* class does not match the element ("
+ << given_dofs_per_cell << " != " << this->data->dofs_per_cell << "). "
+ << "Suggestion: ";
+ if (this->data->dofs_per_cell == dofs_per_cell)
+ str << "FEEvaluation/FEEvaluationGeneral";
+ else if (this->data->dofs_per_cell ==
+ internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value)
+ str << "FEEvaluationDGP";
+ else if (this->data->dofs_per_cell == dofs_per_cell+1)
+ str << "FEEvaluationQ_DG0";
+ else
+ str << "No matching suggestion found";
+
+ Assert(given_dofs_per_cell == this->data->dofs_per_cell,
+ ExcMessage(str.str().c_str()));
+ }
+#endif
+}
+
+
namespace internal
{
out += nn;
break;
case 1:
+ ++in;
+ ++out;
+ // faces 2 and 3 in 3D use local coordinate system zx, which
+ // is the other way around compared to the tensor
+ // product. Need to take that into account.
+ if (dim == 3)
+ {
+ if (dof_to_quad)
+ out += fe_degree;
+ else
+ in += fe_degree;
+ }
+ break;
case 2:
++in;
++out;
Assert (false, ExcNotImplemented());
}
}
- if (face_direction == 1)
+ if (face_direction == 1 && dim == 3)
{
in += mm*(mm-1);
out += nn*(nn-1);
+ // adjust for local coordinate system zx
+ if (dof_to_quad)
+ out -= (fe_degree+1)*(fe_degree+1)-1;
+ else
+ in -= (fe_degree+1)*(fe_degree+1)-1;
}
}
}
const bool evaluate_grad,
const bool evaluate_lapl)
{
+ this->check_dofs_per_cell(dofs_per_cell);
Assert (this->dof_values_initialized == true,
internal::ExcAccessToUninitializedField());
internal::do_evaluate (*this, this->values_dofs, this->values_quad,
::integrate (const bool integrate_val,
const bool integrate_grad)
{
+ this->check_dofs_per_cell(dofs_per_cell);
if (integrate_val == true)
Assert (this->values_quad_submitted == true,
internal::ExcAccessToUninitializedField());
const bool evaluate_grad,
const bool evaluate_lapl)
{
+ this->check_dofs_per_cell(dofs_per_cell);
Assert (this->dof_values_initialized == true,
internal::ExcAccessToUninitializedField());
internal::do_evaluate (*this, this->values_dofs, this->values_quad,
FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
::integrate (bool integrate_val,bool integrate_grad)
{
+ this->check_dofs_per_cell(dofs_per_cell);
if (integrate_val == true)
Assert (this->values_quad_submitted == true,
internal::ExcAccessToUninitializedField());
const bool evaluate_grad,
const bool evaluate_lapl)
{
+ this->check_dofs_per_cell(dofs_per_cell);
Assert (this->cell != numbers::invalid_unsigned_int,
ExcNotInitialized());
Assert (this->dof_values_initialized == true,
FEEvaluationGL<dim,fe_degree,n_components_,Number>
::integrate (const bool integrate_val, const bool integrate_grad)
{
+ this->check_dofs_per_cell(dofs_per_cell);
Assert (this->cell != numbers::invalid_unsigned_int,
ExcNotInitialized());
if (integrate_val == true)
:
BaseClass (data_in, fe_no, quad_no)
{
- // reset values_dofs pointer as it has wider gaps in the base class
- for (unsigned int c=0; c<n_components_; ++c)
- this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
+ AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+ n_components * dofs_per_cell);
}
:
BaseClass (geometry, dof_handler, first_selected_component)
{
- for (unsigned int c=0; c<n_components_; ++c)
- this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
+ AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+ n_components * dofs_per_cell);
}
:
BaseClass (other)
{
- for (unsigned int c=0; c<n_components_; ++c)
- this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
}
const bool evaluate_grad,
const bool evaluate_lapl)
{
+ this->check_dofs_per_cell(dofs_per_cell);
Assert (this->dof_values_initialized == true,
internal::ExcAccessToUninitializedField());
FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
::integrate (bool integrate_val,bool integrate_grad)
{
+ this->check_dofs_per_cell(dofs_per_cell);
if (integrate_val == true)
Assert (this->values_quad_submitted == true,
internal::ExcAccessToUninitializedField());
+/*------------------------- FEEvaluationQ_DG0 -------------------------------*/
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationQ_DG0 (const MatrixFree<dim,Number> &data_in,
+ const unsigned int fe_no,
+ const unsigned int quad_no)
+ :
+ BaseClass (data_in, fe_no, quad_no)
+{
+ AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+ n_components * dofs_per_cell);
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationQ_DG0 (const MappingFEEvaluation<dim,Number> &geometry,
+ const DoFHandler<dim> &dof_handler,
+ const unsigned int first_selected_component)
+ :
+ BaseClass (geometry, dof_handler, first_selected_component)
+{
+ AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+ n_components * dofs_per_cell);
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other)
+ :
+ BaseClass (other)
+{
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+void
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::evaluate (const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl)
+{
+ this->check_dofs_per_cell(dofs_per_cell);
+ Assert (this->dof_values_initialized == true,
+ internal::ExcAccessToUninitializedField());
+ internal::do_evaluate (*this, this->values_dofs, this->values_quad,
+ this->gradients_quad, this->hessians_quad,
+ evaluate_val, evaluate_grad, evaluate_lapl,
+ internal::int2type<dim>());
+ if (evaluate_val)
+ for (unsigned int c=0; c<n_components; ++c)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ this->values_quad[c][q] += this->values_dofs[c][dofs_per_cell-1];
+
+#ifdef DEBUG
+ if (evaluate_val == true)
+ this->values_quad_initialized = true;
+ if (evaluate_grad == true)
+ this->gradients_quad_initialized = true;
+ if (evaluate_lapl == true)
+ this->hessians_quad_initialized = true;
+#endif
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+ typename Number>
+inline
+void
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::integrate (bool integrate_val,bool integrate_grad)
+{
+ this->check_dofs_per_cell(dofs_per_cell);
+ if (integrate_val == true)
+ Assert (this->values_quad_submitted == true,
+ internal::ExcAccessToUninitializedField());
+ if (integrate_grad == true)
+ Assert (this->gradients_quad_submitted == true,
+ internal::ExcAccessToUninitializedField());
+
+ if (integrate_val)
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ this->values_dofs[c][dofs_per_cell-1] = this->values_quad[c][0];
+ for (unsigned int q=1; q<n_q_points; ++q)
+ this->values_dofs[c][dofs_per_cell-1] += this->values_quad[c][q];
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ this->values_dofs[c][dofs_per_cell-1] = VectorizedArray<Number>();
+
+ internal::do_integrate (*this, this->values_dofs, this->values_quad,
+ this->gradients_quad, integrate_val, integrate_grad,
+ internal::int2type<dim>());
+
+#ifdef DEBUG
+ this->dof_values_initialized = true;
+#endif
+}
+
+
+
#endif // ifndef DOXYGEN
const unsigned int base_element = 0);
/**
- * Returns the memory consumption of this
- * class in bytes.
+ * Returns the memory consumption of this class in bytes.
*/
std::size_t memory_consumption () const;
/**
- * Stores the shape values of the 1D finite
- * element evaluated on all 1D quadrature
- * points in vectorized format, i.e., as an
- * array of
- * VectorizedArray<dim>::n_array_elements
- * equal elements. The length of this array is
- * <tt>n_dofs_1d * n_q_points_1d</tt> and
- * quadrature points are the index running
- * fastest.
+ * Stores the shape values of the 1D finite element evaluated on all 1D
+ * quadrature points in vectorized format, i.e., as an array of
+ * VectorizedArray<dim>::n_array_elements equal elements. The length of
+ * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
+ * points are the index running fastest.
*/
AlignedVector<VectorizedArray<Number> > shape_values;
/**
- * Stores the shape gradients of the 1D finite
- * element evaluated on all 1D quadrature
- * points in vectorized format, i.e., as an
- * array of
- * VectorizedArray<dim>::n_array_elements
- * equal elements. The length of this array is
- * <tt>n_dofs_1d * n_q_points_1d</tt> and
- * quadrature points are the index running
- * fastest.
+ * Stores the shape gradients of the 1D finite element evaluated on all
+ * 1D quadrature points in vectorized format, i.e., as an array of
+ * VectorizedArray<dim>::n_array_elements equal elements. The length of
+ * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
+ * points are the index running fastest.
*/
AlignedVector<VectorizedArray<Number> > shape_gradients;
/**
- * Stores the shape Hessians of the 1D finite
- * element evaluated on all 1D quadrature
- * points in vectorized format, i.e., as an
- * array of
- * VectorizedArray<dim>::n_array_elements
- * equal elements. The length of this array is
- * <tt>n_dofs_1d * n_q_points_1d</tt> and
- * quadrature points are the index running
- * fastest.
+ * Stores the shape Hessians of the 1D finite element evaluated on all
+ * 1D quadrature points in vectorized format, i.e., as an array of
+ * VectorizedArray<dim>::n_array_elements equal elements. The length of
+ * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
+ * points are the index running fastest.
*/
AlignedVector<VectorizedArray<Number> > shape_hessians;
/**
- * Stores the indices from cell DoFs to face
- * DoFs. The rows go through the
- * <tt>2*dim</tt> faces, and the columns the
- * DoFs on the faces.
+ * Stores the indices from cell DoFs to face DoFs. The rows go through
+ * the <tt>2*dim</tt> faces, and the columns the DoFs on the faces.
*/
dealii::Table<2,unsigned int> face_indices;
std::vector<Number> subface_value[2];
/**
- * Non-vectorized version of shape
- * values. Needed when evaluating face info.
+ * Non-vectorized version of shape values. Needed when evaluating face
+ * info.
*/
std::vector<Number> shape_values_number;
/**
- * Non-vectorized version of shape
- * gradients. Needed when evaluating face
- * info.
+ * Non-vectorized version of shape gradients. Needed when evaluating
+ * face info.
*/
std::vector<Number> shape_gradient_number;
std::vector<unsigned int> lexicographic_numbering;
/**
- * Stores the number of quadrature points in
- * @p dim dimensions for a cell.
+ * Stores the degree of the element.
+ */
+ unsigned int fe_degree;
+
+ /**
+ * Stores the number of quadrature points in @p dim dimensions for a
+ * cell.
*/
unsigned int n_q_points;
/**
- * Stores the number of DoFs per cell in @p
- * dim dimensions.
+ * Stores the number of DoFs per cell in @p dim dimensions.
*/
unsigned int dofs_per_cell;
/**
- * Stores the number of quadrature points per
- * face in @p dim dimensions.
+ * Stores the number of quadrature points per face in @p dim dimensions.
*/
unsigned int n_q_points_face;
/**
- * Stores the number of DoFs per face in @p
- * dim dimensions.
+ * Stores the number of DoFs per face in @p dim dimensions.
*/
unsigned int dofs_per_face;
};
const FiniteElement<dim> &fe_in,
const unsigned int base_element_number)
:
+ fe_degree (0),
n_q_points (0),
dofs_per_cell (0)
{
#include <deal.II/base/polynomials_piecewise.h>
#include <deal.II/fe/fe_poly.h>
#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/matrix_free/shape_info.h>
Assert (fe->n_components() == 1,
ExcMessage("FEEvaluation only works for scalar finite elements."));
+ fe_degree = fe->degree;
- const unsigned int n_dofs_1d = fe->degree+1,
+ const unsigned int n_dofs_1d = fe_degree+1,
n_q_points_1d = quad.size();
// renumber (this is necessary for FE_Q, for example, since there the
const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim>*>(fe);
+ const FE_Q_DG0<dim> *fe_q_dg0 = dynamic_cast<const FE_Q_DG0<dim>*>(fe);
+
if (fe_poly != 0)
scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse();
else if (fe_poly_piece != 0)
else if (fe_dgp != 0)
{
scalar_lexicographic.resize(fe_dgp->dofs_per_cell);
- for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
- scalar_lexicographic[i] = i;
+ for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
+ scalar_lexicographic[i] = i;
}
+ else if (fe_q_dg0 != 0)
+ scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse();
else
Assert(false, ExcNotImplemented());
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 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 test is similar to matrix_vector_06, but implements the operations for
+// FE_QDG0 instead of FE_Q (where there is an additional discontinuous degree
+// of freedom and different routines need to be used). The data is still not
+// very useful because the matrix does not include face terms actually present
+// in an approximation of the Laplacian. It only contains cell terms.
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.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/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/lac/vector.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename VECTOR>
+void
+helmholtz_operator_qdg0 (const MatrixFree<dim,typename VECTOR::value_type> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range)
+{
+ typedef typename VECTOR::value_type Number;
+ FEEvaluationQ_DG0<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ fe_eval.read_dof_values (src);
+ fe_eval.evaluate (true, true, false);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
+ fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
+ }
+ fe_eval.integrate (true,true);
+ fe_eval.distribute_local_to_global (dst);
+ }
+}
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+ typedef VectorizedArray<Number> vector_t;
+
+ MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+ data (data_in)
+ {};
+
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ const std_cxx1x::function<void(const MatrixFree<dim,typename VECTOR::value_type> &,
+ VECTOR &,
+ const VECTOR &,
+ const std::pair<unsigned int,unsigned int> &)>
+ wrap = helmholtz_operator_qdg0<dim,fe_degree,VECTOR>;
+ data.cell_loop (wrap, dst, src);
+ };
+
+private:
+ const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+ const ConstraintMatrix &constraints)
+{
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ MatrixFree<dim,number> mf_data;
+ {
+ const QGauss<1> quad (fe_degree+1);
+ typename MatrixFree<dim,number>::AdditionalData data;
+ data.tasks_parallel_scheme =
+ MatrixFree<dim,number>::AdditionalData::partition_color;
+ data.tasks_block_size = 3;
+
+ mf_data.reinit (dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+ Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
+ Vector<number> in_dist (dof.n_dofs());
+ Vector<number> out_dist (in_dist);
+
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ if (constraints.is_constrained(i))
+ continue;
+ const double entry = Testing::rand()/(double)RAND_MAX;
+ in(i) = entry;
+ in_dist(i) = entry;
+ }
+
+ mf.vmult (out_dist, in_dist);
+
+
+ // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
+ SparsityPattern sparsity;
+ {
+ CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+ DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
+ sparsity.copy_from(csp);
+ }
+ SparseMatrix<double> sparse_matrix (sparsity);
+ {
+ QGauss<dim> quadrature_formula(fe_degree+1);
+
+ FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
+ update_values | update_gradients |
+ update_JxW_values);
+
+ const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof.begin_active(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix = 0;
+ fe_values.reinit (cell);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point)
+ +
+ 10. *
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point)) *
+ fe_values.JxW(q_point));
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+ constraints.distribute_local_to_global (cell_matrix,
+ local_dof_indices,
+ sparse_matrix);
+ }
+ }
+
+ sparse_matrix.vmult (out, in);
+ out -= out_dist;
+ const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
+
+ deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball (tria);
+ static const HyperBallBoundary<dim> boundary;
+ tria.set_boundary (0, boundary);
+ if (dim < 3 || fe_degree < 2)
+ tria.refine_global(1);
+ tria.begin(tria.n_levels()-1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end();
+ for (; cell!=endc; ++cell)
+ if (cell->center().norm()<1e-8)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+ FE_Q_DG0<dim> fe (fe_degree);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+ ConstraintMatrix constraints;
+
+ do_test<dim, fe_degree, double> (dof, constraints);
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+
+ {
+ deallog.threshold_double(5.e-11);
+ deallog.push("2d");
+ test<2,1>();
+ test<2,2>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3,1>();
+ test<3,2>();
+ deallog.pop();
+ }
+}
+
+
--- /dev/null
+
+DEAL:2d::Testing FE_Q_DG0<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q_DG0<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q_DG0<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q_DG0<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// matrix-vector products by comparing with the result of deal.II sparse
+// matrix. No hanging nodes and no other
+// constraints for a vector-valued problem (stokes equations).
+
+#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/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_q_dg0.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "create_mesh.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;
+ FEEvaluation <dim,degree_p+1,degree_p+2,dim,Number> velocity (data, 0);
+ FEEvaluationQ_DG0<dim,degree_p ,degree_p+2,1, Number> pressure (data, 1);
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ velocity.reinit (cell);
+ velocity.read_dof_values (src, 0);
+ velocity.evaluate (false,true,false);
+ pressure.reinit (cell);
+ pressure.read_dof_values (src, dim);
+ pressure.evaluate (true,false,false);
+
+ for (unsigned int q=0; q<velocity.n_q_points; ++q)
+ {
+ SymmetricTensor<2,dim,vector_t> sym_grad_u =
+ velocity.get_symmetric_gradient (q);
+ vector_t pres = pressure.get_value(q);
+ vector_t div = -velocity.get_divergence(q);
+ pressure.submit_value (div, q);
+
+ // subtract p * I
+ for (unsigned int d=0; d<dim; ++d)
+ sym_grad_u[d][d] -= pres;
+
+ velocity.submit_symmetric_gradient(sym_grad_u, q);
+ }
+
+ velocity.integrate (false,true);
+ velocity.distribute_local_to_global (dst, 0);
+ pressure.integrate (true,false);
+ pressure.distribute_local_to_global (dst, dim);
+ }
+ }
+
+
+ void vmult (VectorType &dst,
+ const VectorType &src) const
+ {
+ AssertDimension (dst.size(), dim+1);
+ for (unsigned int d=0; d<dim+1; ++d)
+ dst[d] = 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;
+ create_mesh (triangulation);
+ if (fe_degree == 1)
+ triangulation.refine_global (4-dim);
+ else
+ triangulation.refine_global (3-dim);
+
+ FE_Q<dim> fe_u (fe_degree+1);
+ FE_Q_DG0<dim> fe_p (fe_degree);
+ FESystem<dim> fe (fe_u, 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;
+
+ BlockSparsityPattern sparsity_pattern;
+ BlockSparseMatrix<double> system_matrix;
+
+ BlockVector<double> solution;
+ BlockVector<double> system_rhs;
+ std::vector<Vector<double> > vec1, vec2;
+
+ dof_handler.distribute_dofs (fe);
+ dof_handler_u.distribute_dofs (fe_u);
+ dof_handler_p.distribute_dofs (fe_p);
+ DoFRenumbering::component_wise (dof_handler);
+
+ constraints.close ();
+
+ std::vector<types::global_dof_index> dofs_per_block (dim+1);
+ DoFTools::count_dofs_per_component (dof_handler, dofs_per_block);
+
+ //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 (dim+1,dim+1);
+
+ for (unsigned int d=0; d<dim+1; ++d)
+ for (unsigned int e=0; e<dim+1; ++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);
+
+ solution.reinit (dim+1);
+ for (unsigned int i=0; i<dim+1; ++i)
+ solution.block(i).reinit (dofs_per_block[i]);
+ solution.collect_sizes ();
+
+ system_rhs.reinit (solution);
+
+ vec1.resize (dim+1);
+ vec2.resize (dim+1);
+ vec1[0].reinit (dofs_per_block[0]);
+ vec2[0].reinit (vec1[0]);
+ for (unsigned int i=1; i<dim; ++i)
+ {
+ vec1[i].reinit (vec1[0]);
+ vec2[i].reinit (vec1[0]);
+ }
+ vec1[dim].reinit (dofs_per_block[dim]);
+ vec2[dim].reinit (vec1[dim]);
+
+ // this is from step-22
+ {
+ QGauss<dim> quadrature_formula(fe_degree+2);
+
+ FEValues<dim> fe_values (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);
+ }
+ }
+
+ // first system_rhs with random numbers
+ for (unsigned int i=0; i<dim+1; ++i)
+ for (unsigned int j=0; j<system_rhs.block(i).size(); ++j)
+ {
+ const double val = -1. + 2.*(double)Testing::rand()/double(RAND_MAX);
+ system_rhs.block(i)(j) = val;
+ vec1[i](j) = val;
+ }
+
+ // setup matrix-free structure
+ {
+ std::vector<const DoFHandler<dim>*> dofs;
+ dofs.push_back(&dof_handler_u);
+ dofs.push_back(&dof_handler_p);
+ ConstraintMatrix dummy_constraints;
+ dummy_constraints.close();
+ std::vector<const ConstraintMatrix *> constraints;
+ constraints.push_back (&dummy_constraints);
+ constraints.push_back (&dummy_constraints);
+ QGauss<1> quad(fe_degree+2);
+ mf_data.reinit (dofs, constraints, quad,
+ typename MatrixFree<dim>::AdditionalData
+ (MPI_COMM_WORLD,
+ MatrixFree<dim>::AdditionalData::none));
+ }
+
+ system_matrix.vmult (solution, system_rhs);
+
+ typedef std::vector<Vector<double> > VectorType;
+ MatrixFreeTest<dim,fe_degree,VectorType> mf (mf_data);
+ mf.vmult (vec2, vec1);
+
+ // Verification
+ double error = 0.;
+ for (unsigned int i=0; i<dim+1; ++i)
+ for (unsigned int j=0; j<system_rhs.block(i).size(); ++j)
+ error += std::fabs (solution.block(i)(j)-vec2[i](j));
+ double relative = solution.block(0).l1_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>();
+ test<3,2>();
+ 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::
+DEAL:3d:: Verification fe degree 2: 0
+DEAL:3d::
projects / dealii.git / commitdiff