Reference documentation for deal.II version GIT 7deb6c54a6 2023-06-09 18:50: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 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
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 
40 #include <deal.II/lac/la_vector.h>
47 #include <deal.II/lac/vector.h>
48 
50 
51 #include <cmath>
52 
54 
55 
56 
57 namespace
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>
80  class Gradient
81  {
82  public:
87  static const UpdateFlags update_flags;
88 
94 
100 
107  template <class InputVector, int spacedim>
108  static ProjectedDerivative
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>
138  inline typename Gradient<dim>::ProjectedDerivative
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:
198  static const UpdateFlags update_flags;
199 
205 
211 
218  template <class InputVector, int spacedim>
219  static ProjectedDerivative
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 
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:
536  static const UpdateFlags update_flags;
537 
544 
550 
557  template <class InputVector, int spacedim>
558  static ProjectedDerivative
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
712 namespace DerivativeApproximation
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 
733 namespace DerivativeApproximation
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
775  Tensor<2, dim> Y;
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
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) =
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 
1014 namespace DerivativeApproximation
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(),
1037  ExcNotImplemented());
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,
1045  derivative_norm);
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(),
1070  ExcNotImplemented());
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,
1078  derivative_norm);
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 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
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3270
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type >> &gradients) const
Definition: fe_values.cc:3409
const FiniteElement< dim, spacedim > & get_fe() const
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type >> &hessians) const
Definition: fe_values.cc:3520
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:109
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:475
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:476
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInsufficientDirections()
static ::ExceptionBase & ExcVectorLengthVsNActiveCells(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
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(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_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_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:259
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 > &)