Reference documentation for deal.II version GIT 35969cdc9b 2023-12-09 01:10:02+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\}}\)
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 
17 #include <deal.II/base/utilities.h>
19 
22 
23 #include <deal.II/fe/fe.h>
24 #include <deal.II/fe/fe_values.h>
25 #include <deal.II/fe/mapping.h>
26 
31 
33 #include <deal.II/hp/fe_values.h>
36 
46 #include <deal.II/lac/vector.h>
47 
49 
50 #include <cmath>
51 
53 
54 
55 
56 namespace
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>
79  class Gradient
80  {
81  public:
86  static const UpdateFlags update_flags;
87 
93 
99 
106  template <class InputVector, int spacedim>
107  static ProjectedDerivative
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>
137  inline typename Gradient<dim>::ProjectedDerivative
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:
197  static const UpdateFlags update_flags;
198 
204 
210 
217  template <class InputVector, int spacedim>
218  static ProjectedDerivative
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 
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]),
484  std::max(std::fabs(EE[1]), std::fabs(EE[2])));
485  }
486 
487 
488 
489  template <int dim>
490  inline double
492  {
493  // computing the spectral norm is
494  // not so simple in general. it is
495  // feasible for dim==3 as shown
496  // above, since then there are
497  // still closed form expressions of
498  // the roots of the characteristic
499  // polynomial, and they can easily
500  // be computed using
501  // maple. however, for higher
502  // dimensions, some other method
503  // needs to be employed. maybe some
504  // steps of the power method would
505  // suffice?
506  Assert(false, ExcNotImplemented());
507  return 0;
508  }
509 
510 
511 
512  template <int dim>
513  inline void
515  {
516  // symmetrize non-diagonal entries
517  for (unsigned int i = 0; i < dim; ++i)
518  for (unsigned int j = i + 1; j < dim; ++j)
519  {
520  const double s = (d[i][j] + d[j][i]) / 2;
521  d[i][j] = d[j][i] = s;
522  }
523  }
524 
525 
526 
527  template <int dim>
529  {
530  public:
535  static const UpdateFlags update_flags;
536 
543 
549 
556  template <class InputVector, int spacedim>
557  static ProjectedDerivative
559  const InputVector &solution,
560  const unsigned int component);
561 
569  static double
570  derivative_norm(const Derivative &d);
571 
581  static void
582  symmetrize(Derivative &derivative_tensor);
583  };
584 
585  template <int dim>
587 
588 
589  template <int dim>
590  template <class InputVector, int spacedim>
593  const FEValues<dim, spacedim> &fe_values,
594  const InputVector &solution,
595  const unsigned int component)
596  {
597  if (fe_values.get_fe().n_components() == 1)
598  {
599  std::vector<Tensor<2, dim, typename InputVector::value_type>> values(
600  1);
601  fe_values.get_function_hessians(solution, values);
602  return ProjectedDerivative(values[0]);
603  }
604  else
605  {
606  std::vector<
607  std::vector<Tensor<2, dim, typename InputVector::value_type>>>
608  values(
609  1,
611  fe_values.get_fe().n_components()));
612  fe_values.get_function_hessians(solution, values);
613  return ProjectedDerivative(values[0][component]);
614  }
615  }
616 
617 
618 
619  template <>
620  inline double
622  {
623  return std::fabs(d[0][0][0]);
624  }
625 
626 
627 
628  template <int dim>
629  inline double
631  {
632  // return the Frobenius-norm. this is a
633  // member function of Tensor<rank_,dim>
634  return d.norm();
635  }
636 
637 
638  template <int dim>
639  inline void
641  {
642  // symmetrize non-diagonal entries
643 
644  // first do it in the case, that i,j,k are
645  // pairwise different (which can onlky happen
646  // in dim >= 3)
647  for (unsigned int i = 0; i < dim; ++i)
648  for (unsigned int j = i + 1; j < dim; ++j)
649  for (unsigned int k = j + 1; k < dim; ++k)
650  {
651  const double s = (d[i][j][k] + d[i][k][j] + d[j][i][k] +
652  d[j][k][i] + d[k][i][j] + d[k][j][i]) /
653  6;
654  d[i][j][k] = d[i][k][j] = d[j][i][k] = d[j][k][i] = d[k][i][j] =
655  d[k][j][i] = s;
656  }
657  // now do the case, where two indices are
658  // equal
659  for (unsigned int i = 0; i < dim; ++i)
660  for (unsigned int j = i + 1; j < dim; ++j)
661  {
662  // case 1: index i (lower one) is
663  // double
664  const double s = (d[i][i][j] + d[i][j][i] + d[j][i][i]) / 3;
665  d[i][i][j] = d[i][j][i] = d[j][i][i] = s;
666 
667  // case 2: index j (higher one) is
668  // double
669  const double t = (d[i][j][j] + d[j][i][j] + d[j][j][i]) / 3;
670  d[i][j][j] = d[j][i][j] = d[j][j][i] = t;
671  }
672  }
673 
674 
675  template <int order, int dim>
677  {
678  public:
684  using DerivDescr = void;
685  };
686 
687  template <int dim>
688  class DerivativeSelector<1, dim>
689  {
690  public:
692  };
693 
694  template <int dim>
695  class DerivativeSelector<2, dim>
696  {
697  public:
699  };
700 
701  template <int dim>
702  class DerivativeSelector<3, dim>
703  {
704  public:
706  };
707  } // namespace internal
708 } // namespace DerivativeApproximation
709 
710 // Dummy structures and dummy function used for WorkStream
711 namespace DerivativeApproximation
712 {
713  namespace internal
714  {
715  namespace Assembler
716  {
717  struct Scratch
718  {
719  Scratch() = default;
720  };
721 
722  struct CopyData
723  {
724  CopyData() = default;
725  };
726  } // namespace Assembler
727  } // namespace internal
728 } // namespace DerivativeApproximation
729 
730 // --------------- now for the functions that do the actual work --------------
731 
732 namespace DerivativeApproximation
733 {
734  namespace internal
735  {
740  template <class DerivativeDescription,
741  int dim,
742  class InputVector,
743  int spacedim>
744  void
746  const Mapping<dim, spacedim> &mapping,
747  const DoFHandler<dim, spacedim> &dof_handler,
748  const InputVector &solution,
749  const unsigned int component,
751  typename DerivativeDescription::Derivative &derivative)
752  {
753  QMidpoint<dim> midpoint_rule;
754 
755  // create collection objects from
756  // single quadratures, mappings,
757  // and finite elements. if we have
758  // an hp-DoFHandler,
759  // dof_handler.get_fe() returns a
760  // collection of which we do a
761  // shallow copy instead
762  const hp::QCollection<dim> q_collection(midpoint_rule);
763  const hp::FECollection<dim> &fe_collection =
764  dof_handler.get_fe_collection();
765  const hp::MappingCollection<dim> mapping_collection(mapping);
766 
767  hp::FEValues<dim> x_fe_midpoint_value(
768  mapping_collection,
769  fe_collection,
770  q_collection,
771  DerivativeDescription::update_flags | update_quadrature_points);
772 
773  // matrix Y=sum_i y_i y_i^T
774  Tensor<2, dim> Y;
775 
776 
777  // vector to hold iterators to all
778  // active neighbors of a cell
779  // reserve the maximal number of
780  // active neighbors
781  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
782  active_neighbors;
783 
784  active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
786 
787  // vector
788  // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
789  // or related type for higher
790  // derivatives
791  typename DerivativeDescription::Derivative projected_derivative;
792 
793  // reinit FE values object...
794  x_fe_midpoint_value.reinit(cell);
795  const FEValues<dim> &fe_midpoint_value =
796  x_fe_midpoint_value.get_present_fe_values();
797 
798  // ...and get the value of the
799  // projected derivative...
800  const typename DerivativeDescription::ProjectedDerivative
801  this_midpoint_value =
802  DerivativeDescription::get_projected_derivative(fe_midpoint_value,
803  solution,
804  component);
805  // ...and the place where it lives
806  // This needs to be a copy. If it was a reference, it would be changed
807  // after the next `reinit` call of the FEValues object. clang-tidy
808  // complains about this not being a reference, so we suppress the warning.
809  const Point<dim> this_center =
810  fe_midpoint_value.quadrature_point(0); // NOLINT
811 
812  // loop over all neighbors and
813  // accumulate the difference
814  // quotients from them. note
815  // that things get a bit more
816  // complicated if the neighbor
817  // is more refined than the
818  // present one
819  //
820  // to make processing simpler,
821  // first collect all neighbor
822  // cells in a vector, and then
823  // collect the data from them
824  GridTools::get_active_neighbors<DoFHandler<dim, spacedim>>(
825  cell, active_neighbors);
826 
827  // now loop over all active
828  // neighbors and collect the
829  // data we need
830  auto neighbor_ptr = active_neighbors.begin();
831  for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
832  {
833  const auto neighbor = *neighbor_ptr;
834 
835  // reinit FE values object...
836  x_fe_midpoint_value.reinit(neighbor);
837  const FEValues<dim> &neighbor_fe_midpoint_value =
838  x_fe_midpoint_value.get_present_fe_values();
839 
840  // ...and get the value of the
841  // solution...
842  const typename DerivativeDescription::ProjectedDerivative
843  neighbor_midpoint_value =
844  DerivativeDescription::get_projected_derivative(
845  neighbor_fe_midpoint_value, solution, component);
846 
847  // ...and the place where it lives
848  const Point<dim> &neighbor_center =
849  neighbor_fe_midpoint_value.quadrature_point(0);
850 
851 
852  // vector for the
853  // normalized
854  // direction between
855  // the centers of two
856  // cells
857 
858  Tensor<1, dim> y = neighbor_center - this_center;
859  const double distance = y.norm();
860  // normalize y
861  y /= distance;
862  // *** note that unlike in
863  // the docs, y denotes the
864  // normalized vector
865  // connecting the centers
866  // of the two cells, rather
867  // than the normal
868  // difference! ***
869 
870  // add up the
871  // contribution of
872  // this cell to Y
873  for (unsigned int i = 0; i < dim; ++i)
874  for (unsigned int j = 0; j < dim; ++j)
875  Y[i][j] += y[i] * y[j];
876 
877  // then update the sum
878  // of difference
879  // quotients
880  typename DerivativeDescription::ProjectedDerivative
881  projected_finite_difference =
882  (neighbor_midpoint_value - this_midpoint_value);
883  projected_finite_difference /= distance;
884 
885  projected_derivative += outer_product(y, projected_finite_difference);
886  }
887 
888  // can we determine an
889  // approximation of the
890  // gradient for the present
891  // cell? if so, then we need to
892  // have passed over vectors y_i
893  // which span the whole space,
894  // otherwise we would not have
895  // all components of the
896  // gradient
898 
899  // compute Y^-1 g
900  const Tensor<2, dim> Y_inverse = invert(Y);
901 
902  derivative = Y_inverse * projected_derivative;
903 
904  // finally symmetrize the derivative
906  }
907 
908 
909 
915  template <class DerivativeDescription,
916  int dim,
917  class InputVector,
918  int spacedim>
919  void
921  const SynchronousIterators<
923  Vector<float>::iterator>> &cell,
924  const Mapping<dim, spacedim> &mapping,
925  const DoFHandler<dim, spacedim> &dof_handler,
926  const InputVector &solution,
927  const unsigned int component)
928  {
929  // if the cell is not locally owned, then there is nothing to do
930  if (std::get<0>(*cell)->is_locally_owned() == false)
931  *std::get<1>(*cell) = 0;
932  else
933  {
934  typename DerivativeDescription::Derivative derivative;
935  // call the function doing the actual
936  // work on this cell
937  approximate_cell<DerivativeDescription, dim, InputVector, spacedim>(
938  mapping,
939  dof_handler,
940  solution,
941  component,
942  std::get<0>(*cell),
943  derivative);
944 
945  // evaluate the norm and fill the vector
946  //*derivative_norm_on_this_cell
947  *std::get<1>(*cell) =
949  }
950  }
951 
952 
963  template <class DerivativeDescription,
964  int dim,
965  class InputVector,
966  int spacedim>
967  void
969  const DoFHandler<dim, spacedim> &dof_handler,
970  const InputVector &solution,
971  const unsigned int component,
973  {
974  Assert(derivative_norm.size() ==
975  dof_handler.get_triangulation().n_active_cells(),
977  derivative_norm.size(),
978  dof_handler.get_triangulation().n_active_cells()));
979  AssertIndexRange(component, dof_handler.get_fe(0).n_components());
980 
981  using Iterators =
982  std::tuple<typename DoFHandler<dim, spacedim>::active_cell_iterator,
985  Iterators(dof_handler.begin_active(), derivative_norm.begin())),
986  end(Iterators(dof_handler.end(), derivative_norm.end()));
987 
988  // There is no need for a copier because there is no conflict between
989  // threads to write in derivative_norm. Scratch and CopyData are also
990  // useless.
992  begin,
993  end,
994  [&mapping, &dof_handler, &solution, component](
996  const Assembler::Scratch &,
998  approximate<DerivativeDescription, dim, InputVector, spacedim>(
999  cell, mapping, dof_handler, solution, component);
1000  },
1001  std::function<void(const internal::Assembler::CopyData &)>(),
1004  }
1005 
1006  } // namespace internal
1007 
1008 } // namespace DerivativeApproximation
1009 
1010 
1011 // ------------------------ finally for the public interface of this namespace
1012 
1013 namespace DerivativeApproximation
1014 {
1015  template <int dim, class InputVector, int spacedim>
1016  void
1018  const DoFHandler<dim, spacedim> &dof_handler,
1019  const InputVector &solution,
1021  const unsigned int component)
1022  {
1023  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1024  mapping, dof_handler, solution, component, derivative_norm);
1025  }
1026 
1027 
1028  template <int dim, class InputVector, int spacedim>
1029  void
1031  const InputVector &solution,
1033  const unsigned int component)
1034  {
1035  Assert(!dof_handler.get_triangulation().is_mixed_mesh(),
1036  ExcNotImplemented());
1037  const auto reference_cell =
1038  dof_handler.get_triangulation().get_reference_cells()[0];
1039  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1040  reference_cell.template get_default_linear_mapping<dim, spacedim>(),
1041  dof_handler,
1042  solution,
1043  component,
1044  derivative_norm);
1045  }
1046 
1047 
1048  template <int dim, class InputVector, int spacedim>
1049  void
1051  const DoFHandler<dim, spacedim> &dof_handler,
1052  const InputVector &solution,
1054  const unsigned int component)
1055  {
1056  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1057  mapping, dof_handler, solution, component, derivative_norm);
1058  }
1059 
1060 
1061  template <int dim, class InputVector, int spacedim>
1062  void
1064  const InputVector &solution,
1066  const unsigned int component)
1067  {
1068  Assert(!dof_handler.get_triangulation().is_mixed_mesh(),
1069  ExcNotImplemented());
1070  const auto reference_cell =
1071  dof_handler.get_triangulation().get_reference_cells()[0];
1072  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1073  reference_cell.template get_default_linear_mapping<dim, spacedim>(),
1074  dof_handler,
1075  solution,
1076  component,
1077  derivative_norm);
1078  }
1079 
1080 
1081  template <int dim, int spacedim, class InputVector, int order>
1082  void
1084  const Mapping<dim, spacedim> &mapping,
1085  const DoFHandler<dim, spacedim> &dof,
1086  const InputVector &solution,
1087 #ifndef _MSC_VER
1089 #else
1091  &cell,
1092 #endif
1093  Tensor<order, dim> &derivative,
1094  const unsigned int component)
1095  {
1098  mapping, dof, solution, component, cell, derivative);
1099  }
1100 
1101 
1102 
1103  template <int dim, int spacedim, class InputVector, int order>
1104  void
1106  const DoFHandler<dim, spacedim> &dof,
1107  const InputVector &solution,
1108 #ifndef _MSC_VER
1110 #else
1112  &cell,
1113 #endif
1114  Tensor<order, dim> &derivative,
1115  const unsigned int component)
1116  {
1117  // just call the respective function with Q1 mapping
1119  cell->reference_cell()
1120  .template get_default_linear_mapping<dim, spacedim>(),
1121  dof,
1122  solution,
1123  cell,
1124  derivative,
1125  component);
1126  }
1127 
1128 
1129 
1130  template <int dim, int order>
1131  double
1133  {
1135  derivative_norm(derivative);
1136  }
1137 
1138 } // namespace DerivativeApproximation
1139 
1140 
1141 // --------------------------- explicit instantiations ---------------------
1142 #include "derivative_approximation.inst"
1143 
1144 
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 Triangulation< dim, spacedim > & get_triangulation() const
active_cell_iterator begin_active(const unsigned int level=0) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() 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_values(const ReadVector< Number > &fe_function, std::vector< Number > &values) 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
Definition: point.h:112
numbers::NumberTraits< Number >::real_type norm() const
unsigned int n_active_cells() const
bool is_mixed_mesh() const
const std::vector< ReferenceCell > & get_reference_cells() const
Definition: vector.h:110
value_type * iterator
Definition: vector.h:132
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInsufficientDirections()
static ::ExceptionBase & ExcVectorLengthVsNActiveCells(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
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(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_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)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
@ eigenvalues
Eigenvalue vector is filled.
static const char T
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(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:1270
static constexpr double PI
Definition: numbers.h:260
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr DEAL_II_HOST Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)