37 std::vector<unsigned int>
38 build_multiplicities(
const std::vector<std::vector<T>> &
functions)
40 std::vector<unsigned int> multiplicities;
41 multiplicities.push_back(1);
42 for (
unsigned int i = 0; i <
functions.size(); i++)
43 multiplicities.push_back(
functions[i].size());
45 return multiplicities;
52 template <
int dim,
int spacedim>
53 std::vector<const FiniteElement<dim, spacedim> *>
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]);
71 template <
int dim,
int spacedim>
75 const std::vector<unsigned int> & multiplicities,
77 const typename ::Triangulation<dim, spacedim>::cell_iterator
83 "FEs and multiplicities should have the same size"));
91 const unsigned int n_comp_base = fes[0]->n_components();
94 for (
unsigned int fe = 1; fe < fes.size(); fe++)
102 "Only dominating FE_Nothing can be used in FE_Enriched"));
107 "All elements must have the same number of components"));
117 template <
int dim,
int spacedim>
123 for (
unsigned int fe = 1; fe < fes.size(); fe++)
135 template <
int dim,
int spacedim>
145 template <
int dim,
int spacedim>
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 &)>>>(
156 std::vector<std::function<const
Function<spacedim> *(
157 const typename
Triangulation<dim, spacedim>::cell_iterator &)>>(
159 [=](const typename
Triangulation<dim, spacedim>::cell_iterator &)
160 -> const
Function<spacedim> * {
return enrichment_function; })))
164 template <
int dim,
int spacedim>
177 template <
int dim,
int spacedim>
180 const std::vector<unsigned int> & multiplicities,
197 Assert(internal::FE_Enriched::consistency_check(fes,
207 for (
unsigned int fe = 1; fe < fes.size(); fe++)
212 this->n_base_elements()));
215 for (
unsigned int system_index = 0; system_index < this->
dofs_per_cell;
218 const unsigned int base_no =
223 const unsigned int base_m =
228 "Size mismatch for base_no_mult_local_enriched_dofs: "
244 for (
unsigned int base_no = 1;
248 for (
unsigned int m = 0;
260 template <
int dim,
int spacedim>
261 const std::vector<std::vector<std::function<const Function<spacedim> *(
269 template <
int dim,
int spacedim>
277 "For enriched finite elements shape_value() can not be defined on the reference element."));
278 return fe_system->shape_value(i, p);
282 template <
int dim,
int spacedim>
283 std::unique_ptr<FiniteElement<dim, spacedim>>
286 std::vector<const FiniteElement<dim, spacedim> *> fes;
287 std::vector<unsigned int> multiplicities;
289 for (
unsigned int i = 0; i < this->n_base_elements(); i++)
291 fes.push_back(&base_element(i));
292 multiplicities.push_back(this->element_multiplicity(i));
295 return std::unique_ptr<FE_Enriched<dim, spacedim>>(
300 template <
int dim,
int spacedim>
304 UpdateFlags out = fe_system->requires_update_flags(flags);
323 template <
int dim,
int spacedim>
325 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
333 auto update_each_flags = fes_data->update_each;
334 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
335 data_ptr = std_cxx14::make_unique<InternalData>(std::move(fes_data));
342 data.enrichment.resize(this->n_base_elements());
344 const unsigned int n_q_points = quadrature.
size();
346 for (
unsigned int base = 0; base < this->n_base_elements(); ++base)
348 data.enrichment[base].resize(this->element_multiplicity(base));
349 for (
unsigned int m = 0; m < this->element_multiplicity(base); ++m)
352 data.enrichment[base][m].values.resize(n_q_points);
355 data.enrichment[base][m].gradients.resize(n_q_points);
358 data.enrichment[base][m].hessians.resize(n_q_points);
366 template <
int dim,
int spacedim>
367 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
376 fe_system->get_face_data(update_flags, mapping, quadrature, output_data);
385 template <
int dim,
int spacedim>
386 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
396 fe_system->get_subface_data(update_flags, mapping, quadrature, output_data);
405 template <
int dim,
int spacedim>
406 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
414 auto data = fe_system->get_data(flags, mapping, quadrature, output_data);
423 template <
int dim,
int spacedim>
427 const std::vector<unsigned int> & multiplicities)
429 Assert(fes.size() == multiplicities.size(),
434 this->base_to_block_indices.reinit(0, 0);
436 for (
unsigned int i = 0; i < fes.size(); i++)
437 if (multiplicities[i] > 0)
438 this->base_to_block_indices.push_back(multiplicities[i]);
443 this->system_to_component_table.resize(this->dofs_per_cell);
444 this->face_system_to_component_table.resize(this->dofs_per_face);
447 this->system_to_component_table,
448 this->component_to_base_table,
453 this->face_system_to_base_table,
454 this->face_system_to_component_table,
463 this->interface_constraints = fe_system->interface_constraints;
472 this->unit_support_points = fe_system->unit_support_points;
473 this->unit_face_support_points = fe_system->unit_face_support_points;
478 this->adjust_line_dof_index_for_line_orientation_table =
479 fe_system->adjust_line_dof_index_for_line_orientation_table;
484 template <
int dim,
int spacedim>
488 std::ostringstream namebuf;
491 for (
unsigned int i = 0; i < this->n_base_elements(); ++i)
493 namebuf << base_element(i).get_name();
494 if (this->element_multiplicity(i) != 1)
495 namebuf <<
'^' << this->element_multiplicity(i);
496 if (i != this->n_base_elements() - 1)
501 return namebuf.str();
505 template <
int dim,
int spacedim>
509 return fe_system->base_element(index);
513 template <
int dim,
int spacedim>
521 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
533 fe_system->fill_fe_values(cell,
543 multiply_by_enrichment(
544 quadrature, fe_data, mapping_data, cell, output_data);
548 template <
int dim,
int spacedim>
552 const unsigned int face_no,
556 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
568 fe_system->fill_fe_face_values(cell,
578 multiply_by_enrichment(
579 quadrature, fe_data, mapping_data, cell, output_data);
583 template <
int dim,
int spacedim>
587 const unsigned int face_no,
588 const unsigned int sub_no,
592 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
604 fe_system->fill_fe_subface_values(cell,
615 multiply_by_enrichment(
616 quadrature, fe_data, mapping_data, cell, output_data);
620 template <
int dim,
int spacedim>
642 const unsigned int n_q_points = quadrature.
size();
652 for (
unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
663 for (
unsigned int system_index = 0; system_index < this->dofs_per_cell;
665 if (this->system_to_base_table[system_index].
first.first == base_no)
667 const unsigned int base_index =
668 this->system_to_base_table[system_index].second;
677 unsigned int out_index = 0;
678 for (
unsigned int i = 0; i < system_index; ++i)
679 out_index += this->n_nonzero_components(i);
680 unsigned int in_index = 0;
681 for (
unsigned int i = 0; i < base_index; ++i)
685 Assert(this->n_nonzero_components(system_index) ==
688 for (
unsigned int s = 0;
689 s < this->n_nonzero_components(system_index);
693 for (
unsigned int q = 0; q < n_q_points; ++q)
695 base_data.shape_values(in_index + s, q);
698 for (
unsigned int q = 0; q < n_q_points; ++q)
700 base_data.shape_gradients[in_index + s][q];
703 for (
unsigned int q = 0; q < n_q_points; ++q)
705 base_data.shape_hessians[in_index + s][q];
710 Assert(base_no_mult_local_enriched_dofs.size() == fe_data.
enrichment.size(),
714 for (
unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
717 base_no_mult_local_enriched_dofs[base_no].size() ==
721 for (
unsigned int m = 0;
722 m < base_no_mult_local_enriched_dofs[base_no].size();
728 if (base_no_mult_local_enriched_dofs[base_no][m].size() == 0)
731 Assert(enrichments[base_no - 1][m](cell) !=
nullptr,
733 "The pointer to the enrichment function is not set"));
737 "Only scalar-valued enrichment functions are allowed"));
744 fe_data.
enrichment[base_no][m].hessians.size(),
746 for (
unsigned int q = 0; q < n_q_points; q++)
748 enrichments[base_no - 1][m](cell)->hessian(
757 fe_data.
enrichment[base_no][m].gradients.size(),
759 for (
unsigned int q = 0; q < n_q_points; q++)
761 enrichments[base_no - 1][m](cell)->gradient(
771 for (
unsigned int q = 0; q < n_q_points; q++)
773 enrichments[base_no - 1][m](cell)->value(
787 for (
unsigned int base_no = 1; base_no < this->n_base_elements();
790 for (
unsigned int m = 0;
791 m < base_no_mult_local_enriched_dofs[base_no].size();
793 for (
unsigned int i = 0;
794 i < base_no_mult_local_enriched_dofs[base_no][m].size();
797 const unsigned int enriched_dof =
798 base_no_mult_local_enriched_dofs[base_no][m][i];
799 for (
unsigned int q = 0; q < n_q_points; ++q)
819 for (
unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
821 for (
unsigned int m = 0;
822 m < base_no_mult_local_enriched_dofs[base_no].size();
824 for (
unsigned int i = 0;
825 i < base_no_mult_local_enriched_dofs[base_no][m].size();
828 const unsigned int enriched_dof =
829 base_no_mult_local_enriched_dofs[base_no][m][i];
830 for (
unsigned int q = 0; q < n_q_points; ++q)
842 for (
unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)
844 for (
unsigned int m = 0;
845 m < base_no_mult_local_enriched_dofs[base_no].size();
847 for (
unsigned int i = 0;
848 i < base_no_mult_local_enriched_dofs[base_no][m].size();
851 const unsigned int enriched_dof =
852 base_no_mult_local_enriched_dofs[base_no][m][i];
853 for (
unsigned int q = 0; q < n_q_points; ++q)
863 template <
int dim,
int spacedim>
871 template <
int dim,
int spacedim>
879 template <
int dim,
int spacedim>
888 fe_system->get_face_interpolation_matrix(fe_enr_other->get_fe_system(),
896 spacedim>::ExcInterpolationNotImplemented()));
901 template <
int dim,
int spacedim>
905 const unsigned int subface,
911 fe_system->get_subface_interpolation_matrix(fe_enr_other->get_fe_system(),
920 spacedim>::ExcInterpolationNotImplemented()));
925 template <
int dim,
int spacedim>
926 std::vector<std::pair<unsigned int, unsigned int>>
933 return fe_system->hp_vertex_dof_identities(fe_enr_other->get_fe_system());
938 return std::vector<std::pair<unsigned int, unsigned int>>();
943 template <
int dim,
int spacedim>
944 std::vector<std::pair<unsigned int, unsigned int>>
951 return fe_system->hp_line_dof_identities(fe_enr_other->get_fe_system());
956 return std::vector<std::pair<unsigned int, unsigned int>>();
961 template <
int dim,
int spacedim>
962 std::vector<std::pair<unsigned int, unsigned int>>
969 return fe_system->hp_quad_dof_identities(fe_enr_other->get_fe_system());
974 return std::vector<std::pair<unsigned int, unsigned int>>();
979 template <
int dim,
int spacedim>
983 const unsigned int codim)
const
1002 return fe_system->compare_for_domination(fe_enr_other->get_fe_system(),
1013 template <
int dim,
int spacedim>
1016 const unsigned int child,
1019 return fe_system->get_prolongation_matrix(child, refinement_case);
1023 template <
int dim,
int spacedim>
1026 const unsigned int child,
1029 return fe_system->get_restriction_matrix(child, refinement_case);
1036 template <
int dim,
int spacedim>
1039 : fesystem_data(std::move(fesystem_data))
1043 template <
int dim,
int spacedim>
1046 const unsigned int base_no)
const
1048 return fesystem_data->get_fe_data(base_no);
1052 template <
int dim,
int spacedim>
1055 const unsigned int base_no)
const
1057 return fesystem_data->get_fe_output_object(base_no);
1065 template <
int dim,
int spacedim>
1073 std::vector<bool> vertices_subdomain_1(
1078 if (predicate_1(cell))
1080 vertices_subdomain_1[cell->vertex_index(v)] =
true;
1084 if (predicate_2(cell))
1086 if (vertices_subdomain_1[cell->vertex_index(v)] ==
true)
1095 template <
int dim,
int spacedim>
1100 std::vector<unsigned int> & predicate_colors)
1102 const unsigned int num_indices = predicates.size();
1107 dsp.
reinit(num_indices, num_indices);
1114 for (
unsigned int i = 0; i < num_indices; ++i)
1115 for (
unsigned int j = i + 1; j < num_indices; ++j)
1127 predicate_colors.resize(num_indices);
1135 template <
int dim,
int spacedim>
1140 const std::vector<unsigned int> & predicate_colors,
1141 std::map<
unsigned int, std::map<unsigned int, unsigned int>>
1142 & cellwise_color_predicate_map,
1143 std::vector<std::set<unsigned int>> &fe_sets)
1147 cellwise_color_predicate_map.clear();
1175 unsigned int map_index = 0;
1179 cell->set_active_fe_index(0);
1188 cell->set_material_id(map_index);
1189 std::set<unsigned int> color_list;
1197 std::map<unsigned int, unsigned int> &cell_map =
1198 cellwise_color_predicate_map[map_index];
1199 for (
unsigned int i = 0; i < predicates.size(); ++i)
1201 if (predicates[i](cell))
1207 auto ret = cell_map.insert(
1208 std::pair<unsigned int, unsigned int>(predicate_colors[i],
1213 "Only one enrichment function per color"));
1215 color_list.insert(predicate_colors[i]);
1232 if (!color_list.empty())
1235 std::find(fe_sets.begin(), fe_sets.end(), color_list);
1237 if (it == fe_sets.end())
1239 fe_sets.push_back(color_list);
1240 cell->set_active_fe_index(fe_sets.size() - 1);
1245 cell->set_active_fe_index(std::distance(fe_sets.begin(), it));
1292 const unsigned int fe_index = cell->active_fe_index();
1293 const std::set<unsigned int> fe_set = fe_sets.at(fe_index);
1300 if (!cell->at_boundary(face) &&
1301 cell->material_id() < cell->neighbor(face)->material_id())
1303 const auto nbr_fe_index =
1304 cell->neighbor(face)->active_fe_index();
1307 const auto nbr_fe_set = fe_sets.at(nbr_fe_index);
1310 std::set<unsigned int> intersection_set;
1311 std::set_intersection(
1316 std::inserter(intersection_set, intersection_set.begin()));
1319 if (!intersection_set.empty())
1321 const auto it = std::find(fe_sets.begin(),
1325 if (it == fe_sets.end())
1327 fe_sets.push_back(intersection_set);
1337 template <
int dim,
int spacedim>
1340 const unsigned int n_colors,
1342 const std::map<
unsigned int, std::map<unsigned int, unsigned int>>
1343 &cellwise_color_predicate_map,
1348 color_enrichments.clear();
1364 color_enrichments.resize(n_colors);
1365 for (
unsigned int i = 0; i < n_colors; ++i)
1367 color_enrichments[i] =
1370 const unsigned int id = cell->material_id();
1377 return enrichments[cellwise_color_predicate_map.at(
id).at(i + 1)]
1385 template <
int dim,
int spacedim>
1388 const unsigned int n_colors,
1389 const std::vector<std::set<unsigned int>> &fe_sets,
1392 & color_enrichments,
1399 const std::function<const Function<spacedim> *(
1405 ExcMessage(
"Called enrichment function for FE_Nothing"));
1413 for (
const auto &fe_set : fe_sets)
1415 std::vector<const FiniteElement<dim, spacedim> *> vec_fe_enriched(
1416 n_colors, &fe_nothing);
1417 std::vector<std::vector<std::function<const Function<spacedim> *(
1421 for (
const unsigned int color_id : fe_set)
1426 const unsigned int ind = color_id - 1;
1436 vec_fe_enriched[ind] = &fe_enriched;
1440 functions[ind][0] = color_enrichments[ind];
1455 template <
int dim,
int spacedim>
1462 , fe_enriched(fe_enriched)
1464 , predicates(predicates)
1471 ExcMessage(
"Number of predicates should be positive"));
1476 template <
int dim,
int spacedim>
1490 cellwise_color_predicate_map,
1494 internal::make_colorwise_enrichment_functions<dim, spacedim>(
1495 n_colors, enrichments, cellwise_color_predicate_map, color_enrichments);
1506 return fe_collection;
1512 #include "fe_enriched.inst"