Fluid structure interaction suite
compute_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 
19 #include <deal.II/base/function_lib.h>
20 #include <deal.II/base/quadrature_lib.h>
21 
22 #include <deal.II/fe/mapping_q1.h>
23 
24 #include <deal.II/grid/grid_tools.h>
25 #include <deal.II/grid/grid_tools_cache.h>
26 
27 #include <deal.II/cgal/intersections.h>
28 
29 #include <set>
30 #include <tuple>
31 #include <vector>
32 
33 #include "moonolith_tools.h"
34 
35 using namespace dealii;
36 
37 #if defined DEAL_II_PREFER_CGAL_OVER_PARMOONOLITH
38 
39 namespace dealii
40 {
41  namespace NonMatching
42  {
47  template <int dim0, int dim1, int spacedim>
48  dealii::Quadrature<spacedim>
49  compute_cell_intersection(
50  const typename Triangulation<dim0, spacedim>::cell_iterator &cell0,
51  const typename Triangulation<dim1, spacedim>::cell_iterator &cell1,
52  const unsigned int degree,
53  const Mapping<dim0, spacedim> &mapping0,
54  const Mapping<dim1, spacedim> &mapping1,
55  const double tol = 1e-9)
56  {
57  if constexpr (dim0 == 1 && dim1 == 1)
58  {
59  (void)cell0;
60  (void)cell1;
61  (void)degree;
62  (void)mapping0;
63  (void)mapping1;
64  (void)tol;
65  AssertThrow(false, ExcNotImplemented());
66  return dealii::Quadrature<spacedim>();
67  }
68  else
69  {
70  const auto &vec_arrays =
71  ::CGALWrappers::compute_intersection_of_cells(
72  cell0, cell1, mapping0, mapping1, tol);
73  return QGaussSimplex<dim1>(degree).mapped_quadrature(vec_arrays);
74  }
75  }
76 
77 
78 
79 #elif defined DEAL_II_WITH_PARMOONOLITH
80 
81 namespace dealii
82 {
83  namespace NonMatching
84  {
85  template <int dim0, int dim1, int spacedim>
86  dealii::Quadrature<spacedim>
87  compute_cell_intersection(
88  const typename Triangulation<dim0, spacedim>::cell_iterator &cell0,
89  const typename Triangulation<dim1, spacedim>::cell_iterator &cell1,
90  const unsigned int degree,
91  const Mapping<dim0, spacedim> &mapping0,
92  const Mapping<dim1, spacedim> &mapping1)
93  {
94  if constexpr ((dim0 == 1 && dim1 == 3) || (dim0 == 3 && dim1 == 1) ||
95  (dim0 == 1 && dim1 == 1))
96  {
97  (void)cell0;
98  (void)cell1;
99  (void)degree;
100  (void)mapping0;
101  (void)mapping1;
102  AssertThrow(false, ExcNotImplemented());
103  return dealii::Quadrature<spacedim>();
104  }
105  else
106  {
107  return moonolith::compute_intersection(
108  cell0, cell1, degree, mapping0, mapping1);
109  }
110  }
111 #else
112 
113 namespace dealii
114 {
115  namespace NonMatching
116  {
117  template <int dim0, int dim1, int spacedim>
118  Quadrature<spacedim>
119  compute_intersection(
120  const typename Triangulation<dim0, spacedim>::cell_iterator &,
121  const typename Triangulation<dim1, spacedim>::cell_iterator &,
122  const unsigned int,
123  const Mapping<dim0, spacedim> &,
124  const Mapping<dim1, spacedim> &)
125  {
126  Assert(false,
127  ExcMessage(
128  "This function needs CGAL or PARMOONOLITH to be installed, "
129  "but cmake could not find any of them."));
130  return Quadrature<spacedim>();
131  }
132 
133 #endif
134 
135  template <int dim0, int dim1, int spacedim>
136  std::vector<
137  std::tuple<typename Triangulation<dim0, spacedim>::cell_iterator,
138  typename Triangulation<dim1, spacedim>::cell_iterator,
139  Quadrature<spacedim>>>
140  compute_intersection(const GridTools::Cache<dim0, spacedim> &space_cache,
141  const GridTools::Cache<dim1, spacedim> &immersed_cache,
142  const unsigned int degree,
143  const double tol)
144  {
145  Assert(degree >= 1, ExcMessage("degree cannot be less than 1"));
146 
147  std::vector<
148  std::tuple<typename Triangulation<dim0, spacedim>::cell_iterator,
149  typename Triangulation<dim1, spacedim>::cell_iterator,
150  Quadrature<spacedim>>>
151  cells_with_quads;
152 
153 
154  const auto &space_tree =
155  space_cache.get_locally_owned_cell_bounding_boxes_rtree();
156 
157  // The immersed tree *must* contain all cells, also the non-locally owned
158  // ones.
159  const auto &immersed_tree =
160  immersed_cache.get_cell_bounding_boxes_rtree();
161 
162  // references to triangulations' info (cp cstrs marked as delete)
163  const auto &mapping0 = space_cache.get_mapping();
164  const auto &mapping1 = immersed_cache.get_mapping();
165  namespace bgi = boost::geometry::index;
166  // Whenever the BB space_cell intersects the BB of an embedded cell,
167  // store the space_cell in the set of intersected_cells
168  for (const auto &[immersed_box, immersed_cell] : immersed_tree)
169  {
170  for (const auto &[space_box, space_cell] :
171  space_tree |
172  bgi::adaptors::queried(bgi::intersects(immersed_box)))
173  {
174  const auto &test_intersection =
175  compute_cell_intersection<dim0, dim1, spacedim>(
176  space_cell, immersed_cell, degree, mapping0, mapping1);
177 
178  // if (test_intersection.get_points().size() !=
179  const auto &weights = test_intersection.get_weights();
180  const double area =
181  std::accumulate(weights.begin(), weights.end(), 0.0);
182  if (area > tol) // non-trivial intersection
183  {
184  cells_with_quads.push_back(std::make_tuple(
185  space_cell, immersed_cell, test_intersection));
186  }
187  }
188  }
189 
190  return cells_with_quads;
191  }
192 
193 
194  template Quadrature<1>
195  compute_cell_intersection(const Triangulation<1, 1>::cell_iterator &,
196  const Triangulation<1, 1>::cell_iterator &,
197  const unsigned int,
198  const Mapping<1, 1> &,
199  const Mapping<1, 1> &,
200  const double);
201 
202 
203  template Quadrature<2>
204  compute_cell_intersection(const Triangulation<2, 2>::cell_iterator &,
205  const Triangulation<1, 2>::cell_iterator &,
206  const unsigned int,
207  const Mapping<2, 2> &,
208  const Mapping<1, 2> &,
209  const double);
210 
211  template Quadrature<2>
212  compute_cell_intersection(const Triangulation<2, 2>::cell_iterator &,
213  const Triangulation<2, 2>::cell_iterator &,
214  const unsigned int,
215  const Mapping<2, 2> &,
216  const Mapping<2, 2> &,
217  const double);
218 
219 
220  template Quadrature<3>
221  compute_cell_intersection(const Triangulation<3, 3>::cell_iterator &,
222  const Triangulation<2, 3>::cell_iterator &,
223  const unsigned int,
224  const Mapping<3, 3> &,
225  const Mapping<2, 3> &,
226  const double);
227 
228  template Quadrature<3>
229  compute_cell_intersection(const Triangulation<3, 3>::cell_iterator &,
230  const Triangulation<3, 3>::cell_iterator &,
231  const unsigned int,
232  const Mapping<3, 3> &,
233  const Mapping<3, 3> &,
234  const double);
235 
236  template std::vector<
237  std::tuple<typename dealii::Triangulation<1, 1>::cell_iterator,
238  typename dealii::Triangulation<1, 1>::cell_iterator,
239  Quadrature<1>>>
240  NonMatching::compute_intersection(
241  const GridTools::Cache<1, 1> &space_cache,
242  const GridTools::Cache<1, 1> &immersed_cache,
243  const unsigned int degree,
244  const double tol);
245 
246 
247 
248  template std::vector<
249  std::tuple<typename dealii::Triangulation<3, 3>::cell_iterator,
250  typename dealii::Triangulation<1, 3>::cell_iterator,
251  Quadrature<3>>>
252  NonMatching::compute_intersection(
253  const GridTools::Cache<3, 3> &space_cache,
254  const GridTools::Cache<1, 3> &immersed_cache,
255  const unsigned int degree,
256  const double tol);
257 
258 
259  template std::vector<
260  std::tuple<typename dealii::Triangulation<2, 2>::cell_iterator,
261  typename dealii::Triangulation<1, 2>::cell_iterator,
262  Quadrature<2>>>
263  NonMatching::compute_intersection(
264  const GridTools::Cache<2, 2> &space_cache,
265  const GridTools::Cache<1, 2> &immersed_cache,
266  const unsigned int degree,
267  const double tol);
268 
269 
270 
271  template std::vector<std::tuple<typename Triangulation<2, 2>::cell_iterator,
272  typename Triangulation<2, 2>::cell_iterator,
273  Quadrature<2>>>
274  NonMatching::compute_intersection(
275  const GridTools::Cache<2, 2> &space_cache,
276  const GridTools::Cache<2, 2> &immersed_cache,
277  const unsigned int degree,
278  const double tol);
279 
280 
281 
282  template std::vector<std::tuple<typename Triangulation<3, 3>::cell_iterator,
283  typename Triangulation<2, 3>::cell_iterator,
284  Quadrature<3>>>
285  NonMatching::compute_intersection(
286  const GridTools::Cache<3, 3> &space_cache,
287  const GridTools::Cache<2, 3> &immersed_cache,
288  const unsigned int degree,
289  const double tol);
290 
291 
292 
293  template std::vector<std::tuple<typename Triangulation<3, 3>::cell_iterator,
294  typename Triangulation<3, 3>::cell_iterator,
295  Quadrature<3>>>
296  NonMatching::compute_intersection(
297  const GridTools::Cache<3, 3> &space_cache,
298  const GridTools::Cache<3, 3> &immersed_cache,
299  const unsigned int degree,
300  const double tol);
301  }
302 }
GridTools::Cache::get_locally_owned_cell_bounding_boxes_rtree
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_locally_owned_cell_bounding_boxes_rtree() const
GridTools::Cache::get_mapping
const Mapping< dim, spacedim > & get_mapping() const
GridTools::Cache::get_cell_bounding_boxes_rtree
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_cell_bounding_boxes_rtree() const
QGaussSimplex::mapped_quadrature
Quadrature< spacedim > mapped_quadrature(const std::vector< std::array< Point< spacedim >, dim+1 > > &simplices) const
function_lib.h
grid_tools.h
grid_tools_cache.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
mapping_q1.h
moonolith_tools.h
dealii::NonMatching::compute_intersection
Quadrature< spacedim > compute_intersection(const typename Triangulation< dim0, spacedim >::cell_iterator &cell0, const typename Triangulation< dim1, spacedim >::cell_iterator &cell1, const unsigned int degree, const Mapping< dim0, spacedim > &mapping0=(ReferenceCells::get_hypercube< dim0 >() .template get_default_linear_mapping< dim0, spacedim >()), const Mapping< dim1, spacedim > &mapping1=(ReferenceCells::get_hypercube< dim1 >() .template get_default_linear_mapping< dim1, spacedim >()))
Intersect cell0 and cell1 and construct a Quadrature<spacedim> of degree degreeover the intersection,...
Definition: compute_intersections.cc:119
dealii::NonMatching::compute_cell_intersection
template Quadrature< 1 > compute_cell_intersection(const Triangulation< 1, 1 >::cell_iterator &, const Triangulation< 1, 1 >::cell_iterator &, const unsigned int, const Mapping< 1, 1 > &, const Mapping< 1, 1 > &, const double)
quadrature_lib.h