Reference documentation for deal.II version 8.4.1
fe_system.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2016 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/memory_consumption.h>
17 #include <deal.II/base/quadrature.h>
18 #include <deal.II/base/qprojector.h>
19 #include <deal.II/grid/tria.h>
20 #include <deal.II/grid/tria_iterator.h>
21 #include <deal.II/dofs/dof_accessor.h>
22 #include <deal.II/fe/mapping.h>
23 #include <deal.II/fe/fe_system.h>
24 #include <deal.II/fe/fe_values.h>
25 
26 #include <sstream>
27 
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 namespace
32 {
33  bool IsNonZero (unsigned int i)
34  {
35  return i>0;
36  }
37 
38  unsigned int count_nonzeros(const std::vector<unsigned int> &vec)
39  {
40  return std::count_if(vec.begin(), vec.end(), IsNonZero);
41  }
42 
43 }
44 /* ----------------------- FESystem::InternalData ------------------- */
45 
46 
47 template <int dim, int spacedim>
48 FESystem<dim,spacedim>::InternalData::InternalData(const unsigned int n_base_elements)
49  :
50  base_fe_datas(n_base_elements),
51  base_fe_output_objects(n_base_elements)
52 {}
53 
54 
55 
56 template <int dim, int spacedim>
58 {
59  // delete pointers and set them to zero to avoid inadvertent use
60  for (unsigned int i=0; i<base_fe_datas.size(); ++i)
61  if (base_fe_datas[i])
62  {
63  delete base_fe_datas[i];
64  base_fe_datas[i] = 0;
65  }
66 }
67 
68 
69 
70 template <int dim, int spacedim>
73 InternalData::get_fe_data (const unsigned int base_no) const
74 {
75  Assert(base_no<base_fe_datas.size(),
76  ExcIndexRange(base_no,0,base_fe_datas.size()));
77  return *base_fe_datas[base_no];
78 }
79 
80 
81 
82 template <int dim, int spacedim>
83 void
85 InternalData::set_fe_data (const unsigned int base_no,
87 {
88  Assert(base_no<base_fe_datas.size(),
89  ExcIndexRange(base_no,0,base_fe_datas.size()));
90  base_fe_datas[base_no]=ptr;
91 }
92 
93 
94 
95 template <int dim, int spacedim>
98 InternalData::get_fe_output_object (const unsigned int base_no) const
99 {
100  Assert(base_no<base_fe_output_objects.size(),
101  ExcIndexRange(base_no,0,base_fe_output_objects.size()));
102  return base_fe_output_objects[base_no];
103 }
104 
105 
106 
107 /* ---------------------------------- FESystem ------------------- */
108 
109 
110 template <int dim, int spacedim>
112 
113 namespace
114 {
115 
121  template <int dim, int spacedim>
123  multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
124  const std::vector<unsigned int> &multiplicities)
125  {
126  AssertDimension(fes.size(), multiplicities.size());
127 
128  unsigned int multiplied_dofs_per_vertex = 0;
129  unsigned int multiplied_dofs_per_line = 0;
130  unsigned int multiplied_dofs_per_quad = 0;
131  unsigned int multiplied_dofs_per_hex = 0;
132 
133  unsigned int multiplied_n_components = 0;
134 
135  unsigned int degree = 0; // degree is the maximal degree of the components
136 
137  for (unsigned int i=0; i<fes.size(); i++)
138  if (multiplicities[i]>0)
139  {
140  multiplied_dofs_per_vertex += fes[i]->dofs_per_vertex * multiplicities[i];
141  multiplied_dofs_per_line += fes[i]->dofs_per_line * multiplicities[i];
142  multiplied_dofs_per_quad += fes[i]->dofs_per_quad * multiplicities[i];
143  multiplied_dofs_per_hex += fes[i]->dofs_per_hex * multiplicities[i];
144 
145  multiplied_n_components+=fes[i]->n_components() * multiplicities[i];
146 
147  degree = std::max(degree, fes[i]->tensor_degree() );
148  }
149 
150  // assume conformity of the first finite element and then take away
151  // bits as indicated by the base elements. if all multiplicities
152  // happen to be zero, then it doesn't matter what we set it to.
153  typename FiniteElementData<dim>::Conformity total_conformity
155  {
156  unsigned int index = 0;
157  for (index=0; index<fes.size(); ++index)
158  if (multiplicities[index]>0)
159  {
160  total_conformity = fes[index]->conforming_space;
161  break;
162  }
163 
164  for (; index<fes.size(); ++index)
165  if (multiplicities[index]>0)
166  total_conformity =
167  typename FiniteElementData<dim>::Conformity(total_conformity
168  &
169  fes[index]->conforming_space);
170  }
171 
172  std::vector<unsigned int> dpo;
173  dpo.push_back(multiplied_dofs_per_vertex);
174  dpo.push_back(multiplied_dofs_per_line);
175  if (dim>1) dpo.push_back(multiplied_dofs_per_quad);
176  if (dim>2) dpo.push_back(multiplied_dofs_per_hex);
177 
179 
180  for (unsigned int base=0; base < fes.size(); ++base)
181  for (unsigned int m = 0; m < multiplicities[base]; ++m)
182  block_indices.push_back(fes[base]->dofs_per_cell);
183 
184  return FiniteElementData<dim> (dpo,
185  multiplied_n_components,
186  degree,
187  total_conformity,
188  block_indices);
189  }
190 
194  template <int dim, int spacedim>
196  multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
197  const unsigned int N1,
198  const FiniteElement<dim,spacedim> *fe2=NULL,
199  const unsigned int N2=0,
200  const FiniteElement<dim,spacedim> *fe3=NULL,
201  const unsigned int N3=0,
202  const FiniteElement<dim,spacedim> *fe4=NULL,
203  const unsigned int N4=0,
204  const FiniteElement<dim,spacedim> *fe5=NULL,
205  const unsigned int N5=0)
206  {
207  std::vector<const FiniteElement<dim,spacedim>*> fes;
208  fes.push_back(fe1);
209  fes.push_back(fe2);
210  fes.push_back(fe3);
211  fes.push_back(fe4);
212  fes.push_back(fe5);
213 
214  std::vector<unsigned int> mult;
215  mult.push_back(N1);
216  mult.push_back(N2);
217  mult.push_back(N3);
218  mult.push_back(N4);
219  mult.push_back(N5);
220  return multiply_dof_numbers(fes, mult);
221  }
222 
223 
224 
230  template <int dim, int spacedim>
231  std::vector<bool>
232  compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
233  const std::vector<unsigned int> &multiplicities)
234  {
235  AssertDimension(fes.size(), multiplicities.size());
236 
237  // first count the number of dofs and components that will emerge from the
238  // given FEs
239  unsigned int n_shape_functions = 0;
240  for (unsigned int i=0; i<fes.size(); ++i)
241  if (multiplicities[i]>0) // check needed as fe might be NULL
242  n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i];
243 
244  // generate the array that will hold the output
245  std::vector<bool> retval (n_shape_functions, false);
246 
247  // finally go through all the shape functions of the base elements, and copy
248  // their flags. this somehow copies the code in build_cell_table, which is
249  // not nice as it uses too much implicit knowledge about the layout of the
250  // individual bases in the composed FE, but there seems no way around...
251  //
252  // for each shape function, copy the flags from the base element to this
253  // one, taking into account multiplicities, and other complications
254  unsigned int total_index = 0;
255  for (unsigned int vertex_number=0;
256  vertex_number<GeometryInfo<dim>::vertices_per_cell;
257  ++vertex_number)
258  {
259  for (unsigned int base=0; base<fes.size(); ++base)
260  for (unsigned int m=0; m<multiplicities[base]; ++m)
261  for (unsigned int local_index = 0;
262  local_index < fes[base]->dofs_per_vertex;
263  ++local_index, ++total_index)
264  {
265  const unsigned int index_in_base
266  = (fes[base]->dofs_per_vertex*vertex_number +
267  local_index);
268 
269  Assert (index_in_base < fes[base]->dofs_per_cell,
270  ExcInternalError());
271  retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
272  }
273  }
274 
275  // 2. Lines
277  for (unsigned int line_number= 0;
278  line_number != GeometryInfo<dim>::lines_per_cell;
279  ++line_number)
280  {
281  for (unsigned int base=0; base<fes.size(); ++base)
282  for (unsigned int m=0; m<multiplicities[base]; ++m)
283  for (unsigned int local_index = 0;
284  local_index < fes[base]->dofs_per_line;
285  ++local_index, ++total_index)
286  {
287  const unsigned int index_in_base
288  = (fes[base]->dofs_per_line*line_number +
289  local_index +
290  fes[base]->first_line_index);
291 
292  Assert (index_in_base < fes[base]->dofs_per_cell,
293  ExcInternalError());
294  retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
295  }
296  }
297 
298  // 3. Quads
300  for (unsigned int quad_number= 0;
301  quad_number != GeometryInfo<dim>::quads_per_cell;
302  ++quad_number)
303  {
304  for (unsigned int base=0; base<fes.size(); ++base)
305  for (unsigned int m=0; m<multiplicities[base]; ++m)
306  for (unsigned int local_index = 0;
307  local_index < fes[base]->dofs_per_quad;
308  ++local_index, ++total_index)
309  {
310  const unsigned int index_in_base
311  = (fes[base]->dofs_per_quad*quad_number +
312  local_index +
313  fes[base]->first_quad_index);
314 
315  Assert (index_in_base < fes[base]->dofs_per_cell,
316  ExcInternalError());
317  retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
318  }
319  }
320 
321  // 4. Hexes
323  for (unsigned int hex_number= 0;
324  hex_number != GeometryInfo<dim>::hexes_per_cell;
325  ++hex_number)
326  {
327  for (unsigned int base=0; base<fes.size(); ++base)
328  for (unsigned int m=0; m<multiplicities[base]; ++m)
329  for (unsigned int local_index = 0;
330  local_index < fes[base]->dofs_per_hex;
331  ++local_index, ++total_index)
332  {
333  const unsigned int index_in_base
334  = (fes[base]->dofs_per_hex*hex_number +
335  local_index +
336  fes[base]->first_hex_index);
337 
338  Assert (index_in_base < fes[base]->dofs_per_cell,
339  ExcInternalError());
340  retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
341  }
342  }
343 
344  Assert (total_index == n_shape_functions, ExcInternalError());
345 
346  return retval;
347  }
348 
349 
350 
357  template <int dim, int spacedim>
358  std::vector<bool>
359  compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> *fe1,
360  const unsigned int N1,
361  const FiniteElement<dim,spacedim> *fe2=NULL,
362  const unsigned int N2=0,
363  const FiniteElement<dim,spacedim> *fe3=NULL,
364  const unsigned int N3=0,
365  const FiniteElement<dim,spacedim> *fe4=NULL,
366  const unsigned int N4=0,
367  const FiniteElement<dim,spacedim> *fe5=NULL,
368  const unsigned int N5=0)
369  {
370  std::vector<const FiniteElement<dim,spacedim>*> fe_list;
371  std::vector<unsigned int> multiplicities;
372 
373  fe_list.push_back (fe1);
374  multiplicities.push_back (N1);
375 
376  fe_list.push_back (fe2);
377  multiplicities.push_back (N2);
378 
379  fe_list.push_back (fe3);
380  multiplicities.push_back (N3);
381 
382  fe_list.push_back (fe4);
383  multiplicities.push_back (N4);
384 
385  fe_list.push_back (fe5);
386  multiplicities.push_back (N5);
387  return compute_restriction_is_additive_flags (fe_list, multiplicities);
388  }
389 
390 
391 
396  template <int dim, int spacedim>
397  std::vector<ComponentMask>
398  compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
399  const std::vector<unsigned int> &multiplicities)
400  {
401  AssertDimension(fes.size(), multiplicities.size());
402 
403  // first count the number of dofs and components that will emerge from the
404  // given FEs
405  unsigned int n_shape_functions = 0;
406  for (unsigned int i=0; i<fes.size(); ++i)
407  if (multiplicities[i]>0) //needed because fe might be NULL
408  n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i];
409 
410  unsigned int n_components = 0;
411  for (unsigned int i=0; i<fes.size(); ++i)
412  if (multiplicities[i]>0) //needed because fe might be NULL
413  n_components += fes[i]->n_components() * multiplicities[i];
414 
415  // generate the array that will hold the output
416  std::vector<std::vector<bool> >
417  retval (n_shape_functions, std::vector<bool> (n_components, false));
418 
419  // finally go through all the shape functions of the base elements, and copy
420  // their flags. this somehow copies the code in build_cell_table, which is
421  // not nice as it uses too much implicit knowledge about the layout of the
422  // individual bases in the composed FE, but there seems no way around...
423  //
424  // for each shape function, copy the non-zero flags from the base element to
425  // this one, taking into account multiplicities, multiple components in base
426  // elements, and other complications
427  unsigned int total_index = 0;
428  for (unsigned int vertex_number=0;
429  vertex_number<GeometryInfo<dim>::vertices_per_cell;
430  ++vertex_number)
431  {
432  unsigned int comp_start = 0;
433  for (unsigned int base=0; base<fes.size(); ++base)
434  for (unsigned int m=0; m<multiplicities[base];
435  ++m, comp_start+=fes[base]->n_components())
436  for (unsigned int local_index = 0;
437  local_index < fes[base]->dofs_per_vertex;
438  ++local_index, ++total_index)
439  {
440  const unsigned int index_in_base
441  = (fes[base]->dofs_per_vertex*vertex_number +
442  local_index);
443 
444  Assert (comp_start+fes[base]->n_components() <=
445  retval[total_index].size(),
446  ExcInternalError());
447  for (unsigned int c=0; c<fes[base]->n_components(); ++c)
448  {
449  Assert (c < fes[base]->get_nonzero_components(index_in_base).size(),
450  ExcInternalError());
451  retval[total_index][comp_start+c]
452  = fes[base]->get_nonzero_components(index_in_base)[c];
453  }
454  }
455  }
456 
457  // 2. Lines
459  for (unsigned int line_number= 0;
460  line_number != GeometryInfo<dim>::lines_per_cell;
461  ++line_number)
462  {
463  unsigned int comp_start = 0;
464  for (unsigned int base=0; base<fes.size(); ++base)
465  for (unsigned int m=0; m<multiplicities[base];
466  ++m, comp_start+=fes[base]->n_components())
467  for (unsigned int local_index = 0;
468  local_index < fes[base]->dofs_per_line;
469  ++local_index, ++total_index)
470  {
471  const unsigned int index_in_base
472  = (fes[base]->dofs_per_line*line_number +
473  local_index +
474  fes[base]->first_line_index);
475 
476  Assert (comp_start+fes[base]->n_components() <=
477  retval[total_index].size(),
478  ExcInternalError());
479  for (unsigned int c=0; c<fes[base]->n_components(); ++c)
480  {
481  Assert (c < fes[base]->get_nonzero_components(index_in_base).size(),
482  ExcInternalError());
483  retval[total_index][comp_start+c]
484  = fes[base]->get_nonzero_components(index_in_base)[c];
485  }
486  }
487  }
488 
489  // 3. Quads
491  for (unsigned int quad_number= 0;
492  quad_number != GeometryInfo<dim>::quads_per_cell;
493  ++quad_number)
494  {
495  unsigned int comp_start = 0;
496  for (unsigned int base=0; base<fes.size(); ++base)
497  for (unsigned int m=0; m<multiplicities[base];
498  ++m, comp_start+=fes[base]->n_components())
499  for (unsigned int local_index = 0;
500  local_index < fes[base]->dofs_per_quad;
501  ++local_index, ++total_index)
502  {
503  const unsigned int index_in_base
504  = (fes[base]->dofs_per_quad*quad_number +
505  local_index +
506  fes[base]->first_quad_index);
507 
508  Assert (comp_start+fes[base]->n_components() <=
509  retval[total_index].size(),
510  ExcInternalError());
511  for (unsigned int c=0; c<fes[base]->n_components(); ++c)
512  {
513  Assert (c < fes[base]->get_nonzero_components(index_in_base).size(),
514  ExcInternalError());
515  retval[total_index][comp_start+c]
516  = fes[base]->get_nonzero_components(index_in_base)[c];
517  }
518  }
519  }
520 
521  // 4. Hexes
523  for (unsigned int hex_number= 0;
524  hex_number != GeometryInfo<dim>::hexes_per_cell;
525  ++hex_number)
526  {
527  unsigned int comp_start = 0;
528  for (unsigned int base=0; base<fes.size(); ++base)
529  for (unsigned int m=0; m<multiplicities[base];
530  ++m, comp_start+=fes[base]->n_components())
531  for (unsigned int local_index = 0;
532  local_index < fes[base]->dofs_per_hex;
533  ++local_index, ++total_index)
534  {
535  const unsigned int index_in_base
536  = (fes[base]->dofs_per_hex*hex_number +
537  local_index +
538  fes[base]->first_hex_index);
539 
540  Assert (comp_start+fes[base]->n_components() <=
541  retval[total_index].size(),
542  ExcInternalError());
543  for (unsigned int c=0; c<fes[base]->n_components(); ++c)
544  {
545  Assert (c < fes[base]->get_nonzero_components(index_in_base).size(),
546  ExcInternalError());
547  retval[total_index][comp_start+c]
548  = fes[base]->get_nonzero_components(index_in_base)[c];
549  }
550  }
551  }
552 
553  Assert (total_index == n_shape_functions, ExcInternalError());
554 
555  // now copy the vector<vector<bool> > into a vector<ComponentMask>.
556  // this appears complicated but we do it this way since it's just
557  // awkward to generate ComponentMasks directly and so we need the
558  // recourse of the inner vector<bool> anyway.
559  std::vector<ComponentMask> xretval (retval.size());
560  for (unsigned int i=0; i<retval.size(); ++i)
561  xretval[i] = ComponentMask(retval[i]);
562  return xretval;
563  }
564 
565 
569  template <int dim, int spacedim>
570  std::vector<ComponentMask>
571  compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
572  const unsigned int N1,
573  const FiniteElement<dim,spacedim> *fe2=NULL,
574  const unsigned int N2=0,
575  const FiniteElement<dim,spacedim> *fe3=NULL,
576  const unsigned int N3=0,
577  const FiniteElement<dim,spacedim> *fe4=NULL,
578  const unsigned int N4=0,
579  const FiniteElement<dim,spacedim> *fe5=NULL,
580  const unsigned int N5=0)
581  {
582  std::vector<const FiniteElement<dim,spacedim>*> fe_list;
583  std::vector<unsigned int> multiplicities;
584 
585  fe_list.push_back (fe1);
586  multiplicities.push_back (N1);
587 
588  fe_list.push_back (fe2);
589  multiplicities.push_back (N2);
590 
591  fe_list.push_back (fe3);
592  multiplicities.push_back (N3);
593 
594  fe_list.push_back (fe4);
595  multiplicities.push_back (N4);
596 
597  fe_list.push_back (fe5);
598  multiplicities.push_back (N5);
599 
600  return compute_nonzero_components (fe_list, multiplicities);
601  }
602 }
603 
604 
605 
606 template <int dim, int spacedim>
608  const unsigned int n_elements) :
609  FiniteElement<dim,spacedim> (multiply_dof_numbers(&fe, n_elements),
610  compute_restriction_is_additive_flags (&fe, n_elements),
611  compute_nonzero_components(&fe, n_elements)),
612  base_elements((n_elements>0))
613 {
614  std::vector<const FiniteElement<dim,spacedim>*> fes;
615  fes.push_back(&fe);
616  std::vector<unsigned int> multiplicities;
617  multiplicities.push_back(n_elements);
618  initialize(fes, multiplicities);
619 }
620 
621 
622 
623 template <int dim, int spacedim>
625  const unsigned int n1,
626  const FiniteElement<dim,spacedim> &fe2,
627  const unsigned int n2) :
628  FiniteElement<dim,spacedim> (multiply_dof_numbers(&fe1, n1, &fe2, n2),
629  compute_restriction_is_additive_flags (&fe1, n1,
630  &fe2, n2),
631  compute_nonzero_components(&fe1, n1,
632  &fe2, n2)),
633  base_elements((n1>0)+(n2>0))
634 {
635  std::vector<const FiniteElement<dim,spacedim>*> fes;
636  fes.push_back(&fe1);
637  fes.push_back(&fe2);
638  std::vector<unsigned int> multiplicities;
639  multiplicities.push_back(n1);
640  multiplicities.push_back(n2);
641  initialize(fes, multiplicities);
642 }
643 
644 
645 
646 template <int dim, int spacedim>
648  const unsigned int n1,
649  const FiniteElement<dim,spacedim> &fe2,
650  const unsigned int n2,
651  const FiniteElement<dim,spacedim> &fe3,
652  const unsigned int n3) :
653  FiniteElement<dim,spacedim> (multiply_dof_numbers(&fe1, n1,
654  &fe2, n2,
655  &fe3, n3),
656  compute_restriction_is_additive_flags (&fe1, n1,
657  &fe2, n2,
658  &fe3, n3),
659  compute_nonzero_components(&fe1, n1,
660  &fe2, n2,
661  &fe3, n3)),
662  base_elements((n1>0)+(n2>0)+(n3>0))
663 {
664  std::vector<const FiniteElement<dim,spacedim>*> fes;
665  fes.push_back(&fe1);
666  fes.push_back(&fe2);
667  fes.push_back(&fe3);
668  std::vector<unsigned int> multiplicities;
669  multiplicities.push_back(n1);
670  multiplicities.push_back(n2);
671  multiplicities.push_back(n3);
672  initialize(fes, multiplicities);
673 }
674 
675 
676 
677 template <int dim, int spacedim>
679  const unsigned int n1,
680  const FiniteElement<dim,spacedim> &fe2,
681  const unsigned int n2,
682  const FiniteElement<dim,spacedim> &fe3,
683  const unsigned int n3,
684  const FiniteElement<dim,spacedim> &fe4,
685  const unsigned int n4) :
686  FiniteElement<dim,spacedim> (multiply_dof_numbers(&fe1, n1,
687  &fe2, n2,
688  &fe3, n3,
689  &fe4, n4),
690  compute_restriction_is_additive_flags (&fe1, n1,
691  &fe2, n2,
692  &fe3, n3,
693  &fe4, n4),
694  compute_nonzero_components(&fe1, n1,
695  &fe2, n2,
696  &fe3, n3,
697  &fe4 ,n4)),
698  base_elements((n1>0)+(n2>0)+(n3>0)+(n4>0))
699 {
700  std::vector<const FiniteElement<dim,spacedim>*> fes;
701  fes.push_back(&fe1);
702  fes.push_back(&fe2);
703  fes.push_back(&fe3);
704  fes.push_back(&fe4);
705  std::vector<unsigned int> multiplicities;
706  multiplicities.push_back(n1);
707  multiplicities.push_back(n2);
708  multiplicities.push_back(n3);
709  multiplicities.push_back(n4);
710  initialize(fes, multiplicities);
711 }
712 
713 
714 
715 template <int dim, int spacedim>
717  const unsigned int n1,
718  const FiniteElement<dim,spacedim> &fe2,
719  const unsigned int n2,
720  const FiniteElement<dim,spacedim> &fe3,
721  const unsigned int n3,
722  const FiniteElement<dim,spacedim> &fe4,
723  const unsigned int n4,
724  const FiniteElement<dim,spacedim> &fe5,
725  const unsigned int n5) :
726  FiniteElement<dim,spacedim> (multiply_dof_numbers(&fe1, n1,
727  &fe2, n2,
728  &fe3, n3,
729  &fe4, n4,
730  &fe5, n5),
731  compute_restriction_is_additive_flags (&fe1, n1,
732  &fe2, n2,
733  &fe3, n3,
734  &fe4, n4,
735  &fe5, n5),
736  compute_nonzero_components(&fe1, n1,
737  &fe2, n2,
738  &fe3, n3,
739  &fe4 ,n4,
740  &fe5, n5)),
741  base_elements((n1>0)+(n2>0)+(n3>0)+(n4>0)+(n5>0))
742 {
743  std::vector<const FiniteElement<dim,spacedim>*> fes;
744  fes.push_back(&fe1);
745  fes.push_back(&fe2);
746  fes.push_back(&fe3);
747  fes.push_back(&fe4);
748  fes.push_back(&fe5);
749  std::vector<unsigned int> multiplicities;
750  multiplicities.push_back(n1);
751  multiplicities.push_back(n2);
752  multiplicities.push_back(n3);
753  multiplicities.push_back(n4);
754  multiplicities.push_back(n5);
755  initialize(fes, multiplicities);
756 }
757 
758 
759 
760 template <int dim, int spacedim>
762  const std::vector<const FiniteElement<dim,spacedim>*> &fes,
763  const std::vector<unsigned int> &multiplicities)
764  :
765  FiniteElement<dim,spacedim> (multiply_dof_numbers(fes, multiplicities),
766  compute_restriction_is_additive_flags (fes, multiplicities),
767  compute_nonzero_components(fes, multiplicities)),
768  base_elements(count_nonzeros(multiplicities))
769 {
770  initialize(fes, multiplicities);
771 }
772 
773 
774 
775 template <int dim, int spacedim>
777 {}
778 
779 
780 
781 template <int dim, int spacedim>
782 std::string
784 {
785  // note that the
786  // FETools::get_fe_from_name
787  // function depends on the
788  // particular format of the string
789  // this function returns, so they
790  // have to be kept in synch
791 
792  std::ostringstream namebuf;
793 
794  namebuf << "FESystem<"
795  << Utilities::dim_string(dim,spacedim)
796  << ">[";
797  for (unsigned int i=0; i< this->n_base_elements(); ++i)
798  {
799  namebuf << base_element(i).get_name();
800  if (this->element_multiplicity(i) != 1)
801  namebuf << '^' << this->element_multiplicity(i);
802  if (i != this->n_base_elements()-1)
803  namebuf << '-';
804  }
805  namebuf << ']';
806 
807  return namebuf.str();
808 }
809 
810 
811 
812 template <int dim, int spacedim>
815 {
816  std::vector< const FiniteElement<dim,spacedim>* > fes;
817  std::vector<unsigned int> multiplicities;
818 
819  for (unsigned int i=0; i<this->n_base_elements(); i++)
820  {
821  fes.push_back( & base_element(i) );
822  multiplicities.push_back(this->element_multiplicity(i) );
823  }
824  return new FESystem<dim,spacedim>(fes, multiplicities);
825 }
826 
827 
828 
829 template <int dim, int spacedim>
830 double
832  const Point<dim> &p) const
833 {
834  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
835  Assert (this->is_primitive(i),
837 
838  return (base_element(this->system_to_base_table[i].first.first)
839  .shape_value(this->system_to_base_table[i].second, p));
840 }
841 
842 
843 
844 template <int dim, int spacedim>
845 double
847  const Point<dim> &p,
848  const unsigned int component) const
849 {
850  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
851  Assert (component < this->n_components(),
852  ExcIndexRange (component, 0, this->n_components()));
853 
854  // if this value is supposed to be
855  // zero, then return right away...
856  if (this->nonzero_components[i][component] == false)
857  return 0;
858 
859  // ...otherwise: first find out to
860  // which of the base elements this
861  // desired component belongs, and
862  // which component within this base
863  // element it is
864  const unsigned int base = this->component_to_base_index(component).first;
865  const unsigned int component_in_base = this->component_to_base_index(component).second;
866 
867  // then get value from base
868  // element. note that that will
869  // throw an error should the
870  // respective shape function not be
871  // primitive; thus, there is no
872  // need to check this here
873  return (base_element(base).
874  shape_value_component(this->system_to_base_table[i].second,
875  p,
876  component_in_base));
877 }
878 
879 
880 
881 template <int dim, int spacedim>
884  const Point<dim> &p) const
885 {
886  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
887  Assert (this->is_primitive(i),
889 
890  return (base_element(this->system_to_base_table[i].first.first)
891  .shape_grad(this->system_to_base_table[i].second, p));
892 }
893 
894 
895 
896 template <int dim, int spacedim>
899  const Point<dim> &p,
900  const unsigned int component) const
901 {
902  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
903  Assert (component < this->n_components(),
904  ExcIndexRange (component, 0, this->n_components()));
905 
906  // if this value is supposed to be zero, then return right away...
907  if (this->nonzero_components[i][component] == false)
908  return Tensor<1,dim>();
909 
910  // ...otherwise: first find out to which of the base elements this desired
911  // component belongs, and which component within this base element it is
912  const unsigned int base = this->component_to_base_index(component).first;
913  const unsigned int component_in_base = this->component_to_base_index(component).second;
914 
915  // then get value from base element. note that that will throw an error
916  // should the respective shape function not be primitive; thus, there is no
917  // need to check this here
918  return (base_element(base).
919  shape_grad_component(this->system_to_base_table[i].second,
920  p,
921  component_in_base));
922 }
923 
924 
925 
926 template <int dim, int spacedim>
929  const Point<dim> &p) const
930 {
931  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
932  Assert (this->is_primitive(i),
934 
935  return (base_element(this->system_to_base_table[i].first.first)
936  .shape_grad_grad(this->system_to_base_table[i].second, p));
937 }
938 
939 
940 
941 template <int dim, int spacedim>
944  const Point<dim> &p,
945  const unsigned int component) const
946 {
947  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
948  Assert (component < this->n_components(),
949  ExcIndexRange (component, 0, this->n_components()));
950 
951  // if this value is supposed to be zero, then return right away...
952  if (this->nonzero_components[i][component] == false)
953  return Tensor<2,dim>();
954 
955  // ...otherwise: first find out to which of the base elements this desired
956  // component belongs, and which component within this base element it is
957  const unsigned int base = this->component_to_base_index(component).first;
958  const unsigned int component_in_base = this->component_to_base_index(component).second;
959 
960  // then get value from base element. note that that will throw an error
961  // should the respective shape function not be primitive; thus, there is no
962  // need to check this here
963  return (base_element(base).
964  shape_grad_grad_component(this->system_to_base_table[i].second,
965  p,
966  component_in_base));
967 }
968 
969 
970 
971 template <int dim, int spacedim>
974  const Point<dim> &p) const
975 {
976  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
977  Assert (this->is_primitive(i),
979 
980  return (base_element(this->system_to_base_table[i].first.first)
981  .shape_3rd_derivative(this->system_to_base_table[i].second, p));
982 }
983 
984 
985 
986 template <int dim, int spacedim>
989  const Point<dim> &p,
990  const unsigned int component) const
991 {
992  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
993  Assert (component < this->n_components(),
994  ExcIndexRange (component, 0, this->n_components()));
995 
996  // if this value is supposed to be zero, then return right away...
997  if (this->nonzero_components[i][component] == false)
998  return Tensor<3,dim>();
999 
1000  // ...otherwise: first find out to which of the base elements this desired
1001  // component belongs, and which component within this base element it is
1002  const unsigned int base = this->component_to_base_index(component).first;
1003  const unsigned int component_in_base = this->component_to_base_index(component).second;
1004 
1005  // then get value from base element. note that that will throw an error
1006  // should the respective shape function not be primitive; thus, there is no
1007  // need to check this here
1008  return (base_element(base).
1009  shape_3rd_derivative_component(this->system_to_base_table[i].second,
1010  p,
1011  component_in_base));
1012 }
1013 
1014 
1015 
1016 template <int dim, int spacedim>
1019  const Point<dim> &p) const
1020 {
1021  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
1022  Assert (this->is_primitive(i),
1024 
1025  return (base_element(this->system_to_base_table[i].first.first)
1026  .shape_4th_derivative(this->system_to_base_table[i].second, p));
1027 }
1028 
1029 
1030 
1031 template <int dim, int spacedim>
1034  const Point<dim> &p,
1035  const unsigned int component) const
1036 {
1037  Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
1038  Assert (component < this->n_components(),
1039  ExcIndexRange (component, 0, this->n_components()));
1040 
1041  // if this value is supposed to be zero, then return right away...
1042  if (this->nonzero_components[i][component] == false)
1043  return Tensor<4,dim>();
1044 
1045  // ...otherwise: first find out to which of the base elements this desired
1046  // component belongs, and which component within this base element it is
1047  const unsigned int base = this->component_to_base_index(component).first;
1048  const unsigned int component_in_base = this->component_to_base_index(component).second;
1049 
1050  // then get value from base element. note that that will throw an error
1051  // should the respective shape function not be primitive; thus, there is no
1052  // need to check this here
1053  return (base_element(base).
1054  shape_4th_derivative_component(this->system_to_base_table[i].second,
1055  p,
1056  component_in_base));
1057 }
1058 
1059 
1060 
1061 template <int dim, int spacedim>
1062 void
1064  const FiniteElement<dim,spacedim> &x_source_fe,
1065  FullMatrix<double> &interpolation_matrix) const
1066 {
1067  // check that the size of the matrices is correct. for historical
1068  // reasons, if you call matrix.reinit(8,0), it sets the sizes
1069  // to m==n==0 internally. this may happen when we use a FE_Nothing,
1070  // so write the test in a more lenient way
1071  Assert ((interpolation_matrix.m() == this->dofs_per_cell)
1072  ||
1073  (x_source_fe.dofs_per_cell == 0),
1074  ExcDimensionMismatch (interpolation_matrix.m(),
1075  this->dofs_per_cell));
1076  Assert ((interpolation_matrix.n() == x_source_fe.dofs_per_cell)
1077  ||
1078  (this->dofs_per_cell == 0),
1079  ExcDimensionMismatch (interpolation_matrix.m(),
1080  x_source_fe.dofs_per_cell));
1081 
1082  // there are certain conditions that the two elements have to satisfy so
1083  // that this can work.
1084  //
1085  // condition 1: the other element must also be a system element
1086 
1087  typedef FiniteElement<dim,spacedim> FEL;
1088  AssertThrow ((x_source_fe.get_name().find ("FESystem<") == 0)
1089  ||
1090  (dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe) != 0),
1091  typename FEL::
1092  ExcInterpolationNotImplemented());
1093 
1094  // ok, source is a system element, so we may be able to do the work
1095  const FESystem<dim,spacedim> &source_fe
1096  = dynamic_cast<const FESystem<dim,spacedim>&>(x_source_fe);
1097 
1098  // condition 2: same number of basis elements
1099  AssertThrow (this->n_base_elements() == source_fe.n_base_elements(),
1100  typename FEL::
1101  ExcInterpolationNotImplemented());
1102 
1103  // condition 3: same number of basis elements
1104  for (unsigned int i=0; i< this->n_base_elements(); ++i)
1105  AssertThrow (this->element_multiplicity(i) ==
1106  source_fe.element_multiplicity(i),
1107  typename FEL::
1108  ExcInterpolationNotImplemented());
1109 
1110  // ok, so let's try whether it works:
1111 
1112  // first let's see whether all the basis elements actually generate their
1113  // interpolation matrices. if we get past the following loop, then
1114  // apparently none of the called base elements threw an exception, so we're
1115  // fine continuing and assembling the one big matrix from the small ones of
1116  // the base elements
1117  std::vector<FullMatrix<double> > base_matrices (this->n_base_elements());
1118  for (unsigned int i=0; i<this->n_base_elements(); ++i)
1119  {
1120  base_matrices[i].reinit (base_element(i).dofs_per_cell,
1121  source_fe.base_element(i).dofs_per_cell);
1122  base_element(i).get_interpolation_matrix (source_fe.base_element(i),
1123  base_matrices[i]);
1124  }
1125 
1126  // first clear big matrix, to make sure that entries that would couple
1127  // different bases (or multiplicity indices) are really zero. then assign
1128  // entries
1129  interpolation_matrix = 0;
1130  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
1131  for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
1132  if (this->system_to_base_table[i].first ==
1133  source_fe.system_to_base_table[j].first)
1134  interpolation_matrix(i,j)
1135  = (base_matrices[this->system_to_base_table[i].first.first]
1136  (this->system_to_base_table[i].second,
1137  source_fe.system_to_base_table[j].second));
1138 }
1139 
1140 
1141 
1142 template <int dim, int spacedim>
1143 const FullMatrix<double> &
1145 ::get_restriction_matrix (const unsigned int child,
1146  const RefinementCase<dim> &refinement_case) const
1147 {
1149  ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
1150  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
1151  ExcMessage("Restriction matrices are only available for refined cells!"));
1152  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
1153  ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
1154 
1155  // initialization upon first request
1156  if (this->restriction[refinement_case-1][child].n() == 0)
1157  {
1158  Threads::Mutex::ScopedLock lock(this->mutex);
1159 
1160  // check if updated while waiting for lock
1161  if (this->restriction[refinement_case-1][child].n() ==
1162  this->dofs_per_cell)
1163  return this->restriction[refinement_case-1][child];
1164 
1165  // Check if some of the matrices of the base elements are void.
1166  bool do_restriction = true;
1167 
1168  // shortcut for accessing local restrictions further down
1169  std::vector<const FullMatrix<double> *>
1170  base_matrices(this->n_base_elements());
1171 
1172  for (unsigned int i=0; i<this->n_base_elements(); ++i)
1173  {
1174  base_matrices[i] =
1175  &base_element(i).get_restriction_matrix(child, refinement_case);
1176  if (base_matrices[i]->n() != base_element(i).dofs_per_cell)
1177  do_restriction = false;
1178  }
1179  Assert(do_restriction,
1181 
1182  // if we did not encounter void matrices, initialize the matrix sizes
1183  if (do_restriction)
1184  {
1185  FullMatrix<double> restriction(this->dofs_per_cell,
1186  this->dofs_per_cell);
1187 
1188  // distribute the matrices of the base finite elements to the
1189  // matrices of this object. for this, loop over all degrees of
1190  // freedom and take the respective entry of the underlying base
1191  // element.
1192  //
1193  // note that we by definition of a base element, they are
1194  // independent, i.e. do not couple. only DoFs that belong to the
1195  // same instance of a base element may couple
1196  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
1197  for (unsigned int j=0; j<this->dofs_per_cell; ++j)
1198  {
1199  // first find out to which base element indices i and j
1200  // belong, and which instance thereof in case the base element
1201  // has a multiplicity greater than one. if they should not
1202  // happen to belong to the same instance of a base element,
1203  // then they cannot couple, so go on with the next index
1204  if (this->system_to_base_table[i].first !=
1205  this->system_to_base_table[j].first)
1206  continue;
1207 
1208  // so get the common base element and the indices therein:
1209  const unsigned int
1210  base = this->system_to_base_table[i].first.first;
1211 
1212  const unsigned int
1213  base_index_i = this->system_to_base_table[i].second,
1214  base_index_j = this->system_to_base_table[j].second;
1215 
1216  // if we are sure that DoFs i and j may couple, then copy
1217  // entries of the matrices:
1218  restriction(i,j) = (*base_matrices[base])(base_index_i,base_index_j);
1219  }
1220 
1221  restriction.swap(const_cast<FullMatrix<double> &>
1222  (this->restriction[refinement_case-1][child]));
1223  }
1224  }
1225 
1226  return this->restriction[refinement_case-1][child];
1227 }
1228 
1229 
1230 
1231 template <int dim, int spacedim>
1232 const FullMatrix<double> &
1234 ::get_prolongation_matrix (const unsigned int child,
1235  const RefinementCase<dim> &refinement_case) const
1236 {
1238  ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
1239  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
1240  ExcMessage("Restriction matrices are only available for refined cells!"));
1241  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
1242  ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
1243 
1244  // initialization upon first request, construction completely analogous to
1245  // restriction matrix
1246  if (this->prolongation[refinement_case-1][child].n() == 0)
1247  {
1248  Threads::Mutex::ScopedLock lock(this->mutex);
1249 
1250  if (this->prolongation[refinement_case-1][child].n() ==
1251  this->dofs_per_cell)
1252  return this->prolongation[refinement_case-1][child];
1253 
1254  bool do_prolongation = true;
1255  std::vector<const FullMatrix<double> *>
1256  base_matrices(this->n_base_elements());
1257  for (unsigned int i=0; i<this->n_base_elements(); ++i)
1258  {
1259  base_matrices[i] =
1260  &base_element(i).get_prolongation_matrix(child, refinement_case);
1261  if (base_matrices[i]->n() != base_element(i).dofs_per_cell)
1262  do_prolongation = false;
1263  }
1264  Assert(do_prolongation,
1266 
1267  if (do_prolongation)
1268  {
1269  FullMatrix<double> prolongate (this->dofs_per_cell,
1270  this->dofs_per_cell);
1271 
1272  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
1273  for (unsigned int j=0; j<this->dofs_per_cell; ++j)
1274  {
1275  if (this->system_to_base_table[i].first !=
1276  this->system_to_base_table[j].first)
1277  continue;
1278  const unsigned int
1279  base = this->system_to_base_table[i].first.first;
1280 
1281  const unsigned int
1282  base_index_i = this->system_to_base_table[i].second,
1283  base_index_j = this->system_to_base_table[j].second;
1284  prolongate(i,j) = (*base_matrices[base])(base_index_i,base_index_j);
1285  }
1286  prolongate.swap(const_cast<FullMatrix<double> &>
1287  (this->prolongation[refinement_case-1][child]));
1288  }
1289  }
1290 
1291  return this->prolongation[refinement_case-1][child];
1292 }
1293 
1294 
1295 template <int dim, int spacedim>
1296 unsigned int
1298 face_to_cell_index (const unsigned int face_dof_index,
1299  const unsigned int face,
1300  const bool face_orientation,
1301  const bool face_flip,
1302  const bool face_rotation) const
1303 {
1304  // we need to ask the base elements how they want to translate
1305  // the DoFs within their own numbering. thus, translate to
1306  // the base element numbering and then back
1307  const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
1308  face_base_index = this->face_system_to_base_index(face_dof_index);
1309 
1310  const unsigned int
1311  base_face_to_cell_index
1312  = this->base_element(face_base_index.first.first).face_to_cell_index (face_base_index.second,
1313  face,
1314  face_orientation,
1315  face_flip,
1316  face_rotation);
1317 
1318  // it would be nice if we had a base_to_system_index function, but
1319  // all that exists is a component_to_system_index function. we can't do
1320  // this here because it won't work for non-primitive elements. consequently,
1321  // simply do a loop over all dofs till we find whether it corresponds
1322  // to the one we're interested in -- crude, maybe, but works for now
1323  const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
1324  target = std::make_pair (face_base_index.first, base_face_to_cell_index);
1325  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
1326  if (this->system_to_base_index(i) == target)
1327  return i;
1328 
1329  Assert (false, ExcInternalError());
1331 }
1332 
1333 
1334 
1335 //---------------------------------------------------------------------------
1336 // Data field initialization
1337 //---------------------------------------------------------------------------
1338 
1339 
1340 
1341 template <int dim, int spacedim>
1344 {
1346  // generate maximal set of flags
1347  // that are necessary
1348  for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
1349  out |= base_element(base_no).requires_update_flags (flags);
1350  return out;
1351 }
1352 
1353 
1354 
1355 template <int dim, int spacedim>
1358 get_data (const UpdateFlags flags,
1359  const Mapping<dim,spacedim> &mapping,
1360  const Quadrature<dim> &quadrature,
1362 {
1363  // create an internal data object and set the update flags we will need
1364  // to deal with. the current object does not make use of these flags,
1365  // but we need to nevertheless set them correctly since we look
1366  // into the update_each flag of base elements in fill_fe_values,
1367  // and so the current object's update_each flag needs to be
1368  // correct in case the current FESystem is a base element for another,
1369  // higher-level FESystem itself.
1370  InternalData *data = new InternalData(this->n_base_elements());
1371  data->update_each = requires_update_flags(flags);
1372 
1373  // get data objects from each of the base elements and store
1374  // them. one might think that doing this in parallel (over the
1375  // base elements) would be a good idea, but this turns out to
1376  // be wrong because we would then run these jobs on different
1377  // threads/processors and this allocates memory in different
1378  // NUMA domains; this has large detrimental effects when later
1379  // writing into these objects in fill_fe_*_values. all of this
1380  // is particularly true when using FEValues objects in
1381  // WorkStream contexts where we explicitly make sure that
1382  // every function only uses objects previously allocated
1383  // in the same NUMA context and on the same thread as the
1384  // function is called
1385  for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
1386  {
1388  = data->get_fe_output_object(base_no);
1389  base_fe_output_object.initialize (quadrature.size(), base_element(base_no),
1390  flags | base_element(base_no).requires_update_flags(flags));
1391 
1392  // let base objects produce their scratch objects. they may
1393  // also at this time write into the output objects we provide
1394  // for them; it would be nice if we could already copy something
1395  // out of the base output object into the system output object,
1396  // but we can't because we can't know what the elements already
1397  // copied and/or will want to update on every cell
1398  typename FiniteElement<dim,spacedim>::InternalDataBase *base_fe_data =
1399  base_element(base_no).get_data (flags, mapping, quadrature,
1400  base_fe_output_object);
1401 
1402  data->set_fe_data(base_no, base_fe_data);
1403  }
1404 
1405  return data;
1406 }
1407 
1408 // The following function is a clone of get_data, with the exception
1409 // that get_face_data of the base elements is called.
1410 
1411 template <int dim, int spacedim>
1415  const Mapping<dim,spacedim> &mapping,
1416  const Quadrature<dim-1> &quadrature,
1418 {
1419  // create an internal data object and set the update flags we will need
1420  // to deal with. the current object does not make use of these flags,
1421  // but we need to nevertheless set them correctly since we look
1422  // into the update_each flag of base elements in fill_fe_values,
1423  // and so the current object's update_each flag needs to be
1424  // correct in case the current FESystem is a base element for another,
1425  // higher-level FESystem itself.
1426  InternalData *data = new InternalData(this->n_base_elements());
1427  data->update_each = requires_update_flags(flags);
1428 
1429  // get data objects from each of the base elements and store
1430  // them. one might think that doing this in parallel (over the
1431  // base elements) would be a good idea, but this turns out to
1432  // be wrong because we would then run these jobs on different
1433  // threads/processors and this allocates memory in different
1434  // NUMA domains; this has large detrimental effects when later
1435  // writing into these objects in fill_fe_*_values. all of this
1436  // is particularly true when using FEValues objects in
1437  // WorkStream contexts where we explicitly make sure that
1438  // every function only uses objects previously allocated
1439  // in the same NUMA context and on the same thread as the
1440  // function is called
1441  for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
1442  {
1444  = data->get_fe_output_object(base_no);
1445  base_fe_output_object.initialize (quadrature.size(), base_element(base_no),
1446  flags | base_element(base_no).requires_update_flags(flags));
1447 
1448  // let base objects produce their scratch objects. they may
1449  // also at this time write into the output objects we provide
1450  // for them; it would be nice if we could already copy something
1451  // out of the base output object into the system output object,
1452  // but we can't because we can't know what the elements already
1453  // copied and/or will want to update on every cell
1454  typename FiniteElement<dim,spacedim>::InternalDataBase *base_fe_data =
1455  base_element(base_no).get_face_data (flags, mapping, quadrature,
1456  base_fe_output_object);
1457 
1458  data->set_fe_data(base_no, base_fe_data);
1459  }
1460 
1461  return data;
1462 }
1463 
1464 
1465 
1466 // The following function is a clone of get_data, with the exception
1467 // that get_subface_data of the base elements is called.
1468 
1469 template <int dim, int spacedim>
1473  const Mapping<dim,spacedim> &mapping,
1474  const Quadrature<dim-1> &quadrature,
1476 {
1477  // create an internal data object and set the update flags we will need
1478  // to deal with. the current object does not make use of these flags,
1479  // but we need to nevertheless set them correctly since we look
1480  // into the update_each flag of base elements in fill_fe_values,
1481  // and so the current object's update_each flag needs to be
1482  // correct in case the current FESystem is a base element for another,
1483  // higher-level FESystem itself.
1484  InternalData *data = new InternalData(this->n_base_elements());
1485  data->update_each = requires_update_flags(flags);
1486 
1487  // get data objects from each of the base elements and store
1488  // them. one might think that doing this in parallel (over the
1489  // base elements) would be a good idea, but this turns out to
1490  // be wrong because we would then run these jobs on different
1491  // threads/processors and this allocates memory in different
1492  // NUMA domains; this has large detrimental effects when later
1493  // writing into these objects in fill_fe_*_values. all of this
1494  // is particularly true when using FEValues objects in
1495  // WorkStream contexts where we explicitly make sure that
1496  // every function only uses objects previously allocated
1497  // in the same NUMA context and on the same thread as the
1498  // function is called
1499  for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
1500  {
1502  = data->get_fe_output_object(base_no);
1503  base_fe_output_object.initialize (quadrature.size(), base_element(base_no),
1504  flags | base_element(base_no).requires_update_flags(flags));
1505 
1506  // let base objects produce their scratch objects. they may
1507  // also at this time write into the output objects we provide
1508  // for them; it would be nice if we could already copy something
1509  // out of the base output object into the system output object,
1510  // but we can't because we can't know what the elements already
1511  // copied and/or will want to update on every cell
1512  typename FiniteElement<dim,spacedim>::InternalDataBase *base_fe_data =
1513  base_element(base_no).get_subface_data (flags, mapping, quadrature,
1514  base_fe_output_object);
1515 
1516  data->set_fe_data(base_no, base_fe_data);
1517  }
1518 
1519  return data;
1520 }
1521 
1522 
1523 
1524 template <int dim, int spacedim>
1525 void
1528  const CellSimilarity::Similarity cell_similarity,
1529  const Quadrature<dim> &quadrature,
1530  const Mapping<dim,spacedim> &mapping,
1531  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
1532  const ::internal::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
1533  const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
1535 {
1536  compute_fill(mapping, cell, invalid_face_number, invalid_face_number,
1537  quadrature, cell_similarity, mapping_internal, fe_internal,
1538  mapping_data, output_data);
1539 }
1540 
1541 
1542 
1543 template <int dim, int spacedim>
1544 void
1547  const unsigned int face_no,
1548  const Quadrature<dim-1> &quadrature,
1549  const Mapping<dim,spacedim> &mapping,
1550  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
1551  const ::internal::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
1552  const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
1554 {
1555  compute_fill (mapping, cell, face_no, invalid_face_number, quadrature,
1556  CellSimilarity::none, mapping_internal, fe_internal,
1557  mapping_data, output_data);
1558 }
1559 
1560 
1561 
1562 
1563 template <int dim, int spacedim>
1564 void
1567  const unsigned int face_no,
1568  const unsigned int sub_no,
1569  const Quadrature<dim-1> &quadrature,
1570  const Mapping<dim,spacedim> &mapping,
1571  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
1572  const ::internal::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
1573  const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
1575 {
1576  compute_fill (mapping, cell, face_no, sub_no, quadrature,
1577  CellSimilarity::none, mapping_internal, fe_internal,
1578  mapping_data, output_data);
1579 }
1580 
1581 
1582 
1583 template <int dim, int spacedim>
1584 template <int dim_1>
1585 void
1588  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
1589  const unsigned int face_no,
1590  const unsigned int sub_no,
1591  const Quadrature<dim_1> &quadrature,
1592  const CellSimilarity::Similarity cell_similarity,
1593  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
1594  const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,
1597 {
1598  // convert data object to internal
1599  // data for this class. fails with
1600  // an exception if that is not
1601  // possible
1602  Assert (dynamic_cast<const InternalData *> (&fe_internal) != 0, ExcInternalError());
1603  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
1604 
1605  // Either dim_1==dim
1606  // (fill_fe_values) or dim_1==dim-1
1607  // (fill_fe_(sub)face_values)
1608  Assert(dim_1==dim || dim_1==dim-1, ExcInternalError());
1609  const UpdateFlags flags = fe_data.update_each;
1610 
1611 
1612  // loop over the base elements, let them compute what they need to compute,
1613  // and then copy what is necessary.
1614  //
1615  // one may think that it would be a good idea to parallelize this over
1616  // base elements, but it turns out to be not worthwhile: doing so lets
1617  // multiple threads access data objects that were created by the current
1618  // thread, leading to many NUMA memory access inefficiencies. we specifically
1619  // want to avoid this if this class is called in a WorkStream context where
1620  // we very carefully allocate objects only on the thread where they
1621  // will actually be used; spawning new tasks here would be counterproductive
1622  if (flags & (update_values | update_gradients
1624  for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
1625  {
1627  base_fe = base_element(base_no);
1629  base_fe_data = fe_data.get_fe_data(base_no);
1631  base_data = fe_data.get_fe_output_object(base_no);
1632 
1633  // fill_fe_face_values needs argument Quadrature<dim-1> for both cases
1634  // dim_1==dim-1 and dim_1=dim. Hence the following workaround
1635  const Quadrature<dim> *cell_quadrature = 0;
1636  const Quadrature<dim-1> *face_quadrature = 0;
1637  const unsigned int n_q_points = quadrature.size();
1638 
1639  // static cast to the common base class of quadrature being either
1640  // Quadrature<dim> or Quadrature<dim-1>:
1641  const Subscriptor *quadrature_base_pointer = &quadrature;
1642 
1643  if (face_no==invalid_face_number)
1644  {
1645  Assert(dim_1==dim, ExcDimensionMismatch(dim_1,dim));
1646  Assert (dynamic_cast<const Quadrature<dim> *>(quadrature_base_pointer) != 0,
1647  ExcInternalError());
1648 
1649  cell_quadrature
1650  = static_cast<const Quadrature<dim> *>(quadrature_base_pointer);
1651  }
1652  else
1653  {
1654  Assert(dim_1==dim-1, ExcDimensionMismatch(dim_1,dim-1));
1655  Assert (dynamic_cast<const Quadrature<dim-1> *>(quadrature_base_pointer) != 0,
1656  ExcInternalError());
1657 
1658  face_quadrature
1659  = static_cast<const Quadrature<dim-1> *>(quadrature_base_pointer);
1660  }
1661 
1662 
1663  // Make sure that in the case of fill_fe_values the data is only
1664  // copied from base_data to data if base_data is changed. therefore
1665  // use fe_fe_data.current_update_flags()
1666  //
1667  // for the case of fill_fe_(sub)face_values the data needs to be
1668  // copied from base_data to data on each face, therefore use
1669  // base_fe_data.update_flags.
1670  if (face_no==invalid_face_number)
1671  base_fe.fill_fe_values(cell, cell_similarity,
1672  *cell_quadrature,
1673  mapping, mapping_internal, mapping_data,
1674  base_fe_data, base_data);
1675  else if (sub_no==invalid_face_number)
1676  base_fe.fill_fe_face_values(cell, face_no,
1677  *face_quadrature,
1678  mapping, mapping_internal, mapping_data,
1679  base_fe_data, base_data);
1680  else
1681  base_fe.fill_fe_subface_values(cell, face_no, sub_no,
1682  *face_quadrature,
1683  mapping, mapping_internal, mapping_data,
1684  base_fe_data, base_data);
1685 
1686  // now data has been generated, so copy it. we used to work by
1687  // looping over all base elements (i.e. this outer loop), then over
1688  // multiplicity, then over the shape functions from that base
1689  // element, but that requires that we can infer the global number of
1690  // a shape function from its number in the base element. for that we
1691  // had the component_to_system_table.
1692  //
1693  // however, this does of course no longer work since we have
1694  // non-primitive elements. so we go the other way round: loop over
1695  // all shape functions of the composed element, and here only treat
1696  // those shape functions that belong to a given base element
1697  //TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style
1698  const UpdateFlags base_flags = base_fe_data.update_each;
1699 
1700  // if the current cell is just a translation of the previous one,
1701  // the underlying data has not changed, and we don't even need to
1702  // enter this section
1703  if (cell_similarity != CellSimilarity::translation)
1704  for (unsigned int system_index=0; system_index<this->dofs_per_cell;
1705  ++system_index)
1706  if (this->system_to_base_table[system_index].first.first == base_no)
1707  {
1708  const unsigned int
1709  base_index = this->system_to_base_table[system_index].second;
1710  Assert (base_index<base_fe.dofs_per_cell, ExcInternalError());
1711 
1712  // now copy. if the shape function is primitive, then there
1713  // is only one value to be copied, but for non-primitive
1714  // elements, there might be more values to be copied
1715  //
1716  // so, find out from which index to take this one value, and
1717  // to which index to put
1718  unsigned int out_index = 0;
1719  for (unsigned int i=0; i<system_index; ++i)
1720  out_index += this->n_nonzero_components(i);
1721  unsigned int in_index = 0;
1722  for (unsigned int i=0; i<base_index; ++i)
1723  in_index += base_fe.n_nonzero_components(i);
1724 
1725  // then loop over the number of components to be copied
1726  Assert (this->n_nonzero_components(system_index) ==
1727  base_fe.n_nonzero_components(base_index),
1728  ExcInternalError());
1729 
1730  if (base_flags & update_values)
1731  for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
1732  for (unsigned int q=0; q<n_q_points; ++q)
1733  output_data.shape_values[out_index+s][q] =
1734  base_data.shape_values(in_index+s,q);
1735 
1736  if (base_flags & update_gradients)
1737  for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
1738  for (unsigned int q=0; q<n_q_points; ++q)
1739  output_data.shape_gradients[out_index+s][q] =
1740  base_data.shape_gradients[in_index+s][q];
1741 
1742  if (base_flags & update_hessians)
1743  for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
1744  for (unsigned int q=0; q<n_q_points; ++q)
1745  output_data.shape_hessians[out_index+s][q] =
1746  base_data.shape_hessians[in_index+s][q];
1747 
1748  if (base_flags & update_3rd_derivatives)
1749  for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
1750  for (unsigned int q=0; q<n_q_points; ++q)
1751  output_data.shape_3rd_derivatives[out_index+s][q] =
1752  base_data.shape_3rd_derivatives[in_index+s][q];
1753 
1754  }
1755  }
1756 }
1757 
1758 
1759 
1760 template <int dim, int spacedim>
1761 void
1763 {
1764  // If the system is not primitive, these have not been initialized by
1765  // FiniteElement
1766  this->system_to_component_table.resize(this->dofs_per_cell);
1767  this->face_system_to_component_table.resize(this->dofs_per_face);
1768 
1769  unsigned int total_index = 0;
1770 
1771  for (unsigned int base=0; base < this->n_base_elements(); ++base)
1772  for (unsigned int m = 0; m < this->element_multiplicity(base); ++m)
1773  {
1774  for (unsigned int k=0; k<base_element(base).n_components(); ++k)
1775  this->component_to_base_table[total_index++]
1776  = std::make_pair(std::make_pair(base,k), m);
1777  }
1778  Assert (total_index == this->component_to_base_table.size(),
1779  ExcInternalError());
1780 
1781  // Initialize index tables. Multi-component base elements have to be
1782  // thought of. For non-primitive shape functions, have a special invalid
1783  // index.
1784  const std::pair<unsigned int, unsigned int>
1785  non_primitive_index (numbers::invalid_unsigned_int,
1787 
1788  // First enumerate vertex indices, where we first enumerate all indices on
1789  // the first vertex in the order of the base elements, then of the second
1790  // vertex, etc
1791  total_index = 0;
1792  for (unsigned int vertex_number=0;
1793  vertex_number<GeometryInfo<dim>::vertices_per_cell;
1794  ++vertex_number)
1795  {
1796  unsigned int comp_start = 0;
1797  for (unsigned int base=0; base<this->n_base_elements(); ++base)
1798  for (unsigned int m=0; m<this->element_multiplicity(base);
1799  ++m, comp_start+=base_element(base).n_components())
1800  for (unsigned int local_index = 0;
1801  local_index < base_element(base).dofs_per_vertex;
1802  ++local_index, ++total_index)
1803  {
1804  const unsigned int index_in_base
1805  = (base_element(base).dofs_per_vertex*vertex_number +
1806  local_index);
1807 
1808  this->system_to_base_table[total_index]
1809  = std::make_pair (std::make_pair(base, m), index_in_base);
1810 
1811  if (base_element(base).is_primitive(index_in_base))
1812  {
1813  const unsigned int comp_in_base
1814  = base_element(base).system_to_component_index(index_in_base).first;
1815  const unsigned int comp
1816  = comp_start + comp_in_base;
1817  const unsigned int index_in_comp
1818  = base_element(base).system_to_component_index(index_in_base).second;
1819  this->system_to_component_table[total_index]
1820  = std::make_pair (comp, index_in_comp);
1821  }
1822  else
1823  this->system_to_component_table[total_index] = non_primitive_index;
1824  }
1825  }
1826 
1827  // 2. Lines
1829  for (unsigned int line_number= 0;
1830  line_number != GeometryInfo<dim>::lines_per_cell;
1831  ++line_number)
1832  {
1833  unsigned int comp_start = 0;
1834  for (unsigned int base=0; base<this->n_base_elements(); ++base)
1835  for (unsigned int m=0; m<this->element_multiplicity(base);
1836  ++m, comp_start+=base_element(base).n_components())
1837  for (unsigned int local_index = 0;
1838  local_index < base_element(base).dofs_per_line;
1839  ++local_index, ++total_index)
1840  {
1841  const unsigned int index_in_base
1842  = (base_element(base).dofs_per_line*line_number +
1843  local_index +
1844  base_element(base).first_line_index);
1845 
1846  this->system_to_base_table[total_index]
1847  = std::make_pair (std::make_pair(base,m), index_in_base);
1848 
1849  if (base_element(base).is_primitive(index_in_base))
1850  {
1851  const unsigned int comp_in_base
1852  = base_element(base).system_to_component_index(index_in_base).first;
1853  const unsigned int comp
1854  = comp_start + comp_in_base;
1855  const unsigned int index_in_comp
1856  = base_element(base).system_to_component_index(index_in_base).second;
1857  this->system_to_component_table[total_index]
1858  = std::make_pair (comp, index_in_comp);
1859  }
1860  else
1861  this->system_to_component_table[total_index] = non_primitive_index;
1862  }
1863  }
1864 
1865  // 3. Quads
1867  for (unsigned int quad_number= 0;
1868  quad_number != GeometryInfo<dim>::quads_per_cell;
1869  ++quad_number)
1870  {
1871  unsigned int comp_start = 0;
1872  for (unsigned int base=0; base<this->n_base_elements(); ++base)
1873  for (unsigned int m=0; m<this->element_multiplicity(base);
1874  ++m, comp_start += base_element(base).n_components())
1875  for (unsigned int local_index = 0;
1876  local_index < base_element(base).dofs_per_quad;
1877  ++local_index, ++total_index)
1878  {
1879  const unsigned int index_in_base
1880  = (base_element(base).dofs_per_quad*quad_number +
1881  local_index +
1882  base_element(base).first_quad_index);
1883 
1884  this->system_to_base_table[total_index]
1885  = std::make_pair (std::make_pair(base,m), index_in_base);
1886 
1887  if (base_element(base).is_primitive(index_in_base))
1888  {
1889  const unsigned int comp_in_base
1890  = base_element(base).system_to_component_index(index_in_base).first;
1891  const unsigned int comp
1892  = comp_start + comp_in_base;
1893  const unsigned int index_in_comp
1894  = base_element(base).system_to_component_index(index_in_base).second;
1895  this->system_to_component_table[total_index]
1896  = std::make_pair (comp, index_in_comp);
1897  }
1898  else
1899  this->system_to_component_table[total_index] = non_primitive_index;
1900  }
1901  }
1902 
1903  // 4. Hexes
1905  for (unsigned int hex_number= 0;
1906  hex_number != GeometryInfo<dim>::hexes_per_cell;
1907  ++hex_number)
1908  {
1909  unsigned int comp_start = 0;
1910  for (unsigned int base=0; base<this->n_base_elements(); ++base)
1911  for (unsigned int m=0; m<this->element_multiplicity(base);
1912  ++m, comp_start+=base_element(base).n_components())
1913  for (unsigned int local_index = 0;
1914  local_index < base_element(base).dofs_per_hex;
1915  ++local_index, ++total_index)
1916  {
1917  const unsigned int index_in_base
1918  = (base_element(base).dofs_per_hex*hex_number +
1919  local_index +
1920  base_element(base).first_hex_index);
1921 
1922  this->system_to_base_table[total_index]
1923  = std::make_pair (std::make_pair(base,m), index_in_base);
1924 
1925  if (base_element(base).is_primitive(index_in_base))
1926  {
1927  const unsigned int comp_in_base
1928  = base_element(base).system_to_component_index(index_in_base).first;
1929  const unsigned int comp
1930  = comp_start + comp_in_base;
1931  const unsigned int index_in_comp
1932  = base_element(base).system_to_component_index(index_in_base).second;
1933  this->system_to_component_table[total_index]
1934  = std::make_pair (comp, index_in_comp);
1935  }
1936  else
1937  this->system_to_component_table[total_index] = non_primitive_index;
1938  }
1939  }
1940 }
1941 
1942 
1943 
1944 template <int dim, int spacedim>
1945 void
1947 {
1948  // Initialize index tables. do this in the same way as done for the cell
1949  // tables, except that we now loop over the objects of faces
1950 
1951  // For non-primitive shape functions, have a special invalid index
1952  const std::pair<unsigned int, unsigned int>
1953  non_primitive_index (numbers::invalid_unsigned_int,
1955 
1956  // 1. Vertices
1957  unsigned int total_index = 0;
1958  for (unsigned int vertex_number=0;
1959  vertex_number<GeometryInfo<dim>::vertices_per_face;
1960  ++vertex_number)
1961  {
1962  unsigned int comp_start = 0;
1963  for (unsigned int base=0; base<this->n_base_elements(); ++base)
1964  for (unsigned int m=0; m<this->element_multiplicity(base);
1965  ++m, comp_start += base_element(base).n_components())
1966  for (unsigned int local_index = 0;
1967  local_index < base_element(base).dofs_per_vertex;
1968  ++local_index, ++total_index)
1969  {
1970  // get (cell) index of this shape function inside the base
1971  // element to see whether the shape function is primitive
1972  // (assume that all shape functions on vertices share the same
1973  // primitivity property; assume likewise for all shape functions
1974  // located on lines, quads, etc. this way, we can ask for
1975  // primitivity of only _one_ shape function, which is taken as
1976  // representative for all others located on the same type of
1977  // object):
1978  const unsigned int index_in_base
1979  = (base_element(base).dofs_per_vertex*vertex_number +
1980  local_index);
1981 
1982  const unsigned int face_index_in_base
1983  = (base_element(base).dofs_per_vertex*vertex_number +
1984  local_index);
1985 
1986  this->face_system_to_base_table[total_index]
1987  = std::make_pair (std::make_pair(base,m), face_index_in_base);
1988 
1989  if (base_element(base).is_primitive(index_in_base))
1990  {
1991  const unsigned int comp_in_base
1992  = base_element(base).face_system_to_component_index(face_index_in_base).first;
1993  const unsigned int comp
1994  = comp_start + comp_in_base;
1995  const unsigned int face_index_in_comp
1996  = base_element(base).face_system_to_component_index(face_index_in_base).second;
1997  this->face_system_to_component_table[total_index]
1998  = std::make_pair (comp, face_index_in_comp);
1999  }
2000  else
2001  this->face_system_to_component_table[total_index] = non_primitive_index;
2002  }
2003  }
2004 
2005  // 2. Lines
2007  for (unsigned int line_number= 0;
2008  line_number != GeometryInfo<dim>::lines_per_face;
2009  ++line_number)
2010  {
2011  unsigned int comp_start = 0;
2012  for (unsigned int base = 0; base < this->n_base_elements(); ++base)
2013  for (unsigned int m=0; m<this->element_multiplicity(base);
2014  ++m, comp_start += base_element(base).n_components())
2015  for (unsigned int local_index = 0;
2016  local_index < base_element(base).dofs_per_line;
2017  ++local_index, ++total_index)
2018  {
2019  // do everything alike for this type of object
2020  const unsigned int index_in_base
2021  = (base_element(base).dofs_per_line*line_number +
2022  local_index +
2023  base_element(base).first_line_index);
2024 
2025  const unsigned int face_index_in_base
2026  = (base_element(base).first_face_line_index +
2027  base_element(base).dofs_per_line * line_number +
2028  local_index);
2029 
2030  this->face_system_to_base_table[total_index]
2031  = std::make_pair (std::make_pair(base,m), face_index_in_base);
2032 
2033  if (base_element(base).is_primitive(index_in_base))
2034  {
2035  const unsigned int comp_in_base
2036  = base_element(base).face_system_to_component_index(face_index_in_base).first;
2037  const unsigned int comp
2038  = comp_start + comp_in_base;
2039  const unsigned int face_index_in_comp
2040  = base_element(base).face_system_to_component_index(face_index_in_base).second;
2041  this->face_system_to_component_table[total_index]
2042  = std::make_pair (comp, face_index_in_comp);
2043  }
2044  else
2045  this->face_system_to_component_table[total_index] = non_primitive_index;
2046  }
2047  }
2048 
2049  // 3. Quads
2051  for (unsigned int quad_number= 0;
2052  quad_number != GeometryInfo<dim>::quads_per_face;
2053  ++quad_number)
2054  {
2055  unsigned int comp_start = 0;
2056  for (unsigned int base=0; base<this->n_base_elements(); ++base)
2057  for (unsigned int m=0; m<this->element_multiplicity(base);
2058  ++m, comp_start += base_element(base).n_components())
2059  for (unsigned int local_index = 0;
2060  local_index < base_element(base).dofs_per_quad;
2061  ++local_index, ++total_index)
2062  {
2063  // do everything alike for this type of object
2064  const unsigned int index_in_base
2065  = (base_element(base).dofs_per_quad*quad_number +
2066  local_index +
2067  base_element(base).first_quad_index);
2068 
2069  const unsigned int face_index_in_base
2070  = (base_element(base).first_face_quad_index +
2071  base_element(base).dofs_per_quad * quad_number +
2072  local_index);
2073 
2074  this->face_system_to_base_table[total_index]
2075  = std::make_pair (std::make_pair(base,m), face_index_in_base);
2076 
2077  if (base_element(base).is_primitive(index_in_base))
2078  {
2079  const unsigned int comp_in_base
2080  = base_element(base).face_system_to_component_index(face_index_in_base).first;
2081  const unsigned int comp
2082  = comp_start + comp_in_base;
2083  const unsigned int face_index_in_comp
2084  = base_element(base).face_system_to_component_index(face_index_in_base).second;
2085  this->face_system_to_component_table[total_index]
2086  = std::make_pair (comp, face_index_in_comp);
2087  }
2088  else
2089  this->face_system_to_component_table[total_index] = non_primitive_index;
2090  }
2091  }
2092  Assert (total_index == this->dofs_per_face, ExcInternalError());
2093  Assert (total_index == this->face_system_to_component_table.size(),
2094  ExcInternalError());
2095  Assert (total_index == this->face_system_to_base_table.size(),
2096  ExcInternalError());
2097 }
2098 
2099 
2100 
2101 template <int dim, int spacedim>
2103 {
2104  // check whether all base elements implement their interface constraint
2105  // matrices. if this is not the case, then leave the interface costraints of
2106  // this composed element empty as well; however, the rest of the element is
2107  // usable
2108  for (unsigned int base=0; base<this->n_base_elements(); ++base)
2109  if (base_element(base).constraints_are_implemented() == false)
2110  return;
2111 
2112  this->interface_constraints.
2113  TableBase<2,double>::reinit (this->interface_constraints_size());
2114 
2115  // the layout of the constraints matrix is described in the FiniteElement
2116  // class. you may want to look there first before trying to understand the
2117  // following, especially the mapping of the @p{m} index.
2118  //
2119  // in order to map it to the fe-system class, we have to know which base
2120  // element a degree of freedom within a vertex, line, etc belongs to. this
2121  // can be accomplished by the system_to_component_index function in
2122  // conjunction with the numbers first_{line,quad,...}_index
2123  for (unsigned int n=0; n<this->interface_constraints.n(); ++n)
2124  for (unsigned int m=0; m<this->interface_constraints.m(); ++m)
2125  {
2126  // for the pair (n,m) find out which base element they belong to and
2127  // the number therein
2128  //
2129  // first for the n index. this is simple since the n indices are in
2130  // the same order as they are usually on a face. note that for the
2131  // data type, first value in pair is (base element,instance of base
2132  // element), second is index within this instance
2133  const std::pair<std::pair<unsigned int,unsigned int>, unsigned int> n_index
2134  = this->face_system_to_base_table[n];
2135 
2136  // likewise for the m index. this is more complicated due to the
2137  // strange ordering we have for the dofs on the refined faces.
2138  std::pair<std::pair<unsigned int,unsigned int>, unsigned int> m_index;
2139  switch (dim)
2140  {
2141  case 1:
2142  {
2143  // we should never get here! (in 1d, the constraints matrix
2144  // should be of size zero)
2145  Assert (false, ExcInternalError());
2146  break;
2147  };
2148 
2149  case 2:
2150  {
2151  // the indices m=0..d_v-1 are from the center vertex. their order
2152  // is the same as for the first vertex of the whole cell, so we
2153  // can use the system_to_base_table variable (using the
2154  // face_s_t_base_t function would yield the same)
2155  if (m < this->dofs_per_vertex)
2156  m_index = this->system_to_base_table[m];
2157  else
2158  // then come the two sets of line indices
2159  {
2160  const unsigned int index_in_line
2161  = (m-this->dofs_per_vertex) % this->dofs_per_line;
2162  const unsigned int sub_line
2163  = (m-this->dofs_per_vertex) / this->dofs_per_line;
2164  Assert (sub_line < 2, ExcInternalError());
2165 
2166  // from this information, try to get base element and instance
2167  // of base element. we do so by constructing the corresponding
2168  // face index of m in the present element, then use
2169  // face_system_to_base_table
2170  const unsigned int tmp1 = 2*this->dofs_per_vertex+index_in_line;
2171  m_index.first = this->face_system_to_base_table[tmp1].first;
2172 
2173  // what we are still missing is the index of m within the base
2174  // elements interface_constraints table
2175  //
2176  // here, the second value of face_system_to_base_table can
2177  // help: it denotes the face index of that shape function
2178  // within the base element. since we know that it is a line
2179  // dof, we can construct the rest: tmp2 will denote the index
2180  // of this shape function among the line shape functions:
2181  Assert (this->face_system_to_base_table[tmp1].second >=
2182  2*base_element(m_index.first.first).dofs_per_vertex,
2183  ExcInternalError());
2184  const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second -
2185  2*base_element(m_index.first.first).dofs_per_vertex;
2186  Assert (tmp2 < base_element(m_index.first.first).dofs_per_line,
2187  ExcInternalError());
2188  m_index.second = base_element(m_index.first.first).dofs_per_vertex +
2189  base_element(m_index.first.first).dofs_per_line*sub_line +
2190  tmp2;
2191  };
2192  break;
2193  };
2194 
2195  case 3:
2196  {
2197  // same way as above, although a little more complicated...
2198 
2199  // the indices m=0..5*d_v-1 are from the center and the four
2200  // subline vertices. their order is the same as for the first
2201  // vertex of the whole cell, so we can use the simple arithmetic
2202  if (m < 5*this->dofs_per_vertex)
2203  m_index = this->system_to_base_table[m];
2204  else
2205  // then come the 12 sets of line indices
2206  if (m < 5*this->dofs_per_vertex + 12*this->dofs_per_line)
2207  {
2208  // for the meaning of all this, see the 2d part
2209  const unsigned int index_in_line
2210  = (m-5*this->dofs_per_vertex) % this->dofs_per_line;
2211  const unsigned int sub_line
2212  = (m-5*this->dofs_per_vertex) / this->dofs_per_line;
2213  Assert (sub_line < 12, ExcInternalError());
2214 
2215  const unsigned int tmp1 = 4*this->dofs_per_vertex+index_in_line;
2216  m_index.first = this->face_system_to_base_table[tmp1].first;
2217 
2218  Assert (this->face_system_to_base_table[tmp1].second >=
2219  4*base_element(m_index.first.first).dofs_per_vertex,
2220  ExcInternalError());
2221  const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second -
2222  4*base_element(m_index.first.first).dofs_per_vertex;
2223  Assert (tmp2 < base_element(m_index.first.first).dofs_per_line,
2224  ExcInternalError());
2225  m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex +
2226  base_element(m_index.first.first).dofs_per_line*sub_line +
2227  tmp2;
2228  }
2229  else
2230  // on one of the four sub-quads
2231  {
2232  // for the meaning of all this, see the 2d part
2233  const unsigned int index_in_quad
2234  = (m-5*this->dofs_per_vertex-12*this->dofs_per_line) %
2235  this->dofs_per_quad;
2236  Assert (index_in_quad < this->dofs_per_quad,
2237  ExcInternalError());
2238  const unsigned int sub_quad
2239  = ((m-5*this->dofs_per_vertex-12*this->dofs_per_line) /
2240  this->dofs_per_quad);
2241  Assert (sub_quad < 4, ExcInternalError());
2242 
2243  const unsigned int tmp1 = 4*this->dofs_per_vertex +
2244  4*this->dofs_per_line +
2245  index_in_quad;
2246  Assert (tmp1 < this->face_system_to_base_table.size(),
2247  ExcInternalError());
2248  m_index.first = this->face_system_to_base_table[tmp1].first;
2249 
2250  Assert (this->face_system_to_base_table[tmp1].second >=
2251  4*base_element(m_index.first.first).dofs_per_vertex +
2252  4*base_element(m_index.first.first).dofs_per_line,
2253  ExcInternalError());
2254  const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second -
2255  4*base_element(m_index.first.first).dofs_per_vertex -
2256  4*base_element(m_index.first.first).dofs_per_line;
2257  Assert (tmp2 < base_element(m_index.first.first).dofs_per_quad,
2258  ExcInternalError());
2259  m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex +
2260  12*base_element(m_index.first.first).dofs_per_line +
2261  base_element(m_index.first.first).dofs_per_quad*sub_quad +
2262  tmp2;
2263  };
2264 
2265  break;
2266  };
2267 
2268  default:
2269  Assert (false, ExcNotImplemented());
2270  };
2271 
2272  // now that we gathered all information: use it to build the
2273  // matrix. note that if n and m belong to different base elements or
2274  // instances, then there definitely will be no coupling
2275  if (n_index.first == m_index.first)
2276  this->interface_constraints(m,n)
2277  = (base_element(n_index.first.first).constraints()(m_index.second,
2278  n_index.second));
2279  };
2280 }
2281 
2282 
2283 
2284 template <int dim, int spacedim>
2286  const std::vector<unsigned int> &multiplicities)
2287 {
2288  Assert (fes.size() == multiplicities.size(),
2289  ExcDimensionMismatch (fes.size(), multiplicities.size()) );
2290  Assert (fes.size() > 0,
2291  ExcMessage ("Need to pass at least one finite element."));
2292  Assert (count_nonzeros(multiplicities) > 0,
2293  ExcMessage("You only passed FiniteElements with multiplicity 0."));
2294 
2295  // Note that we need to skip every fe with multiplicity 0 in the following block of code
2296 
2297  this->base_to_block_indices.reinit(0, 0);
2298 
2299  for (unsigned int i=0; i<fes.size(); i++)
2300  if (multiplicities[i]>0)
2301  this->base_to_block_indices.push_back( multiplicities[i] );
2302 
2303  std::vector<Threads::Task<FiniteElement<dim,spacedim>*> > clone_base_elements;
2304 
2305  for (unsigned int i=0; i<fes.size(); i++)
2306  if (multiplicities[i]>0)
2307  clone_base_elements.push_back (Threads::new_task (&FiniteElement<dim,spacedim>::clone,
2308  *fes[i]));
2309 
2310  unsigned int ind=0;
2311  for (unsigned int i=0; i<fes.size(); i++)
2312  {
2313  if (multiplicities[i]>0)
2314  {
2315  base_elements[ind] =
2316  std::make_pair (std_cxx11::shared_ptr<const FiniteElement<dim,spacedim> >
2317  (clone_base_elements[ind].return_value()),
2318  multiplicities[i]);
2319  ++ind;
2320  }
2321  }
2322 
2323  Assert(ind>0, ExcInternalError());
2324 
2325  build_cell_tables();
2326  build_face_tables();
2327 
2328  // restriction and prolongation matrices are build on demand
2329 
2330  // now set up the interface constraints. this is kind'o hairy, so don't try
2331  // to do it dimension independent
2332  build_interface_constraints ();
2333 
2334  // finally fill in support points on cell and face
2335  initialize_unit_support_points ();
2336  initialize_unit_face_support_points ();
2337 
2338  initialize_quad_dof_index_permutation ();
2339 }
2340 
2341 
2342 
2343 
2344 
2345 
2346 template <int dim, int spacedim>
2347 void
2350 {
2351  // if one of the base elements has no support points, then it makes no sense
2352  // to define support points for the composed element, so return an empty
2353  // array to demonstrate that fact. Note that we ignore FE_Nothing in this logic.
2354  for (unsigned int base_el=0; base_el<this->n_base_elements(); ++base_el)
2355  if (!base_element(base_el).has_support_points() && base_element(base_el).dofs_per_cell!=0)
2356  {
2357  this->unit_support_points.resize(0);
2358  return;
2359  };
2360 
2361  // generate unit support points from unit support points of sub elements
2362  this->unit_support_points.resize(this->dofs_per_cell);
2363 
2364  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
2365  {
2366  const unsigned int
2367  base = this->system_to_base_table[i].first.first,
2368  base_index = this->system_to_base_table[i].second;
2369  Assert (base<this->n_base_elements(), ExcInternalError());
2370  Assert (base_index<base_element(base).unit_support_points.size(),
2371  ExcInternalError());
2372  this->unit_support_points[i] = base_element(base).unit_support_points[base_index];
2373  };
2374 }
2375 
2376 
2377 
2378 template <int dim, int spacedim>
2379 void
2382 {
2383  // Nothing to do in 1D
2384  if (dim == 1)
2385  return;
2386 
2387  // if one of the base elements has no support points, then it makes no sense
2388  // to define support points for the composed element. In that case, return
2389  // an empty array to demonstrate that fact (note that we ask whether the
2390  // base element has no support points at all, not only none on the face!)
2391  //
2392  // on the other hand, if there is an element that simply has no degrees of
2393  // freedom on the face at all, then we don't care whether it has support
2394  // points or not. this is, for example, the case for the stable Stokes
2395  // element Q(p)^dim \times DGP(p-1).
2396  for (unsigned int base_el=0; base_el<this->n_base_elements(); ++base_el)
2397  if (!base_element(base_el).has_support_points()
2398  &&
2399  (base_element(base_el).dofs_per_face > 0))
2400  {
2401  this->unit_face_support_points.resize(0);
2402  return;
2403  }
2404 
2405 
2406  // generate unit face support points from unit support points of sub
2407  // elements
2408  this->unit_face_support_points.resize(this->dofs_per_face);
2409 
2410  for (unsigned int i=0; i<this->dofs_per_face; ++i)
2411  {
2412  const unsigned int base_i = this->face_system_to_base_table[i].first.first;
2413  const unsigned int index_in_base = this->face_system_to_base_table[i].second;
2414 
2415  Assert (index_in_base < base_element(base_i).unit_face_support_points.size(),
2416  ExcInternalError());
2417 
2418  this->unit_face_support_points[i]
2419  = base_element(base_i).unit_face_support_points[index_in_base];
2420  }
2421 }
2422 
2423 
2424 
2425 template <int dim, int spacedim>
2426 void
2428 {
2429  // nothing to do in other dimensions than 3
2430  if (dim < 3)
2431  return;
2432 
2433  // the array into which we want to write should have the correct size
2434  // already.
2435  Assert (this->adjust_quad_dof_index_for_face_orientation_table.n_elements()==
2436  8*this->dofs_per_quad, ExcInternalError());
2437 
2438  // to obtain the shifts for this composed element, copy the shift
2439  // information of the base elements
2440  unsigned int index = 0;
2441  for (unsigned int b=0; b<this->n_base_elements(); ++b)
2442  {
2443  const Table<2,int> &temp
2444  = this->base_element(b).adjust_quad_dof_index_for_face_orientation_table;
2445  for (unsigned int c=0; c<this->element_multiplicity(b); ++c)
2446  {
2447  for (unsigned int i=0; i<temp.size(0); ++i)
2448  for (unsigned int j=0; j<8; ++j)
2449  this->adjust_quad_dof_index_for_face_orientation_table(index+i,j)=
2450  temp(i,j);
2451  index += temp.size(0);
2452  }
2453  }
2454  Assert (index == this->dofs_per_quad,
2455  ExcInternalError());
2456 
2457  // aditionally compose the permutation information for lines
2458  Assert (this->adjust_line_dof_index_for_line_orientation_table.size()==
2459  this->dofs_per_line, ExcInternalError());
2460  index = 0;
2461  for (unsigned int b=0; b<this->n_base_elements(); ++b)
2462  {
2463  const std::vector<int> &temp2
2464  = this->base_element(b).adjust_line_dof_index_for_line_orientation_table;
2465  for (unsigned int c=0; c<this->element_multiplicity(b); ++c)
2466  {
2467  std::copy(temp2.begin(), temp2.end(),
2468  this->adjust_line_dof_index_for_line_orientation_table.begin()
2469  +index);
2470  index += temp2.size();
2471  }
2472  }
2473  Assert (index == this->dofs_per_line,
2474  ExcInternalError());
2475 }
2476 
2477 
2478 
2479 template <int dim, int spacedim>
2480 bool
2483 {
2484  for (unsigned int b=0; b<this->n_base_elements(); ++b)
2485  if (base_element(b).hp_constraints_are_implemented() == false)
2486  return false;
2487 
2488  return true;
2489 }
2490 
2491 
2492 
2493 template <int dim, int spacedim>
2494 void
2497  FullMatrix<double> &interpolation_matrix) const
2498 {
2499  typedef FiniteElement<dim,spacedim> FEL;
2500  AssertThrow ((x_source_fe.get_name().find ("FE_System<") == 0)
2501  ||
2502  (dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe) != 0),
2503  typename FEL::
2504  ExcInterpolationNotImplemented());
2505 
2506  Assert (interpolation_matrix.n() == this->dofs_per_face,
2507  ExcDimensionMismatch (interpolation_matrix.n(),
2508  this->dofs_per_face));
2509  Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
2510  ExcDimensionMismatch (interpolation_matrix.m(),
2511  x_source_fe.dofs_per_face));
2512 
2513  // since dofs for each base are independent, we only have to stack things up
2514  // from base element to base element
2515  //
2516  // the problem is that we have to work with two FEs (this and
2517  // fe_other). only deal with the case that both are FESystems and that they
2518  // both have the same number of bases (counting multiplicity) each of which
2519  // match in their number of components. this covers
2520  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2521  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2522  // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2523  const FESystem<dim,spacedim> *fe_other_system
2524  = dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe);
2525 
2526  // clear matrix, since we will not get to set all elements
2527  interpolation_matrix = 0;
2528 
2529  // loop over all the base elements of this and the other element, counting
2530  // their multiplicities
2531  unsigned int base_index = 0,
2532  base_index_other = 0;
2533  unsigned int multiplicity = 0,
2534  multiplicity_other = 0;
2535 
2536  FullMatrix<double> base_to_base_interpolation;
2537 
2538  while (true)
2539  {
2541  &base = base_element(base_index),
2542  &base_other = fe_other_system->base_element(base_index_other);
2543 
2544  Assert (base.n_components() == base_other.n_components(),
2545  ExcNotImplemented());
2546 
2547  // get the interpolation from the bases
2548  base_to_base_interpolation.reinit (base_other.dofs_per_face,
2549  base.dofs_per_face);
2550  base.get_face_interpolation_matrix (base_other,
2551  base_to_base_interpolation);
2552 
2553  // now translate entries. we'd like to have something like
2554  // face_base_to_system_index, but that doesn't exist. rather, all we
2555  // have is the reverse. well, use that then
2556  for (unsigned int i=0; i<this->dofs_per_face; ++i)
2557  if (this->face_system_to_base_index(i).first
2558  ==
2559  std::make_pair (base_index, multiplicity))
2560  for (unsigned int j=0; j<fe_other_system->dofs_per_face; ++j)
2561  if (fe_other_system->face_system_to_base_index(j).first
2562  ==
2563  std::make_pair (base_index_other, multiplicity_other))
2564  interpolation_matrix(j, i)
2565  = base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second,
2566  this->face_system_to_base_index(i).second);
2567 
2568  // advance to the next base element for this and the other fe_system;
2569  // see if we can simply advance the multiplicity by one, or if have to
2570  // move on to the next base element
2571  ++multiplicity;
2572  if (multiplicity == this->element_multiplicity(base_index))
2573  {
2574  multiplicity = 0;
2575  ++base_index;
2576  }
2577  ++multiplicity_other;
2578  if (multiplicity_other ==
2579  fe_other_system->element_multiplicity(base_index_other))
2580  {
2581  multiplicity_other = 0;
2582  ++base_index_other;
2583  }
2584 
2585  // see if we have reached the end of the present element. if so, we
2586  // should have reached the end of the other one as well
2587  if (base_index == this->n_base_elements())
2588  {
2589  Assert (base_index_other == fe_other_system->n_base_elements(),
2590  ExcInternalError());
2591  break;
2592  }
2593 
2594  // if we haven't reached the end of this element, we shouldn't have
2595  // reached the end of the other one either
2596  Assert (base_index_other != fe_other_system->n_base_elements(),
2597  ExcInternalError());
2598  }
2599 }
2600 
2601 
2602 
2603 template <int dim, int spacedim>
2604 void
2607  const unsigned int subface,
2608  FullMatrix<double> &interpolation_matrix) const
2609 {
2610  typedef FiniteElement<dim,spacedim> FEL;
2611  AssertThrow ((x_source_fe.get_name().find ("FE_System<") == 0)
2612  ||
2613  (dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe) != 0),
2614  typename FEL::
2615  ExcInterpolationNotImplemented());
2616 
2617  Assert (interpolation_matrix.n() == this->dofs_per_face,
2618  ExcDimensionMismatch (interpolation_matrix.n(),
2619  this->dofs_per_face));
2620  Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
2621  ExcDimensionMismatch (interpolation_matrix.m(),
2622  x_source_fe.dofs_per_face));
2623 
2624  // since dofs for each base are independent, we only have to stack things up
2625  // from base element to base element
2626  //
2627  // the problem is that we have to work with two FEs (this and
2628  // fe_other). only deal with the case that both are FESystems and that they
2629  // both have the same number of bases (counting multiplicity) each of which
2630  // match in their number of components. this covers
2631  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2632  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2633  // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2634  const FESystem<dim,spacedim> *fe_other_system
2635  = dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe);
2636 
2637  // clear matrix, since we will not get to set all elements
2638  interpolation_matrix = 0;
2639 
2640  // loop over all the base elements of this and the other element, counting
2641  // their multiplicities
2642  unsigned int base_index = 0,
2643  base_index_other = 0;
2644  unsigned int multiplicity = 0,
2645  multiplicity_other = 0;
2646 
2647  FullMatrix<double> base_to_base_interpolation;
2648 
2649  while (true)
2650  {
2652  &base = base_element(base_index),
2653  &base_other = fe_other_system->base_element(base_index_other);
2654 
2655  Assert (base.n_components() == base_other.n_components(),
2656  ExcNotImplemented());
2657 
2658  // get the interpolation from the bases
2659  base_to_base_interpolation.reinit (base_other.dofs_per_face,
2660  base.dofs_per_face);
2661  base.get_subface_interpolation_matrix (base_other,
2662  subface,
2663  base_to_base_interpolation);
2664 
2665  // now translate entries. we'd like to have something like
2666  // face_base_to_system_index, but that doesn't exist. rather, all we
2667  // have is the reverse. well, use that then
2668  for (unsigned int i=0; i<this->dofs_per_face; ++i)
2669  if (this->face_system_to_base_index(i).first
2670  ==
2671  std::make_pair (base_index, multiplicity))
2672  for (unsigned int j=0; j<fe_other_system->dofs_per_face; ++j)
2673  if (fe_other_system->face_system_to_base_index(j).first
2674  ==
2675  std::make_pair (base_index_other, multiplicity_other))
2676  interpolation_matrix(j, i)
2677  = base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second,
2678  this->face_system_to_base_index(i).second);
2679 
2680  // advance to the next base element for this and the other fe_system;
2681  // see if we can simply advance the multiplicity by one, or if have to
2682  // move on to the next base element
2683  ++multiplicity;
2684  if (multiplicity == this->element_multiplicity(base_index))
2685  {
2686  multiplicity = 0;
2687  ++base_index;
2688  }
2689  ++multiplicity_other;
2690  if (multiplicity_other ==
2691  fe_other_system->element_multiplicity(base_index_other))
2692  {
2693  multiplicity_other = 0;
2694  ++base_index_other;
2695  }
2696 
2697  // see if we have reached the end of the present element. if so, we
2698  // should have reached the end of the other one as well
2699  if (base_index == this->n_base_elements())
2700  {
2701  Assert (base_index_other == fe_other_system->n_base_elements(),
2702  ExcInternalError());
2703  break;
2704  }
2705 
2706  // if we haven't reached the end of this element, we shouldn't have
2707  // reached the end of the other one either
2708  Assert (base_index_other != fe_other_system->n_base_elements(),
2709  ExcInternalError());
2710  }
2711 }
2712 
2713 
2714 
2715 template <int dim, int spacedim>
2716 template <int structdim>
2717 std::vector<std::pair<unsigned int, unsigned int> >
2719 {
2720  // since dofs on each subobject (vertex, line, ...) are ordered such that
2721  // first come all from the first base element all multiplicities, then
2722  // second base element all multiplicities, etc., we simply have to stack all
2723  // the identities after each other
2724  //
2725  // the problem is that we have to work with two FEs (this and
2726  // fe_other). only deal with the case that both are FESystems and that they
2727  // both have the same number of bases (counting multiplicity) each of which
2728  // match in their number of components. this covers
2729  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2730  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2731  // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2732  if (const FESystem<dim,spacedim> *fe_other_system
2733  = dynamic_cast<const FESystem<dim,spacedim>*>(&fe_other))
2734  {
2735  // loop over all the base elements of this and the other element,
2736  // counting their multiplicities
2737  unsigned int base_index = 0,
2738  base_index_other = 0;
2739  unsigned int multiplicity = 0,
2740  multiplicity_other = 0;
2741 
2742  // we also need to keep track of the number of dofs already treated for
2743  // each of the elements
2744  unsigned int dof_offset = 0,
2745  dof_offset_other = 0;
2746 
2747  std::vector<std::pair<unsigned int, unsigned int> > identities;
2748 
2749  while (true)
2750  {
2752  &base = base_element(base_index),
2753  &base_other = fe_other_system->base_element(base_index_other);
2754 
2755  Assert (base.n_components() == base_other.n_components(),
2756  ExcNotImplemented());
2757 
2758  // now translate the identities returned by the base elements to the
2759  // indices of this system element
2760  std::vector<std::pair<unsigned int, unsigned int> > base_identities;
2761  switch (structdim)
2762  {
2763  case 0:
2764  base_identities = base.hp_vertex_dof_identities (base_other);
2765  break;
2766  case 1:
2767  base_identities = base.hp_line_dof_identities (base_other);
2768  break;
2769  case 2:
2770  base_identities = base.hp_quad_dof_identities (base_other);
2771  break;
2772  default:
2773  Assert (false, ExcNotImplemented());
2774  }
2775 
2776  for (unsigned int i=0; i<base_identities.size(); ++i)
2777  identities.push_back
2778  (std::make_pair (base_identities[i].first + dof_offset,
2779  base_identities[i].second + dof_offset_other));
2780 
2781  // record the dofs treated above as already taken care of
2782  dof_offset += base.template n_dofs_per_object<structdim>();
2783  dof_offset_other += base_other.template n_dofs_per_object<structdim>();
2784 
2785  // advance to the next base element for this and the other
2786  // fe_system; see if we can simply advance the multiplicity by one,
2787  // or if have to move on to the next base element
2788  ++multiplicity;
2789  if (multiplicity == this->element_multiplicity(base_index))
2790  {
2791  multiplicity = 0;
2792  ++base_index;
2793  }
2794  ++multiplicity_other;
2795  if (multiplicity_other ==
2796  fe_other_system->element_multiplicity(base_index_other))
2797  {
2798  multiplicity_other = 0;
2799  ++base_index_other;
2800  }
2801 
2802  // see if we have reached the end of the present element. if so, we
2803  // should have reached the end of the other one as well
2804  if (base_index == this->n_base_elements())
2805  {
2806  Assert (base_index_other == fe_other_system->n_base_elements(),
2807  ExcInternalError());
2808  break;
2809  }
2810 
2811  // if we haven't reached the end of this element, we shouldn't have
2812  // reached the end of the other one either
2813  Assert (base_index_other != fe_other_system->n_base_elements(),
2814  ExcInternalError());
2815  }
2816 
2817  return identities;
2818  }
2819  else
2820  {
2821  Assert (false, ExcNotImplemented());
2822  return std::vector<std::pair<unsigned int, unsigned int> >();
2823  }
2824 }
2825 
2826 
2827 
2828 template <int dim, int spacedim>
2829 std::vector<std::pair<unsigned int, unsigned int> >
2831 {
2832  return hp_object_dof_identities<0> (fe_other);
2833 }
2834 
2835 template <int dim, int spacedim>
2836 std::vector<std::pair<unsigned int, unsigned int> >
2838 {
2839  return hp_object_dof_identities<1> (fe_other);
2840 }
2841 
2842 
2843 
2844 template <int dim, int spacedim>
2845 std::vector<std::pair<unsigned int, unsigned int> >
2847 {
2848  return hp_object_dof_identities<2> (fe_other);
2849 }
2850 
2851 
2852 
2853 template <int dim, int spacedim>
2857 {
2858  // at present all we can do is to compare with other FESystems that have the
2859  // same number of components and bases
2860  if (const FESystem<dim,spacedim> *fe_sys_other
2861  = dynamic_cast<const FESystem<dim,spacedim>*>(&fe_other))
2862  {
2863  Assert (this->n_components() == fe_sys_other->n_components(),
2864  ExcNotImplemented());
2865  Assert (this->n_base_elements() == fe_sys_other->n_base_elements(),
2866  ExcNotImplemented());
2867 
2869  domination = FiniteElementDomination::no_requirements;
2870 
2871  // loop over all base elements and do some sanity checks
2872  for (unsigned int b=0; b<this->n_base_elements(); ++b)
2873  {
2874  Assert (this->base_element(b).n_components() ==
2875  fe_sys_other->base_element(b).n_components(),
2876  ExcNotImplemented());
2877  Assert (this->element_multiplicity(b) ==
2878  fe_sys_other->element_multiplicity(b),
2879  ExcNotImplemented());
2880 
2881  // for this pair of base elements, check who dominates and combine
2882  // with previous result
2884  base_domination = (this->base_element(b)
2885  .compare_for_face_domination (fe_sys_other->base_element(b)));
2886  domination = domination & base_domination;
2887  }
2888 
2889  return domination;
2890  }
2891 
2892  Assert (false, ExcNotImplemented());
2893  return FiniteElementDomination::neither_element_dominates;
2894 }
2895 
2896 
2897 
2898 template <int dim, int spacedim>
2900 FESystem<dim,spacedim>::base_element (const unsigned int index) const
2901 {
2902  Assert (index < base_elements.size(),
2903  ExcIndexRange(index, 0, base_elements.size()));
2904  return *base_elements[index].first;
2905 }
2906 
2907 
2908 
2909 template <int dim, int spacedim>
2910 bool
2911 FESystem<dim,spacedim>::has_support_on_face (const unsigned int shape_index,
2912  const unsigned int face_index) const
2913 {
2914  return (base_element(this->system_to_base_index(shape_index).first.first)
2915  .has_support_on_face(this->system_to_base_index(shape_index).second,
2916  face_index));
2917 }
2918 
2919 
2920 
2921 template <int dim, int spacedim>
2922 Point<dim>
2923 FESystem<dim,spacedim>::unit_support_point (const unsigned int index) const
2924 {
2925  Assert (index < this->dofs_per_cell,
2926  ExcIndexRange (index, 0, this->dofs_per_cell));
2927  typedef FiniteElement<dim,spacedim> FEL;
2928  Assert ((this->unit_support_points.size() == this->dofs_per_cell) ||
2929  (this->unit_support_points.size() == 0),
2930  typename FEL::ExcFEHasNoSupportPoints ());
2931 
2932  // let's see whether we have the information pre-computed
2933  if (this->unit_support_points.size() != 0)
2934  return this->unit_support_points[index];
2935  else
2936  // no. ask the base element whether it would like to provide this
2937  // information
2938  return (base_element(this->system_to_base_index(index).first.first)
2939  .unit_support_point(this->system_to_base_index(index).second));
2940 }
2941 
2942 
2943 
2944 template <int dim, int spacedim>
2945 Point<dim-1>
2946 FESystem<dim,spacedim>::unit_face_support_point (const unsigned int index) const
2947 {
2948  Assert (index < this->dofs_per_face,
2949  ExcIndexRange (index, 0, this->dofs_per_face));
2950  typedef FiniteElement<dim,spacedim> FEL;
2951  Assert ((this->unit_face_support_points.size() == this->dofs_per_face) ||
2952  (this->unit_face_support_points.size() == 0),
2953  typename FEL::ExcFEHasNoSupportPoints ());
2954 
2955  // let's see whether we have the information pre-computed
2956  if (this->unit_face_support_points.size() != 0)
2957  return this->unit_face_support_points[index];
2958  else
2959  // no. ask the base element whether it would like to provide this
2960  // information
2961  return (base_element(this->face_system_to_base_index(index).first.first)
2962  .unit_face_support_point(this->face_system_to_base_index(index).second));
2963 }
2964 
2965 
2966 
2967 template <int dim, int spacedim>
2968 std::pair<Table<2,bool>, std::vector<unsigned int> >
2970 {
2971  // Note that this->n_components() is actually only an estimate of how many
2972  // constant modes we will need. There might be more than one such mode
2973  // (e.g. FE_Q_DG0).
2974  Table<2,bool> constant_modes(this->n_components(), this->dofs_per_cell);
2975  std::vector<unsigned int> components;
2976  for (unsigned int i=0; i<base_elements.size(); ++i)
2977  {
2978  std::pair<Table<2,bool>, std::vector<unsigned int> >
2979  base_table = base_elements[i].first->get_constant_modes();
2980  AssertDimension(base_table.first.n_rows(), base_table.second.size());
2981  const unsigned int element_multiplicity = this->element_multiplicity(i);
2982 
2983  // there might be more than one constant mode for some scalar elements,
2984  // so make sure the table actually fits: Create a new table with more
2985  // rows
2986  const unsigned int comp = components.size();
2987  if (constant_modes.n_rows() < comp+base_table.first.n_rows()*element_multiplicity)
2988  {
2989  Table<2,bool> new_constant_modes(comp+base_table.first.n_rows()*
2990  element_multiplicity,
2991  constant_modes.n_cols());
2992  for (unsigned int r=0; r<comp; ++r)
2993  for (unsigned int c=0; c<this->dofs_per_cell; ++c)
2994  new_constant_modes(r,c) = constant_modes(r,c);
2995  constant_modes.swap(new_constant_modes);
2996  }
2997 
2998  // next, fill the constant modes from the individual components as well
2999  // as the component numbers corresponding to the constant mode rows
3000  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
3001  {
3002  std::pair<std::pair<unsigned int,unsigned int>, unsigned int> ind
3003  = this->system_to_base_index(k);
3004  if (ind.first.first == i)
3005  for (unsigned int c=0; c<base_table.first.n_rows(); ++c)
3006  constant_modes(comp+ind.first.second*base_table.first.n_rows()+c,k)
3007  = base_table.first(c,ind.second);
3008  }
3009  for (unsigned int r=0; r<element_multiplicity; ++r)
3010  for (unsigned int c=0; c<base_table.second.size(); ++c)
3011  components.push_back(comp+r*this->base_elements[i].first->n_components()
3012  +base_table.second[c]);
3013  }
3014  AssertDimension(components.size(), constant_modes.n_rows());
3015  return std::pair<Table<2,bool>, std::vector<unsigned int> >(constant_modes,
3016  components);
3017 }
3018 
3019 
3020 
3021 template <int dim, int spacedim>
3022 std::size_t
3024 {
3025  // neglect size of data stored in @p{base_elements} due to some problems
3026  // with the compiler. should be neglectable after all, considering the size
3027  // of the data of the subelements
3028  std::size_t mem = (FiniteElement<dim,spacedim>::memory_consumption () +
3029  sizeof (base_elements));
3030  for (unsigned int i=0; i<base_elements.size(); ++i)
3031  mem += MemoryConsumption::memory_consumption (*base_elements[i].first);
3032  return mem;
3033 }
3034 
3035 
3036 
3037 
3038 // explicit instantiations
3039 #include "fe_system.inst"
3040 
3041 DEAL_II_NAMESPACE_CLOSE
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe_system.cc:943
void initialize_unit_support_points()
Definition: fe_system.cc:2349
Shape function values.
virtual bool hp_constraints_are_implemented() const
Definition: fe_system.cc:2482
virtual void fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal,::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
virtual std::size_t memory_consumption() const
Definition: fe_system.cc:3023
static const unsigned int invalid_unsigned_int
Definition: types.h:164
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe_system.cc:1234
FiniteElement< dim, spacedim >::InternalDataBase & get_fe_data(const unsigned int base_no) const
Definition: fe_system.cc:73
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const
Definition: fe_system.cc:928
void build_interface_constraints()
Definition: fe_system.cc:2102
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe_system.cc:1063
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1052
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const
Definition: fe_system.cc:2969
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:913
void swap(TableBase< N, T > &v)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe_system.cc:2830
std::vector< std::pair< std_cxx11::shared_ptr< const FiniteElement< dim, spacedim > >, unsigned int > > base_elements
Definition: fe_system.h:957
virtual Tensor< 3, dim > shape_3rd_derivative(const unsigned int i, const Point< dim > &p) const
Definition: fe_system.cc:973
virtual FiniteElement< dim, spacedim > * clone() const
Definition: fe_system.cc:814
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const
Definition: fe_system.cc:2606
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe_system.cc:2856
::ExceptionBase & ExcMessage(std::string arg1)
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:54
const unsigned int degree
Definition: fe_base.h:299
void initialize_quad_dof_index_permutation()
Definition: fe_system.cc:2427
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:924
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe_system.cc:2837
void set_fe_data(const unsigned int base_no, typename FiniteElement< dim, spacedim >::InternalDataBase *)
Definition: fe_system.cc:85
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const
Definition: fe_system.cc:883
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe_system.cc:898
virtual std::string get_name() const
Definition: fe_system.cc:783
#define AssertThrow(cond, exc)
Definition: exceptions.h:358
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
Definition: fe_system.cc:1298
InternalData(const unsigned int n_base_elements)
Definition: fe_system.cc:48
void initialize(const unsigned int n_quadrature_points, const FiniteElement< dim, spacedim > &fe, const UpdateFlags flags)
Definition: fe_values.cc:2203
virtual FiniteElement< dim, spacedim >::InternalDataBase * get_subface_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim-1 > &quadrature,::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe_system.cc:1472
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:902
const BlockIndices & block_indices() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
Definition: fe_system.cc:2911
unsigned int tensor_degree() const
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:2766
std::vector< std::pair< unsigned int, unsigned int > > hp_object_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe_system.cc:2718
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe_system.cc:1145
void build_face_tables()
Definition: fe_system.cc:1946
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const
Definition: fe_system.cc:1343
unsigned int size() const
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe_system.cc:846
size_type n() const
No update.
Third derivatives of shape functions.
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:294
UpdateFlags
unsigned int n_nonzero_components(const unsigned int i) const
Definition: fe.h:2926
Abstract base class for mapping classes.
Definition: dof_tools.h:52
Definition: fe.h:34
virtual Point< dim > unit_support_point(const unsigned int index) const
Definition: fe_system.cc:2923
void initialize_unit_face_support_points()
Definition: fe_system.cc:2381
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim-1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal,::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
virtual Tensor< 4, dim > shape_4th_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe_system.cc:1033
virtual Point< dim-1 > unit_face_support_point(const unsigned int index) const
Definition: fe_system.cc:2946
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:2915
virtual std::string get_name() const =0
unsigned int n_base_elements() const
Definition: fe.h:2756
unsigned int n_components() const
size_type m() const
virtual std::size_t memory_consumption() const
Definition: fe.cc:1173
Second derivatives of shape functions.
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe.cc:868
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const
Definition: fe.cc:885
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:158
const unsigned int dofs_per_cell
Definition: fe_base.h:283
virtual Tensor< 3, dim > shape_3rd_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe_system.cc:988
void initialize(const std::vector< const FiniteElement< dim, spacedim > * > &fes, const std::vector< unsigned int > &multiplicities)
Definition: fe_system.cc:2285
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe_system.cc:2900
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
virtual FiniteElement< dim, spacedim >::InternalDataBase * get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature,::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe_system.cc:1358
void compute_fill(const Mapping< dim, spacedim > &mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim_1 > &quadrature, const CellSimilarity::Similarity cell_similarity, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const typename FiniteElement< dim, spacedim >::InternalDataBase &fe_data, const internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data, internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe_system.cc:1587
internal::FEValues::FiniteElementRelatedData< dim, spacedim > & get_fe_output_object(const unsigned int base_no) const
Definition: fe_system.cc:98
std::pair< std::pair< unsigned int, unsigned int >, unsigned int > face_system_to_base_index(const unsigned int index) const
Definition: fe.h:2840
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim-1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal,::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
Shape function gradients.
const Conformity conforming_space
Definition: fe_base.h:304
FESystem(const FiniteElement< dim, spacedim > &fe, const unsigned int n_elements)
Definition: fe_system.cc:607
void push_back(const size_type size)
const unsigned int dofs_per_face
Definition: fe_base.h:276
virtual FiniteElement< dim, spacedim >::InternalDataBase * get_face_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim-1 > &quadrature,::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe_system.cc:1414
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe_system.cc:2846
void build_cell_tables()
Definition: fe_system.cc:1762
virtual Tensor< 4, dim > shape_4th_derivative(const unsigned int i, const Point< dim > &p) const
Definition: fe_system.cc:1018
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe_system.cc:2496
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
unsigned int size(const unsigned int i) const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
Definition: fe_system.cc:831
virtual ~FESystem()
Definition: fe_system.cc:776
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe.cc:1236
UpdateFlags update_each
Definition: fe.h:644
Task< RT > new_task(const std_cxx11::function< RT()> &function)