16 #include <deal.II/base/quadrature_lib.h> 17 #include <deal.II/base/work_stream.h> 19 #include <deal.II/dofs/dof_accessor.h> 20 #include <deal.II/dofs/dof_handler.h> 22 #include <deal.II/fe/fe.h> 23 #include <deal.II/fe/fe_values.h> 24 #include <deal.II/fe/mapping_q1.h> 26 #include <deal.II/grid/filtered_iterator.h> 27 #include <deal.II/grid/grid_tools.h> 28 #include <deal.II/grid/tria_iterator.h> 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> 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_vector.h> 43 #include <deal.II/lac/vector.h> 45 #include <deal.II/numerics/derivative_approximation.h> 49 DEAL_II_NAMESPACE_OPEN
105 template <
class InputVector,
int spacedim>
108 const InputVector & solution,
109 const unsigned int component);
135 template <
class InputVector,
int spacedim>
139 const InputVector & solution,
140 const unsigned int component)
142 if (fe_values.
get_fe().n_components() == 1)
144 std::vector<typename InputVector::value_type> values(1);
150 std::vector<Vector<typename InputVector::value_type>> values(
152 Vector<typename InputVector::value_type>(
153 fe_values.
get_fe().n_components()));
155 return values[0](component);
166 for (
unsigned int i = 0; i < dim; ++i)
218 template <
class InputVector,
int spacedim>
221 const InputVector & solution,
222 const unsigned int component);
252 template <
class InputVector,
int spacedim>
256 const InputVector & solution,
257 const unsigned int component)
259 if (fe_values.
get_fe().n_components() == 1)
261 std::vector<Tensor<1, dim, typename InputVector::value_type>> values(
269 std::vector<Tensor<1, dim, typename InputVector::value_type>>>
273 fe_values.
get_fe().n_components()));
285 return std::fabs(d[0][0]);
302 const double radicand =
303 ::sqr(d[0][0] - d[1][1]) + 4 * ::sqr(d[0][1]);
305 0.5 * (d[0][0] + d[1][1] + std::sqrt(radicand)),
306 0.5 * (d[0][0] + d[1][1] - std::sqrt(radicand))};
308 return std::max(std::fabs(eigenvalues[0]), std::fabs(eigenvalues[1]));
426 const double am =
trace(d) / 3.;
430 for (
unsigned int i = 0; i < 3; ++i)
433 const double ss01 = s[0][1] * s[0][1], ss12 = s[1][2] * s[1][2],
434 ss02 = s[0][2] * s[0][2];
436 const double J2 = (s[0][0] * s[0][0] + s[1][1] * s[1][1] +
437 s[2][2] * s[2][2] + 2 * (ss01 + ss02 + ss12)) /
440 (std::pow(s[0][0], 3) + std::pow(s[1][1], 3) + std::pow(s[2][2], 3) +
441 3. * s[0][0] * (ss01 + ss02) + 3. * s[1][1] * (ss01 + ss12) +
442 3. * s[2][2] * (ss02 + ss12) + 6. * s[0][1] * s[0][2] * s[1][2]) /
445 const double R = std::sqrt(4. * J2 / 3.);
447 double EE[3] = {0, 0, 0};
454 if (R <= 1e-14 * std::fabs(am))
455 EE[0] = EE[1] = EE[2] = am;
460 const double R3 = R * R * R;
461 const double XX = 4. * J3 / R3;
462 const double YY = 1. - std::fabs(XX);
469 const double a = (XX > 0 ? -1. : 1.) * R / 2;
470 EE[0] = EE[1] = am + a;
475 const double theta = std::acos(XX) / 3.;
476 EE[0] = am + R * std::cos(theta);
477 EE[1] = am + R * std::cos(theta + 2. / 3. *
numbers::PI);
478 EE[2] = am + R * std::cos(theta + 4. / 3. *
numbers::PI);
482 return std::max(std::fabs(EE[0]),
483 std::max(std::fabs(EE[1]), std::fabs(EE[2])));
516 for (
unsigned int i = 0; i < dim; ++i)
517 for (
unsigned int j = i + 1; j < dim; ++j)
519 const double s = (d[i][j] + d[j][i]) / 2;
520 d[i][j] = d[j][i] = s;
527 class ThirdDerivative
555 template <
class InputVector,
int spacedim>
558 const InputVector & solution,
559 const unsigned int component);
589 template <
class InputVector,
int spacedim>
591 ThirdDerivative<dim>::get_projected_derivative(
593 const InputVector & solution,
594 const unsigned int component)
596 if (fe_values.
get_fe().n_components() == 1)
598 std::vector<Tensor<2, dim, typename InputVector::value_type>> values(
606 std::vector<Tensor<2, dim, typename InputVector::value_type>>>
610 fe_values.
get_fe().n_components()));
620 ThirdDerivative<1>::derivative_norm(
const Derivative &d)
622 return std::fabs(d[0][0][0]);
629 ThirdDerivative<dim>::derivative_norm(
const Derivative &d)
646 for (
unsigned int i = 0; i < dim; ++i)
647 for (
unsigned int j = i + 1; j < dim; ++j)
648 for (
unsigned int k = j + 1; k < dim; ++k)
650 const double s = (d[i][j][k] + d[i][k][j] + d[j][i][k] +
651 d[j][k][i] + d[k][i][j] + d[k][j][i]) /
653 d[i][j][k] = d[i][k][j] = d[j][i][k] = d[j][k][i] = d[k][i][j] =
658 for (
unsigned int i = 0; i < dim; ++i)
659 for (
unsigned int j = i + 1; j < dim; ++j)
663 const double s = (d[i][i][j] + d[i][j][i] + d[j][i][i]) / 3;
664 d[i][i][j] = d[i][j][i] = d[j][i][i] = s;
668 const double t = (d[i][j][j] + d[j][i][j] + d[j][j][i]) / 3;
669 d[i][j][j] = d[j][i][j] = d[j][j][i] = t;
674 template <
int order,
int dim>
675 class DerivativeSelector
683 using DerivDescr = void;
687 class DerivativeSelector<1, dim>
694 class DerivativeSelector<2, dim>
701 class DerivativeSelector<3, dim>
704 using DerivDescr = ThirdDerivative<dim>;
723 CopyData() =
default;
739 template <
class DerivativeDescription,
741 template <
int,
int>
class DoFHandlerType,
747 const DoFHandlerType<dim, spacedim> &dof_handler,
748 const InputVector & solution,
749 const unsigned int component,
752 typename DerivativeDescription::Derivative &derivative)
765 dof_handler.get_fe_collection();
793 typename DerivativeDescription::Derivative projected_derivative;
796 x_fe_midpoint_value.
reinit(cell);
802 const typename DerivativeDescription::ProjectedDerivative
803 this_midpoint_value =
804 DerivativeDescription::get_projected_derivative(fe_midpoint_value,
822 GridTools::get_active_neighbors<DoFHandlerType<dim, spacedim>>(
823 cell, active_neighbors);
829 ::DoFCellAccessor<DoFHandlerType<dim, spacedim>,
false>>>::
830 const_iterator neighbor_ptr = active_neighbors.begin();
831 for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr)
834 ::DoFCellAccessor<DoFHandlerType<dim, spacedim>,
false>>
835 neighbor = *neighbor_ptr;
838 x_fe_midpoint_value.
reinit(neighbor);
844 const typename DerivativeDescription::ProjectedDerivative
845 neighbor_midpoint_value =
846 DerivativeDescription::get_projected_derivative(
847 neighbor_fe_midpoint_value, solution, component);
860 const double distance = y.
norm();
874 for (
unsigned int i = 0; i < dim; ++i)
875 for (
unsigned int j = 0; j < dim; ++j)
876 Y[i][j] += y[i] * y[j];
881 typename DerivativeDescription::ProjectedDerivative
882 projected_finite_difference =
883 (neighbor_midpoint_value - this_midpoint_value);
884 projected_finite_difference /= distance;
886 projected_derivative +=
outer_product(y, projected_finite_difference);
903 derivative = Y_inverse * projected_derivative;
916 template <
class DerivativeDescription,
918 template <
int,
int>
class DoFHandlerType,
926 Vector<float>::iterator>>
const & cell,
928 const DoFHandlerType<dim, spacedim> &dof_handler,
929 const InputVector & solution,
930 const unsigned int component)
933 if (std::get<0>(*cell)->is_locally_owned() ==
false)
934 *std::get<1>(*cell) = 0;
937 typename DerivativeDescription::Derivative derivative;
940 approximate_cell<DerivativeDescription,
953 *std::get<1>(*cell) =
954 DerivativeDescription::derivative_norm(derivative);
969 template <
class DerivativeDescription,
971 template <
int,
int>
class DoFHandlerType,
976 const DoFHandlerType<dim, spacedim> &dof_handler,
977 const InputVector & solution,
978 const unsigned int component,
981 Assert(derivative_norm.size() ==
982 dof_handler.get_triangulation().n_active_cells(),
984 derivative_norm.size(),
985 dof_handler.get_triangulation().n_active_cells()));
986 Assert(component < dof_handler.get_fe(0).n_components(),
987 ExcIndexRange(component, 0, dof_handler.get_fe(0).n_components()));
989 using Iterators = std::tuple<
992 Vector<float>::iterator>;
994 Iterators(dof_handler.begin_active(), derivative_norm.begin())),
995 end(Iterators(dof_handler.end(), derivative_norm.end()));
1003 static_cast<std::function<void(SynchronousIterators<Iterators>
const &,
1004 Assembler::Scratch
const &,
1005 Assembler::CopyData &)
>>(
1006 std::bind(&approximate<DerivativeDescription,
1011 std::placeholders::_1,
1013 std::cref(dof_handler),
1014 std::cref(solution),
1016 std::function<
void(internal::Assembler::CopyData
const &)>(),
1017 internal::Assembler::Scratch(),
1018 internal::Assembler::CopyData());
1031 template <
int,
int>
class DoFHandlerType,
1036 const DoFHandlerType<dim, spacedim> &dof_handler,
1037 const InputVector & solution,
1039 const unsigned int component)
1041 internal::approximate_derivative<internal::Gradient<dim>, dim>(
1047 template <
int,
int>
class DoFHandlerType,
1052 const InputVector & solution,
1054 const unsigned int component)
1056 internal::approximate_derivative<internal::Gradient<dim>, dim>(
1066 template <
int,
int>
class DoFHandlerType,
1072 const DoFHandlerType<dim, spacedim> &dof_handler,
1073 const InputVector & solution,
1075 const unsigned int component)
1077 internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1083 template <
int,
int>
class DoFHandlerType,
1088 const DoFHandlerType<dim, spacedim> &dof_handler,
1089 const InputVector & solution,
1091 const unsigned int component)
1093 internal::approximate_derivative<internal::SecondDerivative<dim>, dim>(
1102 template <
typename DoFHandlerType,
class InputVector,
int order>
1107 const DoFHandlerType &dof,
1108 const InputVector & solution,
1110 const typename DoFHandlerType::active_cell_iterator &cell,
1116 const unsigned int component)
1118 internal::approximate_cell<
1119 typename internal::DerivativeSelector<order, DoFHandlerType::dimension>::
1120 DerivDescr>(mapping, dof, solution, component, cell, derivative);
1125 template <
typename DoFHandlerType,
class InputVector,
int order>
1128 const DoFHandlerType &dof,
1129 const InputVector & solution,
1131 const typename DoFHandlerType::active_cell_iterator &cell,
1137 const unsigned int component)
1142 DoFHandlerType::space_dimension>::mapping,
1152 template <
int dim,
int order>
1156 return internal::DerivativeSelector<order, dim>::DerivDescr::
1157 derivative_norm(derivative);
1164 #include "derivative_approximation.inst" 1167 DEAL_II_NAMESPACE_CLOSE
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
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)
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type >> &gradients) const
static const UpdateFlags update_flags
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Transformed quadrature points.
numbers::NumberTraits< Number >::real_type norm() const
#define AssertThrow(cond, exc)
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 > &)
static const UpdateFlags update_flags
#define Assert(cond, exc)
static double derivative_norm(const Derivative &d)
Abstract base class for mapping classes.
const FiniteElement< dim, spacedim > & get_fe() const
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)
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)
Second derivatives of shape functions.
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)
Shape function gradients.
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)
const ::FEValues< dim, spacedim > & get_present_fe_values() const
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)
const Point< spacedim > & quadrature_point(const unsigned int q) const
Tensor< 0, dim > ProjectedDerivative
static::ExceptionBase & ExcInternalError()
static::ExceptionBase & ExcVectorLengthVsNActiveCells(int arg1, int arg2)
void get_function_hessians(const InputVector &fe_function, std::vector< Tensor< 2, spacedim, typename InputVector::value_type >> &hessians) const