Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
mapping_info.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: mapping_info.templates.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2011 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
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>
23 
24 #include <deal.II/matrix_free/mapping_info.h>
25 
26 
28 
29 
30 namespace internal
31 {
32  namespace MatrixFreeFunctions
33  {
34  // ----------------- actual MappingInfo functions -------------------------
35 
36  template <int dim, typename Number>
38  :
39  JxW_values_initialized (false),
40  second_derivatives_initialized (false),
41  quadrature_points_initialized (false)
42  {}
43 
44 
45 
46  template <int dim, typename Number>
47  void
49  {
50  JxW_values_initialized = false;
51  quadrature_points_initialized = false;
52  second_derivatives_initialized = false;
53  mapping_data_gen.clear();
54  cell_type.clear();
55  cartesian_data.clear();
56  affine_data.clear();
57  }
58 
59 
60 
61  template <int dim, typename Number>
64  compute_update_flags (const UpdateFlags update_flags,
65  const std::vector<::hp::QCollection<1> > &quad) const
66  {
67  // this class is build around the evaluation this class is build around
68  // the evaluation of inverse gradients, so compute them in any case
70 
71  // if the user requested gradients, need inverse Jacobians
72  if (update_flags & update_gradients || update_flags & update_inverse_jacobians)
73  new_flags |= update_inverse_jacobians;
74 
75  // for JxW, would only need JxW values.
76  if (update_flags & update_JxW_values)
77  new_flags |= update_JxW_values;
78 
79  // for Hessian information, need inverse Jacobians and the derivative of
80  // Jacobians (these two together will give use the gradients of the
81  // inverse Jacobians, which is what we need)
82  if (update_flags & update_hessians || update_flags & update_jacobian_grads)
83  new_flags |= update_jacobian_grads;
84 
85  if (update_flags & update_quadrature_points)
86  new_flags |= update_quadrature_points;
87 
88  // there is one more thing: if we have a quadrature formula with only
89  // one quadrature point on the first component, but more points on later
90  // components, we need to have Jacobian gradients anyway in order to
91  // determine whether the Jacobian is constant throughout a cell
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)
95  {
96  formula_with_one_point = true;
97  break;
98  }
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)
103  {
104  new_flags |= update_jacobian_grads;
105  goto end_set;
106  }
107 end_set:
108  return new_flags;
109  }
110 
111 
112 
113  namespace internal
114  {
115  template <int dim>
116  double get_jacobian_size (const ::Triangulation<dim> &tria)
117  {
118  if (tria.n_cells() == 0)
119  return 1;
120  else return tria.begin()->diameter();
121  }
122  }
123 
124 
125 
126  template <int dim, typename Number>
127  void
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,
132  const Mapping<dim> &mapping,
133  const std::vector<::hp::QCollection<1> > &quad,
134  const UpdateFlags update_flags_input)
135  {
136  clear();
137  const unsigned int n_quads = quad.size();
138  const unsigned int n_cells = cells.size();
139  const unsigned int vectorization_length =
141  Assert (n_cells%vectorization_length == 0, ExcInternalError());
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);
145 
146  // dummy FE that is used to set up an FEValues object. Do not need the
147  // actual finite element because we will only evaluate quantities for
148  // the mapping that are independent of the FE
149  FE_Nothing<dim> dummy_fe;
150  UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
151 
152  if (update_flags & update_JxW_values)
153  JxW_values_initialized = true;
154  if (update_flags & update_jacobian_grads)
155  second_derivatives_initialized = true;
156  if (update_flags & update_quadrature_points)
157  quadrature_points_initialized = true;
158 
159  // when we make comparisons about the size of Jacobians we need to know
160  // the approximate size of typical entries in Jacobians. We need to fix
161  // the Jacobian size once and for all. We choose the diameter of the
162  // first cell (on level zero, which is the best accuracy we can hope
163  // for, since diameters on finer levels are computed by differences of
164  // nearby cells). If the mesh extends over a certain domain, the
165  // precision of double values is essentially limited by this precision.
166  const double jacobian_size = internal::get_jacobian_size(tria);
167 
168  // objects that hold the data for up to vectorization_length cells while
169  // we fill them up. Only after all vectorization_length cells have been
170  // processed, we can insert the data into the data structures of this
171  // class
172  CellData data (jacobian_size);
173 
174  for (unsigned int my_q=0; my_q<n_quads; ++my_q)
175  {
176  MappingInfoDependent &current_data = mapping_data_gen[my_q];
177  const unsigned int n_hp_quads = quad[my_q].size();
178  AssertIndexRange (0, n_hp_quads);
179  current_data.n_q_points.reserve (n_hp_quads);
180  current_data.n_q_points_face.reserve (n_hp_quads);
181  current_data.quadrature_weights.resize (n_hp_quads);
182  std::vector<unsigned int> n_q_points_1d (n_hp_quads),
183  step_size_cartesian (n_hp_quads);
184  if (n_hp_quads > 1)
185  current_data.quad_index_conversion.resize(n_hp_quads);
186  for (unsigned int q=0; q<n_hp_quads; ++q)
187  {
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);
192 
193  current_data.n_q_points_face.push_back
194  (Utilities::fixed_power<dim-1>(n_q_points_1d[q]));
195  current_data.quadrature.push_back
196  (Quadrature<dim>(quad[my_q][q]));
197  current_data.face_quadrature.push_back
198  (Quadrature<dim-1>(quad[my_q][q]));
199 
200  // set quadrature weights in vectorized form
201  current_data.quadrature_weights[q].resize(n_q_points);
202  for (unsigned int i=0; i<n_q_points; ++i)
203  current_data.quadrature_weights[q][i] =
204  current_data.quadrature[q].get_weights()[i];
205 
206  if (n_hp_quads > 1)
207  current_data.quad_index_conversion[q] = n_q_points;
208 
209  // To walk on the diagonal for lexicographic ordering, we have
210  // to jump one index ahead in each direction. For direction 0,
211  // this is just the next point, for direction 1, it means adding
212  // n_q_points_1d, and so on.
213  step_size_cartesian[q] = 0;
214  unsigned int factor = 1;
215  for (unsigned int d=0; d<dim; ++d)
216  {
217  step_size_cartesian[q] += factor;
218  factor *= n_q_points_1d[q];
219  }
220  }
221 
222  // if there are no cells, there is nothing to do
223  if (cells.size() == 0)
224  continue;
225 
226  Tensor<3,dim,VectorizedArray<Number> > jac_grad, grad_jac_inv;
228 
229  // encodes the cell types of the current cell. Since several cells
230  // must be considered together, this variable holds the individual
231  // info of the last chunk of cells
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;
236 
237  // fe_values object that is used to compute the mapping data. for
238  // the hp case there might be more than one finite element. since we
239  // manually select the active FE index and not via a
240  // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
241  // select the correct finite element, so just hold a vector of
242  // FEValues
243  std::vector<std_cxx1x::shared_ptr<FEValues<dim> > >
244  fe_values (current_data.quadrature.size());
245  UpdateFlags update_flags_feval =
247  (update_flags & update_jacobian_grads ? update_jacobian_grads : update_default) |
248  (update_flags & update_quadrature_points ? update_quadrature_points : update_default);
249 
250  // resize the fields that have fixed size or for which we know
251  // something from an earlier loop
252  current_data.rowstart_q_points.resize (n_macro_cells+1);
253  if (my_q > 0)
254  {
255  const unsigned int n_cells_var =
256  mapping_data_gen[0].rowstart_jacobians.size()-1;
257  current_data.rowstart_jacobians.reserve (n_cells_var+1);
258  const unsigned int reserve_size = n_cells_var *
259  current_data.n_q_points[0];
260  if (mapping_data_gen[0].jacobians.size() > 0)
261  current_data.jacobians.reserve (reserve_size);
262  if (mapping_data_gen[0].JxW_values.size() > 0)
263  current_data.jacobians.reserve (reserve_size);
264  if (mapping_data_gen[0].jacobians_grad_diag.size() > 0)
265  current_data.jacobians_grad_diag.reserve (reserve_size);
266  if (mapping_data_gen[0].jacobians_grad_upper.size() > 0)
267  current_data.jacobians_grad_upper.reserve (reserve_size);
268  }
269 
270  // we would like to put a Tensor<1,dim,VectorizedArray<Number> > as
271  // key into the std::map, but std::map allocation does not align the
272  // allocated memory correctly, so put it into a tensor of the
273  // correct length instead
274  FPArrayComparator<Number> comparator(jacobian_size);
275  typedef Tensor<1,VectorizedArray<Number>::n_array_elements,Number> VEC_ARRAY;
276  std::map<Tensor<1,dim,VEC_ARRAY>, unsigned int,
277  FPArrayComparator<Number> > cartesians(comparator);
278  std::map<Tensor<2,dim,VEC_ARRAY>, unsigned int,
279  FPArrayComparator<Number> > affines(comparator);
280 
281  // loop over all cells
282  for (unsigned int cell=0; cell<n_macro_cells; ++cell)
283  {
284  // GENERAL OUTLINE: First generate the data in format "number"
285  // for vectorization_length cells, and then find the most
286  // general type of cell for appropriate vectorized formats. then
287  // fill this data in
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
293  (new FEValues<dim> (mapping, dummy_fe,
294  current_data.quadrature[fe_index],
295  update_flags_feval));
296  FEValues<dim> &fe_val = *fe_values[fe_index];
297  data.resize (n_q_points);
298 
299  // if the fe index has changed from the previous cell, set the
300  // old cell type to invalid (otherwise, we might detect
301  // similarity due to some cells further ahead)
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);
307 
308  // now reorder the data into vectorized types. if we are here
309  // for the first time, we need to find out whether the Jacobian
310  // allows for some simplification (Cartesian, affine) taking
311  // vectorization_length cell together and we have to insert that
312  // data into the respective fields. Also, we have to compress
313  // different cell indicators into one structure.
314 
315  if (my_q == 0)
316  {
317  // find the most general cell type (most general type is 2
318  // (general cell))
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];
323  AssertIndexRange (most_general_type, 3);
324  unsigned int insert_position = numbers::invalid_unsigned_int;
325 
326  // Cartesian cell with diagonal Jacobian: only insert the
327  // diagonal of the inverse and the Jacobian determinant. We
328  // do this by using an std::map that collects pointers to
329  // all Cartesian Jacobians. We need a pointer in the
330  // std::map because it cannot store data based on
331  // VectorizedArray (alignment issue). We circumvent the
332  // problem by temporarily filling the next value into the
333  // cartesian_data field and, in case we did an insertion,
334  // the data is already in the correct place.
335  if (most_general_type == cartesian)
336  {
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];
342 
343  std::pair<typename std::map<Tensor<1,dim,VEC_ARRAY>,
344  unsigned int, FPArrayComparator<Number> >::iterator,
345  bool> it = cartesians.insert(new_entry);
346  if (it.second == false)
347  insert_position = it.first->second;
348  else
349  insert_position = new_entry.second;
350  }
351 
352  // Constant Jacobian case. same strategy as before, but with
353  // other data fields
354  else if (most_general_type == affine)
355  {
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];
362 
363  std::pair<typename std::map<Tensor<2,dim,VEC_ARRAY>,
364  unsigned int, FPArrayComparator<Number> >::iterator,
365  bool> it = affines.insert(new_entry);
366  if (it.second == false)
367  insert_position = it.first->second;
368  else
369  insert_position = new_entry.second;
370  }
371 
372  // general cell case: first resize the data field to fit the
373  // new data. if we are here the first time, assume that
374  // there are many general cells to come, so reserve some
375  // memory in order to not have too many reallocations and
376  // memcpy's. The scheme used here involves at most one
377  // reallocation.
378  else
379  {
380  Assert (most_general_type == general, ExcInternalError());
381  insert_position = current_data.rowstart_jacobians.size();
382  if (current_data.rowstart_jacobians.size() == 0)
383  {
384  unsigned int reserve_size = (n_macro_cells-cell+1)/2;
385  current_data.rowstart_jacobians.reserve
386  (reserve_size);
387  reserve_size *= n_q_points;
388  current_data.jacobians.reserve (reserve_size);
389  if (update_flags & update_JxW_values)
390  current_data.JxW_values.reserve (reserve_size);
391  if (update_flags & update_jacobian_grads)
392  current_data.jacobians_grad_diag.reserve (reserve_size);
393  if (update_flags & update_jacobian_grads)
394  current_data.jacobians_grad_upper.reserve (reserve_size);
395  }
396  }
397 
398  cell_type[cell] = ((insert_position << n_cell_type_bits) +
399  (unsigned int)most_general_type);
400 
401  } // end if (my_q == 0)
402 
403  // general cell case: now go through all quadrature points and
404  // collect the data. done for all different quadrature formulas,
405  // so do it outside the above loop.
406  if (get_cell_type(cell) == general)
407  {
408  const unsigned int previous_size =
409  current_data.jacobians.size();
410  current_data.rowstart_jacobians.push_back (previous_size);
411  if (update_flags & update_JxW_values)
412  {
413  AssertDimension (previous_size,
414  current_data.JxW_values.size());
415  }
416  if (update_flags & update_jacobian_grads)
417  {
418  AssertDimension (previous_size,
419  current_data.jacobians_grad_diag.size());
420  AssertDimension (previous_size,
421  current_data.jacobians_grad_upper.size());
422  }
423  for (unsigned int q=0; q<n_q_points; ++q)
424  {
425  Tensor<2,dim,VectorizedArray<Number> > &jac = data.general_jac[q];
426  Tensor<3,dim,VectorizedArray<Number> > &jacobian_grad = data.general_jac_grad[q];
427  for (unsigned int j=0; j<vectorization_length; ++j)
428  if (cell_t[j] == cartesian || cell_t[j] == affine)
429  {
430  for (unsigned int d=0; d<dim; ++d)
431  for (unsigned int e=0; e<dim; ++e)
432  {
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.;
436  }
437  }
438 
439  const VectorizedArray<Number> det = determinant (jac);
440  current_data.jacobians.push_back (transpose(invert(jac)));
441  const Tensor<2,dim,VectorizedArray<Number> > &inv_jac = current_data.jacobians.back();
442 
443  if (update_flags & update_JxW_values)
444  current_data.JxW_values.push_back
445  (det * current_data.quadrature_weights[fe_index][q]);
446 
447  if (update_flags & update_jacobian_grads)
448  {
449  // for second derivatives on the real cell, need
450  // also the gradient of the inverse Jacobian J. This
451  // involves some calculus and is done
452  // vectorized. This is very cheap compared to what
453  // fe_values does (in early 2011). If L is the
454  // gradient of the jacobian on the unit cell, the
455  // gradient of the inverse is given by
456  // (multidimensional calculus) - J * (J * L) * J
457  // (the third J is because we need to transform the
458  // gradient L from the unit to the real cell, and
459  // then apply the inverse Jacobian). Compare this
460  // with 1D with j(x) = 1/k(phi(x)), where j = phi'
461  // is the inverse of the jacobian and k is the
462  // derivative of the jacobian on the unit cell. Then
463  // j' = phi' k'/k^2 = j k' j^2.
464 
465  // compute: jac_grad = J*grad_unit(J^-1)
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)
469  {
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]);
475  }
476 
477  // compute: transpose (-jac * jac_grad[d] * jac)
478  for (unsigned int d=0; d<dim; ++d)
479  for (unsigned int e=0; e<dim; ++e)
480  {
481  for (unsigned int f=0; f<dim; ++f)
482  {
483  tmp[f] = VectorizedArray<Number>();
484  for (unsigned int g=0; g<dim; ++g)
485  tmp[f] -= jac_grad[d][f][g] * inv_jac[g][e];
486  }
487 
488  // needed for non-diagonal part of Jacobian
489  // grad
490  for (unsigned int f=0; f<dim; ++f)
491  {
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];
495  }
496  }
497 
498  {
499  VectorizedArray<Number> grad_diag[dim][dim];
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];
503  current_data.jacobians_grad_diag.push_back
504  (Tensor<2,dim,VectorizedArray<Number> >(grad_diag));
505  }
506 
507  // sets upper-diagonal part of Jacobian
508  Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > grad_upper;
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];
513  current_data.jacobians_grad_upper.push_back(grad_upper);
514  }
515  }
516  }
517 
518  if (update_flags & update_quadrature_points)
519  {
520  // eventually we turn to the quadrature points that we can
521  // compress in case we have Cartesian cells. we also need to
522  // reorder them into arrays of vectorized data types. first
523  // go through the cells and find out how much memory we need
524  // to allocate for the quadrature points. We store
525  // n_q_points for all cells but Cartesian cells. For
526  // Cartesian cells, only need to store the values on a
527  // diagonal through the cell (n_q_points_1d). This will give
528  // (with some little indexing) the location of all
529  // quadrature points.
530  const unsigned int old_size =
531  current_data.quadrature_points.size();
532  current_data.rowstart_q_points[cell] = old_size;
533 
535 
536  if (get_cell_type(cell) == cartesian)
537  {
538  current_data.quadrature_points.resize (old_size+
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)
542  current_data.quadrature_points[old_size+q][d] =
543  data.quadrature_points[q*step_size_cartesian[fe_index]][d];
544  }
545  else
546  {
547  current_data.quadrature_points.resize (old_size + n_q_points);
548  for (unsigned int q=0; q<n_q_points; ++q)
549  for (unsigned int d=0; d<dim; ++d)
550  current_data.quadrature_points[old_size+q][d] =
551  data.quadrature_points[q][d];
552  }
553  }
554  } // end for ( cell < n_macro_cells )
555  current_data.rowstart_jacobians.push_back
556  (current_data.jacobians.size());
557  current_data.rowstart_q_points[n_macro_cells] =
558  current_data.quadrature_points.size();
559 
560  // finally, fill the accumulated data for Cartesian and affine cells
561  // into cartesian_data and affine_data, invert and transpose the
562  // Jacobians, and compute the JxW value.
563  if (my_q == 0)
564  {
565  cartesian_data.resize(cartesians.size());
566  for (typename std::map<Tensor<1,dim,VEC_ARRAY>,
567  unsigned int, FPArrayComparator<Number> >::iterator
568  it = cartesians.begin(); it != cartesians.end(); ++it)
569  {
570  VectorizedArray<Number> det = make_vectorized_array<Number>(1.);
571  for (unsigned int d=0; d<dim; ++d)
572  {
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;
578  det *= jac_d;
579  }
580  cartesian_data[it->second].second = det;
581  }
582  affine_data.resize(affines.size());
583  for (typename std::map<Tensor<2,dim,VEC_ARRAY>,
584  unsigned int, FPArrayComparator<Number> >::iterator
585  it = affines.begin(); it != affines.end(); ++it)
586  {
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];
593 
594  affine_data[it->second].second = determinant(jac);
595  affine_data[it->second].first = transpose(invert(jac));
596  }
597  }
598  }
599  }
600 
601 
602 
603  template<int dim, typename Number>
604  void
605  MappingInfo<dim,Number>::evaluate_on_cell (const ::Triangulation<dim> &tria,
606  const std::pair<unsigned int,unsigned int> *cells,
607  const unsigned int cell,
608  const unsigned int my_q,
609  CellType (&cell_t_prev)[VectorizedArray<Number>::n_array_elements],
610  CellType (&cell_t)[VectorizedArray<Number>::n_array_elements],
611  FEValues<dim,dim> &fe_val,
612  CellData &data) const
613  {
614  const unsigned int vectorization_length =
615  VectorizedArray<Number>::n_array_elements;
616  const unsigned int n_q_points = fe_val.n_quadrature_points;
617  const UpdateFlags update_flags = fe_val.get_update_flags();
618 
619  // this should be the same value as used in HashValue::scaling (but we
620  // not have that field here)
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)
624  {
625  typename ::Triangulation<dim>::cell_iterator
626  cell_it (&tria, cells[j].first, cells[j].second);
627  fe_val.reinit(cell_it);
628  cell_t[j] = undefined;
629 
630  // extract quadrature points and store them temporarily. if we have
631  // Cartesian cells, we can compress the indices
632  if (update_flags & update_quadrature_points)
633  for (unsigned int q=0; q<n_q_points; ++q)
634  {
635  const Point<dim> &point = fe_val.quadrature_point(q);
636  for (unsigned int d=0; d<dim; ++d)
637  data.quadrature_points[q][d][j] = point[d];
638  }
639 
640  // if this is not the first quadrature formula and we already have
641  // determined that this cell is either Cartesian or with constant
642  // Jacobian, we have nothing more to do.
643  if (my_q > 0 && (get_cell_type(cell) == cartesian
644  || get_cell_type(cell) == affine) )
645  continue;
646 
647  // first round: if the transformation is detected to be the same as
648  // on the old cell, we only need to copy over the data.
649  if (fe_val.get_cell_similarity() == CellSimilarity::translation
650  &&
651  my_q == 0)
652  {
653  if (j==0)
654  {
655  Assert (cell>0, ExcInternalError());
656  cell_t[j] = cell_t_prev[vectorization_length-1];
657  }
658  else
659  cell_t[j] = cell_t[j-1];
660  }
661 
662  const DerivativeForm<1,dim,dim> &jac_0 = fe_val.jacobian(0);
663 
664  if (my_q == 0)
665  {
666  // check whether the Jacobian is constant on this cell the first
667  // time we come around here
668  if (cell_t[j] == undefined)
669  {
670  bool jacobian_constant = true;
671  for (unsigned int q=1; q<n_q_points; ++q)
672  {
673  const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(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)
680  break;
681  }
682 
683  // check whether the Jacobian is diagonal to machine
684  // accuracy
685  bool cell_cartesian = jacobian_constant;
686  for (unsigned int d=0; d<dim; ++d)
687  for (unsigned int e=0; e<dim; ++e)
688  if (d!=e)
689  if (std::fabs(jac_0[d][e]) >
690  zero_tolerance_double)
691  {
692  cell_cartesian=false;
693  break;
694  }
695 
696  // in case we have only one quadrature point, we can have
697  // non-constant Jacobians, but we cannot detect it by
698  // comparison from one quadrature point to the next: in that
699  // case, need to look at second derivatives and see whether
700  // there are some non-zero entries (this is necessary since
701  // we determine the constness of the Jacobian for the first
702  // quadrature formula and might not look at them any more
703  // for the second, third quadrature formula). in any case,
704  // the flag update_jacobian_grads will be set in that case
705  if (cell_cartesian == false && n_q_points == 1 &&
706  update_flags & update_jacobian_grads)
707  {
708  const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(0);
709  const DerivativeForm<2,dim,dim> &jacobian_grad =
710  fe_val.jacobian_grad(0);
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)
714  {
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;
723  }
724  }
725  // set cell type
726  if (cell_cartesian == true)
727  cell_t[j] = cartesian;
728  else if (jacobian_constant == true)
729  cell_t[j] = affine;
730  else
731  cell_t[j] = general;
732  }
733 
734  // Cartesian cell
735  if (cell_t[j] == cartesian)
736  {
737  // set Jacobian into diagonal and clear off-diagonal part
738  for (unsigned int d=0; d<dim; ++d)
739  {
740  data.const_jac[d][d][j] = jac_0[d][d];
741  for (unsigned int e=d+1; e<dim; ++e)
742  {
743  data.const_jac[d][e][j] = 0.;
744  data.const_jac[e][d][j] = 0.;
745  }
746  }
747  continue;
748  }
749 
750  // cell with affine mapping
751  else if (cell_t[j] == affine)
752  {
753  // compress out very small values
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 ?
758  0 : jac_0[d][e];
759  continue;
760  }
761  }
762 
763  // general cell case
764 
765  // go through all quadrature points and fill in the data into the
766  // temporary data structures with slots for the vectorized data
767  // types
768  for (unsigned int q=0; q<n_q_points; ++q)
769  {
770 
771  // compress out very small numbers which are only noise. Then it
772  // is cleaner to use zero straight away (though it does not save
773  // any memory)
774  const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(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];
779 
780  // need to do some calculus based on the gradient of the
781  // Jacobian, in order to find the gradient of the inverse
782  // Jacobian which is needed in user code. however, we would like
783  // to perform that on vectorized data types instead of doubles
784  // or floats. to this end, copy the gradients first
785  if (update_flags & update_jacobian_grads)
786  {
787  const DerivativeForm<2,dim,dim> &jacobian_grad = fe_val.jacobian_grad(q);
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];
792  }
793  }
794  } // end loop over all entries in vectorization (vectorization_length
795  // cells)
796 
797  // set information for next cell
798  for (unsigned int j=0; j<vectorization_length; ++j)
799  cell_t_prev[j] = cell_t[j];
800  }
801 
802 
803  template <int dim, typename Number>
804  MappingInfo<dim,Number>::CellData::CellData (const double jac_size_in)
805  :
806  jac_size (jac_size_in)
807  {}
808 
809 
810 
811  template <int dim, typename Number>
812  void
813  MappingInfo<dim,Number>::CellData::resize (const unsigned int size)
814  {
815  if (general_jac.size() != size)
816  {
817  quadrature_points.resize(size);
818  general_jac.resize(size);
819  general_jac_grad.resize(size);
820  }
821  }
822 
823 
824 
825  template <int dim, typename Number>
827  {
828  std::size_t
829  memory = MemoryConsumption::memory_consumption (jacobians);
830  memory += MemoryConsumption::memory_consumption (JxW_values);
831  memory += MemoryConsumption::memory_consumption (jacobians_grad_diag);
832  memory += MemoryConsumption::memory_consumption (jacobians_grad_upper);
833  memory += MemoryConsumption::memory_consumption (rowstart_q_points);
834  memory += MemoryConsumption::memory_consumption (quadrature_points);
835  memory += MemoryConsumption::memory_consumption (quadrature);
836  memory += MemoryConsumption::memory_consumption (face_quadrature);
837  memory += MemoryConsumption::memory_consumption (quadrature_weights);
838  memory += MemoryConsumption::memory_consumption (n_q_points);
839  memory += MemoryConsumption::memory_consumption (n_q_points_face);
840  memory += MemoryConsumption::memory_consumption (quad_index_conversion);
841  return memory;
842  }
843 
844 
845 
846  template <int dim, typename Number>
848  {
849  std::size_t
854  memory += sizeof (this);
855  return memory;
856  }
857 
858 
859 
860  template <int dim, typename Number>
861  template <typename STREAM>
863  (STREAM &out,
864  const SizeInfo &size_info) const
865  {
866  // print_memory_statistics involves global communication, so we can
867  // disable the check here only if no processor has any such data
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,
871  MPI_MAX, size_info.communicator);
872 #else
873  unsigned int general_size_glob = jacobians.size();
874 #endif
875  if (general_size_glob > 0)
876  {
877  out << " Memory Jacobian data: ";
878  size_info.print_memory_statistics
879  (out, MemoryConsumption::memory_consumption (jacobians) +
881  out << " Memory second derivative data: ";
882  size_info.print_memory_statistics
883  (out,MemoryConsumption::memory_consumption (jacobians_grad_diag) +
884  MemoryConsumption::memory_consumption (jacobians_grad_upper));
885  }
886 
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,
890  MPI_MAX, size_info.communicator);
891 #else
892  unsigned int quad_size_glob = quadrature_points.size();
893 #endif
894  if (quad_size_glob > 0)
895  {
896  out << " Memory quadrature points: ";
897  size_info.print_memory_statistics
898  (out, MemoryConsumption::memory_consumption (rowstart_q_points) +
899  MemoryConsumption::memory_consumption (quadrature_points));
900  }
901  }
902 
903 
904 
905  template <int dim, typename Number>
906  template <typename STREAM>
908  const SizeInfo &size_info) const
909  {
910  out << " Cell types: ";
911  size_info.print_memory_statistics
913  out << " Memory transformations compr: ";
914  size_info.print_memory_statistics
917  for (unsigned int j=0; j<mapping_data_gen.size(); ++j)
918  {
919  out << " Data component " << j << std::endl;
920  mapping_data_gen[j].print_memory_consumption(out, size_info);
921  }
922  }
923 
924  } // end of namespace MatrixFreeFunctions
925 } // end of namespace internal
926 
927 
928 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
static const unsigned int invalid_unsigned_int
Definition: types.h:191
void push_back(const Quadrature< dim > &new_quadrature)
Definition: q_collection.h:251
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
AlignedVector< std::pair< Tensor< 1, dim, VectorizedArray< Number > >, VectorizedArray< Number > > > cartesian_data
Definition: mapping_info.h:139
AlignedVector< VectorizedArray< Number > > JxW_values
Definition: mapping_info.h:180
AlignedVector< std::pair< Tensor< 2, dim, VectorizedArray< Number > >, VectorizedArray< Number > > > affine_data
Definition: mapping_info.h:154
const DerivativeForm< 1, dim, spacedim > & jacobian(const unsigned int quadrature_point) const
Volume element.
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
void reserve(const size_type size_alloc)
Transformed quadrature points.
std::vector< MappingInfoDependent > mapping_data_gen
Definition: mapping_info.h:284
void push_back(const T in_data)
CellSimilarity::Similarity get_cell_similarity() const
size_type size() 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
No update.
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
Definition: mapping_info.h:124
#define Assert(cond, exc)
Definition: exceptions.h:299
UpdateFlags
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
Definition: mapping_info.h:202
void reinit(const TriaIterator< DoFCellAccessor< DH, level_dof_access > > cell)
Definition: tensor.h:26
void print_memory_consumption(STREAM &out, const SizeInfo &size_info) const
void resize(const size_type size_in, const T &init=T())
const Point< spacedim > & quadrature_point(const unsigned int i) const
Second derivatives of shape functions.
Gradient of volume element.
AlignedVector< Point< dim, VectorizedArray< Number > > > quadrature_points
Definition: mapping_info.h:216
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
Definition: mapping_info.h:247
const unsigned int n_quadrature_points
Definition: fe_values.h:1411
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
Definition: mapping_info.h:174
AlignedVector< Tensor< 2, dim, VectorizedArray< Number > > > jacobians_grad_diag
Definition: mapping_info.h:191
reference back()