ASPECT
compat.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2015 - 2023 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef _aspect_compat_h
22 #define _aspect_compat_h
23 
24 #include <aspect/global.h>
25 // C++11 related includes.
26 #include <array>
27 #include <functional>
28 #include <memory>
29 
30 namespace big_mpi
31 {
32 
33  using ::Utilities::MPI::broadcast;
34 
35 }
36 
37 // deal.II 9.6 introduces the new MGTransferMF class as a replacement
38 // for MGTransferMatrixFree. Instead of putting an ifdef in every place,
39 // do this in one central location:
40 #if !DEAL_II_VERSION_GTE(9,6,0)
41 #include <deal.II/multigrid/mg_transfer_matrix_free.h>
42 namespace dealii
43 {
44  template<int dim, class NumberType>
45  using MGTransferMF = MGTransferMatrixFree<dim,NumberType>;
46 }
47 #endif
48 
49 
50 // Implement VectorTools::compute_no_normal_flux_constraints_on_level
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>
55 namespace dealii
56 {
57  namespace VectorTools
58  {
59  namespace internal
60  {
65  template <int dim>
67  {
68  std::array<types::global_dof_index,dim> dof_indices;
69 
71  {
72  for (unsigned int i = 0; i < dim; ++i)
73  dof_indices[i] = numbers::invalid_dof_index;
74  }
75 
76 
77  bool
78  operator<(const VectorDoFTuple<dim> &other) const
79  {
80  return (dof_indices < other.dof_indices);
81  }
82 
83  bool
84  operator==(const VectorDoFTuple<dim> &other) const
85  {
86  return (dof_indices == other.dof_indices);
87  }
88 
89  bool
90  operator!=(const VectorDoFTuple<dim> &other) const
91  {
92  return (dof_indices != other.dof_indices);
93  }
94  };
95 
96 
97  template <int dim>
98  std::ostream &
99  operator<<(std::ostream &out, const VectorDoFTuple<dim> &vdt)
100  {
101  for (unsigned int d = 0; d < dim; ++d)
102  out << vdt.dof_indices[d] << (d < dim - 1 ? " " : "");
103  return out;
104  }
105 
106 
107 
118  template <int dim>
119  void
121  const Tensor<1, dim> &constraining_vector,
122  AffineConstraints<double> &constraints,
123  const double inhomogeneity = 0)
124  {
125  // choose the DoF that has the largest component in the
126  // constraining_vector as the one to be constrained as this makes the
127  // process stable in cases where the constraining_vector has the form
128  // n=(1,0) or n=(0,1)
129  //
130  // we get constraints of the form x0 = a_1*x1 + a2*x2 + ... if one of
131  // the weights is essentially zero then skip this part. the
132  // AffineConstraints can also deal with cases like x0 = 0 if
133  // necessary
134  //
135  // there is a problem if we have a normal vector of the form
136  // (a,a,small) or (a,a,a). Depending on round-off we may choose the
137  // first or second component (or third, in the latter case) as the
138  // largest one, and depending on our choice one or another degree of
139  // freedom will be constrained. On a single processor this is not
140  // much of a problem, but it's a nightmare when we run in parallel
141  // and two processors disagree on which DoF should be constrained.
142  // This led to an incredibly difficult to find bug in step-32 when
143  // running in parallel with 9 or more processors.
144  //
145  // in practice, such normal vectors of the form (a,a,small) or
146  // (a,a,a) happen not infrequently since they lie on the diagonals
147  // where vertices frequently happen to land upon mesh refinement if
148  // one starts from a symmetric and regular body. we work around this
149  // problem in the following way: if we have a normal vector of the
150  // form (a,b) (similarly algorithm in 3d), we choose 'a' as the
151  // largest coefficient not if a>b but if a>b+1e-10. this shifts the
152  // problem away from the frequently visited diagonal to a line that
153  // is off the diagonal. there will of course be problems where the
154  // exact values of a and b differ by exactly 1e-10 and we get into
155  // the same instability, but from a practical viewpoint such problems
156  // should be much rarer. in particular, meshes have to be very fine
157  // for a vertex to land on this line if the original body had a
158  // vertex on the diagonal as well
159  switch (dim)
160  {
161  case 2:
162  {
163  if (std::fabs(constraining_vector[0]) >
164  std::fabs(constraining_vector[1]) + 1e-10)
165  {
166  if (!constraints.is_constrained(dof_indices.dof_indices[0]) &&
167  constraints.can_store_line(dof_indices.dof_indices[0]))
168  {
169  constraints.add_line(dof_indices.dof_indices[0]);
170 
171  if (std::fabs(constraining_vector[1] /
172  constraining_vector[0]) >
173  std::numeric_limits<double>::epsilon())
174  constraints.add_entry(dof_indices.dof_indices[0],
175  dof_indices.dof_indices[1],
176  -constraining_vector[1] /
177  constraining_vector[0]);
178 
179  if (std::fabs(inhomogeneity / constraining_vector[0]) >
180  std::numeric_limits<double>::epsilon())
181  constraints.set_inhomogeneity(
182  dof_indices.dof_indices[0],
183  inhomogeneity / constraining_vector[0]);
184  }
185  }
186  else
187  {
188  if (!constraints.is_constrained(dof_indices.dof_indices[1]) &&
189  constraints.can_store_line(dof_indices.dof_indices[1]))
190  {
191  constraints.add_line(dof_indices.dof_indices[1]);
192 
193  if (std::fabs(constraining_vector[0] /
194  constraining_vector[1]) >
195  std::numeric_limits<double>::epsilon())
196  constraints.add_entry(dof_indices.dof_indices[1],
197  dof_indices.dof_indices[0],
198  -constraining_vector[0] /
199  constraining_vector[1]);
200 
201  if (std::fabs(inhomogeneity / constraining_vector[1]) >
202  std::numeric_limits<double>::epsilon())
203  constraints.set_inhomogeneity(
204  dof_indices.dof_indices[1],
205  inhomogeneity / constraining_vector[1]);
206  }
207  }
208  break;
209  }
210 
211  case 3:
212  {
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))
217  {
218  if (!constraints.is_constrained(dof_indices.dof_indices[0]) &&
219  constraints.can_store_line(dof_indices.dof_indices[0]))
220  {
221  constraints.add_line(dof_indices.dof_indices[0]);
222 
223  if (std::fabs(constraining_vector[1] /
224  constraining_vector[0]) >
225  std::numeric_limits<double>::epsilon())
226  constraints.add_entry(dof_indices.dof_indices[0],
227  dof_indices.dof_indices[1],
228  -constraining_vector[1] /
229  constraining_vector[0]);
230 
231  if (std::fabs(constraining_vector[2] /
232  constraining_vector[0]) >
233  std::numeric_limits<double>::epsilon())
234  constraints.add_entry(dof_indices.dof_indices[0],
235  dof_indices.dof_indices[2],
236  -constraining_vector[2] /
237  constraining_vector[0]);
238 
239  if (std::fabs(inhomogeneity / constraining_vector[0]) >
240  std::numeric_limits<double>::epsilon())
241  constraints.set_inhomogeneity(
242  dof_indices.dof_indices[0],
243  inhomogeneity / constraining_vector[0]);
244  }
245  }
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))
250  {
251  if (!constraints.is_constrained(dof_indices.dof_indices[1]) &&
252  constraints.can_store_line(dof_indices.dof_indices[1]))
253  {
254  constraints.add_line(dof_indices.dof_indices[1]);
255 
256  if (std::fabs(constraining_vector[0] /
257  constraining_vector[1]) >
258  std::numeric_limits<double>::epsilon())
259  constraints.add_entry(dof_indices.dof_indices[1],
260  dof_indices.dof_indices[0],
261  -constraining_vector[0] /
262  constraining_vector[1]);
263 
264  if (std::fabs(constraining_vector[2] /
265  constraining_vector[1]) >
266  std::numeric_limits<double>::epsilon())
267  constraints.add_entry(dof_indices.dof_indices[1],
268  dof_indices.dof_indices[2],
269  -constraining_vector[2] /
270  constraining_vector[1]);
271 
272  if (std::fabs(inhomogeneity / constraining_vector[1]) >
273  std::numeric_limits<double>::epsilon())
274  constraints.set_inhomogeneity(
275  dof_indices.dof_indices[1],
276  inhomogeneity / constraining_vector[1]);
277  }
278  }
279  else
280  {
281  if (!constraints.is_constrained(dof_indices.dof_indices[2]) &&
282  constraints.can_store_line(dof_indices.dof_indices[2]))
283  {
284  constraints.add_line(dof_indices.dof_indices[2]);
285 
286  if (std::fabs(constraining_vector[0] /
287  constraining_vector[2]) >
288  std::numeric_limits<double>::epsilon())
289  constraints.add_entry(dof_indices.dof_indices[2],
290  dof_indices.dof_indices[0],
291  -constraining_vector[0] /
292  constraining_vector[2]);
293 
294  if (std::fabs(constraining_vector[1] /
295  constraining_vector[2]) >
296  std::numeric_limits<double>::epsilon())
297  constraints.add_entry(dof_indices.dof_indices[2],
298  dof_indices.dof_indices[1],
299  -constraining_vector[1] /
300  constraining_vector[2]);
301 
302  if (std::fabs(inhomogeneity / constraining_vector[2]) >
303  std::numeric_limits<double>::epsilon())
304  constraints.set_inhomogeneity(
305  dof_indices.dof_indices[2],
306  inhomogeneity / constraining_vector[2]);
307  }
308  }
309 
310  break;
311  }
312 
313  default:
314  Assert(false, ExcNotImplemented());
315  }
316  }
317 
318 
319 
331  template <int dim>
332  void
334  const VectorDoFTuple<dim> &dof_indices,
335  const Tensor<1, dim> &tangent_vector,
336  AffineConstraints<double> &constraints,
337  const Vector<double> &b_values = Vector<double>(dim))
338  {
339  // choose the DoF that has the
340  // largest component in the
341  // tangent_vector as the
342  // independent component, and
343  // then constrain the others to
344  // it. specifically, if, say,
345  // component 0 of the tangent
346  // vector t is largest by
347  // magnitude, then
348  // x1=(b[1]*t[0]-b[0]*t[1])/t[0]+t[1]/t[0]*x_0, etc.
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;
354 
355  // then constrain all of the
356  // other degrees of freedom in
357  // terms of the one just found
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]))
362  {
363  constraints.add_line(dof_indices.dof_indices[d]);
364 
365  if (std::fabs(tangent_vector[d] /
366  tangent_vector[largest_component]) >
367  std::numeric_limits<double>::epsilon())
368  constraints.add_entry(
369  dof_indices.dof_indices[d],
370  dof_indices.dof_indices[largest_component],
371  tangent_vector[d] / tangent_vector[largest_component]);
372 
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];
377 
378  if (std::fabs(inhomogeneity) >
379  std::numeric_limits<double>::epsilon())
380  constraints.set_inhomogeneity(dof_indices.dof_indices[d],
381  inhomogeneity);
382  }
383  }
384 
385 
386 
391  template <int dim, int spacedim>
392  void
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,
397  const std::map<types::boundary_id, const Function<spacedim> *>
398  &function_map,
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,
403  std::multimap<
405  std::pair<Tensor<1, dim>,
406  typename DoFHandler<dim, spacedim>::cell_iterator>>
407  &dof_to_normals_map,
408  std::map<internal::VectorDoFTuple<dim>, Vector<double>>
409  &dof_vector_to_b_values)
410  {
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())) !=
414  boundary_ids.end())
415  {
416  const FiniteElement<dim> &fe = cell->get_fe();
417  typename DoFHandler<dim, spacedim>::level_face_iterator face =
418  cell->face(face_no);
419 
420  std::vector<types::global_dof_index> face_dofs;
421  // get the indices of the dofs on this cell...
422  face_dofs.resize(fe.n_dofs_per_face(face_no));
423 
424  if (level != numbers::invalid_unsigned_int)
425  face->get_mg_dof_indices(level,
426  face_dofs,
427  cell->active_fe_index());
428  else
429  face->get_dof_indices(face_dofs, cell->active_fe_index());
430 
431  x_fe_face_values.reinit(cell, face_no);
432  const FEFaceValues<dim> &fe_values =
433  x_fe_face_values.get_present_fe_values();
434 
435  const auto &face_quadrature_collection =
436  x_fe_face_values.get_quadrature_collection();
437 
438  // then identify which of them correspond to the selected set of
439  // vector components
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)
443  // Refinement edge indices are going to be constrained to 0
444  // during a multigrid cycle and do not need no-normal-flux
445  // constraints, so skip them:
446  if (level == numbers::invalid_unsigned_int ||
447  !refinement_edge_indices.is_element(face_dofs[i]))
448  {
449  // find corresponding other components of vector
450  internal::VectorDoFTuple<dim> vector_dofs;
451  vector_dofs.dof_indices[0] = face_dofs[i];
452 
453  Assert(
454  first_vector_component + dim <= fe.n_components(),
455  ExcMessage(
456  "Error: the finite element does not have enough components "
457  "to define a normal direction."));
458 
459  for (unsigned int k = 0; k < fe.n_dofs_per_face(face_no);
460  ++k)
461  if ((k != i) &&
462  (face_quadrature_collection[cell->active_fe_index()]
463  .point(k) ==
464  face_quadrature_collection[cell->active_fe_index()]
465  .point(i)) &&
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))
470  vector_dofs.dof_indices
471  [fe.face_system_to_component_index(k, face_no).first -
472  first_vector_component] = face_dofs[k];
473 
474  for (unsigned int d = 0; d < dim; ++d)
475  Assert(vector_dofs.dof_indices[d] < n_dofs,
476  ExcInternalError());
477  (void)n_dofs;
478 
479  // we need the normal vector on this face. we know that
480  // it is a vector of length 1 but at least with higher
481  // order mappings it isn't always possible to guarantee
482  // that each component is exact up to zero tolerance. in
483  // particular, as shown in the deal.II/no_flux_06 test,
484  // if we just take the normal vector as given by the
485  // fe_values object, we can get entries in the normal
486  // vectors of the unit cube that have entries up to
487  // several times 1e-14.
488  //
489  // the problem with this is that this later yields
490  // constraints that are circular (e.g., in the testcase,
491  // we get constraints of the form
492  //
493  // x22 = 2.93099e-14*x21 + 2.93099e-14*x23
494  // x21 = -2.93099e-14*x22 + 2.93099e-14*x21
495  //
496  // in both of these constraints, the small numbers
497  // should be zero and the constraints should simply be
498  // x22 = x21 = 0
499  //
500  // to achieve this, we utilize that we know that the
501  // normal vector has (or should have) length 1 and that
502  // we can simply set small elements to zero (without
503  // having to check that they are small *relative to
504  // something else*). we do this and then normalize the
505  // length of the vector back to one, just to be on the
506  // safe side
507  //
508  // one more point: we would like to use the "real"
509  // normal vector here, as provided by the boundary
510  // description and as opposed to what we get from the
511  // FEValues object. we do this in the immediately next
512  // line, but as is obvious, the boundary only has a
513  // vague idea which side of a cell it is on -- indicated
514  // by the face number. in other words, it may provide
515  // the inner or outer normal. by and large, there is no
516  // harm from this, since the tangential vector we
517  // compute is still the same. however, we do average
518  // over normal vectors from adjacent cells and if they
519  // have recorded normal vectors from the inside once and
520  // from the outside the other time, then this averaging
521  // is going to run into trouble. as a consequence we ask
522  // the mapping after all for its normal vector, but we
523  // only ask it so that we can possibly correct the sign
524  // of the normal vector provided by the boundary if they
525  // should point in different directions. this is the
526  // case in tests/deal.II/no_flux_11.
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)
531  normal_vector *= -1;
532  Assert(std::fabs(normal_vector.norm() - 1) < 1e-14,
533  ExcInternalError());
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();
538 
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);
542 
543  // now enter the (dofs,(normal_vector,cell)) entry into
544  // the map
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));
550 
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()
555  << std::endl
556  << " normal=" << normal_vector << std::endl;
557 #endif
558  }
559  }
560  }
561 
562 
563 
572  template <int dim, int spacedim>
573  void
575  const DoFHandler<dim, spacedim> &dof_handler,
576  const unsigned int first_vector_component,
577  const std::set<types::boundary_id> &boundary_ids,
578  const std::map<types::boundary_id, const Function<spacedim> *>
579  &function_map,
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)
584  {
585  Assert(dim > 1,
586  ExcMessage("This function is not useful in 1d because it amounts "
587  "to imposing Dirichlet values on the vector-valued "
588  "quantity."));
589 
590  // create FE and mapping collections for all elements in use by this
591  // DoFHandler
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);
597 
598  // TODO: the implementation makes the assumption that all faces have the
599  // same number of dofs
600  AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
601  const unsigned int face_no = 0;
602 
603  // now also create a quadrature collection for the faces of a cell. fill
604  // it with a quadrature formula with the support points on faces for each
605  // FE
606  hp::QCollection<dim - 1> face_quadrature_collection;
607  for (unsigned int i = 0; i < fe_collection.size(); ++i)
608  {
609  const std::vector<Point<dim - 1>> &unit_support_points =
610  fe_collection[i].get_unit_face_support_points(face_no);
611 
612  Assert(unit_support_points.size() ==
613  fe_collection[i].n_dofs_per_face(face_no),
614  ExcInternalError());
615 
616  face_quadrature_collection.push_back(
617  Quadrature<dim - 1>(unit_support_points));
618  }
619 
620  // now create the object with which we will generate the normal vectors
621  hp::FEFaceValues<dim, spacedim> x_fe_face_values(
622  mapping_collection,
623  fe_collection,
624  face_quadrature_collection,
625  update_quadrature_points | update_normal_vectors);
626 
627  // have a map that stores normal vectors for each vector-dof tuple we want
628  // to constrain. since we can get at the same vector dof tuple more than
629  // once (for example if it is located at a vertex that we visit from all
630  // adjacent cells), we will want to average later on the normal vectors
631  // computed on different cells as described in the documentation of this
632  // function. however, we can only average if the contributions came from
633  // different cells, whereas we want to constrain twice or more in case the
634  // contributions came from different faces of the same cell
635  // (i.e. constrain not just the *average normal direction* but *all normal
636  // directions* we find). consequently, we also have to store which cell a
637  // normal vector was computed on
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;
644 
645  DoFToNormalsMap dof_to_normals_map;
646 
647  const unsigned int n_dof = dof_handler.n_dofs();
648 
649  if (level == numbers::invalid_unsigned_int)
650  {
651  // active cells
652  for (const auto &cell : dof_handler.active_cell_iterators())
653  if (!cell->is_artificial())
654  {
656  cell,
657  first_vector_component,
658  boundary_ids,
659  function_map,
660  x_fe_face_values,
661  n_dof,
662  refinement_edge_indices,
663  level,
664  dof_to_normals_map,
665  dof_vector_to_b_values);
666  }
667  }
668  else
669  {
670  // level cells
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)
675  {
677  cell,
678  first_vector_component,
679  boundary_ids,
680  function_map,
681  x_fe_face_values,
682  n_dof,
683  refinement_edge_indices,
684  level,
685  dof_to_normals_map,
686  dof_vector_to_b_values);
687  }
688  }
689 
690  // Now do something with the collected information. To this end, loop
691  // through all sets of pairs (dofs,normal_vector) and identify which
692  // entries belong to the same set of dofs and then do as described in the
693  // documentation, i.e. either average the normal vector or don't for this
694  // particular set of dofs
695  typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
696 
697  while (p != dof_to_normals_map.end())
698  {
699  // first find the range of entries in the multimap that corresponds to
700  // the same vector-dof tuple. as usual, we define the range
701  // half-open. the first entry of course is 'p'
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)
705  {
706  same_dof_range[1] = p;
707  break;
708  }
709  if (p == dof_to_normals_map.end())
710  same_dof_range[1] = dof_to_normals_map.end();
711 
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];
717  ++q)
718  std::cout << " " << q->second.first << " from cell "
719  << q->second.second << std::endl;
720 #endif
721 
722 
723  // now compute the reverse mapping: for each of the cells that
724  // contributed to the current set of vector dofs, add up the normal
725  // vectors. the values of the map are pairs of normal vectors and
726  // number of cells that have contributed
727  using CellToNormalsMap =
728  std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
729  std::pair<Tensor<1, dim>, unsigned int>>;
730 
731  CellToNormalsMap cell_to_normals_map;
732  for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
733  q != same_dof_range[1];
734  ++q)
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);
739  else
740  {
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;
745 
746  Assert(old_count > 0, ExcInternalError());
747 
748  // in the same entry, store again the now averaged normal vector
749  // and the new count
750  cell_to_normals_map[q->second.second] =
751  std::make_pair((old_normal * old_count + q->second.first) /
752  (old_count + 1),
753  old_count + 1);
754  }
755  Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
756 
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();
762  ++x)
763  std::cout << " " << x->first << " -> (" << x->second.first
764  << ',' << x->second.second << ')' << std::endl;
765 #endif
766 
767  // count the maximum number of contributions from each cell
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();
772  ++x)
773  max_n_contributions_per_cell =
774  std::max(max_n_contributions_per_cell, x->second.second);
775 
776  // verify that each cell can have only contributed at most dim times,
777  // since that is the maximum number of faces that come together at a
778  // single place
779  Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
780 
781  switch (max_n_contributions_per_cell)
782  {
783  // first deal with the case that a number of cells all have
784  // registered that they have a normal vector defined at the
785  // location of a given vector dof, and that each of them have
786  // encountered this vector dof exactly once while looping over all
787  // their faces. as stated in the documentation, this is the case
788  // where we want to simply average over all normal vectors
789  //
790  // the typical case is in 2d where multiple cells meet at one
791  // vertex sitting on the boundary. same in 3d for a vertex that
792  // is associated with only one of the boundary indicators passed
793  // to this function
794  case 1:
795  {
796  // compute the average normal vector from all the ones that
797  // have the same set of dofs. we could add them up and divide
798  // them by the number of additions, or simply normalize them
799  // right away since we want them to have unit length anyway
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();
804  ++x)
805  normal += x->second.first;
806  normal /= normal.norm();
807 
808  // normalize again
809  for (unsigned int d = 0; d < dim; ++d)
810  if (std::fabs(normal[d]) < 1e-13)
811  normal[d] = 0;
812  normal /= normal.norm();
813 
814  // then construct constraints from this:
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];
822  internal::add_constraint(dof_indices,
823  normal,
824  constraints,
825  normal_value);
826 
827  break;
828  }
829 
830  // this is the slightly more complicated case that a single cell
831  // has contributed with exactly DIM normal vectors to the same set
832  // of vector dofs. this is what happens in a corner in 2d and 3d
833  // (but not on an edge in 3d, where we have only 2, i.e. <DIM,
834  // contributions. Here we do not want to average the normal
835  // vectors. Since we have DIM contributions, let's assume (and
836  // verify) that they are in fact all linearly independent; in that
837  // case, all vector components are constrained and we need to set
838  // all of them to the corresponding boundary values
839  case dim:
840  {
841  // assert that indeed only a single cell has contributed
842  Assert(cell_to_normals_map.size() == 1, ExcInternalError());
843 
844  // check linear independence by computing the determinant of
845  // the matrix created from all the normal vectors. if they are
846  // linearly independent, then the determinant is nonzero. if
847  // they are orthogonal, then the matrix is in fact equal to 1
848  // (since they are all unit vectors); make sure the
849  // determinant is larger than 1e-3 to avoid cases where cells
850  // are degenerate
851  {
852  Tensor<2, dim> t;
853 
854  typename DoFToNormalsMap::const_iterator x =
855  same_dof_range[0];
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];
859 
860  Assert(
861  std::fabs(determinant(t)) > 1e-3,
862  ExcMessage(
863  "Found a set of normal vectors that are nearly collinear."));
864  }
865 
866  // so all components of this vector dof are constrained. enter
867  // this into the AffineConstraints object
868  //
869  // ignore dofs already constrained
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]))
879  {
880  const types::global_dof_index line =
881  dof_indices.dof_indices[i];
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]);
886  // no add_entries here
887  }
888 
889  break;
890  }
891 
892  // this is the case of an edge contribution in 3d, i.e. the vector
893  // is constrained in two directions but not the third.
894  default:
895  {
896  Assert(dim >= 3, ExcNotImplemented());
897  Assert(max_n_contributions_per_cell == 2, ExcInternalError());
898 
899  // as described in the documentation, let us first collect
900  // what each of the cells contributed at the current point. we
901  // use a std::list instead of a std::set (which would be more
902  // natural) because std::set requires that the stored elements
903  // are comparable with operator<
904  using CellContributions =
905  std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
906  std::list<Tensor<1, dim>>>;
907  CellContributions cell_contributions;
908 
909  for (typename DoFToNormalsMap::const_iterator q =
910  same_dof_range[0];
911  q != same_dof_range[1];
912  ++q)
913  cell_contributions[q->second.second].push_back(
914  q->second.first);
915  Assert(cell_contributions.size() >= 1, ExcInternalError());
916 
917  // now for each cell that has contributed determine the number
918  // of normal vectors it has contributed. we currently only
919  // implement if this is dim-1 for all cells (if a single cell
920  // has contributed dim, or if all adjacent cells have
921  // contributed 1 normal vector, this is already handled
922  // above).
923  //
924  // we only implement the case that all cells contribute
925  // dim-1 because we assume that we are following an edge
926  // of the domain (think: we are looking at a vertex
927  // located on one of the edges of a refined cube where the
928  // boundary indicators of the two adjacent faces of the
929  // cube are both listed in the set of boundary indicators
930  // passed to this function). in that case, all cells along
931  // that edge of the domain are assumed to have contributed
932  // dim-1 normal vectors. however, there are cases where
933  // this assumption is not justified (see the lengthy
934  // explanation in test no_flux_12.cc) and in those cases
935  // we simply ignore the cell that contributes only
936  // once. this is also discussed at length in the
937  // documentation of this function.
938  //
939  // for each contributing cell compute the tangential vector
940  // that remains unconstrained
941  std::list<Tensor<1, dim>> tangential_vectors;
942  for (typename CellContributions::const_iterator contribution =
943  cell_contributions.begin();
944  contribution != cell_contributions.end();
945  ++contribution)
946  {
947 #ifdef DEBUG_NO_NORMAL_FLUX
948  std::cout
949  << " Treating edge case with dim-1 contributions."
950  << std::endl
951  << " Looking at cell " << contribution->first
952  << " which has contributed these normal vectors:"
953  << std::endl;
954  for (typename std::list<Tensor<1, dim>>::const_iterator
955  t = contribution->second.begin();
956  t != contribution->second.end();
957  ++t)
958  std::cout << " " << *t << std::endl;
959 #endif
960 
961  // as mentioned above, simply ignore cells that only
962  // contribute once
963  if (contribution->second.size() < dim - 1)
964  continue;
965 
966  Tensor<1, dim> normals[dim - 1];
967  {
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();
972  ++t, ++index)
973  normals[index] = *t;
974  Assert(index == dim - 1, ExcInternalError());
975  }
976 
977  // calculate the tangent as the outer product of the
978  // normal vectors. since these vectors do not need to be
979  // orthogonal (think, for example, the case of the
980  // deal.II/no_flux_07 test: a sheared cube in 3d, with Q2
981  // elements, where we have constraints from the two normal
982  // vectors of two faces of the sheared cube that are not
983  // perpendicular to each other), we have to normalize the
984  // outer product
985  Tensor<1, dim> tangent;
986  switch (dim)
987  {
988  case 3:
989  // take cross product between normals[0] and
990  // normals[1]. write it in the current form (with
991  // [dim-2]) to make sure that compilers don't warn
992  // about out-of-bounds accesses -- the warnings are
993  // bogus since we get here only for dim==3, but at
994  // least one isn't quite smart enough to notice this
995  // and warns when compiling the function in 2d
996  tangent =
997  cross_product_3d(normals[0], normals[dim - 2]);
998  break;
999  default:
1000  Assert(false, ExcNotImplemented());
1001  }
1002 
1003  Assert(
1004  std::fabs(tangent.norm()) > 1e-12,
1005  ExcMessage(
1006  "Two normal vectors from adjacent faces are almost "
1007  "parallel."));
1008  tangent /= tangent.norm();
1009 
1010  tangential_vectors.push_back(tangent);
1011  }
1012 
1013  // go through the list of tangents and make sure that they all
1014  // roughly point in the same direction as the first one (i.e.
1015  // have an angle less than 90 degrees); if they don't then
1016  // flip their sign
1017  {
1018  const Tensor<1, dim> first_tangent =
1019  tangential_vectors.front();
1020  typename std::list<Tensor<1, dim>>::iterator t =
1021  tangential_vectors.begin();
1022  ++t;
1023  for (; t != tangential_vectors.end(); ++t)
1024  if (*t * first_tangent < 0)
1025  *t *= -1;
1026  }
1027 
1028  // now compute the average tangent and normalize it
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();
1033  ++t)
1034  average_tangent += *t;
1035  average_tangent /= average_tangent.norm();
1036 
1037  // now all that is left is that we add the constraints that
1038  // the vector is parallel to the tangent
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];
1044  average_tangent,
1045  constraints,
1046  b_values);
1047  }
1048  }
1049  }
1050  }
1051 
1052  } // namespece internal
1053 
1054  template <int dim, int spacedim>
1055  void
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)
1064  {
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;
1070  dof_handler,
1071  first_vector_component,
1072  boundary_ids,
1073  function_map,
1074  constraints,
1075  mapping,
1076  refinement_edge_indices,
1077  level);
1078  }
1079 
1080  } // namespace VectorTools
1081 } // namespace dealii
1082 
1083 #endif
1084 
1085 
1086 // deal.II versions up to 9.5 had a poorly designed interface of the
1087 // SphericalManifold class that made it impossible for us to use.
1088 // This file thus contains a copy of it.
1089 #if !DEAL_II_VERSION_GTE(9,6,0)
1090 
1091 #include <deal.II/grid/manifold.h>
1092 #include <deal.II/grid/manifold_lib.h>
1093 
1094 namespace aspect
1095 {
1096  using namespace dealii;
1097 
1106  template <int dim, int spacedim = dim>
1107  class SphericalManifold : public Manifold<dim, spacedim>
1108  {
1109  public:
1116  SphericalManifold(const Point<spacedim> center = Point<spacedim>());
1117 
1121  virtual std::unique_ptr<Manifold<dim, spacedim>>
1122  clone() const override;
1123 
1131  virtual Point<spacedim>
1132  get_intermediate_point(const Point<spacedim> &p1,
1133  const Point<spacedim> &p2,
1134  const double w) const override;
1135 
1140  virtual Tensor<1, spacedim>
1141  get_tangent_vector(const Point<spacedim> &x1,
1142  const Point<spacedim> &x2) const override;
1143 
1147  virtual Tensor<1, spacedim>
1148  normal_vector(
1149  const typename Triangulation<dim, spacedim>::face_iterator &face,
1150  const Point<spacedim> &p) const override;
1151 
1155  virtual void
1156  get_normals_at_vertices(
1157  const typename Triangulation<dim, spacedim>::face_iterator &face,
1158  typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
1159  const override;
1160 
1175  virtual void
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;
1179 
1184  virtual Point<spacedim>
1185  get_new_point(const ArrayView<const Point<spacedim>> &vertices,
1186  const ArrayView<const double> &weights) const override;
1187 
1191  const Point<spacedim> center;
1192 
1193  private:
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;
1204 
1222  void
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;
1226 
1230  const PolarManifold<spacedim> polar_manifold;
1231  };
1232 }
1233 
1234 #else
1235 
1236 // For sufficiently new deal.II versions, we can use the deal.II class, but to
1237 // avoid name clashes, we have to import the class into namespace aspect. Once
1238 // we rely on these sufficiently new versions of deal.II, we can not only remove
1239 // the code above, but also the following lines, and in all places where we
1240 // reference 'aspect::SphericalManifold' simply use 'SphericalManifold' instead
1241 // (which then refers to the deal.II class).
1242 
1243 #include <deal.II/grid/manifold_lib.h>
1244 namespace aspect
1245 {
1246  using ::SphericalManifold;
1247 }
1248 
1249 #endif
1250 
1251 #endif
void compute_nonzero_normal_flux_constraints_active_or_level(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, const std::map< types::boundary_id, const Function< spacedim > *> &function_map, AffineConstraints< double > &constraints, const Mapping< dim, spacedim > &mapping, const IndexSet &refinement_edge_indices=IndexSet(), const unsigned int level=numbers::invalid_unsigned_int)
Definition: compat.h:574
bool operator==(const VectorDoFTuple< dim > &other) const
Definition: compat.h:84
void map_dofs_to_normal_vectors_and_normal_fluxes(const typename DoFHandler< dim, spacedim >::cell_iterator &cell, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, const std::map< types::boundary_id, const Function< spacedim > *> &function_map, hp::FEFaceValues< dim, spacedim > &x_fe_face_values, const unsigned int n_dofs, const IndexSet &refinement_edge_indices, const unsigned int level, std::multimap< internal::VectorDoFTuple< dim >, std::pair< Tensor< 1, dim >, typename DoFHandler< dim, spacedim >::cell_iterator >> &dof_to_normals_map, std::map< internal::VectorDoFTuple< dim >, Vector< double >> &dof_vector_to_b_values)
Definition: compat.h:393
void add_tangentiality_constraints(const VectorDoFTuple< dim > &dof_indices, const Tensor< 1, dim > &tangent_vector, AffineConstraints< double > &constraints, const Vector< double > &b_values=Vector< double >(dim))
Definition: compat.h:333
void compute_no_normal_flux_constraints_on_level(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, AffineConstraints< double > &constraints, const Mapping< dim, spacedim > &mapping, const IndexSet &refinement_edge_indices, const unsigned int level)
Definition: compat.h:1056
std::array< types::global_dof_index, dim > dof_indices
Definition: compat.h:68
MGTransferMatrixFree< dim, NumberType > MGTransferMF
Definition: compat.h:45
void add_constraint(const VectorDoFTuple< dim > &dof_indices, const Tensor< 1, dim > &constraining_vector, AffineConstraints< double > &constraints, const double inhomogeneity=0)
Definition: compat.h:120
const PolarManifold< spacedim > polar_manifold
Definition: compat.h:1230
bool operator!=(const VectorDoFTuple< dim > &other) const
Definition: compat.h:90
const Point< spacedim > center
Definition: compat.h:1191
Definition: compat.h:30
Definition: compat.h:42