Reference documentation for deal.II version GIT a3f17f8a20 2023-06-02 10:50:02+00:00
\(\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_enriched.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 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 
17 #include <deal.II/fe/fe_enriched.h>
18 #include <deal.II/fe/fe_tools.h>
19 
22 
23 #include <memory>
24 
26 
27 namespace internal
28 {
29  namespace FE_Enriched
30  {
31  namespace
32  {
37  template <typename T>
38  std::vector<unsigned int>
39  build_multiplicities(const std::vector<std::vector<T>> &functions)
40  {
41  std::vector<unsigned int> multiplicities;
42  multiplicities.push_back(1); // the first one is non-enriched FE
43  for (unsigned int i = 0; i < functions.size(); ++i)
44  multiplicities.push_back(functions[i].size());
45 
46  return multiplicities;
47  }
48 
49 
53  template <int dim, int spacedim>
54  std::vector<const FiniteElement<dim, spacedim> *>
55  build_fes(
56  const FiniteElement<dim, spacedim> * fe_base,
57  const std::vector<const FiniteElement<dim, spacedim> *> &fe_enriched)
58  {
59  std::vector<const FiniteElement<dim, spacedim> *> fes;
60  fes.push_back(fe_base);
61  for (unsigned int i = 0; i < fe_enriched.size(); ++i)
62  fes.push_back(fe_enriched[i]);
63 
64  return fes;
65  }
66 
67 
72  template <int dim, int spacedim>
73  bool
74  consistency_check(
75  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
76  const std::vector<unsigned int> & multiplicities,
77  const std::vector<std::vector<std::function<const Function<spacedim> *(
78  const typename ::Triangulation<dim, spacedim>::cell_iterator
79  &)>>> & functions)
80  {
81  AssertThrow(fes.size() > 0, ExcMessage("FEs size should be >=1"));
82  AssertThrow(fes.size() == multiplicities.size(),
83  ExcMessage(
84  "FEs and multiplicities should have the same size"));
85 
86  AssertThrow(functions.size() == fes.size() - 1,
87  ExcDimensionMismatch(functions.size(), fes.size() - 1));
88 
89  AssertThrow(multiplicities[0] == 1,
90  ExcMessage("First multiplicity should be 1"));
91 
92  const unsigned int n_comp_base = fes[0]->n_components();
93 
94  // start from fe=1 as 0th is always non-enriched FE.
95  for (unsigned int fe = 1; fe < fes.size(); ++fe)
96  {
97  const FE_Nothing<dim> *fe_nothing =
98  dynamic_cast<const FE_Nothing<dim> *>(fes[fe]);
99  if (fe_nothing)
100  AssertThrow(
101  fe_nothing->is_dominating(),
102  ExcMessage(
103  "Only dominating FE_Nothing can be used in FE_Enriched"));
104 
105  AssertThrow(
106  fes[fe]->n_components() == n_comp_base,
107  ExcMessage(
108  "All elements must have the same number of components"));
109  }
110  return true;
111  }
112 
113 
118  template <int dim, int spacedim>
119  bool
120  check_if_enriched(
121  const std::vector<const FiniteElement<dim, spacedim> *> &fes)
122  {
123  // start from fe=1 as 0th is always non-enriched FE.
124  for (unsigned int fe = 1; fe < fes.size(); ++fe)
125  if (dynamic_cast<const FE_Nothing<dim> *>(fes[fe]) == nullptr)
126  // this is not FE_Nothing => there will be enrichment
127  return true;
128 
129  return false;
130  }
131  } // namespace
132  } // namespace FE_Enriched
133 } // namespace internal
134 
135 
136 template <int dim, int spacedim>
138  const FiniteElement<dim, spacedim> &fe_base)
139  : FE_Enriched<dim, spacedim>(fe_base,
140  FE_Nothing<dim, spacedim>(fe_base.n_components(),
141  true),
142  nullptr)
143 {}
144 
145 
146 template <int dim, int spacedim>
148  const FiniteElement<dim, spacedim> &fe_base,
149  const FiniteElement<dim, spacedim> &fe_enriched,
150  const Function<spacedim> * enrichment_function)
151  : FE_Enriched<dim, spacedim>(
152  &fe_base,
153  std::vector<const FiniteElement<dim, spacedim> *>(1, &fe_enriched),
154  std::vector<std::vector<std::function<const Function<spacedim> *(
155  const typename Triangulation<dim, spacedim>::cell_iterator &)>>>(
156  1,
157  std::vector<std::function<const Function<spacedim> *(
158  const typename Triangulation<dim, spacedim>::cell_iterator &)>>(
159  1,
160  [=](const typename Triangulation<dim, spacedim>::cell_iterator &)
161  -> const Function<spacedim> * { return enrichment_function; })))
162 {}
163 
164 
165 template <int dim, int spacedim>
167  const FiniteElement<dim, spacedim> * fe_base,
168  const std::vector<const FiniteElement<dim, spacedim> *> &fe_enriched,
169  const std::vector<std::vector<std::function<const Function<spacedim> *(
171  : FE_Enriched<dim, spacedim>(
172  internal::FE_Enriched::build_fes(fe_base, fe_enriched),
173  internal::FE_Enriched::build_multiplicities(functions),
174  functions)
175 {}
176 
177 
178 template <int dim, int spacedim>
180  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
181  const std::vector<unsigned int> & multiplicities,
182  const std::vector<std::vector<std::function<const Function<spacedim> *(
184  : FiniteElement<dim, spacedim>(
185  FETools::Compositing::multiply_dof_numbers(fes, multiplicities, false),
187  fes,
188  multiplicities),
189  FETools::Compositing::compute_nonzero_components(fes,
190  multiplicities,
191  false))
192  , enrichments(functions)
193  , is_enriched(internal::FE_Enriched::check_if_enriched(fes))
194  , fe_system(std::make_unique<FESystem<dim, spacedim>>(fes, multiplicities))
195 {
196  // descriptive error are thrown within the function.
197  Assert(internal::FE_Enriched::consistency_check(fes,
198  multiplicities,
199  functions),
200  ExcInternalError());
201 
202  initialize(fes, multiplicities);
203 
204  // resize to be consistent with all FEs used to construct the FE_Enriched,
205  // even though we will never use the 0th element.
206  base_no_mult_local_enriched_dofs.resize(fes.size());
207  for (unsigned int fe = 1; fe < fes.size(); ++fe)
208  base_no_mult_local_enriched_dofs[fe].resize(multiplicities[fe]);
209 
210  Assert(base_no_mult_local_enriched_dofs.size() == this->n_base_elements(),
212  this->n_base_elements()));
213 
214  // build the map: (base_no, base_m) -> vector of local element DoFs
215  for (unsigned int system_index = 0; system_index < this->n_dofs_per_cell();
216  ++system_index)
217  {
218  const unsigned int base_no =
219  this->system_to_base_table[system_index].first.first;
220  if (base_no == 0) // 0th is always non-enriched FE
221  continue;
222 
223  const unsigned int base_m =
224  this->system_to_base_table[system_index].first.second;
225 
226  Assert(base_m < base_no_mult_local_enriched_dofs[base_no].size(),
227  ExcMessage(
228  "Size mismatch for base_no_mult_local_enriched_dofs: "
229  "base_index = " +
230  std::to_string(this->system_to_base_table[system_index].second) +
231  "; base_no = " + std::to_string(base_no) +
232  "; base_m = " + std::to_string(base_m) +
233  "; system_index = " + std::to_string(system_index)));
234 
235  Assert(base_m < base_no_mult_local_enriched_dofs[base_no].size(),
237  base_m, base_no_mult_local_enriched_dofs[base_no].size()));
238 
239  base_no_mult_local_enriched_dofs[base_no][base_m].push_back(system_index);
240  }
241 
242  // make sure that local_enriched_dofs.size() is correct, that is equals to
243  // DoFs per cell of the corresponding FE.
244  for (unsigned int base_no = 1;
245  base_no < base_no_mult_local_enriched_dofs.size();
246  base_no++)
247  {
248  for (unsigned int m = 0;
249  m < base_no_mult_local_enriched_dofs[base_no].size();
250  m++)
251  Assert(base_no_mult_local_enriched_dofs[base_no][m].size() ==
252  fes[base_no]->n_dofs_per_cell(),
254  base_no_mult_local_enriched_dofs[base_no][m].size(),
255  fes[base_no]->n_dofs_per_cell()));
256  }
257 }
258 
259 
260 template <int dim, int spacedim>
261 const std::vector<std::vector<std::function<const Function<spacedim> *(
264 {
265  return enrichments;
266 }
267 
268 
269 template <int dim, int spacedim>
270 double
272  const Point<dim> & p) const
273 {
274  Assert(
275  !is_enriched,
276  ExcMessage(
277  "For enriched finite elements shape_value() can not be defined on the reference element."));
278  return fe_system->shape_value(i, p);
279 }
280 
281 
282 template <int dim, int spacedim>
283 std::unique_ptr<FiniteElement<dim, spacedim>>
285 {
286  std::vector<const FiniteElement<dim, spacedim> *> fes;
287  std::vector<unsigned int> multiplicities;
288 
289  for (unsigned int i = 0; i < this->n_base_elements(); ++i)
290  {
291  fes.push_back(&base_element(i));
292  multiplicities.push_back(this->element_multiplicity(i));
293  }
294 
295  return std::unique_ptr<FE_Enriched<dim, spacedim>>(
296  new FE_Enriched<dim, spacedim>(fes, multiplicities, get_enrichments()));
297 }
298 
299 
300 template <int dim, int spacedim>
303 {
304  UpdateFlags out = fe_system->requires_update_flags(flags);
305 
306  if (is_enriched)
307  {
308  // if we ask for values or gradients, then we would need quadrature points
309  if (flags & (update_values | update_gradients))
311 
312  // if need gradients, add update_values due to product rule
313  if (out & update_gradients)
314  out |= update_values;
315  }
316 
318 
319  return out;
320 }
321 
322 
323 template <int dim, int spacedim>
324 template <int dim_1>
325 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
327  std::unique_ptr<typename FESystem<dim, spacedim>::InternalData> fes_data,
328  const UpdateFlags flags,
329  const Quadrature<dim_1> &quadrature) const
330 {
331  // Pass ownership of the FiniteElement::InternalDataBase object
332  // that fes_data points to, to the new InternalData object.
333  auto update_each_flags = fes_data->update_each;
334  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
335  data_ptr = std::make_unique<InternalData>(std::move(fes_data));
336  auto &data = dynamic_cast<InternalData &>(*data_ptr);
337 
338  // copy update_each from FESystem data:
339  data.update_each = update_each_flags;
340 
341  // resize cache array according to requested flags
342  data.enrichment.resize(this->n_base_elements());
343 
344  const unsigned int n_q_points = quadrature.size();
345 
346  for (unsigned int base = 0; base < this->n_base_elements(); ++base)
347  {
348  data.enrichment[base].resize(this->element_multiplicity(base));
349  for (unsigned int m = 0; m < this->element_multiplicity(base); ++m)
350  {
351  if (flags & update_values)
352  data.enrichment[base][m].values.resize(n_q_points);
353 
354  if (flags & update_gradients)
355  data.enrichment[base][m].gradients.resize(n_q_points);
356 
357  if (flags & update_hessians)
358  data.enrichment[base][m].hessians.resize(n_q_points);
359  }
360  }
361 
362  return data_ptr;
363 }
364 
365 
366 template <int dim, int spacedim>
367 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
369  const UpdateFlags update_flags,
370  const Mapping<dim, spacedim> & mapping,
371  const hp::QCollection<dim - 1> &quadrature,
373  &output_data) const
374 {
375  AssertDimension(quadrature.size(), 1);
376 
377  auto data =
378  fe_system->get_face_data(update_flags, mapping, quadrature, output_data);
379  return setup_data(Utilities::dynamic_unique_cast<
381  std::move(data)),
382  update_flags,
383  quadrature[0]);
384 }
385 
386 
387 template <int dim, int spacedim>
388 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
390  const UpdateFlags update_flags,
391  const Mapping<dim, spacedim> &mapping,
392  const Quadrature<dim - 1> & quadrature,
394  spacedim>
395  &output_data) const
396 {
397  auto data =
398  fe_system->get_subface_data(update_flags, mapping, quadrature, output_data);
399  return setup_data(Utilities::dynamic_unique_cast<
401  std::move(data)),
402  update_flags,
403  quadrature);
404 }
405 
406 
407 template <int dim, int spacedim>
408 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
410  const UpdateFlags flags,
411  const Mapping<dim, spacedim> &mapping,
412  const Quadrature<dim> & quadrature,
414  &output_data) const
415 {
416  auto data = fe_system->get_data(flags, mapping, quadrature, output_data);
417  return setup_data(Utilities::dynamic_unique_cast<
419  std::move(data)),
420  flags,
421  quadrature);
422 }
423 
424 
425 template <int dim, int spacedim>
426 void
428  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
429  const std::vector<unsigned int> & multiplicities)
430 {
431  Assert(fes.size() == multiplicities.size(),
432  ExcDimensionMismatch(fes.size(), multiplicities.size()));
433 
434  // Note that we need to skip every FE with multiplicity 0 in the following
435  // block of code
436  this->base_to_block_indices.reinit(0, 0);
437 
438  for (unsigned int i = 0; i < fes.size(); ++i)
439  if (multiplicities[i] > 0)
440  this->base_to_block_indices.push_back(multiplicities[i]);
441 
442  {
443  // If the system is not primitive, these have not been initialized by
444  // FiniteElement
445  this->system_to_component_table.resize(this->n_dofs_per_cell());
446 
447  FETools::Compositing::build_cell_tables(this->system_to_base_table,
448  this->system_to_component_table,
449  this->component_to_base_table,
450  *this,
451  false);
452 
453  this->face_system_to_component_table.resize(this->n_unique_faces());
454 
455  for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
456  {
457  this->face_system_to_component_table[0].resize(
458  this->n_dofs_per_face(face_no));
459 
460 
462  this->face_system_to_base_table[face_no],
463  this->face_system_to_component_table[face_no],
464  *this,
465  false,
466  face_no);
467  }
468  }
469 
470  // restriction and prolongation matrices are built on demand
471 
472  // now set up the interface constraints for h-refinement.
473  // take them from fe_system:
474  this->interface_constraints = fe_system->interface_constraints;
475 
476  // if we just wrap another FE (i.e. use FE_Nothing as a second FE)
477  // then it makes sense to have support points.
478  // However, functions like interpolate_boundary_values() need all FEs inside
479  // FECollection to be able to provide support points irrespectively whether
480  // this FE sits on the boundary or not. Thus for moment just copy support
481  // points from FE system:
482  {
483  this->unit_support_points = fe_system->unit_support_points;
484  this->unit_face_support_points = fe_system->unit_face_support_points;
485  }
486 
487  // take adjust_quad_dof_index_for_face_orientation_table from FESystem:
488  {
489  this->adjust_line_dof_index_for_line_orientation_table =
490  fe_system->adjust_line_dof_index_for_line_orientation_table;
491  }
492 }
493 
494 
495 template <int dim, int spacedim>
496 std::string
498 {
499  std::ostringstream namebuf;
500 
501  namebuf << "FE_Enriched<" << Utilities::dim_string(dim, spacedim) << ">[";
502  for (unsigned int i = 0; i < this->n_base_elements(); ++i)
503  {
504  namebuf << base_element(i).get_name();
505  if (this->element_multiplicity(i) != 1)
506  namebuf << '^' << this->element_multiplicity(i);
507  if (i != this->n_base_elements() - 1)
508  namebuf << '-';
509  }
510  namebuf << ']';
511 
512  return namebuf.str();
513 }
514 
515 
516 template <int dim, int spacedim>
518 FE_Enriched<dim, spacedim>::base_element(const unsigned int index) const
519 {
520  return fe_system->base_element(index);
521 }
522 
523 
524 template <int dim, int spacedim>
525 void
527  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
528  const CellSimilarity::Similarity cell_similarity,
529  const Quadrature<dim> & quadrature,
530  const Mapping<dim, spacedim> & mapping,
531  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
533  & mapping_data,
534  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
536  &output_data) const
537 {
538  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
539  ExcInternalError());
540  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
541 
542  // call FESystem's method to fill everything without enrichment function
543  fe_system->fill_fe_values(cell,
544  cell_similarity,
545  quadrature,
546  mapping,
547  mapping_internal,
548  mapping_data,
549  *fe_data.fesystem_data,
550  output_data);
551 
552  if (is_enriched)
553  multiply_by_enrichment(
554  quadrature, fe_data, mapping_data, cell, output_data);
555 }
556 
557 
558 template <int dim, int spacedim>
559 void
561  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
562  const unsigned int face_no,
563  const hp::QCollection<dim - 1> & quadrature,
564  const Mapping<dim, spacedim> & mapping,
565  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
567  & mapping_data,
568  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
570  &output_data) const
571 {
572  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
573  ExcInternalError());
574  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
575 
576  AssertDimension(quadrature.size(), 1);
577 
578  // call FESystem's method to fill everything without enrichment function
579  fe_system->fill_fe_face_values(cell,
580  face_no,
581  quadrature,
582  mapping,
583  mapping_internal,
584  mapping_data,
585  *fe_data.fesystem_data,
586  output_data);
587 
588  if (is_enriched)
589  multiply_by_enrichment(
590  quadrature[0], fe_data, mapping_data, cell, output_data);
591 }
592 
593 
594 template <int dim, int spacedim>
595 void
597  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
598  const unsigned int face_no,
599  const unsigned int sub_no,
600  const Quadrature<dim - 1> & quadrature,
601  const Mapping<dim, spacedim> & mapping,
602  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
604  & mapping_data,
605  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
607  &output_data) const
608 {
609  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
610  ExcInternalError());
611  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
612 
613  // call FESystem's method to fill everything without enrichment function
614  fe_system->fill_fe_subface_values(cell,
615  face_no,
616  sub_no,
617  quadrature,
618  mapping,
619  mapping_internal,
620  mapping_data,
621  *fe_data.fesystem_data,
622  output_data);
623 
624  if (is_enriched)
625  multiply_by_enrichment(
626  quadrature, fe_data, mapping_data, cell, output_data);
627 }
628 
629 
630 template <int dim, int spacedim>
631 template <int dim_1>
632 void
634  const Quadrature<dim_1> &quadrature,
635  const InternalData & fe_data,
637  & mapping_data,
638  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
640  &output_data) const
641 {
642  // mapping_data will contain quadrature points on the real element.
643  // fe_internal is needed to get update flags
644  // finally, output_data should store all the results we need.
645 
646  // Either dim_1==dim
647  // (fill_fe_values) or dim_1==dim-1
648  // (fill_fe_(sub)face_values)
649  Assert(dim_1 == dim || dim_1 == dim - 1, ExcInternalError());
650  const UpdateFlags flags = fe_data.update_each;
651 
652  const unsigned int n_q_points = quadrature.size();
653 
654  // First, populate output_data object (that shall hold everything requested
655  // such as shape value, gradients, hessians, etc) from each base element. That
656  // is almost identical to FESystem::compute_fill_one_base(), the difference
657  // being that we do it irrespectively of cell_similarity and use
658  // base_fe_data.update_flags
659 
660  // TODO: do we need it only for dim_1 == dim (i.e. fill_fe_values)?
661  if (dim_1 == dim)
662  for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
663  {
664  const FiniteElement<dim, spacedim> &base_fe = base_element(base_no);
665  typename FiniteElement<dim, spacedim>::InternalDataBase &base_fe_data =
666  fe_data.get_fe_data(base_no);
668  spacedim>
669  &base_data = fe_data.get_fe_output_object(base_no);
670 
671  const UpdateFlags base_flags = base_fe_data.update_each;
672 
673  for (unsigned int system_index = 0;
674  system_index < this->n_dofs_per_cell();
675  ++system_index)
676  if (this->system_to_base_table[system_index].first.first == base_no)
677  {
678  const unsigned int base_index =
679  this->system_to_base_table[system_index].second;
680  Assert(base_index < base_fe.n_dofs_per_cell(),
681  ExcInternalError());
682 
683  // now copy. if the shape function is primitive, then there
684  // is only one value to be copied, but for non-primitive
685  // elements, there might be more values to be copied
686  //
687  // so, find out from which index to take this one value, and
688  // to which index to put
689  unsigned int out_index = 0;
690  for (unsigned int i = 0; i < system_index; ++i)
691  out_index += this->n_nonzero_components(i);
692  unsigned int in_index = 0;
693  for (unsigned int i = 0; i < base_index; ++i)
694  in_index += base_fe.n_nonzero_components(i);
695 
696  // then loop over the number of components to be copied
697  Assert(this->n_nonzero_components(system_index) ==
698  base_fe.n_nonzero_components(base_index),
699  ExcInternalError());
700  for (unsigned int s = 0;
701  s < this->n_nonzero_components(system_index);
702  ++s)
703  {
704  if (base_flags & update_values)
705  for (unsigned int q = 0; q < n_q_points; ++q)
706  output_data.shape_values[out_index + s][q] =
707  base_data.shape_values(in_index + s, q);
708 
709  if (base_flags & update_gradients)
710  for (unsigned int q = 0; q < n_q_points; ++q)
711  output_data.shape_gradients[out_index + s][q] =
712  base_data.shape_gradients[in_index + s][q];
713 
714  if (base_flags & update_hessians)
715  for (unsigned int q = 0; q < n_q_points; ++q)
716  output_data.shape_hessians[out_index + s][q] =
717  base_data.shape_hessians[in_index + s][q];
718  }
719  }
720  }
721 
722  Assert(base_no_mult_local_enriched_dofs.size() == fe_data.enrichment.size(),
723  ExcDimensionMismatch(base_no_mult_local_enriched_dofs.size(),
724  fe_data.enrichment.size()));
725  // calculate hessians, gradients and values for each function
726  for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
727  {
728  Assert(
729  base_no_mult_local_enriched_dofs[base_no].size() ==
730  fe_data.enrichment[base_no].size(),
731  ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no].size(),
732  fe_data.enrichment[base_no].size()));
733  for (unsigned int m = 0;
734  m < base_no_mult_local_enriched_dofs[base_no].size();
735  m++)
736  {
737  // Avoid evaluating quadrature points if no dofs are assigned. This
738  // happens when FE_Nothing is used together with other FE (i.e. FE_Q)
739  // as enrichments.
740  if (base_no_mult_local_enriched_dofs[base_no][m].size() == 0)
741  continue;
742 
743  Assert(enrichments[base_no - 1][m](cell) != nullptr,
744  ExcMessage(
745  "The pointer to the enrichment function is not set"));
746 
747  Assert(enrichments[base_no - 1][m](cell)->n_components == 1,
748  ExcMessage(
749  "Only scalar-valued enrichment functions are allowed"));
750 
751  if (flags & update_hessians)
752  {
753  Assert(fe_data.enrichment[base_no][m].hessians.size() ==
754  n_q_points,
756  fe_data.enrichment[base_no][m].hessians.size(),
757  n_q_points));
758  for (unsigned int q = 0; q < n_q_points; ++q)
759  fe_data.enrichment[base_no][m].hessians[q] =
760  enrichments[base_no - 1][m](cell)->hessian(
761  mapping_data.quadrature_points[q]);
762  }
763 
764  if (flags & update_gradients)
765  {
766  Assert(fe_data.enrichment[base_no][m].gradients.size() ==
767  n_q_points,
769  fe_data.enrichment[base_no][m].gradients.size(),
770  n_q_points));
771  for (unsigned int q = 0; q < n_q_points; ++q)
772  fe_data.enrichment[base_no][m].gradients[q] =
773  enrichments[base_no - 1][m](cell)->gradient(
774  mapping_data.quadrature_points[q]);
775  }
776 
777  if (flags & update_values)
778  {
779  Assert(fe_data.enrichment[base_no][m].values.size() == n_q_points,
781  fe_data.enrichment[base_no][m].values.size(),
782  n_q_points));
783  for (unsigned int q = 0; q < n_q_points; ++q)
784  fe_data.enrichment[base_no][m].values[q] =
785  enrichments[base_no - 1][m](cell)->value(
786  mapping_data.quadrature_points[q]);
787  }
788  }
789  }
790 
791  // Finally, update the standard data stored in output_data
792  // by expanding the product rule for enrichment function.
793  // note that the order if important, namely
794  // output_data.shape_XYZ contains values of standard FEM and we overwrite
795  // it with the updated one in the following order: hessians -> gradients ->
796  // values
797  if (flags & update_hessians)
798  {
799  for (unsigned int base_no = 1; base_no < this->n_base_elements();
800  base_no++)
801  {
802  for (unsigned int m = 0;
803  m < base_no_mult_local_enriched_dofs[base_no].size();
804  m++)
805  for (unsigned int i = 0;
806  i < base_no_mult_local_enriched_dofs[base_no][m].size();
807  i++)
808  {
809  const unsigned int enriched_dof =
810  base_no_mult_local_enriched_dofs[base_no][m][i];
811  for (unsigned int q = 0; q < n_q_points; ++q)
812  {
813  const Tensor<2, spacedim> grad_grad = outer_product(
814  output_data.shape_gradients[enriched_dof][q],
815  fe_data.enrichment[base_no][m].gradients[q]);
816  const Tensor<2, spacedim, double> sym_grad_grad =
817  symmetrize(grad_grad) * 2.0; // symmetrize does [s+s^T]/2
818 
819  output_data.shape_hessians[enriched_dof][q] *=
820  fe_data.enrichment[base_no][m].values[q];
821  output_data.shape_hessians[enriched_dof][q] +=
822  sym_grad_grad +
823  output_data.shape_values[enriched_dof][q] *
824  fe_data.enrichment[base_no][m].hessians[q];
825  }
826  }
827  }
828  }
829 
830  if (flags & update_gradients)
831  for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
832  {
833  for (unsigned int m = 0;
834  m < base_no_mult_local_enriched_dofs[base_no].size();
835  m++)
836  for (unsigned int i = 0;
837  i < base_no_mult_local_enriched_dofs[base_no][m].size();
838  i++)
839  {
840  const unsigned int enriched_dof =
841  base_no_mult_local_enriched_dofs[base_no][m][i];
842  for (unsigned int q = 0; q < n_q_points; ++q)
843  {
844  output_data.shape_gradients[enriched_dof][q] *=
845  fe_data.enrichment[base_no][m].values[q];
846  output_data.shape_gradients[enriched_dof][q] +=
847  output_data.shape_values[enriched_dof][q] *
848  fe_data.enrichment[base_no][m].gradients[q];
849  }
850  }
851  }
852 
853  if (flags & update_values)
854  for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
855  {
856  for (unsigned int m = 0;
857  m < base_no_mult_local_enriched_dofs[base_no].size();
858  m++)
859  for (unsigned int i = 0;
860  i < base_no_mult_local_enriched_dofs[base_no][m].size();
861  i++)
862  {
863  const unsigned int enriched_dof =
864  base_no_mult_local_enriched_dofs[base_no][m][i];
865  for (unsigned int q = 0; q < n_q_points; ++q)
866  {
867  output_data.shape_values[enriched_dof][q] *=
868  fe_data.enrichment[base_no][m].values[q];
869  }
870  }
871  }
872 }
873 
874 
875 template <int dim, int spacedim>
878 {
879  return *fe_system;
880 }
881 
882 
883 template <int dim, int spacedim>
884 bool
886 {
887  return true;
888 }
889 
890 
891 template <int dim, int spacedim>
892 void
894  const FiniteElement<dim, spacedim> &source,
896  const unsigned int face_no) const
897 {
898  if (const FE_Enriched<dim, spacedim> *fe_enr_other =
899  dynamic_cast<const FE_Enriched<dim, spacedim> *>(&source))
900  {
901  fe_system->get_face_interpolation_matrix(fe_enr_other->get_fe_system(),
902  matrix,
903  face_no);
904  }
905  else
906  {
907  AssertThrow(
908  false,
909  (typename FiniteElement<dim,
910  spacedim>::ExcInterpolationNotImplemented()));
911  }
912 }
913 
914 
915 template <int dim, int spacedim>
916 void
918  const FiniteElement<dim, spacedim> &source,
919  const unsigned int subface,
921  const unsigned int face_no) const
922 {
923  if (const FE_Enriched<dim, spacedim> *fe_enr_other =
924  dynamic_cast<const FE_Enriched<dim, spacedim> *>(&source))
925  {
926  fe_system->get_subface_interpolation_matrix(fe_enr_other->get_fe_system(),
927  subface,
928  matrix,
929  face_no);
930  }
931  else
932  {
933  AssertThrow(
934  false,
935  (typename FiniteElement<dim,
936  spacedim>::ExcInterpolationNotImplemented()));
937  }
938 }
939 
940 
941 template <int dim, int spacedim>
942 std::vector<std::pair<unsigned int, unsigned int>>
944  const FiniteElement<dim, spacedim> &fe_other) const
945 {
946  if (const FE_Enriched<dim, spacedim> *fe_enr_other =
947  dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
948  {
949  return fe_system->hp_vertex_dof_identities(fe_enr_other->get_fe_system());
950  }
951  else
952  {
953  Assert(false, ExcNotImplemented());
954  return std::vector<std::pair<unsigned int, unsigned int>>();
955  }
956 }
957 
958 
959 template <int dim, int spacedim>
960 std::vector<std::pair<unsigned int, unsigned int>>
962  const FiniteElement<dim, spacedim> &fe_other) const
963 {
964  if (const FE_Enriched<dim, spacedim> *fe_enr_other =
965  dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
966  {
967  return fe_system->hp_line_dof_identities(fe_enr_other->get_fe_system());
968  }
969  else
970  {
971  Assert(false, ExcNotImplemented());
972  return std::vector<std::pair<unsigned int, unsigned int>>();
973  }
974 }
975 
976 
977 template <int dim, int spacedim>
978 std::vector<std::pair<unsigned int, unsigned int>>
980  const FiniteElement<dim, spacedim> &fe_other,
981  const unsigned int face_no) const
982 {
983  if (const FE_Enriched<dim, spacedim> *fe_enr_other =
984  dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
985  {
986  return fe_system->hp_quad_dof_identities(fe_enr_other->get_fe_system(),
987  face_no);
988  }
989  else
990  {
991  Assert(false, ExcNotImplemented());
992  return std::vector<std::pair<unsigned int, unsigned int>>();
993  }
994 }
995 
996 
997 template <int dim, int spacedim>
1000  const FiniteElement<dim, spacedim> &fe_other,
1001  const unsigned int codim) const
1002 {
1003  Assert(codim <= dim, ExcImpossibleInDim(dim));
1004 
1005  // vertex/line/face/cell domination
1006  // --------------------------------
1007  // need to decide which element constrain another.
1008  // for example Q(2) dominate Q(4) and thus some DoFs of Q(4) will be
1009  // constrained. If we have Q(2) and Q(4)+POU, then it's clear that Q(2)
1010  // dominates, namely our DoFs will be constrained to make field continuous.
1011  // However, we need to check for situations like Q(4) vs Q(2)+POU.
1012  // In that case the domination for the underlying FEs should be the other way,
1013  // but this implies that we can't constrain POU dofs to make the field
1014  // continuous. In that case, throw an error
1015 
1016  // if it's also enriched, do domination based on each one's FESystem
1017  if (const FE_Enriched<dim, spacedim> *fe_enr_other =
1018  dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
1019  {
1020  return fe_system->compare_for_domination(fe_enr_other->get_fe_system(),
1021  codim);
1022  }
1023  else
1024  {
1025  Assert(false, ExcNotImplemented());
1027  }
1028 }
1029 
1030 
1031 template <int dim, int spacedim>
1032 const FullMatrix<double> &
1034  const unsigned int child,
1035  const RefinementCase<dim> &refinement_case) const
1036 {
1037  return fe_system->get_prolongation_matrix(child, refinement_case);
1038 }
1039 
1040 
1041 template <int dim, int spacedim>
1042 const FullMatrix<double> &
1044  const unsigned int child,
1045  const RefinementCase<dim> &refinement_case) const
1046 {
1047  return fe_system->get_restriction_matrix(child, refinement_case);
1048 }
1049 
1050 
1051 /* ----------------------- FESystem::InternalData ------------------- */
1052 
1053 
1054 template <int dim, int spacedim>
1056  std::unique_ptr<typename FESystem<dim, spacedim>::InternalData> fesystem_data)
1057  : fesystem_data(std::move(fesystem_data))
1058 {}
1059 
1060 
1061 template <int dim, int spacedim>
1064  const unsigned int base_no) const
1065 {
1066  return fesystem_data->get_fe_data(base_no);
1067 }
1068 
1069 
1070 template <int dim, int spacedim>
1073  const unsigned int base_no) const
1074 {
1075  return fesystem_data->get_fe_output_object(base_no);
1076 }
1077 
1078 
1079 namespace ColorEnriched
1080 {
1081  namespace internal
1082  {
1083  template <int dim, int spacedim>
1084  bool
1086  const DoFHandler<dim, spacedim> & dof_handler,
1087  const predicate_function<dim, spacedim> &predicate_1,
1088  const predicate_function<dim, spacedim> &predicate_2)
1089  {
1090  // Use a vector to mark vertices
1091  std::vector<bool> vertices_subdomain_1(
1092  dof_handler.get_triangulation().n_vertices(), false);
1093 
1094  // Mark vertices that belong to cells in subdomain 1
1095  for (const auto &cell : dof_handler.active_cell_iterators())
1096  if (predicate_1(cell)) // True ==> part of subdomain 1
1097  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1098  vertices_subdomain_1[cell->vertex_index(v)] = true;
1099 
1100  // Find if cells in subdomain 2 and subdomain 1 share vertices.
1101  for (const auto &cell : dof_handler.active_cell_iterators())
1102  if (predicate_2(cell)) // True ==> part of subdomain 2
1103  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1104  if (vertices_subdomain_1[cell->vertex_index(v)] == true)
1105  {
1106  return true;
1107  }
1108  return false;
1109  }
1110 
1111 
1112 
1113  template <int dim, int spacedim>
1114  unsigned int
1116  const DoFHandler<dim, spacedim> & mesh,
1117  const std::vector<predicate_function<dim, spacedim>> &predicates,
1118  std::vector<unsigned int> & predicate_colors)
1119  {
1120  const unsigned int num_indices = predicates.size();
1121 
1122  // Use sparsity pattern to represent connections between subdomains.
1123  // Each predicate (i.e a subdomain) is a node in the graph.
1125  dsp.reinit(num_indices, num_indices);
1126 
1127  /*
1128  * Find connections between subdomains taken two at a time.
1129  * If the connection exists, add it to a graph object represented
1130  * by dynamic sparsity pattern.
1131  */
1132  for (unsigned int i = 0; i < num_indices; ++i)
1133  for (unsigned int j = i + 1; j < num_indices; ++j)
1135  predicates[i],
1136  predicates[j]))
1137  dsp.add(i, j);
1138 
1139  dsp.symmetrize();
1140 
1141  // Copy dynamic sparsity pattern to sparsity pattern needed by
1142  // coloring function
1143  SparsityPattern sp_graph;
1144  sp_graph.copy_from(dsp);
1145  predicate_colors.resize(num_indices);
1146 
1147  // Assign each predicate with a color and return number of colors
1148  return SparsityTools::color_sparsity_pattern(sp_graph, predicate_colors);
1149  }
1150 
1151 
1152 
1153  template <int dim, int spacedim>
1154  void
1156  DoFHandler<dim, spacedim> & dof_handler,
1157  const std::vector<predicate_function<dim, spacedim>> &predicates,
1158  const std::vector<unsigned int> & predicate_colors,
1159  std::map<unsigned int, std::map<unsigned int, unsigned int>>
1160  & cellwise_color_predicate_map,
1161  std::vector<std::set<unsigned int>> &fe_sets)
1162  {
1163  // clear output variables first
1164  fe_sets.clear();
1165  cellwise_color_predicate_map.clear();
1166 
1167  /*
1168  * Add first element of fe_sets which is empty by default. This means that
1169  * the default, FE index = 0 is associated with an empty set, since no
1170  * predicate is active in these regions.
1171  */
1172  fe_sets.resize(1);
1173 
1174  /*
1175  * Loop through cells and find set of predicate colors associated
1176  * with the cell. As an example, a cell with an FE index associated
1177  * with colors {a,b} means that predicates active in the cell have
1178  * colors a or b.
1179  *
1180  * Create new active FE index in case of the color
1181  * set is not already listed in fe_sets. If the set already exists,
1182  * find index of the set in fe_sets. In either case, use the id in
1183  * fe_sets to modify cell->active_fe_index.
1184  *
1185  * Associate each cell_id with a set of pairs. The pair represents
1186  * predicate color and the active predicate with that color.
1187  * Each color can only correspond to a single predicate since
1188  * predicates with the same color correspond to disjoint domains.
1189  * This is what the graph coloring in color_predicates
1190  * function ensures. The number of pairs is equal to the number
1191  * of predicates active in the given cell.
1192  */
1193  unsigned int map_index = 0;
1194  for (const auto &cell : dof_handler.active_cell_iterators())
1195  {
1196  // set default FE index ==> no enrichment and no active predicates
1197  cell->set_active_fe_index(0);
1198 
1199  // Give each cell a unique id, which the cellwise_color_predicate_map
1200  // can later use to associate a colors of active predicates with
1201  // the actual predicate id.
1202  //
1203  // When the grid is refined, material id is inherited to
1204  // children cells. So, the cellwise_color_predicate_map stays
1205  // relevant.
1206  cell->set_material_id(map_index);
1207  std::set<unsigned int> color_list;
1208 
1209  // loop through active predicates for the cell and insert map.
1210  // Eg: if the cell with material id 100 has active
1211  // predicates 4 (color = 1) and 5 (color = 2), the map will insert
1212  // pairs (1, 4) and (2, 5) at key 100 (i.e unique id of cell is
1213  // mapped with a map which associates color with predicate id)
1214  // Note that color list for the cell would be {1,2}.
1215  std::map<unsigned int, unsigned int> &cell_map =
1216  cellwise_color_predicate_map[map_index];
1217  for (unsigned int i = 0; i < predicates.size(); ++i)
1218  {
1219  if (predicates[i](cell))
1220  {
1221  /*
1222  * create a pair predicate color and predicate id and add it
1223  * to a map associated with each enriched cell
1224  */
1225  auto ret = cell_map.insert(
1226  std::pair<unsigned int, unsigned int>(predicate_colors[i],
1227  i));
1228 
1229  AssertThrow(ret.second == 1,
1230  ExcMessage(
1231  "Only one enrichment function per color"));
1232 
1233  color_list.insert(predicate_colors[i]);
1234  }
1235  }
1236 
1237 
1238  /*
1239  * check if color combination is already added.
1240  * If already added, set the active FE index based on
1241  * its index in the fe_sets. If the combination doesn't
1242  * exist, add the set to fe_sets and once again set the
1243  * active FE index as last index in fe_sets.
1244  *
1245  * Eg: if cell has color list {1,2} associated and
1246  * fe_sets = { {}, {1}, {2} } for now, we need to add a new set {1,2}
1247  * to fe_sets and a new active FE index 3 because 0 to 2 FE indices
1248  * are already taken by existing sets in fe_sets.
1249  */
1250  if (!color_list.empty())
1251  {
1252  const auto it =
1253  std::find(fe_sets.begin(), fe_sets.end(), color_list);
1254  // when entry is not found
1255  if (it == fe_sets.end())
1256  {
1257  fe_sets.push_back(color_list);
1258  cell->set_active_fe_index(fe_sets.size() - 1);
1259  }
1260  // when entry is found
1261  else
1262  {
1263  cell->set_active_fe_index(std::distance(fe_sets.begin(), it));
1264  }
1265  }
1266  /*
1267  * map_index is used to identify cells and should be unique. The
1268  * index is stored in the material_id of the cell and hence
1269  * stays relevant even when the cells are refined (which is
1270  * why cell_id is not used).
1271  * Two distant cells can have the same set of colors but different
1272  * enrichment functions can be associated with any given
1273  * color. So, in order to figure which enrichment function
1274  * belongs to a color, we use a map that uses this index.
1275  */
1276  ++map_index;
1277  }
1278 
1279 
1280  /*
1281  * Treat interface between enriched cells specially,
1282  * until #1496 (https://github.com/dealii/dealii/issues/1496) is resolved.
1283  * Each time we build constraints at the interface between two different
1284  * FE_Enriched, we look for the least dominating FE of their common
1285  * subspace via hp::FECollection::find_dominating_fe_extended().
1286  * If we don't take further actions, we may find a dominating FE that is
1287  * too restrictive, i.e. enriched FE consisting of only FE_Nothing. New
1288  * elements needs to be added to FECollection object to help find the
1289  * correct enriched FE underlying the spaces in the adjacent cells. This
1290  * is done by creating an appropriate set in fe_sets and a call to the
1291  * function make_fe_collection_from_colored_enrichments at a later stage.
1292  *
1293  * Consider a domain with three predicates and hence with three different
1294  * enrichment functions. Let the enriched finite element of a cell with
1295  * enrichment functions 1 and 2 be represented by [1 1 0], with the last
1296  * entry as zero since the 3rd enrichment function is not relevant for
1297  * the cell. If the interface has enriched FE [1 0 1] and [0 1 1]
1298  * on adjacent cells, an enriched FE [0 0 1] should exist and is
1299  * found as the least dominating finite element for the two cells by
1300  * DoFTools::make_hanging_node_constraints, using the above mentioned
1301  * hp::FECollection functions. Denoting the FE set in adjacent cells as
1302  * {1,3} and {2,3}, this implies that an FE set {3} needs to be added!
1303  * Based on the predicate configuration, this may not be automatically
1304  * done without the following special treatment.
1305  */
1306 
1307  // loop through faces
1308  for (const auto &cell : dof_handler.active_cell_iterators())
1309  {
1310  const unsigned int fe_index = cell->active_fe_index();
1311  const std::set<unsigned int> fe_set = fe_sets.at(fe_index);
1312  for (const unsigned int face : GeometryInfo<dim>::face_indices())
1313  {
1314  // cell shouldn't be at the boundary and
1315  // neighboring cell is not already visited (to avoid visiting
1316  // same face twice). Note that the cells' material ids are
1317  // labeled according to their order in dof_handler previously.
1318  if (!cell->at_boundary(face) &&
1319  cell->material_id() < cell->neighbor(face)->material_id())
1320  {
1321  const auto nbr_fe_index =
1322  cell->neighbor(face)->active_fe_index();
1323 
1324  // find corresponding FE set
1325  const auto nbr_fe_set = fe_sets.at(nbr_fe_index);
1326 
1327  // find intersection of the FE sets: fe_set and nbr_fe_set
1328  std::set<unsigned int> intersection_set;
1329  std::set_intersection(
1330  fe_set.begin(),
1331  fe_set.end(),
1332  nbr_fe_set.begin(),
1333  nbr_fe_set.end(),
1334  std::inserter(intersection_set, intersection_set.begin()));
1335 
1336  // add only the new sets
1337  if (!intersection_set.empty())
1338  {
1339  const auto it = std::find(fe_sets.begin(),
1340  fe_sets.end(),
1341  intersection_set);
1342  // add the set if it is not found
1343  if (it == fe_sets.end())
1344  {
1345  fe_sets.push_back(intersection_set);
1346  }
1347  }
1348  }
1349  }
1350  }
1351  }
1352 
1353 
1354 
1355  template <int dim, int spacedim>
1356  void
1358  const unsigned int n_colors,
1359  const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments,
1360  const std::map<unsigned int, std::map<unsigned int, unsigned int>>
1361  &cellwise_color_predicate_map,
1362  std::vector<std::function<const Function<spacedim> *(
1364  &color_enrichments)
1365  {
1366  color_enrichments.clear();
1367 
1368  // Each color should be associated with a single enrichment function
1369  // called color enrichment function which calls the correct enrichment
1370  // function for a given cell.
1371  //
1372  // Assume that a cell has a active predicates with ids 4 (color = 1) and
1373  // 5 (color = 2). cellwise_color_predicate_map has this information,
1374  // provided we know the material id.
1375  //
1376  // The constructed color_enrichments is such that
1377  // color_enrichments[1](cell) will return return a pointer to
1378  // function with id=4, i.e. enrichments[4].
1379  // In other words, using the previously collected information in
1380  // this function we translate a vector of user provided enrichment
1381  // functions into a vector of functions suitable for FE_Enriched class.
1382  color_enrichments.resize(n_colors);
1383  for (unsigned int i = 0; i < n_colors; ++i)
1384  {
1385  color_enrichments[i] =
1386  [&, i](const typename Triangulation<dim, spacedim>::cell_iterator
1387  &cell) {
1388  const unsigned int id = cell->material_id();
1389 
1390  /*
1391  * i'th color_enrichment function corresponds to (i+1)'th color.
1392  * Since FE_Enriched takes function pointers, we return a
1393  * function pointer.
1394  */
1395  return enrichments[cellwise_color_predicate_map.at(id).at(i + 1)]
1396  .get();
1397  };
1398  }
1399  }
1400 
1401 
1402 
1403  template <int dim, int spacedim>
1404  void
1406  const unsigned int n_colors,
1407  const std::vector<std::set<unsigned int>> &fe_sets,
1408  const std::vector<std::function<const Function<spacedim> *(
1410  & color_enrichments,
1411  const FiniteElement<dim, spacedim> &fe_base,
1412  const FiniteElement<dim, spacedim> &fe_enriched,
1413  const FE_Nothing<dim, spacedim> & fe_nothing,
1414  hp::FECollection<dim, spacedim> & fe_collection)
1415  {
1416  // define dummy function which is associated with FE_Nothing
1417  const std::function<const Function<spacedim> *(
1419  dummy_function =
1420  [=](const typename Triangulation<dim, spacedim>::cell_iterator &)
1421  -> const Function<spacedim> * {
1422  AssertThrow(false,
1423  ExcMessage("Called enrichment function for FE_Nothing"));
1424  return nullptr;
1425  };
1426 
1427 
1428  // loop through color sets and create FE_enriched element for each
1429  // of them provided before calling this function, we have color
1430  // enrichment function associated with each color.
1431  for (const auto &fe_set : fe_sets)
1432  {
1433  std::vector<const FiniteElement<dim, spacedim> *> vec_fe_enriched(
1434  n_colors, &fe_nothing);
1435  std::vector<std::vector<std::function<const Function<spacedim> *(
1436  const typename Triangulation<dim, spacedim>::cell_iterator &)>>>
1437  functions(n_colors, {dummy_function});
1438 
1439  for (const unsigned int color_id : fe_set)
1440  {
1441  // Given a color id, corresponding color enrichment
1442  // function is at index id-1 because color_enrichments are
1443  // indexed from zero and colors are indexed from 1.
1444  const unsigned int ind = color_id - 1;
1445 
1446  AssertIndexRange(ind, vec_fe_enriched.size());
1447  AssertIndexRange(ind, functions.size());
1448  AssertIndexRange(ind, color_enrichments.size());
1449 
1450  // Assume an active predicate colors {1,2} for a cell.
1451  // We then need to create a vector of FE enriched elements
1452  // with vec_fe_enriched[0] = vec_fe_enriched[1] = &fe_enriched
1453  // which can later be associated with enrichment functions.
1454  vec_fe_enriched[ind] = &fe_enriched;
1455 
1456  // color_set_id'th color function is (color_set_id-1)
1457  // element of color wise enrichments
1458  functions[ind][0] = color_enrichments[ind];
1459  }
1460 
1461  AssertDimension(vec_fe_enriched.size(), functions.size());
1462 
1463  FE_Enriched<dim, spacedim> fe_component(&fe_base,
1464  vec_fe_enriched,
1465  functions);
1466  fe_collection.push_back(fe_component);
1467  }
1468  }
1469  } // namespace internal
1470 
1471 
1472 
1473  template <int dim, int spacedim>
1475  const FiniteElement<dim, spacedim> & fe_base,
1476  const FiniteElement<dim, spacedim> & fe_enriched,
1477  const std::vector<predicate_function<dim, spacedim>> & predicates,
1478  const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments)
1479  : fe_base(fe_base)
1480  , fe_enriched(fe_enriched)
1481  , fe_nothing(fe_base.n_components(), true)
1482  , predicates(predicates)
1484  , n_colors(numbers::invalid_unsigned_int)
1485  {
1486  AssertDimension(predicates.size(), enrichments.size());
1488  AssertThrow(predicates.size() > 0,
1489  ExcMessage("Number of predicates should be positive"));
1490  }
1491 
1492 
1493 
1494  template <int dim, int spacedim>
1497  DoFHandler<dim, spacedim> &dof_handler)
1498  {
1499  // color the predicates based on connections between corresponding
1500  // subdomains
1501  n_colors =
1502  internal::color_predicates(dof_handler, predicates, predicate_colors);
1503 
1504  // create color maps and color list for each cell
1506  predicates,
1507  predicate_colors,
1508  cellwise_color_predicate_map,
1509  fe_sets);
1510  // setup color wise enrichment functions
1511  // i'th function corresponds to (i+1) color!
1512  internal::make_colorwise_enrichment_functions<dim, spacedim>(
1513  n_colors, enrichments, cellwise_color_predicate_map, color_enrichments);
1514 
1515  // make FE_Collection
1517  fe_sets,
1518  color_enrichments,
1519  fe_base,
1520  fe_enriched,
1521  fe_nothing,
1522  fe_collection);
1523 
1524  return fe_collection;
1525  }
1526 } // namespace ColorEnriched
1527 
1528 
1529 // explicit instantiations
1530 #include "fe_enriched.inst"
1531 
const Triangulation< dim, spacedim > & get_triangulation() const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
void add(const size_type i, const size_type j)
std::vector< std::vector< EnrichmentValues > > enrichment
Definition: fe_enriched.h:512
internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > & get_fe_output_object(const unsigned int base_no) const
InternalData(std::unique_ptr< typename FESystem< dim, spacedim >::InternalData > fesystem_data)
FiniteElement< dim, spacedim >::InternalDataBase & get_fe_data(const unsigned int base_no) const
std::unique_ptr< typename FESystem< dim, spacedim >::InternalData > fesystem_data
Definition: fe_enriched.h:495
virtual std::string get_name() const override
Definition: fe_enriched.cc:497
std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > setup_data(std::unique_ptr< typename FESystem< dim, spacedim >::InternalData > fes_data, const UpdateFlags flags, const Quadrature< dim_1 > &quadrature) const
Definition: fe_enriched.cc:326
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_enriched.cc:979
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_enriched.cc:596
std::vector< std::vector< std::vector< unsigned int > > > base_no_mult_local_enriched_dofs
Definition: fe_enriched.h:522
void multiply_by_enrichment(const Quadrature< dim_1 > &quadrature, const InternalData &fe_data, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename Triangulation< dim, spacedim >::cell_iterator &cell, internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe_enriched.cc:633
virtual std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
Definition: fe_enriched.cc:409
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_enriched.cc:368
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_enriched.cc:560
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_enriched.cc:999
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
Definition: fe_enriched.cc:271
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_enriched.cc:284
virtual bool hp_constraints_are_implemented() const override
Definition: fe_enriched.cc:885
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_enriched.cc:389
void initialize(const std::vector< const FiniteElement< dim, spacedim > * > &fes, const std::vector< unsigned int > &multiplicities)
Definition: fe_enriched.cc:427
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_enriched.cc:943
const FESystem< dim, spacedim > & get_fe_system() const
Definition: fe_enriched.cc:877
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const override
Definition: fe_enriched.cc:518
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_enriched.cc:961
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
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_enriched.cc:526
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition: fe_enriched.cc:893
FE_Enriched(const FiniteElement< dim, spacedim > &fe_base, const FiniteElement< dim, spacedim > &fe_enriched, const Function< spacedim > *enrichment_function)
Definition: fe_enriched.cc:147
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: fe_enriched.cc:302
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_enriched.cc:917
const std::vector< std::vector< std::function< const Function< spacedim > *(const typename Triangulation< dim, spacedim >::cell_iterator &)> > > enrichments
Definition: fe_enriched.h:533
bool is_dominating() const
Definition: fe_nothing.cc:197
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > system_to_base_table
Definition: fe.h:2527
unsigned int n_nonzero_components(const unsigned int i) const
Definition: point.h:112
unsigned int size() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
unsigned int n_vertices() const
unsigned int size() const
Definition: collection.h:264
void push_back(const FiniteElement< dim, spacedim > &new_fe)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 2 > second
Definition: grid_out.cc:4616
Point< 2 > first
Definition: grid_out.cc:4615
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void make_fe_collection_from_colored_enrichments(const unsigned int n_colors, const std::vector< std::set< unsigned int >> &fe_sets, const std::vector< std::function< const Function< spacedim > *(const typename Triangulation< dim, spacedim >::cell_iterator &)>> &color_enrichments, const FiniteElement< dim, spacedim > &fe_base, const FiniteElement< dim, spacedim > &fe_enriched, const FE_Nothing< dim, spacedim > &fe_nothing, hp::FECollection< dim, spacedim > &fe_collection)
unsigned int color_predicates(const DoFHandler< dim, spacedim > &mesh, const std::vector< predicate_function< dim, spacedim >> &predicates, std::vector< unsigned int > &predicate_colors)
void set_cellwise_color_set_and_fe_index(DoFHandler< dim, spacedim > &dof_handler, const std::vector< predicate_function< dim, spacedim >> &predicates, const std::vector< unsigned int > &predicate_colors, std::map< unsigned int, std::map< unsigned int, unsigned int >> &cellwise_color_predicate_map, std::vector< std::set< unsigned int >> &fe_sets)
void make_colorwise_enrichment_functions(const unsigned int n_colors, const std::vector< std::shared_ptr< Function< spacedim >>> &enrichments, const std::map< unsigned int, std::map< unsigned int, unsigned int >> &cellwise_color_predicate_map, std::vector< std::function< const Function< spacedim > *(const typename Triangulation< dim, spacedim >::cell_iterator &)>> &color_enrichments)
bool find_connection_between_subdomains(const DoFHandler< dim, spacedim > &dof_handler, const predicate_function< dim, spacedim > &predicate_1, const predicate_function< dim, spacedim > &predicate_2)
std::function< bool(const typename Triangulation< dim, spacedim >::cell_iterator &)> predicate_function
Definition: fe_enriched.h:700
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)
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)
std::vector< bool > compute_restriction_is_additive_flags(const std::vector< const FiniteElement< dim, spacedim > * > &fes, const std::vector< unsigned int > &multiplicities)
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)
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)
@ matrix
Contents is actually a matrix.
std::string to_string(const T &t)
Definition: patterns.h:2392
unsigned int color_sparsity_pattern(const SparsityPattern &sparsity_pattern, std::vector< unsigned int > &color_indices)
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:556
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:1593
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)
int(&) functions(const void *v1, const void *v2)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
unsigned short int fe_index
Definition: types.h:60
const hp::FECollection< dim, spacedim > & build_fe_collection(DoFHandler< dim, spacedim > &dof_handler)
const FiniteElement< dim, spacedim > & fe_enriched
Definition: fe_enriched.h:1108
const std::vector< predicate_function< dim, spacedim > > predicates
Definition: fe_enriched.h:1122
const std::vector< std::shared_ptr< Function< spacedim > > > enrichments
Definition: fe_enriched.h:1129
const FiniteElement< dim, spacedim > & fe_base
Definition: fe_enriched.h:1102