21 #ifndef _aspect_compat_h 22 #define _aspect_compat_h 33 using ::Utilities::MPI::broadcast;
40 #if !DEAL_II_VERSION_GTE(9,6,0) 41 #include <deal.II/multigrid/mg_transfer_matrix_free.h> 44 template<
int dim,
class NumberType>
51 #if !DEAL_II_VERSION_GTE(9,5,0) 52 #include <deal.II/numerics/vector_tools.h> 53 #include <deal.II/hp/fe_values.h> 54 #include <deal.II/grid/manifold.h> 72 for (
unsigned int i = 0; i < dim; ++i)
73 dof_indices[i] = numbers::invalid_dof_index;
78 operator<(const VectorDoFTuple<dim> &other)
const 80 return (dof_indices < other.dof_indices);
99 operator<<(std::ostream &out, const VectorDoFTuple<dim> &vdt)
101 for (
unsigned int d = 0; d < dim; ++d)
102 out << vdt.dof_indices[d] << (d < dim - 1 ?
" " :
"");
121 const Tensor<1, dim> &constraining_vector,
122 AffineConstraints<double> &constraints,
123 const double inhomogeneity = 0)
163 if (std::fabs(constraining_vector[0]) >
164 std::fabs(constraining_vector[1]) + 1e-10)
166 if (!constraints.is_constrained(dof_indices.
dof_indices[0]) &&
167 constraints.can_store_line(dof_indices.
dof_indices[0]))
171 if (std::fabs(constraining_vector[1] /
172 constraining_vector[0]) >
173 std::numeric_limits<double>::epsilon())
176 -constraining_vector[1] /
177 constraining_vector[0]);
179 if (std::fabs(inhomogeneity / constraining_vector[0]) >
180 std::numeric_limits<double>::epsilon())
181 constraints.set_inhomogeneity(
183 inhomogeneity / constraining_vector[0]);
188 if (!constraints.is_constrained(dof_indices.
dof_indices[1]) &&
189 constraints.can_store_line(dof_indices.
dof_indices[1]))
193 if (std::fabs(constraining_vector[0] /
194 constraining_vector[1]) >
195 std::numeric_limits<double>::epsilon())
198 -constraining_vector[0] /
199 constraining_vector[1]);
201 if (std::fabs(inhomogeneity / constraining_vector[1]) >
202 std::numeric_limits<double>::epsilon())
203 constraints.set_inhomogeneity(
205 inhomogeneity / constraining_vector[1]);
213 if ((std::fabs(constraining_vector[0]) >=
214 std::fabs(constraining_vector[1]) + 1e-10) &&
215 (std::fabs(constraining_vector[0]) >=
216 std::fabs(constraining_vector[2]) + 2e-10))
218 if (!constraints.is_constrained(dof_indices.
dof_indices[0]) &&
219 constraints.can_store_line(dof_indices.
dof_indices[0]))
223 if (std::fabs(constraining_vector[1] /
224 constraining_vector[0]) >
225 std::numeric_limits<double>::epsilon())
228 -constraining_vector[1] /
229 constraining_vector[0]);
231 if (std::fabs(constraining_vector[2] /
232 constraining_vector[0]) >
233 std::numeric_limits<double>::epsilon())
236 -constraining_vector[2] /
237 constraining_vector[0]);
239 if (std::fabs(inhomogeneity / constraining_vector[0]) >
240 std::numeric_limits<double>::epsilon())
241 constraints.set_inhomogeneity(
243 inhomogeneity / constraining_vector[0]);
246 else if ((std::fabs(constraining_vector[1]) + 1e-10 >=
247 std::fabs(constraining_vector[0])) &&
248 (std::fabs(constraining_vector[1]) >=
249 std::fabs(constraining_vector[2]) + 1e-10))
251 if (!constraints.is_constrained(dof_indices.
dof_indices[1]) &&
252 constraints.can_store_line(dof_indices.
dof_indices[1]))
256 if (std::fabs(constraining_vector[0] /
257 constraining_vector[1]) >
258 std::numeric_limits<double>::epsilon())
261 -constraining_vector[0] /
262 constraining_vector[1]);
264 if (std::fabs(constraining_vector[2] /
265 constraining_vector[1]) >
266 std::numeric_limits<double>::epsilon())
269 -constraining_vector[2] /
270 constraining_vector[1]);
272 if (std::fabs(inhomogeneity / constraining_vector[1]) >
273 std::numeric_limits<double>::epsilon())
274 constraints.set_inhomogeneity(
276 inhomogeneity / constraining_vector[1]);
281 if (!constraints.is_constrained(dof_indices.
dof_indices[2]) &&
282 constraints.can_store_line(dof_indices.
dof_indices[2]))
286 if (std::fabs(constraining_vector[0] /
287 constraining_vector[2]) >
288 std::numeric_limits<double>::epsilon())
291 -constraining_vector[0] /
292 constraining_vector[2]);
294 if (std::fabs(constraining_vector[1] /
295 constraining_vector[2]) >
296 std::numeric_limits<double>::epsilon())
299 -constraining_vector[1] /
300 constraining_vector[2]);
302 if (std::fabs(inhomogeneity / constraining_vector[2]) >
303 std::numeric_limits<double>::epsilon())
304 constraints.set_inhomogeneity(
306 inhomogeneity / constraining_vector[2]);
314 Assert(
false, ExcNotImplemented());
335 const Tensor<1, dim> &tangent_vector,
336 AffineConstraints<double> &constraints,
337 const Vector<double> &b_values = Vector<double>(dim))
349 unsigned int largest_component = 0;
350 for (
unsigned int d = 1; d < dim; ++d)
351 if (std::fabs(tangent_vector[d]) >
352 std::fabs(tangent_vector[largest_component]) + 1e-10)
353 largest_component = d;
358 for (
unsigned int d = 0; d < dim; ++d)
359 if (d != largest_component)
360 if (!constraints.is_constrained(dof_indices.
dof_indices[d]) &&
361 constraints.can_store_line(dof_indices.
dof_indices[d]))
365 if (std::fabs(tangent_vector[d] /
366 tangent_vector[largest_component]) >
367 std::numeric_limits<double>::epsilon())
368 constraints.add_entry(
371 tangent_vector[d] / tangent_vector[largest_component]);
373 const double inhomogeneity =
374 (b_values(d) * tangent_vector[largest_component] -
375 b_values(largest_component) * tangent_vector[d]) /
376 tangent_vector[largest_component];
378 if (std::fabs(inhomogeneity) >
379 std::numeric_limits<double>::epsilon())
380 constraints.set_inhomogeneity(dof_indices.
dof_indices[d],
391 template <
int dim,
int spacedim>
394 const typename DoFHandler<dim, spacedim>::cell_iterator &cell,
395 const unsigned int first_vector_component,
396 const std::set<types::boundary_id> &boundary_ids,
399 hp::FEFaceValues<dim, spacedim> &x_fe_face_values,
400 const unsigned int n_dofs,
401 const IndexSet &refinement_edge_indices,
402 const unsigned int level,
405 std::pair<Tensor<1, dim>,
406 typename DoFHandler<dim, spacedim>::cell_iterator>>
409 &dof_vector_to_b_values)
411 std::set<types::boundary_id>::iterator b_id;
412 for (
const unsigned int face_no : cell->face_indices())
413 if ((b_id = boundary_ids.find(cell->face(face_no)->boundary_id())) !=
416 const FiniteElement<dim> &fe = cell->get_fe();
417 typename DoFHandler<dim, spacedim>::level_face_iterator face =
420 std::vector<types::global_dof_index> face_dofs;
422 face_dofs.resize(fe.n_dofs_per_face(face_no));
424 if (level != numbers::invalid_unsigned_int)
425 face->get_mg_dof_indices(level,
427 cell->active_fe_index());
429 face->get_dof_indices(face_dofs, cell->active_fe_index());
431 x_fe_face_values.reinit(cell, face_no);
432 const FEFaceValues<dim> &fe_values =
433 x_fe_face_values.get_present_fe_values();
435 const auto &face_quadrature_collection =
436 x_fe_face_values.get_quadrature_collection();
440 for (
unsigned int i = 0; i < face_dofs.size(); ++i)
441 if (fe.face_system_to_component_index(i, face_no).first ==
442 first_vector_component)
446 if (level == numbers::invalid_unsigned_int ||
447 !refinement_edge_indices.is_element(face_dofs[i]))
454 first_vector_component + dim <= fe.n_components(),
456 "Error: the finite element does not have enough components " 457 "to define a normal direction."));
459 for (
unsigned int k = 0; k < fe.n_dofs_per_face(face_no);
462 (face_quadrature_collection[cell->active_fe_index()]
464 face_quadrature_collection[cell->active_fe_index()]
466 (fe.face_system_to_component_index(k, face_no)
467 .first >= first_vector_component) &&
468 (fe.face_system_to_component_index(k, face_no).first <
469 first_vector_component + dim))
471 [fe.face_system_to_component_index(k, face_no).first -
472 first_vector_component] = face_dofs[k];
474 for (
unsigned int d = 0; d < dim; ++d)
527 Tensor<1, dim> normal_vector =
528 (cell->face(face_no)->get_manifold().normal_vector(
529 cell->face(face_no), fe_values.quadrature_point(i)));
530 if (normal_vector * fe_values.normal_vector(i) < 0)
532 Assert(std::fabs(normal_vector.norm() - 1) < 1e-14,
534 for (
unsigned int d = 0; d < dim; ++d)
535 if (std::fabs(normal_vector[d]) < 1e-13)
536 normal_vector[d] = 0;
537 normal_vector /= normal_vector.norm();
539 const Point<dim> &point = fe_values.quadrature_point(i);
540 Vector<double> b_values(dim);
541 function_map.at(*b_id)->vector_value(point, b_values);
545 dof_to_normals_map.insert(
546 std::make_pair(vector_dofs,
547 std::make_pair(normal_vector, cell)));
548 dof_vector_to_b_values.insert(
549 std::make_pair(vector_dofs, b_values));
551 #ifdef DEBUG_NO_NORMAL_FLUX 552 std::cout <<
"Adding normal vector:" << std::endl
553 <<
" dofs=" << vector_dofs << std::endl
554 <<
" cell=" << cell <<
" at " << cell->center()
556 <<
" normal=" << normal_vector << std::endl;
572 template <
int dim,
int spacedim>
575 const DoFHandler<dim, spacedim> &dof_handler,
576 const unsigned int first_vector_component,
577 const std::set<types::boundary_id> &boundary_ids,
580 AffineConstraints<double> &constraints,
581 const Mapping<dim, spacedim> &mapping,
582 const IndexSet &refinement_edge_indices = IndexSet(),
583 const unsigned int level = numbers::invalid_unsigned_int)
586 ExcMessage(
"This function is not useful in 1d because it amounts " 587 "to imposing Dirichlet values on the vector-valued " 592 const hp::FECollection<dim, spacedim> &fe_collection =
593 dof_handler.get_fe_collection();
594 hp::MappingCollection<dim, spacedim> mapping_collection;
595 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
596 mapping_collection.push_back(mapping);
600 AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
601 const unsigned int face_no = 0;
606 hp::QCollection<dim - 1> face_quadrature_collection;
607 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
609 const std::vector<Point<dim - 1>> &unit_support_points =
610 fe_collection[i].get_unit_face_support_points(face_no);
612 Assert(unit_support_points.size() ==
613 fe_collection[i].n_dofs_per_face(face_no),
616 face_quadrature_collection.push_back(
617 Quadrature<dim - 1>(unit_support_points));
621 hp::FEFaceValues<dim, spacedim> x_fe_face_values(
624 face_quadrature_collection,
625 update_quadrature_points | update_normal_vectors);
638 using DoFToNormalsMap = std::multimap<
640 std::pair<Tensor<1, dim>,
641 typename DoFHandler<dim, spacedim>::cell_iterator>>;
642 std::map<internal::VectorDoFTuple<dim>, Vector<double>>
643 dof_vector_to_b_values;
645 DoFToNormalsMap dof_to_normals_map;
647 const unsigned int n_dof = dof_handler.n_dofs();
649 if (level == numbers::invalid_unsigned_int)
652 for (
const auto &cell : dof_handler.active_cell_iterators())
653 if (!cell->is_artificial())
657 first_vector_component,
662 refinement_edge_indices,
665 dof_vector_to_b_values);
671 for (
const auto &cell : dof_handler.cell_iterators_on_level(level))
672 if (cell->level_subdomain_id() !=
673 numbers::artificial_subdomain_id &&
674 cell->level_subdomain_id() != numbers::invalid_subdomain_id)
678 first_vector_component,
683 refinement_edge_indices,
686 dof_vector_to_b_values);
695 typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
697 while (p != dof_to_normals_map.end())
702 typename DoFToNormalsMap::const_iterator same_dof_range[2] = {p};
703 for (++p; p != dof_to_normals_map.end(); ++p)
704 if (p->first != same_dof_range[0]->first)
706 same_dof_range[1] = p;
709 if (p == dof_to_normals_map.end())
710 same_dof_range[1] = dof_to_normals_map.end();
712 #ifdef DEBUG_NO_NORMAL_FLUX 713 std::cout <<
"For dof indices <" << p->first
714 <<
">, found the following normals" << std::endl;
715 for (
typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
716 q != same_dof_range[1];
718 std::cout <<
" " << q->second.first <<
" from cell " 719 << q->second.second << std::endl;
727 using CellToNormalsMap =
728 std::map<
typename DoFHandler<dim, spacedim>::cell_iterator,
729 std::pair<Tensor<1, dim>,
unsigned int>>;
731 CellToNormalsMap cell_to_normals_map;
732 for (
typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
733 q != same_dof_range[1];
735 if (cell_to_normals_map.find(q->second.second) ==
736 cell_to_normals_map.end())
737 cell_to_normals_map[q->second.second] =
738 std::make_pair(q->second.first, 1U);
741 const Tensor<1, dim> old_normal =
742 cell_to_normals_map[q->second.second].first;
743 const unsigned int old_count =
744 cell_to_normals_map[q->second.second].second;
746 Assert(old_count > 0, ExcInternalError());
750 cell_to_normals_map[q->second.second] =
751 std::make_pair((old_normal * old_count + q->second.first) /
755 Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
757 #ifdef DEBUG_NO_NORMAL_FLUX 758 std::cout <<
" cell_to_normals_map:" << std::endl;
759 for (
typename CellToNormalsMap::const_iterator x =
760 cell_to_normals_map.begin();
761 x != cell_to_normals_map.end();
763 std::cout <<
" " << x->first <<
" -> (" << x->second.first
764 <<
',' << x->second.second <<
')' << std::endl;
768 unsigned int max_n_contributions_per_cell = 1;
769 for (
typename CellToNormalsMap::const_iterator x =
770 cell_to_normals_map.begin();
771 x != cell_to_normals_map.end();
773 max_n_contributions_per_cell =
774 std::max(max_n_contributions_per_cell, x->second.second);
779 Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
781 switch (max_n_contributions_per_cell)
800 Tensor<1, dim> normal;
801 for (
typename CellToNormalsMap::const_iterator x =
802 cell_to_normals_map.begin();
803 x != cell_to_normals_map.end();
805 normal += x->second.first;
806 normal /= normal.norm();
809 for (
unsigned int d = 0; d < dim; ++d)
810 if (std::fabs(normal[d]) < 1e-13)
812 normal /= normal.norm();
815 const internal::VectorDoFTuple<dim> &dof_indices =
816 same_dof_range[0]->first;
817 double normal_value = 0.;
818 const Vector<double> b_values =
819 dof_vector_to_b_values[dof_indices];
820 for (
unsigned int i = 0; i < dim; ++i)
821 normal_value += b_values[i] * normal[i];
842 Assert(cell_to_normals_map.size() == 1, ExcInternalError());
854 typename DoFToNormalsMap::const_iterator x =
856 for (
unsigned int i = 0; i < dim; ++i, ++x)
857 for (
unsigned int j = 0; j < dim; ++j)
858 t[i][j] = x->second.first[j];
861 std::fabs(determinant(t)) > 1e-3,
863 "Found a set of normal vectors that are nearly collinear."));
870 const internal::VectorDoFTuple<dim> &dof_indices =
871 same_dof_range[0]->first;
872 const Vector<double> b_values =
873 dof_vector_to_b_values[dof_indices];
874 for (
unsigned int i = 0; i < dim; ++i)
875 if (!constraints.is_constrained(
876 same_dof_range[0]->first.dof_indices[i]) &&
877 constraints.can_store_line(
878 same_dof_range[0]->first.dof_indices[i]))
880 const types::global_dof_index line =
882 constraints.add_line(line);
883 if (std::fabs(b_values[i]) >
884 std::numeric_limits<double>::epsilon())
885 constraints.set_inhomogeneity(line, b_values[i]);
896 Assert(dim >= 3, ExcNotImplemented());
897 Assert(max_n_contributions_per_cell == 2, ExcInternalError());
904 using CellContributions =
905 std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
906 std::list<Tensor<1, dim>>>;
907 CellContributions cell_contributions;
909 for (
typename DoFToNormalsMap::const_iterator q =
911 q != same_dof_range[1];
913 cell_contributions[q->second.second].push_back(
915 Assert(cell_contributions.size() >= 1, ExcInternalError());
941 std::list<Tensor<1, dim>> tangential_vectors;
942 for (
typename CellContributions::const_iterator contribution =
943 cell_contributions.begin();
944 contribution != cell_contributions.end();
947 #ifdef DEBUG_NO_NORMAL_FLUX 949 <<
" Treating edge case with dim-1 contributions." 951 <<
" Looking at cell " << contribution->first
952 <<
" which has contributed these normal vectors:" 954 for (
typename std::list<Tensor<1, dim>>::const_iterator
955 t = contribution->second.begin();
956 t != contribution->second.end();
958 std::cout <<
" " << *t << std::endl;
963 if (contribution->second.size() < dim - 1)
966 Tensor<1, dim> normals[dim - 1];
968 unsigned int index = 0;
969 for (
typename std::list<Tensor<1, dim>>::const_iterator
970 t = contribution->second.begin();
971 t != contribution->second.end();
974 Assert(index == dim - 1, ExcInternalError());
985 Tensor<1, dim> tangent;
997 cross_product_3d(normals[0], normals[dim - 2]);
1000 Assert(
false, ExcNotImplemented());
1004 std::fabs(tangent.norm()) > 1e-12,
1006 "Two normal vectors from adjacent faces are almost " 1008 tangent /= tangent.norm();
1010 tangential_vectors.push_back(tangent);
1018 const Tensor<1, dim> first_tangent =
1019 tangential_vectors.front();
1020 typename std::list<Tensor<1, dim>>::iterator t =
1021 tangential_vectors.begin();
1023 for (; t != tangential_vectors.end(); ++t)
1024 if (*t * first_tangent < 0)
1029 Tensor<1, dim> average_tangent;
1030 for (
typename std::list<Tensor<1, dim>>::const_iterator t =
1031 tangential_vectors.begin();
1032 t != tangential_vectors.end();
1034 average_tangent += *t;
1035 average_tangent /= average_tangent.norm();
1039 const internal::VectorDoFTuple<dim> &dof_indices =
1040 same_dof_range[0]->first;
1041 const Vector<double> b_values =
1042 dof_vector_to_b_values[dof_indices];
1054 template <
int dim,
int spacedim>
1057 const DoFHandler<dim, spacedim> &dof_handler,
1058 const unsigned int first_vector_component,
1059 const std::set<types::boundary_id> &boundary_ids,
1060 AffineConstraints<double> &constraints,
1061 const Mapping<dim, spacedim> &mapping,
1062 const IndexSet &refinement_edge_indices,
1063 const unsigned int level)
1065 Functions::ZeroFunction<dim> zero_function(dim);
1066 std::map<types::boundary_id, const Function<spacedim> *> function_map;
1067 for (
const types::boundary_id boundary_id : boundary_ids)
1068 function_map[boundary_id] = &zero_function;
1071 first_vector_component,
1076 refinement_edge_indices,
1089 #if !DEAL_II_VERSION_GTE(9,6,0) 1091 #include <deal.II/grid/manifold.h> 1092 #include <deal.II/grid/manifold_lib.h> 1106 template <
int dim,
int spacedim = dim>
1121 virtual std::unique_ptr<Manifold<dim, spacedim>>
1122 clone()
const override;
1131 virtual Point<spacedim>
1132 get_intermediate_point(
const Point<spacedim> &p1,
1133 const Point<spacedim> &p2,
1134 const double w)
const override;
1140 virtual Tensor<1, spacedim>
1141 get_tangent_vector(
const Point<spacedim> &x1,
1142 const Point<spacedim> &x2)
const override;
1147 virtual Tensor<1, spacedim>
1149 const typename Triangulation<dim, spacedim>::face_iterator &face,
1150 const Point<spacedim> &p)
const override;
1156 get_normals_at_vertices(
1157 const typename Triangulation<dim, spacedim>::face_iterator &face,
1158 typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
1176 get_new_points(
const ArrayView<
const Point<spacedim>> &surrounding_points,
1177 const Table<2, double> &weights,
1178 ArrayView<Point<spacedim>> new_points)
const override;
1184 virtual Point<spacedim>
1185 get_new_point(
const ArrayView<
const Point<spacedim>> &vertices,
1186 const ArrayView<const double> &weights)
const override;
1200 std::pair<double, Tensor<1, spacedim>>
1201 guess_new_point(
const ArrayView<
const Tensor<1, spacedim>> &directions,
1202 const ArrayView<const double> &distances,
1203 const ArrayView<const double> &weights)
const;
1223 do_get_new_points(
const ArrayView<
const Point<spacedim>> &surrounding_points,
1224 const ArrayView<const double> &weights,
1225 ArrayView<Point<spacedim>> new_points)
const;
1243 #include <deal.II/grid/manifold_lib.h> 1246 using ::SphericalManifold;
MGTransferMatrixFree< dim, NumberType > MGTransferMF
const PolarManifold< spacedim > polar_manifold
const Point< spacedim > center