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