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