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