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