18 #include <deal.II/base/utilities.h>
19 #include <deal.II/base/memory_consumption.h>
20 #include <deal.II/fe/fe_nothing.h>
21 #include <deal.II/fe/fe_values.h>
22 #include <deal.II/fe/mapping_q1.h>
24 #include <deal.II/matrix_free/mapping_info.h>
32 namespace MatrixFreeFunctions
36 template <
int dim,
typename Number>
39 JxW_values_initialized (false),
40 second_derivatives_initialized (false),
41 quadrature_points_initialized (false)
46 template <
int dim,
typename Number>
50 JxW_values_initialized =
false;
51 quadrature_points_initialized =
false;
52 second_derivatives_initialized =
false;
53 mapping_data_gen.clear();
55 cartesian_data.clear();
61 template <
int dim,
typename Number>
92 bool formula_with_one_point =
false;
93 for (
unsigned int i=0; i<quad[0].size(); ++i)
94 if (quad[0][i].size() == 1)
96 formula_with_one_point =
true;
99 if (formula_with_one_point ==
true)
100 for (
unsigned int comp=1; comp<quad.size(); ++comp)
101 for (
unsigned int i=0; i<quad[comp].size(); ++i)
102 if (quad[comp][i].size() > 1)
116 double get_jacobian_size (const ::Triangulation<dim> &tria)
118 if (tria.n_cells() == 0)
120 else return tria.begin()->diameter();
126 template <
int dim,
typename Number>
129 (const ::Triangulation<dim> &tria,
130 const std::vector<std::pair<unsigned int,unsigned int> > &cells,
131 const std::vector<unsigned int> &active_fe_index,
137 const unsigned int n_quads = quad.size();
138 const unsigned int n_cells = cells.size();
139 const unsigned int vectorization_length =
142 const unsigned int n_macro_cells = n_cells/vectorization_length;
143 mapping_data_gen.resize (n_quads);
144 cell_type.resize (n_macro_cells);
150 UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
153 JxW_values_initialized =
true;
155 second_derivatives_initialized =
true;
157 quadrature_points_initialized =
true;
166 const double jacobian_size = internal::get_jacobian_size(tria);
174 for (
unsigned int my_q=0; my_q<n_quads; ++my_q)
177 const unsigned int n_hp_quads = quad[my_q].size();
182 std::vector<unsigned int> n_q_points_1d (n_hp_quads),
183 step_size_cartesian (n_hp_quads);
186 for (
unsigned int q=0; q<n_hp_quads; ++q)
188 n_q_points_1d[q] = quad[my_q][q].size();
189 const unsigned int n_q_points =
190 Utilities::fixed_power<dim>(n_q_points_1d[q]);
191 current_data.
n_q_points.push_back (n_q_points);
194 (Utilities::fixed_power<dim-1>(n_q_points_1d[q]));
202 for (
unsigned int i=0; i<n_q_points; ++i)
213 step_size_cartesian[q] = 0;
214 unsigned int factor = 1;
215 for (
unsigned int d=0; d<dim; ++d)
217 step_size_cartesian[q] += factor;
218 factor *= n_q_points_1d[q];
223 if (cells.size() == 0)
232 CellType cell_t [vectorization_length],
233 cell_t_prev [vectorization_length];
234 for (
unsigned int j=0; j<vectorization_length; ++j)
235 cell_t_prev[j] = undefined;
243 std::vector<std_cxx1x::shared_ptr<FEValues<dim> > >
247 (update_flags & update_jacobian_grads ? update_jacobian_grads :
update_default) |
248 (update_flags & update_quadrature_points ? update_quadrature_points :
update_default);
255 const unsigned int n_cells_var =
256 mapping_data_gen[0].rowstart_jacobians.size()-1;
258 const unsigned int reserve_size = n_cells_var *
260 if (mapping_data_gen[0].jacobians.size() > 0)
262 if (mapping_data_gen[0].JxW_values.size() > 0)
264 if (mapping_data_gen[0].jacobians_grad_diag.size() > 0)
266 if (mapping_data_gen[0].jacobians_grad_upper.size() > 0)
276 std::map<Tensor<1,dim,VEC_ARRAY>,
unsigned int,
278 std::map<Tensor<2,dim,VEC_ARRAY>,
unsigned int,
282 for (
unsigned int cell=0; cell<n_macro_cells; ++cell)
288 const unsigned int fe_index = active_fe_index.size() > 0 ?
289 active_fe_index[cell] : 0;
290 const unsigned int n_q_points = current_data.
n_q_points[fe_index];
291 if (fe_values[fe_index].
get() == 0)
292 fe_values[fe_index].reset
295 update_flags_feval));
297 data.resize (n_q_points);
302 if (cell > 0 && active_fe_index.size() > 0 &&
303 active_fe_index[cell] != active_fe_index[cell-1])
304 cell_t_prev[vectorization_length-1] = undefined;
305 evaluate_on_cell (tria, &cells[cell*vectorization_length],
306 cell, my_q, cell_t_prev, cell_t, fe_val, data);
319 CellType most_general_type = cartesian;
320 for (
unsigned int j=0; j<vectorization_length; ++j)
321 if (cell_t[j] > most_general_type)
322 most_general_type = cell_t[j];
335 if (most_general_type == cartesian)
337 std::pair<Tensor<1,dim,VEC_ARRAY>,
unsigned int> new_entry;
338 new_entry.second = cartesians.size();
339 for (
unsigned int d=0; d<dim; ++d)
340 for (
unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
341 new_entry.first[d][v] = data.const_jac[d][d][v];
345 bool> it = cartesians.insert(new_entry);
346 if (it.second ==
false)
347 insert_position = it.first->second;
349 insert_position = new_entry.second;
354 else if (most_general_type == affine)
356 std::pair<Tensor<2,dim,VEC_ARRAY>,
unsigned int> new_entry;
357 new_entry.second = affines.size();
358 for (
unsigned int d=0; d<dim; ++d)
359 for (
unsigned int e=0; e<dim; ++e)
360 for (
unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
361 new_entry.first[d][e][v] = data.const_jac[d][e][v];
365 bool> it = affines.insert(new_entry);
366 if (it.second ==
false)
367 insert_position = it.first->second;
369 insert_position = new_entry.second;
384 unsigned int reserve_size = (n_macro_cells-cell+1)/2;
387 reserve_size *= n_q_points;
389 if (update_flags & update_JxW_values)
391 if (update_flags & update_jacobian_grads)
393 if (update_flags & update_jacobian_grads)
398 cell_type[cell] = ((insert_position << n_cell_type_bits) +
399 (
unsigned int)most_general_type);
406 if (get_cell_type(cell) == general)
408 const unsigned int previous_size =
411 if (update_flags & update_JxW_values)
416 if (update_flags & update_jacobian_grads)
423 for (
unsigned int q=0; q<n_q_points; ++q)
427 for (
unsigned int j=0; j<vectorization_length; ++j)
428 if (cell_t[j] == cartesian || cell_t[j] == affine)
430 for (
unsigned int d=0; d<dim; ++d)
431 for (
unsigned int e=0; e<dim; ++e)
433 jac[d][e][j] = data.const_jac[d][e][j];
434 for (
unsigned int f=0; f<dim; ++f)
435 jacobian_grad[d][e][f][j] = 0.;
443 if (update_flags & update_JxW_values)
447 if (update_flags & update_jacobian_grads)
466 for (
unsigned int d=0; d<dim; ++d)
467 for (
unsigned int e=0; e<dim; ++e)
468 for (
unsigned int f=0; f<dim; ++f)
470 jac_grad[f][e][d] = (inv_jac[f][0] *
471 jacobian_grad[d][e][0]);
472 for (
unsigned int g=1; g<dim; ++g)
473 jac_grad[f][e][d] += (inv_jac[f][g] *
474 jacobian_grad[d][e][g]);
478 for (
unsigned int d=0; d<dim; ++d)
479 for (
unsigned int e=0; e<dim; ++e)
481 for (
unsigned int f=0; f<dim; ++f)
484 for (
unsigned int g=0; g<dim; ++g)
485 tmp[f] -= jac_grad[d][f][g] * inv_jac[g][e];
490 for (
unsigned int f=0; f<dim; ++f)
492 grad_jac_inv[f][d][e] = inv_jac[f][0] * tmp[0];
493 for (
unsigned int g=1; g<dim; ++g)
494 grad_jac_inv[f][d][e] += inv_jac[f][g] * tmp[g];
500 for (
unsigned int d=0; d<dim; ++d)
501 for (
unsigned int e=0; e<dim; ++e)
502 grad_diag[d][e] = grad_jac_inv[d][d][e];
509 for (
unsigned int d=0, count=0; d<dim; ++d)
510 for (
unsigned int e=d+1; e<dim; ++e, ++count)
511 for (
unsigned int f=0; f<dim; ++f)
512 grad_upper[count][f] = grad_jac_inv[d][e][f];
518 if (update_flags & update_quadrature_points)
530 const unsigned int old_size =
536 if (get_cell_type(cell) == cartesian)
539 n_q_points_1d[fe_index]);
540 for (
unsigned int q=0; q<n_q_points_1d[fe_index]; ++q)
541 for (
unsigned int d=0; d<dim; ++d)
543 data.quadrature_points[q*step_size_cartesian[fe_index]][d];
548 for (
unsigned int q=0; q<n_q_points; ++q)
549 for (
unsigned int d=0; d<dim; ++d)
551 data.quadrature_points[q][d];
565 cartesian_data.resize(cartesians.size());
568 it = cartesians.begin(); it != cartesians.end(); ++it)
571 for (
unsigned int d=0; d<dim; ++d)
574 for (
unsigned int v=0;
575 v<VectorizedArray<Number>::n_array_elements; ++v)
576 jac_d[v] = it->first[d][v];
577 cartesian_data[it->second].first[d] = 1./jac_d;
580 cartesian_data[it->second].second = det;
582 affine_data.resize(affines.size());
585 it = affines.begin(); it != affines.end(); ++it)
588 for (
unsigned int d=0; d<dim; ++d)
589 for (
unsigned int e=0; e<dim; ++e)
590 for (
unsigned int v=0;
591 v<VectorizedArray<Number>::n_array_elements; ++v)
592 jac[d][e][v] = it->first[d][e][v];
594 affine_data[it->second].second = determinant(jac);
595 affine_data[it->second].first = transpose(invert(jac));
603 template<
int dim,
typename Number>
606 const std::pair<unsigned int,unsigned int> *cells,
607 const unsigned int cell,
608 const unsigned int my_q,
610 CellType (&cell_t)[VectorizedArray<Number>::n_array_elements],
614 const unsigned int vectorization_length =
615 VectorizedArray<Number>::n_array_elements;
621 const double zero_tolerance_double = data.jac_size *
622 std::numeric_limits<double>::epsilon() * 1024.;
623 for (
unsigned int j=0; j<vectorization_length; ++j)
625 typename ::Triangulation<dim>::cell_iterator
626 cell_it (&tria, cells[j].first, cells[j].second);
628 cell_t[j] = undefined;
633 for (
unsigned int q=0; q<n_q_points; ++q)
636 for (
unsigned int d=0; d<dim; ++d)
637 data.quadrature_points[q][d][j] = point[d];
643 if (my_q > 0 && (get_cell_type(cell) == cartesian
644 || get_cell_type(cell) == affine) )
656 cell_t[j] = cell_t_prev[vectorization_length-1];
659 cell_t[j] = cell_t[j-1];
668 if (cell_t[j] == undefined)
670 bool jacobian_constant =
true;
671 for (
unsigned int q=1; q<n_q_points; ++q)
674 for (
unsigned int d=0; d<dim; ++d)
675 for (
unsigned int e=0; e<dim; ++e)
676 if (std::fabs(jac_0[d][e]-jac[d][e]) >
677 zero_tolerance_double)
678 jacobian_constant =
false;
679 if (jacobian_constant ==
false)
685 bool cell_cartesian = jacobian_constant;
686 for (
unsigned int d=0; d<dim; ++d)
687 for (
unsigned int e=0; e<dim; ++e)
689 if (std::fabs(jac_0[d][e]) >
690 zero_tolerance_double)
692 cell_cartesian=
false;
705 if (cell_cartesian ==
false && n_q_points == 1 &&
711 for (
unsigned int d=0; d<dim; ++d)
712 for (
unsigned int e=0; e<dim; ++e)
713 for (
unsigned int f=0; f<dim; ++f)
715 double jac_grad_comp = (jac[f][0] *
716 jacobian_grad[d][e][0]);
717 for (
unsigned int g=1; g<dim; ++g)
718 jac_grad_comp += (jac[f][g] *
719 jacobian_grad[d][e][g]);
720 if (std::fabs(jac_grad_comp) >
721 zero_tolerance_double)
722 jacobian_constant =
false;
726 if (cell_cartesian ==
true)
727 cell_t[j] = cartesian;
728 else if (jacobian_constant ==
true)
735 if (cell_t[j] == cartesian)
738 for (
unsigned int d=0; d<dim; ++d)
740 data.const_jac[d][d][j] = jac_0[d][d];
741 for (
unsigned int e=d+1; e<dim; ++e)
743 data.const_jac[d][e][j] = 0.;
744 data.const_jac[e][d][j] = 0.;
751 else if (cell_t[j] == affine)
754 for (
unsigned int d=0; d<dim; ++d)
755 for (
unsigned int e=0; e<dim; ++e)
756 data.const_jac[d][e][j] =
757 std::fabs(jac_0[d][e]) < zero_tolerance_double ?
768 for (
unsigned int q=0; q<n_q_points; ++q)
775 for (
unsigned int d=0; d<dim; ++d)
776 for (
unsigned int e=0; e<dim; ++e)
777 data.general_jac[q][d][e][j] =
778 std::fabs(jac[d][e]) < zero_tolerance_double ? 0. : jac[d][e];
788 for (
unsigned int d=0; d<dim; ++d)
789 for (
unsigned int e=0; e<dim; ++e)
790 for (
unsigned int f=0; f<dim; ++f)
791 data.general_jac_grad[q][d][e][f][j] = jacobian_grad[d][e][f];
798 for (
unsigned int j=0; j<vectorization_length; ++j)
799 cell_t_prev[j] = cell_t[j];
803 template <
int dim,
typename Number>
806 jac_size (jac_size_in)
811 template <
int dim,
typename Number>
813 MappingInfo<dim,Number>::CellData::resize (
const unsigned int size)
815 if (general_jac.
size() != size)
817 quadrature_points.
resize(size);
819 general_jac_grad.
resize(size);
825 template <
int dim,
typename Number>
846 template <
int dim,
typename Number>
854 memory +=
sizeof (
this);
860 template <
int dim,
typename Number>
861 template <
typename STREAM>
868 #ifdef DEAL_II_WITH_MPI
869 unsigned int general_size_glob = 0, general_size_loc = jacobians.size();
870 MPI_Allreduce (&general_size_loc, &general_size_glob, 1, MPI_UNSIGNED,
873 unsigned int general_size_glob = jacobians.size();
875 if (general_size_glob > 0)
877 out <<
" Memory Jacobian data: ";
881 out <<
" Memory second derivative data: ";
887 #ifdef DEAL_II_WITH_MPI
888 unsigned int quad_size_glob = 0, quad_size_loc = quadrature_points.size();
889 MPI_Allreduce (&quad_size_loc, &quad_size_glob, 1, MPI_UNSIGNED,
892 unsigned int quad_size_glob = quadrature_points.size();
894 if (quad_size_glob > 0)
896 out <<
" Memory quadrature points: ";
905 template <
int dim,
typename Number>
906 template <
typename STREAM>
910 out <<
" Cell types: ";
913 out <<
" Memory transformations compr: ";
919 out <<
" Data component " << j << std::endl;
928 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
static const unsigned int invalid_unsigned_int
void push_back(const Quadrature< dim > &new_quadrature)
#define AssertDimension(dim1, dim2)
AlignedVector< std::pair< Tensor< 1, dim, VectorizedArray< Number > >, VectorizedArray< Number > > > cartesian_data
::hp::QCollection< dim > quadrature
AlignedVector< VectorizedArray< Number > > JxW_values
AlignedVector< std::pair< Tensor< 2, dim, VectorizedArray< Number > >, VectorizedArray< Number > > > affine_data
const DerivativeForm< 1, dim, spacedim > & jacobian(const unsigned int quadrature_point) const
std::vector< unsigned int > n_q_points
#define AssertIndexRange(index, range)
void reserve(const size_type size_alloc)
Transformed quadrature points.
std::vector< MappingInfoDependent > mapping_data_gen
std::vector< unsigned int > rowstart_jacobians
std::vector< unsigned int > n_q_points_face
void push_back(const T in_data)
CellSimilarity::Similarity get_cell_similarity() const
void print_memory_consumption(STREAM &out, const SizeInfo &size_info) const
const DerivativeForm< 2, dim, spacedim > & jacobian_grad(const unsigned int quadrature_point) const
UpdateFlags compute_update_flags(const UpdateFlags update_flags, const std::vector<::hp::QCollection< 1 > > &quad) const
void print_memory_statistics(STREAM &out, std::size_t data_length) const
std::vector< unsigned int > cell_type
#define Assert(cond, exc)
std::vector< unsigned int > quad_index_conversion
std::size_t memory_consumption(const T &t)
AlignedVector< Tensor< 1,(dim >1?dim *(dim-1)/2:1), Tensor< 1, dim, VectorizedArray< Number > > > > jacobians_grad_upper
void reinit(const TriaIterator< DoFCellAccessor< DH, level_dof_access > > cell)
void print_memory_consumption(STREAM &out, const SizeInfo &size_info) const
void resize(const size_type size_in, const T &init=T())
std::vector< unsigned int > rowstart_q_points
const Point< spacedim > & quadrature_point(const unsigned int i) const
std::size_t memory_consumption() const
Second derivatives of shape functions.
Gradient of volume element.
std::size_t memory_consumption() const
::hp::QCollection< dim-1 > face_quadrature
AlignedVector< Point< dim, VectorizedArray< Number > > > quadrature_points
Shape function gradients.
void initialize(const ::Triangulation< dim > &tria, const std::vector< std::pair< unsigned int, unsigned int > > &cells, const std::vector< unsigned int > &active_fe_index, const Mapping< dim > &mapping, const std::vector<::hp::QCollection< 1 > > &quad, const UpdateFlags update_flags)
UpdateFlags get_update_flags() const
::ExceptionBase & ExcInternalError()
std::vector< AlignedVector< VectorizedArray< Number > > > quadrature_weights
const unsigned int n_quadrature_points
void evaluate_on_cell(const ::Triangulation< dim > &tria, const std::pair< unsigned int, unsigned int > *cells, const unsigned int cell, const unsigned int my_q, CellType(&cell_t_prev)[VectorizedArray< Number >::n_array_elements], CellType(&cell_t)[VectorizedArray< Number >::n_array_elements], FEValues< dim, dim > &fe_values, CellData &cell_data) const
AlignedVector< Tensor< 2, dim, VectorizedArray< Number > > > jacobians
AlignedVector< Tensor< 2, dim, VectorizedArray< Number > > > jacobians_grad_diag