Reference documentation for deal.II version 9.2.0
\(\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\}}\)
coupling.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2020 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/exceptions.h>
17 #include <deal.II/base/point.h>
18 #include <deal.II/base/tensor.h>
19 
20 #include <deal.II/distributed/shared_tria.h>
21 #include <deal.II/distributed/tria.h>
22 
23 #include <deal.II/fe/fe_values.h>
24 
25 #include <deal.II/grid/grid_tools.h>
26 #include <deal.II/grid/grid_tools_cache.h>
27 
28 #include <deal.II/lac/block_sparse_matrix.h>
29 #include <deal.II/lac/block_sparsity_pattern.h>
30 #include <deal.II/lac/petsc_block_sparse_matrix.h>
31 #include <deal.II/lac/petsc_sparse_matrix.h>
32 #include <deal.II/lac/sparse_matrix.h>
33 #include <deal.II/lac/sparsity_pattern.h>
34 #include <deal.II/lac/trilinos_block_sparse_matrix.h>
35 #include <deal.II/lac/trilinos_sparse_matrix.h>
36 #include <deal.II/lac/trilinos_sparsity_pattern.h>
37 
38 #include <deal.II/non_matching/coupling.h>
39 
40 DEAL_II_NAMESPACE_OPEN
41 namespace NonMatching
42 {
43  namespace internal
44  {
62  template <int dim0, int dim1, int spacedim>
63  std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>
64  qpoints_over_locally_owned_cells(
65  const GridTools::Cache<dim0, spacedim> &cache,
66  const DoFHandler<dim1, spacedim> & immersed_dh,
67  const Quadrature<dim1> & quad,
68  const Mapping<dim1, spacedim> & immersed_mapping,
69  const bool tria_is_parallel)
70  {
71  const auto & immersed_fe = immersed_dh.get_fe();
72  std::vector<Point<spacedim>> points_over_local_cells;
73  // Keep track of which cells we actually used
74  std::vector<unsigned int> used_cells_ids;
75  {
76  FEValues<dim1, spacedim> fe_v(immersed_mapping,
77  immersed_fe,
78  quad,
79  update_quadrature_points);
80  unsigned int cell_id = 0;
81  for (const auto &cell : immersed_dh.active_cell_iterators())
82  {
83  bool use_cell = false;
84  if (tria_is_parallel)
85  {
86  const auto bbox = cell->bounding_box();
87  std::vector<std::pair<
88  BoundingBox<spacedim>,
89  typename Triangulation<dim0, spacedim>::active_cell_iterator>>
90  out_vals;
91  cache.get_cell_bounding_boxes_rtree().query(
92  boost::geometry::index::intersects(bbox),
93  std::back_inserter(out_vals));
94  // Each bounding box corresponds to an active cell
95  // of the embedding triangulation: we now check if
96  // the current cell, of the embedded triangulation,
97  // overlaps a locally owned cell of the embedding one
98  for (const auto &bbox_it : out_vals)
99  if (bbox_it.second->is_locally_owned())
100  {
101  use_cell = true;
102  used_cells_ids.emplace_back(cell_id);
103  break;
104  }
105  }
106  else
107  // for sequential triangulations, simply use all cells
108  use_cell = true;
109 
110  if (use_cell)
111  {
112  // Reinitialize the cell and the fe_values
113  fe_v.reinit(cell);
114  const std::vector<Point<spacedim>> &x_points =
115  fe_v.get_quadrature_points();
116 
117  // Insert the points to the vector
118  points_over_local_cells.insert(points_over_local_cells.end(),
119  x_points.begin(),
120  x_points.end());
121  }
122  ++cell_id;
123  }
124  }
125  return {std::move(points_over_local_cells), std::move(used_cells_ids)};
126  }
127 
128 
135  template <int dim0, int dim1, int spacedim>
136  std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
137  compute_components_coupling(const ComponentMask & comps0,
138  const ComponentMask & comps1,
139  const FiniteElement<dim0, spacedim> &fe0,
140  const FiniteElement<dim1, spacedim> &fe1)
141  {
142  // Take care of components
143  const ComponentMask mask0 =
144  (comps0.size() == 0 ? ComponentMask(fe0.n_components(), true) : comps0);
145 
146  const ComponentMask mask1 =
147  (comps1.size() == 0 ? ComponentMask(fe1.n_components(), true) : comps1);
148 
149  AssertDimension(mask0.size(), fe0.n_components());
150  AssertDimension(mask1.size(), fe1.n_components());
151 
152  // Global to local indices
153  std::vector<unsigned int> gtl0(fe0.n_components(),
154  numbers::invalid_unsigned_int);
155  std::vector<unsigned int> gtl1(fe1.n_components(),
156  numbers::invalid_unsigned_int);
157 
158  for (unsigned int i = 0, j = 0; i < gtl0.size(); ++i)
159  if (mask0[i])
160  gtl0[i] = j++;
161 
162  for (unsigned int i = 0, j = 0; i < gtl1.size(); ++i)
163  if (mask1[i])
164  gtl1[i] = j++;
165  return {gtl0, gtl1};
166  }
167  } // namespace internal
168 
169  template <int dim0,
170  int dim1,
171  int spacedim,
172  typename Sparsity,
173  typename number>
174  void
175  create_coupling_sparsity_pattern(
176  const DoFHandler<dim0, spacedim> &space_dh,
177  const DoFHandler<dim1, spacedim> &immersed_dh,
178  const Quadrature<dim1> & quad,
179  Sparsity & sparsity,
180  const AffineConstraints<number> & constraints,
181  const ComponentMask & space_comps,
182  const ComponentMask & immersed_comps,
183  const Mapping<dim0, spacedim> & space_mapping,
184  const Mapping<dim1, spacedim> & immersed_mapping)
185  {
186  GridTools::Cache<dim0, spacedim> cache(space_dh.get_triangulation(),
187  space_mapping);
188  create_coupling_sparsity_pattern(cache,
189  space_dh,
190  immersed_dh,
191  quad,
192  sparsity,
193  constraints,
194  space_comps,
195  immersed_comps,
196  immersed_mapping);
197  }
198 
199 
200 
201  template <int dim0,
202  int dim1,
203  int spacedim,
204  typename Sparsity,
205  typename number>
206  void
207  create_coupling_sparsity_pattern(
208  const GridTools::Cache<dim0, spacedim> &cache,
209  const DoFHandler<dim0, spacedim> & space_dh,
210  const DoFHandler<dim1, spacedim> & immersed_dh,
211  const Quadrature<dim1> & quad,
212  Sparsity & sparsity,
213  const AffineConstraints<number> & constraints,
214  const ComponentMask & space_comps,
215  const ComponentMask & immersed_comps,
216  const Mapping<dim1, spacedim> & immersed_mapping)
217  {
218  AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
219  AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
220  Assert(dim1 <= dim0,
221  ExcMessage("This function can only work if dim1 <= dim0"));
222  Assert((dynamic_cast<
223  const parallel::distributed::Triangulation<dim1, spacedim> *>(
224  &immersed_dh.get_triangulation()) == nullptr),
225  ExcNotImplemented());
226 
227  const bool tria_is_parallel =
228  (dynamic_cast<const parallel::TriangulationBase<dim1, spacedim> *>(
229  &space_dh.get_triangulation()) != nullptr);
230  const auto &space_fe = space_dh.get_fe();
231  const auto &immersed_fe = immersed_dh.get_fe();
232 
233  // Dof indices
234  std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
235  std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
236 
237  // Take care of components
238  const ComponentMask space_c =
239  (space_comps.size() == 0 ? ComponentMask(space_fe.n_components(), true) :
240  space_comps);
241 
242  const ComponentMask immersed_c =
243  (immersed_comps.size() == 0 ?
244  ComponentMask(immersed_fe.n_components(), true) :
245  immersed_comps);
246 
247  AssertDimension(space_c.size(), space_fe.n_components());
248  AssertDimension(immersed_c.size(), immersed_fe.n_components());
249 
250  // Global to local indices
251  std::vector<unsigned int> space_gtl(space_fe.n_components(),
252  numbers::invalid_unsigned_int);
253  std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(),
254  numbers::invalid_unsigned_int);
255 
256  for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
257  if (space_c[i])
258  space_gtl[i] = j++;
259 
260  for (unsigned int i = 0, j = 0; i < immersed_gtl.size(); ++i)
261  if (immersed_c[i])
262  immersed_gtl[i] = j++;
263 
264  const unsigned int n_q_points = quad.size();
265  const unsigned int n_active_c =
266  immersed_dh.get_triangulation().n_active_cells();
267 
268  const auto qpoints_cells_data = internal::qpoints_over_locally_owned_cells(
269  cache, immersed_dh, quad, immersed_mapping, tria_is_parallel);
270 
271  const auto &points_over_local_cells = std::get<0>(qpoints_cells_data);
272  const auto &used_cells_ids = std::get<1>(qpoints_cells_data);
273 
274  // [TODO]: when the add_entries_local_to_global below will implement
275  // the version with the dof_mask, this should be uncommented.
276  //
277  // // Construct a dof_mask, used to distribute entries to the sparsity
278  // able< 2, bool > dof_mask(space_fe.dofs_per_cell,
279  // immersed_fe.dofs_per_cell);
280  // of_mask.fill(false);
281  // or (unsigned int i=0; i<space_fe.dofs_per_cell; ++i)
282  // {
283  // const auto comp_i = space_fe.system_to_component_index(i).first;
284  // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
285  // for (unsigned int j=0; j<immersed_fe.dofs_per_cell; ++j)
286  // {
287  // const auto comp_j =
288  // immersed_fe.system_to_component_index(j).first; if
289  // (immersed_gtl[comp_j] == space_gtl[comp_i])
290  // dof_mask(i,j) = true;
291  // }
292  // }
293 
294 
295  // Get a list of outer cells, qpoints and maps.
296  const auto cpm =
297  GridTools::compute_point_locations(cache, points_over_local_cells);
298  const auto &all_cells = std::get<0>(cpm);
299  const auto &maps = std::get<2>(cpm);
300 
301  std::vector<
302  std::set<typename Triangulation<dim0, spacedim>::active_cell_iterator>>
303  cell_sets(n_active_c);
304 
305  for (unsigned int i = 0; i < maps.size(); ++i)
306  {
307  // Quadrature points should be reasonably clustered:
308  // the following index keeps track of the last id
309  // where the current cell was inserted
310  unsigned int last_id = std::numeric_limits<unsigned int>::max();
311  unsigned int cell_id;
312  for (const unsigned int idx : maps[i])
313  {
314  // Find in which cell of immersed triangulation the point lies
315  if (tria_is_parallel)
316  cell_id = used_cells_ids[idx / n_q_points];
317  else
318  cell_id = idx / n_q_points;
319 
320  if (last_id != cell_id)
321  {
322  cell_sets[cell_id].insert(all_cells[i]);
323  last_id = cell_id;
324  }
325  }
326  }
327 
328  // Now we run on each cell of the immersed
329  // and build the sparsity
330  unsigned int i = 0;
331  for (const auto &cell : immersed_dh.active_cell_iterators())
332  {
333  // Reinitialize the cell
334  cell->get_dof_indices(dofs);
335 
336  // List of outer cells
337  const auto &cells = cell_sets[i];
338 
339  for (const auto &cell_c : cells)
340  {
341  // Get the ones in the current outer cell
342  typename DoFHandler<dim0, spacedim>::cell_iterator ocell(*cell_c,
343  &space_dh);
344  // Make sure we act only on locally_owned cells
345  if (ocell->is_locally_owned())
346  {
347  ocell->get_dof_indices(odofs);
348  // [TODO]: When the following function will be implemented
349  // for the case of non-trivial dof_mask, we should
350  // uncomment the missing part.
351  constraints.add_entries_local_to_global(
352  odofs, dofs, sparsity); //, true, dof_mask);
353  }
354  }
355  ++i;
356  }
357  }
358 
359 
360 
361  template <int dim0, int dim1, int spacedim, typename Matrix>
362  void
363  create_coupling_mass_matrix(
364  const DoFHandler<dim0, spacedim> & space_dh,
365  const DoFHandler<dim1, spacedim> & immersed_dh,
366  const Quadrature<dim1> & quad,
367  Matrix & matrix,
368  const AffineConstraints<typename Matrix::value_type> &constraints,
369  const ComponentMask & space_comps,
370  const ComponentMask & immersed_comps,
371  const Mapping<dim0, spacedim> & space_mapping,
372  const Mapping<dim1, spacedim> & immersed_mapping)
373  {
374  GridTools::Cache<dim0, spacedim> cache(space_dh.get_triangulation(),
375  space_mapping);
376  create_coupling_mass_matrix(cache,
377  space_dh,
378  immersed_dh,
379  quad,
380  matrix,
381  constraints,
382  space_comps,
383  immersed_comps,
384  immersed_mapping);
385  }
386 
387 
388 
389  template <int dim0, int dim1, int spacedim, typename Matrix>
390  void
391  create_coupling_mass_matrix(
392  const GridTools::Cache<dim0, spacedim> & cache,
393  const DoFHandler<dim0, spacedim> & space_dh,
394  const DoFHandler<dim1, spacedim> & immersed_dh,
395  const Quadrature<dim1> & quad,
396  Matrix & matrix,
397  const AffineConstraints<typename Matrix::value_type> &constraints,
398  const ComponentMask & space_comps,
399  const ComponentMask & immersed_comps,
400  const Mapping<dim1, spacedim> & immersed_mapping)
401  {
402  AssertDimension(matrix.m(), space_dh.n_dofs());
403  AssertDimension(matrix.n(), immersed_dh.n_dofs());
404  Assert(dim1 <= dim0,
405  ExcMessage("This function can only work if dim1 <= dim0"));
406  Assert((dynamic_cast<
407  const parallel::distributed::Triangulation<dim1, spacedim> *>(
408  &immersed_dh.get_triangulation()) == nullptr),
409  ExcNotImplemented());
410 
411  const bool tria_is_parallel =
412  (dynamic_cast<const parallel::TriangulationBase<dim1, spacedim> *>(
413  &space_dh.get_triangulation()) != nullptr);
414 
415  const auto &space_fe = space_dh.get_fe();
416  const auto &immersed_fe = immersed_dh.get_fe();
417 
418  // Dof indices
419  std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
420  std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
421 
422  // Take care of components
423  const ComponentMask space_c =
424  (space_comps.size() == 0 ? ComponentMask(space_fe.n_components(), true) :
425  space_comps);
426 
427  const ComponentMask immersed_c =
428  (immersed_comps.size() == 0 ?
429  ComponentMask(immersed_fe.n_components(), true) :
430  immersed_comps);
431 
432  AssertDimension(space_c.size(), space_fe.n_components());
433  AssertDimension(immersed_c.size(), immersed_fe.n_components());
434 
435  std::vector<unsigned int> space_gtl(space_fe.n_components(),
436  numbers::invalid_unsigned_int);
437  std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(),
438  numbers::invalid_unsigned_int);
439 
440  for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
441  if (space_c[i])
442  space_gtl[i] = j++;
443 
444  for (unsigned int i = 0, j = 0; i < immersed_gtl.size(); ++i)
445  if (immersed_c[i])
446  immersed_gtl[i] = j++;
447 
448  FullMatrix<typename Matrix::value_type> cell_matrix(
449  space_dh.get_fe().dofs_per_cell, immersed_dh.get_fe().dofs_per_cell);
450 
451  FEValues<dim1, spacedim> fe_v(immersed_mapping,
452  immersed_dh.get_fe(),
453  quad,
454  update_JxW_values | update_quadrature_points |
455  update_values);
456 
457  const unsigned int n_q_points = quad.size();
458  const unsigned int n_active_c =
459  immersed_dh.get_triangulation().n_active_cells();
460 
461  const auto used_cells_data = internal::qpoints_over_locally_owned_cells(
462  cache, immersed_dh, quad, immersed_mapping, tria_is_parallel);
463 
464  const auto &points_over_local_cells = std::get<0>(used_cells_data);
465  const auto &used_cells_ids = std::get<1>(used_cells_data);
466 
467  // Get a list of outer cells, qpoints and maps.
468  const auto cpm =
469  GridTools::compute_point_locations(cache, points_over_local_cells);
470  const auto &all_cells = std::get<0>(cpm);
471  const auto &all_qpoints = std::get<1>(cpm);
472  const auto &all_maps = std::get<2>(cpm);
473 
474  std::vector<
475  std::vector<typename Triangulation<dim0, spacedim>::active_cell_iterator>>
476  cell_container(n_active_c);
477  std::vector<std::vector<std::vector<Point<dim0>>>> qpoints_container(
478  n_active_c);
479  std::vector<std::vector<std::vector<unsigned int>>> maps_container(
480  n_active_c);
481 
482  // Cycle over all cells of underling mesh found
483  // call it omesh, elaborating the output
484  for (unsigned int o = 0; o < all_cells.size(); ++o)
485  {
486  for (unsigned int j = 0; j < all_maps[o].size(); ++j)
487  {
488  // Find the index of the "owner" cell and qpoint
489  // with regard to the immersed mesh
490  // Find in which cell of immersed triangulation the point lies
491  unsigned int cell_id;
492  if (tria_is_parallel)
493  cell_id = used_cells_ids[all_maps[o][j] / n_q_points];
494  else
495  cell_id = all_maps[o][j] / n_q_points;
496 
497  const unsigned int n_pt = all_maps[o][j] % n_q_points;
498 
499  // If there are no cells, we just add our data
500  if (cell_container[cell_id].empty())
501  {
502  cell_container[cell_id].emplace_back(all_cells[o]);
503  qpoints_container[cell_id].emplace_back(
504  std::vector<Point<dim0>>{all_qpoints[o][j]});
505  maps_container[cell_id].emplace_back(
506  std::vector<unsigned int>{n_pt});
507  }
508  // If there are already cells, we begin by looking
509  // at the last inserted cell, which is more likely:
510  else if (cell_container[cell_id].back() == all_cells[o])
511  {
512  qpoints_container[cell_id].back().emplace_back(
513  all_qpoints[o][j]);
514  maps_container[cell_id].back().emplace_back(n_pt);
515  }
516  else
517  {
518  // We don't need to check the last element
519  const auto cell_p = std::find(cell_container[cell_id].begin(),
520  cell_container[cell_id].end() - 1,
521  all_cells[o]);
522 
523  if (cell_p == cell_container[cell_id].end() - 1)
524  {
525  cell_container[cell_id].emplace_back(all_cells[o]);
526  qpoints_container[cell_id].emplace_back(
527  std::vector<Point<dim0>>{all_qpoints[o][j]});
528  maps_container[cell_id].emplace_back(
529  std::vector<unsigned int>{n_pt});
530  }
531  else
532  {
533  const unsigned int pos =
534  cell_p - cell_container[cell_id].begin();
535  qpoints_container[cell_id][pos].emplace_back(
536  all_qpoints[o][j]);
537  maps_container[cell_id][pos].emplace_back(n_pt);
538  }
539  }
540  }
541  }
542 
543  typename DoFHandler<dim1, spacedim>::active_cell_iterator
544  cell = immersed_dh.begin_active(),
545  endc = immersed_dh.end();
546 
547  for (unsigned int j = 0; cell != endc; ++cell, ++j)
548  {
549  // Reinitialize the cell and the fe_values
550  fe_v.reinit(cell);
551  cell->get_dof_indices(dofs);
552 
553  // Get a list of outer cells, qpoints and maps.
554  const auto &cells = cell_container[j];
555  const auto &qpoints = qpoints_container[j];
556  const auto &maps = maps_container[j];
557 
558  for (unsigned int c = 0; c < cells.size(); ++c)
559  {
560  // Get the ones in the current outer cell
561  typename DoFHandler<dim0, spacedim>::active_cell_iterator ocell(
562  *cells[c], &space_dh);
563  // Make sure we act only on locally_owned cells
564  if (ocell->is_locally_owned())
565  {
566  const std::vector<Point<dim0>> & qps = qpoints[c];
567  const std::vector<unsigned int> &ids = maps[c];
568 
569  FEValues<dim0, spacedim> o_fe_v(cache.get_mapping(),
570  space_dh.get_fe(),
571  qps,
572  update_values);
573  o_fe_v.reinit(ocell);
574  ocell->get_dof_indices(odofs);
575 
576  // Reset the matrices.
577  cell_matrix = typename Matrix::value_type();
578 
579  for (unsigned int i = 0; i < space_dh.get_fe().dofs_per_cell;
580  ++i)
581  {
582  const auto comp_i =
583  space_dh.get_fe().system_to_component_index(i).first;
584  if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
585  for (unsigned int j = 0;
586  j < immersed_dh.get_fe().dofs_per_cell;
587  ++j)
588  {
589  const auto comp_j = immersed_dh.get_fe()
590  .system_to_component_index(j)
591  .first;
592  if (space_gtl[comp_i] == immersed_gtl[comp_j])
593  for (unsigned int oq = 0;
594  oq < o_fe_v.n_quadrature_points;
595  ++oq)
596  {
597  // Get the corresponding q point
598  const unsigned int q = ids[oq];
599 
600  cell_matrix(i, j) +=
601  (fe_v.shape_value(j, q) *
602  o_fe_v.shape_value(i, oq) * fe_v.JxW(q));
603  }
604  }
605  }
606 
607  // Now assemble the matrices
608  constraints.distribute_local_to_global(cell_matrix,
609  odofs,
610  dofs,
611  matrix);
612  }
613  }
614  }
615  }
616 
617  template <int dim0,
618  int dim1,
619  int spacedim,
620  typename Sparsity,
621  typename Number>
622  void
623  create_coupling_sparsity_pattern(
624  const double & epsilon,
625  const GridTools::Cache<dim0, spacedim> &cache0,
626  const GridTools::Cache<dim1, spacedim> &cache1,
627  const DoFHandler<dim0, spacedim> & dh0,
628  const DoFHandler<dim1, spacedim> & dh1,
629  const Quadrature<dim1> & quad,
630  Sparsity & sparsity,
631  const AffineConstraints<Number> & constraints0,
632  const ComponentMask & comps0,
633  const ComponentMask & comps1)
634  {
635  if (epsilon == 0.0)
636  {
637  Assert(dim1 <= dim0,
638  ExcMessage("When epsilon is zero, you can only "
639  "call this function with dim1 <= dim0."));
640  create_coupling_sparsity_pattern(cache0,
641  dh0,
642  dh1,
643  quad,
644  sparsity,
645  constraints0,
646  comps0,
647  comps1,
648  cache1.get_mapping());
649  return;
650  }
651  AssertDimension(sparsity.n_rows(), dh0.n_dofs());
652  AssertDimension(sparsity.n_cols(), dh1.n_dofs());
653 
654  const bool zero_is_distributed =
655  (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
656  *>(&dh0.get_triangulation()) != nullptr);
657  const bool one_is_distributed =
658  (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
659  *>(&dh1.get_triangulation()) != nullptr);
660 
661  // We bail out if both are distributed triangulations
662  Assert(!zero_is_distributed || !one_is_distributed, ExcNotImplemented());
663 
664  // If we can loop on both, we decide where to make the outer loop according
665  // to the size of the triangulation. The reasoning is the following:
666  // - cost for accessing the tree: log(N)
667  // - cost for computing the intersection for each of the outer loop cells: M
668  // Total cost (besides the setup) is: M log(N)
669  // If we can, make sure M is the smaller number of the two.
670  const bool outer_loop_on_zero =
671  (zero_is_distributed && !one_is_distributed) ||
672  (dh1.get_triangulation().n_active_cells() >
673  dh0.get_triangulation().n_active_cells());
674 
675  const auto &fe0 = dh0.get_fe();
676  const auto &fe1 = dh1.get_fe();
677 
678  // Dof indices
679  std::vector<types::global_dof_index> dofs0(fe0.dofs_per_cell);
680  std::vector<types::global_dof_index> dofs1(fe1.dofs_per_cell);
681 
682  if (outer_loop_on_zero)
683  {
684  Assert(one_is_distributed == false, ExcInternalError());
685 
686  const auto &tree1 = cache1.get_cell_bounding_boxes_rtree();
687 
688  std::vector<std::pair<
689  BoundingBox<spacedim>,
690  typename Triangulation<dim1, spacedim>::active_cell_iterator>>
691  intersection;
692 
693  for (const auto &cell0 : dh0.active_cell_iterators())
694  if (cell0->is_locally_owned())
695  {
696  intersection.resize(0);
697  BoundingBox<spacedim> box0 =
698  cache0.get_mapping().get_bounding_box(cell0);
699  box0.extend(epsilon);
700  boost::geometry::index::query(tree1,
701  boost::geometry::index::intersects(
702  box0),
703  std::back_inserter(intersection));
704  if (!intersection.empty())
705  {
706  cell0->get_dof_indices(dofs0);
707  for (const auto &entry : intersection)
708  {
709  typename DoFHandler<dim1, spacedim>::cell_iterator cell1(
710  *entry.second, &dh1);
711  cell1->get_dof_indices(dofs1);
712  constraints0.add_entries_local_to_global(dofs0,
713  dofs1,
714  sparsity);
715  }
716  }
717  }
718  }
719  else
720  {
721  Assert(zero_is_distributed == false, ExcInternalError());
722  const auto &tree0 = cache0.get_cell_bounding_boxes_rtree();
723 
724  std::vector<std::pair<
725  BoundingBox<spacedim>,
726  typename Triangulation<dim0, spacedim>::active_cell_iterator>>
727  intersection;
728 
729  for (const auto &cell1 : dh1.active_cell_iterators())
730  if (cell1->is_locally_owned())
731  {
732  intersection.resize(0);
733  BoundingBox<spacedim> box1 =
734  cache1.get_mapping().get_bounding_box(cell1);
735  box1.extend(epsilon);
736  boost::geometry::index::query(tree0,
737  boost::geometry::index::intersects(
738  box1),
739  std::back_inserter(intersection));
740  if (!intersection.empty())
741  {
742  cell1->get_dof_indices(dofs1);
743  for (const auto &entry : intersection)
744  {
745  typename DoFHandler<dim0, spacedim>::cell_iterator cell0(
746  *entry.second, &dh0);
747  cell0->get_dof_indices(dofs0);
748  constraints0.add_entries_local_to_global(dofs0,
749  dofs1,
750  sparsity);
751  }
752  }
753  }
754  }
755  }
756 
757 
758 
759  template <int dim0, int dim1, int spacedim, typename Matrix>
760  void
761  create_coupling_mass_matrix(
762  Functions::CutOffFunctionBase<spacedim> & kernel,
763  const double & epsilon,
764  const GridTools::Cache<dim0, spacedim> & cache0,
765  const GridTools::Cache<dim1, spacedim> & cache1,
766  const DoFHandler<dim0, spacedim> & dh0,
767  const DoFHandler<dim1, spacedim> & dh1,
768  const Quadrature<dim0> & quadrature0,
769  const Quadrature<dim1> & quadrature1,
770  Matrix & matrix,
771  const AffineConstraints<typename Matrix::value_type> &constraints0,
772  const ComponentMask & comps0,
773  const ComponentMask & comps1)
774  {
775  if (epsilon == 0)
776  {
777  Assert(dim1 <= dim0,
778  ExcMessage("When epsilon is zero, you can only "
779  "call this function with dim1 <= dim0."));
780  create_coupling_mass_matrix(cache0,
781  dh0,
782  dh1,
783  quadrature1,
784  matrix,
785  constraints0,
786  comps0,
787  comps1,
788  cache1.get_mapping());
789  return;
790  }
791 
792  AssertDimension(matrix.m(), dh0.n_dofs());
793  AssertDimension(matrix.n(), dh1.n_dofs());
794 
795  const bool zero_is_distributed =
796  (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
797  *>(&dh0.get_triangulation()) != nullptr);
798  const bool one_is_distributed =
799  (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
800  *>(&dh1.get_triangulation()) != nullptr);
801 
802  // We bail out if both are distributed triangulations
803  Assert(!zero_is_distributed || !one_is_distributed, ExcNotImplemented());
804 
805  // If we can loop on both, we decide where to make the outer loop according
806  // to the size of the triangulation. The reasoning is the following:
807  // - cost for accessing the tree: log(N)
808  // - cost for computing the intersection for each of the outer loop cells: M
809  // Total cost (besides the setup) is: M log(N)
810  // If we can, make sure M is the smaller number of the two.
811  const bool outer_loop_on_zero =
812  (zero_is_distributed && !one_is_distributed) ||
813  (dh1.get_triangulation().n_active_cells() >
814  dh0.get_triangulation().n_active_cells());
815 
816  const auto &fe0 = dh0.get_fe();
817  const auto &fe1 = dh1.get_fe();
818 
819  FEValues<dim0, spacedim> fev0(cache0.get_mapping(),
820  fe0,
821  quadrature0,
822  update_values | update_JxW_values |
823  update_quadrature_points);
824 
825  FEValues<dim1, spacedim> fev1(cache1.get_mapping(),
826  fe1,
827  quadrature1,
828  update_values | update_JxW_values |
829  update_quadrature_points);
830 
831  // Dof indices
832  std::vector<types::global_dof_index> dofs0(fe0.dofs_per_cell);
833  std::vector<types::global_dof_index> dofs1(fe1.dofs_per_cell);
834 
835  // Local Matrix
836  FullMatrix<typename Matrix::value_type> cell_matrix(fe0.dofs_per_cell,
837  fe1.dofs_per_cell);
838 
839  // Global to local indices
840  const auto p =
841  internal::compute_components_coupling(comps0, comps1, fe0, fe1);
842  const auto &gtl0 = p.first;
843  const auto &gtl1 = p.second;
844 
845  kernel.set_radius(epsilon);
846  std::vector<double> kernel_values(quadrature1.size());
847 
848  auto assemble_one_pair = [&]() {
849  cell_matrix = 0;
850  for (unsigned int q0 = 0; q0 < quadrature0.size(); ++q0)
851  {
852  kernel.set_center(fev0.quadrature_point(q0));
853  kernel.value_list(fev1.get_quadrature_points(), kernel_values);
854  for (unsigned int j = 0; j < fe1.dofs_per_cell; ++j)
855  {
856  const auto comp_j = fe1.system_to_component_index(j).first;
857 
858  // First compute the part of the integral that does not
859  // depend on i
860  typename Matrix::value_type sum_q1 = {};
861  for (unsigned int q1 = 0; q1 < quadrature1.size(); ++q1)
862  sum_q1 +=
863  fev1.shape_value(j, q1) * kernel_values[q1] * fev1.JxW(q1);
864  sum_q1 *= fev0.JxW(q0);
865 
866  // Now compute the main integral with the sum over q1 already
867  // completed - this gives a cubic complexity as usual rather
868  // than a quartic one with naive loops
869  for (unsigned int i = 0; i < fe0.dofs_per_cell; ++i)
870  {
871  const auto comp_i = fe0.system_to_component_index(i).first;
872  if (gtl0[comp_i] != numbers::invalid_unsigned_int &&
873  gtl1[comp_j] == gtl0[comp_i])
874  cell_matrix(i, j) += fev0.shape_value(i, q0) * sum_q1;
875  }
876  }
877  }
878 
879  constraints0.distribute_local_to_global(cell_matrix,
880  dofs0,
881  dofs1,
882  matrix);
883  };
884 
885  if (outer_loop_on_zero)
886  {
887  Assert(one_is_distributed == false, ExcInternalError());
888 
889  const auto &tree1 = cache1.get_cell_bounding_boxes_rtree();
890 
891  std::vector<std::pair<
892  BoundingBox<spacedim>,
893  typename Triangulation<dim1, spacedim>::active_cell_iterator>>
894  intersection;
895 
896  for (const auto &cell0 : dh0.active_cell_iterators())
897  if (cell0->is_locally_owned())
898  {
899  intersection.resize(0);
900  BoundingBox<spacedim> box0 =
901  cache0.get_mapping().get_bounding_box(cell0);
902  box0.extend(epsilon);
903  boost::geometry::index::query(tree1,
904  boost::geometry::index::intersects(
905  box0),
906  std::back_inserter(intersection));
907  if (!intersection.empty())
908  {
909  cell0->get_dof_indices(dofs0);
910  fev0.reinit(cell0);
911  for (const auto &entry : intersection)
912  {
913  typename DoFHandler<dim1, spacedim>::cell_iterator cell1(
914  *entry.second, &dh1);
915  cell1->get_dof_indices(dofs1);
916  fev1.reinit(cell1);
917  assemble_one_pair();
918  }
919  }
920  }
921  }
922  else
923  {
924  Assert(zero_is_distributed == false, ExcInternalError());
925  const auto &tree0 = cache0.get_cell_bounding_boxes_rtree();
926 
927  std::vector<std::pair<
928  BoundingBox<spacedim>,
929  typename Triangulation<dim0, spacedim>::active_cell_iterator>>
930  intersection;
931 
932  for (const auto &cell1 : dh1.active_cell_iterators())
933  if (cell1->is_locally_owned())
934  {
935  intersection.resize(0);
936  BoundingBox<spacedim> box1 =
937  cache1.get_mapping().get_bounding_box(cell1);
938  box1.extend(epsilon);
939  boost::geometry::index::query(tree0,
940  boost::geometry::index::intersects(
941  box1),
942  std::back_inserter(intersection));
943  if (!intersection.empty())
944  {
945  cell1->get_dof_indices(dofs1);
946  fev1.reinit(cell1);
947  for (const auto &entry : intersection)
948  {
949  typename DoFHandler<dim0, spacedim>::cell_iterator cell0(
950  *entry.second, &dh0);
951  cell0->get_dof_indices(dofs0);
952  fev0.reinit(cell0);
953  assemble_one_pair();
954  }
955  }
956  }
957  }
958  }
959 
960 #include "coupling.inst"
961 } // namespace NonMatching
962 
963 DEAL_II_NAMESPACE_CLOSE
NonMatching::internal::qpoints_over_locally_owned_cells
std::pair< std::vector< Point< spacedim > >, std::vector< unsigned int > > qpoints_over_locally_owned_cells(const GridTools::Cache< dim0, spacedim > &cache, const DoFHandler< dim1, spacedim > &immersed_dh, const Quadrature< dim1 > &quad, const Mapping< dim1, spacedim > &immersed_mapping, const bool tria_is_parallel)
Definition: coupling.cc:64
Physics::Elasticity::Kinematics::epsilon
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
fe_values.h
sparse_matrix.h
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:357
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:329
trilinos_sparsity_pattern.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
petsc_sparse_matrix.h
Function< dim >::value_list
virtual void value_list(const std::vector< Point< dim >> &points, std::vector< double > &values, const unsigned int component=0) const
Functions::CutOffFunctionBase
Definition: function_lib.h:924
trilinos_sparse_matrix.h
ComponentMask
Definition: component_mask.h:83
block_sparsity_pattern.h
LocalIntegrators::Advection::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &velocity, const double factor=1.)
Definition: advection.h:80
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
DoFHandler
Definition: dof_handler.h:205
BoundingBox::extend
void extend(const Number &amount)
FEValues
Definition: fe.h:38
BoundingBox
Definition: bounding_box.h:128
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: dof_handler.cc:1030
ComponentMask::size
unsigned int size() const
Mapping
Abstract base class for mapping classes.
Definition: mapping.h:302
tensor.h
AffineConstraints::distribute_local_to_global
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
Definition: affine_constraints.h:1859
trilinos_block_sparse_matrix.h
AffineConstraints::add_entries_local_to_global
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
Definition: affine_constraints.h:2034
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
grid_tools.h
GridTools::compute_point_locations
return_type compute_point_locations(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:4193
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
FEValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
parallel::distributed::Triangulation
Definition: tria.h:242
shared_tria.h
GridTools::Cache::get_cell_bounding_boxes_rtree
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_cell_bounding_boxes_rtree() const
Definition: grid_tools_cache.cc:126
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
sparsity_pattern.h
GridTools::Cache
Definition: grid_tools.h:128
NonMatching
Definition: coupling.h:45
FiniteElement
Definition: fe.h:648
exceptions.h
TriaActiveIterator
Definition: tria_iterator.h:759
DoFHandler::end
cell_iterator end() const
Definition: dof_handler.cc:951
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
tria.h
AffineConstraints
Definition: affine_constraints.h:180
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
block_sparse_matrix.h
grid_tools_cache.h
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
NonMatching::internal::compute_components_coupling
std::pair< std::vector< unsigned int >, std::vector< unsigned int > > compute_components_coupling(const ComponentMask &comps0, const ComponentMask &comps1, const FiniteElement< dim0, spacedim > &fe0, const FiniteElement< dim1, spacedim > &fe1)
Definition: coupling.cc:137
Point
Definition: point.h:111
FiniteElementData::n_components
unsigned int n_components() const
parallel::TriangulationBase
Definition: tria_base.h:78
Quadrature::size
unsigned int size() const
FullMatrix
Definition: full_matrix.h:71
coupling.h
internal
Definition: aligned_vector.h:369
NonMatching::create_coupling_mass_matrix
void create_coupling_mass_matrix(const DoFHandler< dim0, spacedim > &space_dh, const DoFHandler< dim1, spacedim > &immersed_dh, const Quadrature< dim1 > &quad, Matrix &matrix, const AffineConstraints< typename Matrix::value_type > &constraints=AffineConstraints< typename Matrix::value_type >(), const ComponentMask &space_comps=ComponentMask(), const ComponentMask &immersed_comps=ComponentMask(), const Mapping< dim0, spacedim > &space_mapping=StaticMappingQ1< dim0, spacedim >::mapping, const Mapping< dim1, spacedim > &immersed_mapping=StaticMappingQ1< dim1, spacedim >::mapping)
Definition: coupling.cc:363
GridTools::Cache::get_mapping
const Mapping< dim, spacedim > & get_mapping() const
Definition: grid_tools_cache.h:274
Functions::CutOffFunctionBase::set_center
virtual void set_center(const Point< dim > &p)
Definition: function_lib_cutoff.cc:62
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Quadrature
Definition: quadrature.h:85
NonMatching::create_coupling_sparsity_pattern
void create_coupling_sparsity_pattern(const DoFHandler< dim0, spacedim > &space_dh, const DoFHandler< dim1, spacedim > &immersed_dh, const Quadrature< dim1 > &quad, Sparsity &sparsity, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const ComponentMask &space_comps=ComponentMask(), const ComponentMask &immersed_comps=ComponentMask(), const Mapping< dim0, spacedim > &space_mapping=StaticMappingQ1< dim0, spacedim >::mapping, const Mapping< dim1, spacedim > &immersed_mapping=StaticMappingQ1< dim1, spacedim >::mapping)
Definition: coupling.cc:175
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:935
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
petsc_block_sparse_matrix.h
Functions::CutOffFunctionBase::set_radius
virtual void set_radius(const double r)
Definition: function_lib_cutoff.cc:89
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
point.h
Generated by   doxygen 1.8.17