Reference documentation for deal.II version GIT b533f29dc3 2022-09-27 06:15:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
18 
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_values.h>
24 #include <deal.II/fe/mapping.h>
25 
30 
32 #include <deal.II/hp/fe_values.h>
35 
39 #include <deal.II/lac/la_vector.h>
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  (std::pow(s[0][0], 3) + std::pow(s[1][1], 3) + std::pow(s[2][2], 3) +
440  3. * s[0][0] * (ss01 + ss02) + 3. * s[1][1] * (ss01 + ss12) +
441  3. * s[2][2] * (ss02 + ss12) + 6. * s[0][1] * s[0][2] * s[1][2]) /
442  3.;
443 
444  const double R = std::sqrt(4. * J2 / 3.);
445 
446  double EE[3] = {0, 0, 0};
447  // the eigenvalues are away from
448  // @p{am} in the order of R. thus,
449  // if R<<AM, then we have the
450  // degenerate case with three
451  // identical eigenvalues. check
452  // this first
453  if (R <= 1e-14 * std::fabs(am))
454  EE[0] = EE[1] = EE[2] = am;
455  else
456  {
457  // at least two eigenvalues are
458  // distinct
459  const double R3 = R * R * R;
460  const double XX = 4. * J3 / R3;
461  const double YY = 1. - std::fabs(XX);
462 
463  Assert(YY > -1e-14, ExcInternalError());
464 
465  if (YY < 0)
466  {
467  // two roots are equal
468  const double a = (XX > 0 ? -1. : 1.) * R / 2;
469  EE[0] = EE[1] = am + a;
470  EE[2] = am - 2. * a;
471  }
472  else
473  {
474  const double theta = std::acos(XX) / 3.;
475  EE[0] = am + R * std::cos(theta);
476  EE[1] = am + R * std::cos(theta + 2. / 3. * numbers::PI);
477  EE[2] = am + R * std::cos(theta + 4. / 3. * numbers::PI);
478  }
479  }
480 
481  return std::max(std::fabs(EE[0]),
482  std::max(std::fabs(EE[1]), std::fabs(EE[2])));
483  }
484 
485 
486 
487  template <int dim>
488  inline double
490  {
491  // computing the spectral norm is
492  // not so simple in general. it is
493  // feasible for dim==3 as shown
494  // above, since then there are
495  // still closed form expressions of
496  // the roots of the characteristic
497  // polynomial, and they can easily
498  // be computed using
499  // maple. however, for higher
500  // dimensions, some other method
501  // needs to be employed. maybe some
502  // steps of the power method would
503  // suffice?
504  Assert(false, ExcNotImplemented());
505  return 0;
506  }
507 
508 
509 
510  template <int dim>
511  inline void
513  {
514  // symmetrize non-diagonal entries
515  for (unsigned int i = 0; i < dim; ++i)
516  for (unsigned int j = i + 1; j < dim; ++j)
517  {
518  const double s = (d[i][j] + d[j][i]) / 2;
519  d[i][j] = d[j][i] = s;
520  }
521  }
522 
523 
524 
525  template <int dim>
527  {
528  public:
533  static const UpdateFlags update_flags;
534 
541 
547 
554  template <class InputVector, int spacedim>
555  static ProjectedDerivative
557  const InputVector & solution,
558  const unsigned int component);
559 
567  static double
568  derivative_norm(const Derivative &d);
569 
579  static void
580  symmetrize(Derivative &derivative_tensor);
581  };
582 
583  template <int dim>
585 
586 
587  template <int dim>
588  template <class InputVector, int spacedim>
591  const FEValues<dim, spacedim> &fe_values,
592  const InputVector & solution,
593  const unsigned int component)
594  {
595  if (fe_values.get_fe().n_components() == 1)
596  {
597  std::vector<Tensor<2, dim, typename InputVector::value_type>> values(
598  1);
599  fe_values.get_function_hessians(solution, values);
600  return ProjectedDerivative(values[0]);
601  }
602  else
603  {
604  std::vector<
605  std::vector<Tensor<2, dim, typename InputVector::value_type>>>
606  values(
607  1,
609  fe_values.get_fe().n_components()));
610  fe_values.get_function_hessians(solution, values);
611  return ProjectedDerivative(values[0][component]);
612  }
613  }
614 
615 
616 
617  template <>
618  inline double
620  {
621  return std::fabs(d[0][0][0]);
622  }
623 
624 
625 
626  template <int dim>
627  inline double
629  {
630  // return the Frobenius-norm. this is a
631  // member function of Tensor<rank_,dim>
632  return d.norm();
633  }
634 
635 
636  template <int dim>
637  inline void
639  {
640  // symmetrize non-diagonal entries
641 
642  // first do it in the case, that i,j,k are
643  // pairwise different (which can onlky happen
644  // in dim >= 3)
645  for (unsigned int i = 0; i < dim; ++i)
646  for (unsigned int j = i + 1; j < dim; ++j)
647  for (unsigned int k = j + 1; k < dim; ++k)
648  {
649  const double s = (d[i][j][k] + d[i][k][j] + d[j][i][k] +
650  d[j][k][i] + d[k][i][j] + d[k][j][i]) /
651  6;
652  d[i][j][k] = d[i][k][j] = d[j][i][k] = d[j][k][i] = d[k][i][j] =
653  d[k][j][i] = s;
654  }
655  // now do the case, where two indices are
656  // equal
657  for (unsigned int i = 0; i < dim; ++i)
658  for (unsigned int j = i + 1; j < dim; ++j)
659  {
660  // case 1: index i (lower one) is
661  // double
662  const double s = (d[i][i][j] + d[i][j][i] + d[j][i][i]) / 3;
663  d[i][i][j] = d[i][j][i] = d[j][i][i] = s;
664 
665  // case 2: index j (higher one) is
666  // double
667  const double t = (d[i][j][j] + d[j][i][j] + d[j][j][i]) / 3;
668  d[i][j][j] = d[j][i][j] = d[j][j][i] = t;
669  }
670  }
671 
672 
673  template <int order, int dim>
675  {
676  public:
682  using DerivDescr = void;
683  };
684 
685  template <int dim>
686  class DerivativeSelector<1, dim>
687  {
688  public:
690  };
691 
692  template <int dim>
693  class DerivativeSelector<2, dim>
694  {
695  public:
697  };
698 
699  template <int dim>
700  class DerivativeSelector<3, dim>
701  {
702  public:
704  };
705  } // namespace internal
706 } // namespace DerivativeApproximation
707 
708 // Dummy structures and dummy function used for WorkStream
709 namespace DerivativeApproximation
710 {
711  namespace internal
712  {
713  namespace Assembler
714  {
715  struct Scratch
716  {
717  Scratch() = default;
718  };
719 
720  struct CopyData
721  {
722  CopyData() = default;
723  };
724  } // namespace Assembler
725  } // namespace internal
726 } // namespace DerivativeApproximation
727 
728 // --------------- now for the functions that do the actual work --------------
729 
730 namespace DerivativeApproximation
731 {
732  namespace internal
733  {
738  template <class DerivativeDescription,
739  int dim,
740  class InputVector,
741  int spacedim>
742  void
744  const Mapping<dim, spacedim> & mapping,
745  const DoFHandler<dim, spacedim> &dof_handler,
746  const InputVector & solution,
747  const unsigned int component,
749  typename DerivativeDescription::Derivative &derivative)
750  {
751  QMidpoint<dim> midpoint_rule;
752 
753  // create collection objects from
754  // single quadratures, mappings,
755  // and finite elements. if we have
756  // an hp-DoFHandler,
757  // dof_handler.get_fe() returns a
758  // collection of which we do a
759  // shallow copy instead
760  const hp::QCollection<dim> q_collection(midpoint_rule);
761  const hp::FECollection<dim> &fe_collection =
762  dof_handler.get_fe_collection();
763  const hp::MappingCollection<dim> mapping_collection(mapping);
764 
765  hp::FEValues<dim> x_fe_midpoint_value(
766  mapping_collection,
767  fe_collection,
768  q_collection,
769  DerivativeDescription::update_flags | update_quadrature_points);
770 
771  // matrix Y=sum_i y_i y_i^T
772  Tensor<2, dim> Y;
773 
774 
775  // vector to hold iterators to all
776  // active neighbors of a cell
777  // reserve the maximal number of
778  // active neighbors
779  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
780  active_neighbors;
781 
782  active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
784 
785  // vector
786  // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
787  // or related type for higher
788  // derivatives
789  typename DerivativeDescription::Derivative projected_derivative;
790 
791  // reinit FE values object...
792  x_fe_midpoint_value.reinit(cell);
793  const FEValues<dim> &fe_midpoint_value =
794  x_fe_midpoint_value.get_present_fe_values();
795 
796  // ...and get the value of the
797  // projected derivative...
798  const typename DerivativeDescription::ProjectedDerivative
799  this_midpoint_value =
800  DerivativeDescription::get_projected_derivative(fe_midpoint_value,
801  solution,
802  component);
803  // ...and the place where it lives
804  // This needs to be a copy. If it was a reference, it would be changed
805  // after the next `reinit` call of the FEValues object. clang-tidy
806  // complains about this not being a reference, so we suppress the warning.
807  const Point<dim> this_center =
808  fe_midpoint_value.quadrature_point(0); // NOLINT
809 
810  // loop over all neighbors and
811  // accumulate the difference
812  // quotients from them. note
813  // that things get a bit more
814  // complicated if the neighbor
815  // is more refined than the
816  // present one
817  //
818  // to make processing simpler,
819  // first collect all neighbor
820  // cells in a vector, and then
821  // collect the data from them
822  GridTools::get_active_neighbors<DoFHandler<dim, spacedim>>(
823  cell, active_neighbors);
824 
825  // now loop over all active
826  // neighbors and collect the
827  // data we need
828  auto neighbor_ptr = active_neighbors.begin();
829  for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
830  {
831  const auto neighbor = *neighbor_ptr;
832 
833  // reinit FE values object...
834  x_fe_midpoint_value.reinit(neighbor);
835  const FEValues<dim> &neighbor_fe_midpoint_value =
836  x_fe_midpoint_value.get_present_fe_values();
837 
838  // ...and get the value of the
839  // solution...
840  const typename DerivativeDescription::ProjectedDerivative
841  neighbor_midpoint_value =
842  DerivativeDescription::get_projected_derivative(
843  neighbor_fe_midpoint_value, solution, component);
844 
845  // ...and the place where it lives
846  const Point<dim> &neighbor_center =
847  neighbor_fe_midpoint_value.quadrature_point(0);
848 
849 
850  // vector for the
851  // normalized
852  // direction between
853  // the centers of two
854  // cells
855 
856  Tensor<1, dim> y = neighbor_center - this_center;
857  const double distance = y.norm();
858  // normalize y
859  y /= distance;
860  // *** note that unlike in
861  // the docs, y denotes the
862  // normalized vector
863  // connecting the centers
864  // of the two cells, rather
865  // than the normal
866  // difference! ***
867 
868  // add up the
869  // contribution of
870  // this cell to Y
871  for (unsigned int i = 0; i < dim; ++i)
872  for (unsigned int j = 0; j < dim; ++j)
873  Y[i][j] += y[i] * y[j];
874 
875  // then update the sum
876  // of difference
877  // quotients
878  typename DerivativeDescription::ProjectedDerivative
879  projected_finite_difference =
880  (neighbor_midpoint_value - this_midpoint_value);
881  projected_finite_difference /= distance;
882 
883  projected_derivative += outer_product(y, projected_finite_difference);
884  }
885 
886  // can we determine an
887  // approximation of the
888  // gradient for the present
889  // cell? if so, then we need to
890  // have passed over vectors y_i
891  // which span the whole space,
892  // otherwise we would not have
893  // all components of the
894  // gradient
896 
897  // compute Y^-1 g
898  const Tensor<2, dim> Y_inverse = invert(Y);
899 
900  derivative = Y_inverse * projected_derivative;
901 
902  // finally symmetrize the derivative
904  }
905 
906 
907 
913  template <class DerivativeDescription,
914  int dim,
915  class InputVector,
916  int spacedim>
917  void
921  Vector<float>::iterator>> const &cell,
922  const Mapping<dim, spacedim> & mapping,
923  const DoFHandler<dim, spacedim> & dof_handler,
924  const InputVector & solution,
925  const unsigned int component)
926  {
927  // if the cell is not locally owned, then there is nothing to do
928  if (std::get<0>(*cell)->is_locally_owned() == false)
929  *std::get<1>(*cell) = 0;
930  else
931  {
932  typename DerivativeDescription::Derivative derivative;
933  // call the function doing the actual
934  // work on this cell
935  approximate_cell<DerivativeDescription, dim, InputVector, spacedim>(
936  mapping,
937  dof_handler,
938  solution,
939  component,
940  std::get<0>(*cell),
941  derivative);
942 
943  // evaluate the norm and fill the vector
944  //*derivative_norm_on_this_cell
945  *std::get<1>(*cell) =
947  }
948  }
949 
950 
961  template <class DerivativeDescription,
962  int dim,
963  class InputVector,
964  int spacedim>
965  void
967  const DoFHandler<dim, spacedim> &dof_handler,
968  const InputVector & solution,
969  const unsigned int component,
971  {
972  Assert(derivative_norm.size() ==
973  dof_handler.get_triangulation().n_active_cells(),
975  derivative_norm.size(),
976  dof_handler.get_triangulation().n_active_cells()));
977  AssertIndexRange(component, dof_handler.get_fe(0).n_components());
978 
979  using Iterators =
980  std::tuple<typename DoFHandler<dim, spacedim>::active_cell_iterator,
983  Iterators(dof_handler.begin_active(), derivative_norm.begin())),
984  end(Iterators(dof_handler.end(), derivative_norm.end()));
985 
986  // There is no need for a copier because there is no conflict between
987  // threads to write in derivative_norm. Scratch and CopyData are also
988  // useless.
990  begin,
991  end,
992  [&mapping, &dof_handler, &solution, component](
994  Assembler::Scratch const &,
996  approximate<DerivativeDescription, dim, InputVector, spacedim>(
997  cell, mapping, dof_handler, solution, component);
998  },
999  std::function<void(internal::Assembler::CopyData const &)>(),
1002  }
1003 
1004  } // namespace internal
1005 
1006 } // namespace DerivativeApproximation
1007 
1008 
1009 // ------------------------ finally for the public interface of this namespace
1010 
1011 namespace DerivativeApproximation
1012 {
1013  template <int dim, class InputVector, int spacedim>
1014  void
1016  const DoFHandler<dim, spacedim> &dof_handler,
1017  const InputVector & solution,
1019  const unsigned int component)
1020  {
1021  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1022  mapping, dof_handler, solution, component, derivative_norm);
1023  }
1024 
1025 
1026  template <int dim, class InputVector, int spacedim>
1027  void
1029  const InputVector & solution,
1031  const unsigned int component)
1032  {
1033  Assert(!dof_handler.get_triangulation().is_mixed_mesh(),
1034  ExcNotImplemented());
1035  const auto reference_cell =
1036  dof_handler.get_triangulation().get_reference_cells()[0];
1037  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1038  reference_cell.template get_default_linear_mapping<dim, spacedim>(),
1039  dof_handler,
1040  solution,
1041  component,
1042  derivative_norm);
1043  }
1044 
1045 
1046  template <int dim, class InputVector, int spacedim>
1047  void
1049  const DoFHandler<dim, spacedim> &dof_handler,
1050  const InputVector & solution,
1052  const unsigned int component)
1053  {
1054  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1055  mapping, dof_handler, solution, component, derivative_norm);
1056  }
1057 
1058 
1059  template <int dim, class InputVector, int spacedim>
1060  void
1062  const InputVector & solution,
1064  const unsigned int component)
1065  {
1066  Assert(!dof_handler.get_triangulation().is_mixed_mesh(),
1067  ExcNotImplemented());
1068  const auto reference_cell =
1069  dof_handler.get_triangulation().get_reference_cells()[0];
1070  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1071  reference_cell.template get_default_linear_mapping<dim, spacedim>(),
1072  dof_handler,
1073  solution,
1074  component,
1075  derivative_norm);
1076  }
1077 
1078 
1079  template <int dim, int spacedim, class InputVector, int order>
1080  void
1082  const Mapping<dim, spacedim> & mapping,
1083  const DoFHandler<dim, spacedim> &dof,
1084  const InputVector & solution,
1085 #ifndef _MSC_VER
1087 #else
1089  &cell,
1090 #endif
1091  Tensor<order, dim> &derivative,
1092  const unsigned int component)
1093  {
1096  mapping, dof, solution, component, cell, derivative);
1097  }
1098 
1099 
1100 
1101  template <int dim, int spacedim, class InputVector, int order>
1102  void
1104  const DoFHandler<dim, spacedim> &dof,
1105  const InputVector & solution,
1106 #ifndef _MSC_VER
1108 #else
1110  &cell,
1111 #endif
1112  Tensor<order, dim> &derivative,
1113  const unsigned int component)
1114  {
1115  // just call the respective function with Q1 mapping
1117  cell->reference_cell()
1118  .template get_default_linear_mapping<dim, spacedim>(),
1119  dof,
1120  solution,
1121  cell,
1122  derivative,
1123  component);
1124  }
1125 
1126 
1127 
1128  template <int dim, int order>
1129  double
1131  {
1133  derivative_norm(derivative);
1134  }
1135 
1136 } // namespace DerivativeApproximation
1137 
1138 
1139 // --------------------------- explicit instantiations ---------------------
1140 #include "derivative_approximation.inst"
1141 
1142 
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 unsigned int index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3356
const Point< spacedim > & quadrature_point(const unsigned int q) const
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type >> &gradients) const
Definition: fe_values.cc:3495
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:3606
unsigned int n_components() const
Definition: point.h:111
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:692
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:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInsufficientDirections()
static ::ExceptionBase & ExcVectorLengthVsNActiveCells(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
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:1274
static constexpr double PI
Definition: numbers.h:248
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)