*/
std::vector<unsigned int> plain_dof_indices;
+ /**
+ * Stores the dimension of the underlying DoFHandler. Since the indices
+ * are not templated, this is the variable that makes the dimension
+ * accessible in the (rare) cases it is needed inside this class.
+ */
+ unsigned int dimension;
+
/**
* Stores the number of components in the DoFHandler where the indices
* have been read from.
constrained_dofs (dof_info_in.constrained_dofs),
row_starts_plain_indices (dof_info_in.row_starts_plain_indices),
plain_dof_indices (dof_info_in.plain_dof_indices),
+ dimension (dof_info_in.dimension),
n_components (dof_info_in.n_components),
dofs_per_cell (dof_info_in.dofs_per_cell),
dofs_per_face (dof_info_in.dofs_per_face),
ghost_dofs.clear();
dofs_per_cell.clear();
dofs_per_face.clear();
+ dimension = 2;
n_components = 0;
row_starts_plain_indices.clear();
plain_dof_indices.clear();
*/
unsigned int cell_data_number;
- /**
- * If the present cell chunk for vectorization is not completely filled up
- * with data, this field stores how many physical cells are underlying. Is
- * between 1 and VectorizedArray<Number>::n_array_elements-1 (inclusive).
- */
- unsigned int n_irreg_components_filled;
-
- /**
- * Stores whether the present cell chunk used in vectorization is not
- * completely filled up with physical cells. E.g. if vectorization dictates
- * that four cells should be worked with but only three physical cells are
- * left, this flag will be set to true, otherwise to false. Mainly used for
- * internal checking when reading from vectors or writing to vectors.
- */
- bool at_irregular_cell;
-
/**
* Debug information to track whether dof values have been initialized
* before accessed. Used to control exceptions when uninitialized data is
jacobian_grad_upper(0),
cell (numbers::invalid_unsigned_int),
cell_type (internal::MatrixFreeFunctions::undefined),
- cell_data_number (0),
- n_irreg_components_filled (0),
- at_irregular_cell (false)
+ cell_data_number (0)
{
- Assert (matrix_info.indices_initialized() == true,
- ExcNotInitialized());
Assert (matrix_info.mapping_initialized() == true,
ExcNotInitialized());
AssertDimension (matrix_info.get_size_info().vectorization_length,
VectorizedArray<Number>::n_array_elements);
Assert (n_fe_components == 1 ||
+ n_components == 1 ||
n_components == n_fe_components,
ExcMessage ("The underlying FE is vector-valued. In this case, the "
"template argument n_components must be a the same "
}
}
- n_irreg_components_filled =
- std_cxx1x::get<2>(dof_info.row_starts[cell_in]);
- at_irregular_cell = n_irreg_components_filled > 0;
#ifdef DEBUG
dof_values_initialized = false;
values_quad_initialized = false;
// into the local data field or write local data into the vector. Certain
// operations are no-ops for the given use case.
+ Assert (matrix_info.indices_initialized() == true,
+ ExcNotInitialized());
Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
// loop over all local dofs. ind_local holds local number on cell, index
dof_info.end_indicators(cell);
unsigned int ind_local = 0;
+ const unsigned int n_irreg_components_filled =
+ std_cxx1x::get<2>(dof_info.row_starts[cell]);
+ const bool at_irregular_cell = n_irreg_components_filled > 0;
+
// scalar case (or case when all components have the same degrees of freedom
// and sit on a different vector each)
if (n_fe_components == 1)
{
// this is different from the other three operations because we do not use
// constraints here, so this is a separate function.
+ Assert (matrix_info.indices_initialized() == true,
+ ExcNotInitialized());
Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
// points to the global indices stored in index_local_to_global
const unsigned int *dof_indices = dof_info.begin_indices_plain(cell);
+ const unsigned int n_irreg_components_filled =
+ std_cxx1x::get<2>(dof_info.row_starts[cell]);
+ const bool at_irregular_cell = n_irreg_components_filled > 0;
+
// scalar case (or case when all components have the same degrees of freedom
// and sit on a different vector each)
if (n_fe_components == 1)
initialize_indices (constraint, locally_owned_set);
}
+ // initialize bare structures
+ else if (dof_info.size() != dof_handler.size())
+ {
+ initialize_dof_handlers(dof_handler, additional_data.level_mg_handler);
+ std::vector<unsigned int> dummy;
+ size_info.make_layout (cell_level_index.size(),
+ VectorizedArray<Number>::n_array_elements,
+ dummy, dummy);
+ for (unsigned int i=0; i<dof_info.size(); ++i)
+ {
+ dof_info[i].dimension = dim;
+ dof_info[i].n_components = dof_handler[i]->get_fe().element_multiplicity(0);
+ dof_info[i].dofs_per_cell.push_back(dof_handler[i]->get_fe().dofs_per_cell);
+ dof_info[i].row_starts.resize(size_info.n_macro_cells+1);
+ std_cxx1x::get<2>(dof_info[i].row_starts.back()) =
+ cell_level_index.size() % VectorizedArray<Number>::n_array_elements;
+
+ // if indices are not initialized, the cell_level_index might not be
+ // divisible by the vectorization length. But it must be for
+ // mapping_info...
+ while (cell_level_index.size() % VectorizedArray<Number>::n_array_elements
+ != 0)
+ cell_level_index.push_back(cell_level_index.back());
+ }
+ }
+
// Reads out the FE information and stores the shape function values,
// gradients and Hessians for quadrature points.
const unsigned int n_fe = dof_handler.size();
initialize_indices (constraint, locally_owned_set);
}
+ // initialize bare structures
+ else if (dof_info.size() != dof_handler.size())
+ {
+ initialize_dof_handlers(dof_handler, additional_data.level_mg_handler);
+ std::vector<unsigned int> dummy;
+ size_info.make_layout (cell_level_index.size(),
+ VectorizedArray<Number>::n_array_elements,
+ dummy, dummy);
+ for (unsigned int i=0; i<dof_info.size(); ++i)
+ {
+ Assert(dof_handler[i]->get_fe().size() == 1, ExcNotImplemented());
+ dof_info[i].dimension = dim;
+ dof_info[i].n_components = dof_handler[i]->get_fe()[0].element_multiplicity(0);
+ dof_info[i].dofs_per_cell.push_back(dof_handler[i]->get_fe()[0].dofs_per_cell);
+ dof_info[i].row_starts.resize(size_info.n_macro_cells+1);
+ std_cxx1x::get<2>(dof_info[i].row_starts.back()) =
+ cell_level_index.size() % VectorizedArray<Number>::n_array_elements;
+
+ // if indices are not initialized, the cell_level_index might not be
+ // divisible by the vectorization length. But it must be for
+ // mapping_info...
+ while (cell_level_index.size() % VectorizedArray<Number>::n_array_elements
+ != 0)
+ cell_level_index.push_back(cell_level_index.back());
+ }
+ }
+
// Reads out the FE information and stores the shape function values,
// gradients and Hessians for quadrature points.
const unsigned int n_components = dof_handler.size();
// cache number of finite elements and dofs_per_cell
dof_info[no].dofs_per_cell.push_back (fe.dofs_per_cell);
dof_info[no].dofs_per_face.push_back (fe.dofs_per_face);
- dof_info[no].n_components = n_fe_components;
-
+ dof_info[no].dimension = dim;
+ dof_info[no].n_components = n_fe_components;
// get permutation that gives lexicographic renumbering of the cell
// dofs renumber (this is necessary for FE_Q, for example, since
virtual ~MatrixFreeTest ()
{}
-
+
// make function virtual to allow derived
// classes to define a different function
virtual void
deallog.pop();
}
}
-
--- /dev/null
+//------------------ no_index_initialize.cc ------------------------
+// $Id$
+// Version: $Name$
+//
+//------------------ no_index_initialize.cc ------------------------
+
+
+// check that FEEvaluation can be evaluated without indices initialized (and
+// throws an exception when trying to read/write from/to vectors)
+
+#include "../tests.h"
+
+
+std::ofstream logfile("no_index_initialize/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/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/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+
+// forward declare this function
+template <int dim, int fe_degree>
+void test ();
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1, typename Number=double>
+class MatrixFreeTest
+{
+ public:
+ bool read_vector;
+
+ MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+ read_vector(false),
+ data (data_in)
+ {};
+
+ void
+ operator () (const MatrixFree<dim,Number> &data,
+ Vector<Number> &,
+ const Vector<Number> &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval (data);
+
+ for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
+ {
+ fe_eval.reinit (cell);
+
+ // do not read from vector as that is disallowed
+ if (read_vector)
+ fe_eval.read_dof_values(src);
+ else
+ for (unsigned int e=0; e<fe_eval.dofs_per_cell; ++e)
+ fe_eval.submit_dof_value(make_vectorized_array<Number>(1.), e);
+ fe_eval.evaluate (true,true,true);
+
+ // values should evaluate to one, derivatives to zero
+ for (unsigned int q=0; q<fe_eval.n_q_points; ++q)
+ for (unsigned int d=0; d<VectorizedArray<Number>::n_array_elements; ++d)
+ {
+ Assert(fe_eval.get_value(q)[d] == 1., ExcInternalError());
+ for (unsigned int e=0; e<dim; ++e)
+ Assert(fe_eval.get_gradient(q)[e][d] == 0., ExcInternalError());
+ Assert(fe_eval.get_laplacian(q)[d] == 0., ExcInternalError());
+ }
+ }
+ }
+
+
+
+ void test_functions (const Vector<Number> &src) const
+ {
+ Vector<Number> dst_dummy;
+ data.cell_loop (&MatrixFreeTest::operator(), this, dst_dummy, src);
+ deallog << "OK" << std::endl;
+ };
+
+protected:
+ const MatrixFree<dim,Number> &data;
+};
+
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+ const ConstraintMatrix&constraints)
+{
+ // use this for info on problem
+ //std::cout << "Number of cells: " << dof.get_tria().n_active_cells()
+ // << std::endl;
+ //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+ //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+ Vector<number> solution (dof.n_dofs());
+
+ // create vector with random entries
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ if(constraints.is_constrained(i))
+ continue;
+ const double entry = rand()/(double)RAND_MAX;
+ solution(i) = entry;
+ }
+
+ constraints.distribute(solution);
+ 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::none;
+ data.mapping_update_flags = update_gradients | update_second_derivatives;
+ data.initialize_indices = false;
+ mf_data.reinit (dof, constraints, quad, data);
+ }
+
+ deallog << "Testing " << dof.get_fe().get_name() << " without read" << std::endl;
+ MatrixFreeTest<dim,fe_degree,fe_degree+1,number> mf (mf_data);
+ mf.test_functions(solution);
+
+ deallog << "Testing " << dof.get_fe().get_name() << " with read" << std::endl;
+ try
+ {
+ mf.read_vector = true;
+ mf.test_functions(solution);
+ }
+ catch (ExceptionBase &e)
+ {
+ deallog << e.get_exc_name() << std::endl;
+ }
+}
+
+
+int main ()
+{
+ deal_II_exceptions::disable_abort_on_exception();
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog << std::setprecision (3);
+ test<2,1>();
+}
+
+
+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);
+ // refine first and last cell
+ tria.begin(tria.n_levels()-1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ tria.refine_global (4-dim);
+
+ FE_Q<dim> fe (fe_degree);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ constraints.close();
+
+ do_test <dim, fe_degree, double> (dof, constraints);
+}
--- /dev/null
+
+DEAL::Testing FE_Q<2>(1) without read
+DEAL::OK
+DEAL::Testing FE_Q<2>(1) with read
+DEAL::ExcNotInitialized()
--- /dev/null
+//------------------ update_mapping_only.cc ------------------------
+// $Id$
+// Version: $Name$
+//
+//------------------ update_mapping_only.cc ------------------------
+
+
+// correctness of matrix-free initialization when the mapping changes and the
+// indices are not updated, only the mapping info
+
+#include "../tests.h"
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+
+
+std::ofstream logfile("update_mapping_only/output");
+
+#include "get_functions_common.h"
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube (tria);
+ tria.refine_global(1);
+
+ FE_Q<dim> fe (fe_degree);
+ DoFHandler<dim> dof (tria);
+ dof.distribute_dofs(fe);
+
+ ConstraintMatrix constraints;
+ constraints.close();
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+ // use this for info on problem
+ //std::cout << "Number of cells: " << dof.get_tria().n_active_cells()
+ // << std::endl;
+ //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+ //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+ Vector<double> solution (dof.n_dofs());
+
+ // create vector with random entries
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ if(constraints.is_constrained(i))
+ continue;
+ const double entry = rand()/(double)RAND_MAX;
+ solution(i) = entry;
+ }
+
+ constraints.distribute(solution);
+
+ FESystem<dim> fe_sys(dof.get_fe(), dim);
+ DoFHandler<dim> dofh_eulerian(dof.get_tria());
+ dofh_eulerian.distribute_dofs(fe_sys);
+
+ MatrixFree<dim,double> mf_data;
+ Vector<double> shift(dofh_eulerian.n_dofs());
+ for (unsigned int i=0; i<2; ++i)
+ {
+ if (i == 1)
+ {
+ shift(0) = 0.121;
+ shift(1) = -0.005;
+ }
+ MappingQEulerian<dim> mapping(2,shift,dofh_eulerian);
+
+ {
+ const QGauss<1> quad (fe_degree+1);
+ typename MatrixFree<dim,double>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::none;
+ data.mapping_update_flags = update_gradients | update_second_derivatives;
+ if (i==1)
+ data.initialize_indices = false;
+ mf_data.reinit (mapping, dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim,fe_degree,fe_degree+1> mf (mf_data, mapping);
+ mf.test_functions(solution);
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
projects / dealii-svn.git / commitdiff