Reference documentation for deal.II version 9.3.3
fe_raviart_thomas.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2003 - 2021 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
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
16
20#include <deal.II/base/table.h>
22
24
25#include <deal.II/fe/fe.h>
27#include <deal.II/fe/fe_tools.h>
29#include <deal.II/fe/mapping.h>
30
31#include <deal.II/grid/tria.h>
33
34#include <iostream>
35#include <memory>
36#include <sstream>
37
38
40
41
42template <int dim>
44 : FE_PolyTensor<dim>(
47 dim,
48 deg + 1,
49 FiniteElementData<dim>::Hdiv),
50 std::vector<bool>(PolynomialsRaviartThomas<dim>::n_polynomials(deg),
51 true),
53 deg),
54 std::vector<bool>(dim, true)))
55{
56 Assert(dim >= 2, ExcImpossibleInDim(dim));
57 const unsigned int n_dofs = this->n_dofs_per_cell();
58
60 // First, initialize the
61 // generalized support points and
62 // quadrature weights, since they
63 // are required for interpolation.
65
66 // Now compute the inverse node matrix, generating the correct
67 // basis functions from the raw ones. For a discussion of what
68 // exactly happens here, see FETools::compute_node_matrix.
70 this->inverse_node_matrix.reinit(n_dofs, n_dofs);
72 // From now on, the shape functions provided by FiniteElement::shape_value
73 // and similar functions will be the correct ones, not
74 // the raw shape functions from the polynomial space anymore.
75
76 // Reinit the vectors of
77 // restriction and prolongation
78 // matrices to the right sizes.
79 // Restriction only for isotropic
80 // refinement
82 // Fill prolongation matrices with embedding operators
85
86 // TODO: the implementation makes the assumption that all faces have the
87 // same number of dofs
89 const unsigned int face_no = 0;
90
91 // TODO[TL]: for anisotropic refinement we will probably need a table of
92 // submatrices with an array for each refine case
94 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
95 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
96 this->n_dofs_per_face(face_no));
97 FETools::compute_face_embedding_matrices<dim, double>(*this,
98 face_embeddings,
99 0,
100 0);
101 this->interface_constraints.reinit((1 << (dim - 1)) *
102 this->n_dofs_per_face(face_no),
103 this->n_dofs_per_face(face_no));
104 unsigned int target_row = 0;
105 for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
106 for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
107 {
108 for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
109 this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
110 ++target_row;
111 }
112
113 // We need to initialize the dof permutation table and the one for the sign
114 // change.
116}
117
118
119
120template <int dim>
121std::string
123{
124 // note that the
125 // FETools::get_fe_by_name
126 // function depends on the
127 // particular format of the string
128 // this function returns, so they
129 // have to be kept in synch
130
131 // note that this->degree is the maximal
132 // polynomial degree and is thus one higher
133 // than the argument given to the
134 // constructor
135 std::ostringstream namebuf;
136 namebuf << "FE_RaviartThomas<" << dim << ">(" << this->degree - 1 << ")";
137
138 return namebuf.str();
139}
140
141
142template <int dim>
143std::unique_ptr<FiniteElement<dim, dim>>
145{
146 return std::make_unique<FE_RaviartThomas<dim>>(*this);
147}
148
149
150//---------------------------------------------------------------------------
151// Auxiliary and internal functions
152//---------------------------------------------------------------------------
153
154
155template <int dim>
156void
158{
160 const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.size() : 0;
161
162 // TODO: the implementation makes the assumption that all faces have the
163 // same number of dofs
164 AssertDimension(this->n_unique_faces(), 1);
165 const unsigned int face_no = 0;
166
167 unsigned int n_face_points = (dim > 1) ? 1 : 0;
168 // compute (deg+1)^(dim-1)
169 for (unsigned int d = 1; d < dim; ++d)
170 n_face_points *= deg + 1;
171
172
173 this->generalized_support_points.resize(
174 GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
175 this->generalized_face_support_points[face_no].resize(n_face_points);
176
177 // Number of the point being entered
178 unsigned int current = 0;
179
180 if (dim > 1)
181 {
182 QGauss<dim - 1> face_points(deg + 1);
185
186 boundary_weights.reinit(n_face_points, legendre.n());
187
188 for (unsigned int k = 0; k < n_face_points; ++k)
189 {
190 this->generalized_face_support_points[face_no][k] =
191 face_points.point(k);
193 // contribution for each
194 // moment.
195 for (unsigned int i = 0; i < legendre.n(); ++i)
196 {
197 boundary_weights(k, i) =
198 face_points.weight(k) *
199 legendre.compute_value(i, face_points.point(k));
200 }
201 }
202
205 face_points);
206 for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
207 ++current)
208 {
209 // Enter the support point
210 // into the vector
211 this->generalized_support_points[current] = faces.point(
212 current +
214 this->reference_cell(), 0, true, false, false, n_face_points));
215 }
216 }
217
218 if (deg == 0)
219 return;
220
221 // Create Legendre basis for the space D_xi Q_k
222 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
223 for (unsigned int dd = 0; dd < dim; ++dd)
224 {
225 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
226 for (unsigned int d = 0; d < dim; ++d)
229
230 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
231 }
232
233 interior_weights.reinit(
234 TableIndices<3>(n_interior_points, polynomials[0]->n(), dim));
235
236 for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
237 {
239 for (unsigned int i = 0; i < polynomials[0]->n(); ++i)
240 for (unsigned int d = 0; d < dim; ++d)
241 interior_weights(k, i, d) =
244 }
245
246 Assert(current == this->generalized_support_points.size(),
248}
249
250
251template <int dim>
252void
254{
255 // For 1D do nothing.
256 //
257 // TODO: For 2D we simply keep the legacy behavior for now. This should be
258 // changed in the future and can be taken care of by similar means as the 3D
259 // case below. The legacy behavior can be found in fe_poly_tensor.cc in the
260 // function internal::FE_PolyTensor::get_dof_sign_change_h_div(...)
261 if (dim < 3)
262 return;
263
264 // TODO: the implementation makes the assumption that all faces have the
265 // same number of dofs
266 AssertDimension(this->n_unique_faces(), 1);
267 const unsigned int face_no = 0;
268
270 .n_elements() == 8 * this->n_dofs_per_quad(face_no),
272
273 // The 3D RaviartThomas space has tensor_degree*tensor_degree face dofs
274 const unsigned int n = this->tensor_degree();
275 Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
276
278 // contains offsets for local face_dofs and is being filled with zeros in
279 // fe.cc. We need to fill it with the correct values in case of non-standard,
280 // flipped (rotated by +180 degrees) or rotated (rotated by +90 degrees) faces
281
282 // The dofs on a face are connected to a n x n matrix. for example, for
283 // tensor_degree==3 we have the following dofs on a quad:
284
285 // ___________
286 // | |
287 // | 6 7 8 |
288 // | |
289 // | 3 4 5 |
290 // | |
291 // | 0 1 2 |
292 // |___________|
293 //
294 // We have dof_index=i+n*j with index i in x-direction and index j in
295 // y-direction running from 0 to n-1. to extract i and j we can use
296 // i=dof_index%n and j=dof_index/n. The indices i and j can then be used to
297 // compute the offset.
298
299 // Example: if the switches are (true | true | true) that means we rotate the
300 // face first by + 90 degree(counterclockwise) then by another +180
301 // degrees but we do not flip it since the face has standard
302 // orientation. The flip axis is the diagonal from the lower left to the upper
303 // right corner of the face. With these flags the configuration above becomes:
304 // ___________
305 // | |
306 // | 2 5 8 |
307 // | |
308 // | 1 4 7 |
309 // | |
310 // | 0 3 6 |
311 // |___________|
312 //
313 // This is exactly what is realized by the formulas implemented below. Note
314 // that the necessity of a permuattion depends on the three flags.
315
316 // There is also a pattern for the sign change of the mapped face_dofs
317 // depending on the switches. In the above example it would be
318 // ___________
319 // | |
320 // | + - + |
321 // | |
322 // | + - + |
323 // | |
324 // | + - + |
325 // |___________|
326 //
327
328 for (unsigned int local_face_dof = 0;
330 ++local_face_dof)
331 {
332 // Row and column
333 unsigned int i = local_face_dof % n;
334 unsigned int j = local_face_dof / n;
335
336 // We have 8 cases that are all treated the same way. Note that the
337 // corresponding case to case_no is just its binary representation.
338 // The above example of (false | true | true) would be case_no=3
339 for (unsigned int case_no = 0; case_no < 8; ++case_no)
340 {
341 // Get the binary representation of the case
342 const bool face_orientation = Utilities::get_bit(case_no, 2);
343 const bool face_flip = Utilities::get_bit(case_no, 1);
344 const bool face_rotation = Utilities::get_bit(case_no, 0);
345
346 if (((!face_orientation) && (!face_rotation)) ||
347 ((face_orientation) && (face_rotation)))
348 {
349 // We flip across the diagonal
350 // This is the local face dof offset
352 local_face_dof, case_no) = j + i * n - local_face_dof;
353 }
354 else
355 {
356 // Offset is zero
358 local_face_dof, case_no) = 0;
359 } // if face needs dof permutation
360
361 // Get new local face_dof by adding offset
362 const unsigned int new_local_face_dof =
363 local_face_dof +
365 local_face_dof, case_no);
366 // compute new row and column index
367 i = new_local_face_dof % n;
368 j = new_local_face_dof / n;
369
370 /*
371 * Now compute if a sign change is necessary. This is done for the
372 * case of face_orientation==true
373 */
374 // flip = false, rotation=true
375 if (!face_flip && face_rotation)
376 {
378 local_face_dof, case_no) = ((j % 2) == 1);
379 }
380 // flip = true, rotation=false
381 else if (face_flip && !face_rotation)
382 {
383 // This case is symmetric (although row and column may be
384 // switched)
386 local_face_dof, case_no) = ((j % 2) == 1) != ((i % 2) == 1);
387 }
388 // flip = true, rotation=true
389 else if (face_flip && face_rotation)
390 {
392 local_face_dof, case_no) = ((i % 2) == 1);
393 }
394 /*
395 * flip = false, rotation=false => nothing to do
396 */
397
398 /*
399 * If face_orientation==false the sign flip is exactly the opposite.
400 */
401 if (!face_orientation)
403 local_face_dof, case_no) =
405 local_face_dof, case_no);
406 } // case_no
407 } // local_face_dof
408}
409
410
411template <>
412void
414{
415 // there is only one refinement case in 1d,
416 // which is the isotropic one (first index of
417 // the matrix array has to be 0)
418 for (unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
419 this->restriction[0][i].reinit(0, 0);
420}
421
422
423
424// This function is the same Raviart-Thomas interpolation performed by
425// interpolate. Still, we cannot use interpolate, since it was written
426// for smooth functions. The functions interpolated here are not
427// smooth, maybe even not continuous. Therefore, we must double the
428// number of quadrature points in each direction in order to integrate
429// only smooth functions.
430
431// Then again, the interpolated function is chosen such that the
432// moments coincide with the function to be interpolated.
433
434template <int dim>
435void
437{
438 const unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
439
440 QGauss<dim - 1> q_base(this->degree);
441 const unsigned int n_face_points = q_base.size();
442 // First, compute interpolation on
443 // subfaces
444 for (unsigned int face : GeometryInfo<dim>::face_indices())
445 {
446 // The shape functions of the
447 // child cell are evaluated
448 // in the quadrature points
449 // of a full face.
452 // Store shape values, since the
453 // evaluation suffers if not
454 // ordered by point
455 Table<2, double> cached_values_on_face(this->n_dofs_per_cell(),
456 q_face.size());
457 for (unsigned int k = 0; k < q_face.size(); ++k)
458 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
459 cached_values_on_face(i, k) = this->shape_value_component(
461
462 for (unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
463 ++sub)
464 {
465 // The weight functions for
466 // the coarse face are
467 // evaluated on the subface
468 // only.
470 this->reference_cell(), q_base, face, sub);
471 const unsigned int child = GeometryInfo<dim>::child_cell_on_face(
473
474 // On a certain face, we must
475 // compute the moments of ALL
476 // fine level functions with
477 // the coarse level weight
478 // functions belonging to
479 // that face. Due to the
480 // orthogonalization process
481 // when building the shape
482 // functions, these weights
483 // are equal to the
484 // corresponding shape
485 // functions.
486 for (unsigned int k = 0; k < n_face_points; ++k)
487 for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
488 ++i_child)
489 for (unsigned int i_face = 0;
490 i_face < this->n_dofs_per_face(face);
491 ++i_face)
492 {
494 // weights on the
495 // subcell are NOT
496 // transformed, so we
497 // have to do it here.
498 this->restriction[iso][child](
499 face * this->n_dofs_per_face(face) + i_face, i_child) +=
500 Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
501 cached_values_on_face(i_child, k) *
502 this->shape_value_component(
503 face * this->n_dofs_per_face(face) + i_face,
504 q_sub.point(k),
506 }
507 }
508 }
509
510 if (this->degree == 1)
511 return;
512
513 // Create Legendre basis for the space D_xi Q_k. Here, we cannot
514 // use the shape functions
515 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
516 for (unsigned int dd = 0; dd < dim; ++dd)
517 {
518 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
519 for (unsigned int d = 0; d < dim; ++d)
520 poly[d] =
522 poly[dd] =
524
525 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
526 }
527
528 // TODO: the implementation makes the assumption that all faces have the
529 // same number of dofs
530 AssertDimension(this->n_unique_faces(), 1);
531 const unsigned int face_no = 0;
532
533 QGauss<dim> q_cell(this->degree);
534 const unsigned int start_cell_dofs =
535 GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
536
537 // Store shape values, since the
538 // evaluation suffers if not
539 // ordered by point
540 Table<3, double> cached_values_on_cell(this->n_dofs_per_cell(),
541 q_cell.size(),
542 dim);
543 for (unsigned int k = 0; k < q_cell.size(); ++k)
544 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
545 for (unsigned int d = 0; d < dim; ++d)
546 cached_values_on_cell(i, k, d) =
547 this->shape_value_component(i, q_cell.point(k), d);
548
549 for (unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
550 ++child)
551 {
554 q_cell,
555 child);
556
557 for (unsigned int k = 0; k < q_sub.size(); ++k)
558 for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
559 ++i_child)
560 for (unsigned int d = 0; d < dim; ++d)
561 for (unsigned int i_weight = 0; i_weight < polynomials[d]->n();
562 ++i_weight)
563 {
564 this->restriction[iso][child](start_cell_dofs + i_weight * dim +
565 d,
566 i_child) +=
567 q_sub.weight(k) * cached_values_on_cell(i_child, k, d) *
568 polynomials[d]->compute_value(i_weight, q_sub.point(k));
569 }
570 }
571}
572
573
574
575template <int dim>
576std::vector<unsigned int>
578{
579 // the element is face-based and we have
580 // (deg+1)^(dim-1) DoFs per face
581 unsigned int dofs_per_face = 1;
582 for (unsigned int d = 1; d < dim; ++d)
583 dofs_per_face *= deg + 1;
584
585 // and then there are interior dofs
586 const unsigned int interior_dofs = dim * deg * dofs_per_face;
587
588 std::vector<unsigned int> dpo(dim + 1);
589 dpo[dim - 1] = dofs_per_face;
590 dpo[dim] = interior_dofs;
591
592 return dpo;
593}
594
595
596
597template <int dim>
598std::pair<Table<2, bool>, std::vector<unsigned int>>
600{
601 Table<2, bool> constant_modes(dim, this->n_dofs_per_cell());
602 for (unsigned int d = 0; d < dim; ++d)
603 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
604 constant_modes(d, i) = true;
605 std::vector<unsigned int> components;
606 for (unsigned int d = 0; d < dim; ++d)
607 components.push_back(d);
608 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
609 components);
610}
611
612
613
614//---------------------------------------------------------------------------
615// Data field initialization
616//---------------------------------------------------------------------------
617
618
619template <int dim>
620bool
621FE_RaviartThomas<dim>::has_support_on_face(const unsigned int shape_index,
622 const unsigned int face_index) const
623{
624 AssertIndexRange(shape_index, this->n_dofs_per_cell());
626
627 // Return computed values if we
628 // know them easily. Otherwise, it
629 // is always safe to return true.
630 switch (this->degree)
631 {
632 case 1:
633 {
634 switch (dim)
635 {
636 case 2:
637 {
638 // only on the one
640 // are the values
641 // actually zero. list
642 // these in a table
643 return (face_index !=
645 }
646
647 default:
648 return true;
649 }
650 }
651
652 default: // other rt_order
653 return true;
654 }
655
656 return true;
657}
658
659
660
661template <int dim>
662void
664 const std::vector<Vector<double>> &support_point_values,
665 std::vector<double> & nodal_values) const
666{
667 Assert(support_point_values.size() == this->generalized_support_points.size(),
668 ExcDimensionMismatch(support_point_values.size(),
669 this->generalized_support_points.size()));
670 Assert(nodal_values.size() == this->n_dofs_per_cell(),
671 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
672 Assert(support_point_values[0].size() == this->n_components(),
673 ExcDimensionMismatch(support_point_values[0].size(),
674 this->n_components()));
675
676 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
677
678 const unsigned int n_face_points = boundary_weights.size(0);
679 for (unsigned int face : GeometryInfo<dim>::face_indices())
680 for (unsigned int k = 0; k < n_face_points; ++k)
681 for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
682 {
683 nodal_values[i + face * this->n_dofs_per_face(face)] +=
684 boundary_weights(k, i) *
685 support_point_values[face * n_face_points + k](
687 }
688
689 // TODO: the implementation makes the assumption that all faces have the
690 // same number of dofs
691 AssertDimension(this->n_unique_faces(), 1);
692 const unsigned int face_no = 0;
693
694 const unsigned int start_cell_dofs =
695 GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
696 const unsigned int start_cell_points =
697 GeometryInfo<dim>::faces_per_cell * n_face_points;
698
699 for (unsigned int k = 0; k < interior_weights.size(0); ++k)
700 for (unsigned int i = 0; i < interior_weights.size(1); ++i)
701 for (unsigned int d = 0; d < dim; ++d)
702 nodal_values[start_cell_dofs + i * dim + d] +=
703 interior_weights(k, i, d) *
704 support_point_values[k + start_cell_points](d);
705}
706
707
708
709template <int dim>
710std::size_t
712{
713 Assert(false, ExcNotImplemented());
714 return 0;
715}
716
717
718
719// explicit instantiations
720#include "fe_raviart_thomas.inst"
721
722
FullMatrix< double > inverse_node_matrix
std::vector< MappingKind > mapping_kind
virtual std::size_t memory_consumption() const override
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
friend class FE_RaviartThomas
void initialize_support_points(const unsigned int rt_degree)
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::string get_name() const override
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FullMatrix< double > interface_constraints
Definition: fe.h:2430
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2418
size_type n() const
void invert(const FullMatrix< number2 > &M)
size_type m() const
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:744
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
Definition: qprojector.h:579
Definition: qprojector.cc:1238
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
@ mapping_raviart_thomas
Definition: mapping.h:120
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
bool get_bit(const unsigned char number, const unsigned int n)
Definition: utilities.h:1376
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
double legendre(unsigned int l, double x)
Definition: cmath.h:75
STL namespace.
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)