Reference documentation for deal.II version Git f346ebc02b 2020-06-06 16:01:17 +0200
\(\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 - 2020 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_q1.h>
25 
29 
31 #include <deal.II/hp/fe_values.h>
34 
38 #include <deal.II/lac/la_vector.h>
45 #include <deal.II/lac/vector.h>
46 
48 
49 #include <cmath>
50 
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  {
77  template <int dim>
78  class Gradient
79  {
80  public:
85  static const UpdateFlags update_flags;
86 
92 
98 
105  template <class InputVector, int spacedim>
106  static ProjectedDerivative
108  const InputVector & solution,
109  const unsigned int component);
110 
115  static double
116  derivative_norm(const Derivative &d);
117 
125  static void
126  symmetrize(Derivative &derivative_tensor);
127  };
128 
129  // static variables
130  template <int dim>
132 
133 
134  template <int dim>
135  template <class InputVector, int spacedim>
136  inline typename Gradient<dim>::ProjectedDerivative
138  const FEValues<dim, spacedim> &fe_values,
139  const InputVector & solution,
140  const unsigned int component)
141  {
142  if (fe_values.get_fe().n_components() == 1)
143  {
144  std::vector<typename InputVector::value_type> values(1);
145  fe_values.get_function_values(solution, values);
146  return values[0];
147  }
148  else
149  {
150  std::vector<Vector<typename InputVector::value_type>> values(
151  1,
152  Vector<typename InputVector::value_type>(
153  fe_values.get_fe().n_components()));
154  fe_values.get_function_values(solution, values);
155  return values[0](component);
156  }
157  }
158 
159 
160 
161  template <int dim>
162  inline double
164  {
165  double s = 0;
166  for (unsigned int i = 0; i < dim; ++i)
167  s += d[i] * d[i];
168  return std::sqrt(s);
169  }
170 
171 
172 
173  template <int dim>
174  inline void
176  {
177  // nothing to do here
178  }
179 
180 
181 
188  template <int dim>
190  {
191  public:
196  static const UpdateFlags update_flags;
197 
203 
209 
216  template <class InputVector, int spacedim>
217  static ProjectedDerivative
219  const InputVector & solution,
220  const unsigned int component);
221 
229  static double
230  derivative_norm(const Derivative &d);
231 
241  static void
242  symmetrize(Derivative &derivative_tensor);
243  };
244 
245  template <int dim>
247 
248 
249  template <int dim>
250  template <class InputVector, int spacedim>
253  const FEValues<dim, spacedim> &fe_values,
254  const InputVector & solution,
255  const unsigned int component)
256  {
257  if (fe_values.get_fe().n_components() == 1)
258  {
259  std::vector<Tensor<1, dim, typename InputVector::value_type>> values(
260  1);
261  fe_values.get_function_gradients(solution, values);
262  return ProjectedDerivative(values[0]);
263  }
264  else
265  {
266  std::vector<
267  std::vector<Tensor<1, dim, typename InputVector::value_type>>>
268  values(
269  1,
271  fe_values.get_fe().n_components()));
272  fe_values.get_function_gradients(solution, values);
273  return ProjectedDerivative(values[0][component]);
274  }
275  }
276 
277 
278 
279  template <>
280  inline double
282  {
283  return std::fabs(d[0][0]);
284  }
285 
286 
287 
288  template <>
289  inline double
291  {
292  // note that d should be a
293  // symmetric 2x2 tensor, so the
294  // eigenvalues are:
295  //
296  // 1/2(a+b\pm\sqrt((a-b)^2+4c^2))
297  //
298  // if the d_11=a, d_22=b,
299  // d_12=d_21=c
300  const double radicand =
301  ::sqr(d[0][0] - d[1][1]) + 4 * ::sqr(d[0][1]);
302  const double eigenvalues[2] = {
303  0.5 * (d[0][0] + d[1][1] + std::sqrt(radicand)),
304  0.5 * (d[0][0] + d[1][1] - std::sqrt(radicand))};
305 
306  return std::max(std::fabs(eigenvalues[0]), std::fabs(eigenvalues[1]));
307  }
308 
309 
310 
311  template <>
312  inline double
314  {
315  /*
316  compute the three eigenvalues of the tensor @p{d} and take the
317  largest. one could use the following maple script to generate C
318  code:
319 
320  with(linalg);
321  readlib(C);
322  A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]);
323  E:=eigenvals(A);
324  EE:=vector(3,[E[1],E[2],E[3]]);
325  C(EE);
326 
327  Unfortunately, with both optimized and non-optimized output, at some
328  places the code `sqrt(-1.0)' is emitted, and I don't know what
329  Maple intends to do with it. This happens both with Maple4 and
330  Maple5.
331 
332  Fortunately, Roger Young provided the following Fortran code, which
333  is transcribed below to C. The code uses an algorithm that uses the
334  invariants of a symmetric matrix. (The translated algorithm is
335  augmented by a test for R>0, since R==0 indicates that all three
336  eigenvalues are equal.)
337 
338 
339  PROGRAM MAIN
340 
341  C FIND EIGENVALUES OF REAL SYMMETRIC MATRIX
342  C (ROGER YOUNG, 2001)
343 
344  IMPLICIT NONE
345 
346  REAL*8 A11, A12, A13, A22, A23, A33
347  REAL*8 I1, J2, J3, AM
348  REAL*8 S11, S12, S13, S22, S23, S33
349  REAL*8 SS12, SS23, SS13
350  REAL*8 R,R3, XX,YY, THETA
351  REAL*8 A1,A2,A3
352  REAL*8 PI
353  PARAMETER (PI=3.141592653587932384D0)
354  REAL*8 A,B,C, TOL
355  PARAMETER (TOL=1.D-14)
356 
357  C DEFINE A TEST MATRIX
358 
359  A11 = -1.D0
360  A12 = 5.D0
361  A13 = 3.D0
362  A22 = -2.D0
363  A23 = 0.5D0
364  A33 = 4.D0
365 
366 
367  I1 = A11 + A22 + A33
368  AM = I1/3.D0
369 
370  S11 = A11 - AM
371  S22 = A22 - AM
372  S33 = A33 - AM
373  S12 = A12
374  S13 = A13
375  S23 = A23
376 
377  SS12 = S12*S12
378  SS23 = S23*S23
379  SS13 = S13*S13
380 
381  J2 = S11*S11 + S22*S22 + S33*S33
382  J2 = J2 + 2.D0*(SS12 + SS23 + SS13)
383  J2 = J2/2.D0
384 
385  J3 = S11**3 + S22**3 + S33**3
386  J3 = J3 + 3.D0*S11*(SS12 + SS13)
387  J3 = J3 + 3.D0*S22*(SS12 + SS23)
388  J3 = J3 + 3.D0*S33*(SS13 + SS23)
389  J3 = J3 + 6.D0*S12*S23*S13
390  J3 = J3/3.D0
391 
392  R = SQRT(4.D0*J2/3.D0)
393  R3 = R*R*R
394  XX = 4.D0*J3/R3
395 
396  YY = 1.D0 - DABS(XX)
397  IF(YY.LE.0.D0)THEN
398  IF(YY.GT.(-TOL))THEN
399  WRITE(6,*)'Equal roots: XX= ',XX
400  A = -(XX/DABS(XX))*SQRT(J2/3.D0)
401  B = AM + A
402  C = AM - 2.D0*A
403  WRITE(6,*)B,' (twice) ',C
404  STOP
405  ELSE
406  WRITE(6,*)'Error: XX= ',XX
407  STOP
408  ENDIF
409  ENDIF
410 
411  THETA = (ACOS(XX))/3.D0
412 
413  A1 = AM + R*COS(THETA)
414  A2 = AM + R*COS(THETA + 2.D0*PI/3.D0)
415  A3 = AM + R*COS(THETA + 4.D0*PI/3.D0)
416 
417  WRITE(6,*)A1,A2,A3
418 
419  STOP
420  END
421 
422  */
423 
424  const double am = trace(d) / 3.;
425 
426  // s := d - trace(d) I
427  Tensor<2, 3> s = d;
428  for (unsigned int i = 0; i < 3; ++i)
429  s[i][i] -= am;
430 
431  const double ss01 = s[0][1] * s[0][1], ss12 = s[1][2] * s[1][2],
432  ss02 = s[0][2] * s[0][2];
433 
434  const double J2 = (s[0][0] * s[0][0] + s[1][1] * s[1][1] +
435  s[2][2] * s[2][2] + 2 * (ss01 + ss02 + ss12)) /
436  2.;
437  const double J3 =
438  (std::pow(s[0][0], 3) + std::pow(s[1][1], 3) + std::pow(s[2][2], 3) +
439  3. * s[0][0] * (ss01 + ss02) + 3. * s[1][1] * (ss01 + ss12) +
440  3. * s[2][2] * (ss02 + ss12) + 6. * s[0][1] * s[0][2] * s[1][2]) /
441  3.;
442 
443  const double R = std::sqrt(4. * J2 / 3.);
444 
445  double EE[3] = {0, 0, 0};
446  // the eigenvalues are away from
447  // @p{am} in the order of R. thus,
448  // if R<<AM, then we have the
449  // degenerate case with three
450  // identical eigenvalues. check
451  // this first
452  if (R <= 1e-14 * std::fabs(am))
453  EE[0] = EE[1] = EE[2] = am;
454  else
455  {
456  // at least two eigenvalues are
457  // distinct
458  const double R3 = R * R * R;
459  const double XX = 4. * J3 / R3;
460  const double YY = 1. - std::fabs(XX);
461 
462  Assert(YY > -1e-14, ExcInternalError());
463 
464  if (YY < 0)
465  {
466  // two roots are equal
467  const double a = (XX > 0 ? -1. : 1.) * R / 2;
468  EE[0] = EE[1] = am + a;
469  EE[2] = am - 2. * a;
470  }
471  else
472  {
473  const double theta = std::acos(XX) / 3.;
474  EE[0] = am + R * std::cos(theta);
475  EE[1] = am + R * std::cos(theta + 2. / 3. * numbers::PI);
476  EE[2] = am + R * std::cos(theta + 4. / 3. * numbers::PI);
477  }
478  }
479 
480  return std::max(std::fabs(EE[0]),
481  std::max(std::fabs(EE[1]), std::fabs(EE[2])));
482  }
483 
484 
485 
486  template <int dim>
487  inline double
489  {
490  // computing the spectral norm is
491  // not so simple in general. it is
492  // feasible for dim==3 as shown
493  // above, since then there are
494  // still closed form expressions of
495  // the roots of the characteristic
496  // polynomial, and they can easily
497  // be computed using
498  // maple. however, for higher
499  // dimensions, some other method
500  // needs to be employed. maybe some
501  // steps of the power method would
502  // suffice?
503  Assert(false, ExcNotImplemented());
504  return 0;
505  }
506 
507 
508 
509  template <int dim>
510  inline void
512  {
513  // symmetrize non-diagonal entries
514  for (unsigned int i = 0; i < dim; ++i)
515  for (unsigned int j = i + 1; j < dim; ++j)
516  {
517  const double s = (d[i][j] + d[j][i]) / 2;
518  d[i][j] = d[j][i] = s;
519  }
520  }
521 
522 
523 
524  template <int dim>
526  {
527  public:
532  static const UpdateFlags update_flags;
533 
540 
546 
553  template <class InputVector, int spacedim>
554  static ProjectedDerivative
556  const InputVector & solution,
557  const unsigned int component);
558 
566  static double
567  derivative_norm(const Derivative &d);
568 
578  static void
579  symmetrize(Derivative &derivative_tensor);
580  };
581 
582  template <int dim>
584 
585 
586  template <int dim>
587  template <class InputVector, int spacedim>
590  const FEValues<dim, spacedim> &fe_values,
591  const InputVector & solution,
592  const unsigned int component)
593  {
594  if (fe_values.get_fe().n_components() == 1)
595  {
596  std::vector<Tensor<2, dim, typename InputVector::value_type>> values(
597  1);
598  fe_values.get_function_hessians(solution, values);
599  return ProjectedDerivative(values[0]);
600  }
601  else
602  {
603  std::vector<
604  std::vector<Tensor<2, dim, typename InputVector::value_type>>>
605  values(
606  1,
608  fe_values.get_fe().n_components()));
609  fe_values.get_function_hessians(solution, values);
610  return ProjectedDerivative(values[0][component]);
611  }
612  }
613 
614 
615 
616  template <>
617  inline double
619  {
620  return std::fabs(d[0][0][0]);
621  }
622 
623 
624 
625  template <int dim>
626  inline double
628  {
629  // return the Frobenius-norm. this is a
630  // member function of Tensor<rank_,dim>
631  return d.norm();
632  }
633 
634 
635  template <int dim>
636  inline void
638  {
639  // symmetrize non-diagonal entries
640 
641  // first do it in the case, that i,j,k are
642  // pairwise different (which can onlky happen
643  // in dim >= 3)
644  for (unsigned int i = 0; i < dim; ++i)
645  for (unsigned int j = i + 1; j < dim; ++j)
646  for (unsigned int k = j + 1; k < dim; ++k)
647  {
648  const double s = (d[i][j][k] + d[i][k][j] + d[j][i][k] +
649  d[j][k][i] + d[k][i][j] + d[k][j][i]) /
650  6;
651  d[i][j][k] = d[i][k][j] = d[j][i][k] = d[j][k][i] = d[k][i][j] =
652  d[k][j][i] = s;
653  }
654  // now do the case, where two indices are
655  // equal
656  for (unsigned int i = 0; i < dim; ++i)
657  for (unsigned int j = i + 1; j < dim; ++j)
658  {
659  // case 1: index i (lower one) is
660  // double
661  const double s = (d[i][i][j] + d[i][j][i] + d[j][i][i]) / 3;
662  d[i][i][j] = d[i][j][i] = d[j][i][i] = s;
663 
664  // case 2: index j (higher one) is
665  // double
666  const double t = (d[i][j][j] + d[j][i][j] + d[j][j][i]) / 3;
667  d[i][j][j] = d[j][i][j] = d[j][j][i] = t;
668  }
669  }
670 
671 
672  template <int order, int dim>
674  {
675  public:
681  using DerivDescr = void;
682  };
683 
684  template <int dim>
685  class DerivativeSelector<1, dim>
686  {
687  public:
689  };
690 
691  template <int dim>
692  class DerivativeSelector<2, dim>
693  {
694  public:
696  };
697 
698  template <int dim>
699  class DerivativeSelector<3, dim>
700  {
701  public:
703  };
704  } // namespace internal
705 } // namespace DerivativeApproximation
706 
707 // Dummy structures and dummy function used for WorkStream
708 namespace DerivativeApproximation
709 {
710  namespace internal
711  {
712  namespace Assembler
713  {
714  struct Scratch
715  {
716  Scratch() = default;
717  };
718 
719  struct CopyData
720  {
721  CopyData() = default;
722  };
723  } // namespace Assembler
724  } // namespace internal
725 } // namespace DerivativeApproximation
726 
727 // --------------- now for the functions that do the actual work --------------
728 
729 namespace DerivativeApproximation
730 {
731  namespace internal
732  {
737  template <class DerivativeDescription,
738  int dim,
739  template <int, int> class DoFHandlerType,
740  class InputVector,
741  int spacedim>
742  void
744  const Mapping<dim, spacedim> & mapping,
745  const DoFHandlerType<dim, spacedim> &dof_handler,
746  const InputVector & solution,
747  const unsigned int component,
748  const TriaActiveIterator<
749  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>> &cell,
750  typename DerivativeDescription::Derivative &derivative)
751  {
752  QMidpoint<dim> midpoint_rule;
753 
754  // create collection objects from
755  // single quadratures, mappings,
756  // and finite elements. if we have
757  // an hp DoFHandler,
758  // dof_handler.get_fe() returns a
759  // collection of which we do a
760  // shallow copy instead
761  const hp::QCollection<dim> q_collection(midpoint_rule);
762  const hp::FECollection<dim> &fe_collection =
763  dof_handler.get_fe_collection();
764  const hp::MappingCollection<dim> mapping_collection(mapping);
765 
766  hp::FEValues<dim> x_fe_midpoint_value(
767  mapping_collection,
768  fe_collection,
769  q_collection,
770  DerivativeDescription::update_flags | update_quadrature_points);
771 
772  // matrix Y=sum_i y_i y_i^T
773  Tensor<2, dim> Y;
774 
775 
776  // vector to hold iterators to all
777  // active neighbors of a cell
778  // reserve the maximal number of
779  // active neighbors
780  std::vector<TriaActiveIterator<
782  active_neighbors;
783 
784  active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
786 
787  // vector
788  // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
789  // or related type for higher
790  // derivatives
791  typename DerivativeDescription::Derivative projected_derivative;
792 
793  // reinit fe values object...
794  x_fe_midpoint_value.reinit(cell);
795  const FEValues<dim> &fe_midpoint_value =
796  x_fe_midpoint_value.get_present_fe_values();
797 
798  // ...and get the value of the
799  // projected derivative...
800  const typename DerivativeDescription::ProjectedDerivative
801  this_midpoint_value =
802  DerivativeDescription::get_projected_derivative(fe_midpoint_value,
803  solution,
804  component);
805  // ...and the place where it lives
806  const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
807 
808  // loop over all neighbors and
809  // accumulate the difference
810  // quotients from them. note
811  // that things get a bit more
812  // complicated if the neighbor
813  // is more refined than the
814  // present one
815  //
816  // to make processing simpler,
817  // first collect all neighbor
818  // cells in a vector, and then
819  // collect the data from them
820  GridTools::get_active_neighbors<DoFHandlerType<dim, spacedim>>(
821  cell, active_neighbors);
822 
823  // now loop over all active
824  // neighbors and collect the
825  // data we need
826  typename std::vector<TriaActiveIterator<
827  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>>>::
828  const_iterator neighbor_ptr = active_neighbors.begin();
829  for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
830  {
831  const TriaActiveIterator<
832  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>>
833  neighbor = *neighbor_ptr;
834 
835  // reinit fe values object...
836  x_fe_midpoint_value.reinit(neighbor);
837  const FEValues<dim> &neighbor_fe_midpoint_value =
838  x_fe_midpoint_value.get_present_fe_values();
839 
840  // ...and get the value of the
841  // solution...
842  const typename DerivativeDescription::ProjectedDerivative
843  neighbor_midpoint_value =
844  DerivativeDescription::get_projected_derivative(
845  neighbor_fe_midpoint_value, solution, component);
846 
847  // ...and the place where it lives
848  const Point<dim> neighbor_center =
849  neighbor_fe_midpoint_value.quadrature_point(0);
850 
851 
852  // vector for the
853  // normalized
854  // direction between
855  // the centers of two
856  // cells
857  Tensor<1, dim> y = neighbor_center - this_center;
858  const double distance = y.norm();
859  // normalize y
860  y /= distance;
861  // *** note that unlike in
862  // the docs, y denotes the
863  // normalized vector
864  // connecting the centers
865  // of the two cells, rather
866  // than the normal
867  // difference! ***
868 
869  // add up the
870  // contribution of
871  // this cell to Y
872  for (unsigned int i = 0; i < dim; ++i)
873  for (unsigned int j = 0; j < dim; ++j)
874  Y[i][j] += y[i] * y[j];
875 
876  // then update the sum
877  // of difference
878  // quotients
879  typename DerivativeDescription::ProjectedDerivative
880  projected_finite_difference =
881  (neighbor_midpoint_value - this_midpoint_value);
882  projected_finite_difference /= distance;
883 
884  projected_derivative += outer_product(y, projected_finite_difference);
885  }
886 
887  // can we determine an
888  // approximation of the
889  // gradient for the present
890  // cell? if so, then we need to
891  // have passed over vectors y_i
892  // which span the whole space,
893  // otherwise we would not have
894  // all components of the
895  // gradient
897 
898  // compute Y^-1 g
899  const Tensor<2, dim> Y_inverse = invert(Y);
900 
901  derivative = Y_inverse * projected_derivative;
902 
903  // finally symmetrize the derivative
905  }
906 
907 
908 
914  template <class DerivativeDescription,
915  int dim,
916  template <int, int> class DoFHandlerType,
917  class InputVector,
918  int spacedim>
919  void
921  SynchronousIterators<std::tuple<
923  ::DoFCellAccessor<DoFHandlerType<dim, spacedim>, false>>,
924  Vector<float>::iterator>> const & cell,
925  const Mapping<dim, spacedim> & mapping,
926  const DoFHandlerType<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,
939  dim,
940  DoFHandlerType,
941  InputVector,
942  spacedim>(mapping,
943  dof_handler,
944  solution,
945  component,
946  std::get<0>(*cell),
947  derivative);
948 
949  // evaluate the norm and fill the vector
950  //*derivative_norm_on_this_cell
951  *std::get<1>(*cell) =
953  }
954  }
955 
956 
967  template <class DerivativeDescription,
968  int dim,
969  template <int, int> class DoFHandlerType,
970  class InputVector,
971  int spacedim>
972  void
974  const DoFHandlerType<dim, spacedim> &dof_handler,
975  const InputVector & solution,
976  const unsigned int component,
977  Vector<float> & derivative_norm)
978  {
979  Assert(derivative_norm.size() ==
980  dof_handler.get_triangulation().n_active_cells(),
982  derivative_norm.size(),
983  dof_handler.get_triangulation().n_active_cells()));
984  AssertIndexRange(component, dof_handler.get_fe(0).n_components());
985 
986  using Iterators = std::tuple<
991  Iterators(dof_handler.begin_active(), derivative_norm.begin())),
992  end(Iterators(dof_handler.end(), derivative_norm.end()));
993 
994  // There is no need for a copier because there is no conflict between
995  // threads to write in derivative_norm. Scratch and CopyData are also
996  // useless.
998  begin,
999  end,
1000  [&mapping, &dof_handler, &solution, component](
1001  SynchronousIterators<Iterators> const &cell,
1002  Assembler::Scratch const &,
1003  Assembler::CopyData &) {
1004  approximate<DerivativeDescription,
1005  dim,
1006  DoFHandlerType,
1007  InputVector,
1008  spacedim>(
1009  cell, mapping, dof_handler, solution, component);
1010  },
1011  std::function<void(internal::Assembler::CopyData const &)>(),
1014  }
1015 
1016  } // namespace internal
1017 
1018 } // namespace DerivativeApproximation
1019 
1020 
1021 // ------------------------ finally for the public interface of this namespace
1022 
1023 namespace DerivativeApproximation
1024 {
1025  template <int dim,
1026  template <int, int> class DoFHandlerType,
1027  class InputVector,
1028  int spacedim>
1029  void
1031  const DoFHandlerType<dim, spacedim> &dof_handler,
1032  const InputVector & solution,
1033  Vector<float> & derivative_norm,
1034  const unsigned int component)
1035  {
1036  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1037  mapping, dof_handler, solution, component, derivative_norm);
1038  }
1039 
1040 
1041  template <int dim,
1042  template <int, int> class DoFHandlerType,
1043  class InputVector,
1044  int spacedim>
1045  void
1046  approximate_gradient(const DoFHandlerType<dim, spacedim> &dof_handler,
1047  const InputVector & solution,
1048  Vector<float> & derivative_norm,
1049  const unsigned int component)
1050  {
1051  internal::approximate_derivative<internal::Gradient<dim>, dim>(
1053  dof_handler,
1054  solution,
1055  component,
1056  derivative_norm);
1057  }
1058 
1059 
1060  template <int dim,
1061  template <int, int> class DoFHandlerType,
1062  class InputVector,
1063  int spacedim>
1064  void
1066  const Mapping<dim, spacedim> & mapping,
1067  const DoFHandlerType<dim, spacedim> &dof_handler,
1068  const InputVector & solution,
1069  Vector<float> & derivative_norm,
1070  const unsigned int component)
1071  {
1072  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1073  mapping, dof_handler, solution, component, derivative_norm);
1074  }
1075 
1076 
1077  template <int dim,
1078  template <int, int> class DoFHandlerType,
1079  class InputVector,
1080  int spacedim>
1081  void
1083  const DoFHandlerType<dim, spacedim> &dof_handler,
1084  const InputVector & solution,
1085  Vector<float> & derivative_norm,
1086  const unsigned int component)
1087  {
1088  internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1090  dof_handler,
1091  solution,
1092  component,
1093  derivative_norm);
1094  }
1095 
1096 
1097  template <typename DoFHandlerType, class InputVector, int order>
1098  void
1101  & mapping,
1102  const DoFHandlerType &dof,
1103  const InputVector & solution,
1104 #ifndef _MSC_VER
1105  const typename DoFHandlerType::active_cell_iterator &cell,
1106 #else
1108  &cell,
1109 #endif
1111  const unsigned int component)
1112  {
1115  DerivDescr>(mapping, dof, solution, component, cell, derivative);
1116  }
1117 
1118 
1119 
1120  template <typename DoFHandlerType, class InputVector, int order>
1121  void
1123  const DoFHandlerType &dof,
1124  const InputVector & solution,
1125 #ifndef _MSC_VER
1126  const typename DoFHandlerType::active_cell_iterator &cell,
1127 #else
1129  &cell,
1130 #endif
1132  const unsigned int component)
1133  {
1134  // just call the respective function with Q1 mapping
1136  StaticMappingQ1<DoFHandlerType::dimension,
1137  DoFHandlerType::space_dimension>::mapping,
1138  dof,
1139  solution,
1140  cell,
1141  derivative,
1142  component);
1143  }
1144 
1145 
1146 
1147  template <int dim, int order>
1148  double
1150  {
1152  derivative_norm(derivative);
1153  }
1154 
1155 } // namespace DerivativeApproximation
1156 
1157 
1158 // --------------------------- explicit instantiations ---------------------
1159 #include "derivative_approximation.inst"
1160 
1161 
Shape function values.
static ::ExceptionBase & ExcInsufficientDirections()
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
static void symmetrize(Derivative &derivative_tensor)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void approximate_cell(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component, const TriaActiveIterator< ::DoFCellAccessor< DoFHandlerType< dim, spacedim >, false >> &cell, typename DerivativeDescription::Derivative &derivative)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
static ProjectedDerivative get_projected_derivative(const FEValues< dim, spacedim > &fe_values, const InputVector &solution, const unsigned int component)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
const Point< spacedim > & quadrature_point(const unsigned int q) const
static void symmetrize(Derivative &derivative_tensor)
static const char T
Expression acos(const Expression &x)
#define Assert(cond, exc)
Definition: exceptions.h:1403
UpdateFlags
Abstract base class for mapping classes.
Definition: mapping.h:301
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:361
static void symmetrize(Derivative &derivative_tensor)
VectorType::value_type * end(VectorType &V)
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:228
Expression fabs(const Expression &x)
void approximate_derivative(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component, Vector< float > &derivative_norm)
Second derivatives of shape functions.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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:231
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
VectorType::value_type * begin(VectorType &V)
Shape function gradients.
void approximate(SynchronousIterators< std::tuple< TriaActiveIterator< ::DoFCellAccessor< DoFHandlerType< dim, spacedim >, false >>, Vector< float >::iterator >> const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
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:1183
static ProjectedDerivative get_projected_derivative(const FEValues< dim, spacedim > &fe_values, const InputVector &solution, const unsigned int component)
static double derivative_norm(const Derivative &d)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
numbers::NumberTraits< Number >::real_type norm() const
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)
Eigenvalue vector is filled.
T max(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcInternalError()
const ::FEValues< dim, spacedim > & get_present_fe_values() const
Definition: fe_values.h:610
static ::ExceptionBase & ExcVectorLengthVsNActiveCells(int arg1, int arg2)