17 #ifndef dealii_matrix_free_evaluation_kernels_h
18 #define dealii_matrix_free_evaluation_kernels_h
37 template <MatrixFreeFunctions::ElementType element,
bool is_
long>
41 template <
bool is_
long>
59 template <
bool is_
long>
78 template <
bool is_
long>
128 evaluate(
const unsigned int n_components,
130 const Number *values_dofs_actual,
134 integrate(
const unsigned int n_components,
136 Number *values_dofs_actual,
138 const bool add_into_values_array);
143 *univariate_shape_data)
166 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
174 evaluate(
const unsigned int n_components,
176 const Number *values_dofs_actual,
180 integrate(
const unsigned int n_components,
182 Number *values_dofs_actual,
184 const bool add_into_values_array);
196 const unsigned int n_components,
198 const Number *values_dofs_actual,
204 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
205 univariate_shape_data;
209 univariate_shape_data.fill(&shape_data.front());
211 if (shape_data.size() == dim)
212 for (
int i = 1; i < dim; ++i)
213 univariate_shape_data[i] = &shape_data[i];
215 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
216 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
217 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
219 const unsigned int temp_size =
222 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
223 Eval::n_rows_of_product :
224 Eval::n_columns_of_product);
229 temp2 = temp1 +
std::max(Utilities::fixed_power<dim>(
230 shape_data.front().fe_degree + 1),
231 Utilities::fixed_power<dim>(
232 shape_data.front().n_q_points_1d));
236 temp2 = temp1 + temp_size;
239 const std::size_t n_q_points = temp_size == 0 ?
241 Eval::n_columns_of_product;
242 const std::size_t dofs_per_comp =
246 const Number *values_dofs = values_dofs_actual;
249 const std::size_t n_dofs_per_comp =
251 Number *values_dofs_tmp =
252 temp1 + 2 * (
std::max(n_dofs_per_comp, n_q_points));
254 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
255 for (
unsigned int c = 0; c < n_components; ++c)
256 for (
int i = 0, count_p = 0, count_q = 0;
257 i < (dim > 2 ? degree + 1 : 1);
260 for (
int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
262 for (
int k = 0; k < degree + 1 - j - i;
263 ++k, ++count_p, ++count_q)
264 values_dofs_tmp[c * dofs_per_comp + count_q] =
265 values_dofs_actual[c * n_dofs_per_comp + count_p];
266 for (
int k = degree + 1 - j - i; k < degree + 1;
268 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
270 for (
int j = degree + 1 - i; j < degree + 1; ++j)
271 for (
int k = 0; k < degree + 1; ++k, ++count_q)
272 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
274 values_dofs = values_dofs_tmp;
283 for (
unsigned int c = 0; c < n_components; ++c)
286 eval0.template values<0, true, false>(values_dofs, values_quad);
288 eval0.template gradients<0, true, false>(values_dofs,
292 values_dofs += dofs_per_comp;
293 values_quad += n_q_points;
294 gradients_quad += n_q_points;
299 for (
unsigned int c = 0; c < n_components; ++c)
304 eval0.template gradients<0, true, false>(values_dofs, temp1);
305 eval1.template values<1, true, false, 2>(temp1,
310 eval0.template values<0, true, false>(values_dofs, temp1);
312 eval1.template gradients<1, true, false, 2>(temp1,
317 eval1.template values<1, true, false>(temp1, values_quad);
320 values_dofs += dofs_per_comp;
321 values_quad += n_q_points;
322 gradients_quad += 2 * n_q_points;
327 for (
unsigned int c = 0; c < n_components; ++c)
332 eval0.template gradients<0, true, false>(values_dofs, temp1);
333 eval1.template values<1, true, false>(temp1, temp2);
334 eval2.template values<2, true, false, 3>(temp2,
339 eval0.template values<0, true, false>(values_dofs, temp1);
342 eval1.template gradients<1, true, false>(temp1, temp2);
343 eval2.template values<2, true, false, 3>(temp2,
349 eval1.template values<1, true, false>(temp1, temp2);
351 eval2.template gradients<2, true, false, 3>(temp2,
357 eval2.template values<2, true, false>(temp2, values_quad);
360 values_dofs += dofs_per_comp;
361 values_quad += n_q_points;
362 gradients_quad += 3 * n_q_points;
375 values_quad -= n_components * n_q_points;
376 values_dofs -= n_components * dofs_per_comp;
377 for (std::size_t c = 0; c < n_components; ++c)
378 for (std::size_t q = 0; q < n_q_points; ++q)
379 values_quad[c * n_q_points + q] +=
380 values_dofs[(c + 1) * dofs_per_comp - 1];
393 const unsigned int n_components,
395 Number *values_dofs_actual,
397 const bool add_into_values_array)
399 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
400 univariate_shape_data;
403 univariate_shape_data.fill(&shape_data.front());
405 if (shape_data.size() == dim)
406 for (
int i = 1; i < dim; ++i)
407 univariate_shape_data[i] = &shape_data[i];
409 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
410 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
411 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
413 const unsigned int temp_size =
416 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
417 Eval::n_rows_of_product :
418 Eval::n_columns_of_product);
423 temp2 = temp1 +
std::max(Utilities::fixed_power<dim>(
424 shape_data.front().fe_degree + 1),
425 Utilities::fixed_power<dim>(
426 shape_data.front().n_q_points_1d));
430 temp2 = temp1 + temp_size;
433 const std::size_t n_q_points = temp_size == 0 ?
435 Eval::n_columns_of_product;
436 const unsigned int dofs_per_comp =
438 Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
441 Number *values_dofs =
443 temp1 + 2 * (std::max<std::size_t>(
454 for (
unsigned int c = 0; c < n_components; ++c)
458 if (add_into_values_array ==
false)
459 eval0.template values<0, false, false>(values_quad,
462 eval0.template values<0, false, true>(values_quad,
468 add_into_values_array ==
true)
469 eval0.template gradients<0, false, true>(gradients_quad,
472 eval0.template gradients<0, false, false>(gradients_quad,
477 values_dofs += dofs_per_comp;
478 values_quad += n_q_points;
479 gradients_quad += n_q_points;
484 for (
unsigned int c = 0; c < n_components; ++c)
489 eval1.template values<1, false, false>(values_quad, temp1);
490 if (add_into_values_array ==
false)
491 eval0.template values<0, false, false>(temp1, values_dofs);
493 eval0.template values<0, false, true>(temp1, values_dofs);
497 eval1.template gradients<1, false, false, 2>(gradients_quad +
501 eval1.template values<1, false, true>(values_quad, temp1);
502 if (add_into_values_array ==
false)
503 eval0.template values<0, false, false>(temp1, values_dofs);
505 eval0.template values<0, false, true>(temp1, values_dofs);
506 eval1.template values<1, false, false, 2>(gradients_quad,
508 eval0.template gradients<0, false, true>(temp1, values_dofs);
512 values_dofs += dofs_per_comp;
513 values_quad += n_q_points;
514 gradients_quad += 2 * n_q_points;
519 for (
unsigned int c = 0; c < n_components; ++c)
524 eval2.template values<2, false, false>(values_quad, temp1);
525 eval1.template values<1, false, false>(temp1, temp2);
526 if (add_into_values_array ==
false)
527 eval0.template values<0, false, false>(temp2, values_dofs);
529 eval0.template values<0, false, true>(temp2, values_dofs);
533 eval2.template gradients<2, false, false, 3>(gradients_quad +
537 eval2.template values<2, false, true>(values_quad, temp1);
538 eval1.template values<1, false, false>(temp1, temp2);
539 eval2.template values<2, false, false, 3>(gradients_quad + 1,
541 eval1.template gradients<1, false, true>(temp1, temp2);
542 if (add_into_values_array ==
false)
543 eval0.template values<0, false, false>(temp2, values_dofs);
545 eval0.template values<0, false, true>(temp2, values_dofs);
546 eval2.template values<2, false, false, 3>(gradients_quad,
548 eval1.template values<1, false, false>(temp1, temp2);
549 eval0.template gradients<0, false, true>(temp2, values_dofs);
553 values_dofs += dofs_per_comp;
554 values_quad += n_q_points;
555 gradients_quad += 3 * n_q_points;
566 values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
567 values_quad -= n_components * n_q_points;
569 for (
unsigned int c = 0; c < n_components; ++c)
571 values_dofs[0] = values_quad[0];
572 for (
unsigned int q = 1; q < n_q_points; ++q)
573 values_dofs[0] += values_quad[q];
574 values_dofs += dofs_per_comp;
575 values_quad += n_q_points;
579 for (
unsigned int c = 0; c < n_components; ++c)
580 values_dofs[c * dofs_per_comp] = Number();
581 values_dofs += n_components * dofs_per_comp;
587 const std::size_t n_dofs_per_comp =
589 values_dofs -= dofs_per_comp * n_components;
591 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
592 for (
unsigned int c = 0; c < n_components; ++c)
593 for (
int i = 0, count_p = 0, count_q = 0;
594 i < (dim > 2 ? degree + 1 : 1);
597 for (
int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
599 for (
int k = 0; k < degree + 1 - j - i;
600 ++k, ++count_p, ++count_q)
601 values_dofs_actual[c * n_dofs_per_comp + count_p] =
602 values_dofs[c * dofs_per_comp + count_q];
605 count_q += i * (degree + 1);
612 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
619 Number>::evaluate(
const unsigned int n_components,
621 const Number *values_dofs_actual,
626 const std::size_t n_dofs =
639 const auto *
const shape_values = shape_data.front().
shape_values.data();
641 const auto *values_dofs_actual_ptr = values_dofs_actual;
643 Eval eval(shape_values,
nullptr,
nullptr, n_dofs, n_q_points);
644 for (
unsigned int c = 0; c < n_components; ++c)
646 eval.template values<0, true, false>(values_dofs_actual_ptr,
649 values_quad_ptr += n_q_points;
650 values_dofs_actual_ptr += n_dofs;
656 const auto *
const shape_gradients =
657 shape_data.front().shape_gradients.data();
659 const auto *values_dofs_actual_ptr = values_dofs_actual;
661 for (
unsigned int c = 0; c < n_components; ++c)
663 for (
unsigned int d = 0;
d < dim; ++
d)
666 shape_gradients + n_q_points * n_dofs *
d,
671 eval.template gradients<0, true, false, dim>(
672 values_dofs_actual_ptr, gradients_quad_ptr +
d);
674 gradients_quad_ptr += n_q_points * dim;
675 values_dofs_actual_ptr += n_dofs;
682 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
689 Number>::integrate(
const unsigned int n_components,
691 Number *values_dofs_actual,
693 const bool add_into_values_array)
698 const std::size_t n_dofs =
711 const auto *
const shape_values = shape_data.front().
shape_values.data();
713 auto *values_dofs_actual_ptr = values_dofs_actual;
715 Eval eval(shape_values,
nullptr,
nullptr, n_dofs, n_q_points);
716 for (
unsigned int c = 0; c < n_components; ++c)
718 if (add_into_values_array ==
false)
719 eval.template values<0, false, false>(values_quad_ptr,
720 values_dofs_actual_ptr);
722 eval.template values<0, false, true>(values_quad_ptr,
723 values_dofs_actual_ptr);
725 values_quad_ptr += n_q_points;
726 values_dofs_actual_ptr += n_dofs;
732 const auto *
const shape_gradients =
733 shape_data.front().shape_gradients.data();
735 auto *values_dofs_actual_ptr = values_dofs_actual;
737 for (
unsigned int c = 0; c < n_components; ++c)
739 for (
unsigned int d = 0;
d < dim; ++
d)
742 shape_gradients + n_q_points * n_dofs *
d,
747 if ((add_into_values_array ==
false &&
750 eval.template gradients<0, false, false, dim>(
751 gradients_quad_ptr +
d, values_dofs_actual_ptr);
753 eval.template gradients<0, false, true, dim>(
754 gradients_quad_ptr +
d, values_dofs_actual_ptr);
756 gradients_quad_ptr += n_q_points * dim;
757 values_dofs_actual_ptr += n_dofs;
780 static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
781 "The second dimension must not be smaller than the first");
805 template <
typename Number,
typename Number2>
812 const Number *values_in,
814 const unsigned int basis_size_1_variable =
816 const unsigned int basis_size_2_variable =
820 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
821 ExcMessage(
"The second dimension must not be smaller than the first"));
829 constexpr
int next_dim = (dim == 1 || (dim == 2 && basis_size_1 > 0 &&
830 basis_size_1 == basis_size_2)) ?
837 (basis_size_1 == 0 ? 0 : basis_size_2),
840 eval_val(transformation_matrix,
843 basis_size_1_variable,
844 basis_size_2_variable);
845 const unsigned int np_1 =
846 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
847 const unsigned int np_2 =
848 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
850 ExcMessage(
"Cannot transform with 0-point basis"));
852 ExcMessage(
"Cannot transform with 0-point basis"));
856 values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
858 values_out + n_components * Utilities::fixed_power<dim>(np_2);
859 for (
unsigned int c = n_components; c != 0; --c)
861 values_in -= Utilities::fixed_power<dim>(np_1);
862 values_out -= Utilities::fixed_power<dim>(np_2);
864 for (
unsigned int q = np_1; q != 0; --q)
871 transformation_matrix,
873 (q - 1) * Utilities::fixed_power<next_dim>(np_1),
875 (q - 1) * Utilities::fixed_power<next_dim>(np_2),
876 basis_size_1_variable,
877 basis_size_2_variable);
882 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
884 eval_val.template values<0, true, false>(values_in, values_out);
885 eval_val.template values<1, true, false>(values_out, values_out);
888 eval_val.template
values<dim - 1,
true,
false>(values_in,
891 eval_val.template
values<dim - 1,
true,
false>(values_out,
926 template <
typename Number,
typename Number2>
933 const bool add_into_result,
936 const unsigned int basis_size_1_variable =
938 const unsigned int basis_size_2_variable =
942 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
943 ExcMessage(
"The second dimension must not be smaller than the first"));
944 Assert(add_into_result ==
false || values_in != values_out,
946 "Input and output cannot alias with each other when "
947 "adding the result of the basis change to existing data"));
953 constexpr
int next_dim =
955 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
961 (basis_size_1 == 0 ? 0 : basis_size_2),
964 eval_val(transformation_matrix,
965 transformation_matrix,
966 transformation_matrix,
967 basis_size_1_variable,
968 basis_size_2_variable);
969 const unsigned int np_1 =
970 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
971 const unsigned int np_2 =
972 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
974 ExcMessage(
"Cannot transform with 0-point basis"));
976 ExcMessage(
"Cannot transform with 0-point basis"));
978 for (
unsigned int c = 0; c < n_components; ++c)
980 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
983 eval_val.template values<1, false, false>(values_in, values_in);
985 eval_val.template hessians<1, false, false>(values_in,
991 eval_val.template values<0, false, true>(values_in,
994 eval_val.template hessians<0, false, true>(values_in,
1000 eval_val.template values<0, false, false>(values_in,
1003 eval_val.template hessians<0, false, false>(values_in,
1009 if (dim == 1 && add_into_result)
1012 eval_val.template values<0, false, true>(values_in,
1015 eval_val.template hessians<0, false, true>(values_in,
1021 eval_val.template values<0, false, false>(values_in,
1024 eval_val.template hessians<0, false, false>(values_in,
1030 eval_val.template
values<dim - 1,
false,
false>(values_in,
1033 eval_val.template
hessians<dim - 1,
false,
false>(
1034 values_in, values_in);
1038 for (
unsigned int q = 0; q < np_1; ++q)
1045 transformation_matrix,
1048 q * Utilities::fixed_power<next_dim>(np_2),
1050 q * Utilities::fixed_power<next_dim>(np_1),
1051 basis_size_1_variable,
1052 basis_size_2_variable);
1054 values_in += Utilities::fixed_power<dim>(np_2);
1055 values_out += Utilities::fixed_power<dim>(np_1);
1079 template <
typename Number,
typename Number2>
1084 const Number *values_in,
1085 Number *scratch_data,
1088 constexpr
int next_dim = dim > 1 ? dim - 1 : dim;
1089 Number *my_scratch =
1090 basis_size_1 != basis_size_2 ? scratch_data : values_out;
1092 const unsigned int size_per_component =
Utilities::pow(basis_size_2, dim);
1093 Assert(coefficients.
size() == size_per_component ||
1094 coefficients.
size() == n_components * size_per_component,
1096 const unsigned int stride =
1097 coefficients.
size() == size_per_component ? 0 : 1;
1099 for (
unsigned int q = basis_size_1; q != 0; --q)
1106 transformation_matrix,
1119 eval_val(transformation_matrix);
1120 const unsigned int n_inner_blocks =
1121 (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1123 for (
unsigned int ii = 0; ii <
n_blocks; ii += n_inner_blocks)
1124 for (
unsigned int c = 0; c < n_components; ++c)
1126 for (
unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1127 eval_val.template values_one_line<dim - 1, true, false>(
1128 my_scratch + i, my_scratch + i);
1129 for (
unsigned int q = 0; q < basis_size_2; ++q)
1130 for (
unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1131 my_scratch[i + q *
n_blocks + c * size_per_component] *=
1133 c * stride * size_per_component];
1134 for (
unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1135 eval_val.template values_one_line<dim - 1, false, false>(
1136 my_scratch + i, my_scratch + i);
1138 for (
unsigned int q = 0; q < basis_size_1; ++q)
1145 transformation_matrix,
1161 template <
int n_po
ints_1d,
int dim,
typename Number,
typename Number2>
1169 (n_points_1d + 1) / 2 * n_points_1d);
1192 constexpr
unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
1193 constexpr
unsigned int n_points_2d = n_points_1d * n_points_1d;
1194 const Number *in =
values + (loop_bound - 1) * n_points_2d;
1195 Number *out =
gradients + (loop_bound - 1) * dim * n_points_2d;
1196 for (
unsigned int l = 0;
l < loop_bound; ++
l)
1198 eval_2d.template gradients<0, true, false, dim>(in, out);
1199 eval_2d.template gradients<1, true, false, dim>(in, out + 1);
1201 out -= dim * n_points_2d;
1215 template <
int n_po
ints_1d,
int dim,
typename Number,
typename Number2>
1221 const bool add_into_values_array)
1224 (n_points_1d + 1) / 2 * n_points_1d);
1243 if (add_into_values_array)
1250 constexpr
unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
1251 constexpr
unsigned int n_points_2d = n_points_1d * n_points_1d;
1253 const Number *in =
gradients + (loop_bound - 1) * dim * n_points_2d;
1254 Number *out =
values + (loop_bound - 1) * n_points_2d;
1255 for (
unsigned int l = 0;
l < loop_bound; ++
l)
1257 if (add_into_values_array)
1258 eval_2d.template gradients<0, false, true, dim>(in, out);
1260 eval_2d.template gradients<0, false, false, dim>(in, out);
1261 eval_2d.template gradients<1, false, true, dim>(in + 1, out);
1262 in -= dim * n_points_2d;
1278 template <
int n_po
ints_1d,
int dim,
typename Number>
1309 for (
unsigned int comp = 0; comp < n_components; ++comp)
1318 eval.template gradients<0, true, false>(
values, scratch);
1319 eval.template gradients<1, true, false>(scratch,
1322 eval.template hessians<1, true, false>(
values,
hessians + n_points);
1327 eval.template gradients<2, true, false>(scratch,
1330 eval.template gradients<1, true, false>(
values, scratch);
1331 eval.template gradients<2, true, false>(scratch,
1334 eval.template hessians<2, true, false>(
values,
1339 hessians += (dim * (dim + 1)) / 2 * n_points;
1351 template <
int n_q_po
ints_1d,
int dim,
typename Number>
1355 const bool add_into_values_array)
1380 for (
unsigned int comp = 0; comp < n_components; ++comp)
1383 if (add_into_values_array ==
true)
1390 eval.template hessians<1, false, true>(
hessians + n_points,
values);
1394 eval.template hessians<2, false, true>(
hessians + 2 * n_points,
1397 eval.template gradients<2, false, false>(
hessians + 5 * n_points,
1399 eval.template gradients<1, false, true>(scratch,
values);
1402 eval.template gradients<2, false, false>(
hessians + 4 * n_points,
1412 eval.template gradients<0, false, true>(scratch,
values);
1416 hessians += (dim * (dim + 1)) / 2 * n_points;
1428 template <
int dim,
typename Number>
1431 const Number *values_dofs,
1437 using Eval =
typename Impl::Eval;
1439 Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
1440 Eval eval1 = Impl::create_evaluator_tensor_product(
1441 &univariate_shape_data[std::min<int>(1,
1442 univariate_shape_data.size() - 1)]);
1443 Eval eval2 = Impl::create_evaluator_tensor_product(
1444 &univariate_shape_data[std::min<int>(2,
1445 univariate_shape_data.size() - 1)]);
1450 tmp1 +
std::max(Utilities::fixed_power<dim>(
1451 univariate_shape_data.front().fe_degree + 1),
1452 Utilities::fixed_power<dim>(
1453 univariate_shape_data.front().n_q_points_1d));
1456 for (
unsigned int comp = 0; comp < n_components;
1458 hessians += n_points * dim * (dim + 1) / 2,
1464 eval0.template hessians<0, true, false>(values_dofs,
hessians);
1468 eval0.template hessians<0, true, false>(values_dofs, tmp1);
1469 eval1.template values<1, true, false>(tmp1,
hessians);
1471 eval0.template gradients<0, true, false>(values_dofs, tmp1);
1472 eval1.template gradients<1, true, false>(tmp1,
1475 eval0.template values<0, true, false>(values_dofs, tmp1);
1476 eval1.template hessians<1, true, false>(tmp1,
hessians + n_points);
1480 eval0.template hessians<0, true, false>(values_dofs, tmp1);
1481 eval1.template values<1, true, false>(tmp1, tmp2);
1482 eval2.template values<2, true, false>(tmp2,
hessians);
1484 eval0.template gradients<0, true, false>(values_dofs, tmp1);
1485 eval1.template gradients<1, true, false>(tmp1, tmp2);
1486 eval2.template values<2, true, false>(tmp2,
1489 eval1.template values<1, true, false>(tmp1, tmp2);
1490 eval2.template gradients<2, true, false>(tmp2,
1493 eval0.template values<0, true, false>(values_dofs, tmp1);
1494 eval1.template hessians<1, true, false>(tmp1, tmp2);
1495 eval2.template values<2, true, false>(tmp2,
hessians + n_points);
1497 eval1.template gradients<1, true, false>(tmp1, tmp2);
1498 eval2.template gradients<2, true, false>(tmp2,
1501 eval1.template values<1, true, false>(tmp1, tmp2);
1502 eval2.template hessians<2, true, false>(tmp2,
1509 "Only 1d, 2d and 3d implemented for Hessian"));
1522 template <
int dim,
typename Number>
1526 Number *values_dofs,
1527 const bool add_into_values_array)
1532 using Eval =
typename Impl::Eval;
1534 Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
1535 Eval eval1 = Impl::create_evaluator_tensor_product(
1536 &univariate_shape_data[std::min<int>(1,
1537 univariate_shape_data.size() - 1)]);
1538 Eval eval2 = Impl::create_evaluator_tensor_product(
1539 &univariate_shape_data[std::min<int>(2,
1540 univariate_shape_data.size() - 1)]);
1545 tmp1 +
std::max(Utilities::fixed_power<dim>(
1546 univariate_shape_data.front().fe_degree + 1),
1547 Utilities::fixed_power<dim>(
1548 univariate_shape_data.front().n_q_points_1d));
1551 for (
unsigned int comp = 0; comp < n_components;
1553 hessians += n_points * dim * (dim + 1) / 2,
1559 if (add_into_values_array)
1560 eval0.template hessians<0, false, true>(
hessians, values_dofs);
1562 eval0.template hessians<0, false, false>(
hessians, values_dofs);
1566 eval1.template values<1, false, false>(
hessians, tmp1);
1567 if (add_into_values_array)
1568 eval0.template hessians<0, false, true>(tmp1, values_dofs);
1570 eval0.template hessians<0, false, false>(tmp1, values_dofs);
1573 eval1.template gradients<1, false, false>(
hessians + 2 * n_points,
1575 eval0.template gradients<0, false, true>(tmp1, values_dofs);
1577 eval1.template hessians<1, false, false>(
hessians + n_points, tmp1);
1578 eval0.template values<0, false, true>(tmp1, values_dofs);
1582 eval2.template values<2, false, false>(
hessians, tmp1);
1583 eval1.template values<1, false, false>(tmp1, tmp2);
1585 if (add_into_values_array)
1586 eval0.template hessians<0, false, true>(tmp2, values_dofs);
1588 eval0.template hessians<0, false, false>(tmp2, values_dofs);
1591 eval2.template values<2, false, false>(
hessians + 3 * n_points,
1593 eval1.template gradients<1, false, false>(tmp1, tmp2);
1595 eval2.template gradients<2, false, false>(
hessians + 4 * n_points,
1597 eval1.template values<1, false, true>(tmp1, tmp2);
1598 eval1.template values<0, false, true>(tmp2, values_dofs);
1601 eval2.template values<2, false, false>(
hessians + n_points, tmp1);
1602 eval1.template hessians<1, false, false>(tmp1, tmp2);
1605 eval2.template gradients<2, false, false>(
hessians + 5 * n_points,
1607 eval1.template gradients<1, false, true>(tmp1, tmp2);
1610 eval2.template hessians<2, false, false>(
hessians + 2 * n_points,
1612 eval1.template values<1, false, true>(tmp1, tmp2);
1613 eval0.template values<0, false, true>(tmp2, values_dofs);
1619 "Only 1d, 2d and 3d implemented for Hessian"));
1637 template <
int dim,
int fe_degree,
typename Number>
1652 const Number *values_dofs,
1655 constexpr std::size_t n_points =
Utilities::pow(fe_degree + 1, dim);
1657 for (
unsigned int c = 0; c < n_components; ++c)
1660 for (
unsigned int i = 0; i < n_points; ++i)
1662 values_dofs[n_points * c + i];
1665 evaluate_gradients_collocation<fe_degree + 1, dim>(
1667 values_dofs + c * n_points,
1675 Number *values_dofs,
1677 const bool add_into_values_array)
1679 constexpr std::size_t n_points =
Utilities::pow(fe_degree + 1, dim);
1681 for (
unsigned int c = 0; c < n_components; ++c)
1685 if (add_into_values_array)
1686 for (
unsigned int i = 0; i < n_points; ++i)
1687 values_dofs[n_points * c + i] +=
1690 for (
unsigned int i = 0; i < n_points; ++i)
1691 values_dofs[n_points * c + i] =
1696 integrate_gradients_collocation<fe_degree + 1, dim>(
1698 values_dofs + c * n_points,
1700 add_into_values_array ||
1718 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
1724 const Number *values_dofs,
1729 Assert(n_q_points_1d > fe_degree,
1730 ExcMessage(
"You lose information when going to a collocation "
1731 "space of lower degree, so the evaluation results "
1732 "would be wrong. Thus, this class does not permit "
1733 "the chosen operation."));
1734 constexpr std::size_t n_dofs =
Utilities::pow(fe_degree + 1, dim);
1735 constexpr std::size_t n_q_points =
Utilities::pow(n_q_points_1d, dim);
1737 for (
unsigned int c = 0; c < n_components; ++c)
1743 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1744 n_q_points_1d>::do_forward(1,
1745 shape_data.shape_values_eo,
1746 values_dofs + c * n_dofs,
1751 evaluate_gradients_collocation<n_q_points_1d, dim>(
1761 Number *values_dofs,
1763 const bool add_into_values_array)
1767 Assert(n_q_points_1d > fe_degree,
1768 ExcMessage(
"You lose information when going to a collocation "
1769 "space of lower degree, so the evaluation results "
1770 "would be wrong. Thus, this class does not permit "
1771 "the chosen operation."));
1772 constexpr std::size_t n_q_points =
Utilities::pow(n_q_points_1d, dim);
1774 for (
unsigned int c = 0; c < n_components; ++c)
1778 integrate_gradients_collocation<n_q_points_1d, dim>(
1790 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1791 n_q_points_1d>::do_backward(1,
1792 shape_data.shape_values_eo,
1793 add_into_values_array,
1808 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
1818 template <
bool integrate>
1820 evaluate_or_integrate(
1822 Number *values_dofs_actual,
1824 const bool add_into_values_array =
false);
1829 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
1830 template <
bool integrate>
1837 evaluate_or_integrate(
1839 Number *values_dofs,
1843 Assert(dim == 2 || dim == 3,
1844 ExcMessage(
"Only dim = 2,3 implemented for Raviart-Thomas "
1845 "evaluation/integration"));
1859 const unsigned int dofs_per_component =
1861 const unsigned int n_points =
Utilities::pow(n_q_points_1d, dim);
1872 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1876 if constexpr (dim > 2)
1877 eval.template tangential<2, 0>(shape_data[1],
values,
values);
1878 eval.template tangential<1, 0>(shape_data[1],
values,
values);
1879 eval.template normal<0>(shape_data[0],
values, values_dofs, add);
1883 values_dofs += dofs_per_component;
1886 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1890 if constexpr (dim > 2)
1891 eval.template tangential<2, 1>(shape_data[1],
values,
values);
1892 eval.template tangential<0, 1>(shape_data[1],
values,
values);
1893 eval.template normal<1>(shape_data[0],
values, values_dofs, add);
1895 if constexpr (dim > 2)
1899 values_dofs += dofs_per_component;
1902 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1906 eval.template tangential<1, 2>(shape_data[1],
values,
values);
1907 eval.template tangential<0, 2>(shape_data[1],
values,
values);
1908 eval.template normal<0>(shape_data[0],
values, values_dofs, add);
1915 eval.template normal<0>(shape_data[0], values_dofs,
values);
1916 eval.template tangential<1, 0>(shape_data[1],
values,
values);
1917 if constexpr (dim > 2)
1918 eval.template tangential<2, 0>(shape_data[1],
values,
values);
1920 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1926 values_dofs += dofs_per_component;
1928 eval.template normal<1>(shape_data[0], values_dofs,
values);
1929 eval.template tangential<0, 1>(shape_data[1],
values,
values);
1930 if constexpr (dim > 2)
1931 eval.template tangential<2, 1>(shape_data[1],
values,
values);
1933 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1937 if constexpr (dim > 2)
1941 values_dofs += dofs_per_component;
1943 eval.template normal<2>(shape_data[0], values_dofs,
values);
1944 eval.template tangential<0, 2>(shape_data[1],
values,
values);
1945 eval.template tangential<1, 2>(shape_data[1],
values,
values);
1947 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1971 template <
int dim,
typename Number,
bool do_
integrate>
1974 template <
int fe_degree,
int n_q_po
ints_1d,
typename OtherNumber>
1976 run(
const unsigned int n_components,
1978 OtherNumber *values_dofs,
1980 const bool sum_into_values_array_in =
false)
1984 static_assert(std::is_same_v<Number, std::remove_const_t<OtherNumber>>,
1985 "Type of Number and of OtherNumber do not match.");
1997 bool sum_into_values_array = sum_into_values_array_in;
2003 if constexpr (do_integrate)
2007 integrate_hessians_collocation<n_q_points_1d>(
2016 sum_into_values_array);
2017 sum_into_values_array =
true;
2022 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2031 sum_into_values_array);
2035 else if (fe_degree >= 0 &&
2048 sum_into_values_array);
2050 else if (fe_degree >= 0 &&
2062 sum_into_values_array);
2071 Number>>(n_components,
2075 sum_into_values_array);
2088 sum_into_values_array);
2101 sum_into_values_array);
2105 if constexpr (fe_degree > 0 && n_q_points_1d > 0 && dim > 1)
2112 template evaluate_or_integrate<do_integrate>(
2114 const_cast<Number *
>(values_dofs),
2116 sum_into_values_array);
2122 "in 2d/3d and with templated degree"));
2136 sum_into_values_array);
2145 evaluate_hessians_collocation<n_q_points_1d>(n_components, fe_eval);
2154 template <
typename T>
2157 const unsigned int n_components,
2159 const Number *values_dofs,
2161 const bool sum_into_values_array,
2162 std::bool_constant<false>)
2164 (void)sum_into_values_array;
2166 T::evaluate(n_components, evaluation_flag, values_dofs, fe_eval);
2169 template <
typename T>
2172 const unsigned int n_components,
2174 Number *values_dofs,
2176 const bool sum_into_values_array,
2177 std::bool_constant<true>)
2179 T::integrate(n_components,
2183 sum_into_values_array);
2186 template <
typename T,
typename OtherNumber>
2189 const unsigned int n_components,
2191 OtherNumber *values_dofs,
2193 const bool sum_into_values_array)
2195 evaluate_or_integrate<T>(n_components,
2199 sum_into_values_array,
2200 std::bool_constant<do_integrate>());
2210 template <
int dim,
typename Number>
2216 template <
int fe_degree,
int = 0>
2218 run(
const unsigned int n_components,
2220 const Number *in_array,
2223 const unsigned int given_degree =
2224 (fe_degree > -1) ? fe_degree :
2227 const unsigned int dofs_per_component =
2247 for (
unsigned int d = 0;
d < n_components; ++
d)
2249 const Number *in = in_array +
d * dofs_per_component;
2250 Number *out = out_array +
d * dofs_per_component;
2253 evaluator.template hessians<0, true, false>(in, out);
2255 evaluator.template hessians<1, true, false>(out, out);
2257 evaluator.template hessians<2, true, false>(out, out);
2259 for (
unsigned int q = 0; q < dofs_per_component; ++q)
2261 const Number inverse_JxW_q = Number(1.) / fe_eval.
JxW(q);
2262 for (
unsigned int d = 0;
d < n_components; ++
d)
2263 out_array[q +
d * dofs_per_component] *= inverse_JxW_q;
2265 for (
unsigned int d = 0;
d < n_components; ++
d)
2267 Number *out = out_array +
d * dofs_per_component;
2269 evaluator.template hessians<2, false, false>(out, out);
2271 evaluator.template hessians<1, false, false>(out, out);
2272 evaluator.template hessians<0, false, false>(out, out);
2286 template <
int dim,
typename Number>
2292 template <
int fe_degree,
int = 0>
2294 run(
const unsigned int n_desired_components,
2297 const bool dyadic_coefficients,
2298 const Number *in_array,
2301 const unsigned int given_degree =
2302 (fe_degree > -1) ? fe_degree :
2305 const unsigned int dofs_per_component =
2309 inverse_coefficients.
size() % dofs_per_component == 0,
2311 "Expected diagonal to be a multiple of scalar dof per cells"));
2313 if (!dyadic_coefficients)
2315 if (inverse_coefficients.
size() != dofs_per_component)
2317 inverse_coefficients.
size());
2323 inverse_coefficients.
size());
2343 const Number *in = in_array;
2344 Number *out = out_array;
2346 const Number *inv_coefficient = inverse_coefficients.
data();
2348 const unsigned int shift_coefficient =
2349 inverse_coefficients.
size() > dofs_per_component ? dofs_per_component :
2352 const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
2353 const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
2355 for (
unsigned int d = 0;
d < n_comp_outer; ++
d)
2357 for (
unsigned int di = 0; di < n_comp_inner; ++di)
2359 const Number *in_ = in + di * dofs_per_component;
2360 Number *out_ = out + di * dofs_per_component;
2361 evaluator.template hessians<0, true, false>(in_, out_);
2363 evaluator.template hessians<1, true, false>(out_, out_);
2365 evaluator.template hessians<2, true, false>(out_, out_);
2367 if (dyadic_coefficients)
2369 const auto n_coeff_components =
2370 n_desired_components * n_desired_components;
2371 if (n_desired_components == dim)
2373 for (
unsigned int q = 0; q < dofs_per_component; ++q)
2374 vmult<dim>(&inv_coefficient[q * n_coeff_components],
2377 dofs_per_component);
2381 for (
unsigned int q = 0; q < dofs_per_component; ++q)
2382 vmult<-1>(&inv_coefficient[q * n_coeff_components],
2386 n_desired_components);
2390 for (
unsigned int q = 0; q < dofs_per_component; ++q)
2391 out[q] *= inv_coefficient[q];
2393 for (
unsigned int di = 0; di < n_comp_inner; ++di)
2395 Number *out_ = out + di * dofs_per_component;
2397 evaluator.template hessians<2, false, false>(out_, out_);
2399 evaluator.template hessians<1, false, false>(out_, out_);
2400 evaluator.template hessians<0, false, false>(out_, out_);
2403 in += dofs_per_component;
2404 out += dofs_per_component;
2405 inv_coefficient += shift_coefficient;
2412 template <
int n_components>
2414 vmult(
const Number *inverse_coefficients,
2417 const unsigned int dofs_per_component,
2418 const unsigned int n_given_components = 0)
2420 const unsigned int n_desired_components =
2421 (n_components > -1) ? n_components : n_given_components;
2423 std::array<Number, dim + 2> tmp = {};
2424 Assert(n_desired_components <= dim + 2,
2426 "Number of components larger than dim+2 not supported."));
2428 for (
unsigned int d = 0;
d < n_desired_components; ++
d)
2429 tmp[
d] = src[
d * dofs_per_component];
2431 for (
unsigned int d1 = 0; d1 < n_desired_components; ++d1)
2433 const Number *inv_coeff_row =
2434 &inverse_coefficients[d1 * n_desired_components];
2435 Number
sum = inv_coeff_row[0] * tmp[0];
2436 for (
unsigned int d2 = 1; d2 < n_desired_components; ++d2)
2437 sum += inv_coeff_row[d2] * tmp[d2];
2438 dst[d1 * dofs_per_component] =
sum;
2451 template <
int dim,
typename Number>
2454 template <
int fe_degree,
int n_q_po
ints_1d>
2456 run(
const unsigned int n_desired_components,
2458 const Number *in_array,
2461 static const bool do_inplace =
2462 fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
2468 const auto &inverse_shape =
2473 const std::size_t dofs_per_component =
2476 const std::size_t n_q_points = do_inplace ?
2494 for (
unsigned int d = 0;
d < n_desired_components; ++
d)
2496 const Number *in = in_array +
d * n_q_points;
2497 Number *out = out_array +
d * dofs_per_component;
2500 auto *temp_2 = do_inplace ?
2502 (temp_1 +
std::max(n_q_points, dofs_per_component));
2506 evaluator.template hessians<2, false, false>(in, temp_1);
2507 evaluator.template hessians<1, false, false>(temp_1, temp_2);
2508 evaluator.template hessians<0, false, false>(temp_2, out);
2512 evaluator.template hessians<1, false, false>(in, temp_1);
2513 evaluator.template hessians<0, false, false>(temp_1, out);
2516 evaluator.template hessians<0, false, false>(in, out);
value_type * data() const noexcept
Number JxW(const unsigned int q_point) const
const Number * begin_values() const
const Number * begin_hessians() const
Number shape_info_number_type
const Number * begin_gradients() const
ArrayView< Number > get_scratch_data() const
const ShapeInfoType & get_shape_info() const
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ tensor_symmetric_no_collocation
@ tensor_symmetric_collocation
@ tensor_symmetric_plus_dg0
EvaluationFlags
The EvaluationFlags enum.
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
constexpr T pow(const T base, const int iexp)
void evaluate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
void integrate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, Number *values, const Number *gradients, const bool add_into_values_array)
void evaluate_hessians_slow(const unsigned int n_components, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
void integrate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
void evaluate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, const Number *values, Number *gradients)
void integrate_hessians_slow(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, Number *values_dofs, const bool add_into_values_array)
static const unsigned int invalid_unsigned_int
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const ArrayView< const Number > &inverse_coefficients, const bool dyadic_coefficients, const Number *in_array, Number *out_array)
static void vmult(const Number *inverse_coefficients, const Number *src, Number *dst, const unsigned int dofs_per_component, const unsigned int n_given_components=0)
const Number2 * shape_values
static void do_backward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const bool add_into_result, Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static void do_forward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static void do_mass(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, OtherNumber *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array, std::bool_constant< true >)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, OtherNumber *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array_in=false)
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array, std::bool_constant< false >)
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
static const EvaluatorVariant variant
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
EvaluatorTensorProduct< variant, dim, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > *univariate_shape_data)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval)
unsigned int dofs_per_component_on_cell
std::vector< UnivariateShapeData< Number > > data
AlignedVector< Number > shape_values
AlignedVector< Number > shape_hessians_collocation
AlignedVector< Number > shape_values_eo
AlignedVector< Number > shape_hessians_eo
AlignedVector< Number > shape_gradients_collocation_eo
unsigned int n_q_points_1d
AlignedVector< Number > shape_gradients_eo
AlignedVector< Number > shape_hessians
AlignedVector< Number > shape_gradients
AlignedVector< Number > shape_gradients_collocation