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