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