Reference documentation for deal.II version GIT 3a483e2de5 2022-11-28 08:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mg_tools.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/logstream.h>
20 
22 
25 #include <deal.II/dofs/dof_tools.h>
26 
28 #include <deal.II/fe/fe.h>
29 
31 #include <deal.II/grid/tria.h>
33 
38 
42 
44 
45 #include <algorithm>
46 #include <numeric>
47 #include <vector>
48 
50 
51 
52 namespace MGTools
53 {
54  // specializations for 1D
55  template <>
56  void
58  const unsigned int,
59  std::vector<unsigned int> &,
60  const DoFTools::Coupling)
61  {
62  Assert(false, ExcNotImplemented());
63  }
64 
65 
66 
67  template <>
68  void
70  const unsigned int,
71  std::vector<unsigned int> &,
74  {
75  Assert(false, ExcNotImplemented());
76  }
77 
78 
79 
80  template <>
81  void
83  const unsigned int,
84  std::vector<unsigned int> &,
85  const DoFTools::Coupling)
86  {
87  Assert(false, ExcNotImplemented());
88  }
89 
90 
91  template <>
92  void
94  const unsigned int,
95  std::vector<unsigned int> &,
98  {
99  Assert(false, ExcNotImplemented());
100  }
101 
102 
103 
104  // Template for 2D and 3D. For 1D see specialization above
105  template <int dim, int spacedim>
106  void
108  const unsigned int level,
109  std::vector<unsigned int> & row_lengths,
110  const DoFTools::Coupling flux_coupling)
111  {
112  Assert(row_lengths.size() == dofs.n_dofs(),
113  ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));
114 
115  // Function starts here by
116  // resetting the counters.
117  std::fill(row_lengths.begin(), row_lengths.end(), 0);
118 
119  std::vector<bool> face_touched(dim == 2 ?
120  dofs.get_triangulation().n_raw_lines() :
121  dofs.get_triangulation().n_raw_quads());
122 
123  std::vector<types::global_dof_index> cell_indices;
124  std::vector<types::global_dof_index> neighbor_indices;
125 
126  // We loop over cells and go from
127  // cells to lower dimensional
128  // objects. This is the only way to
129  // cope with the fact, that an
130  // unknown number of cells may
131  // share an object of dimension
132  // smaller than dim-1.
133  for (const auto &cell : dofs.cell_iterators_on_level(level))
134  {
135  const FiniteElement<dim> &fe = cell->get_fe();
136  cell_indices.resize(fe.n_dofs_per_cell());
137  cell->get_mg_dof_indices(cell_indices);
138  unsigned int i = 0;
139  // First, dofs on
140  // vertices. We assume that
141  // each vertex dof couples
142  // with all dofs on
143  // adjacent grid cells.
144 
145  // Adding all dofs of the cells
146  // will add dofs of the faces
147  // of the cell adjacent to the
148  // vertex twice. Therefore, we
149  // subtract these here and add
150  // them in a loop over the
151  // faces below.
152 
153  // in 1d, faces and vertices
154  // are identical. Nevertheless,
155  // this will only work if
156  // dofs_per_face is zero and
157  // n_dofs_per_vertex() is
158  // arbitrary, not the other way
159  // round.
160  // TODO: This assumes that the dofs per face on all faces coincide!
161  const unsigned int face_no = 0;
162 
163  Assert(fe.reference_cell() == ReferenceCells::get_hypercube<dim>(),
165 
166  unsigned int increment =
167  fe.n_dofs_per_cell() - dim * fe.n_dofs_per_face(face_no);
168  while (i < fe.get_first_line_index())
169  row_lengths[cell_indices[i++]] += increment;
170  // From now on, if an object is
171  // a cell, its dofs only couple
172  // inside the cell. Since the
173  // faces are handled below, we
174  // have to subtract ALL faces
175  // in this case.
176 
177  // In all other cases we
178  // subtract adjacent faces to be
179  // added in the loop below.
180  increment =
181  (dim > 1) ?
182  fe.n_dofs_per_cell() - (dim - 1) * fe.n_dofs_per_face(face_no) :
183  fe.n_dofs_per_cell() -
185  while (i < fe.get_first_quad_index(face_no))
186  row_lengths[cell_indices[i++]] += increment;
187 
188  // Now quads in 2D and 3D
189  increment =
190  (dim > 2) ?
191  fe.n_dofs_per_cell() - (dim - 2) * fe.n_dofs_per_face(face_no) :
192  fe.n_dofs_per_cell() -
194  while (i < fe.get_first_hex_index())
195  row_lengths[cell_indices[i++]] += increment;
196  // Finally, cells in 3D
198  fe.n_dofs_per_face(face_no);
199  while (i < fe.n_dofs_per_cell())
200  row_lengths[cell_indices[i++]] += increment;
201 
202  // At this point, we have
203  // counted all dofs
204  // contributiong from cells
205  // coupled topologically to the
206  // adjacent cells, but we
207  // subtracted some faces.
208 
209  // Now, let's go by the faces
210  // and add the missing
211  // contribution as well as the
212  // flux contributions.
213  for (const unsigned int iface : GeometryInfo<dim>::face_indices())
214  {
215  bool level_boundary = cell->at_boundary(iface);
217  if (!level_boundary)
218  {
219  neighbor = cell->neighbor(iface);
220  if (static_cast<unsigned int>(neighbor->level()) != level)
221  level_boundary = true;
222  }
223 
224  if (level_boundary)
225  {
226  for (unsigned int local_dof = 0;
227  local_dof < fe.n_dofs_per_cell();
228  ++local_dof)
229  row_lengths[cell_indices[local_dof]] +=
230  fe.n_dofs_per_face(face_no);
231  continue;
232  }
233 
234  const FiniteElement<dim> &nfe = neighbor->get_fe();
236  cell->face(iface);
237 
238  // Flux couplings are
239  // computed from both sides
240  // for simplicity.
241 
242  // The dofs on the common face
243  // will be handled below,
244  // therefore, we subtract them
245  // here.
246  if (flux_coupling != DoFTools::none)
247  {
248  const unsigned int dof_increment =
249  nfe.n_dofs_per_cell() - nfe.n_dofs_per_face(face_no);
250  for (unsigned int local_dof = 0;
251  local_dof < fe.n_dofs_per_cell();
252  ++local_dof)
253  row_lengths[cell_indices[local_dof]] += dof_increment;
254  }
255 
256  // Do this only once per
257  // face.
258  if (face_touched[face->index()])
259  continue;
260  face_touched[face->index()] = true;
261 
262  // At this point, we assume
263  // that each cell added its
264  // dofs minus the face to
265  // the couplings of the
266  // face dofs. Since we
267  // subtracted two faces, we
268  // have to re-add one.
269 
270  // If one side of the face
271  // is refined, all the fine
272  // face dofs couple with
273  // the coarse one.
274  neighbor_indices.resize(nfe.n_dofs_per_cell());
275  neighbor->get_mg_dof_indices(neighbor_indices);
276  for (unsigned int local_dof = 0; local_dof < fe.n_dofs_per_cell();
277  ++local_dof)
278  row_lengths[cell_indices[local_dof]] +=
279  nfe.n_dofs_per_face(face_no);
280  for (unsigned int local_dof = 0; local_dof < nfe.n_dofs_per_cell();
281  ++local_dof)
282  row_lengths[neighbor_indices[local_dof]] +=
283  fe.n_dofs_per_face(face_no);
284  }
285  }
286  }
287 
288 
289  // This is the template for 2D and 3D. See version for 1D above
290  template <int dim, int spacedim>
291  void
293  const unsigned int level,
294  std::vector<unsigned int> & row_lengths,
295  const Table<2, DoFTools::Coupling> &couplings,
296  const Table<2, DoFTools::Coupling> &flux_couplings)
297  {
298  Assert(row_lengths.size() == dofs.n_dofs(),
299  ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));
300 
301  // Function starts here by
302  // resetting the counters.
303  std::fill(row_lengths.begin(), row_lengths.end(), 0);
304 
305  std::vector<bool> face_touched(dim == 2 ?
306  dofs.get_triangulation().n_raw_lines() :
307  dofs.get_triangulation().n_raw_quads());
308 
309  std::vector<types::global_dof_index> cell_indices;
310  std::vector<types::global_dof_index> neighbor_indices;
311 
312  // We have to translate the
313  // couplings from components to
314  // blocks, so this works for
315  // nonprimitive elements as well.
316  std::vector<Table<2, DoFTools::Coupling>> couple_cell;
317  std::vector<Table<2, DoFTools::Coupling>> couple_face;
318  DoFTools::convert_couplings_to_blocks(dofs, couplings, couple_cell);
319  DoFTools::convert_couplings_to_blocks(dofs, flux_couplings, couple_face);
320 
321  // We loop over cells and go from
322  // cells to lower dimensional
323  // objects. This is the only way to
324  // cope with the fact, that an
325  // unknown number of cells may
326  // share an object of dimension
327  // smaller than dim-1.
328  for (const auto &cell : dofs.cell_iterators_on_level(level))
329  {
330  const FiniteElement<dim> &fe = cell->get_fe();
331  const unsigned int fe_index = cell->active_fe_index();
332 
333 
334  // TODO: This assumes that the dofs per face on all faces coincide!
335  const unsigned int face_no = 0;
336  Assert(fe.reference_cell() == ReferenceCells::get_hypercube<dim>(),
338 
339  Assert(couplings.n_rows() == fe.n_components(),
340  ExcDimensionMismatch(couplings.n_rows(), fe.n_components()));
341  Assert(couplings.n_cols() == fe.n_components(),
342  ExcDimensionMismatch(couplings.n_cols(), fe.n_components()));
343  Assert(flux_couplings.n_rows() == fe.n_components(),
344  ExcDimensionMismatch(flux_couplings.n_rows(),
345  fe.n_components()));
346  Assert(flux_couplings.n_cols() == fe.n_components(),
347  ExcDimensionMismatch(flux_couplings.n_cols(),
348  fe.n_components()));
349 
350  cell_indices.resize(fe.n_dofs_per_cell());
351  cell->get_mg_dof_indices(cell_indices);
352  unsigned int i = 0;
353  // First, dofs on
354  // vertices. We assume that
355  // each vertex dof couples
356  // with all dofs on
357  // adjacent grid cells.
358 
359  // Adding all dofs of the cells
360  // will add dofs of the faces
361  // of the cell adjacent to the
362  // vertex twice. Therefore, we
363  // subtract these here and add
364  // them in a loop over the
365  // faces below.
366 
367  // in 1d, faces and vertices
368  // are identical. Nevertheless,
369  // this will only work if
370  // dofs_per_face is zero and
371  // n_dofs_per_vertex() is
372  // arbitrary, not the other way
373  // round.
374  unsigned int increment;
375  while (i < fe.get_first_line_index())
376  {
377  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
378  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
379  ++mult)
380  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
381  fe.first_block_of_base(base) +
382  mult) != DoFTools::none)
383  {
384  increment =
385  fe.base_element(base).n_dofs_per_cell() -
386  dim * fe.base_element(base).n_dofs_per_face(face_no);
387  row_lengths[cell_indices[i]] += increment;
388  }
389  ++i;
390  }
391  // From now on, if an object is
392  // a cell, its dofs only couple
393  // inside the cell. Since the
394  // faces are handled below, we
395  // have to subtract ALL faces
396  // in this case.
397 
398  // In all other cases we
399  // subtract adjacent faces to be
400  // added in the loop below.
401  while (i < fe.get_first_quad_index(face_no))
402  {
403  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
404  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
405  ++mult)
406  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
407  fe.first_block_of_base(base) +
408  mult) != DoFTools::none)
409  {
410  increment =
411  fe.base_element(base).n_dofs_per_cell() -
412  ((dim > 1) ? (dim - 1) :
414  fe.base_element(base).n_dofs_per_face(face_no);
415  row_lengths[cell_indices[i]] += increment;
416  }
417  ++i;
418  }
419 
420  // Now quads in 2D and 3D
421  while (i < fe.get_first_hex_index())
422  {
423  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
424  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
425  ++mult)
426  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
427  fe.first_block_of_base(base) +
428  mult) != DoFTools::none)
429  {
430  increment =
431  fe.base_element(base).n_dofs_per_cell() -
432  ((dim > 2) ? (dim - 2) :
434  fe.base_element(base).n_dofs_per_face(face_no);
435  row_lengths[cell_indices[i]] += increment;
436  }
437  ++i;
438  }
439 
440  // Finally, cells in 3D
441  while (i < fe.n_dofs_per_cell())
442  {
443  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
444  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
445  ++mult)
446  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
447  fe.first_block_of_base(base) +
448  mult) != DoFTools::none)
449  {
450  increment =
451  fe.base_element(base).n_dofs_per_cell() -
453  fe.base_element(base).n_dofs_per_face(face_no);
454  row_lengths[cell_indices[i]] += increment;
455  }
456  ++i;
457  }
458 
459  // At this point, we have
460  // counted all dofs
461  // contributiong from cells
462  // coupled topologically to the
463  // adjacent cells, but we
464  // subtracted some faces.
465 
466  // Now, let's go by the faces
467  // and add the missing
468  // contribution as well as the
469  // flux contributions.
470  for (const unsigned int iface : GeometryInfo<dim>::face_indices())
471  {
472  bool level_boundary = cell->at_boundary(iface);
474  if (!level_boundary)
475  {
476  neighbor = cell->neighbor(iface);
477  if (static_cast<unsigned int>(neighbor->level()) != level)
478  level_boundary = true;
479  }
480 
481  if (level_boundary)
482  {
483  for (unsigned int local_dof = 0;
484  local_dof < fe.n_dofs_per_cell();
485  ++local_dof)
486  row_lengths[cell_indices[local_dof]] +=
487  fe.n_dofs_per_face(face_no);
488  continue;
489  }
490 
491  const FiniteElement<dim> &nfe = neighbor->get_fe();
493  cell->face(iface);
494 
495  // Flux couplings are
496  // computed from both sides
497  // for simplicity.
498 
499  // The dofs on the common face
500  // will be handled below,
501  // therefore, we subtract them
502  // here.
503  for (unsigned int base = 0; base < nfe.n_base_elements(); ++base)
504  for (unsigned int mult = 0; mult < nfe.element_multiplicity(base);
505  ++mult)
506  for (unsigned int local_dof = 0;
507  local_dof < fe.n_dofs_per_cell();
508  ++local_dof)
509  if (couple_face[fe_index](
510  fe.system_to_block_index(local_dof).first,
511  nfe.first_block_of_base(base) + mult) != DoFTools::none)
512  {
513  const unsigned int dof_increment =
514  nfe.base_element(base).n_dofs_per_cell() -
515  nfe.base_element(base).n_dofs_per_face(face_no);
516  row_lengths[cell_indices[local_dof]] += dof_increment;
517  }
518 
519  // Do this only once per face and not on the hanging faces.
520  if (face_touched[face->index()])
521  continue;
522  face_touched[face->index()] = true;
523 
524  // At this point, we assume
525  // that each cell added its
526  // dofs minus the face to
527  // the couplings of the
528  // face dofs. Since we
529  // subtracted two faces, we
530  // have to re-add one.
531 
532  // If one side of the face
533  // is refined, all the fine
534  // face dofs couple with
535  // the coarse one.
536 
537  // Wolfgang, do they couple
538  // with each other by
539  // constraints?
540 
541  // This will not work with
542  // different couplings on
543  // different cells.
544  neighbor_indices.resize(nfe.n_dofs_per_cell());
545  neighbor->get_mg_dof_indices(neighbor_indices);
546  for (unsigned int base = 0; base < nfe.n_base_elements(); ++base)
547  for (unsigned int mult = 0; mult < nfe.element_multiplicity(base);
548  ++mult)
549  for (unsigned int local_dof = 0;
550  local_dof < fe.n_dofs_per_cell();
551  ++local_dof)
552  if (couple_cell[fe_index](
553  fe.system_to_component_index(local_dof).first,
554  nfe.first_block_of_base(base) + mult) != DoFTools::none)
555  row_lengths[cell_indices[local_dof]] +=
556  nfe.base_element(base).n_dofs_per_face(face_no);
557  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
558  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
559  ++mult)
560  for (unsigned int local_dof = 0;
561  local_dof < nfe.n_dofs_per_cell();
562  ++local_dof)
563  if (couple_cell[fe_index](
564  nfe.system_to_component_index(local_dof).first,
565  fe.first_block_of_base(base) + mult) != DoFTools::none)
566  row_lengths[neighbor_indices[local_dof]] +=
567  fe.base_element(base).n_dofs_per_face(face_no);
568  }
569  }
570  }
571 
572 
573 
574  template <int dim, int spacedim, typename number>
575  void
577  SparsityPatternBase & sparsity,
578  const unsigned int level,
579  const AffineConstraints<number> &constraints,
580  const bool keep_constrained_dofs)
581  {
582  const types::global_dof_index n_dofs = dof.n_dofs(level);
583  (void)n_dofs;
584 
585  Assert(sparsity.n_rows() == n_dofs,
586  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
587  Assert(sparsity.n_cols() == n_dofs,
588  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
589 
590  const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
591  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
592  for (const auto &cell : dof.cell_iterators_on_level(level))
593  if (cell->is_locally_owned_on_level())
594  {
595  cell->get_mg_dof_indices(dofs_on_this_cell);
596  constraints.add_entries_local_to_global(dofs_on_this_cell,
597  sparsity,
598  keep_constrained_dofs);
599  }
600  }
601 
602 
603 
604  template <int dim, int spacedim, typename number>
605  void
607  SparsityPatternBase & sparsity,
608  const unsigned int level,
609  const AffineConstraints<number> &constraints,
610  const bool keep_constrained_dofs)
611  {
612  const types::global_dof_index n_dofs = dof.n_dofs(level);
613  (void)n_dofs;
614 
615  Assert(sparsity.n_rows() == n_dofs,
616  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
617  Assert(sparsity.n_cols() == n_dofs,
618  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
619 
620  const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
621  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
622  std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
623  for (const auto &cell : dof.cell_iterators_on_level(level))
624  {
625  if (!cell->is_locally_owned_on_level())
626  continue;
627 
628  cell->get_mg_dof_indices(dofs_on_this_cell);
629  // make sparsity pattern for this cell
630  constraints.add_entries_local_to_global(dofs_on_this_cell,
631  sparsity,
632  keep_constrained_dofs);
633 
634  // Loop over all interior neighbors
635  for (const unsigned int face : GeometryInfo<dim>::face_indices())
636  {
637  bool use_face = false;
638  if ((!cell->at_boundary(face)) &&
639  (static_cast<unsigned int>(cell->neighbor_level(face)) ==
640  level))
641  use_face = true;
642  else if (cell->has_periodic_neighbor(face) &&
643  (static_cast<unsigned int>(
644  cell->periodic_neighbor_level(face)) == level))
645  use_face = true;
646 
647  if (use_face)
648  {
649  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
650  cell->neighbor_or_periodic_neighbor(face);
651  neighbor->get_mg_dof_indices(dofs_on_other_cell);
652  // only add one direction The other is taken care of by
653  // neighbor (except when the neighbor is not owned by the same
654  // processor)
655  constraints.add_entries_local_to_global(dofs_on_this_cell,
656  dofs_on_other_cell,
657  sparsity,
658  keep_constrained_dofs);
659 
660  if (neighbor->is_locally_owned_on_level() == false)
661  {
662  constraints.add_entries_local_to_global(
663  dofs_on_other_cell, sparsity, keep_constrained_dofs);
664 
665  constraints.add_entries_local_to_global(
666  dofs_on_other_cell,
667  dofs_on_this_cell,
668  sparsity,
669  keep_constrained_dofs);
670  }
671  }
672  }
673  }
674  }
675 
676 
677 
678  template <int dim, int spacedim>
679  void
681  SparsityPatternBase & sparsity,
682  const unsigned int level)
683  {
684  Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
686 
687  const types::global_dof_index fine_dofs = dof.n_dofs(level);
688  const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
689  (void)fine_dofs;
690  (void)coarse_dofs;
691 
692  // Matrix maps from fine level to coarse level
693 
694  Assert(sparsity.n_rows() == coarse_dofs,
695  ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
696  Assert(sparsity.n_cols() == fine_dofs,
697  ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));
698 
699  const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
700  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
701  std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
702  for (const auto &cell : dof.cell_iterators_on_level(level))
703  {
704  if (!cell->is_locally_owned_on_level())
705  continue;
706 
707  cell->get_mg_dof_indices(dofs_on_this_cell);
708  // Loop over all interior neighbors
709  for (const unsigned int face : GeometryInfo<dim>::face_indices())
710  {
711  // Neighbor is coarser
712  bool use_face = false;
713  if ((!cell->at_boundary(face)) &&
714  (static_cast<unsigned int>(cell->neighbor_level(face)) !=
715  level))
716  use_face = true;
717  else if (cell->has_periodic_neighbor(face) &&
718  (static_cast<unsigned int>(
719  cell->periodic_neighbor_level(face)) != level))
720  use_face = true;
721 
722  if (use_face)
723  {
724  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
725  cell->neighbor_or_periodic_neighbor(face);
726  neighbor->get_mg_dof_indices(dofs_on_other_cell);
727 
728  std::sort(dofs_on_this_cell.begin(), dofs_on_this_cell.end());
729  for (unsigned int i = 0; i < dofs_per_cell; ++i)
730  sparsity.add_row_entries(dofs_on_other_cell[i],
731  make_array_view(dofs_on_this_cell),
732  true);
733  }
734  }
735  }
736  }
737 
738 
739 
740  template <int dim, int spacedim>
741  void
743  SparsityPatternBase & sparsity,
744  const unsigned int level,
745  const Table<2, DoFTools::Coupling> &int_mask,
746  const Table<2, DoFTools::Coupling> &flux_mask)
747  {
748  const FiniteElement<dim> & fe = dof.get_fe();
749  const types::global_dof_index n_dofs = dof.n_dofs(level);
750  const unsigned int n_comp = fe.n_components();
751  (void)n_dofs;
752  (void)n_comp;
753 
754  Assert(sparsity.n_rows() == n_dofs,
755  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
756  Assert(sparsity.n_cols() == n_dofs,
757  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
758  Assert(int_mask.n_rows() == n_comp,
759  ExcDimensionMismatch(int_mask.n_rows(), n_comp));
760  Assert(int_mask.n_cols() == n_comp,
761  ExcDimensionMismatch(int_mask.n_cols(), n_comp));
762  Assert(flux_mask.n_rows() == n_comp,
763  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
764  Assert(flux_mask.n_cols() == n_comp,
765  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
766 
767  const unsigned int total_dofs = fe.n_dofs_per_cell();
768  std::vector<types::global_dof_index> dofs_on_this_cell(total_dofs);
769  std::vector<types::global_dof_index> dofs_on_other_cell(total_dofs);
770  std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
771  cell_entries;
772 
773  Table<2, bool> support_on_face(total_dofs,
775 
777  int_dof_mask =
779  flux_dof_mask =
781 
782  for (unsigned int i = 0; i < total_dofs; ++i)
783  for (auto f : GeometryInfo<dim>::face_indices())
784  support_on_face(i, f) = fe.has_support_on_face(i, f);
785 
786  std::vector<bool> face_touched(dim == 2 ?
789 
790  for (const auto &cell : dof.cell_iterators_on_level(level))
791  {
792  if (!cell->is_locally_owned_on_level())
793  continue;
794 
795  cell->get_mg_dof_indices(dofs_on_this_cell);
796  // make sparsity pattern for this cell
797  for (unsigned int i = 0; i < total_dofs; ++i)
798  for (unsigned int j = 0; j < total_dofs; ++j)
799  if (int_dof_mask[i][j] != DoFTools::none)
800  cell_entries.emplace_back(dofs_on_this_cell[i],
801  dofs_on_this_cell[j]);
802 
803  // Loop over all interior neighbors
804  for (const unsigned int face : GeometryInfo<dim>::face_indices())
805  {
806  typename DoFHandler<dim, spacedim>::face_iterator cell_face =
807  cell->face(face);
808  if (face_touched[cell_face->index()])
809  continue;
810 
811  if (cell->at_boundary(face) && !cell->has_periodic_neighbor(face))
812  {
813  for (unsigned int i = 0; i < total_dofs; ++i)
814  {
815  const bool i_non_zero_i = support_on_face(i, face);
816  for (unsigned int j = 0; j < total_dofs; ++j)
817  {
818  const bool j_non_zero_i = support_on_face(j, face);
819 
820  if (flux_dof_mask(i, j) == DoFTools::always)
821  cell_entries.emplace_back(dofs_on_this_cell[i],
822  dofs_on_this_cell[j]);
823  if (flux_dof_mask(i, j) == DoFTools::nonzero &&
824  i_non_zero_i && j_non_zero_i)
825  cell_entries.emplace_back(dofs_on_this_cell[i],
826  dofs_on_this_cell[j]);
827  }
828  }
829  }
830  else
831  {
832  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
833  cell->neighbor_or_periodic_neighbor(face);
834 
835  if (neighbor->level() < cell->level())
836  continue;
837 
838  unsigned int neighbor_face =
839  cell->has_periodic_neighbor(face) ?
840  cell->periodic_neighbor_of_periodic_neighbor(face) :
841  cell->neighbor_of_neighbor(face);
842 
843  neighbor->get_mg_dof_indices(dofs_on_other_cell);
844  for (unsigned int i = 0; i < total_dofs; ++i)
845  {
846  const bool i_non_zero_i = support_on_face(i, face);
847  const bool i_non_zero_e = support_on_face(i, neighbor_face);
848  for (unsigned int j = 0; j < total_dofs; ++j)
849  {
850  const bool j_non_zero_i = support_on_face(j, face);
851  const bool j_non_zero_e =
852  support_on_face(j, neighbor_face);
853  if (flux_dof_mask(i, j) == DoFTools::always)
854  {
855  cell_entries.emplace_back(dofs_on_this_cell[i],
856  dofs_on_other_cell[j]);
857  cell_entries.emplace_back(dofs_on_other_cell[i],
858  dofs_on_this_cell[j]);
859  cell_entries.emplace_back(dofs_on_this_cell[i],
860  dofs_on_this_cell[j]);
861  cell_entries.emplace_back(dofs_on_other_cell[i],
862  dofs_on_other_cell[j]);
863  }
864  if (flux_dof_mask(i, j) == DoFTools::nonzero)
865  {
866  if (i_non_zero_i && j_non_zero_e)
867  cell_entries.emplace_back(dofs_on_this_cell[i],
868  dofs_on_other_cell[j]);
869  if (i_non_zero_e && j_non_zero_i)
870  cell_entries.emplace_back(dofs_on_other_cell[i],
871  dofs_on_this_cell[j]);
872  if (i_non_zero_i && j_non_zero_i)
873  cell_entries.emplace_back(dofs_on_this_cell[i],
874  dofs_on_this_cell[j]);
875  if (i_non_zero_e && j_non_zero_e)
876  cell_entries.emplace_back(dofs_on_other_cell[i],
877  dofs_on_other_cell[j]);
878  }
879 
880  if (flux_dof_mask(j, i) == DoFTools::always)
881  {
882  cell_entries.emplace_back(dofs_on_this_cell[j],
883  dofs_on_other_cell[i]);
884  cell_entries.emplace_back(dofs_on_other_cell[j],
885  dofs_on_this_cell[i]);
886  cell_entries.emplace_back(dofs_on_this_cell[j],
887  dofs_on_this_cell[i]);
888  cell_entries.emplace_back(dofs_on_other_cell[j],
889  dofs_on_other_cell[i]);
890  }
891  if (flux_dof_mask(j, i) == DoFTools::nonzero)
892  {
893  if (j_non_zero_i && i_non_zero_e)
894  cell_entries.emplace_back(dofs_on_this_cell[j],
895  dofs_on_other_cell[i]);
896  if (j_non_zero_e && i_non_zero_i)
897  cell_entries.emplace_back(dofs_on_other_cell[j],
898  dofs_on_this_cell[i]);
899  if (j_non_zero_i && i_non_zero_i)
900  cell_entries.emplace_back(dofs_on_this_cell[j],
901  dofs_on_this_cell[i]);
902  if (j_non_zero_e && i_non_zero_e)
903  cell_entries.emplace_back(dofs_on_other_cell[j],
904  dofs_on_other_cell[i]);
905  }
906  }
907  }
908  face_touched[neighbor->face(neighbor_face)->index()] = true;
909  }
910  }
911  sparsity.add_entries(make_array_view(cell_entries));
912  cell_entries.clear();
913  }
914  }
915 
916 
917 
918  template <int dim, int spacedim>
919  void
921  SparsityPatternBase & sparsity,
922  const unsigned int level,
923  const Table<2, DoFTools::Coupling> &flux_mask)
924  {
925  const FiniteElement<dim> &fe = dof.get_fe();
926  const unsigned int n_comp = fe.n_components();
927  (void)n_comp;
928 
929  Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
931 
932  const types::global_dof_index fine_dofs = dof.n_dofs(level);
933  const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
934  (void)fine_dofs;
935  (void)coarse_dofs;
936 
937  // Matrix maps from fine level to coarse level
938 
939  Assert(sparsity.n_rows() == coarse_dofs,
940  ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
941  Assert(sparsity.n_cols() == fine_dofs,
942  ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));
943  Assert(flux_mask.n_rows() == n_comp,
944  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
945  Assert(flux_mask.n_cols() == n_comp,
946  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
947 
948  const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
949  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
950  std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
951  std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
952  cell_entries;
953 
954  Table<2, bool> support_on_face(dofs_per_cell,
956 
957  const Table<2, DoFTools::Coupling> flux_dof_mask =
959 
960  for (unsigned int i = 0; i < dofs_per_cell; ++i)
961  for (auto f : GeometryInfo<dim>::face_indices())
962  support_on_face(i, f) = fe.has_support_on_face(i, f);
963 
964  for (const auto &cell : dof.cell_iterators_on_level(level))
965  {
966  if (!cell->is_locally_owned_on_level())
967  continue;
968 
969  cell->get_mg_dof_indices(dofs_on_this_cell);
970  // Loop over all interior neighbors
971  for (const unsigned int face : GeometryInfo<dim>::face_indices())
972  {
973  // Neighbor is coarser
974  bool use_face = false;
975  if ((!cell->at_boundary(face)) &&
976  (static_cast<unsigned int>(cell->neighbor_level(face)) !=
977  level))
978  use_face = true;
979  else if (cell->has_periodic_neighbor(face) &&
980  (static_cast<unsigned int>(
981  cell->periodic_neighbor_level(face)) != level))
982  use_face = true;
983 
984  if (use_face)
985  {
986  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
987  cell->neighbor_or_periodic_neighbor(face);
988  neighbor->get_mg_dof_indices(dofs_on_other_cell);
989 
990  for (unsigned int i = 0; i < dofs_per_cell; ++i)
991  {
992  for (unsigned int j = 0; j < dofs_per_cell; ++j)
993  {
994  if (flux_dof_mask(i, j) != DoFTools::none)
995  {
996  cell_entries.emplace_back(dofs_on_other_cell[i],
997  dofs_on_this_cell[j]);
998  cell_entries.emplace_back(dofs_on_other_cell[j],
999  dofs_on_this_cell[i]);
1000  }
1001  }
1002  }
1003  }
1004  }
1005  sparsity.add_entries(make_array_view(cell_entries));
1006  cell_entries.clear();
1007  }
1008  }
1009 
1010 
1011 
1012  template <int dim, int spacedim>
1013  void
1015  const MGConstrainedDoFs &mg_constrained_dofs,
1016  SparsityPatternBase & sparsity,
1017  const unsigned int level)
1018  {
1019  const types::global_dof_index n_dofs = dof.n_dofs(level);
1020  (void)n_dofs;
1021 
1022  Assert(sparsity.n_rows() == n_dofs,
1023  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1024  Assert(sparsity.n_cols() == n_dofs,
1025  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1026 
1027  const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
1028  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
1029  std::vector<types::global_dof_index> cols;
1030  cols.reserve(dofs_per_cell);
1031 
1032  for (const auto &cell : dof.cell_iterators_on_level(level))
1033  if (cell->is_locally_owned_on_level())
1034  {
1035  cell->get_mg_dof_indices(dofs_on_this_cell);
1036  std::sort(dofs_on_this_cell.begin(), dofs_on_this_cell.end());
1037  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1038  {
1039  for (unsigned int j = 0; j < dofs_per_cell; ++j)
1040  if (mg_constrained_dofs.is_interface_matrix_entry(
1041  level, dofs_on_this_cell[i], dofs_on_this_cell[j]))
1042  cols.push_back(dofs_on_this_cell[j]);
1043  sparsity.add_row_entries(dofs_on_this_cell[i],
1044  make_array_view(cols),
1045  true);
1046  cols.clear();
1047  }
1048  }
1049  }
1050 
1051 
1052 
1053  template <int dim, int spacedim>
1054  void
1056  const DoFHandler<dim, spacedim> & dof_handler,
1057  std::vector<std::vector<types::global_dof_index>> &result,
1058  bool only_once,
1059  std::vector<unsigned int> target_component)
1060  {
1061  const FiniteElement<dim> &fe = dof_handler.get_fe();
1062  const unsigned int n_components = fe.n_components();
1063  const unsigned int nlevels =
1064  dof_handler.get_triangulation().n_global_levels();
1065 
1066  Assert(result.size() == nlevels,
1067  ExcDimensionMismatch(result.size(), nlevels));
1068 
1069  if (target_component.size() == 0)
1070  {
1071  target_component.resize(n_components);
1072  for (unsigned int i = 0; i < n_components; ++i)
1073  target_component[i] = i;
1074  }
1075 
1076  Assert(target_component.size() == n_components,
1077  ExcDimensionMismatch(target_component.size(), n_components));
1078 
1079  for (unsigned int l = 0; l < nlevels; ++l)
1080  {
1081  result[l].resize(n_components);
1082  std::fill(result[l].begin(), result[l].end(), 0U);
1083 
1084  // special case for only one
1085  // component. treat this first
1086  // since it does not require any
1087  // computations
1088  if (n_components == 1)
1089  {
1090  result[l][0] = dof_handler.n_dofs(l);
1091  }
1092  else
1093  {
1094  // otherwise determine the number
1095  // of dofs in each component
1096  // separately. do so in parallel
1097  std::vector<std::vector<bool>> dofs_in_component(
1098  n_components, std::vector<bool>(dof_handler.n_dofs(l), false));
1099  std::vector<ComponentMask> component_select(n_components);
1100  Threads::TaskGroup<> tasks;
1101  for (unsigned int i = 0; i < n_components; ++i)
1102  {
1103  void (*fun_ptr)(const unsigned int level,
1104  const DoFHandler<dim, spacedim> &,
1105  const ComponentMask &,
1106  std::vector<bool> &) =
1107  &DoFTools::extract_level_dofs<dim, spacedim>;
1108 
1109  std::vector<bool> tmp(n_components, false);
1110  tmp[i] = true;
1111  component_select[i] = ComponentMask(tmp);
1112 
1113  tasks += Threads::new_task(fun_ptr,
1114  l,
1115  dof_handler,
1116  component_select[i],
1117  dofs_in_component[i]);
1118  }
1119  tasks.join_all();
1120 
1121  // next count what we got
1122  unsigned int component = 0;
1123  for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
1124  {
1125  const FiniteElement<dim> &base = fe.base_element(b);
1126  // Dimension of base element
1127  unsigned int d = base.n_components();
1128 
1129  for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
1130  {
1131  for (unsigned int dd = 0; dd < d; ++dd)
1132  {
1133  if (base.is_primitive() || (!only_once || dd == 0))
1134  result[l][target_component[component]] +=
1135  std::count(dofs_in_component[component].begin(),
1136  dofs_in_component[component].end(),
1137  true);
1138  ++component;
1139  }
1140  }
1141  }
1142  // finally sanity check
1143  Assert(!dof_handler.get_fe().is_primitive() ||
1144  std::accumulate(result[l].begin(),
1145  result[l].end(),
1147  dof_handler.n_dofs(l),
1148  ExcInternalError());
1149  }
1150  }
1151  }
1152 
1153 
1154 
1155  template <int dim, int spacedim>
1156  void
1158  const DoFHandler<dim, spacedim> & dof_handler,
1159  std::vector<std::vector<types::global_dof_index>> &dofs_per_block,
1160  std::vector<unsigned int> target_block)
1161  {
1162  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
1163  const unsigned int n_blocks = fe.n_blocks();
1164  const unsigned int n_levels =
1165  dof_handler.get_triangulation().n_global_levels();
1166 
1167  AssertDimension(dofs_per_block.size(), n_levels);
1168 
1169  for (unsigned int l = 0; l < n_levels; ++l)
1170  std::fill(dofs_per_block[l].begin(), dofs_per_block[l].end(), 0U);
1171  // If the empty vector was given as
1172  // default argument, set up this
1173  // vector as identity.
1174  if (target_block.size() == 0)
1175  {
1176  target_block.resize(n_blocks);
1177  for (unsigned int i = 0; i < n_blocks; ++i)
1178  target_block[i] = i;
1179  }
1180  Assert(target_block.size() == n_blocks,
1181  ExcDimensionMismatch(target_block.size(), n_blocks));
1182 
1183  const unsigned int max_block =
1184  *std::max_element(target_block.begin(), target_block.end());
1185  const unsigned int n_target_blocks = max_block + 1;
1186  (void)n_target_blocks;
1187 
1188  for (unsigned int l = 0; l < n_levels; ++l)
1189  AssertDimension(dofs_per_block[l].size(), n_target_blocks);
1190 
1191  // special case for only one
1192  // block. treat this first
1193  // since it does not require any
1194  // computations
1195  if (n_blocks == 1)
1196  {
1197  for (unsigned int l = 0; l < n_levels; ++l)
1198  dofs_per_block[l][0] = dof_handler.n_dofs(l);
1199  return;
1200  }
1201  // otherwise determine the number
1202  // of dofs in each block
1203  // separately. do so in parallel
1204  for (unsigned int l = 0; l < n_levels; ++l)
1205  {
1206  std::vector<std::vector<bool>> dofs_in_block(
1207  n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
1208  std::vector<BlockMask> block_select(n_blocks);
1209  Threads::TaskGroup<> tasks;
1210  for (unsigned int i = 0; i < n_blocks; ++i)
1211  {
1212  void (*fun_ptr)(const unsigned int level,
1213  const DoFHandler<dim, spacedim> &,
1214  const BlockMask &,
1215  std::vector<bool> &) =
1216  &DoFTools::extract_level_dofs<dim, spacedim>;
1217 
1218  std::vector<bool> tmp(n_blocks, false);
1219  tmp[i] = true;
1220  block_select[i] = tmp;
1221 
1222  tasks += Threads::new_task(
1223  fun_ptr, l, dof_handler, block_select[i], dofs_in_block[i]);
1224  }
1225  tasks.join_all();
1226 
1227  // next count what we got
1228  for (unsigned int block = 0; block < fe.n_blocks(); ++block)
1229  dofs_per_block[l][target_block[block]] +=
1230  std::count(dofs_in_block[block].begin(),
1231  dofs_in_block[block].end(),
1232  true);
1233  }
1234  }
1235 
1236 
1237 
1238  template <int dim, int spacedim>
1239  void
1241  const DoFHandler<dim, spacedim> &dof,
1242  const std::map<types::boundary_id, const Function<spacedim> *>
1243  & function_map,
1244  std::vector<std::set<types::global_dof_index>> &boundary_indices,
1245  const ComponentMask & component_mask)
1246  {
1247  Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1248  ExcDimensionMismatch(boundary_indices.size(),
1250 
1251  std::set<types::boundary_id> boundary_ids;
1252  for (const auto &boundary_function : function_map)
1253  boundary_ids.insert(boundary_function.first);
1254 
1255  std::vector<IndexSet> boundary_indexset;
1256  make_boundary_list(dof, boundary_ids, boundary_indexset, component_mask);
1257  for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1258  boundary_indices[i].insert(boundary_indexset[i].begin(),
1259  boundary_indexset[i].end());
1260  }
1261 
1262 
1263  template <int dim, int spacedim>
1264  void
1266  const std::map<types::boundary_id,
1267  const Function<spacedim> *> &function_map,
1268  std::vector<IndexSet> &boundary_indices,
1269  const ComponentMask & component_mask)
1270  {
1271  Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1272  ExcDimensionMismatch(boundary_indices.size(),
1274 
1275  std::set<types::boundary_id> boundary_ids;
1276  for (const auto &boundary_function : function_map)
1277  boundary_ids.insert(boundary_function.first);
1278 
1279  make_boundary_list(dof, boundary_ids, boundary_indices, component_mask);
1280  }
1281 
1282 
1283 
1284  template <int dim, int spacedim>
1285  void
1287  const std::set<types::boundary_id> &boundary_ids,
1288  std::vector<IndexSet> & boundary_indices,
1289  const ComponentMask & component_mask)
1290  {
1291  boundary_indices.resize(dof.get_triangulation().n_global_levels());
1292 
1293  // if for whatever reason we were passed an empty set, return immediately
1294  if (boundary_ids.size() == 0)
1295  return;
1296 
1297  for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1298  if (boundary_indices[i].size() == 0)
1299  boundary_indices[i] = IndexSet(dof.n_dofs(i));
1300 
1301  const unsigned int n_components = dof.get_fe_collection().n_components();
1302  const bool fe_is_system = (n_components != 1);
1303 
1304  std::vector<types::global_dof_index> local_dofs;
1305  local_dofs.reserve(dof.get_fe_collection().max_dofs_per_face());
1306  std::fill(local_dofs.begin(), local_dofs.end(), numbers::invalid_dof_index);
1307 
1308  std::vector<std::vector<types::global_dof_index>> dofs_by_level(
1309  dof.get_triangulation().n_levels());
1310 
1311  // First, deal with the simpler case when we have to identify all boundary
1312  // dofs
1313  if (component_mask.n_selected_components(n_components) == n_components)
1314  {
1315  for (const auto &cell : dof.cell_iterators())
1316  {
1317  if (cell->is_artificial_on_level())
1318  continue;
1319  const FiniteElement<dim> &fe = cell->get_fe();
1320  const unsigned int level = cell->level();
1321 
1322  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
1323  if (cell->at_boundary(face_no) == true)
1324  {
1325  const typename DoFHandler<dim, spacedim>::face_iterator face =
1326  cell->face(face_no);
1327  const types::boundary_id bi = face->boundary_id();
1328  // Face is listed in boundary map
1329  if (boundary_ids.find(bi) != boundary_ids.end())
1330  {
1331  local_dofs.resize(fe.n_dofs_per_face(face_no));
1332  face->get_mg_dof_indices(level, local_dofs);
1333  dofs_by_level[level].insert(dofs_by_level[level].end(),
1334  local_dofs.begin(),
1335  local_dofs.end());
1336  }
1337  }
1338  }
1339  }
1340  else
1341  {
1342  Assert(component_mask.n_selected_components(n_components) > 0,
1343  ExcMessage(
1344  "It's probably worthwhile to select at least one component."));
1345 
1346  for (const auto &cell : dof.cell_iterators())
1347  if (!cell->is_artificial_on_level())
1348  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
1349  {
1350  if (cell->at_boundary(face_no) == false)
1351  continue;
1352 
1353  const FiniteElement<dim> &fe = cell->get_fe();
1354  const unsigned int level = cell->level();
1355 
1357  cell->face(face_no);
1358  const types::boundary_id boundary_component =
1359  face->boundary_id();
1360  if (boundary_ids.find(boundary_component) != boundary_ids.end())
1361  // we want to constrain this boundary
1362  {
1363  for (unsigned int i = 0;
1364  i < cell->get_fe().n_dofs_per_cell();
1365  ++i)
1366  {
1367  const ComponentMask &nonzero_component_array =
1368  cell->get_fe().get_nonzero_components(i);
1369  // if we want to constrain one of the nonzero
1370  // components, we have to constrain all of them
1371 
1372  bool selected = false;
1373  for (unsigned int c = 0; c < n_components; ++c)
1374  if (nonzero_component_array[c] == true &&
1375  component_mask[c] == true)
1376  {
1377  selected = true;
1378  break;
1379  }
1380  if (selected)
1381  for (unsigned int c = 0; c < n_components; ++c)
1382  Assert(
1383  nonzero_component_array[c] == false ||
1384  component_mask[c] == true,
1385  ExcMessage(
1386  "You are using a non-primitive FiniteElement "
1387  "and try to constrain just some of its components!"));
1388  }
1389 
1390  // get indices, physical location and boundary values of
1391  // dofs on this face
1392  local_dofs.resize(fe.n_dofs_per_face(face_no));
1393  face->get_mg_dof_indices(level, local_dofs);
1394  if (fe_is_system)
1395  {
1396  for (unsigned int i = 0; i < local_dofs.size(); ++i)
1397  {
1398  unsigned int component =
1400  if (fe.is_primitive())
1401  component =
1402  fe.face_system_to_component_index(i, face_no)
1403  .first;
1404  else
1405  {
1406  // Just pick the first of the components
1407  // We already know that either all or none
1408  // of the components are selected
1409  const ComponentMask &nonzero_component_array =
1410  cell->get_fe().get_nonzero_components(i);
1411  for (unsigned int c = 0; c < n_components; ++c)
1412  if (nonzero_component_array[c] == true)
1413  {
1414  component = c;
1415  break;
1416  }
1417  }
1419  ExcInternalError());
1420  if (component_mask[component] == true)
1421  dofs_by_level[level].push_back(local_dofs[i]);
1422  }
1423  }
1424  else
1425  dofs_by_level[level].insert(dofs_by_level[level].end(),
1426  local_dofs.begin(),
1427  local_dofs.end());
1428  }
1429  }
1430  }
1431  for (unsigned int level = 0; level < dof.get_triangulation().n_levels();
1432  ++level)
1433  {
1434  std::sort(dofs_by_level[level].begin(), dofs_by_level[level].end());
1435  boundary_indices[level].add_indices(
1436  dofs_by_level[level].begin(),
1437  std::unique(dofs_by_level[level].begin(),
1438  dofs_by_level[level].end()));
1439  }
1440  }
1441 
1442 
1443 
1444  template <int dim, int spacedim>
1445  void
1447  std::vector<IndexSet> & interface_dofs)
1448  {
1449  Assert(interface_dofs.size() ==
1450  mg_dof_handler.get_triangulation().n_global_levels(),
1452  interface_dofs.size(),
1453  mg_dof_handler.get_triangulation().n_global_levels()));
1454 
1455  std::vector<std::vector<types::global_dof_index>> tmp_interface_dofs(
1456  interface_dofs.size());
1457 
1458  const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
1459 
1460  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1461 
1462  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1463 
1464  std::vector<bool> cell_dofs(dofs_per_cell, false);
1465 
1466  for (const auto &cell : mg_dof_handler.cell_iterators())
1467  {
1468  // Do not look at artificial level cells (in a serial computation we
1469  // need to ignore the level_subdomain_id() because it is never set).
1470  if (cell->is_artificial_on_level())
1471  continue;
1472 
1473  bool has_coarser_neighbor = false;
1474 
1475  std::fill(cell_dofs.begin(), cell_dofs.end(), false);
1476 
1477  for (const unsigned int face_nr : GeometryInfo<dim>::face_indices())
1478  {
1479  const typename DoFHandler<dim, spacedim>::face_iterator face =
1480  cell->face(face_nr);
1481  if (!face->at_boundary() || cell->has_periodic_neighbor(face_nr))
1482  {
1483  // interior face
1484  const typename DoFHandler<dim>::cell_iterator neighbor =
1485  cell->neighbor_or_periodic_neighbor(face_nr);
1486 
1487  // only process cell pairs if one or both of them are owned by
1488  // me (ignore if running in serial)
1489  if (neighbor->is_artificial_on_level())
1490  continue;
1491 
1492  // Do refinement face from the coarse side
1493  if (neighbor->level() < cell->level())
1494  {
1495  for (unsigned int j = 0; j < fe.n_dofs_per_face(face_nr);
1496  ++j)
1497  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1498 
1499  has_coarser_neighbor = true;
1500  }
1501  }
1502  }
1503 
1504  if (has_coarser_neighbor == false)
1505  continue;
1506 
1507  const unsigned int level = cell->level();
1508  cell->get_mg_dof_indices(local_dof_indices);
1509 
1510  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1511  {
1512  if (cell_dofs[i])
1513  tmp_interface_dofs[level].push_back(local_dof_indices[i]);
1514  }
1515  }
1516 
1517  for (unsigned int l = 0;
1518  l < mg_dof_handler.get_triangulation().n_global_levels();
1519  ++l)
1520  {
1521  interface_dofs[l].clear();
1522  std::sort(tmp_interface_dofs[l].begin(), tmp_interface_dofs[l].end());
1523  interface_dofs[l].add_indices(tmp_interface_dofs[l].begin(),
1524  std::unique(tmp_interface_dofs[l].begin(),
1525  tmp_interface_dofs[l].end()));
1526  interface_dofs[l].compress();
1527  }
1528  }
1529 
1530 
1531 
1532  template <int dim, int spacedim>
1533  unsigned int
1535  {
1536  // Find minimum level for an active cell in
1537  // this locally owned subdomain
1538  // Note: with the way active cells are traversed,
1539  // the first locally owned cell we find will have
1540  // the lowest level in the particular subdomain.
1541  unsigned int min_level = tria.n_global_levels();
1542  for (const auto &cell : tria.active_cell_iterators())
1543  if (cell->is_locally_owned())
1544  {
1545  min_level = cell->level();
1546  break;
1547  }
1548 
1549  unsigned int global_min = min_level;
1550  // If necessary, communicate to find minimum
1551  // level for an active cell over all subdomains
1553  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1554  &tria))
1555  global_min = Utilities::MPI::min(min_level, tr->get_communicator());
1556 
1557  AssertIndexRange(global_min, tria.n_global_levels());
1558 
1559  return global_min;
1560  }
1561 
1562 
1563 
1564  namespace internal
1565  {
1566  double
1568  const std::vector<types::global_dof_index> &n_cells_on_levels,
1569  const MPI_Comm comm)
1570  {
1571  std::vector<types::global_dof_index> n_cells_on_leves_max(
1572  n_cells_on_levels.size());
1573  std::vector<types::global_dof_index> n_cells_on_leves_sum(
1574  n_cells_on_levels.size());
1575 
1576  Utilities::MPI::max(n_cells_on_levels, comm, n_cells_on_leves_max);
1577  Utilities::MPI::sum(n_cells_on_levels, comm, n_cells_on_leves_sum);
1578 
1579  const unsigned int n_proc = Utilities::MPI::n_mpi_processes(comm);
1580 
1581  const double ideal_work = std::accumulate(n_cells_on_leves_sum.begin(),
1582  n_cells_on_leves_sum.end(),
1583  0) /
1584  static_cast<double>(n_proc);
1585  const double workload_imbalance =
1586  std::accumulate(n_cells_on_leves_max.begin(),
1587  n_cells_on_leves_max.end(),
1588  0) /
1589  ideal_work;
1590 
1591  return workload_imbalance;
1592  }
1593 
1594 
1595 
1596  double
1598  const std::vector<
1599  std::pair<types::global_dof_index, types::global_dof_index>> &cells,
1600  const MPI_Comm comm)
1601  {
1602  std::vector<types::global_dof_index> cells_local(cells.size());
1603  std::vector<types::global_dof_index> cells_remote(cells.size());
1604 
1605  for (unsigned int i = 0; i < cells.size(); ++i)
1606  {
1607  cells_local[i] = cells[i].first;
1608  cells_remote[i] = cells[i].second;
1609  }
1610 
1611  std::vector<types::global_dof_index> cells_local_sum(cells_local.size());
1612  Utilities::MPI::sum(cells_local, comm, cells_local_sum);
1613 
1614  std::vector<types::global_dof_index> cells_remote_sum(
1615  cells_remote.size());
1616  Utilities::MPI::sum(cells_remote, comm, cells_remote_sum);
1617 
1618  const auto n_cells_local =
1619  std::accumulate(cells_local_sum.begin(), cells_local_sum.end(), 0);
1620  const auto n_cells_remote =
1621  std::accumulate(cells_remote_sum.begin(), cells_remote_sum.end(), 0);
1622 
1623  return static_cast<double>(n_cells_local) /
1624  (n_cells_local + n_cells_remote);
1625  }
1626  } // namespace internal
1627 
1628 
1629 
1630  template <int dim, int spacedim>
1631  std::vector<types::global_dof_index>
1633  {
1635  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1636  &tria))
1637  Assert(
1638  tr->is_multilevel_hierarchy_constructed(),
1639  ExcMessage(
1640  "We can only compute the workload imbalance if the multilevel hierarchy has been constructed!"));
1641 
1642  const unsigned int n_global_levels = tria.n_global_levels();
1643 
1644  std::vector<types::global_dof_index> n_cells_on_levels(n_global_levels);
1645 
1646  for (unsigned int lvl = 0; lvl < n_global_levels; ++lvl)
1647  for (const auto &cell : tria.cell_iterators_on_level(lvl))
1648  if (cell->is_locally_owned_on_level())
1649  ++n_cells_on_levels[lvl];
1650 
1651  return n_cells_on_levels;
1652  }
1653 
1654 
1655 
1656  template <int dim, int spacedim>
1657  std::vector<types::global_dof_index>
1659  const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1660  &trias)
1661  {
1662  const unsigned int n_global_levels = trias.size();
1663 
1664  std::vector<types::global_dof_index> n_cells_on_levels(n_global_levels);
1665 
1666  for (unsigned int lvl = 0; lvl < n_global_levels; ++lvl)
1667  for (const auto &cell : trias[lvl]->active_cell_iterators())
1668  if (cell->is_locally_owned())
1669  ++n_cells_on_levels[lvl];
1670 
1671  return n_cells_on_levels;
1672  }
1673 
1674 
1675 
1676  template <int dim, int spacedim>
1677  double
1679  {
1682  }
1683 
1684 
1685 
1686  template <int dim, int spacedim>
1687  double
1689  const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1690  &trias)
1691  {
1693  trias.back()->get_communicator());
1694  }
1695 
1696 
1697 
1698  template <int dim, int spacedim>
1699  std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1701  {
1702  const unsigned int n_global_levels = tria.n_global_levels();
1703 
1704  std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1705  cells(n_global_levels);
1706 
1707  const MPI_Comm communicator = tria.get_communicator();
1708 
1709  const unsigned int my_rank = Utilities::MPI::this_mpi_process(communicator);
1710 
1711  for (unsigned int lvl = 0; lvl < n_global_levels - 1; ++lvl)
1712  for (const auto &cell : tria.cell_iterators_on_level(lvl))
1713  if (cell->is_locally_owned_on_level() && cell->has_children())
1714  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1715  ++i)
1716  {
1717  const auto level_subdomain_id =
1718  cell->child(i)->level_subdomain_id();
1719  if (level_subdomain_id == my_rank)
1720  ++cells[lvl + 1].first;
1721  else if (level_subdomain_id != numbers::invalid_unsigned_int)
1722  ++cells[lvl + 1].second;
1723  else
1724  AssertThrow(false, ExcNotImplemented());
1725  }
1726 
1727  return cells;
1728  }
1729 
1730 
1731 
1732  template <int dim, int spacedim>
1733  std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1735  const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1736  &trias)
1737  {
1738  const unsigned int n_global_levels = trias.size();
1739 
1740  std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1741  cells(n_global_levels);
1742 
1743  const MPI_Comm communicator = trias.back()->get_communicator();
1744 
1745  const unsigned int my_rank = Utilities::MPI::this_mpi_process(communicator);
1746 
1747  for (unsigned int lvl = 0; lvl < n_global_levels - 1; ++lvl)
1748  {
1749  const auto &tria_coarse = *trias[lvl];
1750  const auto &tria_fine = *trias[lvl + 1];
1751 
1752  const unsigned int n_coarse_cells = tria_fine.n_global_coarse_cells();
1753  const unsigned int n_global_levels = tria_fine.n_global_levels();
1754 
1755  const ::internal::CellIDTranslator<dim> cell_id_translator(
1756  n_coarse_cells, n_global_levels);
1757 
1758  IndexSet is_fine_owned(cell_id_translator.size());
1759  IndexSet is_fine_required(cell_id_translator.size());
1760 
1761  for (const auto &cell : tria_fine.active_cell_iterators())
1762  if (!cell->is_artificial() && cell->is_locally_owned())
1763  is_fine_owned.add_index(cell_id_translator.translate(cell));
1764 
1765  for (const auto &cell : tria_coarse.active_cell_iterators())
1766  if (!cell->is_artificial() && cell->is_locally_owned())
1767  {
1768  if (cell->level() + 1u == tria_fine.n_global_levels())
1769  continue;
1770 
1771  for (unsigned int i = 0;
1772  i < GeometryInfo<dim>::max_children_per_cell;
1773  ++i)
1774  is_fine_required.add_index(
1775  cell_id_translator.translate(cell, i));
1776  }
1777 
1778  std::vector<unsigned int> is_fine_required_ranks(
1779  is_fine_required.n_elements());
1780 
1782  process(is_fine_owned,
1783  is_fine_required,
1784  communicator,
1785  is_fine_required_ranks,
1786  false);
1787 
1789  std::vector<
1790  std::pair<types::global_cell_index, types::global_cell_index>>,
1791  std::vector<unsigned int>>
1792  consensus_algorithm;
1793  consensus_algorithm.run(process, communicator);
1794 
1795  for (unsigned i = 0; i < is_fine_required.n_elements(); ++i)
1796  if (is_fine_required_ranks[i] == my_rank)
1797  ++cells[lvl + 1].first;
1798  else if (is_fine_required_ranks[i] != numbers::invalid_unsigned_int)
1799  ++cells[lvl + 1].second;
1800  }
1801 
1802  return cells;
1803  }
1804 
1805 
1806 
1807  template <int dim, int spacedim>
1808  double
1810  {
1813  }
1814 
1815 
1816 
1817  template <int dim, int spacedim>
1818  double
1820  const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1821  &trias)
1822  {
1825  trias.back()->get_communicator());
1826  }
1827 
1828 } // namespace MGTools
1829 
1830 
1831 // explicit instantiations
1832 #include "mg_tools.inst"
1833 
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:699
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
unsigned int get_first_line_index() const
unsigned int n_dofs_per_cell() const
unsigned int get_first_quad_index(const unsigned int quad_no=0) const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_blocks() const
unsigned int n_components() const
ReferenceCell reference_cell() const
unsigned int get_first_hex_index() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
bool is_primitive() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
unsigned int element_multiplicity(const unsigned int index) const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
unsigned int n_base_elements() const
types::global_dof_index first_block_of_base(const unsigned int b) const
size_type n_elements() const
Definition: index_set.h:1774
void add_index(const size_type index)
Definition: index_set.h:1645
bool is_interface_matrix_entry(const unsigned int level, const types::global_dof_index i, const types::global_dof_index j) const
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false)=0
size_type n_cols() const
virtual MPI_Comm get_communicator() const
unsigned int n_raw_lines() const
unsigned int n_levels() const
virtual unsigned int n_global_levels() const
unsigned int n_raw_quads() const
virtual std::vector< unsigned int > run(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm &comm) override
unsigned int max_dofs_per_face() const
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
unsigned int level
Definition: grid_out.cc:4608
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:465
Task< RT > new_task(const std::function< RT()> &function)
void convert_couplings_to_blocks(const DoFHandler< dim, spacedim > &dof_handler, const Table< 2, Coupling > &table_by_component, std::vector< Table< 2, Coupling >> &tables_by_block)
Definition: dof_tools.cc:2377
Table< 2, Coupling > dof_couplings_from_component_couplings(const FiniteElement< dim, spacedim > &fe, const Table< 2, Coupling > &component_couplings)
static const char U
double workload_imbalance(const std::vector< types::global_dof_index > &n_cells_on_levels, const MPI_Comm comm)
Definition: mg_tools.cc:1567
double vertical_communication_efficiency(const std::vector< std::pair< types::global_dof_index, types::global_dof_index >> &cells, const MPI_Comm comm)
Definition: mg_tools.cc:1597
std::vector< std::pair< types::global_dof_index, types::global_dof_index > > local_vertical_communication_cost(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1700
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true)
Definition: mg_tools.cc:576
void make_interface_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, const MGConstrainedDoFs &mg_constrained_dofs, SparsityPatternBase &sparsity, const unsigned int level)
Definition: mg_tools.cc:1014
std::vector< types::global_dof_index > local_workload(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1632
void make_flux_sparsity_pattern_edge(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level)
Definition: mg_tools.cc:680
void compute_row_length_vector(const DoFHandler< dim, spacedim > &dofs, const unsigned int level, std::vector< unsigned int > &row_lengths, const DoFTools::Coupling flux_couplings=DoFTools::none)
Definition: mg_tools.cc:107
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1446
void make_boundary_list(const DoFHandler< dim, spacedim > &mg_dof, const std::map< types::boundary_id, const Function< spacedim > * > &function_map, std::vector< std::set< types::global_dof_index >> &boundary_indices, const ComponentMask &component_mask=ComponentMask())
Definition: mg_tools.cc:1240
double workload_imbalance(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1678
void count_dofs_per_block(const DoFHandler< dim, spacedim > &dof_handler, std::vector< std::vector< types::global_dof_index >> &dofs_per_block, std::vector< unsigned int > target_block={})
Definition: mg_tools.cc:1157
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true)
Definition: mg_tools.cc:606
double vertical_communication_efficiency(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1809
void count_dofs_per_component(const DoFHandler< dim, spacedim > &mg_dof, std::vector< std::vector< types::global_dof_index >> &result, const bool only_once=false, std::vector< unsigned int > target_component={})
Definition: mg_tools.cc:1055
unsigned int max_level_for_coarse_mesh(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1534
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:157
T min(const T &t, const MPI_Comm &mpi_communicator)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:146
T max(const T &t, const MPI_Comm &mpi_communicator)
static const unsigned int invalid_unsigned_int
Definition: types.h:206
const types::global_dof_index invalid_dof_index
Definition: types.h:226
unsigned int global_dof_index
Definition: types.h:81
unsigned short int fe_index
Definition: types.h:59
unsigned int boundary_id
Definition: types.h:134
const MPI_Comm & comm
const ::Triangulation< dim, spacedim > & tria