Fluid structure interaction suite
assemble_coupling_mass_matrix_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 
17 #include <deal.II/base/config.h>
18 
20 
23 
26 
28 
29 #include <deal.II/fe/mapping.h>
30 #include <deal.II/fe/mapping_q1.h>
31 
34 #include <deal.II/grid/tria.h>
35 
36 #include "lac.h"
37 
38 using namespace dealii;
39 
40 namespace dealii
41 {
42  namespace NonMatching
43  {
44 #if defined DEAL_II_WITH_CGAL || defined DEAL_II_WITH_PARMOONOLITH
45 
46  template <int dim0, int dim1, int spacedim, typename Matrix>
47  void
49  const DoFHandler<dim0, spacedim> &space_dh,
50  const DoFHandler<dim1, spacedim> &immersed_dh,
51  const std::vector<
54  Quadrature<spacedim>>> &cells_and_quads,
55  Matrix &matrix,
56  const AffineConstraints<typename Matrix::value_type> &space_constraints,
57  const ComponentMask &space_comps,
58  const ComponentMask &immersed_comps,
59  const Mapping<dim0, spacedim> &space_mapping,
60  const Mapping<dim1, spacedim> &immersed_mapping,
62  &immersed_constraints)
63  {
64  AssertDimension(matrix.m(), space_dh.n_dofs());
65  AssertDimension(matrix.n(), immersed_dh.n_dofs());
66  Assert(dim1 <= dim0,
67  ExcMessage("This function can only work if dim1<=dim0"));
68  Assert((dynamic_cast<
70  &immersed_dh.get_triangulation()) == nullptr),
71  ExcMessage("The immersed triangulation can only be a "
72  "parallel::shared::triangulation"));
73 
74  const auto &space_fe = space_dh.get_fe();
75  const auto &immersed_fe = immersed_dh.get_fe();
76 
77  const unsigned int n_dofs_per_space_cell = space_fe.n_dofs_per_cell();
78  const unsigned int n_dofs_per_immersed_cell =
79  immersed_fe.n_dofs_per_cell();
80 
81  const unsigned int n_space_fe_components = space_fe.n_components();
82  const unsigned int n_immersed_fe_components = immersed_fe.n_components();
83 
84  FullMatrix<double> local_cell_matrix(n_dofs_per_space_cell,
85  n_dofs_per_immersed_cell);
86  // DoF indices
87  std::vector<types::global_dof_index> local_space_dof_indices(
88  n_dofs_per_space_cell);
89  std::vector<types::global_dof_index> local_immersed_dof_indices(
90  n_dofs_per_immersed_cell);
91 
92  const ComponentMask space_c =
93  (space_comps.size() == 0 ? ComponentMask(n_space_fe_components, true) :
94  space_comps);
95  const ComponentMask immersed_c =
96  (immersed_comps.size() == 0 ?
97  ComponentMask(n_immersed_fe_components, true) :
98  immersed_comps);
99 
100  AssertDimension(space_c.size(), n_space_fe_components);
101  AssertDimension(immersed_c.size(), n_immersed_fe_components);
102 
103  std::vector<unsigned int> space_gtl(n_space_fe_components,
105  std::vector<unsigned int> immersed_gtl(n_immersed_fe_components,
107  for (unsigned int i = 0, j = 0; i < n_space_fe_components; ++i)
108  {
109  if (space_c[i])
110  space_gtl[i] = j++;
111  }
112 
113 
114  for (unsigned int i = 0, j = 0; i < n_immersed_fe_components; ++i)
115  {
116  if (immersed_c[i])
117  immersed_gtl[i] = j++;
118  }
119 
120 
121 
122  Table<2, bool> dof_mask(n_dofs_per_space_cell, n_dofs_per_immersed_cell);
123  dof_mask.fill(false); // start off by assuming they don't couple
124 
125  for (unsigned int i = 0; i < n_dofs_per_space_cell; ++i)
126  {
127  const auto comp_i = space_fe.system_to_component_index(i).first;
128  if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
129  {
130  for (unsigned int j = 0; j < n_dofs_per_immersed_cell; ++j)
131  {
132  const auto comp_j =
133  immersed_fe.system_to_component_index(j).first;
134  if (immersed_gtl[comp_j] == space_gtl[comp_i])
135  {
136  dof_mask(i, j) = true;
137  }
138  }
139  }
140  }
141 
142  // const bool dof_mask_is_active =
143  // dof_mask.n_rows() == n_dofs_per_space_cell &&
144  // dof_mask.n_cols() == n_dofs_per_immersed_cell;
145 
146 
147  // Loop over vector of tuples, and gather everything together
148  for (const auto &infos : cells_and_quads)
149  {
150  const auto &[first_cell, second_cell, quad_formula] = infos;
151 
152  local_cell_matrix = typename Matrix::value_type();
153 
154  const unsigned int n_quad_pts = quad_formula.size();
155  const auto &real_qpts = quad_formula.get_points();
156  std::vector<Point<dim0>> ref_pts_space(n_quad_pts);
157  std::vector<Point<dim1>> ref_pts_immersed(n_quad_pts);
158 
159  space_mapping.transform_points_real_to_unit_cell(first_cell,
160  real_qpts,
161  ref_pts_space);
162  immersed_mapping.transform_points_real_to_unit_cell(second_cell,
163  real_qpts,
164  ref_pts_immersed);
165  const auto &JxW = quad_formula.get_weights();
166  for (unsigned int q = 0; q < n_quad_pts; ++q)
167  {
168  for (unsigned int i = 0; i < n_dofs_per_space_cell; ++i)
169  {
170  const unsigned int comp_i =
171  space_dh.get_fe().system_to_component_index(i).first;
172  if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
173  {
174  for (unsigned int j = 0; j < n_dofs_per_immersed_cell;
175  ++j)
176  {
177  const unsigned int comp_j =
178  immersed_dh.get_fe()
180  .first;
181  if (space_gtl[comp_i] == immersed_gtl[comp_j])
182  {
183  local_cell_matrix(i, j) +=
184  space_fe.shape_value(i, ref_pts_space[q]) *
185  immersed_fe.shape_value(j,
186  ref_pts_immersed[q]) *
187  JxW[q];
188  }
189  }
190  }
191  }
192  }
193  typename DoFHandler<dim0, spacedim>::cell_iterator space_cell_dh(
194  *first_cell, &space_dh);
195  typename DoFHandler<dim1, spacedim>::cell_iterator immersed_cell_dh(
196  *second_cell, &immersed_dh);
197  space_cell_dh->get_dof_indices(local_space_dof_indices);
198  immersed_cell_dh->get_dof_indices(local_immersed_dof_indices);
199 
200 
201 
202  // if (dof_mask_is_active)
203  // {
204  // for (unsigned int i = 0; i < n_dofs_per_space_cell; ++i)
205  // {
206  // const unsigned int comp_i =
207  // space_dh.get_fe().system_to_component_index(i).first;
208  // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
209  // {
210  // for (unsigned int j = 0; j < n_dofs_per_immersed_cell;
211  // ++j)
212  // {
213  // const unsigned int comp_j =
214  // immersed_dh.get_fe()
215  // .system_to_component_index(j)
216  // .first;
217  // if (space_gtl[comp_i] == immersed_gtl[comp_j])
218  // {
219  // space_constraints.distribute_local_to_global(
220  // local_cell_matrix,
221  // {local_space_dof_indices[i]},
222  // immersed_constraints,
223  // {local_immersed_dof_indices[j]},
224  // matrix);
225  // }
226  // }
227  // }
228  // }
229  // }
230  // else
231  // {
232  space_constraints.distribute_local_to_global(
233  local_cell_matrix,
234  local_space_dof_indices,
235  immersed_constraints,
236  local_immersed_dof_indices,
237  matrix);
238  // }
239  }
240  matrix.compress(VectorOperation::add);
241  }
242 
243 
244 
245 #else
246 
247 
248 
249  template <int dim0, int dim1, int spacedim, typename Matrix>
250  void
254  const std::vector<
258  Matrix &,
260  const ComponentMask &,
261  const ComponentMask &,
262  const Mapping<dim0, spacedim> &,
263  const Mapping<dim1, spacedim> &,
265  {
266  Assert(false,
267  ExcMessage(
268  "This function needs CGAL or PARMOONOLITH to be installed, "
269  "but cmake could not either."));
270  }
271 
272 
273 
274 #endif
275 
276 
277 
278  template void
280  const DoFHandler<1, 1> &,
281  const DoFHandler<1, 1> &,
282  const std::vector<
283  std::tuple<typename dealii::Triangulation<1, 1>::cell_iterator,
284  typename dealii::Triangulation<1, 1>::cell_iterator,
285  dealii::Quadrature<1>>> &,
286  dealii::SparseMatrix<double> &,
287  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &,
288  const ComponentMask &,
289  const ComponentMask &,
290  const Mapping<1, 1> &space_mapping,
291  const Mapping<1, 1> &immersed_mapping,
292  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &);
293 
294  template void
296  const DoFHandler<2, 2> &,
297  const DoFHandler<1, 2> &,
298  const std::vector<
299  std::tuple<typename dealii::Triangulation<2, 2>::cell_iterator,
300  typename dealii::Triangulation<1, 2>::cell_iterator,
301  dealii::Quadrature<2>>> &,
302  dealii::SparseMatrix<double> &,
303  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &,
304  const ComponentMask &,
305  const ComponentMask &,
306  const Mapping<2, 2> &space_mapping,
307  const Mapping<1, 2> &immersed_mapping,
308  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &);
309 
310  template void
312  const DoFHandler<2, 2> &,
313  const DoFHandler<2, 2> &,
314  const std::vector<
315  std::tuple<typename dealii::Triangulation<2, 2>::cell_iterator,
316  typename dealii::Triangulation<2, 2>::cell_iterator,
317  dealii::Quadrature<2>>> &,
318  dealii::SparseMatrix<double> &,
319  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &,
320  const ComponentMask &,
321  const ComponentMask &,
322  const Mapping<2, 2> &space_mapping,
323  const Mapping<2, 2> &immersed_mapping,
324  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &);
325 
326 
327 
328  template void
330  const DoFHandler<3, 3> &,
331  const DoFHandler<1, 3> &,
332  const std::vector<
333  std::tuple<typename dealii::Triangulation<3, 3>::cell_iterator,
334  typename dealii::Triangulation<1, 3>::cell_iterator,
335  dealii::Quadrature<3>>> &,
336  dealii::SparseMatrix<double> &,
337  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &,
338  const ComponentMask &,
339  const ComponentMask &,
340  const Mapping<3, 3> &space_mapping,
341  const Mapping<1, 3> &immersed_mapping,
342  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &);
343 
344  template void
346  const DoFHandler<3, 3> &,
347  const DoFHandler<2, 3> &,
348  const std::vector<
349  std::tuple<typename dealii::Triangulation<3, 3>::cell_iterator,
350  typename dealii::Triangulation<2, 3>::cell_iterator,
351  dealii::Quadrature<3>>> &,
352  dealii::SparseMatrix<double> &,
353  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &,
354  const ComponentMask &,
355  const ComponentMask &,
356  const Mapping<3, 3> &space_mapping,
357  const Mapping<2, 3> &immersed_mapping,
358  const AffineConstraints<dealii::SparseMatrix<double>::value_type> &);
359 
360  template void
362  const DoFHandler<3, 3> &,
363  const DoFHandler<3, 3> &,
364  const std::vector<std::tuple<typename Triangulation<3, 3>::cell_iterator,
366  Quadrature<3>>> &,
369  const ComponentMask &,
370  const ComponentMask &,
371  const Mapping<3, 3> &space_mapping,
372  const Mapping<3, 3> &immersed_mapping,
374 
375  template void
377  2,
378  1,
379  2,
380  dealii::TrilinosWrappers::SparseMatrix>(
381  dealii::DoFHandler<2, 2> const &,
382  dealii::DoFHandler<1, 2> const &,
383  std::vector<
384  std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
385  dealii::Triangulation<1, 2>::cell_iterator,
386  dealii::Quadrature<2>>,
387  std::allocator<std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
388  dealii::Triangulation<1, 2>::cell_iterator,
389  dealii::Quadrature<2>>>> const &,
390  dealii::TrilinosWrappers::SparseMatrix &,
391  dealii::AffineConstraints<
392  dealii::TrilinosWrappers::SparseMatrix::value_type> const &,
393  dealii::ComponentMask const &,
394  dealii::ComponentMask const &,
395  dealii::Mapping<2, 2> const &,
396  dealii::Mapping<1, 2> const &,
397  dealii::AffineConstraints<
398  dealii::TrilinosWrappers::SparseMatrix::value_type> const &);
399 
400  template void
402  2,
403  2,
404  2,
405  dealii::TrilinosWrappers::SparseMatrix>(
406  dealii::DoFHandler<2, 2> const &,
407  dealii::DoFHandler<2, 2> const &,
408  std::vector<
409  std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
410  dealii::Triangulation<2, 2>::cell_iterator,
411  dealii::Quadrature<2>>,
412  std::allocator<std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
413  dealii::Triangulation<2, 2>::cell_iterator,
414  dealii::Quadrature<2>>>> const &,
415  dealii::TrilinosWrappers::SparseMatrix &,
416  dealii::AffineConstraints<
417  dealii::TrilinosWrappers::SparseMatrix::value_type> const &,
418  dealii::ComponentMask const &,
419  dealii::ComponentMask const &,
420  dealii::Mapping<2, 2> const &,
421  dealii::Mapping<2, 2> const &,
422  dealii::AffineConstraints<
423  dealii::TrilinosWrappers::SparseMatrix::value_type> const &);
424 
425  template void
427  3,
428  2,
429  3,
430  dealii::TrilinosWrappers::SparseMatrix>(
431  dealii::DoFHandler<3, 3> const &,
432  dealii::DoFHandler<2, 3> const &,
433  std::vector<
434  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
435  dealii::Triangulation<2, 3>::cell_iterator,
436  dealii::Quadrature<3>>,
437  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
438  dealii::Triangulation<2, 3>::cell_iterator,
439  dealii::Quadrature<3>>>> const &,
440  dealii::TrilinosWrappers::SparseMatrix &,
441  dealii::AffineConstraints<
442  dealii::TrilinosWrappers::SparseMatrix::value_type> const &,
443  dealii::ComponentMask const &,
444  dealii::ComponentMask const &,
445  dealii::Mapping<3, 3> const &,
446  dealii::Mapping<2, 3> const &,
447  dealii::AffineConstraints<
448  dealii::TrilinosWrappers::SparseMatrix::value_type> const &);
449 
450  template void
452  3,
453  3,
454  3,
455  dealii::TrilinosWrappers::SparseMatrix>(
456  dealii::DoFHandler<3, 3> const &,
457  dealii::DoFHandler<3, 3> const &,
458  std::vector<
459  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
460  dealii::Triangulation<3, 3>::cell_iterator,
461  dealii::Quadrature<3>>,
462  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
463  dealii::Triangulation<3, 3>::cell_iterator,
464  dealii::Quadrature<3>>>> const &,
465  dealii::TrilinosWrappers::SparseMatrix &,
466  dealii::AffineConstraints<
467  dealii::TrilinosWrappers::SparseMatrix::value_type> const &,
468  dealii::ComponentMask const &,
469  dealii::ComponentMask const &,
470  dealii::Mapping<3, 3> const &,
471  dealii::Mapping<3, 3> const &,
472  dealii::AffineConstraints<
473  dealii::TrilinosWrappers::SparseMatrix::value_type> const &);
474 
475  template void
477  2,
478  1,
479  2,
480  dealii::PETScWrappers::MPI::SparseMatrix>(
481  dealii::DoFHandler<2, 2> const &,
482  dealii::DoFHandler<1, 2> const &,
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::PETScWrappers::MPI::SparseMatrix &,
491  dealii::AffineConstraints<
492  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &,
493  dealii::ComponentMask const &,
494  dealii::ComponentMask const &,
495  dealii::Mapping<2, 2> const &,
496  dealii::Mapping<1, 2> const &,
497  dealii::AffineConstraints<
498  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &);
499 
500  template void
502  2,
503  2,
504  2,
505  dealii::PETScWrappers::MPI::SparseMatrix>(
506  dealii::DoFHandler<2, 2> const &,
507  dealii::DoFHandler<2, 2> const &,
508  std::vector<
509  std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
510  dealii::Triangulation<2, 2>::cell_iterator,
511  dealii::Quadrature<2>>,
512  std::allocator<std::tuple<dealii::Triangulation<2, 2>::cell_iterator,
513  dealii::Triangulation<2, 2>::cell_iterator,
514  dealii::Quadrature<2>>>> const &,
515  dealii::PETScWrappers::MPI::SparseMatrix &,
516  dealii::AffineConstraints<
517  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &,
518  dealii::ComponentMask const &,
519  dealii::ComponentMask const &,
520  dealii::Mapping<2, 2> const &,
521  dealii::Mapping<2, 2> const &,
522  dealii::AffineConstraints<
523  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &);
524 
525  template void
527  3,
528  2,
529  3,
530  dealii::PETScWrappers::MPI::SparseMatrix>(
531  dealii::DoFHandler<3, 3> const &,
532  dealii::DoFHandler<2, 3> const &,
533  std::vector<
534  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
535  dealii::Triangulation<2, 3>::cell_iterator,
536  dealii::Quadrature<3>>,
537  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
538  dealii::Triangulation<2, 3>::cell_iterator,
539  dealii::Quadrature<3>>>> const &,
540  dealii::PETScWrappers::MPI::SparseMatrix &,
541  dealii::AffineConstraints<
542  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &,
543  dealii::ComponentMask const &,
544  dealii::ComponentMask const &,
545  dealii::Mapping<3, 3> const &,
546  dealii::Mapping<2, 3> const &,
547  dealii::AffineConstraints<
548  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &);
549 
550  template void
552  3,
553  3,
554  3,
555  dealii::PETScWrappers::MPI::SparseMatrix>(
556  dealii::DoFHandler<3, 3> const &,
557  dealii::DoFHandler<3, 3> const &,
558  std::vector<
559  std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
560  dealii::Triangulation<3, 3>::cell_iterator,
561  dealii::Quadrature<3>>,
562  std::allocator<std::tuple<dealii::Triangulation<3, 3>::cell_iterator,
563  dealii::Triangulation<3, 3>::cell_iterator,
564  dealii::Quadrature<3>>>> const &,
565  dealii::PETScWrappers::MPI::SparseMatrix &,
566  dealii::AffineConstraints<
567  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &,
568  dealii::ComponentMask const &,
569  dealii::ComponentMask const &,
570  dealii::Mapping<3, 3> const &,
571  dealii::Mapping<3, 3> const &,
572  dealii::AffineConstraints<
573  dealii::PETScWrappers::MPI::SparseMatrix::value_type> const &);
574  } // namespace NonMatching
575 } // namespace dealii
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
unsigned int size() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
unsigned int n_dofs_per_cell() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
virtual void transform_points_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< spacedim > > &real_points, const ArrayView< Point< dim > > &unit_points) const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
void assemble_coupling_mass_matrix_with_exact_intersections(const DoFHandler< dim0, spacedim > &, const DoFHandler< dim1, spacedim > &, const std::vector< std::tuple< typename Triangulation< dim0, spacedim >::cell_iterator, typename Triangulation< dim1, spacedim >::cell_iterator, Quadrature< spacedim >>> &, Matrix &matrix, const AffineConstraints< typename Matrix::value_type > &, const ComponentMask &, const ComponentMask &, const Mapping< dim0, spacedim > &, const Mapping< dim1, spacedim > &, const AffineConstraints< typename Matrix::value_type > &)
Create the coupling mass matrix for non-matching, overlapping grids in an "exact" way,...
static const unsigned int invalid_unsigned_int