1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2020 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
16
17 #include <deal.II/base/qprojector.h>
18 #include <deal.II/base/quadrature.h>
19 #include <deal.II/base/quadrature_lib.h>
20 #include <deal.II/base/table.h>
21 #include <deal.II/base/utilities.h>
22
23 #include <deal.II/dofs/dof_accessor.h>
24
25 #include <deal.II/fe/fe.h>
26 #include <deal.II/fe/fe_abf.h>
27 #include <deal.II/fe/fe_tools.h>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping.h>
30
31 #include <deal.II/grid/tria.h>
32 #include <deal.II/grid/tria_iterator.h>
33
34 #include <iostream>
35 #include <memory>
36 #include <sstream>
37
38
39 // TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
40 // adjust_line_dof_index_for_line_orientation_table fields, and write tests
41 // similar to bits/face_orientation_and_fe_q_*
42
43
44 DEAL_II_NAMESPACE_OPEN
45
46
47 template <int dim>
FE_ABF(const unsigned int deg)48 FE_ABF<dim>::FE_ABF(const unsigned int deg)
49 : FE_PolyTensor<dim>(
50 PolynomialsABF<dim>(deg),
51 FiniteElementData<dim>(get_dpo_vector(deg),
52 dim,
53 deg + 2,
54 FiniteElementData<dim>::Hdiv),
55 std::vector<bool>(PolynomialsABF<dim>::n_polynomials(deg), true),
56 std::vector<ComponentMask>(PolynomialsABF<dim>::n_polynomials(deg),
57 std::vector<bool>(dim, true)))
58 , rt_order(deg)
59 {
60 Assert(dim >= 2, ExcImpossibleInDim(dim));
61 const unsigned int n_dofs = this->n_dofs_per_cell();
62
63 this->mapping_kind = {mapping_raviart_thomas};
64 // First, initialize the
65 // generalized support points and
66 // quadrature weights, since they
67 // are required for interpolation.
68 initialize_support_points(deg);
69
70 // Now compute the inverse node matrix, generating the correct
71 // basis functions from the raw ones. For a discussion of what
72 // exactly happens here, see FETools::compute_node_matrix.
73 const FullMatrix<double> M = FETools::compute_node_matrix(*this);
74 this->inverse_node_matrix.reinit(n_dofs, n_dofs);
75 this->inverse_node_matrix.invert(M);
76 // From now on, the shape functions provided by FiniteElement::shape_value
77 // and similar functions will be the correct ones, not
78 // the raw shape functions from the polynomial space anymore.
79
80 // Reinit the vectors of
81 // restriction and prolongation
82 // matrices to the right sizes.
83 // Restriction only for isotropic
84 // refinement
85 this->reinit_restriction_and_prolongation_matrices(true);
86 // Fill prolongation matrices with embedding operators
87 FETools::compute_embedding_matrices(*this, this->prolongation, false, 1.e-10);
88
89 initialize_restriction();
90
91 // TODO: the implementation makes the assumption that all faces have the
92 // same number of dofs
93 AssertDimension(this->n_unique_faces(), 1);
94 const unsigned int face_no = 0;
95
96 // TODO[TL]: for anisotropic refinement we will probably need a table of
97 // submatrices with an array for each refine case
98 std::vector<FullMatrix<double>> face_embeddings(
99 1 << (dim - 1),
100 FullMatrix<double>(this->n_dofs_per_face(face_no),
101 this->n_dofs_per_face(face_no)));
102 // TODO: Something goes wrong there. The error of the least squares fit
103 // is to large ...
104 // FETools::compute_face_embedding_matrices(*this, face_embeddings.data(), 0,
105 // 0);
106 this->interface_constraints.reinit((1 << (dim - 1)) *
107 this->n_dofs_per_face(face_no),
108 this->n_dofs_per_face(face_no));
109 unsigned int target_row = 0;
110 for (const auto &face_embedding : face_embeddings)
111 for (unsigned int i = 0; i < face_embedding.m(); ++i)
112 {
113 for (unsigned int j = 0; j < face_embedding.n(); ++j)
114 this->interface_constraints(target_row, j) = face_embedding(i, j);
115 ++target_row;
116 }
117 }
118
119
120
121 template <int dim>
122 std::string
get_name() const123 FE_ABF<dim>::get_name() const
124 {
125 // note that the
126 // FETools::get_fe_by_name
127 // function depends on the
128 // particular format of the string
129 // this function returns, so they
130 // have to be kept in synch
131
132 std::ostringstream namebuf;
133
134 namebuf << "FE_ABF<" << dim << ">(" << rt_order << ")";
135
136 return namebuf.str();
137 }
138
139
140
141 template <int dim>
142 std::unique_ptr<FiniteElement<dim, dim>>
clone() const143 FE_ABF<dim>::clone() const
144 {
145 return std::make_unique<FE_ABF<dim>>(rt_order);
146 }
147
148
149 //---------------------------------------------------------------------------
150 // Auxiliary and internal functions
151 //---------------------------------------------------------------------------
152
153
154
155 // Version for 2d and higher. See above for 1d version
156 template <int dim>
157 void
initialize_support_points(const unsigned int deg)158 FE_ABF<dim>::initialize_support_points(const unsigned int deg)
159 {
160 QGauss<dim> cell_quadrature(deg + 2);
161 const unsigned int n_interior_points = cell_quadrature.size();
162
163 // TODO: the implementation makes the assumption that all faces have the
164 // same number of dofs
165 AssertDimension(this->n_unique_faces(), 1);
166 const unsigned int face_no = 0;
167
168 unsigned int n_face_points = (dim > 1) ? 1 : 0;
169 // compute (deg+1)^(dim-1)
170 for (unsigned int d = 1; d < dim; ++d)
171 n_face_points *= deg + 1;
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
178 // These might be required when the faces contribution is computed
179 // Therefore they will be initialized at this point.
180 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials_abf;
181
182 // Generate x_1^{i} x_2^{r+1} ...
183 for (unsigned int dd = 0; dd < dim; ++dd)
184 {
185 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
186 for (unsigned int d = 0; d < dim; ++d)
187 poly[d].push_back(Polynomials::Monomial<double>(deg + 1));
188 poly[dd] = Polynomials::Monomial<double>::generate_complete_basis(deg);
189
190 polynomials_abf[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
191 }
192
193 // Number of the point being entered
194 unsigned int current = 0;
195
196 if (dim > 1)
197 {
198 QGauss<dim - 1> face_points(deg + 1);
199 TensorProductPolynomials<dim - 1> legendre =
200 Polynomials::Legendre::generate_complete_basis(deg);
201
202 boundary_weights.reinit(n_face_points, legendre.n());
203
204 // Assert (face_points.size() == this->n_dofs_per_face(),
205 // ExcInternalError());
206
207 for (unsigned int k = 0; k < n_face_points; ++k)
208 {
209 this->generalized_face_support_points[face_no][k] =
210 face_points.point(k);
211 // Compute its quadrature
212 // contribution for each
213 // moment.
214 for (unsigned int i = 0; i < legendre.n(); ++i)
215 {
216 boundary_weights(k, i) =
217 face_points.weight(k) *
218 legendre.compute_value(i, face_points.point(k));
219 }
220 }
221
222 Quadrature<dim> faces =
223 QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
224 face_points);
225 for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
226 ++current)
227 {
228 // Enter the support point
229 // into the vector
230 this->generalized_support_points[current] = faces.point(current);
231 }
232
233
234 // Now initialize edge interior weights for the ABF elements.
235 // These are completely independent from the usual edge moments. They
236 // stem from applying the Gauss theorem to the nodal values, which
237 // was necessary to cast the ABF elements into the deal.II framework
238 // for vector valued elements.
239 boundary_weights_abf.reinit(faces.size(), polynomials_abf[0]->n() * dim);
240 for (unsigned int k = 0; k < faces.size(); ++k)
241 {
242 for (unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
243 {
244 boundary_weights_abf(k, i) =
245 polynomials_abf[i % dim]->compute_value(i / dim,
246 faces.point(k)) *
247 faces.weight(k);
248 }
249 }
250 }
251
252 // Create Legendre basis for the
253 // space D_xi Q_k
254 if (deg > 0)
255 {
256 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
257
258 for (unsigned int dd = 0; dd < dim; ++dd)
259 {
260 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
261 for (unsigned int d = 0; d < dim; ++d)
262 poly[d] = Polynomials::Legendre::generate_complete_basis(deg);
263 poly[dd] = Polynomials::Legendre::generate_complete_basis(deg - 1);
264
265 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
266 }
267
268 interior_weights.reinit(
269 TableIndices<3>(n_interior_points, polynomials[0]->n(), dim));
270
271 for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
272 {
273 for (unsigned int i = 0; i < polynomials[0]->n(); ++i)
274 for (unsigned int d = 0; d < dim; ++d)
275 interior_weights(k, i, d) =
276 cell_quadrature.weight(k) *
277 polynomials[d]->compute_value(i, cell_quadrature.point(k));
278 }
279 }
280
281
282 // Decouple the creation of the generalized support points
283 // from computation of interior weights.
284 for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
285 this->generalized_support_points[current++] = cell_quadrature.point(k);
286
287 // Additional functionality for the ABF elements
288 // TODO: Here the canonical extension of the principle
289 // behind the ABF elements is implemented. It is unclear,
290 // if this really leads to the ABF spaces in 3D!
291 interior_weights_abf.reinit(TableIndices<3>(cell_quadrature.size(),
292 polynomials_abf[0]->n() * dim,
293 dim));
294 Tensor<1, dim> poly_grad;
295
296 for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
297 {
298 for (unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
299 {
300 poly_grad =
301 polynomials_abf[i % dim]->compute_grad(i / dim,
302 cell_quadrature.point(k)) *
303 cell_quadrature.weight(k);
304 // The minus sign comes from the use of the Gauss theorem to replace
305 // the divergence.
306 for (unsigned int d = 0; d < dim; ++d)
307 interior_weights_abf(k, i, d) = -poly_grad[d];
308 }
309 }
310
311 Assert(current == this->generalized_support_points.size(),
312 ExcInternalError());
313 }
314
315
316
317 // This function is the same Raviart-Thomas interpolation performed by
318 // interpolate. Still, we cannot use interpolate, since it was written
319 // for smooth functions. The functions interpolated here are not
320 // smooth, maybe even not continuous. Therefore, we must double the
321 // number of quadrature points in each direction in order to integrate
322 // only smooth functions.
323
324 // Then again, the interpolated function is chosen such that the
325 // moments coincide with the function to be interpolated.
326
327 template <int dim>
328 void
initialize_restriction()329 FE_ABF<dim>::initialize_restriction()
330 {
331 if (dim == 1)
332 {
333 unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
334 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
335 ++i)
336 this->restriction[iso][i].reinit(0, 0);
337 return;
338 }
339 unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
340 QGauss<dim - 1> q_base(rt_order + 1);
341 const unsigned int n_face_points = q_base.size();
342 // First, compute interpolation on
343 // subfaces
344 for (unsigned int face : GeometryInfo<dim>::face_indices())
345 {
346 // The shape functions of the
347 // child cell are evaluated
348 // in the quadrature points
349 // of a full face.
350 Quadrature<dim> q_face =
351 QProjector<dim>::project_to_face(this->reference_cell_type(),
352 q_base,
353 face);
354 // Store shape values, since the
355 // evaluation suffers if not
356 // ordered by point
357 Table<2, double> cached_values_face(this->n_dofs_per_cell(),
358 q_face.size());
359 for (unsigned int k = 0; k < q_face.size(); ++k)
360 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
361 cached_values_face(i, k) = this->shape_value_component(
362 i, q_face.point(k), GeometryInfo<dim>::unit_normal_direction[face]);
363
364 for (unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
365 ++sub)
366 {
367 // The weight functions for
368 // the coarse face are
369 // evaluated on the subface
370 // only.
371 Quadrature<dim> q_sub = QProjector<dim>::project_to_subface(
372 this->reference_cell_type(), q_base, face, sub);
373 const unsigned int child = GeometryInfo<dim>::child_cell_on_face(
374 RefinementCase<dim>::isotropic_refinement, face, sub);
375
376 // On a certain face, we must
377 // compute the moments of ALL
378 // fine level functions with
379 // the coarse level weight
380 // functions belonging to
381 // that face. Due to the
382 // orthogonalization process
383 // when building the shape
384 // functions, these weights
385 // are equal to the
386 // corresponding shape
387 // functions.
388 for (unsigned int k = 0; k < n_face_points; ++k)
389 for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
390 ++i_child)
391 for (unsigned int i_face = 0;
392 i_face < this->n_dofs_per_face(face);
393 ++i_face)
394 {
395 // The quadrature
396 // weights on the
397 // subcell are NOT
398 // transformed, so we
399 // have to do it here.
400 this->restriction[iso][child](
401 face * this->n_dofs_per_face(face) + i_face, i_child) +=
402 Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
403 cached_values_face(i_child, k) *
404 this->shape_value_component(
405 face * this->n_dofs_per_face(face) + i_face,
406 q_sub.point(k),
407 GeometryInfo<dim>::unit_normal_direction[face]);
408 }
409 }
410 }
411
412 if (rt_order == 0)
413 return;
414
415 // Create Legendre basis for the
416 // space D_xi Q_k. Here, we cannot
417 // use the shape functions
418 std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
419 for (unsigned int dd = 0; dd < dim; ++dd)
420 {
421 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
422 for (unsigned int d = 0; d < dim; ++d)
423 poly[d] = Polynomials::Legendre::generate_complete_basis(rt_order);
424 poly[dd] = Polynomials::Legendre::generate_complete_basis(rt_order - 1);
425
426 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
427 }
428
429 // TODO: the implementation makes the assumption that all faces have the
430 // same number of dofs
431 AssertDimension(this->n_unique_faces(), 1);
432 const unsigned int face_no = 0;
433
434 QGauss<dim> q_cell(rt_order + 1);
435 const unsigned int start_cell_dofs =
436 GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
437
438 // Store shape values, since the
439 // evaluation suffers if not
440 // ordered by point
441 Table<3, double> cached_values_cell(this->n_dofs_per_cell(),
442 q_cell.size(),
443 dim);
444 for (unsigned int k = 0; k < q_cell.size(); ++k)
445 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
446 for (unsigned int d = 0; d < dim; ++d)
447 cached_values_cell(i, k, d) =
448 this->shape_value_component(i, q_cell.point(k), d);
449
450 for (unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
451 ++child)
452 {
453 Quadrature<dim> q_sub =
454 QProjector<dim>::project_to_child(this->reference_cell_type(),
455 q_cell,
456 child);
457
458 for (unsigned int k = 0; k < q_sub.size(); ++k)
459 for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
460 ++i_child)
461 for (unsigned int d = 0; d < dim; ++d)
462 for (unsigned int i_weight = 0; i_weight < polynomials[d]->n();
463 ++i_weight)
464 {
465 this->restriction[iso][child](start_cell_dofs + i_weight * dim +
466 d,
467 i_child) +=
468 q_sub.weight(k) * cached_values_cell(i_child, k, d) *
469 polynomials[d]->compute_value(i_weight, q_sub.point(k));
470 }
471 }
472 }
473
474
475
476 template <int dim>
477 std::vector<unsigned int>
get_dpo_vector(const unsigned int rt_order)478 FE_ABF<dim>::get_dpo_vector(const unsigned int rt_order)
479 {
480 if (dim == 1)
481 {
482 Assert(false, ExcImpossibleInDim(1));
483 return std::vector<unsigned int>();
484 }
485
486 // the element is face-based (not
487 // to be confused with George
488 // W. Bush's Faith Based
489 // Initiative...), and we have
490 // (rt_order+1)^(dim-1) DoFs per face
491 unsigned int dofs_per_face = 1;
492 for (int d = 0; d < dim - 1; ++d)
493 dofs_per_face *= rt_order + 1;
494
495 // and then there are interior dofs
496 const unsigned int interior_dofs = dim * (rt_order + 1) * dofs_per_face;
497
498 std::vector<unsigned int> dpo(dim + 1);
499 dpo[dim - 1] = dofs_per_face;
500 dpo[dim] = interior_dofs;
501
502 return dpo;
503 }
504
505 //---------------------------------------------------------------------------
506 // Data field initialization
507 //---------------------------------------------------------------------------
508
509 template <int dim>
510 bool
has_support_on_face(const unsigned int shape_index,const unsigned int face_index) const511 FE_ABF<dim>::has_support_on_face(const unsigned int shape_index,
512 const unsigned int face_index) const
513 {
514 AssertIndexRange(shape_index, this->n_dofs_per_cell());
515 AssertIndexRange(face_index, GeometryInfo<dim>::faces_per_cell);
516
517 // Return computed values if we
518 // know them easily. Otherwise, it
519 // is always safe to return true.
520 switch (rt_order)
521 {
522 case 0:
523 {
524 switch (dim)
525 {
526 case 2:
527 {
528 // only on the one
529 // non-adjacent face
530 // are the values
531 // actually zero. list
532 // these in a table
533 return (face_index !=
534 GeometryInfo<dim>::opposite_face[shape_index]);
535 }
536
537 default:
538 return true;
539 }
540 }
541
542 default: // other rt_order
543 return true;
544 }
545
546 return true;
547 }
548
549
550
551 template <int dim>
552 void
convert_generalized_support_point_values_to_dof_values(const std::vector<Vector<double>> & support_point_values,std::vector<double> & nodal_values) const553 FE_ABF<dim>::convert_generalized_support_point_values_to_dof_values(
554 const std::vector<Vector<double>> &support_point_values,
555 std::vector<double> & nodal_values) const
556 {
557 Assert(support_point_values.size() == this->generalized_support_points.size(),
558 ExcDimensionMismatch(support_point_values.size(),
559 this->generalized_support_points.size()));
560 Assert(support_point_values[0].size() == this->n_components(),
561 ExcDimensionMismatch(support_point_values[0].size(),
562 this->n_components()));
563 Assert(nodal_values.size() == this->n_dofs_per_cell(),
564 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
565
566 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
567
568 const unsigned int n_face_points = boundary_weights.size(0);
569 for (unsigned int face : GeometryInfo<dim>::face_indices())
570 for (unsigned int k = 0; k < n_face_points; ++k)
571 for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
572 {
573 nodal_values[i + face * this->n_dofs_per_face(face)] +=
574 boundary_weights(k, i) *
575 support_point_values[face * n_face_points + k][GeometryInfo<
576 dim>::unit_normal_direction[face]];
577 }
578
579 // TODO: the implementation makes the assumption that all faces have the
580 // same number of dofs
581 AssertDimension(this->n_unique_faces(), 1);
582 const unsigned int face_no = 0;
583
584 const unsigned int start_cell_dofs =
585 GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
586 const unsigned int start_cell_points =
587 GeometryInfo<dim>::faces_per_cell * n_face_points;
588
589 for (unsigned int k = 0; k < interior_weights.size(0); ++k)
590 for (unsigned int i = 0; i < interior_weights.size(1); ++i)
591 for (unsigned int d = 0; d < dim; ++d)
592 nodal_values[start_cell_dofs + i * dim + d] +=
593 interior_weights(k, i, d) *
594 support_point_values[k + start_cell_points][d];
595
596 const unsigned int start_abf_dofs =
597 start_cell_dofs + interior_weights.size(1) * dim;
598
599 // Cell integral of ABF terms
600 for (unsigned int k = 0; k < interior_weights_abf.size(0); ++k)
601 for (unsigned int i = 0; i < interior_weights_abf.size(1); ++i)
602 for (unsigned int d = 0; d < dim; ++d)
603 nodal_values[start_abf_dofs + i] +=
604 interior_weights_abf(k, i, d) *
605 support_point_values[k + start_cell_points][d];
606
607 // Face integral of ABF terms
608 for (unsigned int face : GeometryInfo<dim>::face_indices())
609 {
610 const double n_orient = GeometryInfo<dim>::unit_normal_orientation[face];
611 for (unsigned int fp = 0; fp < n_face_points; ++fp)
612 {
613 // TODO: Check what the face_orientation, face_flip and face_rotation
614 // have to be in 3D
615 unsigned int k = QProjector<dim>::DataSetDescriptor::face(
616 this->reference_cell_type(),
617 face,
618 false,
619 false,
620 false,
621 n_face_points);
622 for (unsigned int i = 0; i < boundary_weights_abf.size(1); ++i)
623 nodal_values[start_abf_dofs + i] +=
624 n_orient * boundary_weights_abf(k + fp, i) *
625 support_point_values[face * n_face_points + fp][GeometryInfo<
626 dim>::unit_normal_direction[face]];
627 }
628 }
629
630 // TODO: Check if this "correction" can be removed.
631 for (unsigned int i = 0; i < boundary_weights_abf.size(1); ++i)
632 if (std::fabs(nodal_values[start_abf_dofs + i]) < 1.0e-16)
633 nodal_values[start_abf_dofs + i] = 0.0;
634 }
635
636
637
638 template <int dim>
639 std::size_t
memory_consumption() const640 FE_ABF<dim>::memory_consumption() const
641 {
642 Assert(false, ExcNotImplemented());
643 return 0;
644 }
645
646
647
648 /*-------------- Explicit Instantiations -------------------------------*/
649 #include "fe_abf.inst"
650
651 DEAL_II_NAMESPACE_CLOSE
652