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