Reference documentation for deal.II version 9.3.3
\(\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\}}\)
derivative_approximation.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 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
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
18
21
22#include <deal.II/fe/fe.h>
25
30
35
46#include <deal.II/lac/vector.h>
47
49
50#include <cmath>
51
53
54
55
56namespace
57{
58 template <typename T>
59 inline T
60 sqr(const T t)
61 {
62 return t * t;
63 }
64} // namespace
65
66// --- First the classes and functions that describe individual derivatives ---
67
69{
70 namespace internal
71 {
78 template <int dim>
80 {
81 public:
87
93
99
106 template <class InputVector, int spacedim>
109 const InputVector & solution,
110 const unsigned int component);
111
116 static double
118
126 static void
127 symmetrize(Derivative &derivative_tensor);
128 };
129
130 // static variables
131 template <int dim>
133
134
135 template <int dim>
136 template <class InputVector, int spacedim>
139 const FEValues<dim, spacedim> &fe_values,
140 const InputVector & solution,
141 const unsigned int component)
142 {
143 if (fe_values.get_fe().n_components() == 1)
144 {
145 std::vector<typename InputVector::value_type> values(1);
146 fe_values.get_function_values(solution, values);
147 return values[0];
148 }
149 else
150 {
151 std::vector<Vector<typename InputVector::value_type>> values(
152 1,
154 fe_values.get_fe().n_components()));
155 fe_values.get_function_values(solution, values);
156 return values[0](component);
157 }
158 }
159
160
161
162 template <int dim>
163 inline double
165 {
166 double s = 0;
167 for (unsigned int i = 0; i < dim; ++i)
168 s += d[i] * d[i];
169 return std::sqrt(s);
170 }
171
172
173
174 template <int dim>
175 inline void
177 {
178 // nothing to do here
179 }
180
181
182
189 template <int dim>
191 {
192 public:
198
204
210
217 template <class InputVector, int spacedim>
220 const InputVector & solution,
221 const unsigned int component);
222
230 static double
232
242 static void
243 symmetrize(Derivative &derivative_tensor);
244 };
245
246 template <int dim>
248
249
250 template <int dim>
251 template <class InputVector, int spacedim>
254 const FEValues<dim, spacedim> &fe_values,
255 const InputVector & solution,
256 const unsigned int component)
257 {
258 if (fe_values.get_fe().n_components() == 1)
259 {
260 std::vector<Tensor<1, dim, typename InputVector::value_type>> values(
261 1);
262 fe_values.get_function_gradients(solution, values);
263 return ProjectedDerivative(values[0]);
264 }
265 else
266 {
267 std::vector<
268 std::vector<Tensor<1, dim, typename InputVector::value_type>>>
269 values(
270 1,
272 fe_values.get_fe().n_components()));
273 fe_values.get_function_gradients(solution, values);
274 return ProjectedDerivative(values[0][component]);
275 }
276 }
277
278
279
280 template <>
281 inline double
283 {
284 return std::fabs(d[0][0]);
285 }
286
287
288
289 template <>
290 inline double
292 {
293 // note that d should be a
294 // symmetric 2x2 tensor, so the
295 // eigenvalues are:
296 //
297 // 1/2(a+b\pm\sqrt((a-b)^2+4c^2))
298 //
299 // if the d_11=a, d_22=b,
300 // d_12=d_21=c
301 const double radicand =
302 ::sqr(d[0][0] - d[1][1]) + 4 * ::sqr(d[0][1]);
303 const double eigenvalues[2] = {
304 0.5 * (d[0][0] + d[1][1] + std::sqrt(radicand)),
305 0.5 * (d[0][0] + d[1][1] - std::sqrt(radicand))};
306
308 }
309
310
311
312 template <>
313 inline double
315 {
316 /*
317 compute the three eigenvalues of the tensor @p{d} and take the
318 largest. one could use the following maple script to generate C
319 code:
320
321 with(linalg);
322 readlib(C);
323 A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]);
324 E:=eigenvals(A);
325 EE:=vector(3,[E[1],E[2],E[3]]);
326 C(EE);
327
328 Unfortunately, with both optimized and non-optimized output, at some
329 places the code `sqrt(-1.0)' is emitted, and I don't know what
330 Maple intends to do with it. This happens both with Maple4 and
331 Maple5.
332
333 Fortunately, Roger Young provided the following Fortran code, which
334 is transcribed below to C. The code uses an algorithm that uses the
335 invariants of a symmetric matrix. (The translated algorithm is
336 augmented by a test for R>0, since R==0 indicates that all three
337 eigenvalues are equal.)
338
339
340 PROGRAM MAIN
341
342 C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX
343 C (ROGER YOUNG, 2001)
344
345 IMPLICIT NONE
346
347 REAL*8 A11, A12, A13, A22, A23, A33
348 REAL*8 I1, J2, J3, AM
349 REAL*8 S11, S12, S13, S22, S23, S33
350 REAL*8 SS12, SS23, SS13
351 REAL*8 R,R3, XX,YY, THETA
352 REAL*8 A1,A2,A3
353 REAL*8 PI
354 PARAMETER (PI=3.141592653587932384D0)
355 REAL*8 A,B,C, TOL
356 PARAMETER (TOL=1.D-14)
357
358 C DEFINE A TEST MATRIX
359
360 A11 = -1.D0
361 A12 = 5.D0
362 A13 = 3.D0
363 A22 = -2.D0
364 A23 = 0.5D0
365 A33 = 4.D0
366
367
368 I1 = A11 + A22 + A33
369 AM = I1/3.D0
370
371 S11 = A11 - AM
372 S22 = A22 - AM
373 S33 = A33 - AM
374 S12 = A12
375 S13 = A13
376 S23 = A23
377
378 SS12 = S12*S12
379 SS23 = S23*S23
380 SS13 = S13*S13
381
382 J2 = S11*S11 + S22*S22 + S33*S33
383 J2 = J2 + 2.D0*(SS12 + SS23 + SS13)
384 J2 = J2/2.D0
385
386 J3 = S11**3 + S22**3 + S33**3
387 J3 = J3 + 3.D0*S11*(SS12 + SS13)
388 J3 = J3 + 3.D0*S22*(SS12 + SS23)
389 J3 = J3 + 3.D0*S33*(SS13 + SS23)
390 J3 = J3 + 6.D0*S12*S23*S13
391 J3 = J3/3.D0
392
393 R = SQRT(4.D0*J2/3.D0)
394 R3 = R*R*R
395 XX = 4.D0*J3/R3
396
397 YY = 1.D0 - DABS(XX)
398 IF(YY.LE.0.D0)THEN
399 IF(YY.GT.(-TOL))THEN
400 WRITE(6,*)'Equal roots: XX= ',XX
401 A = -(XX/DABS(XX))*SQRT(J2/3.D0)
402 B = AM + A
403 C = AM - 2.D0*A
404 WRITE(6,*)B,' (twice) ',C
405 STOP
406 ELSE
407 WRITE(6,*)'Error: XX= ',XX
408 STOP
409 ENDIF
410 ENDIF
411
412 THETA = (ACOS(XX))/3.D0
413
414 A1 = AM + R*COS(THETA)
415 A2 = AM + R*COS(THETA + 2.D0*PI/3.D0)
416 A3 = AM + R*COS(THETA + 4.D0*PI/3.D0)
417
418 WRITE(6,*)A1,A2,A3
419
420 STOP
421 END
422
423 */
424
425 const double am = trace(d) / 3.;
426
427 // s := d - trace(d) I
428 Tensor<2, 3> s = d;
429 for (unsigned int i = 0; i < 3; ++i)
430 s[i][i] -= am;
431
432 const double ss01 = s[0][1] * s[0][1], ss12 = s[1][2] * s[1][2],
433 ss02 = s[0][2] * s[0][2];
434
435 const double J2 = (s[0][0] * s[0][0] + s[1][1] * s[1][1] +
436 s[2][2] * s[2][2] + 2 * (ss01 + ss02 + ss12)) /
437 2.;
438 const double J3 =
439 (std::pow(s[0][0], 3) + std::pow(s[1][1], 3) + std::pow(s[2][2], 3) +
440 3. * s[0][0] * (ss01 + ss02) + 3. * s[1][1] * (ss01 + ss12) +
441 3. * s[2][2] * (ss02 + ss12) + 6. * s[0][1] * s[0][2] * s[1][2]) /
442 3.;
443
444 const double R = std::sqrt(4. * J2 / 3.);
445
446 double EE[3] = {0, 0, 0};
447 // the eigenvalues are away from
448 // @p{am} in the order of R. thus,
449 // if R<<AM, then we have the
450 // degenerate case with three
451 // identical eigenvalues. check
452 // this first
453 if (R <= 1e-14 * std::fabs(am))
454 EE[0] = EE[1] = EE[2] = am;
455 else
456 {
457 // at least two eigenvalues are
458 // distinct
459 const double R3 = R * R * R;
460 const double XX = 4. * J3 / R3;
461 const double YY = 1. - std::fabs(XX);
462
463 Assert(YY > -1e-14, ExcInternalError());
464
465 if (YY < 0)
466 {
467 // two roots are equal
468 const double a = (XX > 0 ? -1. : 1.) * R / 2;
469 EE[0] = EE[1] = am + a;
470 EE[2] = am - 2. * a;
471 }
472 else
473 {
474 const double theta = std::acos(XX) / 3.;
475 EE[0] = am + R * std::cos(theta);
476 EE[1] = am + R * std::cos(theta + 2. / 3. * numbers::PI);
477 EE[2] = am + R * std::cos(theta + 4. / 3. * numbers::PI);
478 }
479 }
480
481 return std::max(std::fabs(EE[0]),
482 std::max(std::fabs(EE[1]), std::fabs(EE[2])));
483 }
484
485
486
487 template <int dim>
488 inline double
490 {
491 // computing the spectral norm is
492 // not so simple in general. it is
493 // feasible for dim==3 as shown
494 // above, since then there are
495 // still closed form expressions of
496 // the roots of the characteristic
497 // polynomial, and they can easily
498 // be computed using
499 // maple. however, for higher
500 // dimensions, some other method
501 // needs to be employed. maybe some
502 // steps of the power method would
503 // suffice?
504 Assert(false, ExcNotImplemented());
505 return 0;
506 }
507
508
509
510 template <int dim>
511 inline void
513 {
514 // symmetrize non-diagonal entries
515 for (unsigned int i = 0; i < dim; ++i)
516 for (unsigned int j = i + 1; j < dim; ++j)
517 {
518 const double s = (d[i][j] + d[j][i]) / 2;
519 d[i][j] = d[j][i] = s;
520 }
521 }
522
523
524
525 template <int dim>
527 {
528 public:
534
541
547
554 template <class InputVector, int spacedim>
557 const InputVector & solution,
558 const unsigned int component);
559
567 static double
569
579 static void
580 symmetrize(Derivative &derivative_tensor);
581 };
582
583 template <int dim>
585
586
587 template <int dim>
588 template <class InputVector, int spacedim>
591 const FEValues<dim, spacedim> &fe_values,
592 const InputVector & solution,
593 const unsigned int component)
594 {
595 if (fe_values.get_fe().n_components() == 1)
596 {
597 std::vector<Tensor<2, dim, typename InputVector::value_type>> values(
598 1);
599 fe_values.get_function_hessians(solution, values);
600 return ProjectedDerivative(values[0]);
601 }
602 else
603 {
604 std::vector<
605 std::vector<Tensor<2, dim, typename InputVector::value_type>>>
606 values(
607 1,
609 fe_values.get_fe().n_components()));
610 fe_values.get_function_hessians(solution, values);
611 return ProjectedDerivative(values[0][component]);
612 }
613 }
614
615
616
617 template <>
618 inline double
620 {
621 return std::fabs(d[0][0][0]);
622 }
623
624
625
626 template <int dim>
627 inline double
629 {
630 // return the Frobenius-norm. this is a
631 // member function of Tensor<rank_,dim>
632 return d.norm();
633 }
634
635
636 template <int dim>
637 inline void
639 {
640 // symmetrize non-diagonal entries
641
642 // first do it in the case, that i,j,k are
643 // pairwise different (which can onlky happen
644 // in dim >= 3)
645 for (unsigned int i = 0; i < dim; ++i)
646 for (unsigned int j = i + 1; j < dim; ++j)
647 for (unsigned int k = j + 1; k < dim; ++k)
648 {
649 const double s = (d[i][j][k] + d[i][k][j] + d[j][i][k] +
650 d[j][k][i] + d[k][i][j] + d[k][j][i]) /
651 6;
652 d[i][j][k] = d[i][k][j] = d[j][i][k] = d[j][k][i] = d[k][i][j] =
653 d[k][j][i] = s;
654 }
655 // now do the case, where two indices are
656 // equal
657 for (unsigned int i = 0; i < dim; ++i)
658 for (unsigned int j = i + 1; j < dim; ++j)
659 {
660 // case 1: index i (lower one) is
661 // double
662 const double s = (d[i][i][j] + d[i][j][i] + d[j][i][i]) / 3;
663 d[i][i][j] = d[i][j][i] = d[j][i][i] = s;
664
665 // case 2: index j (higher one) is
666 // double
667 const double t = (d[i][j][j] + d[j][i][j] + d[j][j][i]) / 3;
668 d[i][j][j] = d[j][i][j] = d[j][j][i] = t;
669 }
670 }
671
672
673 template <int order, int dim>
675 {
676 public:
682 using DerivDescr = void;
683 };
684
685 template <int dim>
686 class DerivativeSelector<1, dim>
687 {
688 public:
690 };
691
692 template <int dim>
693 class DerivativeSelector<2, dim>
694 {
695 public:
697 };
698
699 template <int dim>
700 class DerivativeSelector<3, dim>
701 {
702 public:
704 };
705 } // namespace internal
706} // namespace DerivativeApproximation
707
708// Dummy structures and dummy function used for WorkStream
710{
711 namespace internal
712 {
713 namespace Assembler
714 {
715 struct Scratch
716 {
717 Scratch() = default;
718 };
719
720 struct CopyData
721 {
722 CopyData() = default;
723 };
724 } // namespace Assembler
725 } // namespace internal
726} // namespace DerivativeApproximation
727
728// --------------- now for the functions that do the actual work --------------
729
731{
732 namespace internal
733 {
738 template <class DerivativeDescription,
739 int dim,
740 class InputVector,
741 int spacedim>
742 void
744 const Mapping<dim, spacedim> & mapping,
745 const DoFHandler<dim, spacedim> &dof_handler,
746 const InputVector & solution,
747 const unsigned int component,
749 typename DerivativeDescription::Derivative &derivative)
750 {
751 QMidpoint<dim> midpoint_rule;
752
753 // create collection objects from
754 // single quadratures, mappings,
755 // and finite elements. if we have
756 // an hp-DoFHandler,
757 // dof_handler.get_fe() returns a
758 // collection of which we do a
759 // shallow copy instead
760 const hp::QCollection<dim> q_collection(midpoint_rule);
761 const hp::FECollection<dim> &fe_collection =
762 dof_handler.get_fe_collection();
763 const hp::MappingCollection<dim> mapping_collection(mapping);
764
765 hp::FEValues<dim> x_fe_midpoint_value(
766 mapping_collection,
767 fe_collection,
768 q_collection,
769 DerivativeDescription::update_flags | update_quadrature_points);
770
771 // matrix Y=sum_i y_i y_i^T
773
774
775 // vector to hold iterators to all
776 // active neighbors of a cell
777 // reserve the maximal number of
778 // active neighbors
779 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
780 active_neighbors;
781
782 active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
784
785 // vector
786 // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
787 // or related type for higher
788 // derivatives
789 typename DerivativeDescription::Derivative projected_derivative;
790
791 // reinit FE values object...
792 x_fe_midpoint_value.reinit(cell);
793 const FEValues<dim> &fe_midpoint_value =
794 x_fe_midpoint_value.get_present_fe_values();
795
796 // ...and get the value of the
797 // projected derivative...
798 const typename DerivativeDescription::ProjectedDerivative
799 this_midpoint_value =
800 DerivativeDescription::get_projected_derivative(fe_midpoint_value,
801 solution,
802 component);
803 // ...and the place where it lives
804 // This needs to be a copy. If it was a reference, it would be changed
805 // after the next `reinit` call of the FEValues object. clang-tidy
806 // complains about this not being a reference, so we suppress the warning.
807 const Point<dim> this_center =
808 fe_midpoint_value.quadrature_point(0); // NOLINT
809
810 // loop over all neighbors and
811 // accumulate the difference
812 // quotients from them. note
813 // that things get a bit more
814 // complicated if the neighbor
815 // is more refined than the
816 // present one
817 //
818 // to make processing simpler,
819 // first collect all neighbor
820 // cells in a vector, and then
821 // collect the data from them
822 GridTools::get_active_neighbors<DoFHandler<dim, spacedim>>(
823 cell, active_neighbors);
824
825 // now loop over all active
826 // neighbors and collect the
827 // data we need
828 auto neighbor_ptr = active_neighbors.begin();
829 for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
830 {
831 const auto neighbor = *neighbor_ptr;
832
833 // reinit FE values object...
834 x_fe_midpoint_value.reinit(neighbor);
835 const FEValues<dim> &neighbor_fe_midpoint_value =
836 x_fe_midpoint_value.get_present_fe_values();
837
838 // ...and get the value of the
839 // solution...
840 const typename DerivativeDescription::ProjectedDerivative
841 neighbor_midpoint_value =
842 DerivativeDescription::get_projected_derivative(
843 neighbor_fe_midpoint_value, solution, component);
844
845 // ...and the place where it lives
846 const Point<dim> &neighbor_center =
847 neighbor_fe_midpoint_value.quadrature_point(0);
848
849
850 // vector for the
851 // normalized
852 // direction between
853 // the centers of two
854 // cells
855
856 Tensor<1, dim> y = neighbor_center - this_center;
857 const double distance = y.norm();
858 // normalize y
859 y /= distance;
860 // *** note that unlike in
861 // the docs, y denotes the
862 // normalized vector
863 // connecting the centers
864 // of the two cells, rather
865 // than the normal
866 // difference! ***
867
868 // add up the
869 // contribution of
870 // this cell to Y
871 for (unsigned int i = 0; i < dim; ++i)
872 for (unsigned int j = 0; j < dim; ++j)
873 Y[i][j] += y[i] * y[j];
874
875 // then update the sum
876 // of difference
877 // quotients
878 typename DerivativeDescription::ProjectedDerivative
879 projected_finite_difference =
880 (neighbor_midpoint_value - this_midpoint_value);
881 projected_finite_difference /= distance;
882
883 projected_derivative += outer_product(y, projected_finite_difference);
884 }
885
886 // can we determine an
887 // approximation of the
888 // gradient for the present
889 // cell? if so, then we need to
890 // have passed over vectors y_i
891 // which span the whole space,
892 // otherwise we would not have
893 // all components of the
894 // gradient
896
897 // compute Y^-1 g
898 const Tensor<2, dim> Y_inverse = invert(Y);
899
900 derivative = Y_inverse * projected_derivative;
901
902 // finally symmetrize the derivative
904 }
905
906
907
913 template <class DerivativeDescription,
914 int dim,
915 class InputVector,
916 int spacedim>
917 void
921 Vector<float>::iterator>> const &cell,
922 const Mapping<dim, spacedim> & mapping,
923 const DoFHandler<dim, spacedim> & dof_handler,
924 const InputVector & solution,
925 const unsigned int component)
926 {
927 // if the cell is not locally owned, then there is nothing to do
928 if (std::get<0>(*cell)->is_locally_owned() == false)
929 *std::get<1>(*cell) = 0;
930 else
931 {
932 typename DerivativeDescription::Derivative derivative;
933 // call the function doing the actual
934 // work on this cell
935 approximate_cell<DerivativeDescription, dim, InputVector, spacedim>(
936 mapping,
937 dof_handler,
938 solution,
939 component,
940 std::get<0>(*cell),
941 derivative);
942
943 // evaluate the norm and fill the vector
944 //*derivative_norm_on_this_cell
945 *std::get<1>(*cell) =
947 }
948 }
949
950
961 template <class DerivativeDescription,
962 int dim,
963 class InputVector,
964 int spacedim>
965 void
967 const DoFHandler<dim, spacedim> &dof_handler,
968 const InputVector & solution,
969 const unsigned int component,
971 {
972 Assert(derivative_norm.size() ==
973 dof_handler.get_triangulation().n_active_cells(),
975 derivative_norm.size(),
976 dof_handler.get_triangulation().n_active_cells()));
977 AssertIndexRange(component, dof_handler.get_fe(0).n_components());
978
979 using Iterators =
980 std::tuple<typename DoFHandler<dim, spacedim>::active_cell_iterator,
983 Iterators(dof_handler.begin_active(), derivative_norm.begin())),
984 end(Iterators(dof_handler.end(), derivative_norm.end()));
985
986 // There is no need for a copier because there is no conflict between
987 // threads to write in derivative_norm. Scratch and CopyData are also
988 // useless.
990 begin,
991 end,
992 [&mapping, &dof_handler, &solution, component](
994 Assembler::Scratch const &,
996 approximate<DerivativeDescription, dim, InputVector, spacedim>(
997 cell, mapping, dof_handler, solution, component);
998 },
999 std::function<void(internal::Assembler::CopyData const &)>(),
1002 }
1003
1004 } // namespace internal
1005
1006} // namespace DerivativeApproximation
1007
1008
1009// ------------------------ finally for the public interface of this namespace
1010
1012{
1013 template <int dim, class InputVector, int spacedim>
1014 void
1016 const DoFHandler<dim, spacedim> &dof_handler,
1017 const InputVector & solution,
1019 const unsigned int component)
1020 {
1021 internal::approximate_derivative<internal::Gradient<dim>, dim>(
1022 mapping, dof_handler, solution, component, derivative_norm);
1023 }
1024
1025
1026 template <int dim, class InputVector, int spacedim>
1027 void
1029 const InputVector & solution,
1031 const unsigned int component)
1032 {
1033 internal::approximate_derivative<internal::Gradient<dim>, dim>(
1035 dof_handler,
1036 solution,
1037 component,
1039 }
1040
1041
1042 template <int dim, class InputVector, int spacedim>
1043 void
1045 const DoFHandler<dim, spacedim> &dof_handler,
1046 const InputVector & solution,
1048 const unsigned int component)
1049 {
1050 internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1051 mapping, dof_handler, solution, component, derivative_norm);
1052 }
1053
1054
1055 template <int dim, class InputVector, int spacedim>
1056 void
1058 const InputVector & solution,
1060 const unsigned int component)
1061 {
1062 internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1064 dof_handler,
1065 solution,
1066 component,
1068 }
1069
1070
1071 template <int dim, int spacedim, class InputVector, int order>
1072 void
1074 const Mapping<dim, spacedim> & mapping,
1075 const DoFHandler<dim, spacedim> &dof,
1076 const InputVector & solution,
1077#ifndef _MSC_VER
1079#else
1081 &cell,
1082#endif
1083 Tensor<order, dim> &derivative,
1084 const unsigned int component)
1085 {
1088 mapping, dof, solution, component, cell, derivative);
1089 }
1090
1091
1092
1093 template <int dim, int spacedim, class InputVector, int order>
1094 void
1096 const DoFHandler<dim, spacedim> &dof,
1097 const InputVector & solution,
1098#ifndef _MSC_VER
1100#else
1102 &cell,
1103#endif
1104 Tensor<order, dim> &derivative,
1105 const unsigned int component)
1106 {
1107 // just call the respective function with Q1 mapping
1109 cell->reference_cell()
1110 .template get_default_linear_mapping<dim, spacedim>(),
1111 dof,
1112 solution,
1113 cell,
1114 derivative,
1115 component);
1116 }
1117
1118
1119
1120 template <int dim, int order>
1121 double
1123 {
1125 derivative_norm(derivative);
1126 }
1127
1128} // namespace DerivativeApproximation
1129
1130
1131// --------------------------- explicit instantiations ---------------------
1132#include "derivative_approximation.inst"
1133
1134
static ProjectedDerivative get_projected_derivative(const FEValues< dim, spacedim > &fe_values, const InputVector &solution, const unsigned int component)
static void symmetrize(Derivative &derivative_tensor)
static double derivative_norm(const Derivative &d)
static ProjectedDerivative get_projected_derivative(const FEValues< dim, spacedim > &fe_values, const InputVector &solution, const unsigned int component)
static void symmetrize(Derivative &derivative_tensor)
static ProjectedDerivative get_projected_derivative(const FEValues< dim, spacedim > &fe_values, const InputVector &solution, const unsigned int component)
static void symmetrize(Derivative &derivative_tensor)
cell_iterator end() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
active_cell_iterator begin_active(const unsigned int level=0) const
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3523
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type > > &gradients) const
Definition: fe_values.cc:3664
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type > > &hessians) const
Definition: fe_values.cc:3777
const FiniteElement< dim, spacedim > & get_fe() const
unsigned int n_components() const
Definition: point.h:111
numbers::NumberTraits< Number >::real_type norm() const
unsigned int n_active_cells() const
Definition: vector.h:110
const FEValuesType & get_present_fe_values() const
Definition: fe_values.h:693
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda > > &cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition: fe_values.cc:293
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
static ::ExceptionBase & ExcInsufficientDirections()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcVectorLengthVsNActiveCells(int arg1, int arg2)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
value_type * iterator
Definition: vector.h:125
void approximate_derivative(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component, Vector< float > &derivative_norm)
void approximate_cell(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component, const typename DoFHandler< dim, spacedim >::active_cell_iterator &cell, typename DerivativeDescription::Derivative &derivative)
void approximate(SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void approximate_gradient(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InputVector &solution, Vector< float > &derivative_norm, const unsigned int component=0)
double derivative_norm(const Tensor< order, dim > &derivative)
void approximate_derivative_tensor(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InputVector &solution, const typename DoFHandler< dim, spacedim >::active_cell_iterator &cell, Tensor< order, dim > &derivative, const unsigned int component=0)
void approximate_second_derivative(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InputVector &solution, Vector< float > &derivative_norm, const unsigned int component=0)
Expression fabs(const Expression &x)
Expression acos(const Expression &x)
@ eigenvalues
Eigenvalue vector is filled.
static const char T
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1252
static constexpr double PI
Definition: numbers.h:231
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)