Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 <algorithm>
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
117 derivative_norm(const Derivative &d);
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
231 derivative_norm(const Derivative &d);
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
307 return std::max(std::fabs(eigenvalues[0]), std::fabs(eigenvalues[1]));
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 (Utilities::fixed_power<3>(s[0][0]) +
440 Utilities::fixed_power<3>(s[1][1]) +
441 Utilities::fixed_power<3>(s[2][2]) + 3. * s[0][0] * (ss01 + ss02) +
442 3. * s[1][1] * (ss01 + ss12) + 3. * s[2][2] * (ss02 + ss12) +
443 6. * s[0][1] * s[0][2] * s[1][2]) /
444 3.;
445
446 const double R = std::sqrt(4. * J2 / 3.);
447
448 double EE[3] = {0, 0, 0};
449 // the eigenvalues are away from
450 // @p{am} in the order of R. thus,
451 // if R<<AM, then we have the
452 // degenerate case with three
453 // identical eigenvalues. check
454 // this first
455 if (R <= 1e-14 * std::fabs(am))
456 EE[0] = EE[1] = EE[2] = am;
457 else
458 {
459 // at least two eigenvalues are
460 // distinct
461 const double R3 = R * R * R;
462 const double XX = 4. * J3 / R3;
463 const double YY = 1. - std::fabs(XX);
464
465 Assert(YY > -1e-14, ExcInternalError());
466
467 if (YY < 0)
468 {
469 // two roots are equal
470 const double a = (XX > 0 ? -1. : 1.) * R / 2;
471 EE[0] = EE[1] = am + a;
472 EE[2] = am - 2. * a;
473 }
474 else
475 {
476 const double theta = std::acos(XX) / 3.;
477 EE[0] = am + R * std::cos(theta);
478 EE[1] = am + R * std::cos(theta + 2. / 3. * numbers::PI);
479 EE[2] = am + R * std::cos(theta + 4. / 3. * numbers::PI);
480 }
481 }
482
483 return std::max({std::fabs(EE[0]), 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 "numerics/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:320
Definition point.h:113
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:140
const FEValuesType & get_present_fe_values() const
Definition fe_values.h:695
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:296
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
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.
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)
constexpr char T
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)
constexpr double PI
Definition numbers.h:240
::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 > &)