Reference documentation for deal.II version GIT 3e82abc508 2023-06-09 03:50:01+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 
21 #include <deal.II/fe/fe_q.h>
22 #include <deal.II/fe/fe_values.h>
23 
27 #include <deal.II/lac/la_vector.h>
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 <class 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, class 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, class 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, class VectorType>
146  const hp::FECollection<dim> &
148  {
149  return dof_handler->get_fe_collection();
150  }
151 
152 
153 
154  template <int dim, class 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  typename DoFHandler<dim>::active_cell_iterator cell_with_dofs(
162  &dof_handler->get_triangulation(),
163  cell->level(),
164  cell->index(),
165  dof_handler);
166 
167  const unsigned int n_dofs_per_face =
168  dof_handler->get_fe().n_dofs_per_face();
169  std::vector<types::global_dof_index> dof_indices(n_dofs_per_face);
170  cell_with_dofs->face(face_index)->get_dof_indices(dof_indices);
171 
172  local_levelset_values.reinit(dof_indices.size());
173 
174  for (unsigned int i = 0; i < dof_indices.size(); i++)
175  local_levelset_values[i] =
177  dof_indices[i]);
178  }
179 
180 
181 
182  template <int dim, class VectorType>
183  unsigned int
185  const typename Triangulation<dim>::active_cell_iterator &cell) const
186  {
187  typename DoFHandler<dim>::active_cell_iterator cell_with_dofs(
188  &dof_handler->get_triangulation(),
189  cell->level(),
190  cell->index(),
191  dof_handler);
192 
193  return cell_with_dofs->active_fe_index();
194  }
195 
196 
201  template <int dim>
203  {
204  public:
210  const FiniteElement<dim> &element);
211 
216  const hp::FECollection<dim> &
217  get_fe_collection() const override;
218 
223  unsigned int
225  &cell) const override;
226 
231  void
233  const typename Triangulation<dim>::active_cell_iterator &cell,
234  const unsigned int face_index,
235  Vector<double> &local_levelset_values) override;
236 
237  private:
242 
248 
254  };
255 
256 
257 
258  template <int dim>
260  const Function<dim> & level_set,
261  const FiniteElement<dim> &element)
262  : level_set(&level_set)
263  , fe_collection(element)
264  , fe_face_values(element,
265  Quadrature<dim - 1>(
266  element.get_unit_face_support_points()),
268  {}
269 
270 
271 
272  template <int dim>
273  void
275  const typename Triangulation<dim>::active_cell_iterator &cell,
276  const unsigned int face_index,
277  Vector<double> &local_levelset_values)
278  {
279  AssertDimension(local_levelset_values.size(),
280  fe_face_values.n_quadrature_points);
281 
282  fe_face_values.reinit(cell, face_index);
283  const std::vector<Point<dim>> &points =
284  fe_face_values.get_quadrature_points();
285 
286  for (unsigned int i = 0; i < points.size(); i++)
287  local_levelset_values[i] = level_set->value(points[i]);
288  }
289 
290 
291 
292  template <int dim>
293  const hp::FECollection<dim> &
295  {
296  return fe_collection;
297  }
298 
299 
300 
301  template <int dim>
302  unsigned int
304  const typename Triangulation<dim>::active_cell_iterator &) const
305  {
306  return 0;
307  }
308  } // namespace MeshClassifierImplementation
309  } // namespace internal
310 
311 
312 
313  template <int dim>
314  template <class VectorType>
316  const VectorType & level_set)
317  : triangulation(&dof_handler.get_triangulation())
318  , level_set_description(
319  std::make_unique<internal::MeshClassifierImplementation::
320  DiscreteLevelSetDescription<dim, VectorType>>(
321  dof_handler,
322  level_set))
323  {
324 #ifdef DEAL_II_WITH_LAPACK
325  const hp::FECollection<dim> &fe_collection =
326  dof_handler.get_fe_collection();
327  for (unsigned int i = 0; i < fe_collection.size(); i++)
328  {
329  // The level set function must be scalar.
330  AssertDimension(fe_collection[i].n_components(), 1);
331 
332  Assert(fe_collection[i].has_face_support_points(),
333  ExcMessage(
334  "The elements in the FECollection of the incoming DoFHandler "
335  "must have face support points."));
336  }
337 #else
338  AssertThrow(false, ExcNeedsLAPACK());
339 #endif
340  }
341 
342 
343 
344  template <int dim>
346  const Function<dim> & level_set,
347  const FiniteElement<dim> &element)
349  , level_set_description(
350  std::make_unique<internal::MeshClassifierImplementation::
351  AnalyticLevelSetDescription<dim>>(level_set,
352  element))
353  {
354  // The level set function must be scalar.
355  AssertDimension(level_set.n_components, 1);
356  AssertDimension(element.n_components(), 1);
357  }
358 
359 
360 
361  template <int dim>
362  void
364  {
365  initialize();
366  cell_locations.assign(triangulation->n_active_cells(),
368  face_locations.assign(triangulation->n_raw_faces(),
370 
371  // Loop over all cells and determine the location of all non artificial
372  // cells and faces.
373  for (const auto &cell : triangulation->active_cell_iterators())
374  if (!cell->is_artificial())
375  {
376  const LocationToLevelSet face0_location =
377  determine_face_location_to_levelset(cell, 0);
378 
379  face_locations[cell->face(0)->index()] = face0_location;
380  LocationToLevelSet cell_location = face0_location;
381 
382  for (unsigned int f = 1; f < GeometryInfo<dim>::faces_per_cell; ++f)
383  {
384  const LocationToLevelSet face_location =
385  determine_face_location_to_levelset(cell, f);
386 
387  face_locations[cell->face(f)->index()] = face_location;
388 
389  if (face_location != face0_location)
390  cell_location = LocationToLevelSet::intersected;
391  }
392  cell_locations[cell->active_cell_index()] = cell_location;
393  }
394  }
395 
396 
397 
398  template <int dim>
401  const typename Triangulation<dim>::active_cell_iterator &cell,
402  const unsigned int face_index)
403  {
404  // The location of the face might already be computed on the neighboring
405  // cell. If this is the case we just return the value.
406  const LocationToLevelSet location =
407  face_locations.at(cell->face(face_index)->index());
408  if (location != LocationToLevelSet::unassigned)
409  return location;
410 
411  // Determine the location by changing basis to FE_Bernstein and checking
412  // the signs of the dofs.
413  const unsigned int fe_index = level_set_description->active_fe_index(cell);
414  const unsigned int n_local_dofs =
415  lagrange_to_bernstein_face[fe_index][face_index].m();
416 
417  Vector<double> local_levelset_values(n_local_dofs);
418  level_set_description->get_local_level_set_values(cell,
419  face_index,
420  local_levelset_values);
421 
422  lagrange_to_bernstein_face[fe_index][face_index].solve(
423  local_levelset_values);
424 
426  local_levelset_values);
427  }
428 
429 
430 
431  template <int dim>
434  const typename Triangulation<dim>::cell_iterator &cell) const
435  {
436  Assert(cell_locations.size() == triangulation->n_active_cells(),
438  Assert(&cell->get_triangulation() == triangulation,
440 
441  return cell_locations.at(cell->active_cell_index());
442  }
443 
444 
445 
446  template <int dim>
449  const typename Triangulation<dim>::cell_iterator &cell,
450  const unsigned int face_index) const
451  {
453  Assert(face_locations.size() == triangulation->n_raw_faces(),
455  Assert(&cell->get_triangulation() == triangulation,
457 
458  return face_locations.at(cell->face(face_index)->index());
459  }
460 
461 
462 
463  template <int dim>
464  void
466  {
467  const hp::FECollection<dim> &fe_collection =
468  level_set_description->get_fe_collection();
469 
470  // The level set function must be scalar.
471  AssertDimension(fe_collection.n_components(), 1);
472 
473  lagrange_to_bernstein_face.resize(fe_collection.size());
474 
475  for (unsigned int i = 0; i < fe_collection.size(); i++)
476  {
477  const FiniteElement<dim> &element = fe_collection[i];
478  const FE_Q<dim> *fe_q = dynamic_cast<const FE_Q<dim> *>(&element);
479  Assert(fe_q != nullptr, ExcNotImplemented());
480 
481  const FE_Bernstein<dim> fe_bernstein(fe_q->get_degree());
482 
483  const unsigned int dofs_per_face = fe_q->dofs_per_face;
484  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; f++)
485  {
486  FullMatrix<double> face_interpolation_matrix(dofs_per_face,
487  dofs_per_face);
488 
489  fe_bernstein.get_face_interpolation_matrix(
490  *fe_q, face_interpolation_matrix, f);
491  lagrange_to_bernstein_face[i][f].reinit(dofs_per_face);
492  lagrange_to_bernstein_face[i][f] = face_interpolation_matrix;
493  lagrange_to_bernstein_face[i][f].compute_lu_factorization();
494  }
495  }
496  }
497 
498 } // namespace NonMatching
499 
500 #include "mesh_classifier.inst"
501 
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
Definition: fe_q.h:551
const unsigned int dofs_per_face
Definition: fe_data.h:419
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:109
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
unsigned int size() const
Definition: collection.h:264
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:475
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:476
static ::ExceptionBase & ExcNeedsLAPACK()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:488
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:441
@ 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