Reference documentation for deal.II version Git a39983f358 2020-01-20 20:47:10 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
derivative_approximation.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2019 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 
16 #include <deal.II/base/quadrature_lib.h>
17 #include <deal.II/base/work_stream.h>
18 
19 #include <deal.II/dofs/dof_accessor.h>
20 #include <deal.II/dofs/dof_handler.h>
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_values.h>
24 #include <deal.II/fe/mapping_q1.h>
25 
26 #include <deal.II/grid/filtered_iterator.h>
27 #include <deal.II/grid/grid_tools.h>
28 #include <deal.II/grid/tria_iterator.h>
29 
30 #include <deal.II/hp/fe_collection.h>
31 #include <deal.II/hp/fe_values.h>
32 #include <deal.II/hp/mapping_collection.h>
33 #include <deal.II/hp/q_collection.h>
34 
35 #include <deal.II/lac/block_vector.h>
36 #include <deal.II/lac/la_parallel_block_vector.h>
37 #include <deal.II/lac/la_parallel_vector.h>
38 #include <deal.II/lac/la_vector.h>
39 #include <deal.II/lac/petsc_block_vector.h>
40 #include <deal.II/lac/petsc_vector.h>
41 #include <deal.II/lac/trilinos_epetra_vector.h>
42 #include <deal.II/lac/trilinos_parallel_block_vector.h>
43 #include <deal.II/lac/trilinos_tpetra_vector.h>
44 #include <deal.II/lac/trilinos_vector.h>
45 #include <deal.II/lac/vector.h>
46 
47 #include <deal.II/numerics/derivative_approximation.h>
48 
49 #include <cmath>
50 
51 DEAL_II_NAMESPACE_OPEN
52 
53 
54 
55 namespace
56 {
57  template <typename T>
58  inline T
59  sqr(const T t)
60  {
61  return t * t;
62  }
63 } // namespace
64 
65 // --- First the classes and functions that describe individual derivatives ---
66 
68 {
69  namespace internal
70  {
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,
154  Vector<typename InputVector::value_type>(
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 
192  template <int dim>
194  {
195  public:
200  static const UpdateFlags update_flags;
201 
207 
213 
220  template <class InputVector, int spacedim>
221  static ProjectedDerivative
223  const InputVector & solution,
224  const unsigned int component);
225 
233  static double
234  derivative_norm(const Derivative &d);
235 
245  static void
246  symmetrize(Derivative &derivative_tensor);
247  };
248 
249  template <int dim>
251 
252 
253  template <int dim>
254  template <class InputVector, int spacedim>
257  const FEValues<dim, spacedim> &fe_values,
258  const InputVector & solution,
259  const unsigned int component)
260  {
261  if (fe_values.get_fe().n_components() == 1)
262  {
263  std::vector<Tensor<1, dim, typename InputVector::value_type>> values(
264  1);
265  fe_values.get_function_gradients(solution, values);
266  return ProjectedDerivative(values[0]);
267  }
268  else
269  {
270  std::vector<
271  std::vector<Tensor<1, dim, typename InputVector::value_type>>>
272  values(
273  1,
275  fe_values.get_fe().n_components()));
276  fe_values.get_function_gradients(solution, values);
277  return ProjectedDerivative(values[0][component]);
278  }
279  }
280 
281 
282 
283  template <>
284  inline double
286  {
287  return std::fabs(d[0][0]);
288  }
289 
290 
291 
292  template <>
293  inline double
295  {
296  // note that d should be a
297  // symmetric 2x2 tensor, so the
298  // eigenvalues are:
299  //
300  // 1/2(a+b\pm\sqrt((a-b)^2+4c^2))
301  //
302  // if the d_11=a, d_22=b,
303  // d_12=d_21=c
304  const double radicand =
305  ::sqr(d[0][0] - d[1][1]) + 4 * ::sqr(d[0][1]);
306  const double eigenvalues[2] = {
307  0.5 * (d[0][0] + d[1][1] + std::sqrt(radicand)),
308  0.5 * (d[0][0] + d[1][1] - std::sqrt(radicand))};
309 
310  return std::max(std::fabs(eigenvalues[0]), std::fabs(eigenvalues[1]));
311  }
312 
313 
314 
315  template <>
316  inline double
318  {
319  /*
320  compute the three eigenvalues of the tensor @p{d} and take the
321  largest. one could use the following maple script to generate C
322  code:
323 
324  with(linalg);
325  readlib(C);
326  A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]);
327  E:=eigenvals(A);
328  EE:=vector(3,[E[1],E[2],E[3]]);
329  C(EE);
330 
331  Unfortunately, with both optimized and non-optimized output, at some
332  places the code `sqrt(-1.0)' is emitted, and I don't know what
333  Maple intends to do with it. This happens both with Maple4 and
334  Maple5.
335 
336  Fortunately, Roger Young provided the following Fortran code, which
337  is transcribed below to C. The code uses an algorithm that uses the
338  invariants of a symmetric matrix. (The translated algorithm is
339  augmented by a test for R>0, since R==0 indicates that all three
340  eigenvalues are equal.)
341 
342 
343  PROGRAM MAIN
344 
345  C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX
346  C (ROGER YOUNG, 2001)
347 
348  IMPLICIT NONE
349 
350  REAL*8 A11, A12, A13, A22, A23, A33
351  REAL*8 I1, J2, J3, AM
352  REAL*8 S11, S12, S13, S22, S23, S33
353  REAL*8 SS12, SS23, SS13
354  REAL*8 R,R3, XX,YY, THETA
355  REAL*8 A1,A2,A3
356  REAL*8 PI
357  PARAMETER (PI=3.141592653587932384D0)
358  REAL*8 A,B,C, TOL
359  PARAMETER (TOL=1.D-14)
360 
361  C DEFINE A TEST MATRIX
362 
363  A11 = -1.D0
364  A12 = 5.D0
365  A13 = 3.D0
366  A22 = -2.D0
367  A23 = 0.5D0
368  A33 = 4.D0
369 
370 
371  I1 = A11 + A22 + A33
372  AM = I1/3.D0
373 
374  S11 = A11 - AM
375  S22 = A22 - AM
376  S33 = A33 - AM
377  S12 = A12
378  S13 = A13
379  S23 = A23
380 
381  SS12 = S12*S12
382  SS23 = S23*S23
383  SS13 = S13*S13
384 
385  J2 = S11*S11 + S22*S22 + S33*S33
386  J2 = J2 + 2.D0*(SS12 + SS23 + SS13)
387  J2 = J2/2.D0
388 
389  J3 = S11**3 + S22**3 + S33**3
390  J3 = J3 + 3.D0*S11*(SS12 + SS13)
391  J3 = J3 + 3.D0*S22*(SS12 + SS23)
392  J3 = J3 + 3.D0*S33*(SS13 + SS23)
393  J3 = J3 + 6.D0*S12*S23*S13
394  J3 = J3/3.D0
395 
396  R = SQRT(4.D0*J2/3.D0)
397  R3 = R*R*R
398  XX = 4.D0*J3/R3
399 
400  YY = 1.D0 - DABS(XX)
401  IF(YY.LE.0.D0)THEN
402  IF(YY.GT.(-TOL))THEN
403  WRITE(6,*)'Equal roots: XX= ',XX
404  A = -(XX/DABS(XX))*SQRT(J2/3.D0)
405  B = AM + A
406  C = AM - 2.D0*A
407  WRITE(6,*)B,' (twice) ',C
408  STOP
409  ELSE
410  WRITE(6,*)'Error: XX= ',XX
411  STOP
412  ENDIF
413  ENDIF
414 
415  THETA = (ACOS(XX))/3.D0
416 
417  A1 = AM + R*COS(THETA)
418  A2 = AM + R*COS(THETA + 2.D0*PI/3.D0)
419  A3 = AM + R*COS(THETA + 4.D0*PI/3.D0)
420 
421  WRITE(6,*)A1,A2,A3
422 
423  STOP
424  END
425 
426  */
427 
428  const double am = trace(d) / 3.;
429 
430  // s := d - trace(d) I
431  Tensor<2, 3> s = d;
432  for (unsigned int i = 0; i < 3; ++i)
433  s[i][i] -= am;
434 
435  const double ss01 = s[0][1] * s[0][1], ss12 = s[1][2] * s[1][2],
436  ss02 = s[0][2] * s[0][2];
437 
438  const double J2 = (s[0][0] * s[0][0] + s[1][1] * s[1][1] +
439  s[2][2] * s[2][2] + 2 * (ss01 + ss02 + ss12)) /
440  2.;
441  const double J3 =
442  (std::pow(s[0][0], 3) + std::pow(s[1][1], 3) + std::pow(s[2][2], 3) +
443  3. * s[0][0] * (ss01 + ss02) + 3. * s[1][1] * (ss01 + ss12) +
444  3. * s[2][2] * (ss02 + ss12) + 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>
529  class ThirdDerivative
530  {
531  public:
536  static const UpdateFlags update_flags;
537 
543  using Derivative = Tensor<3, dim>;
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>
587  const UpdateFlags ThirdDerivative<dim>::update_flags = update_hessians;
588 
589 
590  template <int dim>
591  template <class InputVector, int spacedim>
593  ThirdDerivative<dim>::get_projected_derivative(
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
622  ThirdDerivative<1>::derivative_norm(const Derivative &d)
623  {
624  return std::fabs(d[0][0][0]);
625  }
626 
627 
628 
629  template <int dim>
630  inline double
631  ThirdDerivative<dim>::derivative_norm(const Derivative &d)
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>
677  class DerivativeSelector
678  {
679  public:
685  using DerivDescr = void;
686  };
687 
688  template <int dim>
689  class DerivativeSelector<1, dim>
690  {
691  public:
692  using DerivDescr = Gradient<dim>;
693  };
694 
695  template <int dim>
696  class DerivativeSelector<2, dim>
697  {
698  public:
699  using DerivDescr = SecondDerivative<dim>;
700  };
701 
702  template <int dim>
703  class DerivativeSelector<3, dim>
704  {
705  public:
706  using DerivDescr = ThirdDerivative<dim>;
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  template <int, int> class DoFHandlerType,
744  class InputVector,
745  int spacedim>
746  void
747  approximate_cell(
748  const Mapping<dim, spacedim> & mapping,
749  const DoFHandlerType<dim, spacedim> &dof_handler,
750  const InputVector & solution,
751  const unsigned int component,
752  const TriaActiveIterator<
753  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>> &cell,
754  typename DerivativeDescription::Derivative &derivative)
755  {
756  QMidpoint<dim> midpoint_rule;
757 
758  // create collection objects from
759  // single quadratures, mappings,
760  // and finite elements. if we have
761  // an hp DoFHandler,
762  // dof_handler.get_fe() returns a
763  // collection of which we do a
764  // shallow copy instead
765  const hp::QCollection<dim> q_collection(midpoint_rule);
766  const hp::FECollection<dim> &fe_collection =
767  dof_handler.get_fe_collection();
768  const hp::MappingCollection<dim> mapping_collection(mapping);
769 
770  hp::FEValues<dim> x_fe_midpoint_value(
771  mapping_collection,
772  fe_collection,
773  q_collection,
774  DerivativeDescription::update_flags | update_quadrature_points);
775 
776  // matrix Y=sum_i y_i y_i^T
777  Tensor<2, dim> Y;
778 
779 
780  // vector to hold iterators to all
781  // active neighbors of a cell
782  // reserve the maximal number of
783  // active neighbors
784  std::vector<TriaActiveIterator<
786  active_neighbors;
787 
788  active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
790 
791  // vector
792  // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
793  // or related type for higher
794  // derivatives
795  typename DerivativeDescription::Derivative projected_derivative;
796 
797  // reinit fe values object...
798  x_fe_midpoint_value.reinit(cell);
799  const FEValues<dim> &fe_midpoint_value =
800  x_fe_midpoint_value.get_present_fe_values();
801 
802  // ...and get the value of the
803  // projected derivative...
804  const typename DerivativeDescription::ProjectedDerivative
805  this_midpoint_value =
806  DerivativeDescription::get_projected_derivative(fe_midpoint_value,
807  solution,
808  component);
809  // ...and the place where it lives
810  const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
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<DoFHandlerType<dim, spacedim>>(
825  cell, active_neighbors);
826 
827  // now loop over all active
828  // neighbors and collect the
829  // data we need
830  typename std::vector<TriaActiveIterator<
831  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>>>::
832  const_iterator neighbor_ptr = active_neighbors.begin();
833  for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
834  {
835  const TriaActiveIterator<
836  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>>
837  neighbor = *neighbor_ptr;
838 
839  // reinit fe values object...
840  x_fe_midpoint_value.reinit(neighbor);
841  const FEValues<dim> &neighbor_fe_midpoint_value =
842  x_fe_midpoint_value.get_present_fe_values();
843 
844  // ...and get the value of the
845  // solution...
846  const typename DerivativeDescription::ProjectedDerivative
847  neighbor_midpoint_value =
848  DerivativeDescription::get_projected_derivative(
849  neighbor_fe_midpoint_value, solution, component);
850 
851  // ...and the place where it lives
852  const Point<dim> neighbor_center =
853  neighbor_fe_midpoint_value.quadrature_point(0);
854 
855 
856  // vector for the
857  // normalized
858  // direction between
859  // the centers of two
860  // cells
861  Tensor<1, dim> y = neighbor_center - this_center;
862  const double distance = y.norm();
863  // normalize y
864  y /= distance;
865  // *** note that unlike in
866  // the docs, y denotes the
867  // normalized vector
868  // connecting the centers
869  // of the two cells, rather
870  // than the normal
871  // difference! ***
872 
873  // add up the
874  // contribution of
875  // this cell to Y
876  for (unsigned int i = 0; i < dim; ++i)
877  for (unsigned int j = 0; j < dim; ++j)
878  Y[i][j] += y[i] * y[j];
879 
880  // then update the sum
881  // of difference
882  // quotients
883  typename DerivativeDescription::ProjectedDerivative
884  projected_finite_difference =
885  (neighbor_midpoint_value - this_midpoint_value);
886  projected_finite_difference /= distance;
887 
888  projected_derivative += outer_product(y, projected_finite_difference);
889  }
890 
891  // can we determine an
892  // approximation of the
893  // gradient for the present
894  // cell? if so, then we need to
895  // have passed over vectors y_i
896  // which span the whole space,
897  // otherwise we would not have
898  // all components of the
899  // gradient
901 
902  // compute Y^-1 g
903  const Tensor<2, dim> Y_inverse = invert(Y);
904 
905  derivative = Y_inverse * projected_derivative;
906 
907  // finally symmetrize the derivative
909  }
910 
911 
912 
918  template <class DerivativeDescription,
919  int dim,
920  template <int, int> class DoFHandlerType,
921  class InputVector,
922  int spacedim>
923  void
924  approximate(
925  SynchronousIterators<std::tuple<
927  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>>,
928  Vector<float>::iterator>> const & cell,
929  const Mapping<dim, spacedim> & mapping,
930  const DoFHandlerType<dim, spacedim> &dof_handler,
931  const InputVector & solution,
932  const unsigned int component)
933  {
934  // if the cell is not locally owned, then there is nothing to do
935  if (std::get<0>(*cell)->is_locally_owned() == false)
936  *std::get<1>(*cell) = 0;
937  else
938  {
939  typename DerivativeDescription::Derivative derivative;
940  // call the function doing the actual
941  // work on this cell
942  approximate_cell<DerivativeDescription,
943  dim,
944  DoFHandlerType,
945  InputVector,
946  spacedim>(mapping,
947  dof_handler,
948  solution,
949  component,
950  std::get<0>(*cell),
951  derivative);
952 
953  // evaluate the norm and fill the vector
954  //*derivative_norm_on_this_cell
955  *std::get<1>(*cell) =
956  DerivativeDescription::derivative_norm(derivative);
957  }
958  }
959 
960 
971  template <class DerivativeDescription,
972  int dim,
973  template <int, int> class DoFHandlerType,
974  class InputVector,
975  int spacedim>
976  void
977  approximate_derivative(const Mapping<dim, spacedim> & mapping,
978  const DoFHandlerType<dim, spacedim> &dof_handler,
979  const InputVector & solution,
980  const unsigned int component,
981  Vector<float> & derivative_norm)
982  {
983  Assert(derivative_norm.size() ==
984  dof_handler.get_triangulation().n_active_cells(),
986  derivative_norm.size(),
987  dof_handler.get_triangulation().n_active_cells()));
988  AssertIndexRange(component, dof_handler.get_fe(0).n_components());
989 
990  using Iterators = std::tuple<
993  Vector<float>::iterator>;
995  Iterators(dof_handler.begin_active(), derivative_norm.begin())),
996  end(Iterators(dof_handler.end(), derivative_norm.end()));
997 
998  // There is no need for a copier because there is no conflict between
999  // threads to write in derivative_norm. Scratch and CopyData are also
1000  // useless.
1002  begin,
1003  end,
1004  [&mapping, &dof_handler, &solution, component](
1005  SynchronousIterators<Iterators> const &cell,
1006  Assembler::Scratch const &,
1007  Assembler::CopyData &) {
1008  approximate<DerivativeDescription,
1009  dim,
1010  DoFHandlerType,
1011  InputVector,
1012  spacedim>(
1013  cell, mapping, dof_handler, solution, component);
1014  },
1015  std::function<void(internal::Assembler::CopyData const &)>(),
1016  internal::Assembler::Scratch(),
1017  internal::Assembler::CopyData());
1018  }
1019 
1020  } // namespace internal
1021 
1022 } // namespace DerivativeApproximation
1023 
1024 
1025 // ------------------------ finally for the public interface of this namespace
1026 
1027 namespace DerivativeApproximation
1028 {
1029  template <int dim,
1030  template <int, int> class DoFHandlerType,
1031  class InputVector,
1032  int spacedim>
1033  void
1035  const DoFHandlerType<dim, spacedim> &dof_handler,
1036  const InputVector & solution,
1037  Vector<float> & derivative_norm,
1038  const unsigned int component)
1039  {
1040  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1041  mapping, dof_handler, solution, component, derivative_norm);
1042  }
1043 
1044 
1045  template <int dim,
1046  template <int, int> class DoFHandlerType,
1047  class InputVector,
1048  int spacedim>
1049  void
1050  approximate_gradient(const DoFHandlerType<dim, spacedim> &dof_handler,
1051  const InputVector & solution,
1052  Vector<float> & derivative_norm,
1053  const unsigned int component)
1054  {
1055  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1057  dof_handler,
1058  solution,
1059  component,
1060  derivative_norm);
1061  }
1062 
1063 
1064  template <int dim,
1065  template <int, int> class DoFHandlerType,
1066  class InputVector,
1067  int spacedim>
1068  void
1070  const Mapping<dim, spacedim> & mapping,
1071  const DoFHandlerType<dim, spacedim> &dof_handler,
1072  const InputVector & solution,
1073  Vector<float> & derivative_norm,
1074  const unsigned int component)
1075  {
1076  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1077  mapping, dof_handler, solution, component, derivative_norm);
1078  }
1079 
1080 
1081  template <int dim,
1082  template <int, int> class DoFHandlerType,
1083  class InputVector,
1084  int spacedim>
1085  void
1087  const DoFHandlerType<dim, spacedim> &dof_handler,
1088  const InputVector & solution,
1089  Vector<float> & derivative_norm,
1090  const unsigned int component)
1091  {
1092  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1094  dof_handler,
1095  solution,
1096  component,
1097  derivative_norm);
1098  }
1099 
1100 
1101  template <typename DoFHandlerType, class InputVector, int order>
1102  void
1105  & mapping,
1106  const DoFHandlerType &dof,
1107  const InputVector & solution,
1108 #ifndef _MSC_VER
1109  const typename DoFHandlerType::active_cell_iterator &cell,
1110 #else
1112  &cell,
1113 #endif
1115  const unsigned int component)
1116  {
1117  internal::approximate_cell<
1118  typename internal::DerivativeSelector<order, DoFHandlerType::dimension>::
1119  DerivDescr>(mapping, dof, solution, component, cell, derivative);
1120  }
1121 
1122 
1123 
1124  template <typename DoFHandlerType, class InputVector, int order>
1125  void
1127  const DoFHandlerType &dof,
1128  const InputVector & solution,
1129 #ifndef _MSC_VER
1130  const typename DoFHandlerType::active_cell_iterator &cell,
1131 #else
1133  &cell,
1134 #endif
1136  const unsigned int component)
1137  {
1138  // just call the respective function with Q1 mapping
1140  StaticMappingQ1<DoFHandlerType::dimension,
1141  DoFHandlerType::space_dimension>::mapping,
1142  dof,
1143  solution,
1144  cell,
1145  derivative,
1146  component);
1147  }
1148 
1149 
1150 
1151  template <int dim, int order>
1152  double
1154  {
1155  return internal::DerivativeSelector<order, dim>::DerivDescr::
1156  derivative_norm(derivative);
1157  }
1158 
1159 } // namespace DerivativeApproximation
1160 
1161 
1162 // --------------------------- explicit instantiations ---------------------
1163 #include "derivative_approximation.inst"
1164 
1165 
1166 DEAL_II_NAMESPACE_CLOSE
Shape function values.
static ::ExceptionBase & ExcInsufficientDirections()
static void symmetrize(Derivative &derivative_tensor)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3572
const FiniteElement< dim, spacedim > & get_fe() const
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1429
static ProjectedDerivative get_projected_derivative(const FEValues< dim, spacedim > &fe_values, const InputVector &solution, const unsigned int component)
const Point< spacedim > & quadrature_point(const unsigned int q) const
#define Assert(cond, exc)
Definition: exceptions.h:1411
UpdateFlags
Abstract base class for mapping classes.
Definition: mapping.h:302
DEAL_II_CONSTEXPR Number trace(const SymmetricTensor< 2, dim, Number > &d)
static void symmetrize(Derivative &derivative_tensor)
DEAL_II_CONSTEXPR SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static ProjectedDerivative get_projected_derivative(const FEValues< dim, spacedim > &fe_values, const InputVector &solution, const unsigned int component)
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType, 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:132
Second derivatives of shape functions.
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type >> &hessians) const
Definition: fe_values.cc:3826
void approximate_derivative_tensor(const Mapping< DoFHandlerType::dimension, DoFHandlerType::space_dimension > &mapping, const DoFHandlerType &dof, const InputVector &solution, const typename DoFHandlerType::active_cell_iterator &cell, Tensor< order, DoFHandlerType::dimension > &derivative, const unsigned int component=0)
static constexpr double PI
Definition: numbers.h:238
Shape function gradients.
DEAL_II_CONSTEXPR Number determinant(const SymmetricTensor< 2, dim, Number > &)
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type >> &gradients) const
Definition: fe_values.cc:3713
Definition: fe.h:38
static ::ExceptionBase & ExcNotImplemented()
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:1187
static double derivative_norm(const Derivative &d)
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
void approximate_second_derivative(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const InputVector &solution, Vector< float > &derivative_norm, const unsigned int component=0)
void approximate_gradient(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const InputVector &solution, Vector< float > &derivative_norm, const unsigned int component=0)
static ::ExceptionBase & ExcInternalError()
const ::FEValues< dim, spacedim > & get_present_fe_values() const
Definition: fe_values.h:569
static ::ExceptionBase & ExcVectorLengthVsNActiveCells(int arg1, int arg2)