deal.II version GIT relicensing-1941-ga9013bea74 2024-10-08 06:50:00+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\}}\)
Loading...
Searching...
No Matches
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
25
26namespace 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(),
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(),
102 "Only dominating FE_Nothing can be used in FE_Enriched"));
103
105 fes[fe]->n_components() == n_comp_base,
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
135template <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
145template <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
164template <int dim, int spacedim>
166 const FiniteElement<dim, spacedim> *fe_base,
167 const std::vector<const FiniteElement<dim, spacedim> *> &fe_enriched,
168 const std::vector<std::vector<std::function<const Function<spacedim> *(
169 const typename Triangulation<dim, spacedim>::cell_iterator &)>>> &functions)
170 : FE_Enriched<dim, spacedim>(
171 internal::FE_Enriched::build_fes(fe_base, fe_enriched),
172 internal::FE_Enriched::build_multiplicities(functions),
173 functions)
174{}
175
176
177template <int dim, int spacedim>
179 const std::vector<const FiniteElement<dim, spacedim> *> &fes,
180 const std::vector<unsigned int> &multiplicities,
181 const std::vector<std::vector<std::function<const Function<spacedim> *(
182 const typename Triangulation<dim, spacedim>::cell_iterator &)>>> &functions)
183 : FiniteElement<dim, spacedim>(
184 FETools::Compositing::multiply_dof_numbers(fes, multiplicities, false),
185 FETools::Compositing::compute_restriction_is_additive_flags(
186 fes,
187 multiplicities),
188 FETools::Compositing::compute_nonzero_components(fes,
189 multiplicities,
190 false))
191 , enrichments(functions)
192 , is_enriched(internal::FE_Enriched::check_if_enriched(fes))
193 , fe_system(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),
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
209 Assert(base_no_mult_local_enriched_dofs.size() == this->n_base_elements(),
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(),
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
259template <int dim, int spacedim>
260std::vector<std::vector<std::function<const Function<spacedim> *(
263{
264 return enrichments;
265}
266
267
268template <int dim, int spacedim>
269double
271 const Point<dim> &p) const
272{
273 Assert(
274 !is_enriched,
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
281template <int dim, int spacedim>
282std::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
299template <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
322template <int dim, int spacedim>
323template <int dim_1>
324std::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
365template <int dim, int spacedim>
366std::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);
378 return setup_data(Utilities::dynamic_unique_cast<
380 std::move(data)),
381 update_flags,
382 quadrature[0]);
383}
384
385
386template <int dim, int spacedim>
387std::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);
398 return setup_data(Utilities::dynamic_unique_cast<
400 std::move(data)),
401 update_flags,
402 quadrature);
403}
404
405
406template <int dim, int spacedim>
407std::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);
416 return setup_data(Utilities::dynamic_unique_cast<
418 std::move(data)),
419 flags,
420 quadrature);
421}
422
423
424template <int dim, int spacedim>
425void
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
446 FETools::Compositing::build_cell_tables(this->system_to_base_table,
447 this->system_to_component_table,
448 this->component_to_base_table,
449 *this,
450 false);
451
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 {
488 this->adjust_line_dof_index_for_line_orientation_table =
489 fe_system->adjust_line_dof_index_for_line_orientation_table;
490 }
491}
492
493
494template <int dim, int spacedim>
495std::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
515template <int dim, int spacedim>
517FE_Enriched<dim, spacedim>::base_element(const unsigned int index) const
518{
519 return fe_system->base_element(index);
520}
521
522
523template <int dim, int spacedim>
524void
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,
532 &mapping_data,
533 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
535 &output_data) const
536{
537 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
539 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
540
541 // call FESystem's method to fill everything without enrichment function
542 fe_system->fill_fe_values(cell,
543 cell_similarity,
544 quadrature,
545 mapping,
546 mapping_internal,
547 mapping_data,
548 *fe_data.fesystem_data,
549 output_data);
550
551 if (is_enriched)
552 multiply_by_enrichment(
553 quadrature, fe_data, mapping_data, cell, output_data);
554}
555
556
557template <int dim, int spacedim>
558void
561 const unsigned int face_no,
562 const hp::QCollection<dim - 1> &quadrature,
563 const Mapping<dim, spacedim> &mapping,
564 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
566 &mapping_data,
567 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
569 &output_data) const
570{
571 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
573 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
574
575 AssertDimension(quadrature.size(), 1);
576
577 // call FESystem's method to fill everything without enrichment function
578 fe_system->fill_fe_face_values(cell,
579 face_no,
580 quadrature,
581 mapping,
582 mapping_internal,
583 mapping_data,
584 *fe_data.fesystem_data,
585 output_data);
586
587 if (is_enriched)
588 multiply_by_enrichment(
589 quadrature[0], fe_data, mapping_data, cell, output_data);
590}
591
592
593template <int dim, int spacedim>
594void
597 const unsigned int face_no,
598 const unsigned int sub_no,
599 const Quadrature<dim - 1> &quadrature,
600 const Mapping<dim, spacedim> &mapping,
601 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
603 &mapping_data,
604 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
606 &output_data) const
607{
608 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
610 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
611
612 // call FESystem's method to fill everything without enrichment function
613 fe_system->fill_fe_subface_values(cell,
614 face_no,
615 sub_no,
616 quadrature,
617 mapping,
618 mapping_internal,
619 mapping_data,
620 *fe_data.fesystem_data,
621 output_data);
622
623 if (is_enriched)
624 multiply_by_enrichment(
625 quadrature, fe_data, mapping_data, cell, output_data);
626}
627
628
629template <int dim, int spacedim>
630template <int dim_1>
631void
633 const Quadrature<dim_1> &quadrature,
634 const InternalData &fe_data,
636 &mapping_data,
639 &output_data) const
640{
641 // mapping_data will contain quadrature points on the real element.
642 // fe_internal is needed to get update flags
643 // finally, output_data should store all the results we need.
644
645 // Either dim_1==dim
646 // (fill_fe_values) or dim_1==dim-1
647 // (fill_fe_(sub)face_values)
648 Assert(dim_1 == dim || dim_1 == dim - 1, ExcInternalError());
649 const UpdateFlags flags = fe_data.update_each;
650
651 const unsigned int n_q_points = quadrature.size();
652
653 // First, populate output_data object (that shall hold everything requested
654 // such as shape value, gradients, hessians, etc) from each base element. That
655 // is almost identical to FESystem::compute_fill_one_base(), the difference
656 // being that we do it irrespectively of cell_similarity and use
657 // base_fe_data.update_flags
658
659 // TODO: do we need it only for dim_1 == dim (i.e. fill_fe_values)?
660 if (dim_1 == dim)
661 for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
662 {
663 const FiniteElement<dim, spacedim> &base_fe = base_element(base_no);
665 fe_data.get_fe_data(base_no);
667 spacedim>
668 &base_data = fe_data.get_fe_output_object(base_no);
669
670 const UpdateFlags base_flags = base_fe_data.update_each;
671
672 for (unsigned int system_index = 0;
673 system_index < this->n_dofs_per_cell();
674 ++system_index)
675 if (this->system_to_base_table[system_index].first.first == base_no)
676 {
677 const unsigned int base_index =
678 this->system_to_base_table[system_index].second;
679 Assert(base_index < base_fe.n_dofs_per_cell(),
681
682 // now copy. if the shape function is primitive, then there
683 // is only one value to be copied, but for non-primitive
684 // elements, there might be more values to be copied
685 //
686 // so, find out from which index to take this one value, and
687 // to which index to put
688 unsigned int out_index = 0;
689 for (unsigned int i = 0; i < system_index; ++i)
690 out_index += this->n_nonzero_components(i);
691 unsigned int in_index = 0;
692 for (unsigned int i = 0; i < base_index; ++i)
693 in_index += base_fe.n_nonzero_components(i);
694
695 // then loop over the number of components to be copied
696 Assert(this->n_nonzero_components(system_index) ==
697 base_fe.n_nonzero_components(base_index),
699 for (unsigned int s = 0;
700 s < this->n_nonzero_components(system_index);
701 ++s)
702 {
703 if (base_flags & update_values)
704 for (unsigned int q = 0; q < n_q_points; ++q)
705 output_data.shape_values[out_index + s][q] =
706 base_data.shape_values(in_index + s, q);
707
708 if (base_flags & update_gradients)
709 for (unsigned int q = 0; q < n_q_points; ++q)
710 output_data.shape_gradients[out_index + s][q] =
711 base_data.shape_gradients[in_index + s][q];
712
713 if (base_flags & update_hessians)
714 for (unsigned int q = 0; q < n_q_points; ++q)
715 output_data.shape_hessians[out_index + s][q] =
716 base_data.shape_hessians[in_index + s][q];
717 }
718 }
719 }
720
721 Assert(base_no_mult_local_enriched_dofs.size() == fe_data.enrichment.size(),
722 ExcDimensionMismatch(base_no_mult_local_enriched_dofs.size(),
723 fe_data.enrichment.size()));
724 // calculate hessians, gradients and values for each function
725 for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
726 {
727 Assert(
728 base_no_mult_local_enriched_dofs[base_no].size() ==
729 fe_data.enrichment[base_no].size(),
730 ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no].size(),
731 fe_data.enrichment[base_no].size()));
732 for (unsigned int m = 0;
733 m < base_no_mult_local_enriched_dofs[base_no].size();
734 m++)
735 {
736 // Avoid evaluating quadrature points if no dofs are assigned. This
737 // happens when FE_Nothing is used together with other FE (i.e. FE_Q)
738 // as enrichments.
739 if (base_no_mult_local_enriched_dofs[base_no][m].empty())
740 continue;
741
742 Assert(enrichments[base_no - 1][m](cell) != nullptr,
744 "The pointer to the enrichment function is not set"));
745
746 Assert(enrichments[base_no - 1][m](cell)->n_components == 1,
748 "Only scalar-valued enrichment functions are allowed"));
749
750 if (flags & update_hessians)
751 {
752 Assert(fe_data.enrichment[base_no][m].hessians.size() ==
753 n_q_points,
755 fe_data.enrichment[base_no][m].hessians.size(),
756 n_q_points));
757 for (unsigned int q = 0; q < n_q_points; ++q)
758 fe_data.enrichment[base_no][m].hessians[q] =
759 enrichments[base_no - 1][m](cell)->hessian(
760 mapping_data.quadrature_points[q]);
761 }
762
763 if (flags & update_gradients)
764 {
765 Assert(fe_data.enrichment[base_no][m].gradients.size() ==
766 n_q_points,
768 fe_data.enrichment[base_no][m].gradients.size(),
769 n_q_points));
770 for (unsigned int q = 0; q < n_q_points; ++q)
771 fe_data.enrichment[base_no][m].gradients[q] =
772 enrichments[base_no - 1][m](cell)->gradient(
773 mapping_data.quadrature_points[q]);
774 }
775
776 if (flags & update_values)
777 {
778 Assert(fe_data.enrichment[base_no][m].values.size() == n_q_points,
780 fe_data.enrichment[base_no][m].values.size(),
781 n_q_points));
782 for (unsigned int q = 0; q < n_q_points; ++q)
783 fe_data.enrichment[base_no][m].values[q] =
784 enrichments[base_no - 1][m](cell)->value(
785 mapping_data.quadrature_points[q]);
786 }
787 }
788 }
789
790 // Finally, update the standard data stored in output_data
791 // by expanding the product rule for enrichment function.
792 // note that the order if important, namely
793 // output_data.shape_XYZ contains values of standard FEM and we overwrite
794 // it with the updated one in the following order: hessians -> gradients ->
795 // values
796 if (flags & update_hessians)
797 {
798 for (unsigned int base_no = 1; base_no < this->n_base_elements();
799 base_no++)
800 {
801 for (unsigned int m = 0;
802 m < base_no_mult_local_enriched_dofs[base_no].size();
803 m++)
804 for (unsigned int i = 0;
805 i < base_no_mult_local_enriched_dofs[base_no][m].size();
806 i++)
807 {
808 const unsigned int enriched_dof =
809 base_no_mult_local_enriched_dofs[base_no][m][i];
810 for (unsigned int q = 0; q < n_q_points; ++q)
811 {
812 const Tensor<2, spacedim> grad_grad = outer_product(
813 output_data.shape_gradients[enriched_dof][q],
814 fe_data.enrichment[base_no][m].gradients[q]);
815 const Tensor<2, spacedim, double> sym_grad_grad =
816 symmetrize(grad_grad) * 2.0; // symmetrize does [s+s^T]/2
817
818 output_data.shape_hessians[enriched_dof][q] *=
819 fe_data.enrichment[base_no][m].values[q];
820 output_data.shape_hessians[enriched_dof][q] +=
821 sym_grad_grad +
822 output_data.shape_values[enriched_dof][q] *
823 fe_data.enrichment[base_no][m].hessians[q];
824 }
825 }
826 }
827 }
828
829 if (flags & update_gradients)
830 for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
831 {
832 for (unsigned int m = 0;
833 m < base_no_mult_local_enriched_dofs[base_no].size();
834 m++)
835 for (unsigned int i = 0;
836 i < base_no_mult_local_enriched_dofs[base_no][m].size();
837 i++)
838 {
839 const unsigned int enriched_dof =
840 base_no_mult_local_enriched_dofs[base_no][m][i];
841 for (unsigned int q = 0; q < n_q_points; ++q)
842 {
843 output_data.shape_gradients[enriched_dof][q] *=
844 fe_data.enrichment[base_no][m].values[q];
845 output_data.shape_gradients[enriched_dof][q] +=
846 output_data.shape_values[enriched_dof][q] *
847 fe_data.enrichment[base_no][m].gradients[q];
848 }
849 }
850 }
851
852 if (flags & update_values)
853 for (unsigned int base_no = 1; base_no < this->n_base_elements(); ++base_no)
854 {
855 for (unsigned int m = 0;
856 m < base_no_mult_local_enriched_dofs[base_no].size();
857 m++)
858 for (unsigned int i = 0;
859 i < base_no_mult_local_enriched_dofs[base_no][m].size();
860 i++)
861 {
862 const unsigned int enriched_dof =
863 base_no_mult_local_enriched_dofs[base_no][m][i];
864 for (unsigned int q = 0; q < n_q_points; ++q)
865 {
866 output_data.shape_values[enriched_dof][q] *=
867 fe_data.enrichment[base_no][m].values[q];
868 }
869 }
870 }
871}
872
873
874template <int dim, int spacedim>
877{
878 return *fe_system;
879}
880
881
882template <int dim, int spacedim>
883bool
888
889
890template <int dim, int spacedim>
891void
893 const FiniteElement<dim, spacedim> &source,
894 FullMatrix<double> &matrix,
895 const unsigned int face_no) const
896{
897 if (const FE_Enriched<dim, spacedim> *fe_enr_other =
898 dynamic_cast<const FE_Enriched<dim, spacedim> *>(&source))
899 {
900 fe_system->get_face_interpolation_matrix(fe_enr_other->get_fe_system(),
901 matrix,
902 face_no);
903 }
904 else
905 {
907 false,
908 (typename FiniteElement<dim,
909 spacedim>::ExcInterpolationNotImplemented()));
910 }
911}
912
913
914template <int dim, int spacedim>
915void
917 const FiniteElement<dim, spacedim> &source,
918 const unsigned int subface,
919 FullMatrix<double> &matrix,
920 const unsigned int face_no) const
921{
922 if (const FE_Enriched<dim, spacedim> *fe_enr_other =
923 dynamic_cast<const FE_Enriched<dim, spacedim> *>(&source))
924 {
925 fe_system->get_subface_interpolation_matrix(fe_enr_other->get_fe_system(),
926 subface,
927 matrix,
928 face_no);
929 }
930 else
931 {
933 false,
934 (typename FiniteElement<dim,
935 spacedim>::ExcInterpolationNotImplemented()));
936 }
937}
938
939
940template <int dim, int spacedim>
941std::vector<std::pair<unsigned int, unsigned int>>
943 const FiniteElement<dim, spacedim> &fe_other) const
944{
945 if (const FE_Enriched<dim, spacedim> *fe_enr_other =
946 dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
947 {
948 return fe_system->hp_vertex_dof_identities(fe_enr_other->get_fe_system());
949 }
950 else
951 {
953 return std::vector<std::pair<unsigned int, unsigned int>>();
954 }
955}
956
957
958template <int dim, int spacedim>
959std::vector<std::pair<unsigned int, unsigned int>>
961 const FiniteElement<dim, spacedim> &fe_other) const
962{
963 if (const FE_Enriched<dim, spacedim> *fe_enr_other =
964 dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
965 {
966 return fe_system->hp_line_dof_identities(fe_enr_other->get_fe_system());
967 }
968 else
969 {
971 return std::vector<std::pair<unsigned int, unsigned int>>();
972 }
973}
974
975
976template <int dim, int spacedim>
977std::vector<std::pair<unsigned int, unsigned int>>
979 const FiniteElement<dim, spacedim> &fe_other,
980 const unsigned int face_no) const
981{
982 if (const FE_Enriched<dim, spacedim> *fe_enr_other =
983 dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
984 {
985 return fe_system->hp_quad_dof_identities(fe_enr_other->get_fe_system(),
986 face_no);
987 }
988 else
989 {
991 return std::vector<std::pair<unsigned int, unsigned int>>();
992 }
993}
994
995
996template <int dim, int spacedim>
999 const FiniteElement<dim, spacedim> &fe_other,
1000 const unsigned int codim) const
1001{
1002 Assert(codim <= dim, ExcImpossibleInDim(dim));
1003
1004 // vertex/line/face/cell domination
1005 // --------------------------------
1006 // need to decide which element constrain another.
1007 // for example Q(2) dominate Q(4) and thus some DoFs of Q(4) will be
1008 // constrained. If we have Q(2) and Q(4)+POU, then it's clear that Q(2)
1009 // dominates, namely our DoFs will be constrained to make field continuous.
1010 // However, we need to check for situations like Q(4) vs Q(2)+POU.
1011 // In that case the domination for the underlying FEs should be the other way,
1012 // but this implies that we can't constrain POU dofs to make the field
1013 // continuous. In that case, throw an error
1014
1015 // if it's also enriched, do domination based on each one's FESystem
1016 if (const FE_Enriched<dim, spacedim> *fe_enr_other =
1017 dynamic_cast<const FE_Enriched<dim, spacedim> *>(&fe_other))
1018 {
1019 return fe_system->compare_for_domination(fe_enr_other->get_fe_system(),
1020 codim);
1021 }
1022 else
1023 {
1026 }
1027}
1028
1029
1030template <int dim, int spacedim>
1031const FullMatrix<double> &
1033 const unsigned int child,
1034 const RefinementCase<dim> &refinement_case) const
1035{
1036 return fe_system->get_prolongation_matrix(child, refinement_case);
1037}
1038
1039
1040template <int dim, int spacedim>
1041const FullMatrix<double> &
1043 const unsigned int child,
1044 const RefinementCase<dim> &refinement_case) const
1045{
1046 return fe_system->get_restriction_matrix(child, refinement_case);
1047}
1048
1049
1050/* ----------------------- FESystem::InternalData ------------------- */
1051
1052
1053template <int dim, int spacedim>
1055 std::unique_ptr<typename FESystem<dim, spacedim>::InternalData> fesystem_data)
1056 : fesystem_data(std::move(fesystem_data))
1057{}
1058
1059
1060template <int dim, int spacedim>
1063 const unsigned int base_no) const
1064{
1065 return fesystem_data->get_fe_data(base_no);
1066}
1067
1068
1069template <int dim, int spacedim>
1072 const unsigned int base_no) const
1073{
1074 return fesystem_data->get_fe_output_object(base_no);
1075}
1076
1077
1078namespace ColorEnriched
1079{
1080 namespace internal
1081 {
1082 template <int dim, int spacedim>
1083 bool
1085 const DoFHandler<dim, spacedim> &dof_handler,
1086 const predicate_function<dim, spacedim> &predicate_1,
1087 const predicate_function<dim, spacedim> &predicate_2)
1088 {
1089 // Use a vector to mark vertices
1090 std::vector<bool> vertices_subdomain_1(
1091 dof_handler.get_triangulation().n_vertices(), false);
1092
1093 // Mark vertices that belong to cells in subdomain 1
1094 for (const auto &cell : dof_handler.active_cell_iterators())
1095 if (predicate_1(cell)) // True ==> part of subdomain 1
1096 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1097 vertices_subdomain_1[cell->vertex_index(v)] = true;
1098
1099 // Find if cells in subdomain 2 and subdomain 1 share vertices.
1100 for (const auto &cell : dof_handler.active_cell_iterators())
1101 if (predicate_2(cell)) // True ==> part of subdomain 2
1102 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1103 if (vertices_subdomain_1[cell->vertex_index(v)] == true)
1104 {
1105 return true;
1106 }
1107 return false;
1108 }
1109
1110
1111
1112 template <int dim, int spacedim>
1113 unsigned int
1115 const DoFHandler<dim, spacedim> &mesh,
1116 const std::vector<predicate_function<dim, spacedim>> &predicates,
1117 std::vector<unsigned int> &predicate_colors)
1118 {
1119 const unsigned int num_indices = predicates.size();
1120
1121 // Use sparsity pattern to represent connections between subdomains.
1122 // Each predicate (i.e a subdomain) is a node in the graph.
1124 dsp.reinit(num_indices, num_indices);
1125
1126 /*
1127 * Find connections between subdomains taken two at a time.
1128 * If the connection exists, add it to a graph object represented
1129 * by dynamic sparsity pattern.
1130 */
1131 for (unsigned int i = 0; i < num_indices; ++i)
1132 for (unsigned int j = i + 1; j < num_indices; ++j)
1133 if (internal::find_connection_between_subdomains(mesh,
1134 predicates[i],
1135 predicates[j]))
1136 dsp.add(i, j);
1137
1138 dsp.symmetrize();
1139
1140 // Copy dynamic sparsity pattern to sparsity pattern needed by
1141 // coloring function
1142 SparsityPattern sp_graph;
1143 sp_graph.copy_from(dsp);
1144 predicate_colors.resize(num_indices);
1145
1146 // Assign each predicate with a color and return number of colors
1147 return SparsityTools::color_sparsity_pattern(sp_graph, predicate_colors);
1148 }
1149
1150
1151
1152 template <int dim, int spacedim>
1153 void
1155 DoFHandler<dim, spacedim> &dof_handler,
1156 const std::vector<predicate_function<dim, spacedim>> &predicates,
1157 const std::vector<unsigned int> &predicate_colors,
1158 std::map<unsigned int, std::map<unsigned int, unsigned int>>
1159 &cellwise_color_predicate_map,
1160 std::vector<std::set<unsigned int>> &fe_sets)
1161 {
1162 // clear output variables first
1163 fe_sets.clear();
1164 cellwise_color_predicate_map.clear();
1165
1166 /*
1167 * Add first element of fe_sets which is empty by default. This means that
1168 * the default, FE index = 0 is associated with an empty set, since no
1169 * predicate is active in these regions.
1170 */
1171 fe_sets.resize(1);
1172
1173 /*
1174 * Loop through cells and find set of predicate colors associated
1175 * with the cell. As an example, a cell with an FE index associated
1176 * with colors {a,b} means that predicates active in the cell have
1177 * colors a or b.
1178 *
1179 * Create new active FE index in case of the color
1180 * set is not already listed in fe_sets. If the set already exists,
1181 * find index of the set in fe_sets. In either case, use the id in
1182 * fe_sets to modify cell->active_fe_index.
1183 *
1184 * Associate each cell_id with a set of pairs. The pair represents
1185 * predicate color and the active predicate with that color.
1186 * Each color can only correspond to a single predicate since
1187 * predicates with the same color correspond to disjoint domains.
1188 * This is what the graph coloring in color_predicates
1189 * function ensures. The number of pairs is equal to the number
1190 * of predicates active in the given cell.
1191 */
1192 unsigned int map_index = 0;
1193 for (const auto &cell : dof_handler.active_cell_iterators())
1194 {
1195 // set default FE index ==> no enrichment and no active predicates
1196 cell->set_active_fe_index(0);
1197
1198 // Give each cell a unique id, which the cellwise_color_predicate_map
1199 // can later use to associate a colors of active predicates with
1200 // the actual predicate id.
1201 //
1202 // When the grid is refined, material id is inherited to
1203 // children cells. So, the cellwise_color_predicate_map stays
1204 // relevant.
1205 cell->set_material_id(map_index);
1206 std::set<unsigned int> color_list;
1207
1208 // loop through active predicates for the cell and insert map.
1209 // Eg: if the cell with material id 100 has active
1210 // predicates 4 (color = 1) and 5 (color = 2), the map will insert
1211 // pairs (1, 4) and (2, 5) at key 100 (i.e unique id of cell is
1212 // mapped with a map which associates color with predicate id)
1213 // Note that color list for the cell would be {1,2}.
1214 std::map<unsigned int, unsigned int> &cell_map =
1215 cellwise_color_predicate_map[map_index];
1216 for (unsigned int i = 0; i < predicates.size(); ++i)
1217 {
1218 if (predicates[i](cell))
1219 {
1220 /*
1221 * create a pair predicate color and predicate id and add it
1222 * to a map associated with each enriched cell
1223 */
1224 auto ret = cell_map.insert(
1225 std::pair<unsigned int, unsigned int>(predicate_colors[i],
1226 i));
1227
1228 AssertThrow(ret.second == 1,
1229 ExcMessage(
1230 "Only one enrichment function per color"));
1231
1232 color_list.insert(predicate_colors[i]);
1233 }
1234 }
1235
1236
1237 /*
1238 * check if color combination is already added.
1239 * If already added, set the active FE index based on
1240 * its index in the fe_sets. If the combination doesn't
1241 * exist, add the set to fe_sets and once again set the
1242 * active FE index as last index in fe_sets.
1243 *
1244 * Eg: if cell has color list {1,2} associated and
1245 * fe_sets = { {}, {1}, {2} } for now, we need to add a new set {1,2}
1246 * to fe_sets and a new active FE index 3 because 0 to 2 FE indices
1247 * are already taken by existing sets in fe_sets.
1248 */
1249 if (!color_list.empty())
1250 {
1251 const auto it =
1252 std::find(fe_sets.begin(), fe_sets.end(), color_list);
1253 // when entry is not found
1254 if (it == fe_sets.end())
1255 {
1256 fe_sets.push_back(color_list);
1257 cell->set_active_fe_index(fe_sets.size() - 1);
1258 }
1259 // when entry is found
1260 else
1261 {
1262 cell->set_active_fe_index(std::distance(fe_sets.begin(), it));
1263 }
1264 }
1265 /*
1266 * map_index is used to identify cells and should be unique. The
1267 * index is stored in the material_id of the cell and hence
1268 * stays relevant even when the cells are refined (which is
1269 * why cell_id is not used).
1270 * Two distant cells can have the same set of colors but different
1271 * enrichment functions can be associated with any given
1272 * color. So, in order to figure which enrichment function
1273 * belongs to a color, we use a map that uses this index.
1274 */
1275 ++map_index;
1276 }
1277
1278
1279 /*
1280 * Treat interface between enriched cells specially,
1281 * until #1496 (https://github.com/dealii/dealii/issues/1496) is resolved.
1282 * Each time we build constraints at the interface between two different
1283 * FE_Enriched, we look for the least dominating FE of their common
1284 * subspace via hp::FECollection::find_dominating_fe_extended().
1285 * If we don't take further actions, we may find a dominating FE that is
1286 * too restrictive, i.e. enriched FE consisting of only FE_Nothing. New
1287 * elements needs to be added to FECollection object to help find the
1288 * correct enriched FE underlying the spaces in the adjacent cells. This
1289 * is done by creating an appropriate set in fe_sets and a call to the
1290 * function make_fe_collection_from_colored_enrichments at a later stage.
1291 *
1292 * Consider a domain with three predicates and hence with three different
1293 * enrichment functions. Let the enriched finite element of a cell with
1294 * enrichment functions 1 and 2 be represented by [1 1 0], with the last
1295 * entry as zero since the 3rd enrichment function is not relevant for
1296 * the cell. If the interface has enriched FE [1 0 1] and [0 1 1]
1297 * on adjacent cells, an enriched FE [0 0 1] should exist and is
1298 * found as the least dominating finite element for the two cells by
1299 * DoFTools::make_hanging_node_constraints, using the above mentioned
1300 * hp::FECollection functions. Denoting the FE set in adjacent cells as
1301 * {1,3} and {2,3}, this implies that an FE set {3} needs to be added!
1302 * Based on the predicate configuration, this may not be automatically
1303 * done without the following special treatment.
1304 */
1305
1306 // loop through faces
1307 for (const auto &cell : dof_handler.active_cell_iterators())
1308 {
1309 const unsigned int fe_index = cell->active_fe_index();
1310 const std::set<unsigned int> fe_set = fe_sets.at(fe_index);
1311 for (const unsigned int face : GeometryInfo<dim>::face_indices())
1312 {
1313 // cell shouldn't be at the boundary and
1314 // neighboring cell is not already visited (to avoid visiting
1315 // same face twice). Note that the cells' material ids are
1316 // labeled according to their order in dof_handler previously.
1317 if (!cell->at_boundary(face) &&
1318 cell->material_id() < cell->neighbor(face)->material_id())
1319 {
1320 const auto nbr_fe_index =
1321 cell->neighbor(face)->active_fe_index();
1322
1323 // find corresponding FE set
1324 const auto nbr_fe_set = fe_sets.at(nbr_fe_index);
1325
1326 // find intersection of the FE sets: fe_set and nbr_fe_set
1327 std::set<unsigned int> intersection_set;
1328 std::set_intersection(
1329 fe_set.begin(),
1330 fe_set.end(),
1331 nbr_fe_set.begin(),
1332 nbr_fe_set.end(),
1333 std::inserter(intersection_set, intersection_set.begin()));
1334
1335 // add only the new sets
1336 if (!intersection_set.empty())
1337 {
1338 const auto it = std::find(fe_sets.begin(),
1339 fe_sets.end(),
1340 intersection_set);
1341 // add the set if it is not found
1342 if (it == fe_sets.end())
1343 {
1344 fe_sets.push_back(intersection_set);
1345 }
1346 }
1347 }
1348 }
1349 }
1350 }
1351
1352
1353
1354 template <int dim, int spacedim>
1355 void
1357 const unsigned int n_colors,
1358 const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments,
1359 const std::map<unsigned int, std::map<unsigned int, unsigned int>>
1360 &cellwise_color_predicate_map,
1361 std::vector<std::function<const Function<spacedim> *(
1363 &color_enrichments)
1364 {
1365 color_enrichments.clear();
1366
1367 // Each color should be associated with a single enrichment function
1368 // called color enrichment function which calls the correct enrichment
1369 // function for a given cell.
1370 //
1371 // Assume that a cell has a active predicates with ids 4 (color = 1) and
1372 // 5 (color = 2). cellwise_color_predicate_map has this information,
1373 // provided we know the material id.
1374 //
1375 // The constructed color_enrichments is such that
1376 // color_enrichments[1](cell) will return return a pointer to
1377 // function with id=4, i.e. enrichments[4].
1378 // In other words, using the previously collected information in
1379 // this function we translate a vector of user provided enrichment
1380 // functions into a vector of functions suitable for FE_Enriched class.
1381 color_enrichments.resize(n_colors);
1382 for (unsigned int i = 0; i < n_colors; ++i)
1383 {
1384 color_enrichments[i] =
1386 &cell) {
1387 const unsigned int id = cell->material_id();
1388
1389 /*
1390 * i'th color_enrichment function corresponds to (i+1)'th color.
1391 * Since FE_Enriched takes function pointers, we return a
1392 * function pointer.
1393 */
1394 return enrichments[cellwise_color_predicate_map.at(id).at(i + 1)]
1395 .get();
1396 };
1397 }
1398 }
1399
1400
1401
1402 template <int dim, int spacedim>
1403 void
1405 const unsigned int n_colors,
1406 const std::vector<std::set<unsigned int>> &fe_sets,
1407 const std::vector<std::function<const Function<spacedim> *(
1409 &color_enrichments,
1410 const FiniteElement<dim, spacedim> &fe_base,
1411 const FiniteElement<dim, spacedim> &fe_enriched,
1412 const FE_Nothing<dim, spacedim> &fe_nothing,
1413 hp::FECollection<dim, spacedim> &fe_collection)
1414 {
1415 // define dummy function which is associated with FE_Nothing
1416 const std::function<const Function<spacedim> *(
1418 dummy_function =
1420 -> const Function<spacedim> * {
1421 AssertThrow(false,
1422 ExcMessage("Called enrichment function for FE_Nothing"));
1423 return nullptr;
1424 };
1425
1426
1427 // loop through color sets and create FE_enriched element for each
1428 // of them provided before calling this function, we have color
1429 // enrichment function associated with each color.
1430 for (const auto &fe_set : fe_sets)
1431 {
1432 std::vector<const FiniteElement<dim, spacedim> *> vec_fe_enriched(
1433 n_colors, &fe_nothing);
1434 std::vector<std::vector<std::function<const Function<spacedim> *(
1436 functions(n_colors, {dummy_function});
1437
1438 for (const unsigned int color_id : fe_set)
1439 {
1440 // Given a color id, corresponding color enrichment
1441 // function is at index id-1 because color_enrichments are
1442 // indexed from zero and colors are indexed from 1.
1443 const unsigned int ind = color_id - 1;
1444
1445 AssertIndexRange(ind, vec_fe_enriched.size());
1446 AssertIndexRange(ind, functions.size());
1447 AssertIndexRange(ind, color_enrichments.size());
1448
1449 // Assume an active predicate colors {1,2} for a cell.
1450 // We then need to create a vector of FE enriched elements
1451 // with vec_fe_enriched[0] = vec_fe_enriched[1] = &fe_enriched
1452 // which can later be associated with enrichment functions.
1453 vec_fe_enriched[ind] = &fe_enriched;
1454
1455 // color_set_id'th color function is (color_set_id-1)
1456 // element of color wise enrichments
1457 functions[ind][0] = color_enrichments[ind];
1458 }
1459
1460 AssertDimension(vec_fe_enriched.size(), functions.size());
1461
1462 FE_Enriched<dim, spacedim> fe_component(&fe_base,
1463 vec_fe_enriched,
1464 functions);
1465 fe_collection.push_back(fe_component);
1466 }
1467 }
1468 } // namespace internal
1469
1470
1471
1472 template <int dim, int spacedim>
1474 const FiniteElement<dim, spacedim> &fe_base,
1475 const FiniteElement<dim, spacedim> &fe_enriched,
1476 const std::vector<predicate_function<dim, spacedim>> &predicates,
1477 const std::vector<std::shared_ptr<Function<spacedim>>> &enrichments)
1478 : fe_base(fe_base)
1479 , fe_enriched(fe_enriched)
1480 , fe_nothing(fe_base.n_components(), true)
1481 , predicates(predicates)
1482 , enrichments(enrichments)
1483 , n_colors(numbers::invalid_unsigned_int)
1484 {
1485 AssertDimension(predicates.size(), enrichments.size());
1487 AssertThrow(predicates.size() > 0,
1488 ExcMessage("Number of predicates should be positive"));
1489 }
1490
1491
1492
1493 template <int dim, int spacedim>
1496 DoFHandler<dim, spacedim> &dof_handler)
1497 {
1498 // color the predicates based on connections between corresponding
1499 // subdomains
1500 n_colors =
1501 internal::color_predicates(dof_handler, predicates, predicate_colors);
1502
1503 // create color maps and color list for each cell
1505 predicates,
1506 predicate_colors,
1507 cellwise_color_predicate_map,
1508 fe_sets);
1509 // setup color wise enrichment functions
1510 // i'th function corresponds to (i+1) color!
1511 internal::make_colorwise_enrichment_functions<dim, spacedim>(
1512 n_colors, enrichments, cellwise_color_predicate_map, color_enrichments);
1513
1514 // make FE_Collection
1516 fe_sets,
1517 color_enrichments,
1518 fe_base,
1519 fe_enriched,
1520 fe_nothing,
1521 fe_collection);
1522
1523 return fe_collection;
1524 }
1525} // namespace ColorEnriched
1526
1527
1528// explicit instantiations
1529#include "fe_enriched.inst"
1530
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:2553
unsigned int n_nonzero_components(const unsigned int i) const
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
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:308
void push_back(const FiniteElement< dim, spacedim > &new_fe)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
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.
#define DEAL_II_NOT_IMPLEMENTED()
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:555
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition utilities.h:1620
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)