Fluid structure interaction suite
create_coupling_sparsity_pattern_with_exact_intersections.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2022 by Luca Heltai
4 //
5 // This file is part of the FSI-suite platform, based on the deal.II library.
6 //
7 // The FSI-suite platform is free software; you can use it, redistribute it,
8 // and/or modify it under the terms of the GNU Lesser General Public License as
9 // published by the Free Software Foundation; either version 3.0 of the License,
10 // or (at your option) any later version. The full text of the license can be
11 // found in the file LICENSE at the top level of the FSI-suite platform
12 // distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include "create_coupling_sparsity_pattern_with_exact_intersections.h"
17 
18 #include "compute_intersections.h"
19 
20 #ifdef DEAL_II_WITH_CGAL
21 
22 using namespace dealii;
23 namespace dealii
24 {
25  namespace NonMatching
26  {
27  template <int dim0,
28  int dim1,
29  int spacedim,
30  typename Sparsity,
31  typename number>
32  void
33  create_coupling_sparsity_pattern_with_exact_intersections(
34  const std::vector<std::tuple<
35  typename dealii::Triangulation<dim0, spacedim>::cell_iterator,
36  typename dealii::Triangulation<dim1, spacedim>::cell_iterator,
37  dealii::Quadrature<spacedim>>> &intersections_info,
38  const DoFHandler<dim0, spacedim> &space_dh,
39  const DoFHandler<dim1, spacedim> &immersed_dh,
40  Sparsity &sparsity,
41  const AffineConstraints<number> &constraints,
42  const ComponentMask &space_comps,
43  const ComponentMask &immersed_comps,
44  const AffineConstraints<number> &immersed_constraints)
45  {
46  AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
47  AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
48  Assert(dim1 <= dim0,
49  ExcMessage("This function can only work if dim1 <= dim0"));
50  Assert((dynamic_cast<
51  const parallel::distributed::Triangulation<dim1, spacedim> *>(
52  &immersed_dh.get_triangulation()) == nullptr),
53  ExcNotImplemented());
54 
55 
56 
57  const auto &space_fe = space_dh.get_fe();
58  const auto &immersed_fe = immersed_dh.get_fe();
59  const unsigned int n_dofs_per_space_cell = space_fe.n_dofs_per_cell();
60  const unsigned int n_dofs_per_immersed_cell =
61  immersed_fe.n_dofs_per_cell();
62  const unsigned int n_space_fe_components = space_fe.n_components();
63  const unsigned int n_immersed_fe_components = immersed_fe.n_components();
64  std::vector<types::global_dof_index> space_dofs(n_dofs_per_space_cell);
65  std::vector<types::global_dof_index> immersed_dofs(
66  n_dofs_per_immersed_cell);
67 
68 
69  const ComponentMask space_c =
70  (space_comps.size() == 0 ? ComponentMask(n_space_fe_components, true) :
71  space_comps);
72 
73 
74  const ComponentMask immersed_c =
75  (immersed_comps.size() == 0 ?
76  ComponentMask(n_immersed_fe_components, true) :
77  immersed_comps);
78 
79  AssertDimension(space_c.size(), n_space_fe_components);
80  AssertDimension(immersed_c.size(), n_immersed_fe_components);
81 
82 
83  // Global 2 Local indices
84  std::vector<unsigned int> space_gtl(n_space_fe_components);
85  std::vector<unsigned int> immersed_gtl(n_immersed_fe_components);
86  for (unsigned int i = 0, j = 0; i < n_space_fe_components; ++i)
87  {
88  if (space_c[i])
89  space_gtl[i] = j++;
90  }
91 
92 
93  for (unsigned int i = 0, j = 0; i < n_immersed_fe_components; ++i)
94  {
95  if (immersed_c[i])
96  immersed_gtl[i] = j++;
97  }
98 
99 
100 
101  Table<2, bool> dof_mask(n_dofs_per_space_cell, n_dofs_per_immersed_cell);
102  dof_mask.fill(false); // start off by assuming they don't couple
103 
104  for (unsigned int i = 0; i < n_dofs_per_space_cell; ++i)
105  {
106  const auto comp_i = space_fe.system_to_component_index(i).first;
107  if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
108  {
109  for (unsigned int j = 0; j < n_dofs_per_immersed_cell; ++j)
110  {
111  const auto comp_j =
112  immersed_fe.system_to_component_index(j).first;
113  if (immersed_gtl[comp_j] == space_gtl[comp_i])
114  {
115  dof_mask(i, j) = true;
116  }
117  }
118  }
119  }
120 
121  const bool dof_mask_is_active =
122  dof_mask.n_rows() == n_dofs_per_space_cell &&
123  dof_mask.n_cols() == n_dofs_per_immersed_cell;
124 
125  // Whenever the BB space_cell intersects the BB of an embedded cell, those
126  // DoFs have to be recorded
127 
128  for (const auto &it : intersections_info)
129  {
130  const auto &space_cell = std::get<0>(it);
131  const auto &immersed_cell = std::get<1>(it);
132  typename DoFHandler<dim0, spacedim>::cell_iterator space_cell_dh(
133  *space_cell, &space_dh);
134  typename DoFHandler<dim1, spacedim>::cell_iterator immersed_cell_dh(
135  *immersed_cell, &immersed_dh);
136 
137  space_cell_dh->get_dof_indices(space_dofs);
138  immersed_cell_dh->get_dof_indices(immersed_dofs);
139 
140  if (dof_mask_is_active)
141  {
142  for (unsigned int i = 0; i < n_dofs_per_space_cell; ++i)
143  {
144  const unsigned int comp_i =
145  space_dh.get_fe().system_to_component_index(i).first;
146  if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
147  {
148  for (unsigned int j = 0; j < n_dofs_per_immersed_cell;
149  ++j)
150  {
151  const unsigned int comp_j =
152  immersed_dh.get_fe()
153  .system_to_component_index(j)
154  .first;
155  if (space_gtl[comp_i] == immersed_gtl[comp_j])
156  {
157  // local_cell_matrix(i, j) +=
158  constraints.add_entries_local_to_global(
159  {space_dofs[i]},
160  immersed_constraints,
161  {immersed_dofs[j]},
162  sparsity,
163  true);
164  }
165  }
166  }
167  }
168  }
169  else
170  {
171  constraints.add_entries_local_to_global(space_dofs,
172  immersed_constraints,
173  immersed_dofs,
174  sparsity,
175  true,
176  dof_mask);
177  }
178  }
179  }
180 
181 #else
182 
183 namespace dealii
184 {
185  namespace NonMatching
186  {
187  template <int dim0,
188  int dim1,
189  int spacedim,
190  typename Sparsity,
191  typename number>
192  void
193 
194  create_coupling_sparsity_pattern_with_exact_intersections(
195  const std::vector<std::tuple<
196  typename dealii::Triangulation<dim0, spacedim>::cell_iterator,
197  typename dealii::Triangulation<dim1, spacedim>::cell_iterator,
198  dealii::Quadrature<spacedim>>> &,
199  const DoFHandler<dim0, spacedim> &,
200  const DoFHandler<dim1, spacedim> &,
201  Sparsity &,
202  const AffineConstraints<number> &,
203  const ComponentMask &,
204  const ComponentMask &,
205  const AffineConstraints<number> &)
206  {
207  Assert(false,
208  ExcMessage("This function needs CGAL to be installed, "
209  "but cmake could not find it."));
210  }
211 
212 #endif
213 
214 
215 
216  template void
217  create_coupling_sparsity_pattern_with_exact_intersections(
218  const std::vector<
219  std::tuple<typename dealii::Triangulation<1, 1>::cell_iterator,
220  typename dealii::Triangulation<1, 1>::cell_iterator,
221  dealii::Quadrature<1>>> &intersections_info,
222  const DoFHandler<1, 1> &space_dh,
223  const DoFHandler<1, 1> &immersed_dh,
224  DynamicSparsityPattern &sparsity,
225  const AffineConstraints<double> &constraints,
226  const ComponentMask &space_comps,
227  const ComponentMask &immersed_comps,
228  const AffineConstraints<double> &immersed_constraint);
229 
230  template void
231  create_coupling_sparsity_pattern_with_exact_intersections(
232  const std::vector<
233  std::tuple<typename dealii::Triangulation<2, 2>::cell_iterator,
234  typename dealii::Triangulation<1, 2>::cell_iterator,
235  dealii::Quadrature<2>>> &intersections_info,
236  const DoFHandler<2, 2> &space_dh,
237  const DoFHandler<1, 2> &immersed_dh,
238  DynamicSparsityPattern &sparsity,
239  const AffineConstraints<double> &constraints,
240  const ComponentMask &space_comps,
241  const ComponentMask &immersed_comps,
242  const AffineConstraints<double> &immersed_constraint);
243 
244 
245  template void
246  create_coupling_sparsity_pattern_with_exact_intersections(
247  const std::vector<
248  std::tuple<typename dealii::Triangulation<3, 3>::cell_iterator,
249  typename dealii::Triangulation<1, 3>::cell_iterator,
250  dealii::Quadrature<3>>> &intersections_info,
251  const DoFHandler<3, 3> &space_dh,
252  const DoFHandler<1, 3> &immersed_dh,
253  DynamicSparsityPattern &sparsity,
254  const AffineConstraints<double> &constraints,
255  const ComponentMask &space_comps,
256  const ComponentMask &immersed_comps,
257  const AffineConstraints<double> &immersed_constraint);
258 
259 
260  template void
261  create_coupling_sparsity_pattern_with_exact_intersections(
262  const std::vector<
263  std::tuple<typename dealii::Triangulation<2, 2>::cell_iterator,
264  typename dealii::Triangulation<2, 2>::cell_iterator,
265  dealii::Quadrature<2>>> &intersections_info,
266  const DoFHandler<2, 2> &space_dh,
267  const DoFHandler<2, 2> &immersed_dh,
268  DynamicSparsityPattern &sparsity,
269  const AffineConstraints<double> &constraints,
270  const ComponentMask &space_comps,
271  const ComponentMask &immersed_comps,
272  const AffineConstraints<double> &immersed_constraint);
273 
274  template void
275  create_coupling_sparsity_pattern_with_exact_intersections(
276  const std::vector<
277  std::tuple<typename dealii::Triangulation<3, 3>::cell_iterator,
278  typename dealii::Triangulation<2, 3>::cell_iterator,
279  dealii::Quadrature<3>>> &intersections_info,
280  const DoFHandler<3, 3> &space_dh,
281  const DoFHandler<2, 3> &immersed_dh,
282  DynamicSparsityPattern &sparsity,
283  const AffineConstraints<double> &constraints,
284  const ComponentMask &space_comps,
285  const ComponentMask &immersed_comps,
286  const AffineConstraints<double> &immersed_constraint);
287 
288  template void
289  create_coupling_sparsity_pattern_with_exact_intersections(
290  const std::vector<
291  std::tuple<typename dealii::Triangulation<3, 3>::cell_iterator,
292  typename dealii::Triangulation<3, 3>::cell_iterator,
293  dealii::Quadrature<3>>> &intersections_info,
294  const DoFHandler<3, 3> &space_dh,
295  const DoFHandler<3, 3> &immersed_dh,
296  DynamicSparsityPattern &sparsity,
297  const AffineConstraints<double> &constraints,
298  const ComponentMask &space_comps,
299  const ComponentMask &immersed_comps,
300  const AffineConstraints<double> &immersed_constraint);
301 
302  template void
303  create_coupling_sparsity_pattern_with_exact_intersections<
304  2,
305  1,
306  2,
307  dealii::TrilinosWrappers::SparsityPattern,
308  double>(
309  std::vector<
310  std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
311  dealii::Triangulation<1, 2>::cell_iterator,
312  dealii::Quadrature<2>>,
313  std::allocator<std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
314  dealii::Triangulation<1, 2>::cell_iterator,
315  dealii::Quadrature<2>>>> const &,
316  dealii::DoFHandler<2, 2> const &,
317  dealii::DoFHandler<1, 2> const &,
318  dealii::TrilinosWrappers::SparsityPattern &,
319  dealii::AffineConstraints<double> const &,
320  dealii::ComponentMask const &,
321  dealii::ComponentMask const &,
322  dealii::AffineConstraints<double> const &);
323  template void
324  create_coupling_sparsity_pattern_with_exact_intersections<
325  2,
326  2,
327  2,
328  dealii::TrilinosWrappers::SparsityPattern,
329  double>(
330  std::vector<
331  std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
332  dealii::Triangulation<2, 2>::cell_iterator,
333  dealii::Quadrature<2>>,
334  std::allocator<std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
335  dealii::Triangulation<2, 2>::cell_iterator,
336  dealii::Quadrature<2>>>> const &,
337  dealii::DoFHandler<2, 2> const &,
338  dealii::DoFHandler<2, 2> const &,
339  dealii::TrilinosWrappers::SparsityPattern &,
340  dealii::AffineConstraints<double> const &,
341  dealii::ComponentMask const &,
342  dealii::ComponentMask const &,
343  dealii::AffineConstraints<double> const &);
344  template void
345  create_coupling_sparsity_pattern_with_exact_intersections<
346  3,
347  2,
348  3,
349  dealii::TrilinosWrappers::SparsityPattern,
350  double>(
351  std::vector<
352  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
353  dealii::Triangulation<2, 3>::cell_iterator,
354  dealii::Quadrature<3>>,
355  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
356  dealii::Triangulation<2, 3>::cell_iterator,
357  dealii::Quadrature<3>>>> const &,
358  dealii::DoFHandler<3, 3> const &,
359  dealii::DoFHandler<2, 3> const &,
360  dealii::TrilinosWrappers::SparsityPattern &,
361  dealii::AffineConstraints<double> const &,
362  dealii::ComponentMask const &,
363  dealii::ComponentMask const &,
364  dealii::AffineConstraints<double> const &);
365  template void
366  create_coupling_sparsity_pattern_with_exact_intersections<
367  3,
368  3,
369  3,
370  dealii::TrilinosWrappers::SparsityPattern,
371  double>(
372  std::vector<
373  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
374  dealii::Triangulation<3, 3>::cell_iterator,
375  dealii::Quadrature<3>>,
376  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
377  dealii::Triangulation<3, 3>::cell_iterator,
378  dealii::Quadrature<3>>>> const &,
379  dealii::DoFHandler<3, 3> const &,
380  dealii::DoFHandler<3, 3> const &,
381  dealii::TrilinosWrappers::SparsityPattern &,
382  dealii::AffineConstraints<double> const &,
383  dealii::ComponentMask const &,
384  dealii::ComponentMask const &,
385  dealii::AffineConstraints<double> const &);
386 
387 
388  template void
389  create_coupling_sparsity_pattern_with_exact_intersections<
390  3,
391  2,
392  3,
393  dealii::SparsityPattern,
394  double>(
395  std::vector<
396  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
397  dealii::Triangulation<2, 3>::cell_iterator,
398  dealii::Quadrature<3>>,
399  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
400  dealii::Triangulation<2, 3>::cell_iterator,
401  dealii::Quadrature<3>>>> const &,
402  dealii::DoFHandler<3, 3> const &,
403  dealii::DoFHandler<2, 3> const &,
404  dealii::SparsityPattern &,
405  dealii::AffineConstraints<double> const &,
406  dealii::ComponentMask const &,
407  dealii::ComponentMask const &,
408  dealii::AffineConstraints<double> const &);
409 
410  template void
411  create_coupling_sparsity_pattern_with_exact_intersections<
412  2,
413  2,
414  2,
415  dealii::SparsityPattern,
416  double>(
417  std::vector<
418  std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
419  dealii::Triangulation<2, 2>::cell_iterator,
420  dealii::Quadrature<2>>,
421  std::allocator<std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
422  dealii::Triangulation<2, 2>::cell_iterator,
423  dealii::Quadrature<2>>>> const &,
424  dealii::DoFHandler<2, 2> const &,
425  dealii::DoFHandler<2, 2> const &,
426  dealii::SparsityPattern &,
427  dealii::AffineConstraints<double> const &,
428  dealii::ComponentMask const &,
429  dealii::ComponentMask const &,
430  dealii::AffineConstraints<double> const &);
431 
432  template void
433  create_coupling_sparsity_pattern_with_exact_intersections<
434  1,
435  1,
436  1,
437  dealii::SparsityPattern,
438  double>(
439  std::vector<
440  std::tuple<dealii::Triangulation<1, 1>::cell_iterator,
441  dealii::Triangulation<1, 1>::cell_iterator,
442  dealii::Quadrature<1>>,
443  std::allocator<std::tuple<dealii::Triangulation<1, 1>::cell_iterator,
444  dealii::Triangulation<1, 1>::cell_iterator,
445  dealii::Quadrature<1>>>> const &,
446  dealii::DoFHandler<1, 1> const &,
447  dealii::DoFHandler<1, 1> const &,
448  dealii::SparsityPattern &,
449  dealii::AffineConstraints<double> const &,
450  dealii::ComponentMask const &,
451  dealii::ComponentMask const &,
452  dealii::AffineConstraints<double> const &);
453 
454  template void
455  create_coupling_sparsity_pattern_with_exact_intersections<
456  3,
457  1,
458  3,
459  dealii::SparsityPattern,
460  double>(
461  std::vector<
462  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
463  dealii::Triangulation<1, 3>::cell_iterator,
464  dealii::Quadrature<3>>,
465  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
466  dealii::Triangulation<1, 3>::cell_iterator,
467  dealii::Quadrature<3>>>> const &,
468  dealii::DoFHandler<3, 3> const &,
469  dealii::DoFHandler<1, 3> const &,
470  dealii::SparsityPattern &,
471  dealii::AffineConstraints<double> const &,
472  dealii::ComponentMask const &,
473  dealii::ComponentMask const &,
474  dealii::AffineConstraints<double> const &);
475 
476  template void
477  create_coupling_sparsity_pattern_with_exact_intersections<
478  2,
479  1,
480  2,
481  dealii::SparsityPattern,
482  double>(
483  std::vector<
484  std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
485  dealii::Triangulation<1, 2>::cell_iterator,
486  dealii::Quadrature<2>>,
487  std::allocator<std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
488  dealii::Triangulation<1, 2>::cell_iterator,
489  dealii::Quadrature<2>>>> const &,
490  dealii::DoFHandler<2, 2> const &,
491  dealii::DoFHandler<1, 2> const &,
492  dealii::SparsityPattern &,
493  dealii::AffineConstraints<double> const &,
494  dealii::ComponentMask const &,
495  dealii::ComponentMask const &,
496  dealii::AffineConstraints<double> const &);
497 
498  template void
499  create_coupling_sparsity_pattern_with_exact_intersections<
500  3,
501  3,
502  3,
503  dealii::SparsityPattern,
504  double>(
505  std::vector<
506  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
507  dealii::Triangulation<3, 3>::cell_iterator,
508  dealii::Quadrature<3>>,
509  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
510  dealii::Triangulation<3, 3>::cell_iterator,
511  dealii::Quadrature<3>>>> const &,
512  dealii::DoFHandler<3, 3> const &,
513  dealii::DoFHandler<3, 3> const &,
514  dealii::SparsityPattern &,
515  dealii::AffineConstraints<double> const &,
516  dealii::ComponentMask const &,
517  dealii::ComponentMask const &,
518  dealii::AffineConstraints<double> const &);
519  } // namespace NonMatching
520 } // namespace dealii
AffineConstraints::add_entries_local_to_global
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
ComponentMask::size
unsigned int size() const
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
FiniteElement::n_dofs_per_cell
unsigned int n_dofs_per_cell() const
FiniteElement::system_to_component_index
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
compute_intersections.h
create_coupling_sparsity_pattern_with_exact_intersections.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
AssertDimension
#define AssertDimension(dim1, dim2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
dealii::NonMatching::create_coupling_sparsity_pattern_with_exact_intersections
void create_coupling_sparsity_pattern_with_exact_intersections(const std::vector< std::tuple< typename Triangulation< dim0, spacedim >::cell_iterator, typename Triangulation< dim1, spacedim >::cell_iterator, Quadrature< spacedim >>> &intersections_info, const DoFHandler< dim0, spacedim > &space_dh, const DoFHandler< dim1, spacedim > &immersed_dh, Sparsity &sparsity, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const ComponentMask &space_comps=ComponentMask(), const ComponentMask &immersed_comps=ComponentMask(), const AffineConstraints< number > &immersed_constraints=AffineConstraints< number >())
Create a coupling sparsity pattern of two non-matching, overlapped grids.
Definition: create_coupling_sparsity_pattern_with_exact_intersections.cc:194
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int