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