21 #ifndef _aspect_compat_h 22 #define _aspect_compat_h 33 #if DEAL_II_VERSION_GTE(9,4,0) 53 const unsigned int root,
56 Assert(root < ::Utilities::MPI::n_mpi_processes(comm),
57 ::ExcMessage(
"Invalid root rank specified."));
61 const size_t max_send_count = std::numeric_limits<signed int>::max();
63 size_t total_sent_count = 0;
64 while (total_sent_count < count)
66 const size_t current_count =
67 std::min(count - total_sent_count, max_send_count);
69 const int ierr = MPI_Bcast(buffer + total_sent_count,
75 total_sent_count += current_count;
84 #if !DEAL_II_VERSION_GTE(9,5,0) 85 #include <deal.II/numerics/vector_tools.h> 86 #include <deal.II/hp/fe_values.h> 87 #include <deal.II/grid/manifold.h> 101 types::global_dof_index dof_indices[dim];
105 for (
unsigned int i = 0; i < dim; ++i)
106 dof_indices[i] = numbers::invalid_dof_index;
111 operator<(const VectorDoFTuple<dim> &other)
const 113 for (
unsigned int i = 0; i < dim; ++i)
114 if (dof_indices[i] < other.dof_indices[i])
116 else if (dof_indices[i] > other.dof_indices[i])
124 for (
unsigned int i = 0; i < dim; ++i)
134 return !(*
this == other);
141 operator<<(std::ostream &out, const VectorDoFTuple<dim> &vdt)
143 for (
unsigned int d = 0; d < dim; ++d)
144 out << vdt.dof_indices[d] << (d < dim - 1 ?
" " :
"");
163 const Tensor<1, dim> &constraining_vector,
164 AffineConstraints<double> &constraints,
165 const double inhomogeneity = 0)
205 if (std::fabs(constraining_vector[0]) >
206 std::fabs(constraining_vector[1]) + 1e-10)
208 if (!constraints.is_constrained(dof_indices.
dof_indices[0]) &&
209 constraints.can_store_line(dof_indices.
dof_indices[0]))
213 if (std::fabs(constraining_vector[1] /
214 constraining_vector[0]) >
215 std::numeric_limits<double>::epsilon())
218 -constraining_vector[1] /
219 constraining_vector[0]);
221 if (std::fabs(inhomogeneity / constraining_vector[0]) >
222 std::numeric_limits<double>::epsilon())
223 constraints.set_inhomogeneity(
225 inhomogeneity / constraining_vector[0]);
230 if (!constraints.is_constrained(dof_indices.
dof_indices[1]) &&
231 constraints.can_store_line(dof_indices.
dof_indices[1]))
235 if (std::fabs(constraining_vector[0] /
236 constraining_vector[1]) >
237 std::numeric_limits<double>::epsilon())
240 -constraining_vector[0] /
241 constraining_vector[1]);
243 if (std::fabs(inhomogeneity / constraining_vector[1]) >
244 std::numeric_limits<double>::epsilon())
245 constraints.set_inhomogeneity(
247 inhomogeneity / constraining_vector[1]);
255 if ((std::fabs(constraining_vector[0]) >=
256 std::fabs(constraining_vector[1]) + 1e-10) &&
257 (std::fabs(constraining_vector[0]) >=
258 std::fabs(constraining_vector[2]) + 2e-10))
260 if (!constraints.is_constrained(dof_indices.
dof_indices[0]) &&
261 constraints.can_store_line(dof_indices.
dof_indices[0]))
265 if (std::fabs(constraining_vector[1] /
266 constraining_vector[0]) >
267 std::numeric_limits<double>::epsilon())
270 -constraining_vector[1] /
271 constraining_vector[0]);
273 if (std::fabs(constraining_vector[2] /
274 constraining_vector[0]) >
275 std::numeric_limits<double>::epsilon())
278 -constraining_vector[2] /
279 constraining_vector[0]);
281 if (std::fabs(inhomogeneity / constraining_vector[0]) >
282 std::numeric_limits<double>::epsilon())
283 constraints.set_inhomogeneity(
285 inhomogeneity / constraining_vector[0]);
288 else if ((std::fabs(constraining_vector[1]) + 1e-10 >=
289 std::fabs(constraining_vector[0])) &&
290 (std::fabs(constraining_vector[1]) >=
291 std::fabs(constraining_vector[2]) + 1e-10))
293 if (!constraints.is_constrained(dof_indices.
dof_indices[1]) &&
294 constraints.can_store_line(dof_indices.
dof_indices[1]))
298 if (std::fabs(constraining_vector[0] /
299 constraining_vector[1]) >
300 std::numeric_limits<double>::epsilon())
303 -constraining_vector[0] /
304 constraining_vector[1]);
306 if (std::fabs(constraining_vector[2] /
307 constraining_vector[1]) >
308 std::numeric_limits<double>::epsilon())
311 -constraining_vector[2] /
312 constraining_vector[1]);
314 if (std::fabs(inhomogeneity / constraining_vector[1]) >
315 std::numeric_limits<double>::epsilon())
316 constraints.set_inhomogeneity(
318 inhomogeneity / constraining_vector[1]);
323 if (!constraints.is_constrained(dof_indices.
dof_indices[2]) &&
324 constraints.can_store_line(dof_indices.
dof_indices[2]))
328 if (std::fabs(constraining_vector[0] /
329 constraining_vector[2]) >
330 std::numeric_limits<double>::epsilon())
333 -constraining_vector[0] /
334 constraining_vector[2]);
336 if (std::fabs(constraining_vector[1] /
337 constraining_vector[2]) >
338 std::numeric_limits<double>::epsilon())
341 -constraining_vector[1] /
342 constraining_vector[2]);
344 if (std::fabs(inhomogeneity / constraining_vector[2]) >
345 std::numeric_limits<double>::epsilon())
346 constraints.set_inhomogeneity(
348 inhomogeneity / constraining_vector[2]);
356 Assert(
false, ExcNotImplemented());
377 const Tensor<1, dim> &tangent_vector,
378 AffineConstraints<double> &constraints,
379 const Vector<double> &b_values = Vector<double>(dim))
391 unsigned int largest_component = 0;
392 for (
unsigned int d = 1; d < dim; ++d)
393 if (std::fabs(tangent_vector[d]) >
394 std::fabs(tangent_vector[largest_component]) + 1e-10)
395 largest_component = d;
400 for (
unsigned int d = 0; d < dim; ++d)
401 if (d != largest_component)
402 if (!constraints.is_constrained(dof_indices.
dof_indices[d]) &&
403 constraints.can_store_line(dof_indices.
dof_indices[d]))
407 if (std::fabs(tangent_vector[d] /
408 tangent_vector[largest_component]) >
409 std::numeric_limits<double>::epsilon())
410 constraints.add_entry(
413 tangent_vector[d] / tangent_vector[largest_component]);
415 const double inhomogeneity =
416 (b_values(d) * tangent_vector[largest_component] -
417 b_values(largest_component) * tangent_vector[d]) /
418 tangent_vector[largest_component];
420 if (std::fabs(inhomogeneity) >
421 std::numeric_limits<double>::epsilon())
422 constraints.set_inhomogeneity(dof_indices.
dof_indices[d],
433 template <
int dim,
int spacedim>
436 const typename DoFHandler<dim, spacedim>::cell_iterator &cell,
437 const unsigned int first_vector_component,
438 const std::set<types::boundary_id> &boundary_ids,
441 hp::FEFaceValues<dim, spacedim> &x_fe_face_values,
442 const unsigned int n_dofs,
443 const IndexSet &refinement_edge_indices,
444 const unsigned int level,
447 std::pair<Tensor<1, dim>,
448 typename DoFHandler<dim, spacedim>::cell_iterator>>
451 &dof_vector_to_b_values)
453 std::set<types::boundary_id>::iterator b_id;
454 for (
const unsigned int face_no : cell->face_indices())
455 if ((b_id = boundary_ids.find(cell->face(face_no)->boundary_id())) !=
458 const FiniteElement<dim> &fe = cell->get_fe();
459 typename DoFHandler<dim, spacedim>::level_face_iterator face =
462 std::vector<types::global_dof_index> face_dofs;
464 face_dofs.resize(fe.n_dofs_per_face(face_no));
466 if (level != numbers::invalid_unsigned_int)
467 face->get_mg_dof_indices(level,
469 cell->active_fe_index());
471 face->get_dof_indices(face_dofs, cell->active_fe_index());
473 x_fe_face_values.reinit(cell, face_no);
474 const FEFaceValues<dim> &fe_values =
475 x_fe_face_values.get_present_fe_values();
477 const auto &face_quadrature_collection =
478 x_fe_face_values.get_quadrature_collection();
482 for (
unsigned int i = 0; i < face_dofs.size(); ++i)
483 if (fe.face_system_to_component_index(i, face_no).first ==
484 first_vector_component)
488 if (level == numbers::invalid_unsigned_int ||
489 !refinement_edge_indices.is_element(face_dofs[i]))
496 first_vector_component + dim <= fe.n_components(),
498 "Error: the finite element does not have enough components " 499 "to define a normal direction."));
501 for (
unsigned int k = 0; k < fe.n_dofs_per_face(face_no);
504 (face_quadrature_collection[cell->active_fe_index()]
506 face_quadrature_collection[cell->active_fe_index()]
508 (fe.face_system_to_component_index(k, face_no)
509 .first >= first_vector_component) &&
510 (fe.face_system_to_component_index(k, face_no).first <
511 first_vector_component + dim))
513 [fe.face_system_to_component_index(k, face_no).first -
514 first_vector_component] = face_dofs[k];
516 for (
unsigned int d = 0; d < dim; ++d)
569 Tensor<1, dim> normal_vector =
570 (cell->face(face_no)->get_manifold().normal_vector(
571 cell->face(face_no), fe_values.quadrature_point(i)));
572 if (normal_vector * fe_values.normal_vector(i) < 0)
574 Assert(std::fabs(normal_vector.norm() - 1) < 1e-14,
576 for (
unsigned int d = 0; d < dim; ++d)
577 if (std::fabs(normal_vector[d]) < 1e-13)
578 normal_vector[d] = 0;
579 normal_vector /= normal_vector.norm();
581 const Point<dim> &point = fe_values.quadrature_point(i);
582 Vector<double> b_values(dim);
583 function_map.at(*b_id)->vector_value(point, b_values);
587 dof_to_normals_map.insert(
588 std::make_pair(vector_dofs,
589 std::make_pair(normal_vector, cell)));
590 dof_vector_to_b_values.insert(
591 std::make_pair(vector_dofs, b_values));
593 #ifdef DEBUG_NO_NORMAL_FLUX 594 std::cout <<
"Adding normal vector:" << std::endl
595 <<
" dofs=" << vector_dofs << std::endl
596 <<
" cell=" << cell <<
" at " << cell->center()
598 <<
" normal=" << normal_vector << std::endl;
614 template <
int dim,
int spacedim>
617 const DoFHandler<dim, spacedim> &dof_handler,
618 const unsigned int first_vector_component,
619 const std::set<types::boundary_id> &boundary_ids,
622 AffineConstraints<double> &constraints,
623 const Mapping<dim, spacedim> &mapping,
624 const IndexSet &refinement_edge_indices = IndexSet(),
625 const unsigned int level = numbers::invalid_unsigned_int)
628 ExcMessage(
"This function is not useful in 1d because it amounts " 629 "to imposing Dirichlet values on the vector-valued " 634 const hp::FECollection<dim, spacedim> &fe_collection =
635 dof_handler.get_fe_collection();
636 hp::MappingCollection<dim, spacedim> mapping_collection;
637 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
638 mapping_collection.push_back(mapping);
642 AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
643 const unsigned int face_no = 0;
648 hp::QCollection<dim - 1> face_quadrature_collection;
649 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
651 const std::vector<Point<dim - 1>> &unit_support_points =
652 fe_collection[i].get_unit_face_support_points(face_no);
654 Assert(unit_support_points.size() ==
655 fe_collection[i].n_dofs_per_face(face_no),
658 face_quadrature_collection.push_back(
659 Quadrature<dim - 1>(unit_support_points));
663 hp::FEFaceValues<dim, spacedim> x_fe_face_values(
666 face_quadrature_collection,
667 update_quadrature_points | update_normal_vectors);
680 using DoFToNormalsMap = std::multimap<
682 std::pair<Tensor<1, dim>,
683 typename DoFHandler<dim, spacedim>::cell_iterator>>;
684 std::map<internal::VectorDoFTuple<dim>, Vector<double>>
685 dof_vector_to_b_values;
687 DoFToNormalsMap dof_to_normals_map;
689 const unsigned int n_dof = dof_handler.n_dofs();
691 if (level == numbers::invalid_unsigned_int)
694 for (
const auto &cell : dof_handler.active_cell_iterators())
695 if (!cell->is_artificial())
699 first_vector_component,
704 refinement_edge_indices,
707 dof_vector_to_b_values);
713 for (
const auto &cell : dof_handler.cell_iterators_on_level(level))
714 if (cell->level_subdomain_id() !=
715 numbers::artificial_subdomain_id &&
716 cell->level_subdomain_id() != numbers::invalid_subdomain_id)
720 first_vector_component,
725 refinement_edge_indices,
728 dof_vector_to_b_values);
737 typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
739 while (p != dof_to_normals_map.end())
744 typename DoFToNormalsMap::const_iterator same_dof_range[2] = {p};
745 for (++p; p != dof_to_normals_map.end(); ++p)
746 if (p->first != same_dof_range[0]->first)
748 same_dof_range[1] = p;
751 if (p == dof_to_normals_map.end())
752 same_dof_range[1] = dof_to_normals_map.end();
754 #ifdef DEBUG_NO_NORMAL_FLUX 755 std::cout <<
"For dof indices <" << p->first
756 <<
">, found the following normals" << std::endl;
757 for (
typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
758 q != same_dof_range[1];
760 std::cout <<
" " << q->second.first <<
" from cell " 761 << q->second.second << std::endl;
769 using CellToNormalsMap =
770 std::map<
typename DoFHandler<dim, spacedim>::cell_iterator,
771 std::pair<Tensor<1, dim>,
unsigned int>>;
773 CellToNormalsMap cell_to_normals_map;
774 for (
typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
775 q != same_dof_range[1];
777 if (cell_to_normals_map.find(q->second.second) ==
778 cell_to_normals_map.end())
779 cell_to_normals_map[q->second.second] =
780 std::make_pair(q->second.first, 1U);
783 const Tensor<1, dim> old_normal =
784 cell_to_normals_map[q->second.second].first;
785 const unsigned int old_count =
786 cell_to_normals_map[q->second.second].second;
788 Assert(old_count > 0, ExcInternalError());
792 cell_to_normals_map[q->second.second] =
793 std::make_pair((old_normal * old_count + q->second.first) /
797 Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
799 #ifdef DEBUG_NO_NORMAL_FLUX 800 std::cout <<
" cell_to_normals_map:" << std::endl;
801 for (
typename CellToNormalsMap::const_iterator x =
802 cell_to_normals_map.begin();
803 x != cell_to_normals_map.end();
805 std::cout <<
" " << x->first <<
" -> (" << x->second.first
806 <<
',' << x->second.second <<
')' << std::endl;
810 unsigned int max_n_contributions_per_cell = 1;
811 for (
typename CellToNormalsMap::const_iterator x =
812 cell_to_normals_map.begin();
813 x != cell_to_normals_map.end();
815 max_n_contributions_per_cell =
816 std::max(max_n_contributions_per_cell, x->second.second);
821 Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
823 switch (max_n_contributions_per_cell)
842 Tensor<1, dim> normal;
843 for (
typename CellToNormalsMap::const_iterator x =
844 cell_to_normals_map.begin();
845 x != cell_to_normals_map.end();
847 normal += x->second.first;
848 normal /= normal.norm();
851 for (
unsigned int d = 0; d < dim; ++d)
852 if (std::fabs(normal[d]) < 1e-13)
854 normal /= normal.norm();
857 const internal::VectorDoFTuple<dim> &dof_indices =
858 same_dof_range[0]->first;
859 double normal_value = 0.;
860 const Vector<double> b_values =
861 dof_vector_to_b_values[dof_indices];
862 for (
unsigned int i = 0; i < dim; ++i)
863 normal_value += b_values[i] * normal[i];
884 Assert(cell_to_normals_map.size() == 1, ExcInternalError());
896 typename DoFToNormalsMap::const_iterator x =
898 for (
unsigned int i = 0; i < dim; ++i, ++x)
899 for (
unsigned int j = 0; j < dim; ++j)
900 t[i][j] = x->second.first[j];
903 std::fabs(determinant(t)) > 1e-3,
905 "Found a set of normal vectors that are nearly collinear."));
912 const internal::VectorDoFTuple<dim> &dof_indices =
913 same_dof_range[0]->first;
914 const Vector<double> b_values =
915 dof_vector_to_b_values[dof_indices];
916 for (
unsigned int i = 0; i < dim; ++i)
917 if (!constraints.is_constrained(
918 same_dof_range[0]->first.dof_indices[i]) &&
919 constraints.can_store_line(
920 same_dof_range[0]->first.dof_indices[i]))
922 const types::global_dof_index line =
924 constraints.add_line(line);
925 if (std::fabs(b_values[i]) >
926 std::numeric_limits<double>::epsilon())
927 constraints.set_inhomogeneity(line, b_values[i]);
938 Assert(dim >= 3, ExcNotImplemented());
939 Assert(max_n_contributions_per_cell == 2, ExcInternalError());
946 using CellContributions =
947 std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
948 std::list<Tensor<1, dim>>>;
949 CellContributions cell_contributions;
951 for (
typename DoFToNormalsMap::const_iterator q =
953 q != same_dof_range[1];
955 cell_contributions[q->second.second].push_back(
957 Assert(cell_contributions.size() >= 1, ExcInternalError());
983 std::list<Tensor<1, dim>> tangential_vectors;
984 for (
typename CellContributions::const_iterator contribution =
985 cell_contributions.begin();
986 contribution != cell_contributions.end();
989 #ifdef DEBUG_NO_NORMAL_FLUX 991 <<
" Treating edge case with dim-1 contributions." 993 <<
" Looking at cell " << contribution->first
994 <<
" which has contributed these normal vectors:" 996 for (
typename std::list<Tensor<1, dim>>::const_iterator
997 t = contribution->second.begin();
998 t != contribution->second.end();
1000 std::cout <<
" " << *t << std::endl;
1005 if (contribution->second.size() < dim - 1)
1008 Tensor<1, dim> normals[dim - 1];
1010 unsigned int index = 0;
1011 for (
typename std::list<Tensor<1, dim>>::const_iterator
1012 t = contribution->second.begin();
1013 t != contribution->second.end();
1015 normals[index] = *t;
1016 Assert(index == dim - 1, ExcInternalError());
1027 Tensor<1, dim> tangent;
1039 cross_product_3d(normals[0], normals[dim - 2]);
1042 Assert(
false, ExcNotImplemented());
1046 std::fabs(tangent.norm()) > 1e-12,
1048 "Two normal vectors from adjacent faces are almost " 1050 tangent /= tangent.norm();
1052 tangential_vectors.push_back(tangent);
1060 const Tensor<1, dim> first_tangent =
1061 tangential_vectors.front();
1062 typename std::list<Tensor<1, dim>>::iterator t =
1063 tangential_vectors.begin();
1065 for (; t != tangential_vectors.end(); ++t)
1066 if (*t * first_tangent < 0)
1071 Tensor<1, dim> average_tangent;
1072 for (
typename std::list<Tensor<1, dim>>::const_iterator t =
1073 tangential_vectors.begin();
1074 t != tangential_vectors.end();
1076 average_tangent += *t;
1077 average_tangent /= average_tangent.norm();
1081 const internal::VectorDoFTuple<dim> &dof_indices =
1082 same_dof_range[0]->first;
1083 const Vector<double> b_values =
1084 dof_vector_to_b_values[dof_indices];
1096 template <
int dim,
int spacedim>
1099 const DoFHandler<dim, spacedim> &dof_handler,
1100 const unsigned int first_vector_component,
1101 const std::set<types::boundary_id> &boundary_ids,
1102 AffineConstraints<double> &constraints,
1103 const Mapping<dim, spacedim> &mapping,
1104 const IndexSet &refinement_edge_indices,
1105 const unsigned int level)
1107 Functions::ZeroFunction<dim> zero_function(dim);
1108 std::map<types::boundary_id, const Function<spacedim> *> function_map;
1109 for (
const types::boundary_id boundary_id : boundary_ids)
1110 function_map[boundary_id] = &zero_function;
1113 first_vector_component,
1118 refinement_edge_indices,
void broadcast(T *buffer, const size_t count, const unsigned int root, const MPI_Comm &comm)
MPI_Datatype mpi_type_id(const char *)