ASPECT
compat.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2015 - 2022 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 #if DEAL_II_VERSION_GTE(9,4,0)
34 
36 
37 #else
38 
39  inline MPI_Datatype
40  mpi_type_id(const char *)
41  {
42  return MPI_CHAR;
43  }
44 
49  template <typename T>
50  void
51  broadcast(T *buffer,
52  const size_t count,
53  const unsigned int root,
54  const MPI_Comm &comm)
55  {
56  Assert(root < ::Utilities::MPI::n_mpi_processes(comm),
57  ::ExcMessage("Invalid root rank specified."));
58 
59  // MPI_Bcast's count is a signed int, so send at most 2^31 in each
60  // iteration:
61  const size_t max_send_count = std::numeric_limits<signed int>::max();
62 
63  size_t total_sent_count = 0;
64  while (total_sent_count < count)
65  {
66  const size_t current_count =
67  std::min(count - total_sent_count, max_send_count);
68 
69  const int ierr = MPI_Bcast(buffer + total_sent_count,
70  current_count,
71  mpi_type_id(buffer),
72  root,
73  comm);
74  AssertThrowMPI(ierr);
75  total_sent_count += current_count;
76  }
77 
78  }
79 #endif
80 
81 }
82 
83 
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>
88 namespace dealii
89 {
90  namespace VectorTools
91  {
92  namespace internal
93  {
98  template <int dim>
100  {
101  types::global_dof_index dof_indices[dim];
102 
104  {
105  for (unsigned int i = 0; i < dim; ++i)
106  dof_indices[i] = numbers::invalid_dof_index;
107  }
108 
109 
110  bool
111  operator<(const VectorDoFTuple<dim> &other) const
112  {
113  for (unsigned int i = 0; i < dim; ++i)
114  if (dof_indices[i] < other.dof_indices[i])
115  return true;
116  else if (dof_indices[i] > other.dof_indices[i])
117  return false;
118  return false;
119  }
120 
121  bool
122  operator==(const VectorDoFTuple<dim> &other) const
123  {
124  for (unsigned int i = 0; i < dim; ++i)
125  if (dof_indices[i] != other.dof_indices[i])
126  return false;
127 
128  return true;
129  }
130 
131  bool
132  operator!=(const VectorDoFTuple<dim> &other) const
133  {
134  return !(*this == other);
135  }
136  };
137 
138 
139  template <int dim>
140  std::ostream &
141  operator<<(std::ostream &out, const VectorDoFTuple<dim> &vdt)
142  {
143  for (unsigned int d = 0; d < dim; ++d)
144  out << vdt.dof_indices[d] << (d < dim - 1 ? " " : "");
145  return out;
146  }
147 
148 
149 
160  template <int dim>
161  void
163  const Tensor<1, dim> &constraining_vector,
164  AffineConstraints<double> &constraints,
165  const double inhomogeneity = 0)
166  {
167  // choose the DoF that has the largest component in the
168  // constraining_vector as the one to be constrained as this makes the
169  // process stable in cases where the constraining_vector has the form
170  // n=(1,0) or n=(0,1)
171  //
172  // we get constraints of the form x0 = a_1*x1 + a2*x2 + ... if one of
173  // the weights is essentially zero then skip this part. the
174  // AffineConstraints can also deal with cases like x0 = 0 if
175  // necessary
176  //
177  // there is a problem if we have a normal vector of the form
178  // (a,a,small) or (a,a,a). Depending on round-off we may choose the
179  // first or second component (or third, in the latter case) as the
180  // largest one, and depending on our choice one or another degree of
181  // freedom will be constrained. On a single processor this is not
182  // much of a problem, but it's a nightmare when we run in parallel
183  // and two processors disagree on which DoF should be constrained.
184  // This led to an incredibly difficult to find bug in step-32 when
185  // running in parallel with 9 or more processors.
186  //
187  // in practice, such normal vectors of the form (a,a,small) or
188  // (a,a,a) happen not infrequently since they lie on the diagonals
189  // where vertices frequently happen to land upon mesh refinement if
190  // one starts from a symmetric and regular body. we work around this
191  // problem in the following way: if we have a normal vector of the
192  // form (a,b) (similarly algorithm in 3d), we choose 'a' as the
193  // largest coefficient not if a>b but if a>b+1e-10. this shifts the
194  // problem away from the frequently visited diagonal to a line that
195  // is off the diagonal. there will of course be problems where the
196  // exact values of a and b differ by exactly 1e-10 and we get into
197  // the same instability, but from a practical viewpoint such problems
198  // should be much rarer. in particular, meshes have to be very fine
199  // for a vertex to land on this line if the original body had a
200  // vertex on the diagonal as well
201  switch (dim)
202  {
203  case 2:
204  {
205  if (std::fabs(constraining_vector[0]) >
206  std::fabs(constraining_vector[1]) + 1e-10)
207  {
208  if (!constraints.is_constrained(dof_indices.dof_indices[0]) &&
209  constraints.can_store_line(dof_indices.dof_indices[0]))
210  {
211  constraints.add_line(dof_indices.dof_indices[0]);
212 
213  if (std::fabs(constraining_vector[1] /
214  constraining_vector[0]) >
215  std::numeric_limits<double>::epsilon())
216  constraints.add_entry(dof_indices.dof_indices[0],
217  dof_indices.dof_indices[1],
218  -constraining_vector[1] /
219  constraining_vector[0]);
220 
221  if (std::fabs(inhomogeneity / constraining_vector[0]) >
222  std::numeric_limits<double>::epsilon())
223  constraints.set_inhomogeneity(
224  dof_indices.dof_indices[0],
225  inhomogeneity / constraining_vector[0]);
226  }
227  }
228  else
229  {
230  if (!constraints.is_constrained(dof_indices.dof_indices[1]) &&
231  constraints.can_store_line(dof_indices.dof_indices[1]))
232  {
233  constraints.add_line(dof_indices.dof_indices[1]);
234 
235  if (std::fabs(constraining_vector[0] /
236  constraining_vector[1]) >
237  std::numeric_limits<double>::epsilon())
238  constraints.add_entry(dof_indices.dof_indices[1],
239  dof_indices.dof_indices[0],
240  -constraining_vector[0] /
241  constraining_vector[1]);
242 
243  if (std::fabs(inhomogeneity / constraining_vector[1]) >
244  std::numeric_limits<double>::epsilon())
245  constraints.set_inhomogeneity(
246  dof_indices.dof_indices[1],
247  inhomogeneity / constraining_vector[1]);
248  }
249  }
250  break;
251  }
252 
253  case 3:
254  {
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))
259  {
260  if (!constraints.is_constrained(dof_indices.dof_indices[0]) &&
261  constraints.can_store_line(dof_indices.dof_indices[0]))
262  {
263  constraints.add_line(dof_indices.dof_indices[0]);
264 
265  if (std::fabs(constraining_vector[1] /
266  constraining_vector[0]) >
267  std::numeric_limits<double>::epsilon())
268  constraints.add_entry(dof_indices.dof_indices[0],
269  dof_indices.dof_indices[1],
270  -constraining_vector[1] /
271  constraining_vector[0]);
272 
273  if (std::fabs(constraining_vector[2] /
274  constraining_vector[0]) >
275  std::numeric_limits<double>::epsilon())
276  constraints.add_entry(dof_indices.dof_indices[0],
277  dof_indices.dof_indices[2],
278  -constraining_vector[2] /
279  constraining_vector[0]);
280 
281  if (std::fabs(inhomogeneity / constraining_vector[0]) >
282  std::numeric_limits<double>::epsilon())
283  constraints.set_inhomogeneity(
284  dof_indices.dof_indices[0],
285  inhomogeneity / constraining_vector[0]);
286  }
287  }
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))
292  {
293  if (!constraints.is_constrained(dof_indices.dof_indices[1]) &&
294  constraints.can_store_line(dof_indices.dof_indices[1]))
295  {
296  constraints.add_line(dof_indices.dof_indices[1]);
297 
298  if (std::fabs(constraining_vector[0] /
299  constraining_vector[1]) >
300  std::numeric_limits<double>::epsilon())
301  constraints.add_entry(dof_indices.dof_indices[1],
302  dof_indices.dof_indices[0],
303  -constraining_vector[0] /
304  constraining_vector[1]);
305 
306  if (std::fabs(constraining_vector[2] /
307  constraining_vector[1]) >
308  std::numeric_limits<double>::epsilon())
309  constraints.add_entry(dof_indices.dof_indices[1],
310  dof_indices.dof_indices[2],
311  -constraining_vector[2] /
312  constraining_vector[1]);
313 
314  if (std::fabs(inhomogeneity / constraining_vector[1]) >
315  std::numeric_limits<double>::epsilon())
316  constraints.set_inhomogeneity(
317  dof_indices.dof_indices[1],
318  inhomogeneity / constraining_vector[1]);
319  }
320  }
321  else
322  {
323  if (!constraints.is_constrained(dof_indices.dof_indices[2]) &&
324  constraints.can_store_line(dof_indices.dof_indices[2]))
325  {
326  constraints.add_line(dof_indices.dof_indices[2]);
327 
328  if (std::fabs(constraining_vector[0] /
329  constraining_vector[2]) >
330  std::numeric_limits<double>::epsilon())
331  constraints.add_entry(dof_indices.dof_indices[2],
332  dof_indices.dof_indices[0],
333  -constraining_vector[0] /
334  constraining_vector[2]);
335 
336  if (std::fabs(constraining_vector[1] /
337  constraining_vector[2]) >
338  std::numeric_limits<double>::epsilon())
339  constraints.add_entry(dof_indices.dof_indices[2],
340  dof_indices.dof_indices[1],
341  -constraining_vector[1] /
342  constraining_vector[2]);
343 
344  if (std::fabs(inhomogeneity / constraining_vector[2]) >
345  std::numeric_limits<double>::epsilon())
346  constraints.set_inhomogeneity(
347  dof_indices.dof_indices[2],
348  inhomogeneity / constraining_vector[2]);
349  }
350  }
351 
352  break;
353  }
354 
355  default:
356  Assert(false, ExcNotImplemented());
357  }
358  }
359 
360 
361 
373  template <int dim>
374  void
376  const VectorDoFTuple<dim> &dof_indices,
377  const Tensor<1, dim> &tangent_vector,
378  AffineConstraints<double> &constraints,
379  const Vector<double> &b_values = Vector<double>(dim))
380  {
381  // choose the DoF that has the
382  // largest component in the
383  // tangent_vector as the
384  // independent component, and
385  // then constrain the others to
386  // it. specifically, if, say,
387  // component 0 of the tangent
388  // vector t is largest by
389  // magnitude, then
390  // x1=(b[1]*t[0]-b[0]*t[1])/t[0]+t[1]/t[0]*x_0, etc.
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;
396 
397  // then constrain all of the
398  // other degrees of freedom in
399  // terms of the one just found
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]))
404  {
405  constraints.add_line(dof_indices.dof_indices[d]);
406 
407  if (std::fabs(tangent_vector[d] /
408  tangent_vector[largest_component]) >
409  std::numeric_limits<double>::epsilon())
410  constraints.add_entry(
411  dof_indices.dof_indices[d],
412  dof_indices.dof_indices[largest_component],
413  tangent_vector[d] / tangent_vector[largest_component]);
414 
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];
419 
420  if (std::fabs(inhomogeneity) >
421  std::numeric_limits<double>::epsilon())
422  constraints.set_inhomogeneity(dof_indices.dof_indices[d],
423  inhomogeneity);
424  }
425  }
426 
427 
428 
433  template <int dim, int spacedim>
434  void
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,
439  const std::map<types::boundary_id, const Function<spacedim> *>
440  &function_map,
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,
445  std::multimap<
447  std::pair<Tensor<1, dim>,
448  typename DoFHandler<dim, spacedim>::cell_iterator>>
449  &dof_to_normals_map,
450  std::map<internal::VectorDoFTuple<dim>, Vector<double>>
451  &dof_vector_to_b_values)
452  {
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())) !=
456  boundary_ids.end())
457  {
458  const FiniteElement<dim> &fe = cell->get_fe();
459  typename DoFHandler<dim, spacedim>::level_face_iterator face =
460  cell->face(face_no);
461 
462  std::vector<types::global_dof_index> face_dofs;
463  // get the indices of the dofs on this cell...
464  face_dofs.resize(fe.n_dofs_per_face(face_no));
465 
466  if (level != numbers::invalid_unsigned_int)
467  face->get_mg_dof_indices(level,
468  face_dofs,
469  cell->active_fe_index());
470  else
471  face->get_dof_indices(face_dofs, cell->active_fe_index());
472 
473  x_fe_face_values.reinit(cell, face_no);
474  const FEFaceValues<dim> &fe_values =
475  x_fe_face_values.get_present_fe_values();
476 
477  const auto &face_quadrature_collection =
478  x_fe_face_values.get_quadrature_collection();
479 
480  // then identify which of them correspond to the selected set of
481  // vector components
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)
485  // Refinement edge indices are going to be constrained to 0
486  // during a multigrid cycle and do not need no-normal-flux
487  // constraints, so skip them:
488  if (level == numbers::invalid_unsigned_int ||
489  !refinement_edge_indices.is_element(face_dofs[i]))
490  {
491  // find corresponding other components of vector
492  internal::VectorDoFTuple<dim> vector_dofs;
493  vector_dofs.dof_indices[0] = face_dofs[i];
494 
495  Assert(
496  first_vector_component + dim <= fe.n_components(),
497  ExcMessage(
498  "Error: the finite element does not have enough components "
499  "to define a normal direction."));
500 
501  for (unsigned int k = 0; k < fe.n_dofs_per_face(face_no);
502  ++k)
503  if ((k != i) &&
504  (face_quadrature_collection[cell->active_fe_index()]
505  .point(k) ==
506  face_quadrature_collection[cell->active_fe_index()]
507  .point(i)) &&
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))
512  vector_dofs.dof_indices
513  [fe.face_system_to_component_index(k, face_no).first -
514  first_vector_component] = face_dofs[k];
515 
516  for (unsigned int d = 0; d < dim; ++d)
517  Assert(vector_dofs.dof_indices[d] < n_dofs,
518  ExcInternalError());
519  (void)n_dofs;
520 
521  // we need the normal vector on this face. we know that
522  // it is a vector of length 1 but at least with higher
523  // order mappings it isn't always possible to guarantee
524  // that each component is exact up to zero tolerance. in
525  // particular, as shown in the deal.II/no_flux_06 test,
526  // if we just take the normal vector as given by the
527  // fe_values object, we can get entries in the normal
528  // vectors of the unit cube that have entries up to
529  // several times 1e-14.
530  //
531  // the problem with this is that this later yields
532  // constraints that are circular (e.g., in the testcase,
533  // we get constraints of the form
534  //
535  // x22 = 2.93099e-14*x21 + 2.93099e-14*x23
536  // x21 = -2.93099e-14*x22 + 2.93099e-14*x21
537  //
538  // in both of these constraints, the small numbers
539  // should be zero and the constraints should simply be
540  // x22 = x21 = 0
541  //
542  // to achieve this, we utilize that we know that the
543  // normal vector has (or should have) length 1 and that
544  // we can simply set small elements to zero (without
545  // having to check that they are small *relative to
546  // something else*). we do this and then normalize the
547  // length of the vector back to one, just to be on the
548  // safe side
549  //
550  // one more point: we would like to use the "real"
551  // normal vector here, as provided by the boundary
552  // description and as opposed to what we get from the
553  // FEValues object. we do this in the immediately next
554  // line, but as is obvious, the boundary only has a
555  // vague idea which side of a cell it is on -- indicated
556  // by the face number. in other words, it may provide
557  // the inner or outer normal. by and large, there is no
558  // harm from this, since the tangential vector we
559  // compute is still the same. however, we do average
560  // over normal vectors from adjacent cells and if they
561  // have recorded normal vectors from the inside once and
562  // from the outside the other time, then this averaging
563  // is going to run into trouble. as a consequence we ask
564  // the mapping after all for its normal vector, but we
565  // only ask it so that we can possibly correct the sign
566  // of the normal vector provided by the boundary if they
567  // should point in different directions. this is the
568  // case in tests/deal.II/no_flux_11.
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)
573  normal_vector *= -1;
574  Assert(std::fabs(normal_vector.norm() - 1) < 1e-14,
575  ExcInternalError());
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();
580 
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);
584 
585  // now enter the (dofs,(normal_vector,cell)) entry into
586  // the map
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));
592 
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()
597  << std::endl
598  << " normal=" << normal_vector << std::endl;
599 #endif
600  }
601  }
602  }
603 
604 
605 
614  template <int dim, int spacedim>
615  void
617  const DoFHandler<dim, spacedim> &dof_handler,
618  const unsigned int first_vector_component,
619  const std::set<types::boundary_id> &boundary_ids,
620  const std::map<types::boundary_id, const Function<spacedim> *>
621  &function_map,
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)
626  {
627  Assert(dim > 1,
628  ExcMessage("This function is not useful in 1d because it amounts "
629  "to imposing Dirichlet values on the vector-valued "
630  "quantity."));
631 
632  // create FE and mapping collections for all elements in use by this
633  // DoFHandler
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);
639 
640  // TODO: the implementation makes the assumption that all faces have the
641  // same number of dofs
642  AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
643  const unsigned int face_no = 0;
644 
645  // now also create a quadrature collection for the faces of a cell. fill
646  // it with a quadrature formula with the support points on faces for each
647  // FE
648  hp::QCollection<dim - 1> face_quadrature_collection;
649  for (unsigned int i = 0; i < fe_collection.size(); ++i)
650  {
651  const std::vector<Point<dim - 1>> &unit_support_points =
652  fe_collection[i].get_unit_face_support_points(face_no);
653 
654  Assert(unit_support_points.size() ==
655  fe_collection[i].n_dofs_per_face(face_no),
656  ExcInternalError());
657 
658  face_quadrature_collection.push_back(
659  Quadrature<dim - 1>(unit_support_points));
660  }
661 
662  // now create the object with which we will generate the normal vectors
663  hp::FEFaceValues<dim, spacedim> x_fe_face_values(
664  mapping_collection,
665  fe_collection,
666  face_quadrature_collection,
667  update_quadrature_points | update_normal_vectors);
668 
669  // have a map that stores normal vectors for each vector-dof tuple we want
670  // to constrain. since we can get at the same vector dof tuple more than
671  // once (for example if it is located at a vertex that we visit from all
672  // adjacent cells), we will want to average later on the normal vectors
673  // computed on different cells as described in the documentation of this
674  // function. however, we can only average if the contributions came from
675  // different cells, whereas we want to constrain twice or more in case the
676  // contributions came from different faces of the same cell
677  // (i.e. constrain not just the *average normal direction* but *all normal
678  // directions* we find). consequently, we also have to store which cell a
679  // normal vector was computed on
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;
686 
687  DoFToNormalsMap dof_to_normals_map;
688 
689  const unsigned int n_dof = dof_handler.n_dofs();
690 
691  if (level == numbers::invalid_unsigned_int)
692  {
693  // active cells
694  for (const auto &cell : dof_handler.active_cell_iterators())
695  if (!cell->is_artificial())
696  {
698  cell,
699  first_vector_component,
700  boundary_ids,
701  function_map,
702  x_fe_face_values,
703  n_dof,
704  refinement_edge_indices,
705  level,
706  dof_to_normals_map,
707  dof_vector_to_b_values);
708  }
709  }
710  else
711  {
712  // level cells
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)
717  {
719  cell,
720  first_vector_component,
721  boundary_ids,
722  function_map,
723  x_fe_face_values,
724  n_dof,
725  refinement_edge_indices,
726  level,
727  dof_to_normals_map,
728  dof_vector_to_b_values);
729  }
730  }
731 
732  // Now do something with the collected information. To this end, loop
733  // through all sets of pairs (dofs,normal_vector) and identify which
734  // entries belong to the same set of dofs and then do as described in the
735  // documentation, i.e. either average the normal vector or don't for this
736  // particular set of dofs
737  typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
738 
739  while (p != dof_to_normals_map.end())
740  {
741  // first find the range of entries in the multimap that corresponds to
742  // the same vector-dof tuple. as usual, we define the range
743  // half-open. the first entry of course is 'p'
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)
747  {
748  same_dof_range[1] = p;
749  break;
750  }
751  if (p == dof_to_normals_map.end())
752  same_dof_range[1] = dof_to_normals_map.end();
753 
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];
759  ++q)
760  std::cout << " " << q->second.first << " from cell "
761  << q->second.second << std::endl;
762 #endif
763 
764 
765  // now compute the reverse mapping: for each of the cells that
766  // contributed to the current set of vector dofs, add up the normal
767  // vectors. the values of the map are pairs of normal vectors and
768  // number of cells that have contributed
769  using CellToNormalsMap =
770  std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
771  std::pair<Tensor<1, dim>, unsigned int>>;
772 
773  CellToNormalsMap cell_to_normals_map;
774  for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
775  q != same_dof_range[1];
776  ++q)
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);
781  else
782  {
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;
787 
788  Assert(old_count > 0, ExcInternalError());
789 
790  // in the same entry, store again the now averaged normal vector
791  // and the new count
792  cell_to_normals_map[q->second.second] =
793  std::make_pair((old_normal * old_count + q->second.first) /
794  (old_count + 1),
795  old_count + 1);
796  }
797  Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
798 
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();
804  ++x)
805  std::cout << " " << x->first << " -> (" << x->second.first
806  << ',' << x->second.second << ')' << std::endl;
807 #endif
808 
809  // count the maximum number of contributions from each cell
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();
814  ++x)
815  max_n_contributions_per_cell =
816  std::max(max_n_contributions_per_cell, x->second.second);
817 
818  // verify that each cell can have only contributed at most dim times,
819  // since that is the maximum number of faces that come together at a
820  // single place
821  Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
822 
823  switch (max_n_contributions_per_cell)
824  {
825  // first deal with the case that a number of cells all have
826  // registered that they have a normal vector defined at the
827  // location of a given vector dof, and that each of them have
828  // encountered this vector dof exactly once while looping over all
829  // their faces. as stated in the documentation, this is the case
830  // where we want to simply average over all normal vectors
831  //
832  // the typical case is in 2d where multiple cells meet at one
833  // vertex sitting on the boundary. same in 3d for a vertex that
834  // is associated with only one of the boundary indicators passed
835  // to this function
836  case 1:
837  {
838  // compute the average normal vector from all the ones that
839  // have the same set of dofs. we could add them up and divide
840  // them by the number of additions, or simply normalize them
841  // right away since we want them to have unit length anyway
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();
846  ++x)
847  normal += x->second.first;
848  normal /= normal.norm();
849 
850  // normalize again
851  for (unsigned int d = 0; d < dim; ++d)
852  if (std::fabs(normal[d]) < 1e-13)
853  normal[d] = 0;
854  normal /= normal.norm();
855 
856  // then construct constraints from this:
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];
864  internal::add_constraint(dof_indices,
865  normal,
866  constraints,
867  normal_value);
868 
869  break;
870  }
871 
872  // this is the slightly more complicated case that a single cell
873  // has contributed with exactly DIM normal vectors to the same set
874  // of vector dofs. this is what happens in a corner in 2d and 3d
875  // (but not on an edge in 3d, where we have only 2, i.e. <DIM,
876  // contributions. Here we do not want to average the normal
877  // vectors. Since we have DIM contributions, let's assume (and
878  // verify) that they are in fact all linearly independent; in that
879  // case, all vector components are constrained and we need to set
880  // all of them to the corresponding boundary values
881  case dim:
882  {
883  // assert that indeed only a single cell has contributed
884  Assert(cell_to_normals_map.size() == 1, ExcInternalError());
885 
886  // check linear independence by computing the determinant of
887  // the matrix created from all the normal vectors. if they are
888  // linearly independent, then the determinant is nonzero. if
889  // they are orthogonal, then the matrix is in fact equal to 1
890  // (since they are all unit vectors); make sure the
891  // determinant is larger than 1e-3 to avoid cases where cells
892  // are degenerate
893  {
894  Tensor<2, dim> t;
895 
896  typename DoFToNormalsMap::const_iterator x =
897  same_dof_range[0];
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];
901 
902  Assert(
903  std::fabs(determinant(t)) > 1e-3,
904  ExcMessage(
905  "Found a set of normal vectors that are nearly collinear."));
906  }
907 
908  // so all components of this vector dof are constrained. enter
909  // this into the AffineConstraints object
910  //
911  // ignore dofs already constrained
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]))
921  {
922  const types::global_dof_index line =
923  dof_indices.dof_indices[i];
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]);
928  // no add_entries here
929  }
930 
931  break;
932  }
933 
934  // this is the case of an edge contribution in 3d, i.e. the vector
935  // is constrained in two directions but not the third.
936  default:
937  {
938  Assert(dim >= 3, ExcNotImplemented());
939  Assert(max_n_contributions_per_cell == 2, ExcInternalError());
940 
941  // as described in the documentation, let us first collect
942  // what each of the cells contributed at the current point. we
943  // use a std::list instead of a std::set (which would be more
944  // natural) because std::set requires that the stored elements
945  // are comparable with operator<
946  using CellContributions =
947  std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
948  std::list<Tensor<1, dim>>>;
949  CellContributions cell_contributions;
950 
951  for (typename DoFToNormalsMap::const_iterator q =
952  same_dof_range[0];
953  q != same_dof_range[1];
954  ++q)
955  cell_contributions[q->second.second].push_back(
956  q->second.first);
957  Assert(cell_contributions.size() >= 1, ExcInternalError());
958 
959  // now for each cell that has contributed determine the number
960  // of normal vectors it has contributed. we currently only
961  // implement if this is dim-1 for all cells (if a single cell
962  // has contributed dim, or if all adjacent cells have
963  // contributed 1 normal vector, this is already handled
964  // above).
965  //
966  // we only implement the case that all cells contribute
967  // dim-1 because we assume that we are following an edge
968  // of the domain (think: we are looking at a vertex
969  // located on one of the edges of a refined cube where the
970  // boundary indicators of the two adjacent faces of the
971  // cube are both listed in the set of boundary indicators
972  // passed to this function). in that case, all cells along
973  // that edge of the domain are assumed to have contributed
974  // dim-1 normal vectors. however, there are cases where
975  // this assumption is not justified (see the lengthy
976  // explanation in test no_flux_12.cc) and in those cases
977  // we simply ignore the cell that contributes only
978  // once. this is also discussed at length in the
979  // documentation of this function.
980  //
981  // for each contributing cell compute the tangential vector
982  // that remains unconstrained
983  std::list<Tensor<1, dim>> tangential_vectors;
984  for (typename CellContributions::const_iterator contribution =
985  cell_contributions.begin();
986  contribution != cell_contributions.end();
987  ++contribution)
988  {
989 #ifdef DEBUG_NO_NORMAL_FLUX
990  std::cout
991  << " Treating edge case with dim-1 contributions."
992  << std::endl
993  << " Looking at cell " << contribution->first
994  << " which has contributed these normal vectors:"
995  << std::endl;
996  for (typename std::list<Tensor<1, dim>>::const_iterator
997  t = contribution->second.begin();
998  t != contribution->second.end();
999  ++t)
1000  std::cout << " " << *t << std::endl;
1001 #endif
1002 
1003  // as mentioned above, simply ignore cells that only
1004  // contribute once
1005  if (contribution->second.size() < dim - 1)
1006  continue;
1007 
1008  Tensor<1, dim> normals[dim - 1];
1009  {
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();
1014  ++t, ++index)
1015  normals[index] = *t;
1016  Assert(index == dim - 1, ExcInternalError());
1017  }
1018 
1019  // calculate the tangent as the outer product of the
1020  // normal vectors. since these vectors do not need to be
1021  // orthogonal (think, for example, the case of the
1022  // deal.II/no_flux_07 test: a sheared cube in 3d, with Q2
1023  // elements, where we have constraints from the two normal
1024  // vectors of two faces of the sheared cube that are not
1025  // perpendicular to each other), we have to normalize the
1026  // outer product
1027  Tensor<1, dim> tangent;
1028  switch (dim)
1029  {
1030  case 3:
1031  // take cross product between normals[0] and
1032  // normals[1]. write it in the current form (with
1033  // [dim-2]) to make sure that compilers don't warn
1034  // about out-of-bounds accesses -- the warnings are
1035  // bogus since we get here only for dim==3, but at
1036  // least one isn't quite smart enough to notice this
1037  // and warns when compiling the function in 2d
1038  tangent =
1039  cross_product_3d(normals[0], normals[dim - 2]);
1040  break;
1041  default:
1042  Assert(false, ExcNotImplemented());
1043  }
1044 
1045  Assert(
1046  std::fabs(tangent.norm()) > 1e-12,
1047  ExcMessage(
1048  "Two normal vectors from adjacent faces are almost "
1049  "parallel."));
1050  tangent /= tangent.norm();
1051 
1052  tangential_vectors.push_back(tangent);
1053  }
1054 
1055  // go through the list of tangents and make sure that they all
1056  // roughly point in the same direction as the first one (i.e.
1057  // have an angle less than 90 degrees); if they don't then
1058  // flip their sign
1059  {
1060  const Tensor<1, dim> first_tangent =
1061  tangential_vectors.front();
1062  typename std::list<Tensor<1, dim>>::iterator t =
1063  tangential_vectors.begin();
1064  ++t;
1065  for (; t != tangential_vectors.end(); ++t)
1066  if (*t * first_tangent < 0)
1067  *t *= -1;
1068  }
1069 
1070  // now compute the average tangent and normalize it
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();
1075  ++t)
1076  average_tangent += *t;
1077  average_tangent /= average_tangent.norm();
1078 
1079  // now all that is left is that we add the constraints that
1080  // the vector is parallel to the tangent
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];
1086  average_tangent,
1087  constraints,
1088  b_values);
1089  }
1090  }
1091  }
1092  }
1093 
1094  } // namespece internal
1095 
1096  template <int dim, int spacedim>
1097  void
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)
1106  {
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;
1112  dof_handler,
1113  first_vector_component,
1114  boundary_ids,
1115  function_map,
1116  constraints,
1117  mapping,
1118  refinement_edge_indices,
1119  level);
1120  }
1121 
1122  } // namespace VectorTools
1123 } // namespace dealii
1124 
1125 #endif
1126 
1127 
1128 #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:616
bool operator==(const VectorDoFTuple< dim > &other) const
Definition: compat.h:122
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:435
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:375
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:1098
void broadcast(T *buffer, const size_t count, const unsigned int root, const MPI_Comm &comm)
Definition: compat.h:51
void add_constraint(const VectorDoFTuple< dim > &dof_indices, const Tensor< 1, dim > &constraining_vector, AffineConstraints< double > &constraints, const double inhomogeneity=0)
Definition: compat.h:162
bool operator!=(const VectorDoFTuple< dim > &other) const
Definition: compat.h:132
types::global_dof_index dof_indices[dim]
Definition: compat.h:101
Definition: compat.h:30
Definition: compat.h:88
MPI_Datatype mpi_type_id(const char *)
Definition: compat.h:40