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