Reference documentation for deal.II version GIT relicensing-767-g85e8698264 2024-05-26 14:00: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\}}\)
Loading...
Searching...
No Matches
mesh_classifier.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
18
21#include <deal.II/fe/fe_q.h>
23
33#include <deal.II/lac/vector.h>
35
37
38#include <algorithm>
39
41
42namespace NonMatching
43{
44 namespace internal
45 {
46 namespace MeshClassifierImplementation
47 {
50 "The Triangulation has not been classified. You need to call the "
51 "reclassify()-function before using this function.");
52
55 "The incoming cell does not belong to the triangulation passed to "
56 "the constructor.");
57
63 template <typename VectorType>
65 location_from_dof_signs(const VectorType &local_levelset_values)
66 {
67 const auto min_max_element =
68 std::minmax_element(local_levelset_values.begin(),
69 local_levelset_values.end());
70
71 if (*min_max_element.second < 0)
73 if (0 < *min_max_element.first)
75
77 }
78
79
80
85 template <int dim, typename VectorType>
87 {
88 public:
93 const VectorType &level_set);
94
99 get_fe_collection() const override;
100
105 unsigned int
107 &cell) const override;
108
113 void
115 const typename Triangulation<dim>::active_cell_iterator &cell,
116 const unsigned int face_index,
117 Vector<double> &local_levelset_values) override;
118
119 private:
124
130 };
131
132
133
134 template <int dim, typename VectorType>
136 const DoFHandler<dim> &dof_handler,
137 const VectorType &level_set)
138 : dof_handler(&dof_handler)
139 , level_set(&level_set)
140 {}
141
142
143
144 template <int dim, typename VectorType>
147 {
148 return dof_handler->get_fe_collection();
149 }
150
151
152
153 template <int dim, typename VectorType>
154 void
156 const typename Triangulation<dim>::active_cell_iterator &cell,
157 const unsigned int face_index,
158 Vector<double> &local_levelset_values)
159 {
160 const auto cell_with_dofs = cell->as_dof_handler_iterator(*dof_handler);
161
162 const unsigned int n_dofs_per_face =
163 dof_handler->get_fe().n_dofs_per_face();
164 std::vector<types::global_dof_index> dof_indices(n_dofs_per_face);
165 cell_with_dofs->face(face_index)->get_dof_indices(dof_indices);
166
167 local_levelset_values.reinit(dof_indices.size());
168
169 for (unsigned int i = 0; i < dof_indices.size(); i++)
170 local_levelset_values[i] =
172 dof_indices[i]);
173 }
174
175
176
177 template <int dim, typename VectorType>
178 unsigned int
180 const typename Triangulation<dim>::active_cell_iterator &cell) const
181 {
182 const auto cell_with_dofs = cell->as_dof_handler_iterator(*dof_handler);
183
184 return cell_with_dofs->active_fe_index();
185 }
186
187
192 template <int dim>
194 {
195 public:
201 const FiniteElement<dim> &element);
202
208 get_fe_collection() const override;
209
214 unsigned int
216 &cell) const override;
217
222 void
224 const typename Triangulation<dim>::active_cell_iterator &cell,
225 const unsigned int face_index,
226 Vector<double> &local_levelset_values) override;
227
228 private:
233
239
245 };
246
247
248
249 template <int dim>
251 const Function<dim> &level_set,
252 const FiniteElement<dim> &element)
253 : level_set(&level_set)
254 , fe_collection(element)
255 , fe_face_values(element,
256 Quadrature<dim - 1>(
257 element.get_unit_face_support_points()),
259 {}
260
261
262
263 template <int dim>
264 void
266 const typename Triangulation<dim>::active_cell_iterator &cell,
267 const unsigned int face_index,
268 Vector<double> &local_levelset_values)
269 {
270 AssertDimension(local_levelset_values.size(),
271 fe_face_values.n_quadrature_points);
272
273 fe_face_values.reinit(cell, face_index);
274 const std::vector<Point<dim>> &points =
275 fe_face_values.get_quadrature_points();
276
277 for (unsigned int i = 0; i < points.size(); i++)
278 local_levelset_values[i] = level_set->value(points[i]);
279 }
280
281
282
283 template <int dim>
286 {
287 return fe_collection;
288 }
289
290
291
292 template <int dim>
293 unsigned int
299 } // namespace MeshClassifierImplementation
300 } // namespace internal
301
302
303
304 template <int dim>
305 template <typename VectorType>
307 const VectorType &level_set)
308 : triangulation(&dof_handler.get_triangulation())
309 , level_set_description(
310 std::make_unique<internal::MeshClassifierImplementation::
311 DiscreteLevelSetDescription<dim, VectorType>>(
312 dof_handler,
313 level_set))
314 {
315#ifdef DEAL_II_WITH_LAPACK
316 const hp::FECollection<dim> &fe_collection =
317 dof_handler.get_fe_collection();
318 for (unsigned int i = 0; i < fe_collection.size(); i++)
319 {
320 // The level set function must be scalar.
321 AssertDimension(fe_collection[i].n_components(), 1);
322
323 Assert(fe_collection[i].has_face_support_points(),
325 "The elements in the FECollection of the incoming DoFHandler "
326 "must have face support points."));
327 }
328#else
329 AssertThrow(false, ExcNeedsLAPACK());
330#endif
331 }
332
333
334
335 template <int dim>
337 const Function<dim> &level_set,
338 const FiniteElement<dim> &element)
340 , level_set_description(
341 std::make_unique<internal::MeshClassifierImplementation::
342 AnalyticLevelSetDescription<dim>>(level_set,
343 element))
344 {
345 // The level set function must be scalar.
346 AssertDimension(level_set.n_components, 1);
347 AssertDimension(element.n_components(), 1);
348 }
349
350
351
352 template <int dim>
353 void
355 {
356 initialize();
357 cell_locations.assign(triangulation->n_active_cells(),
359 face_locations.assign(triangulation->n_raw_faces(),
361
362 // Loop over all cells and determine the location of all non artificial
363 // cells and faces.
364 for (const auto &cell : triangulation->active_cell_iterators())
365 if (!cell->is_artificial())
366 {
367 const LocationToLevelSet face0_location =
368 determine_face_location_to_levelset(cell, 0);
369
370 face_locations[cell->face(0)->index()] = face0_location;
371 LocationToLevelSet cell_location = face0_location;
372
373 for (unsigned int f = 1; f < GeometryInfo<dim>::faces_per_cell; ++f)
374 {
375 const LocationToLevelSet face_location =
376 determine_face_location_to_levelset(cell, f);
377
378 face_locations[cell->face(f)->index()] = face_location;
379
380 if (face_location != face0_location)
381 cell_location = LocationToLevelSet::intersected;
382 }
383 cell_locations[cell->active_cell_index()] = cell_location;
384 }
385 }
386
387
388
389 template <int dim>
392 const typename Triangulation<dim>::active_cell_iterator &cell,
393 const unsigned int face_index)
394 {
395 // The location of the face might already be computed on the neighboring
396 // cell. If this is the case we just return the value.
397 const LocationToLevelSet location =
398 face_locations.at(cell->face(face_index)->index());
399 if (location != LocationToLevelSet::unassigned)
400 return location;
401
402 // Determine the location by changing basis to FE_Bernstein and checking
403 // the signs of the dofs.
404 const unsigned int fe_index = level_set_description->active_fe_index(cell);
405 const unsigned int n_local_dofs =
406 lagrange_to_bernstein_face[fe_index][face_index].m();
407
408 Vector<double> local_levelset_values(n_local_dofs);
409 level_set_description->get_local_level_set_values(cell,
410 face_index,
411 local_levelset_values);
412
413 const FiniteElement<dim> &fe =
414 level_set_description->get_fe_collection()[fe_index];
415
416 const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
417 dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe);
418
419 const FE_Poly<dim> *fe_poly = dynamic_cast<const FE_Poly<dim> *>(&fe);
420
421 const bool is_linear = fe_q_iso_q1 != nullptr ||
422 (fe_poly != nullptr && fe_poly->get_degree() == 1);
423
424 // shortcut for linear elements
425 if (is_linear)
426 {
428 local_levelset_values);
429 }
430
431 lagrange_to_bernstein_face[fe_index][face_index].solve(
432 local_levelset_values);
433
435 local_levelset_values);
436 }
437
438
439
440 template <int dim>
443 const typename Triangulation<dim>::cell_iterator &cell) const
444 {
445 Assert(cell_locations.size() == triangulation->n_active_cells(),
447 Assert(&cell->get_triangulation() == triangulation,
449
450 return cell_locations.at(cell->active_cell_index());
451 }
452
453
454
455 template <int dim>
458 const typename Triangulation<dim>::cell_iterator &cell,
459 const unsigned int face_index) const
460 {
462 Assert(face_locations.size() == triangulation->n_raw_faces(),
464 Assert(&cell->get_triangulation() == triangulation,
466
467 return face_locations.at(cell->face(face_index)->index());
468 }
469
470
471
472 template <int dim>
473 void
475 {
476 const hp::FECollection<dim> &fe_collection =
477 level_set_description->get_fe_collection();
478
479 // The level set function must be scalar.
480 AssertDimension(fe_collection.n_components(), 1);
481
482 lagrange_to_bernstein_face.resize(fe_collection.size());
483
484 for (unsigned int i = 0; i < fe_collection.size(); i++)
485 {
486 const FiniteElement<dim> &element = fe_collection[i];
487 const FE_Q_Base<dim> *fe_q =
488 dynamic_cast<const FE_Q_Base<dim> *>(&element);
489 Assert(fe_q != nullptr, ExcNotImplemented());
490
491 const FE_Bernstein<dim> fe_bernstein(fe_q->get_degree());
492
493 const unsigned int dofs_per_face = fe_q->dofs_per_face;
494 for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; f++)
495 {
496 FullMatrix<double> face_interpolation_matrix(dofs_per_face,
497 dofs_per_face);
498
500 *fe_q, face_interpolation_matrix, f);
501 lagrange_to_bernstein_face[i][f].reinit(dofs_per_face);
502 lagrange_to_bernstein_face[i][f] = face_interpolation_matrix;
503 lagrange_to_bernstein_face[i][f].compute_lu_factorization();
504 }
505 }
506 }
507
508} // namespace NonMatching
509
510#include "mesh_classifier.inst"
511
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
unsigned int get_degree() const
const unsigned int dofs_per_face
Definition fe_data.h:422
unsigned int n_components() const
const unsigned int n_components
Definition function.h:163
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)
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:264
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcNeedsLAPACK()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_quadrature_points
Transformed quadrature points.
LocationToLevelSet location_from_dof_signs(const VectorType &local_levelset_values)
STL namespace.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation