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