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