Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
coupling.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/exceptions.h>
17 #include <deal.II/base/point.h>
18 #include <deal.II/base/tensor.h>
19 
20 #include <deal.II/distributed/shared_tria.h>
21 #include <deal.II/distributed/tria.h>
22 
23 #include <deal.II/fe/fe_values.h>
24 
25 #include <deal.II/grid/grid_tools.h>
26 #include <deal.II/grid/grid_tools_cache.h>
27 
28 #include <deal.II/lac/block_sparse_matrix.h>
29 #include <deal.II/lac/block_sparsity_pattern.h>
30 #include <deal.II/lac/petsc_block_sparse_matrix.h>
31 #include <deal.II/lac/petsc_sparse_matrix.h>
32 #include <deal.II/lac/sparse_matrix.h>
33 #include <deal.II/lac/sparsity_pattern.h>
34 #include <deal.II/lac/trilinos_block_sparse_matrix.h>
35 #include <deal.II/lac/trilinos_sparse_matrix.h>
36 #include <deal.II/lac/trilinos_sparsity_pattern.h>
37 
38 #include <deal.II/non_matching/coupling.h>
39 
40 DEAL_II_NAMESPACE_OPEN
41 namespace NonMatching
42 {
43  namespace internal
44  {
62  template <int dim0, int dim1, int spacedim>
63  std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>
64  qpoints_over_locally_owned_cells(
66  const DoFHandler<dim1, spacedim> & immersed_dh,
67  const Quadrature<dim1> & quad,
68  const Mapping<dim1, spacedim> & immersed_mapping,
69  const bool tria_is_parallel)
70  {
71  const auto & immersed_fe = immersed_dh.get_fe();
72  std::vector<Point<spacedim>> points_over_local_cells;
73  // Keep track of which cells we actually used
74  std::vector<unsigned int> used_cells_ids;
75  {
76  FEValues<dim1, spacedim> fe_v(immersed_mapping,
77  immersed_fe,
78  quad,
80  unsigned int cell_id = 0;
81  for (const auto &cell : immersed_dh.active_cell_iterators())
82  {
83  bool use_cell = false;
84  if (tria_is_parallel)
85  {
86  const auto bbox = cell->bounding_box();
87  std::vector<std::pair<
90  out_vals;
91  cache.get_cell_bounding_boxes_rtree().query(
92  boost::geometry::index::intersects(bbox),
93  std::back_inserter(out_vals));
94  // Each bounding box corresponds to an active cell
95  // of the embedding triangulation: we now check if
96  // the current cell, of the embedded triangulation,
97  // overlaps a locally owned cell of the embedding one
98  for (const auto &bbox_it : out_vals)
99  if (bbox_it.second->is_locally_owned())
100  {
101  use_cell = true;
102  used_cells_ids.emplace_back(cell_id);
103  break;
104  }
105  }
106  else
107  // for sequential triangulations, simply use all cells
108  use_cell = true;
109 
110  if (use_cell)
111  {
112  // Reinitialize the cell and the fe_values
113  fe_v.reinit(cell);
114  const std::vector<Point<spacedim>> &x_points =
115  fe_v.get_quadrature_points();
116 
117  // Insert the points to the vector
118  points_over_local_cells.insert(points_over_local_cells.end(),
119  x_points.begin(),
120  x_points.end());
121  }
122  ++cell_id;
123  }
124  }
125  return {std::move(points_over_local_cells), std::move(used_cells_ids)};
126  }
127  } // namespace internal
128 
129  template <int dim0,
130  int dim1,
131  int spacedim,
132  typename Sparsity,
133  typename number>
134  void
136  const DoFHandler<dim0, spacedim> &space_dh,
137  const DoFHandler<dim1, spacedim> &immersed_dh,
138  const Quadrature<dim1> & quad,
139  Sparsity & sparsity,
140  const AffineConstraints<number> & constraints,
141  const ComponentMask & space_comps,
142  const ComponentMask & immersed_comps,
143  const Mapping<dim0, spacedim> & space_mapping,
144  const Mapping<dim1, spacedim> & immersed_mapping)
145  {
147  space_mapping);
149  space_dh,
150  immersed_dh,
151  quad,
152  sparsity,
153  constraints,
154  space_comps,
155  immersed_comps,
156  immersed_mapping);
157  }
158 
159 
160 
161  template <int dim0,
162  int dim1,
163  int spacedim,
164  typename Sparsity,
165  typename number>
166  void
169  const DoFHandler<dim0, spacedim> & space_dh,
170  const DoFHandler<dim1, spacedim> & immersed_dh,
171  const Quadrature<dim1> & quad,
172  Sparsity & sparsity,
173  const AffineConstraints<number> & constraints,
174  const ComponentMask & space_comps,
175  const ComponentMask & immersed_comps,
176  const Mapping<dim1, spacedim> & immersed_mapping)
177  {
178  AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
179  AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
180  static_assert(dim1 <= dim0, "This function can only work if dim1 <= dim0");
181  Assert((dynamic_cast<
183  &immersed_dh.get_triangulation()) == nullptr),
185 
186  const bool tria_is_parallel =
187  (dynamic_cast<const parallel::Triangulation<dim1, spacedim> *>(
188  &space_dh.get_triangulation()) != nullptr);
189  const auto &space_fe = space_dh.get_fe();
190  const auto &immersed_fe = immersed_dh.get_fe();
191 
192  // Dof indices
193  std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
194  std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
195 
196  // Take care of components
197  const ComponentMask space_c =
198  (space_comps.size() == 0 ? ComponentMask(space_fe.n_components(), true) :
199  space_comps);
200 
201  const ComponentMask immersed_c =
202  (immersed_comps.size() == 0 ?
203  ComponentMask(immersed_fe.n_components(), true) :
204  immersed_comps);
205 
206  AssertDimension(space_c.size(), space_fe.n_components());
207  AssertDimension(immersed_c.size(), immersed_fe.n_components());
208 
209  // Global to local indices
210  std::vector<unsigned int> space_gtl(space_fe.n_components(),
212  std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(),
214 
215  for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
216  if (space_c[i])
217  space_gtl[i] = j++;
218 
219  for (unsigned int i = 0, j = 0; i < immersed_gtl.size(); ++i)
220  if (immersed_c[i])
221  immersed_gtl[i] = j++;
222 
223  const unsigned int n_q_points = quad.size();
224  const unsigned int n_active_c =
225  immersed_dh.get_triangulation().n_active_cells();
226 
227  const auto qpoints_cells_data = internal::qpoints_over_locally_owned_cells(
228  cache, immersed_dh, quad, immersed_mapping, tria_is_parallel);
229 
230  const auto &points_over_local_cells = std::get<0>(qpoints_cells_data);
231  const auto &used_cells_ids = std::get<1>(qpoints_cells_data);
232 
233  // [TODO]: when the add_entries_local_to_global below will implement
234  // the version with the dof_mask, this should be uncommented.
235  //
236  // // Construct a dof_mask, used to distribute entries to the sparsity
237  // able< 2, bool > dof_mask(space_fe.dofs_per_cell,
238  // immersed_fe.dofs_per_cell);
239  // of_mask.fill(false);
240  // or (unsigned int i=0; i<space_fe.dofs_per_cell; ++i)
241  // {
242  // const auto comp_i = space_fe.system_to_component_index(i).first;
243  // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
244  // for (unsigned int j=0; j<immersed_fe.dofs_per_cell; ++j)
245  // {
246  // const auto comp_j =
247  // immersed_fe.system_to_component_index(j).first; if
248  // (immersed_gtl[comp_j] == space_gtl[comp_i])
249  // dof_mask(i,j) = true;
250  // }
251  // }
252 
253 
254  // Get a list of outer cells, qpoints and maps.
255  const auto cpm =
256  GridTools::compute_point_locations(cache, points_over_local_cells);
257  const auto &all_cells = std::get<0>(cpm);
258  const auto &maps = std::get<2>(cpm);
259 
260  std::vector<
261  std::set<typename Triangulation<dim0, spacedim>::active_cell_iterator>>
262  cell_sets(n_active_c);
263 
264  for (unsigned int i = 0; i < maps.size(); ++i)
265  {
266  // Quadrature points should be reasonably clustered:
267  // the following index keeps track of the last id
268  // where the current cell was inserted
269  unsigned int last_id = std::numeric_limits<unsigned int>::max();
270  unsigned int cell_id;
271  for (const unsigned int idx : maps[i])
272  {
273  // Find in which cell of immersed triangulation the point lies
274  if (tria_is_parallel)
275  cell_id = used_cells_ids[idx / n_q_points];
276  else
277  cell_id = idx / n_q_points;
278 
279  if (last_id != cell_id)
280  {
281  cell_sets[cell_id].insert(all_cells[i]);
282  last_id = cell_id;
283  }
284  }
285  }
286 
287  // Now we run on each cell of the immersed
288  // and build the sparsity
289  unsigned int i = 0;
290  for (const auto &cell : immersed_dh.active_cell_iterators())
291  {
292  // Reinitialize the cell
293  cell->get_dof_indices(dofs);
294 
295  // List of outer cells
296  const auto &cells = cell_sets[i];
297 
298  for (const auto &cell_c : cells)
299  {
300  // Get the ones in the current outer cell
301  typename DoFHandler<dim0, spacedim>::cell_iterator ocell(*cell_c,
302  &space_dh);
303  // Make sure we act only on locally_owned cells
304  if (ocell->is_locally_owned())
305  {
306  ocell->get_dof_indices(odofs);
307  // [TODO]: When the following function will be implemented
308  // for the case of non-trivial dof_mask, we should
309  // uncomment the missing part.
310  constraints.add_entries_local_to_global(
311  odofs, dofs, sparsity); //, true, dof_mask);
312  }
313  }
314  ++i;
315  }
316  }
317 
318 
319 
320  template <int dim0, int dim1, int spacedim, typename Matrix>
321  void
323  const DoFHandler<dim0, spacedim> & space_dh,
324  const DoFHandler<dim1, spacedim> & immersed_dh,
325  const Quadrature<dim1> & quad,
326  Matrix & matrix,
328  const ComponentMask & space_comps,
329  const ComponentMask & immersed_comps,
330  const Mapping<dim0, spacedim> & space_mapping,
331  const Mapping<dim1, spacedim> & immersed_mapping)
332  {
334  space_mapping);
336  space_dh,
337  immersed_dh,
338  quad,
339  matrix,
340  constraints,
341  space_comps,
342  immersed_comps,
343  immersed_mapping);
344  }
345 
346 
347 
348  template <int dim0, int dim1, int spacedim, typename Matrix>
349  void
351  const GridTools::Cache<dim0, spacedim> & cache,
352  const DoFHandler<dim0, spacedim> & space_dh,
353  const DoFHandler<dim1, spacedim> & immersed_dh,
354  const Quadrature<dim1> & quad,
355  Matrix & matrix,
357  const ComponentMask & space_comps,
358  const ComponentMask & immersed_comps,
359  const Mapping<dim1, spacedim> & immersed_mapping)
360  {
361  AssertDimension(matrix.m(), space_dh.n_dofs());
362  AssertDimension(matrix.n(), immersed_dh.n_dofs());
363  static_assert(dim1 <= dim0, "This function can only work if dim1 <= dim0");
364  Assert((dynamic_cast<
366  &immersed_dh.get_triangulation()) == nullptr),
368 
369  const bool tria_is_parallel =
370  (dynamic_cast<const parallel::Triangulation<dim1, spacedim> *>(
371  &space_dh.get_triangulation()) != nullptr);
372 
373  const auto &space_fe = space_dh.get_fe();
374  const auto &immersed_fe = immersed_dh.get_fe();
375 
376  // Dof indices
377  std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
378  std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
379 
380  // Take care of components
381  const ComponentMask space_c =
382  (space_comps.size() == 0 ? ComponentMask(space_fe.n_components(), true) :
383  space_comps);
384 
385  const ComponentMask immersed_c =
386  (immersed_comps.size() == 0 ?
387  ComponentMask(immersed_fe.n_components(), true) :
388  immersed_comps);
389 
390  AssertDimension(space_c.size(), space_fe.n_components());
391  AssertDimension(immersed_c.size(), immersed_fe.n_components());
392 
393  std::vector<unsigned int> space_gtl(space_fe.n_components(),
395  std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(),
397 
398  for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
399  if (space_c[i])
400  space_gtl[i] = j++;
401 
402  for (unsigned int i = 0, j = 0; i < immersed_gtl.size(); ++i)
403  if (immersed_c[i])
404  immersed_gtl[i] = j++;
405 
407  space_dh.get_fe().dofs_per_cell, immersed_dh.get_fe().dofs_per_cell);
408 
409  FEValues<dim1, spacedim> fe_v(immersed_mapping,
410  immersed_dh.get_fe(),
411  quad,
413  update_values);
414 
415  const unsigned int n_q_points = quad.size();
416  const unsigned int n_active_c =
417  immersed_dh.get_triangulation().n_active_cells();
418 
419  const auto used_cells_data = internal::qpoints_over_locally_owned_cells(
420  cache, immersed_dh, quad, immersed_mapping, tria_is_parallel);
421 
422  const auto &points_over_local_cells = std::get<0>(used_cells_data);
423  const auto &used_cells_ids = std::get<1>(used_cells_data);
424 
425  // Get a list of outer cells, qpoints and maps.
426  const auto cpm =
427  GridTools::compute_point_locations(cache, points_over_local_cells);
428  const auto &all_cells = std::get<0>(cpm);
429  const auto &all_qpoints = std::get<1>(cpm);
430  const auto &all_maps = std::get<2>(cpm);
431 
432  std::vector<
433  std::vector<typename Triangulation<dim0, spacedim>::active_cell_iterator>>
434  cell_container(n_active_c);
435  std::vector<std::vector<std::vector<Point<dim0>>>> qpoints_container(
436  n_active_c);
437  std::vector<std::vector<std::vector<unsigned int>>> maps_container(
438  n_active_c);
439 
440  // Cycle over all cells of underling mesh found
441  // call it omesh, elaborating the output
442  for (unsigned int o = 0; o < all_cells.size(); ++o)
443  {
444  for (unsigned int j = 0; j < all_maps[o].size(); ++j)
445  {
446  // Find the index of the "owner" cell and qpoint
447  // with regard to the immersed mesh
448  // Find in which cell of immersed triangulation the point lies
449  unsigned int cell_id;
450  if (tria_is_parallel)
451  cell_id = used_cells_ids[all_maps[o][j] / n_q_points];
452  else
453  cell_id = all_maps[o][j] / n_q_points;
454 
455  const unsigned int n_pt = all_maps[o][j] % n_q_points;
456 
457  // If there are no cells, we just add our data
458  if (cell_container[cell_id].empty())
459  {
460  cell_container[cell_id].emplace_back(all_cells[o]);
461  qpoints_container[cell_id].emplace_back(
462  std::vector<Point<dim0>>{all_qpoints[o][j]});
463  maps_container[cell_id].emplace_back(
464  std::vector<unsigned int>{n_pt});
465  }
466  // If there are already cells, we begin by looking
467  // at the last inserted cell, which is more likely:
468  else if (cell_container[cell_id].back() == all_cells[o])
469  {
470  qpoints_container[cell_id].back().emplace_back(
471  all_qpoints[o][j]);
472  maps_container[cell_id].back().emplace_back(n_pt);
473  }
474  else
475  {
476  // We don't need to check the last element
477  const auto cell_p = std::find(cell_container[cell_id].begin(),
478  cell_container[cell_id].end() - 1,
479  all_cells[o]);
480 
481  if (cell_p == cell_container[cell_id].end() - 1)
482  {
483  cell_container[cell_id].emplace_back(all_cells[o]);
484  qpoints_container[cell_id].emplace_back(
485  std::vector<Point<dim0>>{all_qpoints[o][j]});
486  maps_container[cell_id].emplace_back(
487  std::vector<unsigned int>{n_pt});
488  }
489  else
490  {
491  const unsigned int pos =
492  cell_p - cell_container[cell_id].begin();
493  qpoints_container[cell_id][pos].emplace_back(
494  all_qpoints[o][j]);
495  maps_container[cell_id][pos].emplace_back(n_pt);
496  }
497  }
498  }
499  }
500 
502  cell = immersed_dh.begin_active(),
503  endc = immersed_dh.end();
504 
505  for (unsigned int j = 0; cell != endc; ++cell, ++j)
506  {
507  // Reinitialize the cell and the fe_values
508  fe_v.reinit(cell);
509  cell->get_dof_indices(dofs);
510 
511  // Get a list of outer cells, qpoints and maps.
512  const auto &cells = cell_container[j];
513  const auto &qpoints = qpoints_container[j];
514  const auto &maps = maps_container[j];
515 
516  for (unsigned int c = 0; c < cells.size(); ++c)
517  {
518  // Get the ones in the current outer cell
520  *cells[c], &space_dh);
521  // Make sure we act only on locally_owned cells
522  if (ocell->is_locally_owned())
523  {
524  const std::vector<Point<dim0>> & qps = qpoints[c];
525  const std::vector<unsigned int> &ids = maps[c];
526 
527  FEValues<dim0, spacedim> o_fe_v(cache.get_mapping(),
528  space_dh.get_fe(),
529  qps,
530  update_values);
531  o_fe_v.reinit(ocell);
532  ocell->get_dof_indices(odofs);
533 
534  // Reset the matrices.
535  cell_matrix = typename Matrix::value_type();
536 
537  for (unsigned int i = 0; i < space_dh.get_fe().dofs_per_cell;
538  ++i)
539  {
540  const auto comp_i =
541  space_dh.get_fe().system_to_component_index(i).first;
542  if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
543  for (unsigned int j = 0;
544  j < immersed_dh.get_fe().dofs_per_cell;
545  ++j)
546  {
547  const auto comp_j = immersed_dh.get_fe()
548  .system_to_component_index(j)
549  .first;
550  if (space_gtl[comp_i] == immersed_gtl[comp_j])
551  for (unsigned int oq = 0;
552  oq < o_fe_v.n_quadrature_points;
553  ++oq)
554  {
555  // Get the corresponding q point
556  const unsigned int q = ids[oq];
557 
558  cell_matrix(i, j) +=
559  (fe_v.shape_value(j, q) *
560  o_fe_v.shape_value(i, oq) * fe_v.JxW(q));
561  }
562  }
563  }
564 
565  // Now assemble the matrices
566  constraints.distribute_local_to_global(cell_matrix,
567  odofs,
568  dofs,
569  matrix);
570  }
571  }
572  }
573  }
574 
575 #include "coupling.inst"
576 } // namespace NonMatching
577 
578 
579 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
Shape function values.
const Triangulation< dim, spacedim > & get_triangulation() const
static const unsigned int invalid_unsigned_int
Definition: types.h:173
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
const Mapping< dim, spacedim > & get_mapping() const
cell_iterator end() const
Definition: dof_handler.cc:959
Transformed quadrature points.
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void create_coupling_sparsity_pattern(const DoFHandler< dim0, spacedim > &space_dh, const DoFHandler< dim1, spacedim > &immersed_dh, const Quadrature< dim1 > &quad, Sparsity &sparsity, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const ComponentMask &space_comps=ComponentMask(), const ComponentMask &immersed_comps=ComponentMask(), const Mapping< dim0, spacedim > &space_mapping=StaticMappingQ1< dim0, spacedim >::mapping, const Mapping< dim1, spacedim > &immersed_mapping=StaticMappingQ1< dim1, spacedim >::mapping)
Definition: coupling.cc:135
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_cell_bounding_boxes_rtree() const
Definition: point.h:110
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:943
#define Assert(cond, exc)
Definition: exceptions.h:1407
IteratorRange< active_cell_iterator > active_cell_iterators() const
Abstract base class for mapping classes.
Definition: dof_tools.h:57
types::global_dof_index n_dofs() const
unsigned int size() const
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell)
Definition: fe_values.cc:4555
unsigned int size() const
void create_coupling_mass_matrix(const DoFHandler< dim0, spacedim > &space_dh, const DoFHandler< dim1, spacedim > &immersed_dh, const Quadrature< dim1 > &quad, Matrix &matrix, const AffineConstraints< typename Matrix::value_type > &constraints=AffineConstraints< typename Matrix::value_type >(), const ComponentMask &space_comps=ComponentMask(), const ComponentMask &immersed_comps=ComponentMask(), const Mapping< dim0, spacedim > &space_mapping=StaticMappingQ1< dim0, spacedim >::mapping, const Mapping< dim1, spacedim > &immersed_mapping=StaticMappingQ1< dim1, spacedim >::mapping)
Definition: coupling.cc:322
Definition: fe.h:38
static ::ExceptionBase & ExcNotImplemented()
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:4328
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:352
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:324