Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
24
28
37
39
40#include <limits>
41
43namespace 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;
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
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,
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
574 space_dh.get_fe(),
575 qps,
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,
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);
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);
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,
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
835 FullMatrix<typename Matrix::value_type> cell_matrix(fe0.n_dofs_per_cell(),
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
878 constraints0.distribute_local_to_global(cell_matrix,
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);
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);
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 FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) 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
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:317
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:112
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ 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())
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
static const unsigned int invalid_unsigned_int
Definition types.h:213