Reference documentation for deal.II version GIT 35969cdc9b 2023-12-09 01:10:02+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\}}\)
mesh_classifier.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 - 2022 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 
17 
19 
20 #include "deal.II/fe/fe_q_iso_q1.h"
22 #include <deal.II/fe/fe_q.h>
23 #include <deal.II/fe/fe_values.h>
24 
34 #include <deal.II/lac/vector.h>
36 
38 
39 #include <algorithm>
40 
42 
43 namespace NonMatching
44 {
45  namespace internal
46  {
47  namespace MeshClassifierImplementation
48  {
51  "The Triangulation has not been classified. You need to call the "
52  "reclassify()-function before using this function.");
53 
56  "The incoming cell does not belong to the triangulation passed to "
57  "the constructor.");
58 
64  template <typename VectorType>
66  location_from_dof_signs(const VectorType &local_levelset_values)
67  {
68  const auto min_max_element =
69  std::minmax_element(local_levelset_values.begin(),
70  local_levelset_values.end());
71 
72  if (*min_max_element.second < 0)
74  if (0 < *min_max_element.first)
76 
78  }
79 
80 
81 
86  template <int dim, typename VectorType>
88  {
89  public:
94  const VectorType &level_set);
95 
99  const hp::FECollection<dim> &
100  get_fe_collection() const override;
101 
106  unsigned int
108  &cell) const override;
109 
114  void
116  const typename Triangulation<dim>::active_cell_iterator &cell,
117  const unsigned int face_index,
118  Vector<double> &local_levelset_values) override;
119 
120  private:
125 
131  };
132 
133 
134 
135  template <int dim, typename VectorType>
137  const DoFHandler<dim> &dof_handler,
138  const VectorType &level_set)
139  : dof_handler(&dof_handler)
140  , level_set(&level_set)
141  {}
142 
143 
144 
145  template <int dim, typename VectorType>
146  const hp::FECollection<dim> &
148  {
149  return dof_handler->get_fe_collection();
150  }
151 
152 
153 
154  template <int dim, typename VectorType>
155  void
157  const typename Triangulation<dim>::active_cell_iterator &cell,
158  const unsigned int face_index,
159  Vector<double> &local_levelset_values)
160  {
161  const auto cell_with_dofs = cell->as_dof_handler_iterator(*dof_handler);
162 
163  const unsigned int n_dofs_per_face =
164  dof_handler->get_fe().n_dofs_per_face();
165  std::vector<types::global_dof_index> dof_indices(n_dofs_per_face);
166  cell_with_dofs->face(face_index)->get_dof_indices(dof_indices);
167 
168  local_levelset_values.reinit(dof_indices.size());
169 
170  for (unsigned int i = 0; i < dof_indices.size(); i++)
171  local_levelset_values[i] =
173  dof_indices[i]);
174  }
175 
176 
177 
178  template <int dim, typename VectorType>
179  unsigned int
181  const typename Triangulation<dim>::active_cell_iterator &cell) const
182  {
183  const auto cell_with_dofs = cell->as_dof_handler_iterator(*dof_handler);
184 
185  return cell_with_dofs->active_fe_index();
186  }
187 
188 
193  template <int dim>
195  {
196  public:
202  const FiniteElement<dim> &element);
203 
208  const hp::FECollection<dim> &
209  get_fe_collection() const override;
210 
215  unsigned int
217  &cell) const override;
218 
223  void
225  const typename Triangulation<dim>::active_cell_iterator &cell,
226  const unsigned int face_index,
227  Vector<double> &local_levelset_values) override;
228 
229  private:
234 
240 
246  };
247 
248 
249 
250  template <int dim>
252  const Function<dim> &level_set,
253  const FiniteElement<dim> &element)
254  : level_set(&level_set)
255  , fe_collection(element)
256  , fe_face_values(element,
257  Quadrature<dim - 1>(
258  element.get_unit_face_support_points()),
260  {}
261 
262 
263 
264  template <int dim>
265  void
267  const typename Triangulation<dim>::active_cell_iterator &cell,
268  const unsigned int face_index,
269  Vector<double> &local_levelset_values)
270  {
271  AssertDimension(local_levelset_values.size(),
272  fe_face_values.n_quadrature_points);
273 
274  fe_face_values.reinit(cell, face_index);
275  const std::vector<Point<dim>> &points =
276  fe_face_values.get_quadrature_points();
277 
278  for (unsigned int i = 0; i < points.size(); i++)
279  local_levelset_values[i] = level_set->value(points[i]);
280  }
281 
282 
283 
284  template <int dim>
285  const hp::FECollection<dim> &
287  {
288  return fe_collection;
289  }
290 
291 
292 
293  template <int dim>
294  unsigned int
296  const typename Triangulation<dim>::active_cell_iterator &) const
297  {
298  return 0;
299  }
300  } // namespace MeshClassifierImplementation
301  } // namespace internal
302 
303 
304 
305  template <int dim>
306  template <typename VectorType>
308  const VectorType &level_set)
309  : triangulation(&dof_handler.get_triangulation())
310  , level_set_description(
311  std::make_unique<internal::MeshClassifierImplementation::
312  DiscreteLevelSetDescription<dim, VectorType>>(
313  dof_handler,
314  level_set))
315  {
316 #ifdef DEAL_II_WITH_LAPACK
317  const hp::FECollection<dim> &fe_collection =
318  dof_handler.get_fe_collection();
319  for (unsigned int i = 0; i < fe_collection.size(); i++)
320  {
321  // The level set function must be scalar.
322  AssertDimension(fe_collection[i].n_components(), 1);
323 
324  Assert(fe_collection[i].has_face_support_points(),
325  ExcMessage(
326  "The elements in the FECollection of the incoming DoFHandler "
327  "must have face support points."));
328  }
329 #else
330  AssertThrow(false, ExcNeedsLAPACK());
331 #endif
332  }
333 
334 
335 
336  template <int dim>
338  const Function<dim> &level_set,
339  const FiniteElement<dim> &element)
341  , level_set_description(
342  std::make_unique<internal::MeshClassifierImplementation::
343  AnalyticLevelSetDescription<dim>>(level_set,
344  element))
345  {
346  // The level set function must be scalar.
347  AssertDimension(level_set.n_components, 1);
348  AssertDimension(element.n_components(), 1);
349  }
350 
351 
352 
353  template <int dim>
354  void
356  {
357  initialize();
358  cell_locations.assign(triangulation->n_active_cells(),
360  face_locations.assign(triangulation->n_raw_faces(),
362 
363  // Loop over all cells and determine the location of all non artificial
364  // cells and faces.
365  for (const auto &cell : triangulation->active_cell_iterators())
366  if (!cell->is_artificial())
367  {
368  const LocationToLevelSet face0_location =
369  determine_face_location_to_levelset(cell, 0);
370 
371  face_locations[cell->face(0)->index()] = face0_location;
372  LocationToLevelSet cell_location = face0_location;
373 
374  for (unsigned int f = 1; f < GeometryInfo<dim>::faces_per_cell; ++f)
375  {
376  const LocationToLevelSet face_location =
377  determine_face_location_to_levelset(cell, f);
378 
379  face_locations[cell->face(f)->index()] = face_location;
380 
381  if (face_location != face0_location)
382  cell_location = LocationToLevelSet::intersected;
383  }
384  cell_locations[cell->active_cell_index()] = cell_location;
385  }
386  }
387 
388 
389 
390  template <int dim>
393  const typename Triangulation<dim>::active_cell_iterator &cell,
394  const unsigned int face_index)
395  {
396  // The location of the face might already be computed on the neighboring
397  // cell. If this is the case we just return the value.
398  const LocationToLevelSet location =
399  face_locations.at(cell->face(face_index)->index());
400  if (location != LocationToLevelSet::unassigned)
401  return location;
402 
403  // Determine the location by changing basis to FE_Bernstein and checking
404  // the signs of the dofs.
405  const unsigned int fe_index = level_set_description->active_fe_index(cell);
406  const unsigned int n_local_dofs =
407  lagrange_to_bernstein_face[fe_index][face_index].m();
408 
409  Vector<double> local_levelset_values(n_local_dofs);
410  level_set_description->get_local_level_set_values(cell,
411  face_index,
412  local_levelset_values);
413 
414  const FiniteElement<dim> &fe =
415  level_set_description->get_fe_collection()[fe_index];
416 
417  const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
418  dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe);
419 
420  const FE_Poly<dim> *fe_poly = dynamic_cast<const FE_Poly<dim> *>(&fe);
421 
422  const bool is_linear = fe_q_iso_q1 != nullptr ||
423  (fe_poly != nullptr && fe_poly->get_degree() == 1);
424 
425  // shortcut for linear elements
426  if (is_linear)
427  {
429  local_levelset_values);
430  }
431 
432  lagrange_to_bernstein_face[fe_index][face_index].solve(
433  local_levelset_values);
434 
436  local_levelset_values);
437  }
438 
439 
440 
441  template <int dim>
444  const typename Triangulation<dim>::cell_iterator &cell) const
445  {
446  Assert(cell_locations.size() == triangulation->n_active_cells(),
448  Assert(&cell->get_triangulation() == triangulation,
450 
451  return cell_locations.at(cell->active_cell_index());
452  }
453 
454 
455 
456  template <int dim>
459  const typename Triangulation<dim>::cell_iterator &cell,
460  const unsigned int face_index) const
461  {
463  Assert(face_locations.size() == triangulation->n_raw_faces(),
465  Assert(&cell->get_triangulation() == triangulation,
467 
468  return face_locations.at(cell->face(face_index)->index());
469  }
470 
471 
472 
473  template <int dim>
474  void
476  {
477  const hp::FECollection<dim> &fe_collection =
478  level_set_description->get_fe_collection();
479 
480  // The level set function must be scalar.
481  AssertDimension(fe_collection.n_components(), 1);
482 
483  lagrange_to_bernstein_face.resize(fe_collection.size());
484 
485  for (unsigned int i = 0; i < fe_collection.size(); i++)
486  {
487  const FiniteElement<dim> &element = fe_collection[i];
488  const FE_Q_Base<dim> *fe_q =
489  dynamic_cast<const FE_Q_Base<dim> *>(&element);
490  Assert(fe_q != nullptr, ExcNotImplemented());
491 
492  const FE_Bernstein<dim> fe_bernstein(fe_q->get_degree());
493 
494  const unsigned int dofs_per_face = fe_q->dofs_per_face;
495  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; f++)
496  {
497  FullMatrix<double> face_interpolation_matrix(dofs_per_face,
498  dofs_per_face);
499 
500  fe_bernstein.get_face_interpolation_matrix(
501  *fe_q, face_interpolation_matrix, f);
502  lagrange_to_bernstein_face[i][f].reinit(dofs_per_face);
503  lagrange_to_bernstein_face[i][f] = face_interpolation_matrix;
504  lagrange_to_bernstein_face[i][f].compute_lu_factorization();
505  }
506  }
507  }
508 
509 } // namespace NonMatching
510 
511 #include "mesh_classifier.inst"
512 
const hp::FECollection< dim, spacedim > & get_fe_collection() const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition: fe_bernstein.cc:96
unsigned int get_degree() const
const unsigned int dofs_per_face
Definition: fe_data.h:423
unsigned int n_components() const
const unsigned int n_components
Definition: function.h:164
LocationToLevelSet location_to_level_set(const typename Triangulation< dim >::cell_iterator &cell) const
LocationToLevelSet determine_face_location_to_levelset(const typename Triangulation< dim >::active_cell_iterator &cell, const unsigned int face_index)
MeshClassifier(const DoFHandler< dim > &level_set_dof_handler, const VectorType &level_set)
AnalyticLevelSetDescription(const Function< dim > &level_set, const FiniteElement< dim > &element)
unsigned int active_fe_index(const typename Triangulation< dim >::active_cell_iterator &cell) const override
void get_local_level_set_values(const typename Triangulation< dim >::active_cell_iterator &cell, const unsigned int face_index, Vector< double > &local_levelset_values) override
void get_local_level_set_values(const typename Triangulation< dim >::active_cell_iterator &cell, const unsigned int face_index, Vector< double > &local_levelset_values) override
unsigned int active_fe_index(const typename Triangulation< dim >::active_cell_iterator &cell) const override
DiscreteLevelSetDescription(const DoFHandler< dim > &dof_handler, const VectorType &level_set)
Definition: vector.h:110
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
unsigned int size() const
Definition: collection.h:265
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcNeedsLAPACK()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:495
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
@ update_quadrature_points
Transformed quadrature points.
LocationToLevelSet location_from_dof_signs(const VectorType &local_levelset_values)
unsigned short int fe_index
Definition: types.h:60
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation