Reference documentation for deal.II version 9.5.0
\(\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// 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>
23
34#include <deal.II/lac/vector.h>
36
38
39#include <algorithm>
40
42
43namespace 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
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>
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
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>
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(),
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
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
unsigned int get_degree() const
Definition fe_q.h:551
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)
size_type size() const
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ExcNeedsLAPACK()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ 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