From cbd918fa6cd6e09f2d5cd0a62e9a616d561b5405 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 17 Jan 2024 17:18:20 -0700 Subject: [PATCH] Use p[index] instead of p(index) for class Point. --- examples/step-12/step-12.cc | 2 +- examples/step-12b/step-12b.cc | 2 +- examples/step-13/step-13.cc | 20 +- examples/step-16b/step-16b.cc | 2 +- examples/step-17/step-17.cc | 4 +- examples/step-30/step-30.cc | 17 +- examples/step-32/step-32.cc | 4 +- examples/step-35/step-35.cc | 6 +- examples/step-38/step-38.cc | 42 +- examples/step-4/step-4.cc | 2 +- examples/step-41/step-41.cc | 6 +- examples/step-42/step-42.cc | 22 +- examples/step-49/step-49.cc | 6 +- examples/step-51/step-51.cc | 4 +- examples/step-52/step-52.cc | 2 +- examples/step-53/step-53.cc | 8 +- examples/step-56/step-56.cc | 30 +- examples/step-57/step-57.cc | 8 +- examples/step-7/step-7.cc | 4 +- examples/step-78/step-78.cc | 18 +- examples/step-79/step-79.cc | 14 +- examples/step-8/step-8.cc | 4 +- examples/step-82/step-82.cc | 130 ++--- include/deal.II/base/data_out_base.h | 6 +- include/deal.II/base/point.h | 2 +- include/deal.II/base/polynomial_space.h | 2 +- .../base/polynomials_rannacher_turek.h | 16 +- .../deal.II/base/tensor_product_polynomials.h | 4 +- .../base/tensor_product_polynomials_bubbles.h | 18 +- include/deal.II/fe/fe_tools.templates.h | 4 +- include/deal.II/fe/mapping_q_internal.h | 28 +- .../vector_tools_boundary.templates.h | 4 +- include/deal.II/particles/particle_accessor.h | 4 +- source/base/data_out_base.cc | 170 +++--- source/base/flow_function.cc | 44 +- source/base/function_lib.cc | 532 +++++++++--------- source/base/function_restriction.cc | 10 +- source/base/geometric_utilities.cc | 4 +- source/base/mu_parser_internal.cc | 4 +- source/base/polynomial.cc | 6 +- source/base/polynomial_space.cc | 8 +- source/base/polynomials_abf.cc | 2 +- source/base/polynomials_adini.cc | 12 +- source/base/polynomials_bdm.cc | 60 +- source/base/polynomials_nedelec.cc | 16 +- source/base/polynomials_rannacher_turek.cc | 32 +- source/base/polynomials_raviart_thomas.cc | 2 +- source/base/polynomials_rt_bubbles.cc | 14 +- source/base/qprojector.cc | 52 +- source/base/quadrature.cc | 30 +- source/base/tensor_product_polynomials.cc | 12 +- .../tensor_product_polynomials_bubbles.cc | 22 +- source/fe/fe_dgp_monomial.cc | 12 +- source/fe/fe_dgq.cc | 6 +- source/fe/fe_face.cc | 6 +- source/fe/fe_nedelec.cc | 404 ++++++------- source/fe/fe_p1nc.cc | 40 +- source/fe/fe_poly_tensor.cc | 2 +- source/fe/fe_q.cc | 2 +- source/fe/fe_q_base.cc | 12 +- source/fe/fe_q_bubbles.cc | 14 +- source/fe/fe_q_dg0.cc | 12 +- source/fe/fe_q_hierarchical.cc | 36 +- source/fe/fe_rannacher_turek.cc | 8 +- source/fe/mapping.cc | 10 +- source/fe/mapping_cartesian.cc | 42 +- source/fe/mapping_q.cc | 4 +- source/fe/mapping_q_eulerian.cc | 2 +- source/grid/grid_generator.cc | 194 +++---- source/grid/grid_generator_pipe_junction.cc | 6 +- source/grid/grid_in.cc | 62 +- source/grid/grid_out.cc | 52 +- source/grid/grid_tools.cc | 12 +- source/grid/grid_tools_dof_handlers.cc | 4 +- source/grid/grid_tools_nontemplates.cc | 76 +-- source/grid/manifold_lib.cc | 28 +- source/grid/tria_accessor.cc | 64 +-- source/non_matching/quadrature_generator.cc | 2 +- source/numerics/data_out_rotation.cc | 14 +- source/numerics/data_out_stack.cc | 32 +- source/particles/particle.cc | 8 +- 81 files changed, 1321 insertions(+), 1322 deletions(-) diff --git a/examples/step-12/step-12.cc b/examples/step-12/step-12.cc index b365f5b0bd..ac35aad427 100644 --- a/examples/step-12/step-12.cc +++ b/examples/step-12/step-12.cc @@ -101,7 +101,7 @@ namespace Step12 for (unsigned int i = 0; i < values.size(); ++i) { - if (points[i](0) < 0.5) + if (points[i][0] < 0.5) values[i] = 1.; else values[i] = 0.; diff --git a/examples/step-12b/step-12b.cc b/examples/step-12b/step-12b.cc index 9f5d66b513..4f5c9e6d67 100644 --- a/examples/step-12b/step-12b.cc +++ b/examples/step-12b/step-12b.cc @@ -106,7 +106,7 @@ namespace Step12 for (unsigned int i = 0; i < values.size(); ++i) { - if (points[i](0) < 0.5) + if (points[i][0] < 0.5) values[i] = 1.; else values[i] = 0.; diff --git a/examples/step-13/step-13.cc b/examples/step-13/step-13.cc index 805ade33da..4dcc4731ac 100644 --- a/examples/step-13/step-13.cc +++ b/examples/step-13/step-13.cc @@ -1224,9 +1224,9 @@ namespace Step13 { (void)component; AssertIndexRange(component, 1); - double q = p(0); + double q = p[0]; for (unsigned int i = 1; i < dim; ++i) - q += std::sin(10 * p(i) + 5 * p(0) * p(0)); + q += std::sin(10 * p[i] + 5 * p[0] * p[0]); const double exponential = std::exp(q); return exponential; } @@ -1248,19 +1248,19 @@ namespace Step13 { (void)component; AssertIndexRange(component, 1); - double q = p(0); + double q = p[0]; for (unsigned int i = 1; i < dim; ++i) - q += std::sin(10 * p(i) + 5 * p(0) * p(0)); + q += std::sin(10 * p[i] + 5 * p[0] * p[0]); const double u = std::exp(q); double t1 = 1, t2 = 0, t3 = 0; for (unsigned int i = 1; i < dim; ++i) { - t1 += std::cos(10 * p(i) + 5 * p(0) * p(0)) * 10 * p(0); - t2 += 10 * std::cos(10 * p(i) + 5 * p(0) * p(0)) - - 100 * std::sin(10 * p(i) + 5 * p(0) * p(0)) * p(0) * p(0); - t3 += 100 * std::cos(10 * p(i) + 5 * p(0) * p(0)) * - std::cos(10 * p(i) + 5 * p(0) * p(0)) - - 100 * std::sin(10 * p(i) + 5 * p(0) * p(0)); + t1 += std::cos(10 * p[i] + 5 * p[0] * p[0]) * 10 * p[0]; + t2 += 10 * std::cos(10 * p[i] + 5 * p[0] * p[0]) - + 100 * std::sin(10 * p[i] + 5 * p[0] * p[0]) * p[0] * p[0]; + t3 += 100 * std::cos(10 * p[i] + 5 * p[0] * p[0]) * + std::cos(10 * p[i] + 5 * p[0] * p[0]) - + 100 * std::sin(10 * p[i] + 5 * p[0] * p[0]); }; t1 = t1 * t1; diff --git a/examples/step-16b/step-16b.cc b/examples/step-16b/step-16b.cc index 858a539120..2755514570 100644 --- a/examples/step-16b/step-16b.cc +++ b/examples/step-16b/step-16b.cc @@ -146,7 +146,7 @@ namespace Step16 MeshWorker::IntegrationInfo &info) const { AssertDimension(dinfo.n_matrices(), 1); - const double coefficient = (dinfo.cell->center()(0) > 0.) ? .1 : 1.; + const double coefficient = (dinfo.cell->center()[0] > 0.) ? .1 : 1.; LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0, false).matrix, info.fe_values(0), diff --git a/examples/step-17/step-17.cc b/examples/step-17/step-17.cc index ed620ba785..8c549b6afe 100644 --- a/examples/step-17/step-17.cc +++ b/examples/step-17/step-17.cc @@ -187,8 +187,8 @@ namespace Step17 Assert(dim >= 2, ExcInternalError()); Point point_1, point_2; - point_1(0) = 0.5; - point_2(0) = -0.5; + point_1[0] = 0.5; + point_2[0] = -0.5; if (((p - point_1).norm_square() < 0.2 * 0.2) || ((p - point_2).norm_square() < 0.2 * 0.2)) diff --git a/examples/step-30/step-30.cc b/examples/step-30/step-30.cc index d364147c22..7cb4b79a94 100644 --- a/examples/step-30/step-30.cc +++ b/examples/step-30/step-30.cc @@ -81,7 +81,7 @@ namespace Step30 for (unsigned int i = 0; i < values.size(); ++i) { - if (points[i](0) < 0.5) + if (points[i][0] < 0.5) values[i] = 1.; else values[i] = 0.; @@ -111,15 +111,15 @@ namespace Step30 for (unsigned int i = 0; i < points.size(); ++i) { - if (points[i](0) > 0) + if (points[i][0] > 0) { - values[i](0) = -points[i](1); - values[i](1) = points[i](0); + values[i][0] = -points[i][1]; + values[i][1] = points[i][0]; } else { values[i] = Point(); - values[i](0) = -points[i](1); + values[i][0] = -points[i][1]; } } } @@ -873,7 +873,7 @@ namespace Step30 double sum_of_average_jumps = 0.; for (unsigned int i = 0; i < dim; ++i) { - average_jumps[i] = jump(i) / area(i); + average_jumps[i] = jump[i] / area[i]; sum_of_average_jumps += average_jumps[i]; } @@ -946,10 +946,9 @@ namespace Step30 { // Create the rectangular domain. Point p1, p2; - p1(0) = 0; - p1(0) = -1; + p1[0] = -1.; for (unsigned int i = 0; i < dim; ++i) - p2(i) = 1.; + p2[i] = 1.; // Adjust the number of cells in different directions to obtain // completely isotropic cells for the original mesh. std::vector repetitions(dim, 1); diff --git a/examples/step-32/step-32.cc b/examples/step-32/step-32.cc index 32aeccbaa9..d5be30c8b6 100644 --- a/examples/step-32/step-32.cc +++ b/examples/step-32/step-32.cc @@ -169,8 +169,8 @@ namespace Step32 const double s = (r - R0) / h; const double q = - (dim == 3) ? std::max(0.0, cos(numbers::PI * abs(p(2) / R1))) : 1.0; - const double phi = std::atan2(p(0), p(1)); + (dim == 3) ? std::max(0.0, cos(numbers::PI * abs(p[2] / R1))) : 1.0; + const double phi = std::atan2(p[0], p[1]); const double tau = s + 0.2 * s * (1 - s) * std::sin(6 * phi) * q; return T0 * (1.0 - tau) + T1 * tau; diff --git a/examples/step-35/step-35.cc b/examples/step-35/step-35.cc index 3b2a8b6a6b..5b2e4755bc 100644 --- a/examples/step-35/step-35.cc +++ b/examples/step-35/step-35.cc @@ -364,7 +364,7 @@ namespace Step35 { const double Um = 1.5; const double H = 4.1; - return 4. * Um * p(1) * (H - p(1)) / (H * H); + return 4. * Um * p[1] * (H - p[1]) / (H * H); } else return 0.; @@ -398,7 +398,7 @@ namespace Step35 { (void)component; AssertIndexRange(component, 1); - return 25. - p(0); + return 25. - p[0]; } template @@ -1156,7 +1156,7 @@ namespace Step35 { scratch.fe_val.get_function_values(u_star[d], scratch.u_star_tmp); for (unsigned int q = 0; q < scratch.nqp; ++q) - scratch.u_star_local[q](d) = scratch.u_star_tmp[q]; + scratch.u_star_local[q][d] = scratch.u_star_tmp[q]; } for (unsigned int d = 0; d < dim; ++d) diff --git a/examples/step-38/step-38.cc b/examples/step-38/step-38.cc index 20905ef4aa..7e038b0b82 100644 --- a/examples/step-38/step-38.cc +++ b/examples/step-38/step-38.cc @@ -153,7 +153,7 @@ namespace Step38 template <> double Solution<2>::value(const Point<2> &p, const unsigned int) const { - return (-2. * p(0) * p(1)); + return (-2. * p[0] * p[1]); } @@ -162,8 +162,8 @@ namespace Step38 const unsigned int) const { Tensor<1, 2> return_value; - return_value[0] = -2. * p(1) * (1 - 2. * p(0) * p(0)); - return_value[1] = -2. * p(0) * (1 - 2. * p(1) * p(1)); + return_value[0] = -2. * p[1] * (1 - 2. * p[0] * p[0]); + return_value[1] = -2. * p[0] * (1 - 2. * p[1] * p[1]); return return_value; } @@ -172,8 +172,8 @@ namespace Step38 template <> double Solution<3>::value(const Point<3> &p, const unsigned int) const { - return (std::sin(numbers::PI * p(0)) * std::cos(numbers::PI * p(1)) * - exp(p(2))); + return (std::sin(numbers::PI * p[0]) * std::cos(numbers::PI * p[1]) * + exp(p[2])); } @@ -185,9 +185,9 @@ namespace Step38 Tensor<1, 3> return_value; - return_value[0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); - return_value[1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); - return_value[2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + return_value[0] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); + return_value[1] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]); + return_value[2] = sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); return return_value; } @@ -206,7 +206,7 @@ namespace Step38 double RightHandSide<2>::value(const Point<2> &p, const unsigned int /*component*/) const { - return (-8. * p(0) * p(1)); + return (-8. * p[0] * p[1]); } @@ -218,23 +218,23 @@ namespace Step38 Tensor<2, 3> hessian; - hessian[0][0] = -PI * PI * sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); - hessian[1][1] = -PI * PI * sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); - hessian[2][2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + hessian[0][0] = -PI * PI * sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); + hessian[1][1] = -PI * PI * sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); + hessian[2][2] = sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); - hessian[0][1] = -PI * PI * cos(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); - hessian[1][0] = -PI * PI * cos(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + hessian[0][1] = -PI * PI * cos(PI * p[0]) * sin(PI * p[1]) * exp(p[2]); + hessian[1][0] = -PI * PI * cos(PI * p[0]) * sin(PI * p[1]) * exp(p[2]); - hessian[0][2] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); - hessian[2][0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + hessian[0][2] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); + hessian[2][0] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); - hessian[1][2] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); - hessian[2][1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + hessian[1][2] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]); + hessian[2][1] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]); Tensor<1, 3> gradient; - gradient[0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); - gradient[1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); - gradient[2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + gradient[0] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); + gradient[1] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]); + gradient[2] = sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]); Point<3> normal = p; normal /= p.norm(); diff --git a/examples/step-4/step-4.cc b/examples/step-4/step-4.cc index cdfbcd846b..e24dd1f593 100644 --- a/examples/step-4/step-4.cc +++ b/examples/step-4/step-4.cc @@ -190,7 +190,7 @@ double RightHandSide::value(const Point &p, { double return_value = 0.0; for (unsigned int i = 0; i < dim; ++i) - return_value += 4.0 * std::pow(p(i), 4.0); + return_value += 4.0 * std::pow(p[i], 4.0); return return_value; } diff --git a/examples/step-41/step-41.cc b/examples/step-41/step-41.cc index d3416259cf..6d61c69044 100644 --- a/examples/step-41/step-41.cc +++ b/examples/step-41/step-41.cc @@ -157,11 +157,11 @@ namespace Step41 (void)component; Assert(component == 0, ExcIndexRange(component, 0, 1)); - if (p(0) < -0.5) + if (p[0] < -0.5) return -0.2; - else if (p(0) >= -0.5 && p(0) < 0.0) + else if (p[0] >= -0.5 && p[0] < 0.0) return -0.4; - else if (p(0) >= 0.0 && p(0) < 0.5) + else if (p[0] >= 0.0 && p[0] < 0.5) return -0.6; else return -0.8; diff --git a/examples/step-42/step-42.cc b/examples/step-42/step-42.cc index 585f72a243..4b5c465f33 100644 --- a/examples/step-42/step-42.cc +++ b/examples/step-42/step-42.cc @@ -372,14 +372,14 @@ namespace Step42 const unsigned int component) const { if (component == 0) - return p(0); + return p[0]; else if (component == 1) - return p(1); + return p[1]; else if (component == 2) { - if ((p(0) - 0.5) * (p(0) - 0.5) + (p(1) - 0.5) * (p(1) - 0.5) < 0.36) - return (-std::sqrt(0.36 - (p(0) - 0.5) * (p(0) - 0.5) - - (p(1) - 0.5) * (p(1) - 0.5)) + + if ((p[0] - 0.5) * (p[0] - 0.5) + (p[1] - 0.5) * (p[1] - 0.5) < 0.36) + return (-std::sqrt(0.36 - (p[0] - 0.5) * (p[0] - 0.5) - + (p[1] - 0.5) * (p[1] - 0.5)) + z_surface + 0.59); else return 1000; @@ -550,13 +550,13 @@ namespace Step42 const unsigned int component) const { if (component == 0) - return p(0); + return p[0]; if (component == 1) - return p(1); + return p[1]; else if (component == 2) { - if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0) - return z_surface + 0.999 - input_obstacle.get_value(p(0), p(1)); + if (p[0] >= 0.0 && p[0] <= 1.0 && p[1] >= 0.0 && p[1] <= 1.0) + return z_surface + 0.999 - input_obstacle.get_value(p[0], p[1]); } Assert(false, ExcNotImplemented()); @@ -881,7 +881,7 @@ namespace Step42 // indicator one. Point<3> rotate_half_sphere(const Point<3> &in) { - return {in(2), in(1), -in(0)}; + return {in[2], in[1], -in[0]}; } template @@ -1267,7 +1267,7 @@ namespace Step42 obstacle->value(this_support_point, 2); const double solution_here = solution(index_z); const double undeformed_gap = - obstacle_value - this_support_point(2); + obstacle_value - this_support_point[2]; const double c = 100.0 * e_modulus; if ((lambda(index_z) / diff --git a/examples/step-49/step-49.cc b/examples/step-49/step-49.cc index e2cb21552a..db30d67822 100644 --- a/examples/step-49/step-49.cc +++ b/examples/step-49/step-49.cc @@ -209,8 +209,8 @@ void grid_3() for (const auto i : cell->vertex_indices()) { Point<2> &v = cell->vertex(i); - if (std::abs(v(1) - 1.0) < 1e-5) - v(1) += 0.5; + if (std::abs(v[1] - 1.0) < 1e-5) + v[1] += 0.5; } } @@ -315,7 +315,7 @@ struct Grid6Func Point<2> operator()(const Point<2> &in) const { - return {in(0), trans(in(1))}; + return {in[0], trans(in[1])}; } }; diff --git a/examples/step-51/step-51.cc b/examples/step-51/step-51.cc index a631c8c443..2dbdf85b14 100644 --- a/examples/step-51/step-51.cc +++ b/examples/step-51/step-51.cc @@ -1344,8 +1344,8 @@ namespace Step51 for (const auto &cell : triangulation.cell_iterators()) for (const auto &face : cell->face_iterators()) if (face->at_boundary()) - if ((std::fabs(face->center()(0) - (-1)) < 1e-12) || - (std::fabs(face->center()(1) - (-1)) < 1e-12)) + if ((std::fabs(face->center()[0] - (-1)) < 1e-12) || + (std::fabs(face->center()[1] - (-1)) < 1e-12)) face->set_boundary_id(1); } diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc index 94998d7d50..96d247d29b 100644 --- a/examples/step-52/step-52.cc +++ b/examples/step-52/step-52.cc @@ -260,7 +260,7 @@ namespace Step52 const double intensity = 10.; const double frequency = numbers::PI / 10.; const double b = 5.; - const double x = point(0); + const double x = point[0]; return intensity * (frequency * std::cos(frequency * time) * (b * x - x * x) + diff --git a/examples/step-53/step-53.cc b/examples/step-53/step-53.cc index d24da54d13..163fbfe2f5 100644 --- a/examples/step-53/step-53.cc +++ b/examples/step-53/step-53.cc @@ -259,11 +259,11 @@ namespace Step53 { const double b = std::sqrt(R * R * (1 - ellipticity * ellipticity)); const double ep = std::sqrt((R * R - b * b) / (b * b)); - const double p = std::sqrt(x(0) * x(0) + x(1) * x(1)); - const double th = std::atan2(R * x(2), b * p); - const double phi = std::atan2(x(1), x(0)); + const double p = std::sqrt(x[0] * x[0] + x[1] * x[1]); + const double th = std::atan2(R * x[2], b * p); + const double phi = std::atan2(x[1], x[0]); const double theta = - std::atan2(x(2) + ep * ep * b * Utilities::fixed_power<3>(std::sin(th)), + std::atan2(x[2] + ep * ep * b * Utilities::fixed_power<3>(std::sin(th)), (p - (ellipticity * ellipticity * R * Utilities::fixed_power<3>(std::cos(th))))); const double R_bar = diff --git a/examples/step-56/step-56.cc b/examples/step-56/step-56.cc index 188f297e8d..2ae71956ef 100644 --- a/examples/step-56/step-56.cc +++ b/examples/step-56/step-56.cc @@ -116,8 +116,8 @@ namespace Step56 Assert(component <= 2 + 1, ExcIndexRange(component, 0, 2 + 1)); using numbers::PI; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; if (component == 0) return sin(PI * x); @@ -136,9 +136,9 @@ namespace Step56 Assert(component <= 3 + 1, ExcIndexRange(component, 0, 3 + 1)); using numbers::PI; - const double x = p(0); - const double y = p(1); - const double z = p(2); + const double x = p[0]; + const double y = p[1]; + const double z = p[2]; if (component == 0) return 2.0 * sin(PI * x); @@ -160,8 +160,8 @@ namespace Step56 Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); using numbers::PI; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; Tensor<1, 2> return_value; if (component == 0) @@ -190,9 +190,9 @@ namespace Step56 Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); using numbers::PI; - const double x = p(0); - const double y = p(1); - const double z = p(2); + const double x = p[0]; + const double y = p[1]; + const double z = p[2]; Tensor<1, 3> return_value; if (component == 0) @@ -243,8 +243,8 @@ namespace Step56 Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); using numbers::PI; - double x = p(0); - double y = p(1); + const double x = p[0]; + const double y = p[1]; if (component == 0) return PI * PI * sin(PI * x) + PI * cos(PI * x) * cos(PI * y); if (component == 1) @@ -262,9 +262,9 @@ namespace Step56 Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); using numbers::PI; - double x = p(0); - double y = p(1); - double z = p(2); + const double x = p[0]; + const double y = p[1]; + const double z = p[2]; if (component == 0) return 2 * PI * PI * sin(PI * x) + PI * cos(PI * x) * cos(PI * y) * sin(PI * z); diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc index e154a178a5..58086ee94d 100644 --- a/examples/step-57/step-57.cc +++ b/examples/step-57/step-57.cc @@ -788,18 +788,18 @@ namespace Step57 f << "# y u_x u_y" << std::endl; Point p; - p(0) = 0.5; - p(1) = 0.5; + p[0] = 0.5; + p[1] = 0.5; f << std::scientific; for (unsigned int i = 0; i <= 100; ++i) { - p(dim - 1) = i / 100.0; + p[dim - 1] = i / 100.0; Vector tmp_vector(dim + 1); VectorTools::point_value(dof_handler, present_solution, p, tmp_vector); - f << p(dim - 1); + f << p[dim - 1]; for (int j = 0; j < dim; ++j) f << ' ' << tmp_vector(j); diff --git a/examples/step-7/step-7.cc b/examples/step-7/step-7.cc index 55cf0498a7..0e485cf8dd 100644 --- a/examples/step-7/step-7.cc +++ b/examples/step-7/step-7.cc @@ -967,8 +967,8 @@ namespace Step7 for (const auto &face : cell->face_iterators()) { const auto center = face->center(); - if ((std::fabs(center(0) - (-1.0)) < 1e-12) || - (std::fabs(center(1) - (-1.0)) < 1e-12)) + if ((std::fabs(center[0] - (-1.0)) < 1e-12) || + (std::fabs(center[1] - (-1.0)) < 1e-12)) face->set_boundary_id(1); } } diff --git a/examples/step-78/step-78.cc b/examples/step-78/step-78.cc index efd6b6f0a6..383d1433cf 100644 --- a/examples/step-78/step-78.cc +++ b/examples/step-78/step-78.cc @@ -101,7 +101,7 @@ namespace BlackScholesSolver double Solution::value(const Point &p, const unsigned int component) const { - return -Utilities::fixed_power<2, double>(p(component)) - + return -Utilities::fixed_power<2, double>(p[component]) - Utilities::fixed_power<2, double>(this->get_time()) + 6; } @@ -110,7 +110,7 @@ namespace BlackScholesSolver Tensor<1, dim> Solution::gradient(const Point &p, const unsigned int component) const { - return Point(-2 * p(component)); + return Point(-2 * p[component]); } @@ -148,9 +148,9 @@ namespace BlackScholesSolver const unsigned int component) const { #ifdef MMS - return -Utilities::fixed_power<2, double>(p(component)) + 6; + return -Utilities::fixed_power<2, double>(p[component]) + 6; #else - return std::max(p(component) - strike_price, 0.); + return std::max(p[component] - strike_price, 0.); #endif } @@ -208,10 +208,10 @@ namespace BlackScholesSolver const unsigned int component) const { #ifdef MMS - return -Utilities::fixed_power<2, double>(p(component)) - + return -Utilities::fixed_power<2, double>(p[component]) - Utilities::fixed_power<2, double>(this->get_time()) + 6; #else - return (p(component) - strike_price) * + return (p[component] - strike_price) * exp((-interest_rate) * (this->get_time())); #endif } @@ -248,10 +248,10 @@ namespace BlackScholesSolver { #ifdef MMS return 2 * (this->get_time()) - - Utilities::fixed_power<2, double>(asset_volatility * p(component)) - - 2 * interest_rate * Utilities::fixed_power<2, double>(p(component)) - + Utilities::fixed_power<2, double>(asset_volatility * p[component]) - + 2 * interest_rate * Utilities::fixed_power<2, double>(p[component]) - interest_rate * - (-Utilities::fixed_power<2, double>(p(component)) - + (-Utilities::fixed_power<2, double>(p[component]) - Utilities::fixed_power<2, double>(this->get_time()) + 6); #else (void)p; diff --git a/examples/step-79/step-79.cc b/examples/step-79/step-79.cc index 7d4a69c333..b014fec61c 100644 --- a/examples/step-79/step-79.cc +++ b/examples/step-79/step-79.cc @@ -332,9 +332,9 @@ namespace SAND if (face->at_boundary()) { const auto center = face->center(); - if (std::fabs(center(1) - 1) < 1e-12) + if (std::fabs(center[1] - 1) < 1e-12) { - if ((std::fabs(center(0) - 3) < .3)) + if ((std::fabs(center[0] - 3) < .3)) face->set_boundary_id(BoundaryIds::down_force); else face->set_boundary_id(BoundaryIds::no_force); @@ -374,14 +374,14 @@ namespace SAND // boundary, and if it is whether one of its // vertices might be the bottom left or bottom // right vertex: - if (std::fabs(center(1) - 0) < 1e-12) + if (std::fabs(center[1] - 0) < 1e-12) { for (const auto vertex_number : cell->vertex_indices()) { const auto vert = cell->vertex(vertex_number); - if (std::fabs(vert(0) - 0) < 1e-12 && - std::fabs(vert(1) - 0) < 1e-12) + if (std::fabs(vert[0] - 0) < 1e-12 && + std::fabs(vert[1] - 0) < 1e-12) { types::global_dof_index x_displacement = cell->vertex_dof_index(vertex_number, 0); @@ -398,8 +398,8 @@ namespace SAND boundary_values[y_displacement_multiplier] = 0; } - else if (std::fabs(vert(0) - 6) < 1e-12 && - std::fabs(vert(1) - 0) < 1e-12) + else if (std::fabs(vert[0] - 6) < 1e-12 && + std::fabs(vert[1] - 0) < 1e-12) { types::global_dof_index y_displacement = cell->vertex_dof_index(vertex_number, 1); diff --git a/examples/step-8/step-8.cc b/examples/step-8/step-8.cc index cb3d06b5fa..d7b109a80a 100644 --- a/examples/step-8/step-8.cc +++ b/examples/step-8/step-8.cc @@ -150,8 +150,8 @@ namespace Step8 // these areas. Note that upon construction of the Point objects, all // components are set to zero. Point point_1, point_2; - point_1(0) = 0.5; - point_2(0) = -0.5; + point_1[0] = 0.5; + point_2[0] = -0.5; for (unsigned int point_n = 0; point_n < points.size(); ++point_n) { diff --git a/examples/step-82/step-82.cc b/examples/step-82/step-82.cc index e05666259a..7b663d03bf 100644 --- a/examples/step-82/step-82.cc +++ b/examples/step-82/step-82.cc @@ -169,28 +169,28 @@ namespace Step82 if (dim == 2) { - return_value = 24.0 * Utilities::fixed_power<2>(p(1) * (1.0 - p(1))) + - +24.0 * Utilities::fixed_power<2>(p(0) * (1.0 - p(0))) + - 2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * - (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)); + return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) + + +24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0])) + + 2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) * + (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]); } else if (dim == 3) { - return_value = 24.0 * Utilities::fixed_power<2>(p(1) * (1.0 - p(1)) * - p(2) * (1.0 - p(2))) + - 24.0 * Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * - p(2) * (1.0 - p(2))) + - 24.0 * Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * - p(1) * (1.0 - p(1))) + - 2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * - (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * - Utilities::fixed_power<2>(p(2) * (1.0 - p(2))) + - 2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * - (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) * - Utilities::fixed_power<2>(p(1) * (1.0 - p(1))) + - 2.0 * (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * - (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0))); + return_value = 24.0 * Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * + p[2] * (1.0 - p[2])) + + 24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * + p[2] * (1.0 - p[2])) + + 24.0 * Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * + p[1] * (1.0 - p[1])) + + 2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) * + (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) * + Utilities::fixed_power<2>(p[2] * (1.0 - p[2])) + + 2.0 * (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) * + (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) * + Utilities::fixed_power<2>(p[1] * (1.0 - p[1])) + + 2.0 * (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) * + (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0])); } else Assert(false, ExcNotImplemented()); @@ -233,12 +233,12 @@ namespace Step82 if (dim == 2) { return_value = - Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1))); + Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1])); } else if (dim == 3) { return_value = Utilities::fixed_power<2>( - p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2))); + p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2])); } else Assert(false, ExcNotImplemented()); @@ -258,28 +258,28 @@ namespace Step82 if (dim == 2) { return_gradient[0] = - (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) + - 4.0 * Utilities::fixed_power<3>(p(0))) * - Utilities::fixed_power<2>(p(1) * (1.0 - p(1))); + (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) + + 4.0 * Utilities::fixed_power<3>(p[0])) * + Utilities::fixed_power<2>(p[1] * (1.0 - p[1])); return_gradient[1] = - (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) + - 4.0 * Utilities::fixed_power<3>(p(1))) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0))); + (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) + + 4.0 * Utilities::fixed_power<3>(p[1])) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0])); } else if (dim == 3) { return_gradient[0] = - (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) + - 4.0 * Utilities::fixed_power<3>(p(0))) * - Utilities::fixed_power<2>(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2))); + (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) + + 4.0 * Utilities::fixed_power<3>(p[0])) * + Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2])); return_gradient[1] = - (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) + - 4.0 * Utilities::fixed_power<3>(p(1))) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2))); + (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) + + 4.0 * Utilities::fixed_power<3>(p[1])) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2])); return_gradient[2] = - (2.0 * p(2) - 6.0 * Utilities::fixed_power<2>(p(2)) + - 4.0 * Utilities::fixed_power<3>(p(2))) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1))); + (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) + + 4.0 * Utilities::fixed_power<3>(p[2])) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1])); } else Assert(false, ExcNotImplemented()); @@ -298,45 +298,45 @@ namespace Step82 if (dim == 2) { - return_hessian[0][0] = (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * - Utilities::fixed_power<2>(p(1) * (1.0 - p(1))); + return_hessian[0][0] = (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) * + Utilities::fixed_power<2>(p[1] * (1.0 - p[1])); return_hessian[0][1] = - (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) + - 4.0 * Utilities::fixed_power<3>(p(0))) * - (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) + - 4.0 * Utilities::fixed_power<3>(p(1))); - return_hessian[1][1] = (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0))); + (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) + + 4.0 * Utilities::fixed_power<3>(p[0])) * + (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) + + 4.0 * Utilities::fixed_power<3>(p[1])); + return_hessian[1][1] = (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0])); } else if (dim == 3) { return_hessian[0][0] = - (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) * - Utilities::fixed_power<2>(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2))); + (2.0 - 12.0 * p[0] + 12.0 * p[0] * p[0]) * + Utilities::fixed_power<2>(p[1] * (1.0 - p[1]) * p[2] * (1.0 - p[2])); return_hessian[0][1] = - (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) + - 4.0 * Utilities::fixed_power<3>(p(0))) * - (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) + - 4.0 * Utilities::fixed_power<3>(p(1))) * - Utilities::fixed_power<2>(p(2) * (1.0 - p(2))); + (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) + + 4.0 * Utilities::fixed_power<3>(p[0])) * + (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) + + 4.0 * Utilities::fixed_power<3>(p[1])) * + Utilities::fixed_power<2>(p[2] * (1.0 - p[2])); return_hessian[0][2] = - (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) + - 4.0 * Utilities::fixed_power<3>(p(0))) * - (2.0 * p(2) - 6.0 * Utilities::fixed_power<2>(p(2)) + - 4.0 * Utilities::fixed_power<3>(p(2))) * - Utilities::fixed_power<2>(p(1) * (1.0 - p(1))); + (2.0 * p[0] - 6.0 * Utilities::fixed_power<2>(p[0]) + + 4.0 * Utilities::fixed_power<3>(p[0])) * + (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) + + 4.0 * Utilities::fixed_power<3>(p[2])) * + Utilities::fixed_power<2>(p[1] * (1.0 - p[1])); return_hessian[1][1] = - (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2))); + (2.0 - 12.0 * p[1] + 12.0 * p[1] * p[1]) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[2] * (1.0 - p[2])); return_hessian[1][2] = - (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) + - 4.0 * Utilities::fixed_power<3>(p(1))) * - (2.0 * p(2) - 6.0 * Utilities::fixed_power<2>(p(2)) + - 4.0 * Utilities::fixed_power<3>(p(2))) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0))); + (2.0 * p[1] - 6.0 * Utilities::fixed_power<2>(p[1]) + + 4.0 * Utilities::fixed_power<3>(p[1])) * + (2.0 * p[2] - 6.0 * Utilities::fixed_power<2>(p[2]) + + 4.0 * Utilities::fixed_power<3>(p[2])) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0])); return_hessian[2][2] = - (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) * - Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1))); + (2.0 - 12.0 * p[2] + 12.0 * p[2] * p[2]) * + Utilities::fixed_power<2>(p[0] * (1.0 - p[0]) * p[1] * (1.0 - p[1])); } else Assert(false, ExcNotImplemented()); diff --git a/include/deal.II/base/data_out_base.h b/include/deal.II/base/data_out_base.h index e3a7fdf2ba..1d488e0693 100644 --- a/include/deal.II/base/data_out_base.h +++ b/include/deal.II/base/data_out_base.h @@ -1541,9 +1541,9 @@ namespace DataOutBase * return false; */ - return (one(0) < two(0) || - (!(two(0) < one(0)) && - (one(1) < two(1) || (!(two(1) < one(1)) && one(2) < two(2))))); + return (one[0] < two[0] || + (!(two[0] < one[0]) && + (one[1] < two[1] || (!(two[1] < one[1]) && one[2] < two[2])))); } }; diff --git a/include/deal.II/base/point.h b/include/deal.II/base/point.h index 768040ccba..ec31662b7f 100644 --- a/include/deal.II/base/point.h +++ b/include/deal.II/base/point.h @@ -626,7 +626,7 @@ constexpr DEAL_II_HOST_DEVICE typename numbers::NumberTraits::real_type Number sum = internal::NumberType::value(0.0); for (unsigned int i = 0; i < dim; ++i) { - const Number diff = static_cast(this->values[i]) - p(i); + const Number diff = static_cast(this->values[i]) - p[i]; sum += numbers::NumberTraits::abs_square(diff); } diff --git a/include/deal.II/base/polynomial_space.h b/include/deal.II/base/polynomial_space.h index 2739f70662..20fa41c7b1 100644 --- a/include/deal.II/base/polynomial_space.h +++ b/include/deal.II/base/polynomial_space.h @@ -340,7 +340,7 @@ PolynomialSpace::compute_derivative(const unsigned int i, std::vector tmp(order + 1); for (unsigned int d = 0; d < dim; ++d) { - polynomials[indices[d]].value(p(d), tmp); + polynomials[indices[d]].value(p[d], tmp); for (unsigned int j = 0; j < order + 1; ++j) v[d][j] = tmp[j]; } diff --git a/include/deal.II/base/polynomials_rannacher_turek.h b/include/deal.II/base/polynomials_rannacher_turek.h index 14934a1799..e6e47cf750 100644 --- a/include/deal.II/base/polynomials_rannacher_turek.h +++ b/include/deal.II/base/polynomials_rannacher_turek.h @@ -163,23 +163,23 @@ namespace internal *reinterpret_cast *>(&derivative); if (i == 0) { - grad[0] = -2.5 + 3 * p(0); - grad[1] = 1.5 - 3 * p(1); + grad[0] = -2.5 + 3 * p[0]; + grad[1] = 1.5 - 3 * p[1]; } else if (i == 1) { - grad[0] = -0.5 + 3.0 * p(0); - grad[1] = 1.5 - 3.0 * p(1); + grad[0] = -0.5 + 3.0 * p[0]; + grad[1] = 1.5 - 3.0 * p[1]; } else if (i == 2) { - grad[0] = 1.5 - 3.0 * p(0); - grad[1] = -2.5 + 3.0 * p(1); + grad[0] = 1.5 - 3.0 * p[0]; + grad[1] = -2.5 + 3.0 * p[1]; } else if (i == 3) { - grad[0] = 1.5 - 3.0 * p(0); - grad[1] = -0.5 + 3.0 * p(1); + grad[0] = 1.5 - 3.0 * p[0]; + grad[1] = -0.5 + 3.0 * p[1]; } else { diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h index 179a1aa9d1..268ceedd44 100644 --- a/include/deal.II/base/tensor_product_polynomials.h +++ b/include/deal.II/base/tensor_product_polynomials.h @@ -555,7 +555,7 @@ TensorProductPolynomials::compute_derivative( std::vector tmp(5); for (unsigned int d = 0; d < dim; ++d) { - polynomials[indices[d]].value(p(d), tmp); + polynomials[indices[d]].value(p[d], tmp); v[d][0] = tmp[0]; v[d][1] = tmp[1]; v[d][2] = tmp[2]; @@ -739,7 +739,7 @@ AnisotropicPolynomials::compute_derivative(const unsigned int i, std::vector> v(dim, std::vector(order + 1)); for (unsigned int d = 0; d < dim; ++d) - polynomials[d][indices[d]].value(p(d), v[d]); + polynomials[d][indices[d]].value(p[d], v[d]); Tensor derivative; switch (order) diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h index 3753971c44..17f2728233 100644 --- a/include/deal.II/base/tensor_product_polynomials_bubbles.h +++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h @@ -343,10 +343,10 @@ TensorProductPolynomialsBubbles::compute_derivative( // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p) for (unsigned j = 0; j < dim; ++j) derivative_1[d] *= - (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j))); + (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j])); // and multiply with (2*x_i-1)^{r-1} for (unsigned int i = 0; i < q_degree - 1; ++i) - derivative_1[d] *= 2 * p(comp) - 1; + derivative_1[d] *= 2 * p[comp] - 1; } if (q_degree >= 2) @@ -354,11 +354,11 @@ TensorProductPolynomialsBubbles::compute_derivative( // add \prod_{i=1}^d 4*(x_i(1-x_i))(p) double value = 1.; for (unsigned int j = 0; j < dim; ++j) - value *= 4 * p(j) * (1 - p(j)); + value *= 4 * p[j] * (1 - p[j]); // and multiply with grad(2*x_i-1)^{r-1} double tmp = value * 2 * (q_degree - 1); for (unsigned int i = 0; i < q_degree - 2; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; derivative_1[comp] += tmp; } @@ -373,21 +373,21 @@ TensorProductPolynomialsBubbles::compute_derivative( { for (unsigned int c = 0; c < dim; ++c) { - v[c][0] = 4 * p(c) * (1 - p(c)); - v[c][1] = 4 * (1 - 2 * p(c)); + v[c][0] = 4 * p[c] * (1 - p[c]); + v[c][1] = 4 * (1 - 2 * p[c]); v[c][2] = -8; } double tmp = 1.; for (unsigned int i = 0; i < q_degree - 1; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; v[dim][0] = tmp; if (q_degree >= 2) { double tmp = 2 * (q_degree - 1); for (unsigned int i = 0; i < q_degree - 2; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; v[dim][1] = tmp; } else @@ -397,7 +397,7 @@ TensorProductPolynomialsBubbles::compute_derivative( { double tmp = 4 * (q_degree - 2) * (q_degree - 1); for (unsigned int i = 0; i < q_degree - 3; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; v[dim][2] = tmp; } else diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index b0200ce7df..feb89dedf3 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -1792,7 +1792,7 @@ namespace FETools std::vector> q_points_coarse(q_points_fine.size()); for (unsigned int i = 0; i < q_points_fine.size(); ++i) for (unsigned int j = 0; j < dim; ++j) - q_points_coarse[i](j) = q_points_fine[i](j); + q_points_coarse[i][j] = q_points_fine[i][j]; const Quadrature q_coarse(q_points_coarse, fine.get_JxW_values()); FEValues coarse(mapping, @@ -2210,7 +2210,7 @@ namespace FETools std::vector> q_points_coarse(q_points_fine.size()); for (unsigned int q = 0; q < q_points_fine.size(); ++q) for (unsigned int j = 0; j < dim; ++j) - q_points_coarse[q](j) = q_points_fine[q](j); + q_points_coarse[q][j] = q_points_fine[q][j]; Quadrature q_coarse(q_points_coarse, fine.get_JxW_values()); FEValues coarse( coarse_cell->reference_cell() diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index dd56d422f0..7cc9858b2b 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -68,8 +68,8 @@ namespace internal const Point &p) { Assert(spacedim == 1, ExcInternalError()); - return Point<1>((p[0] - vertices[0](0)) / - (vertices[1](0) - vertices[0](0))); + return Point<1>((p[0] - vertices[0][0]) / + (vertices[1][0] - vertices[0][0])); } @@ -86,18 +86,18 @@ namespace internal // For accuracy reasons, we do all arithmetic in extended precision // (long double). This has a noticeable effect on the hit rate for // borderline cases and thus makes the algorithm more robust. - const long double x = p(0); - const long double y = p(1); - - const long double x0 = vertices[0](0); - const long double x1 = vertices[1](0); - const long double x2 = vertices[2](0); - const long double x3 = vertices[3](0); - - const long double y0 = vertices[0](1); - const long double y1 = vertices[1](1); - const long double y2 = vertices[2](1); - const long double y3 = vertices[3](1); + const long double x = p[0]; + const long double y = p[1]; + + const long double x0 = vertices[0][0]; + const long double x1 = vertices[1][0]; + const long double x2 = vertices[2][0]; + const long double x3 = vertices[3][0]; + + const long double y0 = vertices[0][1]; + const long double y1 = vertices[1][1]; + const long double y2 = vertices[2][1]; + const long double y3 = vertices[3][1]; const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3); const long double b = -(x0 - x1 - x2 + x3) * y + (x - 2 * x1 + x3) * y0 - diff --git a/include/deal.II/numerics/vector_tools_boundary.templates.h b/include/deal.II/numerics/vector_tools_boundary.templates.h index c2cfa6bb07..18ae6ed111 100644 --- a/include/deal.II/numerics/vector_tools_boundary.templates.h +++ b/include/deal.II/numerics/vector_tools_boundary.templates.h @@ -1208,9 +1208,9 @@ namespace VectorTools Point shifted_reference_point_2 = reference_quadrature_points[q_point]; - shifted_reference_point_1(edge_coordinate_direction[face][line]) += + shifted_reference_point_1[edge_coordinate_direction[face][line]] += tol; - shifted_reference_point_2(edge_coordinate_direction[face][line]) -= + shifted_reference_point_2[edge_coordinate_direction[face][line]] -= tol; Tensor<1, dim> tangential = (0.5 * diff --git a/include/deal.II/particles/particle_accessor.h b/include/deal.II/particles/particle_accessor.h index 2c51fc98a7..67f3097842 100644 --- a/include/deal.II/particles/particle_accessor.h +++ b/include/deal.II/particles/particle_accessor.h @@ -571,12 +571,12 @@ namespace Particles Point location; for (unsigned int i = 0; i < spacedim; ++i) - location(i) = *pdata++; + location[i] = *pdata++; set_location(location); Point reference_location; for (unsigned int i = 0; i < dim; ++i) - reference_location(i) = *pdata++; + reference_location[i] = *pdata++; set_reference_location(reference_location); // See if there are properties to load diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index a4e1b0d661..098b08fd29 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -443,7 +443,7 @@ namespace DataOutBase Point<3> int_pt; for (unsigned int d = 0; d < dim; ++d) - int_pt(d) = p(d); + int_pt[d] = p[d]; const Map3DPoint::const_iterator it = existing_points.find(int_pt); unsigned int internal_ind; @@ -487,7 +487,7 @@ namespace DataOutBase { for (unsigned int d = 0; d < node_dim; ++d) node_data[node_dim * existing_point.second + d] = - existing_point.first(d); + existing_point.first[d]; } } @@ -1316,13 +1316,13 @@ namespace { float data[dim]; for (unsigned int d = 0; d < dim; ++d) - data[d] = p(d); + data[d] = p[d]; stream.write(reinterpret_cast(data), dim * sizeof(*data)); } else { for (unsigned int d = 0; d < dim; ++d) - stream << p(d) << '\t'; + stream << p[d] << '\t'; stream << '\n'; } } @@ -1459,7 +1459,7 @@ namespace { Assert(selected_component != numbers::invalid_unsigned_int, ExcNotInitialized()); - stream << p(selected_component) << ' '; + stream << p[selected_component] << ' '; } @@ -1533,7 +1533,7 @@ namespace { Assert(selected_component != numbers::invalid_unsigned_int, ExcNotInitialized()); - stream << p(selected_component) << '\n'; + stream << p[selected_component] << '\n'; } @@ -1605,7 +1605,7 @@ namespace stream << index + 1 << " "; // write out coordinates for (unsigned int i = 0; i < dim; ++i) - stream << p(i) << ' '; + stream << p[i] << ' '; // fill with zeroes for (unsigned int i = dim; i < 3; ++i) stream << "0 "; @@ -2197,8 +2197,8 @@ namespace DataOutBase RgbValues rgb_values = {0, 0, 0}; // A difficult color scale: - // xmin = black (1) - // 3/4*xmin+1/4*xmax = blue (2) + // xmin = black [1] + // 3/4*xmin+1/4*xmax = blue [2] // 1/2*xmin+1/2*xmax = green (3) // 1/4*xmin+3/4*xmax = red (4) // xmax = white (5) @@ -2209,17 +2209,17 @@ namespace DataOutBase // / /\ / /\ / // ____/ __/ \/ / \__/ - // { 0 (1) - (3) + // { 0 [1] - (3) // r = { ( 4*x-2*xmin+2*xmax)/(xmax-xmin) (3) - (4) // { 1 (4) - (5) // - // { 0 (1) - (2) - // g = { ( 4*x-3*xmin- xmax)/(xmax-xmin) (2) - (3) + // { 0 [1] - [2] + // g = { ( 4*x-3*xmin- xmax)/(xmax-xmin) [2] - (3) // { (-4*x+ xmin+3*xmax)/(xmax-xmin) (3) - (4) // { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5) // - // { ( 4*x-4*xmin )/(xmax-xmin) (1) - (2) - // b = { (-4*x+2*xmin+2*xmax)/(xmax-xmin) (2) - (3) + // { ( 4*x-4*xmin )/(xmax-xmin) [1] - [2] + // b = { (-4*x+2*xmin+2*xmax)/(xmax-xmin) [2] - (3) // { 0 (3) - (4) // { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5) @@ -4144,37 +4144,37 @@ namespace DataOutBase const unsigned int jr = (j == n_subdivisions) ? j : (j + 1); - h1(0) = - ver[ir * d1 + j * d2](0) - ver[il * d1 + j * d2](0); - h1(1) = patch.data(0, ir * d1 + j * d2) - + h1[0] = + ver[ir * d1 + j * d2][0] - ver[il * d1 + j * d2][0]; + h1[1] = patch.data(0, ir * d1 + j * d2) - patch.data(0, il * d1 + j * d2); - h1(2) = - ver[ir * d1 + j * d2](1) - ver[il * d1 + j * d2](1); + h1[2] = + ver[ir * d1 + j * d2][1] - ver[il * d1 + j * d2][1]; - h2(0) = - ver[i * d1 + jr * d2](0) - ver[i * d1 + jl * d2](0); - h2(1) = patch.data(0, i * d1 + jr * d2) - + h2[0] = + ver[i * d1 + jr * d2][0] - ver[i * d1 + jl * d2][0]; + h2[1] = patch.data(0, i * d1 + jr * d2) - patch.data(0, i * d1 + jl * d2); - h2(2) = - ver[i * d1 + jr * d2](1) - ver[i * d1 + jl * d2](1); + h2[2] = + ver[i * d1 + jr * d2][1] - ver[i * d1 + jl * d2][1]; - nrml[i * d1 + j * d2](0) = - h1(1) * h2(2) - h1(2) * h2(1); - nrml[i * d1 + j * d2](1) = - h1(2) * h2(0) - h1(0) * h2(2); - nrml[i * d1 + j * d2](2) = - h1(0) * h2(1) - h1(1) * h2(0); + nrml[i * d1 + j * d2][0] = + h1[1] * h2[2] - h1[2] * h2[1]; + nrml[i * d1 + j * d2][1] = + h1[2] * h2[0] - h1[0] * h2[2]; + nrml[i * d1 + j * d2][2] = + h1[0] * h2[1] - h1[1] * h2[0]; // normalize Vector - double norm = std::hypot(nrml[i * d1 + j * d2](0), - nrml[i * d1 + j * d2](1), - nrml[i * d1 + j * d2](2)); + double norm = std::hypot(nrml[i * d1 + j * d2][0], + nrml[i * d1 + j * d2][1], + nrml[i * d1 + j * d2][2]); - if (nrml[i * d1 + j * d2](1) < 0) + if (nrml[i * d1 + j * d2][1] < 0) norm *= -1.; for (unsigned int k = 0; k < 3; ++k) - nrml[i * d1 + j * d2](k) /= norm; + nrml[i * d1 + j * d2][k] /= norm; } } @@ -4190,62 +4190,62 @@ namespace DataOutBase // down/right triangle out << "smooth_triangle {" << '\n' - << "\t<" << ver[dl](0) << "," << patch.data(0, dl) - << "," << ver[dl](1) << ">, <" << nrml[dl](0) - << ", " << nrml[dl](1) << ", " << nrml[dl](2) + << "\t<" << ver[dl][0] << "," << patch.data(0, dl) + << "," << ver[dl][1] << ">, <" << nrml[dl][0] + << ", " << nrml[dl][1] << ", " << nrml[dl][2] << ">," << '\n'; - out << " \t<" << ver[dl + d1](0) << "," - << patch.data(0, dl + d1) << "," << ver[dl + d1](1) - << ">, <" << nrml[dl + d1](0) << ", " - << nrml[dl + d1](1) << ", " << nrml[dl + d1](2) + out << " \t<" << ver[dl + d1][0] << "," + << patch.data(0, dl + d1) << "," << ver[dl + d1][1] + << ">, <" << nrml[dl + d1][0] << ", " + << nrml[dl + d1][1] << ", " << nrml[dl + d1][2] << ">," << '\n'; - out << "\t<" << ver[dl + d1 + d2](0) << "," + out << "\t<" << ver[dl + d1 + d2][0] << "," << patch.data(0, dl + d1 + d2) << "," - << ver[dl + d1 + d2](1) << ">, <" - << nrml[dl + d1 + d2](0) << ", " - << nrml[dl + d1 + d2](1) << ", " - << nrml[dl + d1 + d2](2) << ">}" << '\n'; + << ver[dl + d1 + d2][1] << ">, <" + << nrml[dl + d1 + d2][0] << ", " + << nrml[dl + d1 + d2][1] << ", " + << nrml[dl + d1 + d2][2] << ">}" << '\n'; // upper/left triangle out << "smooth_triangle {" << '\n' - << "\t<" << ver[dl](0) << "," << patch.data(0, dl) - << "," << ver[dl](1) << ">, <" << nrml[dl](0) - << ", " << nrml[dl](1) << ", " << nrml[dl](2) + << "\t<" << ver[dl][0] << "," << patch.data(0, dl) + << "," << ver[dl][1] << ">, <" << nrml[dl][0] + << ", " << nrml[dl][1] << ", " << nrml[dl][2] << ">," << '\n'; - out << "\t<" << ver[dl + d1 + d2](0) << "," + out << "\t<" << ver[dl + d1 + d2][0] << "," << patch.data(0, dl + d1 + d2) << "," - << ver[dl + d1 + d2](1) << ">, <" - << nrml[dl + d1 + d2](0) << ", " - << nrml[dl + d1 + d2](1) << ", " - << nrml[dl + d1 + d2](2) << ">," << '\n'; - out << "\t<" << ver[dl + d2](0) << "," - << patch.data(0, dl + d2) << "," << ver[dl + d2](1) - << ">, <" << nrml[dl + d2](0) << ", " - << nrml[dl + d2](1) << ", " << nrml[dl + d2](2) + << ver[dl + d1 + d2][1] << ">, <" + << nrml[dl + d1 + d2][0] << ", " + << nrml[dl + d1 + d2][1] << ", " + << nrml[dl + d1 + d2][2] << ">," << '\n'; + out << "\t<" << ver[dl + d2][0] << "," + << patch.data(0, dl + d2) << "," << ver[dl + d2][1] + << ">, <" << nrml[dl + d2][0] << ", " + << nrml[dl + d2][1] << ", " << nrml[dl + d2][2] << ">}" << '\n'; } else { // writing standard triangles down/right triangle out << "triangle {" << '\n' - << "\t<" << ver[dl](0) << "," << patch.data(0, dl) - << "," << ver[dl](1) << ">," << '\n'; - out << "\t<" << ver[dl + d1](0) << "," - << patch.data(0, dl + d1) << "," << ver[dl + d1](1) + << "\t<" << ver[dl][0] << "," << patch.data(0, dl) + << "," << ver[dl][1] << ">," << '\n'; + out << "\t<" << ver[dl + d1][0] << "," + << patch.data(0, dl + d1) << "," << ver[dl + d1][1] << ">," << '\n'; - out << "\t<" << ver[dl + d1 + d2](0) << "," + out << "\t<" << ver[dl + d1 + d2][0] << "," << patch.data(0, dl + d1 + d2) << "," - << ver[dl + d1 + d2](1) << ">}" << '\n'; + << ver[dl + d1 + d2][1] << ">}" << '\n'; // upper/left triangle out << "triangle {" << '\n' - << "\t<" << ver[dl](0) << "," << patch.data(0, dl) - << "," << ver[dl](1) << ">," << '\n'; - out << "\t<" << ver[dl + d1 + d2](0) << "," + << "\t<" << ver[dl][0] << "," << patch.data(0, dl) + << "," << ver[dl][1] << ">," << '\n'; + out << "\t<" << ver[dl + d1 + d2][0] << "," << patch.data(0, dl + d1 + d2) << "," - << ver[dl + d1 + d2](1) << ">," << '\n'; - out << "\t<" << ver[dl + d2](0) << "," - << patch.data(0, dl + d2) << "," << ver[dl + d2](1) + << ver[dl + d1 + d2][1] << ">," << '\n'; + out << "\t<" << ver[dl + d2][0] << "," + << patch.data(0, dl + d2) << "," << ver[dl + d2][1] << ">}" << '\n'; } } @@ -4263,8 +4263,8 @@ namespace DataOutBase << " v_steps 0" << '\n'; for (int i = 0; i < 16; ++i) { - out << "\t<" << ver[i](0) << "," << patch.data(0, i) << "," - << ver[i](1) << ">"; + out << "\t<" << ver[i][0] << "," << patch.data(0, i) << "," + << ver[i][1] << ">"; if (i != 15) out << ","; out << '\n'; @@ -4425,7 +4425,7 @@ namespace DataOutBase case 3: // Copy z-coordinates into the height vector for (unsigned int i = 0; i < 4; ++i) - heights[i] = points[i](2); + heights[i] = points[i][2]; break; default: Assert(false, ExcNotImplemented()); @@ -4452,11 +4452,11 @@ namespace DataOutBase sz = std::sin(flags.turn_angle * 2 * pi / 360.); for (unsigned int vertex = 0; vertex < 4; ++vertex) { - const double x = points[vertex](0), y = points[vertex](1), + const double x = points[vertex][0], y = points[vertex][1], z = -heights[vertex]; - eps_cell.vertices[vertex](0) = -cz * x + sz * y; - eps_cell.vertices[vertex](1) = + eps_cell.vertices[vertex][0] = -cz * x + sz * y; + eps_cell.vertices[vertex][1] = -cx * sz * x - cx * cz * y - sx * z; // ( 1 0 0 ) @@ -4486,8 +4486,8 @@ namespace DataOutBase -(heights[0] + heights[1] + heights[2] + heights[3]) / 4; // compute the depth into the picture - eps_cell.depth = -sx * sz * center_point(0) - - sx * cz * center_point(1) + cx * center_height; + eps_cell.depth = -sx * sz * center_point[0] - + sx * cz * center_point[1] + cx * center_height; if (flags.draw_cells && flags.shade_cells) { @@ -4533,18 +4533,18 @@ namespace DataOutBase // find out minimum and maximum x and y coordinates to compute offsets and // scaling factors - double x_min = cells.begin()->vertices[0](0); + double x_min = cells.begin()->vertices[0][0]; double x_max = x_min; - double y_min = cells.begin()->vertices[0](1); + double y_min = cells.begin()->vertices[0][1]; double y_max = y_min; for (const auto &cell : cells) for (const auto &vertex : cell.vertices) { - x_min = std::min(x_min, vertex(0)); - x_max = std::max(x_max, vertex(0)); - y_min = std::min(y_min, vertex(1)); - y_max = std::max(y_max, vertex(1)); + x_min = std::min(x_min, vertex[0]); + x_max = std::max(x_max, vertex[0]); + y_min = std::min(y_min, vertex[1]); + y_max = std::max(y_max, vertex[1]); } // scale in x-direction such that in the output 0 <= x <= 300. don't scale diff --git a/source/base/flow_function.cc b/source/base/flow_function.cc index 55c13d91d6..0814474acb 100644 --- a/source/base/flow_function.cc +++ b/source/base/flow_function.cc @@ -217,7 +217,7 @@ namespace Functions // the radius. double r2 = 0; for (unsigned int d = 1; d < dim; ++d) - r2 += p(d) * p(d); + r2 += p[d] * p[d]; r2 *= inv_sqr_radius; // x-velocity @@ -226,7 +226,7 @@ namespace Functions for (unsigned int d = 1; d < dim; ++d) values[d][k] = 0.; // pressure - values[dim][k] = -2 * (dim - 1) * inv_sqr_radius * p(0) / Reynolds + + values[dim][k] = -2 * (dim - 1) * inv_sqr_radius * p[0] / Reynolds + this->mean_pressure; } } @@ -252,7 +252,7 @@ namespace Functions // x-velocity values[0][k][0] = 0.; for (unsigned int d = 1; d < dim; ++d) - values[0][k][d] = -2. * p(d) * inv_sqr_radius; + values[0][k][d] = -2. * p[d] * inv_sqr_radius; // other velocities for (unsigned int d = 1; d < dim; ++d) values[d][k] = 0.; @@ -317,8 +317,8 @@ namespace Functions for (unsigned int k = 0; k < n; ++k) { const Point &p = points[k]; - const double x = numbers::PI / 2. * p(0); - const double y = numbers::PI / 2. * p(1); + const double x = numbers::PI / 2. * p[0]; + const double y = numbers::PI / 2. * p[1]; const double cx = std::cos(x); const double cy = std::cos(y); const double sx = std::sin(x); @@ -332,7 +332,7 @@ namespace Functions } else if (dim == 3) { - const double z = numbers::PI / 2. * p(2); + const double z = numbers::PI / 2. * p[2]; const double cz = std::cos(z); const double sz = std::sin(z); @@ -366,8 +366,8 @@ namespace Functions for (unsigned int k = 0; k < n; ++k) { const Point &p = points[k]; - const double x = numbers::PI / 2. * p(0); - const double y = numbers::PI / 2. * p(1); + const double x = numbers::PI / 2. * p[0]; + const double y = numbers::PI / 2. * p[1]; const double c2x = std::cos(2 * x); const double c2y = std::cos(2 * y); const double s2x = std::sin(2 * x); @@ -386,7 +386,7 @@ namespace Functions } else if (dim == 3) { - const double z = numbers::PI / 2. * p(2); + const double z = numbers::PI / 2. * p[2]; const double c2z = std::cos(2 * z); const double s2z = std::sin(2 * z); const double cz2 = .5 + .5 * c2z; // cos^2 z @@ -446,8 +446,8 @@ namespace Functions for (unsigned int k = 0; k < n; ++k) { const Point &p = points[k]; - const double x = numbers::PI / 2. * p(0); - const double y = numbers::PI / 2. * p(1); + const double x = numbers::PI / 2. * p[0]; + const double y = numbers::PI / 2. * p[1]; const double c2x = std::cos(2 * x); const double c2y = std::cos(2 * y); const double s2x = std::sin(2 * x); @@ -464,7 +464,7 @@ namespace Functions } else if (dim == 3) { - const double z = numbers::PI * p(2); + const double z = numbers::PI * p[2]; const double c2z = std::cos(2 * z); const double s2z = std::sin(2 * z); @@ -556,8 +556,8 @@ namespace Functions for (unsigned int k = 0; k < n; ++k) { const Point<2> &p = points[k]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; if ((x < 0) || (y < 0)) { @@ -596,8 +596,8 @@ namespace Functions for (unsigned int k = 0; k < n; ++k) { const Point<2> &p = points[k]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; if ((x < 0) || (y < 0)) { @@ -689,8 +689,8 @@ namespace Functions for (unsigned int k = 0; k < n; ++k) { const Point<2> &p = points[k]; - const double x = p(0); - const double y = 2. * numbers::PI * p(1); + const double x = p[0]; + const double y = 2. * numbers::PI * p[1]; const double elx = std::exp(lbda * x); values[0][k] = 1. - elx * std::cos(y); @@ -713,8 +713,8 @@ namespace Functions for (unsigned int i = 0; i < n; ++i) { - const double x = points[i](0); - const double y = points[i](1); + const double x = points[i][0]; + const double y = points[i][1]; const double elx = std::exp(lbda * x); const double cy = std::cos(2 * numbers::PI * y); @@ -748,8 +748,8 @@ namespace Functions for (unsigned int k = 0; k < n; ++k) { const Point<2> &p = points[k]; - const double x = p(0); - const double y = zp * p(1); + const double x = p[0]; + const double y = zp * p[1]; const double elx = std::exp(lbda * x); const double u = 1. - elx * std::cos(y); const double ux = -lbda * elx * std::cos(y); diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc index 7e11ff326a..31e5c4d092 100644 --- a/source/base/function_lib.cc +++ b/source/base/function_lib.cc @@ -129,7 +129,7 @@ namespace Functions Q1WedgeFunction::value(const Point &p, const unsigned int) const { Assert(dim >= 2, ExcInternalError()); - return p(0) * p(1); + return p[0] * p[1]; } @@ -147,7 +147,7 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { const Point &p = points[i]; - values[i] = p(0) * p(1); + values[i] = p[0] * p[1]; } } @@ -166,7 +166,7 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { const Point &p = points[i]; - values[i](0) = p(0) * p(1); + values[i](0) = p[0] * p[1]; } } @@ -202,8 +202,8 @@ namespace Functions { Assert(dim >= 2, ExcInternalError()); Tensor<1, dim> erg; - erg[0] = p(1); - erg[1] = p(0); + erg[0] = p[1]; + erg[1] = p[0]; return erg; } @@ -221,8 +221,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { - gradients[i][0] = points[i](1); - gradients[i][1] = points[i](0); + gradients[i][0] = points[i][1]; + gradients[i][1] = points[i][0]; } } @@ -241,8 +241,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { - gradients[i][0][0] = points[i](1); - gradients[i][0][1] = points[i](0); + gradients[i][0][0] = points[i][1]; + gradients[i][0][1] = points[i][0]; } } @@ -263,11 +263,11 @@ namespace Functions switch (dim) { case 1: - return 1. - p(0) * p(0) + offset; + return 1. - p[0] * p[0] + offset; case 2: - return (1. - p(0) * p(0)) * (1. - p(1) * p(1)) + offset; + return (1. - p[0] * p[0]) * (1. - p[1] * p[1]) + offset; case 3: - return (1. - p(0) * p(0)) * (1. - p(1) * p(1)) * (1. - p(2) * p(2)) + + return (1. - p[0] * p[0]) * (1. - p[1] * p[1]) * (1. - p[2] * p[2]) + offset; default: Assert(false, ExcNotImplemented()); @@ -290,14 +290,14 @@ namespace Functions switch (dim) { case 1: - values[i] = 1. - p(0) * p(0) + offset; + values[i] = 1. - p[0] * p[0] + offset; break; case 2: - values[i] = (1. - p(0) * p(0)) * (1. - p(1) * p(1)) + offset; + values[i] = (1. - p[0] * p[0]) * (1. - p[1] * p[1]) + offset; break; case 3: values[i] = - (1. - p(0) * p(0)) * (1. - p(1) * p(1)) * (1. - p(2) * p(2)) + + (1. - p[0] * p[0]) * (1. - p[1] * p[1]) * (1. - p[2] * p[2]) + offset; break; default: @@ -317,11 +317,11 @@ namespace Functions case 1: return -2.; case 2: - return -2. * ((1. - p(0) * p(0)) + (1. - p(1) * p(1))); + return -2. * ((1. - p[0] * p[0]) + (1. - p[1] * p[1])); case 3: - return -2. * ((1. - p(0) * p(0)) * (1. - p(1) * p(1)) + - (1. - p(1) * p(1)) * (1. - p(2) * p(2)) + - (1. - p(2) * p(2)) * (1. - p(0) * p(0))); + return -2. * ((1. - p[0] * p[0]) * (1. - p[1] * p[1]) + + (1. - p[1] * p[1]) * (1. - p[2] * p[2]) + + (1. - p[2] * p[2]) * (1. - p[0] * p[0])); default: Assert(false, ExcNotImplemented()); } @@ -346,12 +346,12 @@ namespace Functions values[i] = -2.; break; case 2: - values[i] = -2. * ((1. - p(0) * p(0)) + (1. - p(1) * p(1))); + values[i] = -2. * ((1. - p[0] * p[0]) + (1. - p[1] * p[1])); break; case 3: - values[i] = -2. * ((1. - p(0) * p(0)) * (1. - p(1) * p(1)) + - (1. - p(1) * p(1)) * (1. - p(2) * p(2)) + - (1. - p(2) * p(2)) * (1. - p(0) * p(0))); + values[i] = -2. * ((1. - p[0] * p[0]) * (1. - p[1] * p[1]) + + (1. - p[1] * p[1]) * (1. - p[2] * p[2]) + + (1. - p[2] * p[2]) * (1. - p[0] * p[0])); break; default: Assert(false, ExcNotImplemented()); @@ -367,16 +367,16 @@ namespace Functions switch (dim) { case 1: - result[0] = -2. * p(0); + result[0] = -2. * p[0]; break; case 2: - result[0] = -2. * p(0) * (1. - p(1) * p(1)); - result[1] = -2. * p(1) * (1. - p(0) * p(0)); + result[0] = -2. * p[0] * (1. - p[1] * p[1]); + result[1] = -2. * p[1] * (1. - p[0] * p[0]); break; case 3: - result[0] = -2. * p(0) * (1. - p(1) * p(1)) * (1. - p(2) * p(2)); - result[1] = -2. * p(1) * (1. - p(0) * p(0)) * (1. - p(2) * p(2)); - result[2] = -2. * p(2) * (1. - p(0) * p(0)) * (1. - p(1) * p(1)); + result[0] = -2. * p[0] * (1. - p[1] * p[1]) * (1. - p[2] * p[2]); + result[1] = -2. * p[1] * (1. - p[0] * p[0]) * (1. - p[2] * p[2]); + result[2] = -2. * p[2] * (1. - p[0] * p[0]) * (1. - p[1] * p[1]); break; default: Assert(false, ExcNotImplemented()); @@ -399,19 +399,19 @@ namespace Functions switch (dim) { case 1: - gradients[i][0] = -2. * p(0); + gradients[i][0] = -2. * p[0]; break; case 2: - gradients[i][0] = -2. * p(0) * (1. - p(1) * p(1)); - gradients[i][1] = -2. * p(1) * (1. - p(0) * p(0)); + gradients[i][0] = -2. * p[0] * (1. - p[1] * p[1]); + gradients[i][1] = -2. * p[1] * (1. - p[0] * p[0]); break; case 3: gradients[i][0] = - -2. * p(0) * (1. - p(1) * p(1)) * (1. - p(2) * p(2)); + -2. * p[0] * (1. - p[1] * p[1]) * (1. - p[2] * p[2]); gradients[i][1] = - -2. * p(1) * (1. - p(0) * p(0)) * (1. - p(2) * p(2)); + -2. * p[1] * (1. - p[0] * p[0]) * (1. - p[2] * p[2]); gradients[i][2] = - -2. * p(2) * (1. - p(0) * p(0)) * (1. - p(1) * p(1)); + -2. * p[2] * (1. - p[0] * p[0]) * (1. - p[1] * p[1]); break; default: Assert(false, ExcNotImplemented()); @@ -435,14 +435,14 @@ namespace Functions switch (dim) { case 1: - return std::cos(numbers::PI_2 * p(0)); + return std::cos(numbers::PI_2 * p[0]); case 2: - return std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); + return std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); case 3: - return std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); + return std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); default: Assert(false, ExcNotImplemented()); } @@ -489,16 +489,16 @@ namespace Functions { case 1: return -numbers::PI_2 * numbers::PI_2 * - std::cos(numbers::PI_2 * p(0)); + std::cos(numbers::PI_2 * p[0]); case 2: return -2 * numbers::PI_2 * numbers::PI_2 * - std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); + std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); case 3: return -3 * numbers::PI_2 * numbers::PI_2 * - std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); + std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); default: Assert(false, ExcNotImplemented()); } @@ -526,24 +526,24 @@ namespace Functions switch (dim) { case 1: - result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)); + result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]); break; case 2: - result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); - result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)); + result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); + result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]); break; case 3: - result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - result[2] = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); + result[0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + result[1] = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + result[2] = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); break; default: Assert(false, ExcNotImplemented()); @@ -566,26 +566,26 @@ namespace Functions switch (dim) { case 1: - gradients[i][0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)); + gradients[i][0] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]); break; case 2: gradients[i][0] = -numbers::PI_2 * - std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); + std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); gradients[i][1] = -numbers::PI_2 * - std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)); + std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]); break; case 3: gradients[i][0] = - -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * std::cos(numbers::PI_2 * p(2)); + -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * std::cos(numbers::PI_2 * p[2]); gradients[i][1] = - -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * std::cos(numbers::PI_2 * p(2)); + -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * std::cos(numbers::PI_2 * p[2]); gradients[i][2] = - -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * std::sin(numbers::PI_2 * p(2)); + -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * std::sin(numbers::PI_2 * p[2]); break; default: Assert(false, ExcNotImplemented()); @@ -603,14 +603,14 @@ namespace Functions switch (dim) { case 1: - result[0][0] = -pi2 * std::cos(numbers::PI_2 * p(0)); + result[0][0] = -pi2 * std::cos(numbers::PI_2 * p[0]); break; case 2: { - const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); - const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)); + const double coco = -pi2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); + const double sisi = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]); result[0][0] = coco; result[1][1] = coco; // for SymmetricTensor we assign [ij] and [ji] simultaneously: @@ -619,18 +619,18 @@ namespace Functions break; case 3: { - const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); - const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); + const double cococo = -pi2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + const double sisico = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + const double sicosi = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); + const double cosisi = pi2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); result[0][0] = cococo; result[1][1] = cococo; @@ -665,14 +665,14 @@ namespace Functions switch (dim) { case 1: - hessians[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p(0)); + hessians[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p[0]); break; case 2: { - const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); - const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)); + const double coco = -pi2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); + const double sisi = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]); hessians[i][0][0] = coco; hessians[i][1][1] = coco; // for SymmetricTensor we assign [ij] and [ji] simultaneously: @@ -681,18 +681,18 @@ namespace Functions break; case 3: { - const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); - const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); + const double cococo = -pi2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + const double sisico = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + const double sicosi = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); + const double cosisi = pi2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); hessians[i][0][0] = cococo; hessians[i][1][1] = cococo; @@ -728,14 +728,14 @@ namespace Functions switch (dim) { case 1: - return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(0))); + return (-numbers::PI_2 * std::sin(numbers::PI_2 * p[0])); case 2: - return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1))); + return (-numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1])); case 3: - return (-numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)) * - std::cos(numbers::PI_2 * p(d2))); + return (-numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]) * + std::cos(numbers::PI_2 * p[d2])); default: Assert(false, ExcNotImplemented()); } @@ -752,24 +752,24 @@ namespace Functions switch (dim) { case 1: - result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)); + result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]); break; case 2: - result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); - result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)); + result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); + result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]); break; case 3: - result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - result(2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); + result(0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + result(1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + result(2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); break; default: Assert(false, ExcNotImplemented()); @@ -795,16 +795,16 @@ namespace Functions switch (dim) { case 1: - values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d)); + values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[d]); break; case 2: - values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)); + values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]); break; case 3: - values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)) * - std::cos(numbers::PI_2 * p(d2)); + values[i] = -numbers::PI_2 * std::sin(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]) * + std::cos(numbers::PI_2 * p[d2]); break; default: Assert(false, ExcNotImplemented()); @@ -828,24 +828,24 @@ namespace Functions switch (dim) { case 1: - values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)); + values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]); break; case 2: - values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); - values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)); + values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); + values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]); break; case 3: - values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - values[i](2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); + values[i](0) = -numbers::PI_2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + values[i](1) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + values[i](2) = -numbers::PI_2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); break; default: Assert(false, ExcNotImplemented()); @@ -877,24 +877,24 @@ namespace Functions switch (dim) { case 1: - result[0] = -pi2 * std::cos(numbers::PI_2 * p(0)); + result[0] = -pi2 * std::cos(numbers::PI_2 * p[0]); break; case 2: - result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)); - result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) * - std::sin(numbers::PI_2 * p(d1)); + result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]); + result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) * + std::sin(numbers::PI_2 * p[d1]); break; case 3: - result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)) * - std::cos(numbers::PI_2 * p(d2)); - result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) * - std::sin(numbers::PI_2 * p(d1)) * - std::cos(numbers::PI_2 * p(d2)); - result[d2] = pi2 * std::sin(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)) * - std::sin(numbers::PI_2 * p(d2)); + result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]) * + std::cos(numbers::PI_2 * p[d2]); + result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) * + std::sin(numbers::PI_2 * p[d1]) * + std::cos(numbers::PI_2 * p[d2]); + result[d2] = pi2 * std::sin(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]) * + std::sin(numbers::PI_2 * p[d2]); break; default: Assert(false, ExcNotImplemented()); @@ -924,24 +924,24 @@ namespace Functions switch (dim) { case 1: - result[0] = -pi2 * std::cos(numbers::PI_2 * p(0)); + result[0] = -pi2 * std::cos(numbers::PI_2 * p[0]); break; case 2: - result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)); - result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) * - std::sin(numbers::PI_2 * p(d1)); + result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]); + result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) * + std::sin(numbers::PI_2 * p[d1]); break; case 3: - result[d] = -pi2 * std::cos(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)) * - std::cos(numbers::PI_2 * p(d2)); - result[d1] = pi2 * std::sin(numbers::PI_2 * p(d)) * - std::sin(numbers::PI_2 * p(d1)) * - std::cos(numbers::PI_2 * p(d2)); - result[d2] = pi2 * std::sin(numbers::PI_2 * p(d)) * - std::cos(numbers::PI_2 * p(d1)) * - std::sin(numbers::PI_2 * p(d2)); + result[d] = -pi2 * std::cos(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]) * + std::cos(numbers::PI_2 * p[d2]); + result[d1] = pi2 * std::sin(numbers::PI_2 * p[d]) * + std::sin(numbers::PI_2 * p[d1]) * + std::cos(numbers::PI_2 * p[d2]); + result[d2] = pi2 * std::sin(numbers::PI_2 * p[d]) * + std::cos(numbers::PI_2 * p[d1]) * + std::sin(numbers::PI_2 * p[d2]); break; default: Assert(false, ExcNotImplemented()); @@ -965,14 +965,14 @@ namespace Functions switch (dim) { case 1: - gradients[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p(0)); + gradients[i][0][0] = -pi2 * std::cos(numbers::PI_2 * p[0]); break; case 2: { - const double coco = -pi2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)); - const double sisi = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)); + const double coco = -pi2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]); + const double sisi = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]); gradients[i][0][0] = coco; gradients[i][1][1] = coco; gradients[i][0][1] = sisi; @@ -981,18 +981,18 @@ namespace Functions break; case 3: { - const double cococo = -pi2 * std::cos(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - const double sisico = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::cos(numbers::PI_2 * p(2)); - const double sicosi = pi2 * std::sin(numbers::PI_2 * p(0)) * - std::cos(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); - const double cosisi = pi2 * std::cos(numbers::PI_2 * p(0)) * - std::sin(numbers::PI_2 * p(1)) * - std::sin(numbers::PI_2 * p(2)); + const double cococo = -pi2 * std::cos(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + const double sisico = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::cos(numbers::PI_2 * p[2]); + const double sicosi = pi2 * std::sin(numbers::PI_2 * p[0]) * + std::cos(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); + const double cosisi = pi2 * std::cos(numbers::PI_2 * p[0]) * + std::sin(numbers::PI_2 * p[1]) * + std::sin(numbers::PI_2 * p[2]); gradients[i][0][0] = cococo; gradients[i][1][1] = cococo; @@ -1021,11 +1021,11 @@ namespace Functions switch (dim) { case 1: - return std::exp(p(0)); + return std::exp(p[0]); case 2: - return std::exp(p(0)) * std::exp(p(1)); + return std::exp(p[0]) * std::exp(p[1]); case 3: - return std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2)); + return std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]); default: Assert(false, ExcNotImplemented()); } @@ -1047,13 +1047,13 @@ namespace Functions switch (dim) { case 1: - values[i] = std::exp(p(0)); + values[i] = std::exp(p[0]); break; case 2: - values[i] = std::exp(p(0)) * std::exp(p(1)); + values[i] = std::exp(p[0]) * std::exp(p[1]); break; case 3: - values[i] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2)); + values[i] = std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]); break; default: Assert(false, ExcNotImplemented()); @@ -1068,11 +1068,11 @@ namespace Functions switch (dim) { case 1: - return std::exp(p(0)); + return std::exp(p[0]); case 2: - return 2 * std::exp(p(0)) * std::exp(p(1)); + return 2 * std::exp(p[0]) * std::exp(p[1]); case 3: - return 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2)); + return 3 * std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]); default: Assert(false, ExcNotImplemented()); } @@ -1094,13 +1094,13 @@ namespace Functions switch (dim) { case 1: - values[i] = std::exp(p(0)); + values[i] = std::exp(p[0]); break; case 2: - values[i] = 2 * std::exp(p(0)) * std::exp(p(1)); + values[i] = 2 * std::exp(p[0]) * std::exp(p[1]); break; case 3: - values[i] = 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2)); + values[i] = 3 * std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]); break; default: Assert(false, ExcNotImplemented()); @@ -1116,14 +1116,14 @@ namespace Functions switch (dim) { case 1: - result[0] = std::exp(p(0)); + result[0] = std::exp(p[0]); break; case 2: - result[0] = std::exp(p(0)) * std::exp(p(1)); + result[0] = std::exp(p[0]) * std::exp(p[1]); result[1] = result[0]; break; case 3: - result[0] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2)); + result[0] = std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]); result[1] = result[0]; result[2] = result[0]; break; @@ -1148,15 +1148,15 @@ namespace Functions switch (dim) { case 1: - gradients[i][0] = std::exp(p(0)); + gradients[i][0] = std::exp(p[0]); break; case 2: - gradients[i][0] = std::exp(p(0)) * std::exp(p(1)); + gradients[i][0] = std::exp(p[0]) * std::exp(p[1]); gradients[i][1] = gradients[i][0]; break; case 3: gradients[i][0] = - std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2)); + std::exp(p[0]) * std::exp(p[1]) * std::exp(p[2]); gradients[i][1] = gradients[i][0]; gradients[i][2] = gradients[i][0]; break; @@ -1172,8 +1172,8 @@ namespace Functions double LSingularityFunction::value(const Point<2> &p, const unsigned int) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; if ((x >= 0) && (y >= 0)) return 0.; @@ -1196,8 +1196,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { - const double x = points[i](0); - const double y = points[i](1); + const double x = points[i][0]; + const double y = points[i][1]; if ((x >= 0) && (y >= 0)) values[i] = 0.; @@ -1225,8 +1225,8 @@ namespace Functions { Assert(values[i].size() == 1, ExcDimensionMismatch(values[i].size(), 1)); - const double x = points[i](0); - const double y = points[i](1); + const double x = points[i][0]; + const double y = points[i][1]; if ((x >= 0) && (y >= 0)) values[i](0) = 0.; @@ -1268,8 +1268,8 @@ namespace Functions Tensor<1, 2> LSingularityFunction::gradient(const Point<2> &p, const unsigned int) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(y, -x) + numbers::PI; const double r43 = std::pow(x * x + y * y, 2. / 3.); @@ -1296,8 +1296,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { const Point<2> &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(y, -x) + numbers::PI; const double r43 = std::pow(x * x + y * y, 2. / 3.); @@ -1325,8 +1325,8 @@ namespace Functions Assert(gradients[i].size() == 1, ExcDimensionMismatch(gradients[i].size(), 1)); const Point<2> &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(y, -x) + numbers::PI; const double r43 = std::pow(x * x + y * y, 2. / 3.); @@ -1352,13 +1352,13 @@ namespace Functions { AssertIndexRange(d, 2); - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(y, -x) + numbers::PI; const double r43 = std::pow(x * x + y * y, 2. / 3.); return 2. / 3. * - (std::sin(2. / 3. * phi) * p(d) + + (std::sin(2. / 3. * phi) * p[d] + (d == 0 ? (std::cos(2. / 3. * phi) * y) : (-std::cos(2. / 3. * phi) * x))) / r43; @@ -1376,13 +1376,13 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { const Point<2> &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(y, -x) + numbers::PI; const double r43 = std::pow(x * x + y * y, 2. / 3.); values[i] = 2. / 3. * - (std::sin(2. / 3. * phi) * p(d) + + (std::sin(2. / 3. * phi) * p[d] + (d == 0 ? (std::cos(2. / 3. * phi) * y) : (-std::cos(2. / 3. * phi) * x))) / r43; @@ -1402,8 +1402,8 @@ namespace Functions { AssertDimension(values[i].size(), 2); const Point<2> &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(y, -x) + numbers::PI; const double r43 = std::pow(x * x + y * y, 2. / 3.); @@ -1473,8 +1473,8 @@ namespace Functions SlitSingularityFunction::value(const Point &p, const unsigned int) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r_squared = x * x + y * y; @@ -1495,8 +1495,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { - const double x = points[i](0); - const double y = points[i](1); + const double x = points[i][0]; + const double y = points[i][1]; const double phi = std::atan2(x, y) + numbers::PI; const double r_squared = x * x + y * y; @@ -1520,8 +1520,8 @@ namespace Functions Assert(values[i].size() == 1, ExcDimensionMismatch(values[i].size(), 1)); - const double x = points[i](0); - const double y = points[i](1); + const double x = points[i][0]; + const double y = points[i][1]; const double phi = std::atan2(x, y) + numbers::PI; const double r_squared = x * x + y * y; @@ -1560,8 +1560,8 @@ namespace Functions SlitSingularityFunction::gradient(const Point &p, const unsigned int) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r64 = std::pow(x * x + y * y, 3. / 4.); @@ -1589,8 +1589,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { const Point &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r64 = std::pow(x * x + y * y, 3. / 4.); @@ -1620,8 +1620,8 @@ namespace Functions ExcDimensionMismatch(gradients[i].size(), 1)); const Point &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r64 = std::pow(x * x + y * y, 3. / 4.); @@ -1643,8 +1643,8 @@ namespace Functions SlitHyperSingularityFunction::value(const Point<2> &p, const unsigned int) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r_squared = x * x + y * y; @@ -1663,8 +1663,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { - const double x = points[i](0); - const double y = points[i](1); + const double x = points[i][0]; + const double y = points[i][1]; const double phi = std::atan2(x, y) + numbers::PI; const double r_squared = x * x + y * y; @@ -1687,8 +1687,8 @@ namespace Functions Assert(values[i].size() == 1, ExcDimensionMismatch(values[i].size(), 1)); - const double x = points[i](0); - const double y = points[i](1); + const double x = points[i][0]; + const double y = points[i][1]; const double phi = std::atan2(x, y) + numbers::PI; const double r_squared = x * x + y * y; @@ -1724,8 +1724,8 @@ namespace Functions SlitHyperSingularityFunction::gradient(const Point<2> &p, const unsigned int) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r78 = std::pow(x * x + y * y, 7. / 8.); @@ -1753,8 +1753,8 @@ namespace Functions for (unsigned int i = 0; i < points.size(); ++i) { const Point<2> &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r78 = std::pow(x * x + y * y, 7. / 8.); @@ -1782,8 +1782,8 @@ namespace Functions ExcDimensionMismatch(gradients[i].size(), 1)); const Point<2> &p = points[i]; - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; const double phi = std::atan2(x, y) + numbers::PI; const double r78 = std::pow(x * x + y * y, 7. / 8.); @@ -1810,7 +1810,7 @@ namespace Functions angle = 0; break; case 2: - angle = std::atan2(direction(0), direction(1)); + angle = std::atan2(direction[0], direction[1]); break; case 3: Assert(false, ExcNotImplemented()); @@ -1825,7 +1825,7 @@ namespace Functions double JumpFunction::value(const Point &p, const unsigned int) const { - const double x = steepness * (-cosine * p(0) + sine * p(1)); + const double x = steepness * (-cosine * p[0] + sine * p[1]); return -std::atan(x); } @@ -1842,7 +1842,7 @@ namespace Functions for (unsigned int i = 0; i < p.size(); ++i) { - const double x = steepness * (-cosine * p[i](0) + sine * p[i](1)); + const double x = steepness * (-cosine * p[i][0] + sine * p[i][1]); values[i] = -std::atan(x); } } @@ -1852,7 +1852,7 @@ namespace Functions double JumpFunction::laplacian(const Point &p, const unsigned int) const { - const double x = steepness * (-cosine * p(0) + sine * p(1)); + const double x = steepness * (-cosine * p[0] + sine * p[1]); const double r = 1 + x * x; return 2 * steepness * steepness * x / (r * r); } @@ -1871,7 +1871,7 @@ namespace Functions for (unsigned int i = 0; i < p.size(); ++i) { - const double x = steepness * (-cosine * p[i](0) + sine * p[i](1)); + const double x = steepness * (-cosine * p[i][0] + sine * p[i][1]); const double r = 1 + x * x; values[i] = f * x / (r * r); } @@ -1883,7 +1883,7 @@ namespace Functions Tensor<1, dim> JumpFunction::gradient(const Point &p, const unsigned int) const { - const double x = steepness * (-cosine * p(0) + sine * p(1)); + const double x = steepness * (-cosine * p[0] + sine * p[1]); const double r = -steepness * (1 + x * x); Tensor<1, dim> erg; erg[0] = cosine * r; @@ -1904,7 +1904,7 @@ namespace Functions for (unsigned int i = 0; i < p.size(); ++i) { - const double x = steepness * (cosine * p[i](0) + sine * p[i](1)); + const double x = steepness * (cosine * p[i][0] + sine * p[i][1]); const double r = -steepness * (1 + x * x); gradients[i][0] = cosine * r; gradients[i][1] = sine * r; @@ -2316,8 +2316,8 @@ namespace Functions { Assert(dim == 2, ExcNotImplemented()); const double r = p.distance(center); - const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r; - const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r; + const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r; + const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r; const double dJn = (order == 0) ? @@ -2344,8 +2344,8 @@ namespace Functions { const Point &p = points[k]; const double r = p.distance(center); - const double co = (r == 0.) ? 0. : (p(0) - center(0)) / r; - const double si = (r == 0.) ? 0. : (p(1) - center(1)) / r; + const double co = (r == 0.) ? 0. : (p[0] - center[0]) / r; + const double si = (r == 0.) ? 0. : (p[1] - center[1]) / r; const double dJn = (order == 0) ? @@ -2960,8 +2960,8 @@ namespace Functions RayleighKotheVortex::vector_value(const Point &point, Vector &values) const { - const double pi_x = numbers::PI * point(0); - const double pi_y = numbers::PI * point(1); + const double pi_x = numbers::PI * point[0]; + const double pi_y = numbers::PI * point[1]; const double pi_t = numbers::PI / T * this->get_time(); values[0] = -2 * std::cos(pi_t) * diff --git a/source/base/function_restriction.cc b/source/base/function_restriction.cc index dfcd2cd87b..9542e1856e 100644 --- a/source/base/function_restriction.cc +++ b/source/base/function_restriction.cc @@ -28,13 +28,13 @@ namespace internal AssertIndexRange(component_in_dim_plus_1, dim + 1); Point output; - output(component_in_dim_plus_1) = coordinate_value; + output[component_in_dim_plus_1] = coordinate_value; for (int d = 0; d < dim; ++d) { const unsigned int component_to_write_to = dealii::internal::coordinate_to_one_dim_higher( component_in_dim_plus_1, d); - output(component_to_write_to) = point(d); + output[component_to_write_to] = point[d]; } return output; @@ -155,7 +155,7 @@ namespace Functions const unsigned int component) const { const Point full_point = - internal::create_higher_dim_point(point, open_direction, point_1D(0)); + internal::create_higher_dim_point(point, open_direction, point_1D[0]); return function->value(full_point, component); } @@ -167,7 +167,7 @@ namespace Functions const unsigned int component) const { const Point full_point = - internal::create_higher_dim_point(point, open_direction, point_1D(0)); + internal::create_higher_dim_point(point, open_direction, point_1D[0]); const Tensor<1, dim + 1> full_gradient = function->gradient(full_point, component); @@ -186,7 +186,7 @@ namespace Functions const unsigned int component) const { const Point full_point = - internal::create_higher_dim_point(point, open_direction, point_1D(0)); + internal::create_higher_dim_point(point, open_direction, point_1D[0]); const Tensor<2, dim + 1> full_hessian = function->hessian(full_point, component); diff --git a/source/base/geometric_utilities.cc b/source/base/geometric_utilities.cc index ee041f2581..0360b7d7c4 100644 --- a/source/base/geometric_utilities.cc +++ b/source/base/geometric_utilities.cc @@ -54,7 +54,7 @@ namespace GeometricUtilities { scoord[0] = position.norm(); // azimuth angle \theta: - scoord[1] = std::atan2(position(1), position(0)); + scoord[1] = std::atan2(position[1], position[0]); // correct to [0,2*pi) if (scoord[1] < 0.0) scoord[1] += 2.0 * numbers::PI; @@ -65,7 +65,7 @@ namespace GeometricUtilities { // acos returns the angle in the range [0,\pi] if (scoord[0] > std::numeric_limits::min()) - scoord[2] = std::acos(position(2) / scoord[0]); + scoord[2] = std::acos(position[2] / scoord[0]); else scoord[2] = 0.0; } diff --git a/source/base/mu_parser_internal.cc b/source/base/mu_parser_internal.cc index 7c26d31cda..04d92e4fe4 100644 --- a/source/base/mu_parser_internal.cc +++ b/source/base/mu_parser_internal.cc @@ -424,7 +424,7 @@ namespace internal init_muparser(); for (unsigned int i = 0; i < dim; ++i) - data.vars[i] = p(i); + data.vars[i] = p[i]; if (dim != this->n_vars) data.vars[dim] = time; @@ -472,7 +472,7 @@ namespace internal init_muparser(); for (unsigned int i = 0; i < dim; ++i) - data.vars[i] = p(i); + data.vars[i] = p[i]; if (dim != this->n_vars) data.vars[dim] = time; diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index 02b4004a5b..68261db079 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -66,8 +66,8 @@ namespace Polynomials for (unsigned int i = 0; i < supp.size(); ++i) if (i != center) { - lagrange_support_points.push_back(supp[i](0)); - tmp_lagrange_weight *= supp[center](0) - supp[i](0); + lagrange_support_points.push_back(supp[i][0]); + tmp_lagrange_weight *= supp[center][0] - supp[i][0]; } // check for underflow and overflow @@ -571,7 +571,7 @@ namespace Polynomials std::vector> points(n + 1); const double one_over_n = 1. / n; for (unsigned int k = 0; k <= n; ++k) - points[k](0) = static_cast(k) * one_over_n; + points[k][0] = static_cast(k) * one_over_n; return points; } } // namespace LagrangeEquidistantImplementation diff --git a/source/base/polynomial_space.cc b/source/base/polynomial_space.cc index 40f57e27c4..d99890486b 100644 --- a/source/base/polynomial_space.cc +++ b/source/base/polynomial_space.cc @@ -134,7 +134,7 @@ PolynomialSpace::compute_value(const unsigned int i, // directions double result = 1.; for (unsigned int d = 0; d < dim; ++d) - result *= polynomials[ix[d]].value(p(d)); + result *= polynomials[ix[d]].value(p[d]); return result; } @@ -155,7 +155,7 @@ PolynomialSpace::compute_grad(const unsigned int i, std::vector v(2); for (unsigned int d = 0; d < dim; ++d) { - polynomials[ix[d]].value(p(d), v); + polynomials[ix[d]].value(p[d], v); result[d] *= v[1]; for (unsigned int d1 = 0; d1 < dim; ++d1) if (d1 != d) @@ -182,7 +182,7 @@ PolynomialSpace::compute_grad_grad(const unsigned int i, std::vector v(3); for (unsigned int d = 0; d < dim; ++d) { - polynomials[ix[d]].value(p(d), v); + polynomials[ix[d]].value(p[d], v); result[d][d] *= v[2]; for (unsigned int d1 = 0; d1 < dim; ++d1) { @@ -263,7 +263,7 @@ PolynomialSpace::evaluate( for (unsigned int i = 0; i < v.size()[1]; ++i) { v(d, i).resize(v_size, 0.); - polynomials[i].value(p(d), v(d, i)); + polynomials[i].value(p[d], v(d, i)); } if (update_values) diff --git a/source/base/polynomials_abf.cc b/source/base/polynomials_abf.cc index 224cdee950..5906df4516 100644 --- a/source/base/polynomials_abf.cc +++ b/source/base/polynomials_abf.cc @@ -109,7 +109,7 @@ PolynomialsABF::evaluate( // directions. Point p; for (unsigned int c = 0; c < dim; ++c) - p(c) = unit_point((c + d) % dim); + p[c] = unit_point[(c + d) % dim]; polynomial_space.evaluate(p, p_values, diff --git a/source/base/polynomials_adini.cc b/source/base/polynomials_adini.cc index 7a1505e01c..1780e56ce8 100644 --- a/source/base/polynomials_adini.cc +++ b/source/base/polynomials_adini.cc @@ -193,8 +193,8 @@ double PolynomialsAdini::compute_value(const unsigned int i, const Point &p) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; return coef(0, i) + coef(1, i) * x + coef(2, i) * y + coef(3, i) * x * x + coef(4, i) * y * y + coef(5, i) * x * y + coef(6, i) * x * x * x + coef(7, i) * y * y * y + coef(8, i) * x * y * y + @@ -209,8 +209,8 @@ Tensor<1, dim> PolynomialsAdini::compute_grad(const unsigned int i, const Point &p) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; Tensor<1, dim> tensor; tensor[0] = dx(0, i) + dx(1, i) * x + dx(2, i) * y + dx(3, i) * x * x + dx(4, i) * y * y + dx(5, i) * x * y + dx(6, i) * x * x * x + @@ -233,8 +233,8 @@ Tensor<2, dim> PolynomialsAdini::compute_grad_grad(const unsigned int i, const Point &p) const { - const double x = p(0); - const double y = p(1); + const double x = p[0]; + const double y = p[1]; Tensor<2, dim> tensor; tensor[0][0] = dxx(0, i) + dxx(1, i) * x + dxx(2, i) * y + dxx(3, i) * x * x + dxx(4, i) * y * y + dxx(5, i) * x * y + dxx(6, i) * x * x * x + diff --git a/source/base/polynomials_bdm.cc b/source/base/polynomials_bdm.cc index b431d9718b..6a2ec18212 100644 --- a/source/base/polynomials_bdm.cc +++ b/source/base/polynomials_bdm.cc @@ -132,22 +132,22 @@ PolynomialsBDM::evaluate( else if (dim == 2) { for (unsigned int d = 0; d < dim; ++d) - monomials[0].value(unit_point(d), monovali[d]); + monomials[0].value(unit_point[d], monovali[d]); if (values.size() != 0) { values[start][0] = monovali[0][0]; - values[start][1] = -unit_point(1) * monovali[0][1]; - values[start + 1][0] = unit_point(0) * monovali[1][1]; + values[start][1] = -unit_point[1] * monovali[0][1]; + values[start + 1][0] = unit_point[0] * monovali[1][1]; values[start + 1][1] = -monovali[1][0]; } if (grads.size() != 0) { grads[start][0][0] = monovali[0][1]; grads[start][0][1] = 0.; - grads[start][1][0] = -unit_point(1) * monovali[0][2]; + grads[start][1][0] = -unit_point[1] * monovali[0][2]; grads[start][1][1] = -monovali[0][1]; grads[start + 1][0][0] = monovali[1][1]; - grads[start + 1][0][1] = unit_point(0) * monovali[1][2]; + grads[start + 1][0][1] = unit_point[0] * monovali[1][2]; grads[start + 1][1][0] = 0.; grads[start + 1][1][1] = -monovali[1][1]; } @@ -157,14 +157,14 @@ PolynomialsBDM::evaluate( grad_grads[start][0][0][1] = 0.; grad_grads[start][0][1][0] = 0.; grad_grads[start][0][1][1] = 0.; - grad_grads[start][1][0][0] = -unit_point(1) * monovali[0][3]; + grad_grads[start][1][0][0] = -unit_point[1] * monovali[0][3]; grad_grads[start][1][0][1] = -monovali[0][2]; grad_grads[start][1][1][0] = -monovali[0][2]; grad_grads[start][1][1][1] = 0.; grad_grads[start + 1][0][0][0] = 0; grad_grads[start + 1][0][0][1] = monovali[1][2]; grad_grads[start + 1][0][1][0] = monovali[1][2]; - grad_grads[start + 1][0][1][1] = unit_point(0) * monovali[1][3]; + grad_grads[start + 1][0][1][1] = unit_point[0] * monovali[1][3]; grad_grads[start + 1][1][0][0] = 0.; grad_grads[start + 1][1][0][1] = 0.; grad_grads[start + 1][1][1][0] = 0.; @@ -184,9 +184,9 @@ PolynomialsBDM::evaluate( for (unsigned int d = 0; d < dim; ++d) { // p(t) = t^(i+1) - monomials[i + 1].value(unit_point(d), monovali[d]); + monomials[i + 1].value(unit_point[d], monovali[d]); // q(t) = t^(k-i) - monomials[this->degree() - 1 - i].value(unit_point(d), + monomials[this->degree() - 1 - i].value(unit_point[d], monovalk[d]); } @@ -194,21 +194,21 @@ PolynomialsBDM::evaluate( { // x p'(y) q(z) values[start][0] = - unit_point(0) * monovali[1][1] * monovalk[2][0]; + unit_point[0] * monovali[1][1] * monovalk[2][0]; // - p(y) q(z) values[start][1] = -monovali[1][0] * monovalk[2][0]; values[start][2] = 0.; // y p'(z) q(x) values[start + 1][1] = - unit_point(1) * monovali[2][1] * monovalk[0][0]; + unit_point[1] * monovali[2][1] * monovalk[0][0]; // - p(z) q(x) values[start + 1][2] = -monovali[2][0] * monovalk[0][0]; values[start + 1][0] = 0.; // z p'(x) q(y) values[start + 2][2] = - unit_point(2) * monovali[0][1] * monovalk[1][0]; + unit_point[2] * monovali[0][1] * monovalk[1][0]; // -p(x) q(y) values[start + 2][0] = -monovali[0][0] * monovalk[1][0]; values[start + 2][1] = 0.; @@ -218,9 +218,9 @@ PolynomialsBDM::evaluate( { grads[start][0][0] = monovali[1][1] * monovalk[2][0]; grads[start][0][1] = - unit_point(0) * monovali[1][2] * monovalk[2][0]; + unit_point[0] * monovali[1][2] * monovalk[2][0]; grads[start][0][2] = - unit_point(0) * monovali[1][1] * monovalk[2][1]; + unit_point[0] * monovali[1][1] * monovalk[2][1]; grads[start][1][0] = 0.; grads[start][1][1] = -monovali[1][1] * monovalk[2][0]; grads[start][1][2] = -monovali[1][0] * monovalk[2][1]; @@ -230,9 +230,9 @@ PolynomialsBDM::evaluate( grads[start + 1][1][1] = monovali[2][1] * monovalk[0][0]; grads[start + 1][1][2] = - unit_point(1) * monovali[2][2] * monovalk[0][0]; + unit_point[1] * monovali[2][2] * monovalk[0][0]; grads[start + 1][1][0] = - unit_point(1) * monovali[2][1] * monovalk[0][1]; + unit_point[1] * monovali[2][1] * monovalk[0][1]; grads[start + 1][2][1] = 0.; grads[start + 1][2][2] = -monovali[2][1] * monovalk[0][0]; grads[start + 1][2][0] = -monovali[2][0] * monovalk[0][1]; @@ -242,9 +242,9 @@ PolynomialsBDM::evaluate( grads[start + 2][2][2] = monovali[0][1] * monovalk[1][0]; grads[start + 2][2][0] = - unit_point(2) * monovali[0][2] * monovalk[1][0]; + unit_point[2] * monovali[0][2] * monovalk[1][0]; grads[start + 2][2][1] = - unit_point(2) * monovali[0][1] * monovalk[1][1]; + unit_point[2] * monovali[0][1] * monovalk[1][1]; grads[start + 2][0][2] = 0.; grads[start + 2][0][0] = -monovali[0][1] * monovalk[1][0]; grads[start + 2][0][1] = -monovali[0][0] * monovalk[1][1]; @@ -260,14 +260,14 @@ PolynomialsBDM::evaluate( grad_grads[start][0][0][2] = monovali[1][1] * monovalk[2][1]; grad_grads[start][0][1][0] = monovali[1][2] * monovalk[2][0]; grad_grads[start][0][1][1] = - unit_point(0) * monovali[1][3] * monovalk[2][0]; + unit_point[0] * monovali[1][3] * monovalk[2][0]; grad_grads[start][0][1][2] = - unit_point(0) * monovali[1][2] * monovalk[2][1]; + unit_point[0] * monovali[1][2] * monovalk[2][1]; grad_grads[start][0][2][0] = monovali[1][1] * monovalk[2][1]; grad_grads[start][0][2][1] = - unit_point(0) * monovali[1][2] * monovalk[2][1]; + unit_point[0] * monovali[1][2] * monovalk[2][1]; grad_grads[start][0][2][2] = - unit_point(0) * monovali[1][1] * monovalk[2][2]; + unit_point[0] * monovali[1][1] * monovalk[2][2]; grad_grads[start][1][0][0] = 0.; grad_grads[start][1][0][1] = 0.; grad_grads[start][1][0][2] = 0.; @@ -297,18 +297,18 @@ PolynomialsBDM::evaluate( grad_grads[start + 1][0][2][1] = 0.; grad_grads[start + 1][0][2][2] = 0.; grad_grads[start + 1][1][0][0] = - unit_point(1) * monovali[2][1] * monovalk[0][2]; + unit_point[1] * monovali[2][1] * monovalk[0][2]; grad_grads[start + 1][1][0][1] = monovali[2][1] * monovalk[0][1]; grad_grads[start + 1][1][0][2] = - unit_point(1) * monovali[2][2] * monovalk[0][1]; + unit_point[1] * monovali[2][2] * monovalk[0][1]; grad_grads[start + 1][1][1][0] = monovalk[0][1] * monovali[2][1]; grad_grads[start + 1][1][1][1] = 0.; grad_grads[start + 1][1][1][2] = monovalk[0][0] * monovali[2][2]; grad_grads[start + 1][1][2][0] = - unit_point(1) * monovalk[0][1] * monovali[2][2]; + unit_point[1] * monovalk[0][1] * monovali[2][2]; grad_grads[start + 1][1][2][1] = monovalk[0][0] * monovali[2][2]; grad_grads[start + 1][1][2][2] = - unit_point(1) * monovalk[0][0] * monovali[2][3]; + unit_point[1] * monovalk[0][0] * monovali[2][3]; grad_grads[start + 1][2][0][0] = -monovalk[0][2] * monovali[2][0]; grad_grads[start + 1][2][0][1] = 0.; grad_grads[start + 1][2][0][2] = -monovalk[0][1] * monovali[2][1]; @@ -338,14 +338,14 @@ PolynomialsBDM::evaluate( grad_grads[start + 2][1][2][1] = 0.; grad_grads[start + 2][1][2][2] = 0.; grad_grads[start + 2][2][0][0] = - unit_point(2) * monovali[0][3] * monovalk[1][0]; + unit_point[2] * monovali[0][3] * monovalk[1][0]; grad_grads[start + 2][2][0][1] = - unit_point(2) * monovali[0][2] * monovalk[1][1]; + unit_point[2] * monovali[0][2] * monovalk[1][1]; grad_grads[start + 2][2][0][2] = monovali[0][2] * monovalk[1][0]; grad_grads[start + 2][2][1][0] = - unit_point(2) * monovali[0][2] * monovalk[1][1]; + unit_point[2] * monovali[0][2] * monovalk[1][1]; grad_grads[start + 2][2][1][1] = - unit_point(2) * monovali[0][1] * monovalk[1][2]; + unit_point[2] * monovali[0][1] * monovalk[1][2]; grad_grads[start + 2][2][1][2] = monovali[0][1] * monovalk[1][1]; grad_grads[start + 2][2][2][0] = monovali[0][2] * monovalk[1][0]; grad_grads[start + 2][2][2][1] = monovali[0][1] * monovalk[1][1]; diff --git a/source/base/polynomials_nedelec.cc b/source/base/polynomials_nedelec.cc index 253eb7f8fd..d5a824921b 100644 --- a/source/base/polynomials_nedelec.cc +++ b/source/base/polynomials_nedelec.cc @@ -133,8 +133,8 @@ PolynomialsNedelec::evaluate( // shifted one step in positive direction Point p; - p(0) = unit_point(1); - p(1) = unit_point(0); + p[0] = unit_point[1]; + p[1] = unit_point[0]; std::vector p_values((values.empty()) ? 0 : n_basis); std::vector> p_grads((grads.empty()) ? 0 : n_basis); @@ -327,18 +327,18 @@ PolynomialsNedelec::evaluate( std::vector> p2_grad_grads( (grad_grads.empty()) ? 0 : n_basis); - p1(0) = unit_point(1); - p1(1) = unit_point(2); - p1(2) = unit_point(0); + p1[0] = unit_point[1]; + p1[1] = unit_point[2]; + p1[2] = unit_point[0]; polynomial_space.evaluate(p1, p1_values, p1_grads, p1_grad_grads, empty_vector_of_3rd_order_tensors, empty_vector_of_4th_order_tensors); - p2(0) = unit_point(2); - p2(1) = unit_point(0); - p2(2) = unit_point(1); + p2[0] = unit_point[2]; + p2[1] = unit_point[0]; + p2[2] = unit_point[1]; polynomial_space.evaluate(p2, p2_values, p2_grads, diff --git a/source/base/polynomials_rannacher_turek.cc b/source/base/polynomials_rannacher_turek.cc index d85986634d..860732821a 100644 --- a/source/base/polynomials_rannacher_turek.cc +++ b/source/base/polynomials_rannacher_turek.cc @@ -39,23 +39,23 @@ PolynomialsRannacherTurek::compute_value(const unsigned int i, Assert(dim == 2, ExcNotImplemented()); if (i == 0) { - return (0.75 - 2.5 * p(0) + 1.5 * p(1) + - 1.5 * (p(0) * p(0) - p(1) * p(1))); + return (0.75 - 2.5 * p[0] + 1.5 * p[1] + + 1.5 * (p[0] * p[0] - p[1] * p[1])); } else if (i == 1) { - return (-0.25 - 0.5 * p(0) + 1.5 * p(1) + - 1.5 * (p(0) * p(0) - p(1) * p(1))); + return (-0.25 - 0.5 * p[0] + 1.5 * p[1] + + 1.5 * (p[0] * p[0] - p[1] * p[1])); } else if (i == 2) { - return (0.75 + 1.5 * p(0) - 2.5 * p(1) - - 1.5 * (p(0) * p(0) - p(1) * p(1))); + return (0.75 + 1.5 * p[0] - 2.5 * p[1] - + 1.5 * (p[0] * p[0] - p[1] * p[1])); } else if (i == 3) { - return (-0.25 + 1.5 * p(0) - 0.5 * p(1) - - 1.5 * (p(0) * p(0) - p(1) * p(1))); + return (-0.25 + 1.5 * p[0] - 0.5 * p[1] - + 1.5 * (p[0] * p[0] - p[1] * p[1])); } Assert(false, ExcNotImplemented()); @@ -74,23 +74,23 @@ PolynomialsRannacherTurek::compute_grad(const unsigned int i, Tensor<1, dim> grad; if (i == 0) { - grad[0] = -2.5 + 3 * p(0); - grad[1] = 1.5 - 3 * p(1); + grad[0] = -2.5 + 3 * p[0]; + grad[1] = 1.5 - 3 * p[1]; } else if (i == 1) { - grad[0] = -0.5 + 3.0 * p(0); - grad[1] = 1.5 - 3.0 * p(1); + grad[0] = -0.5 + 3.0 * p[0]; + grad[1] = 1.5 - 3.0 * p[1]; } else if (i == 2) { - grad[0] = 1.5 - 3.0 * p(0); - grad[1] = -2.5 + 3.0 * p(1); + grad[0] = 1.5 - 3.0 * p[0]; + grad[1] = -2.5 + 3.0 * p[1]; } else if (i == 3) { - grad[0] = 1.5 - 3.0 * p(0); - grad[1] = -0.5 + 3.0 * p(1); + grad[0] = 1.5 - 3.0 * p[0]; + grad[1] = -0.5 + 3.0 * p[1]; } else { diff --git a/source/base/polynomials_raviart_thomas.cc b/source/base/polynomials_raviart_thomas.cc index 56c8cc56ae..2952983a81 100644 --- a/source/base/polynomials_raviart_thomas.cc +++ b/source/base/polynomials_raviart_thomas.cc @@ -172,7 +172,7 @@ PolynomialsRaviartThomas::evaluate( // through all directions. Point p; for (unsigned int c = 0; c < dim; ++c) - p(c) = unit_point((c + d) % dim); + p[c] = unit_point[(c + d) % dim]; polynomial_space.evaluate(p, p_values, diff --git a/source/base/polynomials_rt_bubbles.cc b/source/base/polynomials_rt_bubbles.cc index 1a18565738..3d27437bec 100644 --- a/source/base/polynomials_rt_bubbles.cc +++ b/source/base/polynomials_rt_bubbles.cc @@ -126,14 +126,14 @@ PolynomialsRT_Bubbles::evaluate( // monoval_i = x^i, // monoval_plus = x^(k+1) for (unsigned int d = 0; d < dim; ++d) - monomials[my_degree + 1].value(unit_point(d), + monomials[my_degree + 1].value(unit_point[d], n_derivatives, monoval_plus[d]); for (unsigned int i = 0; i <= my_degree; ++i, ++start) { for (unsigned int d = 0; d < dim; ++d) - monomials[i].value(unit_point(d), n_derivatives, monoval_i[d]); + monomials[i].value(unit_point[d], n_derivatives, monoval_i[d]); if (values.size() != 0) { @@ -216,10 +216,10 @@ PolynomialsRT_Bubbles::evaluate( // monoval_* = x^*, monoval_jplus = x^(j+1) for (unsigned int d = 0; d < dim; ++d) { - monomials[my_degree + 1].value(unit_point(d), + monomials[my_degree + 1].value(unit_point[d], n_derivatives, monoval_plus[d]); - monomials[my_degree].value(unit_point(d), n_derivatives, monoval[d]); + monomials[my_degree].value(unit_point[d], n_derivatives, monoval[d]); } const unsigned int n_curls = (my_degree + 1) * (2 * my_degree + 1); @@ -227,16 +227,16 @@ PolynomialsRT_Bubbles::evaluate( for (unsigned int i = 0; i <= my_degree; ++i) { for (unsigned int d = 0; d < dim; ++d) - monomials[i].value(unit_point(d), n_derivatives, monoval_i[d]); + monomials[i].value(unit_point[d], n_derivatives, monoval_i[d]); for (unsigned int j = 0; j <= my_degree; ++j) { for (unsigned int d = 0; d < dim; ++d) { - monomials[j].value(unit_point(d), + monomials[j].value(unit_point[d], n_derivatives, monoval_j[d]); - monomials[j + 1].value(unit_point(d), + monomials[j + 1].value(unit_point[d], n_derivatives, monoval_jplus[d]); } diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 5203a7eb85..9114743a6a 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -138,8 +138,8 @@ namespace internal for (const Point<2> &p : points) { Point<3> cell_point; - cell_point[xi_index] = xi_scale * p(0) + xi_translation; - cell_point[eta_index] = eta_scale * p(1) + eta_translation; + cell_point[xi_index] = xi_scale * p[0] + xi_translation; + cell_point[eta_index] = eta_scale * p[1] + eta_translation; cell_point[const_index] = const_value; q_points.push_back(cell_point); } @@ -376,14 +376,14 @@ QProjector<2>::project_to_face(const ReferenceCell &reference_cell, switch (face_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0), 0); + q_points[p] = Point(quadrature.point(p)[0], 0); break; case 1: q_points[p] = - Point(1 - quadrature.point(p)(0), quadrature.point(p)(0)); + Point(1 - quadrature.point(p)[0], quadrature.point(p)[0]); break; case 2: - q_points[p] = Point(0, 1 - quadrature.point(p)(0)); + q_points[p] = Point(0, 1 - quadrature.point(p)[0]); break; default: Assert(false, ExcInternalError()); @@ -395,16 +395,16 @@ QProjector<2>::project_to_face(const ReferenceCell &reference_cell, switch (face_no) { case 0: - q_points[p] = Point(0, quadrature.point(p)(0)); + q_points[p] = Point(0, quadrature.point(p)[0]); break; case 1: - q_points[p] = Point(1, quadrature.point(p)(0)); + q_points[p] = Point(1, quadrature.point(p)[0]); break; case 2: - q_points[p] = Point(quadrature.point(p)(0), 0); + q_points[p] = Point(quadrature.point(p)[0], 0); break; case 3: - q_points[p] = Point(quadrature.point(p)(0), 1); + q_points[p] = Point(quadrature.point(p)[0], 1); break; default: Assert(false, ExcInternalError()); @@ -519,11 +519,11 @@ QProjector<2>::project_to_subface(const ReferenceCell &reference_cell, switch (subface_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0) / 2, 0); + q_points[p] = Point(quadrature.point(p)[0] / 2, 0); break; case 1: q_points[p] = - Point(0.5 + quadrature.point(p)(0) / 2, 0); + Point(0.5 + quadrature.point(p)[0] / 2, 0); break; default: Assert(false, ExcInternalError()); @@ -533,12 +533,12 @@ QProjector<2>::project_to_subface(const ReferenceCell &reference_cell, switch (subface_no) { case 0: - q_points[p] = Point(1 - quadrature.point(p)(0) / 2, - quadrature.point(p)(0) / 2); + q_points[p] = Point(1 - quadrature.point(p)[0] / 2, + quadrature.point(p)[0] / 2); break; case 1: - q_points[p] = Point(0.5 - quadrature.point(p)(0) / 2, - 0.5 + quadrature.point(p)(0) / 2); + q_points[p] = Point(0.5 - quadrature.point(p)[0] / 2, + 0.5 + quadrature.point(p)[0] / 2); break; default: Assert(false, ExcInternalError()); @@ -548,11 +548,11 @@ QProjector<2>::project_to_subface(const ReferenceCell &reference_cell, switch (subface_no) { case 0: - q_points[p] = Point(0, 1 - quadrature.point(p)(0) / 2); + q_points[p] = Point(0, 1 - quadrature.point(p)[0] / 2); break; case 1: q_points[p] = - Point(0, 0.5 - quadrature.point(p)(0) / 2); + Point(0, 0.5 - quadrature.point(p)[0] / 2); break; default: Assert(false, ExcInternalError()); @@ -571,11 +571,11 @@ QProjector<2>::project_to_subface(const ReferenceCell &reference_cell, switch (subface_no) { case 0: - q_points[p] = Point(0, quadrature.point(p)(0) / 2); + q_points[p] = Point(0, quadrature.point(p)[0] / 2); break; case 1: q_points[p] = - Point(0, quadrature.point(p)(0) / 2 + 0.5); + Point(0, quadrature.point(p)[0] / 2 + 0.5); break; default: Assert(false, ExcInternalError()); @@ -585,11 +585,11 @@ QProjector<2>::project_to_subface(const ReferenceCell &reference_cell, switch (subface_no) { case 0: - q_points[p] = Point(1, quadrature.point(p)(0) / 2); + q_points[p] = Point(1, quadrature.point(p)[0] / 2); break; case 1: q_points[p] = - Point(1, quadrature.point(p)(0) / 2 + 0.5); + Point(1, quadrature.point(p)[0] / 2 + 0.5); break; default: Assert(false, ExcInternalError()); @@ -599,11 +599,11 @@ QProjector<2>::project_to_subface(const ReferenceCell &reference_cell, switch (subface_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0) / 2, 0); + q_points[p] = Point(quadrature.point(p)[0] / 2, 0); break; case 1: q_points[p] = - Point(quadrature.point(p)(0) / 2 + 0.5, 0); + Point(quadrature.point(p)[0] / 2 + 0.5, 0); break; default: Assert(false, ExcInternalError()); @@ -613,11 +613,11 @@ QProjector<2>::project_to_subface(const ReferenceCell &reference_cell, switch (subface_no) { case 0: - q_points[p] = Point(quadrature.point(p)(0) / 2, 1); + q_points[p] = Point(quadrature.point(p)[0] / 2, 1); break; case 1: q_points[p] = - Point(quadrature.point(p)(0) / 2 + 0.5, 1); + Point(quadrature.point(p)[0] / 2 + 0.5, 1); break; default: Assert(false, ExcInternalError()); @@ -1304,7 +1304,7 @@ QProjector::project_to_line(const ReferenceCell &reference_cell, for (unsigned int k = 0; k < n; ++k) { - const double alpha = quadrature.point(k)(0); + const double alpha = quadrature.point(k)[0]; points[k] = alpha * p2; points[k] += (1. - alpha) * p1; weights[k] = length * quadrature.weight(k); diff --git a/source/base/quadrature.cc b/source/base/quadrature.cc index a3daba9958..d09401f0f6 100644 --- a/source/base/quadrature.cc +++ b/source/base/quadrature.cc @@ -155,8 +155,8 @@ Quadrature::Quadrature(const SubQuadrature &q1, const Quadrature<1> &q2) // compose coordinates of new quadrature point by tensor product in the // last component for (unsigned int d = 0; d < dim - 1; ++d) - quadrature_points[present_index](d) = q1.point(i1)(d); - quadrature_points[present_index](dim - 1) = q2.point(i2)(0); + quadrature_points[present_index][d] = q1.point(i1)[d]; + quadrature_points[present_index][dim - 1] = q2.point(i2)[0]; weights[present_index] = q1.weight(i1) * q2.weight(i2); @@ -197,7 +197,7 @@ Quadrature<1>::Quadrature(const SubQuadrature &, const Quadrature<1> &q2) { // compose coordinates of new quadrature point by tensor product in the // last component - quadrature_points[present_index](0) = q2.point(i2)(0); + quadrature_points[present_index][0] = q2.point(i2)[0]; weights[present_index] = q2.weight(i2); @@ -258,11 +258,11 @@ Quadrature::Quadrature(const Quadrature &q) for (unsigned int i1 = 0; i1 < n1; ++i1) for (unsigned int i0 = 0; i0 < n0; ++i0) { - quadrature_points[k](0) = q.point(i0)(0); + quadrature_points[k][0] = q.point(i0)[0]; if (dim > 1) - quadrature_points[k](1) = q.point(i1)(0); + quadrature_points[k][1] = q.point(i1)[0]; if (dim > 2) - quadrature_points[k](2) = q.point(i2)(0); + quadrature_points[k][2] = q.point(i2)[0]; weights[k] = q.weight(i0); if (dim > 1) weights[k] *= q.weight(i1); @@ -370,7 +370,7 @@ QAnisotropic::QAnisotropic(const Quadrature<1> &qx) unsigned int k = 0; for (unsigned int k1 = 0; k1 < qx.size(); ++k1) { - this->quadrature_points[k](0) = qx.point(k1)(0); + this->quadrature_points[k][0] = qx.point(k1)[0]; this->weights[k++] = qx.weight(k1); } Assert(k == this->size(), ExcInternalError()); @@ -393,8 +393,8 @@ QAnisotropic::QAnisotropic(const Quadrature<1> &qx, for (unsigned int k2 = 0; k2 < qy.size(); ++k2) for (unsigned int k1 = 0; k1 < qx.size(); ++k1) { - this->quadrature_points[k](0) = qx.point(k1)(0); - this->quadrature_points[k](dim_1) = qy.point(k2)(0); + this->quadrature_points[k][0] = qx.point(k1)[0]; + this->quadrature_points[k][dim_1] = qy.point(k2)[0]; this->weights[k++] = qx.weight(k1) * qy.weight(k2); } Assert(k == this->size(), ExcInternalError()); @@ -424,9 +424,9 @@ QAnisotropic::QAnisotropic(const Quadrature<1> &qx, for (unsigned int k2 = 0; k2 < qy.size(); ++k2) for (unsigned int k1 = 0; k1 < qx.size(); ++k1) { - this->quadrature_points[k](0) = qx.point(k1)(0); - this->quadrature_points[k](dim_1) = qy.point(k2)(0); - this->quadrature_points[k](dim_2) = qz.point(k3)(0); + this->quadrature_points[k][0] = qx.point(k1)[0]; + this->quadrature_points[k][dim_1] = qy.point(k2)[0]; + this->quadrature_points[k][dim_2] = qz.point(k3)[0]; this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3); } Assert(k == this->size(), ExcInternalError()); @@ -522,7 +522,7 @@ QIterated<1>::QIterated(const Quadrature<1> &base_quadrature, ++q_point) { this->quadrature_points[next_point] = - Point<1>(base_quadrature.point(q_point)(0) * + Point<1>(base_quadrature.point(q_point)[0] * (intervals[copy + 1][0] - intervals[copy][0]) + intervals[copy][0]); this->weights[next_point] = @@ -564,7 +564,7 @@ QIterated<1>::QIterated(const Quadrature<1> &base_quadrature, { Assert(this->quadrature_points[next_point - double_point_offset] .distance(Point<1>( - base_quadrature.point(q_point)(0) * + base_quadrature.point(q_point)[0] * (intervals[copy + 1][0] - intervals[copy][0]) + intervals[copy][0])) < 1e-10 /*tolerance*/, ExcInternalError()); @@ -577,7 +577,7 @@ QIterated<1>::QIterated(const Quadrature<1> &base_quadrature, } this->quadrature_points[next_point] = - Point<1>(base_quadrature.point(q_point)(0) * + Point<1>(base_quadrature.point(q_point)[0] * (intervals[copy + 1][0] - intervals[copy][0]) + intervals[copy][0]); diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc index 51eabd0f6d..d5ebf56490 100644 --- a/source/base/tensor_product_polynomials.cc +++ b/source/base/tensor_product_polynomials.cc @@ -172,7 +172,7 @@ TensorProductPolynomials::compute_value( double value = 1.; for (unsigned int d = 0; d < dim; ++d) - value *= polynomials[indices[d]].value(p(d)); + value *= polynomials[indices[d]].value(p[d]); return value; } @@ -210,7 +210,7 @@ TensorProductPolynomials::compute_grad( std::vector tmp(2); for (unsigned int d = 0; d < dim; ++d) { - polynomials[indices[d]].value(p(d), tmp); + polynomials[indices[d]].value(p[d], tmp); v[d][0] = tmp[0]; v[d][1] = tmp[1]; } @@ -256,7 +256,7 @@ TensorProductPolynomials::compute_grad_grad( std::vector tmp(3); for (unsigned int d = 0; d < dim; ++d) { - polynomials[indices[d]].value(p(d), tmp); + polynomials[indices[d]].value(p[d], tmp); v[d][0] = tmp[0]; v[d][1] = tmp[1]; v[d][2] = tmp[2]; @@ -707,7 +707,7 @@ AnisotropicPolynomials::compute_value(const unsigned int i, double value = 1.; for (unsigned int d = 0; d < dim; ++d) - value *= polynomials[d][indices[d]].value(p(d)); + value *= polynomials[d][indices[d]].value(p[d]); return value; } @@ -740,7 +740,7 @@ AnisotropicPolynomials::compute_grad(const unsigned int i, // coordinate direction ndarray v; for (unsigned int d = 0; d < dim; ++d) - polynomials[d][indices[d]].value(p(d), 1, v[d].data()); + polynomials[d][indices[d]].value(p[d], 1, v[d].data()); Tensor<1, dim> grad; for (unsigned int d = 0; d < dim; ++d) @@ -777,7 +777,7 @@ AnisotropicPolynomials::compute_grad_grad(const unsigned int i, ndarray v; for (unsigned int d = 0; d < dim; ++d) - polynomials[d][indices[d]].value(p(d), 2, v[d].data()); + polynomials[d][indices[d]].value(p[d], 2, v[d].data()); Tensor<2, dim> grad_grad; for (unsigned int d1 = 0; d1 < dim; ++d1) diff --git a/source/base/tensor_product_polynomials_bubbles.cc b/source/base/tensor_product_polynomials_bubbles.cc index 46fe84f8b6..84209e2990 100644 --- a/source/base/tensor_product_polynomials_bubbles.cc +++ b/source/base/tensor_product_polynomials_bubbles.cc @@ -83,12 +83,12 @@ TensorProductPolynomialsBubbles::compute_value(const unsigned int i, // Compute \prod_{i=1}^d 4*x_i*(1-x_i) double value = 1.; for (unsigned int j = 0; j < dim; ++j) - value *= 4 * p(j) * (1 - p(j)); + value *= 4 * p[j] * (1 - p[j]); // Then multiply with (2x_i-1)^{r-1}. Since q_degree is generally a // small integer, using a loop is likely faster than using std::pow. for (unsigned int i = 0; i < q_degree - 1; ++i) - value *= (2 * p(comp) - 1); + value *= (2 * p[comp] - 1); return value; } @@ -127,10 +127,10 @@ TensorProductPolynomialsBubbles::compute_grad(const unsigned int i, grad[d] = 1.; // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p) for (unsigned j = 0; j < dim; ++j) - grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j))); + grad[d] *= (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j])); // and multiply with (2*x_i-1)^{r-1} for (unsigned int i = 0; i < q_degree - 1; ++i) - grad[d] *= 2 * p(comp) - 1; + grad[d] *= 2 * p[comp] - 1; } if (q_degree >= 2) @@ -138,11 +138,11 @@ TensorProductPolynomialsBubbles::compute_grad(const unsigned int i, // add \prod_{i=1}^d 4*(x_i(1-x_i))(p) double value = 1.; for (unsigned int j = 0; j < dim; ++j) - value *= 4 * p(j) * (1 - p(j)); + value *= 4 * p[j] * (1 - p[j]); // and multiply with grad(2*x_i-1)^{r-1} double tmp = value * 2 * (q_degree - 1); for (unsigned int i = 0; i < q_degree - 2; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; grad[comp] += tmp; } @@ -172,21 +172,21 @@ TensorProductPolynomialsBubbles::compute_grad_grad( { for (unsigned int c = 0; c < dim; ++c) { - v[c][0] = 4 * p(c) * (1 - p(c)); - v[c][1] = 4 * (1 - 2 * p(c)); + v[c][0] = 4 * p[c] * (1 - p[c]); + v[c][1] = 4 * (1 - 2 * p[c]); v[c][2] = -8; } double tmp = 1.; for (unsigned int i = 0; i < q_degree - 1; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; v[dim][0] = tmp; if (q_degree >= 2) { double tmp = 2 * (q_degree - 1); for (unsigned int i = 0; i < q_degree - 2; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; v[dim][1] = tmp; } else @@ -196,7 +196,7 @@ TensorProductPolynomialsBubbles::compute_grad_grad( { double tmp = 4 * (q_degree - 2) * (q_degree - 1); for (unsigned int i = 0; i < q_degree - 3; ++i) - tmp *= 2 * p(comp) - 1; + tmp *= 2 * p[comp] - 1; v[dim][2] = tmp; } else diff --git a/source/fe/fe_dgp_monomial.cc b/source/fe/fe_dgp_monomial.cc index db48efebd4..66fd1f5776 100644 --- a/source/fe/fe_dgp_monomial.cc +++ b/source/fe/fe_dgp_monomial.cc @@ -91,7 +91,7 @@ namespace internal Assert(p.size() == k + 1, ExcDimensionMismatch(p.size(), k + 1)); const double h = 1. / k; for (unsigned int i = 0; i < p.size(); ++i) - p[i](0) = i * h; + p[i][0] = i * h; } template <> @@ -103,8 +103,8 @@ namespace internal ExcInternalError()); for (unsigned int i = 0; i < p.size(); ++i) { - p[i](0) = points2d[start_index2d[k] + i][0]; - p[i](1) = points2d[start_index2d[k] + i][1]; + p[i][0] = points2d[start_index2d[k] + i][0]; + p[i][1] = points2d[start_index2d[k] + i][1]; } } @@ -117,9 +117,9 @@ namespace internal ExcInternalError()); for (unsigned int i = 0; i < p.size(); ++i) { - p[i](0) = points3d[start_index3d[k] + i][0]; - p[i](1) = points3d[start_index3d[k] + i][1]; - p[i](2) = points3d[start_index3d[k] + i][2]; + p[i][0] = points3d[start_index3d[k] + i][0]; + p[i][1] = points3d[start_index3d[k] + i][1]; + p[i][2] = points3d[start_index3d[k] + i][2]; } } } // namespace diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index 7c7db61e47..3f24eb3875 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -890,7 +890,7 @@ FE_DGQArbitraryNodes::get_name() const const QGaussLobatto<1> points_gl(this->degree + 1); bool gauss_lobatto = true; for (unsigned int j = 0; j <= this->degree; ++j) - if (points[j] != points_gl.point(j)(0)) + if (points[j] != points_gl.point(j)[0]) { gauss_lobatto = false; break; @@ -907,7 +907,7 @@ FE_DGQArbitraryNodes::get_name() const const QGauss<1> points_g(this->degree + 1); bool gauss = true; for (unsigned int j = 0; j <= this->degree; ++j) - if (points[j] != points_g.point(j)(0)) + if (points[j] != points_g.point(j)[0]) { gauss = false; break; @@ -924,7 +924,7 @@ FE_DGQArbitraryNodes::get_name() const const QGaussLog<1> points_glog(this->degree + 1); bool gauss_log = true; for (unsigned int j = 0; j <= this->degree; ++j) - if (points[j] != points_glog.point(j)(0)) + if (points[j] != points_glog.point(j)[0]) { gauss_log = false; break; diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc index a0667bac5b..96d4070266 100644 --- a/source/fe/fe_face.cc +++ b/source/fe/fe_face.cc @@ -79,11 +79,11 @@ FE_FaceQ::FE_FaceQ(const unsigned int degree) { Point p; - p(0) = points[ix][0]; + p[0] = points[ix][0]; if (codim > 1) - p(1) = points[iy][0]; + p[1] = points[iy][0]; if (codim > 2) - p(2) = points[iz][0]; + p[2] = points[iz][0]; this->unit_face_support_points[0][k++] = p; } diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index 674082578b..3e5db68626 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -338,7 +338,7 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order) boundary_weights(q_point, i) = reference_edge_quadrature.weight(q_point) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[face_no][q_point](0)); + this->generalized_face_support_points[face_no][q_point][0]); } for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point) @@ -441,20 +441,20 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) { this->generalized_support_points[q_point + (i + 4 * j) * n_edge_points] = - Point(i, reference_edge_quadrature.point(q_point)(0), j); + Point(i, reference_edge_quadrature.point(q_point)[0], j); this->generalized_support_points[q_point + (i + 4 * j + 2) * n_edge_points] = - Point(reference_edge_quadrature.point(q_point)(0), i, j); + Point(reference_edge_quadrature.point(q_point)[0], i, j); this->generalized_support_points[q_point + (i + 2 * (j + 4)) * n_edge_points] = - Point(i, j, reference_edge_quadrature.point(q_point)(0)); + Point(i, j, reference_edge_quadrature.point(q_point)[0]); } for (unsigned int i = 0; i < order; ++i) boundary_weights(q_point, i) = reference_edge_quadrature.weight(q_point) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[face_no][q_point](1)); + this->generalized_face_support_points[face_no][q_point][1]); } // Create support points on faces. @@ -471,19 +471,19 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) reference_face_quadrature.weight(q_point) * lobatto_polynomials_grad[i].value( this->generalized_face_support_points - [face_no][q_point + 4 * n_edge_points](0)) * + [face_no][q_point + 4 * n_edge_points][0]) * lobatto_polynomials[j + 2].value( this->generalized_face_support_points - [face_no][q_point + 4 * n_edge_points](1)); + [face_no][q_point + 4 * n_edge_points][1]); boundary_weights(q_point + n_edge_points, 2 * (i * order + j) + 1) = reference_face_quadrature.weight(q_point) * lobatto_polynomials_grad[i].value( this->generalized_face_support_points - [face_no][q_point + 4 * n_edge_points](1)) * + [face_no][q_point + 4 * n_edge_points][1]) * lobatto_polynomials[j + 2].value( this->generalized_face_support_points - [face_no][q_point + 4 * n_edge_points](0)); + [face_no][q_point + 4 * n_edge_points][0]); } } @@ -542,13 +542,13 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) { this->generalized_support_points[q_point + (i + 4 * j) * n_edge_points] = - Point(i, reference_edge_quadrature.point(q_point)(0), j); + Point(i, reference_edge_quadrature.point(q_point)[0], j); this->generalized_support_points[q_point + (i + 4 * j + 2) * n_edge_points] = - Point(reference_edge_quadrature.point(q_point)(0), i, j); + Point(reference_edge_quadrature.point(q_point)[0], i, j); this->generalized_support_points[q_point + (i + 2 * (j + 4)) * n_edge_points] = - Point(i, j, reference_edge_quadrature.point(q_point)(0)); + Point(i, j, reference_edge_quadrature.point(q_point)[0]); } } } @@ -600,24 +600,24 @@ FE_Nedelec::initialize_restriction() { const double weight = 2.0 * edge_quadrature.weight(q_point); - if (edge_quadrature_points[q_point](0) < 0.5) + if (edge_quadrature_points[q_point][0] < 0.5) { Point quadrature_point( - 0.0, 2.0 * edge_quadrature_points[q_point](0)); + 0.0, 2.0 * edge_quadrature_points[q_point][0]); this->restriction[index][0](0, dof) += weight * this->shape_value_component(dof, quadrature_point, 1); - quadrature_point(0) = 1.0; + quadrature_point[0] = 1.0; this->restriction[index][1](this->degree, dof) += weight * this->shape_value_component(dof, quadrature_point, 1); - quadrature_point(0) = quadrature_point(1); - quadrature_point(1) = 0.0; + quadrature_point[0] = quadrature_point[1]; + quadrature_point[1] = 0.0; this->restriction[index][0](2 * this->degree, dof) += weight * this->shape_value_component(dof, quadrature_point, 0); - quadrature_point(1) = 1.0; + quadrature_point[1] = 1.0; this->restriction[index][2](3 * this->degree, dof) += weight * this->shape_value_component(dof, quadrature_point, 0); @@ -626,21 +626,21 @@ FE_Nedelec::initialize_restriction() else { Point quadrature_point( - 0.0, 2.0 * edge_quadrature_points[q_point](0) - 1.0); + 0.0, 2.0 * edge_quadrature_points[q_point][0] - 1.0); this->restriction[index][2](0, dof) += weight * this->shape_value_component(dof, quadrature_point, 1); - quadrature_point(0) = 1.0; + quadrature_point[0] = 1.0; this->restriction[index][3](this->degree, dof) += weight * this->shape_value_component(dof, quadrature_point, 1); - quadrature_point(0) = quadrature_point(1); - quadrature_point(1) = 0.0; + quadrature_point[0] = quadrature_point[1]; + quadrature_point[1] = 0.0; this->restriction[index][1](2 * this->degree, dof) += weight * this->shape_value_component(dof, quadrature_point, 0); - quadrature_point(1) = 1.0; + quadrature_point[1] = 1.0; this->restriction[index][3](3 * this->degree, dof) += weight * this->shape_value_component(dof, quadrature_point, 0); @@ -673,7 +673,7 @@ FE_Nedelec::initialize_restriction() for (unsigned int i = 0; i < deg; ++i) assembling_matrix(i, q_point) = weight * legendre_polynomials[i + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); } FullMatrix system_matrix(deg, deg); @@ -697,14 +697,14 @@ FE_Nedelec::initialize_restriction() { const double weight = edge_quadrature.weight(q_point); const Point quadrature_point_0( - i, edge_quadrature_points[q_point](0)); + i, edge_quadrature_points[q_point][0]); const Point quadrature_point_1( - edge_quadrature_points[q_point](0), i); + edge_quadrature_points[q_point][0], i); - if (edge_quadrature_points[q_point](0) < 0.5) + if (edge_quadrature_points[q_point][0] < 0.5) { Point quadrature_point_2( - i, 2.0 * edge_quadrature_points[q_point](0)); + i, 2.0 * edge_quadrature_points[q_point][0]); tmp(0) = weight * @@ -723,7 +723,7 @@ FE_Nedelec::initialize_restriction() quadrature_point_0, 1); quadrature_point_2 = Point( - 2.0 * edge_quadrature_points[q_point](0), i); + 2.0 * edge_quadrature_points[q_point][0], i); tmp(2) = weight * (2.0 * this->shape_value_component( @@ -755,7 +755,7 @@ FE_Nedelec::initialize_restriction() Point quadrature_point_2( i, - 2.0 * edge_quadrature_points[q_point](0) - 1.0); + 2.0 * edge_quadrature_points[q_point][0] - 1.0); tmp(1) = weight * @@ -774,7 +774,7 @@ FE_Nedelec::initialize_restriction() this->shape_value_component( (i + 2) * this->degree, quadrature_point_1, 0); quadrature_point_2 = Point( - 2.0 * edge_quadrature_points[q_point](0) - 1.0, + 2.0 * edge_quadrature_points[q_point][0] - 1.0, i); tmp(3) = weight * @@ -792,7 +792,7 @@ FE_Nedelec::initialize_restriction() { const double L_j = legendre_polynomials[j + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); for (unsigned int k = 0; k < tmp.size(); ++k) system_rhs(j, k) += tmp(k) * L_j; @@ -839,13 +839,13 @@ FE_Nedelec::initialize_restriction() { const double L_i = weight * legendre_polynomials[i].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); for (unsigned int j = 0; j < this->degree - 1; ++j) assembling_matrix(i * (this->degree - 1) + j, q_point) = L_i * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); } } @@ -870,13 +870,13 @@ FE_Nedelec::initialize_restriction() { tmp = 0.0; - if (quadrature_points[q_point](0) < 0.5) + if (quadrature_points[q_point][0] < 0.5) { - if (quadrature_points[q_point](1) < 0.5) + if (quadrature_points[q_point][1] < 0.5) { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0), - 2.0 * quadrature_points[q_point](1)); + 2.0 * quadrature_points[q_point][0], + 2.0 * quadrature_points[q_point][1]); tmp(0) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -887,8 +887,8 @@ FE_Nedelec::initialize_restriction() else { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0), - 2.0 * quadrature_points[q_point](1) - 1.0); + 2.0 * quadrature_points[q_point][0], + 2.0 * quadrature_points[q_point][1] - 1.0); tmp(4) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -897,11 +897,11 @@ FE_Nedelec::initialize_restriction() } } - else if (quadrature_points[q_point](1) < 0.5) + else if (quadrature_points[q_point][1] < 0.5) { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0) - 1.0, - 2.0 * quadrature_points[q_point](1)); + 2.0 * quadrature_points[q_point][0] - 1.0, + 2.0 * quadrature_points[q_point][1]); tmp(2) += 2.0 * this->shape_value_component(dof, @@ -916,8 +916,8 @@ FE_Nedelec::initialize_restriction() else { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0) - 1.0, - 2.0 * quadrature_points[q_point](1) - 1.0); + 2.0 * quadrature_points[q_point][0] - 1.0, + 2.0 * quadrature_points[q_point][1] - 1.0); tmp(6) += 2.0 * this->shape_value_component(dof, @@ -965,18 +965,18 @@ FE_Nedelec::initialize_restriction() for (unsigned int i = 0; i < this->degree; ++i) { const double L_i_0 = legendre_polynomials[i].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); const double L_i_1 = legendre_polynomials[i].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); for (unsigned int j = 0; j < this->degree - 1; ++j) { const double l_j_0 = L_i_0 * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); const double l_j_1 = L_i_1 * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); for (unsigned int k = 0; k < 4; ++k) { @@ -1029,12 +1029,12 @@ FE_Nedelec::initialize_restriction() { const double weight = 2.0 * edge_quadrature.weight(q_point); - if (edge_quadrature_points[q_point](0) < 0.5) + if (edge_quadrature_points[q_point][0] < 0.5) for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j < 2; ++j) { Point quadrature_point( - i, 2.0 * edge_quadrature_points[q_point](0), j); + i, 2.0 * edge_quadrature_points[q_point][0], j); this->restriction[index][i + 4 * j]((i + 4 * j) * this->degree, @@ -1042,7 +1042,7 @@ FE_Nedelec::initialize_restriction() weight * this->shape_value_component(dof, quadrature_point, 1); quadrature_point = - Point(2.0 * edge_quadrature_points[q_point](0), + Point(2.0 * edge_quadrature_points[q_point][0], i, j); this->restriction[index][2 * (i + 2 * j)]( @@ -1052,7 +1052,7 @@ FE_Nedelec::initialize_restriction() quadrature_point = Point(i, j, - 2.0 * edge_quadrature_points[q_point](0)); + 2.0 * edge_quadrature_points[q_point][0]); this->restriction[index][i + 2 * j]((i + 2 * (j + 4)) * this->degree, dof) += @@ -1065,7 +1065,7 @@ FE_Nedelec::initialize_restriction() for (unsigned int j = 0; j < 2; ++j) { Point quadrature_point( - i, 2.0 * edge_quadrature_points[q_point](0) - 1.0, j); + i, 2.0 * edge_quadrature_points[q_point][0] - 1.0, j); this->restriction[index][i + 4 * j + 2]((i + 4 * j) * this->degree, @@ -1073,13 +1073,13 @@ FE_Nedelec::initialize_restriction() weight * this->shape_value_component(dof, quadrature_point, 1); quadrature_point = Point( - 2.0 * edge_quadrature_points[q_point](0) - 1.0, i, j); + 2.0 * edge_quadrature_points[q_point][0] - 1.0, i, j); this->restriction[index][2 * (i + 2 * j) + 1]( (i + 4 * j + 2) * this->degree, dof) += weight * this->shape_value_component(dof, quadrature_point, 0); quadrature_point = Point( - i, j, 2.0 * edge_quadrature_points[q_point](0) - 1.0); + i, j, 2.0 * edge_quadrature_points[q_point][0] - 1.0); this->restriction[index][i + 2 * (j + 2)]( (i + 2 * (j + 4)) * this->degree, dof) += weight * @@ -1113,7 +1113,7 @@ FE_Nedelec::initialize_restriction() for (unsigned int i = 0; i < deg; ++i) assembling_matrix(i, q_point) = weight * legendre_polynomials[i + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); } FullMatrix system_matrix(deg, deg); @@ -1139,16 +1139,16 @@ FE_Nedelec::initialize_restriction() { const double weight = edge_quadrature.weight(q_point); const Point quadrature_point_0( - i, edge_quadrature_points[q_point](0), j); + i, edge_quadrature_points[q_point][0], j); const Point quadrature_point_1( - edge_quadrature_points[q_point](0), i, j); + edge_quadrature_points[q_point][0], i, j); const Point quadrature_point_2( - i, j, edge_quadrature_points[q_point](0)); + i, j, edge_quadrature_points[q_point][0]); - if (edge_quadrature_points[q_point](0) < 0.5) + if (edge_quadrature_points[q_point][0] < 0.5) { Point quadrature_point_3( - i, 2.0 * edge_quadrature_points[q_point](0), j); + i, 2.0 * edge_quadrature_points[q_point][0], j); tmp(0) = weight * (2.0 * this->shape_value_component( @@ -1168,7 +1168,7 @@ FE_Nedelec::initialize_restriction() quadrature_point_0, 1); quadrature_point_3 = Point( - 2.0 * edge_quadrature_points[q_point](0), i, j); + 2.0 * edge_quadrature_points[q_point][0], i, j); tmp(2) = weight * (2.0 * this->shape_value_component( @@ -1188,7 +1188,7 @@ FE_Nedelec::initialize_restriction() quadrature_point_1, 0); quadrature_point_3 = Point( - i, j, 2.0 * edge_quadrature_points[q_point](0)); + i, j, 2.0 * edge_quadrature_points[q_point][0]); tmp(4) = weight * (2.0 * this->shape_value_component( @@ -1222,7 +1222,7 @@ FE_Nedelec::initialize_restriction() Point quadrature_point_3( i, - 2.0 * edge_quadrature_points[q_point](0) - 1.0, + 2.0 * edge_quadrature_points[q_point][0] - 1.0, j); tmp(1) = weight * @@ -1243,7 +1243,7 @@ FE_Nedelec::initialize_restriction() quadrature_point_1, 0); quadrature_point_3 = Point( - 2.0 * edge_quadrature_points[q_point](0) - 1.0, + 2.0 * edge_quadrature_points[q_point][0] - 1.0, i, j); tmp(3) = @@ -1267,7 +1267,7 @@ FE_Nedelec::initialize_restriction() quadrature_point_3 = Point( i, j, - 2.0 * edge_quadrature_points[q_point](0) - 1.0); + 2.0 * edge_quadrature_points[q_point][0] - 1.0); tmp(5) = weight * (2.0 * this->shape_value_component( @@ -1284,7 +1284,7 @@ FE_Nedelec::initialize_restriction() { const double L_k = legendre_polynomials[k + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); for (unsigned int l = 0; l < tmp.size(); ++l) system_rhs(k, l) += tmp(l) * L_k; @@ -1339,12 +1339,12 @@ FE_Nedelec::initialize_restriction() { const double L_i = weight * legendre_polynomials[i].value( - face_quadrature_points[q_point](0)); + face_quadrature_points[q_point][0]); for (unsigned int j = 0; j < deg; ++j) assembling_matrix(i * deg + j, q_point) = L_i * lobatto_polynomials[j + 2].value( - face_quadrature_points[q_point](1)); + face_quadrature_points[q_point][1]); } } @@ -1371,30 +1371,30 @@ FE_Nedelec::initialize_restriction() { tmp = 0.0; - if (face_quadrature_points[q_point](0) < 0.5) + if (face_quadrature_points[q_point][0] < 0.5) { - if (face_quadrature_points[q_point](1) < 0.5) + if (face_quadrature_points[q_point][1] < 0.5) { Point quadrature_point_0( i, - 2.0 * face_quadrature_points[q_point](0), - 2.0 * face_quadrature_points[q_point](1)); + 2.0 * face_quadrature_points[q_point][0], + 2.0 * face_quadrature_points[q_point][1]); tmp(0) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 1); tmp(1) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0), + 2.0 * face_quadrature_points[q_point][0], i, - 2.0 * face_quadrature_points[q_point](1)); + 2.0 * face_quadrature_points[q_point][1]); tmp(8) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); tmp(9) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 0); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0), - 2.0 * face_quadrature_points[q_point](1), + 2.0 * face_quadrature_points[q_point][0], + 2.0 * face_quadrature_points[q_point][1], i); tmp(16) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 0); @@ -1406,8 +1406,8 @@ FE_Nedelec::initialize_restriction() { Point quadrature_point_0( i, - 2.0 * face_quadrature_points[q_point](0), - 2.0 * face_quadrature_points[q_point](1) - + 2.0 * face_quadrature_points[q_point][0], + 2.0 * face_quadrature_points[q_point][1] - 1.0); tmp(2) += 2.0 * this->shape_value_component( @@ -1415,17 +1415,17 @@ FE_Nedelec::initialize_restriction() tmp(3) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0), + 2.0 * face_quadrature_points[q_point][0], i, - 2.0 * face_quadrature_points[q_point](1) - + 2.0 * face_quadrature_points[q_point][1] - 1.0); tmp(10) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); tmp(11) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 0); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0), - 2.0 * face_quadrature_points[q_point](1) - + 2.0 * face_quadrature_points[q_point][0], + 2.0 * face_quadrature_points[q_point][1] - 1.0, i); tmp(18) += 2.0 * this->shape_value_component( @@ -1435,28 +1435,28 @@ FE_Nedelec::initialize_restriction() } } - else if (face_quadrature_points[q_point](1) < 0.5) + else if (face_quadrature_points[q_point][1] < 0.5) { Point quadrature_point_0( i, - 2.0 * face_quadrature_points[q_point](0) - 1.0, - 2.0 * face_quadrature_points[q_point](1)); + 2.0 * face_quadrature_points[q_point][0] - 1.0, + 2.0 * face_quadrature_points[q_point][1]); tmp(4) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 1); tmp(5) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0) - 1.0, + 2.0 * face_quadrature_points[q_point][0] - 1.0, i, - 2.0 * face_quadrature_points[q_point](1)); + 2.0 * face_quadrature_points[q_point][1]); tmp(12) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); tmp(13) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 0); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0) - 1.0, - 2.0 * face_quadrature_points[q_point](1), + 2.0 * face_quadrature_points[q_point][0] - 1.0, + 2.0 * face_quadrature_points[q_point][1], i); tmp(20) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 0); @@ -1468,24 +1468,24 @@ FE_Nedelec::initialize_restriction() { Point quadrature_point_0( i, - 2.0 * face_quadrature_points[q_point](0) - 1.0, - 2.0 * face_quadrature_points[q_point](1) - 1.0); + 2.0 * face_quadrature_points[q_point][0] - 1.0, + 2.0 * face_quadrature_points[q_point][1] - 1.0); tmp(6) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 1); tmp(7) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0) - 1.0, + 2.0 * face_quadrature_points[q_point][0] - 1.0, i, - 2.0 * face_quadrature_points[q_point](1) - 1.0); + 2.0 * face_quadrature_points[q_point][1] - 1.0); tmp(14) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 2); tmp(15) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 0); quadrature_point_0 = Point( - 2.0 * face_quadrature_points[q_point](0) - 1.0, - 2.0 * face_quadrature_points[q_point](1) - 1.0, + 2.0 * face_quadrature_points[q_point][0] - 1.0, + 2.0 * face_quadrature_points[q_point][1] - 1.0, i); tmp(22) += 2.0 * this->shape_value_component( dof, quadrature_point_0, 0); @@ -1495,15 +1495,15 @@ FE_Nedelec::initialize_restriction() const Point quadrature_point_0( i, - face_quadrature_points[q_point](0), - face_quadrature_points[q_point](1)); + face_quadrature_points[q_point][0], + face_quadrature_points[q_point][1]); const Point quadrature_point_1( - face_quadrature_points[q_point](0), + face_quadrature_points[q_point][0], i, - face_quadrature_points[q_point](1)); + face_quadrature_points[q_point][1]); const Point quadrature_point_2( - face_quadrature_points[q_point](0), - face_quadrature_points[q_point](1), + face_quadrature_points[q_point][0], + face_quadrature_points[q_point][1], i); for (unsigned int j = 0; j < 2; ++j) @@ -1559,18 +1559,18 @@ FE_Nedelec::initialize_restriction() for (unsigned int j = 0; j <= deg; ++j) { const double L_j_0 = legendre_polynomials[j].value( - face_quadrature_points[q_point](0)); + face_quadrature_points[q_point][0]); const double L_j_1 = legendre_polynomials[j].value( - face_quadrature_points[q_point](1)); + face_quadrature_points[q_point][1]); for (unsigned int k = 0; k < deg; ++k) { const double l_k_0 = L_j_0 * lobatto_polynomials[k + 2].value( - face_quadrature_points[q_point](1)); + face_quadrature_points[q_point][1]); const double l_k_1 = L_j_1 * lobatto_polynomials[k + 2].value( - face_quadrature_points[q_point](0)); + face_quadrature_points[q_point][0]); for (unsigned int l = 0; l < 4; ++l) { @@ -1673,19 +1673,19 @@ FE_Nedelec::initialize_restriction() { const double L_i = weight * legendre_polynomials[i].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); for (unsigned int j = 0; j < deg; ++j) { const double l_j = L_i * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); for (unsigned int k = 0; k < deg; ++k) assembling_matrix((i * deg + j) * deg + k, q_point) = l_j * lobatto_polynomials[k + 2].value( - quadrature_points[q_point](2)); + quadrature_points[q_point][2]); } } } @@ -1711,16 +1711,16 @@ FE_Nedelec::initialize_restriction() { tmp = 0.0; - if (quadrature_points[q_point](0) < 0.5) + if (quadrature_points[q_point][0] < 0.5) { - if (quadrature_points[q_point](1) < 0.5) + if (quadrature_points[q_point][1] < 0.5) { - if (quadrature_points[q_point](2) < 0.5) + if (quadrature_points[q_point][2] < 0.5) { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0), - 2.0 * quadrature_points[q_point](1), - 2.0 * quadrature_points[q_point](2)); + 2.0 * quadrature_points[q_point][0], + 2.0 * quadrature_points[q_point][1], + 2.0 * quadrature_points[q_point][2]); tmp(0) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -1733,9 +1733,9 @@ FE_Nedelec::initialize_restriction() else { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0), - 2.0 * quadrature_points[q_point](1), - 2.0 * quadrature_points[q_point](2) - 1.0); + 2.0 * quadrature_points[q_point][0], + 2.0 * quadrature_points[q_point][1], + 2.0 * quadrature_points[q_point][2] - 1.0); tmp(3) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -1746,12 +1746,12 @@ FE_Nedelec::initialize_restriction() } } - else if (quadrature_points[q_point](2) < 0.5) + else if (quadrature_points[q_point][2] < 0.5) { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0), - 2.0 * quadrature_points[q_point](1) - 1.0, - 2.0 * quadrature_points[q_point](2)); + 2.0 * quadrature_points[q_point][0], + 2.0 * quadrature_points[q_point][1] - 1.0, + 2.0 * quadrature_points[q_point][2]); tmp(6) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -1764,9 +1764,9 @@ FE_Nedelec::initialize_restriction() else { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0), - 2.0 * quadrature_points[q_point](1) - 1.0, - 2.0 * quadrature_points[q_point](2) - 1.0); + 2.0 * quadrature_points[q_point][0], + 2.0 * quadrature_points[q_point][1] - 1.0, + 2.0 * quadrature_points[q_point][2] - 1.0); tmp(9) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -1777,14 +1777,14 @@ FE_Nedelec::initialize_restriction() } } - else if (quadrature_points[q_point](1) < 0.5) + else if (quadrature_points[q_point][1] < 0.5) { - if (quadrature_points[q_point](2) < 0.5) + if (quadrature_points[q_point][2] < 0.5) { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0) - 1.0, - 2.0 * quadrature_points[q_point](1), - 2.0 * quadrature_points[q_point](2)); + 2.0 * quadrature_points[q_point][0] - 1.0, + 2.0 * quadrature_points[q_point][1], + 2.0 * quadrature_points[q_point][2]); tmp(12) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -1797,9 +1797,9 @@ FE_Nedelec::initialize_restriction() else { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0) - 1.0, - 2.0 * quadrature_points[q_point](1), - 2.0 * quadrature_points[q_point](2) - 1.0); + 2.0 * quadrature_points[q_point][0] - 1.0, + 2.0 * quadrature_points[q_point][1], + 2.0 * quadrature_points[q_point][2] - 1.0); tmp(15) += 2.0 * this->shape_value_component( dof, quadrature_point, 0); @@ -1810,12 +1810,12 @@ FE_Nedelec::initialize_restriction() } } - else if (quadrature_points[q_point](2) < 0.5) + else if (quadrature_points[q_point][2] < 0.5) { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0) - 1.0, - 2.0 * quadrature_points[q_point](1) - 1.0, - 2.0 * quadrature_points[q_point](2)); + 2.0 * quadrature_points[q_point][0] - 1.0, + 2.0 * quadrature_points[q_point][1] - 1.0, + 2.0 * quadrature_points[q_point][2]); tmp(18) += 2.0 * this->shape_value_component(dof, @@ -1834,9 +1834,9 @@ FE_Nedelec::initialize_restriction() else { const Point quadrature_point( - 2.0 * quadrature_points[q_point](0) - 1.0, - 2.0 * quadrature_points[q_point](1) - 1.0, - 2.0 * quadrature_points[q_point](2) - 1.0); + 2.0 * quadrature_points[q_point][0] - 1.0, + 2.0 * quadrature_points[q_point][1] - 1.0, + 2.0 * quadrature_points[q_point][2] - 1.0); tmp(21) += 2.0 * this->shape_value_component(dof, @@ -1955,35 +1955,35 @@ FE_Nedelec::initialize_restriction() for (unsigned int i = 0; i <= deg; ++i) { const double L_i_0 = legendre_polynomials[i].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); const double L_i_1 = legendre_polynomials[i].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); const double L_i_2 = legendre_polynomials[i].value( - quadrature_points[q_point](2)); + quadrature_points[q_point][2]); for (unsigned int j = 0; j < deg; ++j) { const double l_j_0 = L_i_0 * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); const double l_j_1 = L_i_1 * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); const double l_j_2 = L_i_2 * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); for (unsigned int k = 0; k < deg; ++k) { const double l_k_0 = l_j_0 * lobatto_polynomials[k + 2].value( - quadrature_points[q_point](2)); + quadrature_points[q_point][2]); const double l_k_1 = l_j_1 * lobatto_polynomials[k + 2].value( - quadrature_points[q_point](2)); + quadrature_points[q_point][2]); const double l_k_2 = l_j_2 * lobatto_polynomials[k + 2].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); for (unsigned int l = 0; l < 8; ++l) { @@ -2643,7 +2643,7 @@ FE_Nedelec::get_subface_interpolation_matrix( ++q_point) { const Point quadrature_point( - 0.0, 0.5 * (edge_quadrature_points[q_point](0) + subface)); + 0.0, 0.5 * (edge_quadrature_points[q_point][0] + subface)); interpolation_matrix(0, dof) += 0.5 * edge_quadrature.weight(q_point) * @@ -2673,7 +2673,7 @@ FE_Nedelec::get_subface_interpolation_matrix( for (unsigned int i = 0; i < source_fe.degree - 1; ++i) assembling_matrix(i, q_point) = weight * legendre_polynomials[i + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); } FullMatrix system_matrix(source_fe.degree - 1, @@ -2697,9 +2697,9 @@ FE_Nedelec::get_subface_interpolation_matrix( { const Point quadrature_point_0( 0.0, - 0.5 * (edge_quadrature_points[q_point](0) + subface)); + 0.5 * (edge_quadrature_points[q_point][0] + subface)); const Point quadrature_point_1( - 0.0, edge_quadrature_points[q_point](0)); + 0.0, edge_quadrature_points[q_point][0]); const double tmp = edge_quadrature.weight(q_point) * (0.5 * this->shape_value_component(dof, @@ -2713,7 +2713,7 @@ FE_Nedelec::get_subface_interpolation_matrix( for (unsigned int i = 0; i < source_fe.degree - 1; ++i) system_rhs(i) += tmp * legendre_polynomials[i + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); } system_matrix_inv.vmult(solution, system_rhs); @@ -2745,7 +2745,7 @@ FE_Nedelec::get_subface_interpolation_matrix( { Point quadrature_point( 0.5 * (i + shifts[subface][0]), - 0.5 * (edge_quadrature_points[q_point](0) + + 0.5 * (edge_quadrature_points[q_point][0] + shifts[subface][1]), 0.0); @@ -2754,7 +2754,7 @@ FE_Nedelec::get_subface_interpolation_matrix( this->shape_value_component( this->face_to_cell_index(dof, 4), quadrature_point, 1); quadrature_point = - Point(0.5 * (edge_quadrature_points[q_point](0) + + Point(0.5 * (edge_quadrature_points[q_point][0] + shifts[subface][0]), 0.5 * (i + shifts[subface][1]), 0.0); @@ -2788,7 +2788,7 @@ FE_Nedelec::get_subface_interpolation_matrix( for (unsigned int i = 0; i < source_fe.degree - 1; ++i) assembling_matrix(i, q_point) = weight * legendre_polynomials[i + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); } FullMatrix system_matrix(source_fe.degree - 1, @@ -2819,11 +2819,11 @@ FE_Nedelec::get_subface_interpolation_matrix( { Point quadrature_point_0( 0.5 * (i + shifts[subface][0]), - 0.5 * (edge_quadrature_points[q_point](0) + + 0.5 * (edge_quadrature_points[q_point][0] + shifts[subface][1]), 0.0); Point quadrature_point_1( - i, edge_quadrature_points[q_point](0), 0.0); + i, edge_quadrature_points[q_point][0], 0.0); tmp(i) = weight * @@ -2836,12 +2836,12 @@ FE_Nedelec::get_subface_interpolation_matrix( i * source_fe.degree, quadrature_point_1, 1)); quadrature_point_0 = Point(0.5 * - (edge_quadrature_points[q_point](0) + + (edge_quadrature_points[q_point][0] + shifts[subface][0]), 0.5 * (i + shifts[subface][1]), 0.0); quadrature_point_1 = - Point(edge_quadrature_points[q_point](0), + Point(edge_quadrature_points[q_point][0], i, 0.0); tmp(i + 2) = @@ -2861,7 +2861,7 @@ FE_Nedelec::get_subface_interpolation_matrix( for (unsigned int i = 0; i < source_fe.degree - 1; ++i) { const double L_i = legendre_polynomials[i + 1].value( - edge_quadrature_points[q_point](0)); + edge_quadrature_points[q_point][0]); for (unsigned int j = 0; j < GeometryInfo::lines_per_face; @@ -2906,13 +2906,13 @@ FE_Nedelec::get_subface_interpolation_matrix( { const double L_i = weight * legendre_polynomials[i].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); for (unsigned int j = 0; j < source_fe.degree - 1; ++j) assembling_matrix(i * (source_fe.degree - 1) + j, q_point) = L_i * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); } } @@ -2938,9 +2938,9 @@ FE_Nedelec::get_subface_interpolation_matrix( { Point quadrature_point( 0.5 * - (quadrature_points[q_point](0) + shifts[subface][0]), + (quadrature_points[q_point][0] + shifts[subface][0]), 0.5 * - (quadrature_points[q_point](1) + shifts[subface][1]), + (quadrature_points[q_point][1] + shifts[subface][1]), 0.0); tmp(0) = 0.5 * this->shape_value_component( this->face_to_cell_index(dof, 4), @@ -2951,8 +2951,8 @@ FE_Nedelec::get_subface_interpolation_matrix( quadrature_point, 1); quadrature_point = - Point(quadrature_points[q_point](0), - quadrature_points[q_point](1), + Point(quadrature_points[q_point][0], + quadrature_points[q_point][1], 0.0); for (unsigned int i = 0; i < 2; ++i) @@ -2976,9 +2976,9 @@ FE_Nedelec::get_subface_interpolation_matrix( for (unsigned int i = 0; i < source_fe.degree; ++i) { const double L_i_0 = legendre_polynomials[i].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); const double L_i_1 = legendre_polynomials[i].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); for (unsigned int j = 0; j < source_fe.degree - 1; ++j) @@ -2986,11 +2986,11 @@ FE_Nedelec::get_subface_interpolation_matrix( system_rhs(i * (source_fe.degree - 1) + j, 0) += tmp(0) * L_i_0 * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](1)); + quadrature_points[q_point][1]); system_rhs(i * (source_fe.degree - 1) + j, 1) += tmp(1) * L_i_1 * lobatto_polynomials[j + 2].value( - quadrature_points[q_point](0)); + quadrature_points[q_point][0]); } } } @@ -3228,8 +3228,8 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( system_matrix(i, j) += boundary_weights(q_point, j) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[face_no][q_point]( - 0)); + this->generalized_face_support_points[face_no][q_point] + [0]); FullMatrix system_matrix_inv(this->degree - 1, this->degree - 1); @@ -3311,19 +3311,19 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( legendre_polynomials[i].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](0)) * + n_edge_points][0]) * lobatto_polynomials[j + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](1)) * + n_edge_points][1]) * lobatto_polynomials_grad[k].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](0)) * + n_edge_points][0]) * lobatto_polynomials[l + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](1)); + n_edge_points][1]); system_matrix_inv.reinit(system_matrix.m(), system_matrix.m()); system_matrix_inv.invert(system_matrix); @@ -3358,11 +3358,11 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( lobatto_polynomials_grad[i].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](0)) * + n_edge_points][0]) * lobatto_polynomials[j + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](1)); + n_edge_points][1]); } solution.reinit(system_matrix.m()); @@ -3409,11 +3409,11 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( lobatto_polynomials_grad[i].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](1)) * + n_edge_points][1]) * lobatto_polynomials[j + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * - n_edge_points](0)); + n_edge_points][0]); } system_matrix_inv.vmult(solution, system_rhs); @@ -3508,8 +3508,8 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( system_matrix(i, j) += boundary_weights(q_point, j) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[face_no][q_point]( - 1)); + this->generalized_face_support_points[face_no][q_point] + [1]); FullMatrix system_matrix_inv(this->degree - 1, this->degree - 1); @@ -3588,10 +3588,10 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( 2 * (k * (this->degree - 1) + l)) * legendre_polynomials[i].value( this->generalized_face_support_points - [face_no][q_point + 4 * n_edge_points](0)) * + [face_no][q_point + 4 * n_edge_points][0]) * lobatto_polynomials[j + 2].value( this->generalized_face_support_points - [face_no][q_point + 4 * n_edge_points](1)); + [face_no][q_point + 4 * n_edge_points][1]); system_matrix_inv.reinit(system_matrix.m(), system_matrix.m()); system_matrix_inv.invert(system_matrix); @@ -3745,42 +3745,42 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](0)) * + n_face_points][0]) * lobatto_polynomials[j + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](1)) * + n_face_points][1]) * lobatto_polynomials[k + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](2)) * + n_face_points][2]) * lobatto_polynomials_grad[l].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](0)) * + n_face_points][0]) * lobatto_polynomials[m + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](1)) * + n_face_points][1]) * lobatto_polynomials[n + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](2)); + n_face_points][2]); system_matrix_inv.reinit(system_matrix.m(), system_matrix.m()); system_matrix_inv.invert(system_matrix); @@ -3847,21 +3847,21 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](0)) * + n_face_points][0]) * lobatto_polynomials[j + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](1)) * + n_face_points][1]) * lobatto_polynomials[k + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](2)); + n_face_points][2]); } solution.reinit(system_rhs.size()); @@ -3958,21 +3958,21 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](1)) * + n_face_points][1]) * lobatto_polynomials[j + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](0)) * + n_face_points][0]) * lobatto_polynomials[k + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](2)); + n_face_points][2]); } system_matrix_inv.vmult(solution, system_rhs); @@ -4051,21 +4051,21 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](2)) * + n_face_points][2]) * lobatto_polynomials[j + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](0)) * + n_face_points][0]) * lobatto_polynomials[k + 2].value( this->generalized_support_points [q_point + GeometryInfo::lines_per_cell * n_edge_points + GeometryInfo::faces_per_cell * - n_face_points](1)); + n_face_points][1]); } system_matrix_inv.vmult(solution, system_rhs); diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc index 07a85344a0..680279c38f 100644 --- a/source/fe/fe_p1nc.cc +++ b/source/fe/fe_p1nc.cc @@ -99,33 +99,33 @@ FE_P1NC::get_linear_shape_coefficients( const Point<2> cpt = (cell->vertex(0) + cell->vertex(1) + cell->vertex(2) + cell->vertex(3)) / 4; - const double det = (mpt[0](0) - mpt[1](0)) * (mpt[2](1) - mpt[3](1)) - - (mpt[2](0) - mpt[3](0)) * (mpt[0](1) - mpt[1](1)); + const double det = (mpt[0][0] - mpt[1][0]) * (mpt[2][1] - mpt[3][1]) - + (mpt[2][0] - mpt[3][0]) * (mpt[0][1] - mpt[1][1]); ndarray coeffs; coeffs[0][0] = - ((mpt[2](1) - mpt[3](1)) * (0.5) - (mpt[0](1) - mpt[1](1)) * (0.5)) / det; + ((mpt[2][1] - mpt[3][1]) * (0.5) - (mpt[0][1] - mpt[1][1]) * (0.5)) / det; coeffs[1][0] = - ((mpt[2](1) - mpt[3](1)) * (-0.5) - (mpt[0](1) - mpt[1](1)) * (0.5)) / det; + ((mpt[2][1] - mpt[3][1]) * (-0.5) - (mpt[0][1] - mpt[1][1]) * (0.5)) / det; coeffs[2][0] = - ((mpt[2](1) - mpt[3](1)) * (0.5) - (mpt[0](1) - mpt[1](1)) * (-0.5)) / det; + ((mpt[2][1] - mpt[3][1]) * (0.5) - (mpt[0][1] - mpt[1][1]) * (-0.5)) / det; coeffs[3][0] = - ((mpt[2](1) - mpt[3](1)) * (-0.5) - (mpt[0](1) - mpt[1](1)) * (-0.5)) / det; + ((mpt[2][1] - mpt[3][1]) * (-0.5) - (mpt[0][1] - mpt[1][1]) * (-0.5)) / det; coeffs[0][1] = - (-(mpt[2](0) - mpt[3](0)) * (0.5) + (mpt[0](0) - mpt[1](0)) * (0.5)) / det; + (-(mpt[2][0] - mpt[3][0]) * (0.5) + (mpt[0][0] - mpt[1][0]) * (0.5)) / det; coeffs[1][1] = - (-(mpt[2](0) - mpt[3](0)) * (-0.5) + (mpt[0](0) - mpt[1](0)) * (0.5)) / det; + (-(mpt[2][0] - mpt[3][0]) * (-0.5) + (mpt[0][0] - mpt[1][0]) * (0.5)) / det; coeffs[2][1] = - (-(mpt[2](0) - mpt[3](0)) * (0.5) + (mpt[0](0) - mpt[1](0)) * (-0.5)) / det; + (-(mpt[2][0] - mpt[3][0]) * (0.5) + (mpt[0][0] - mpt[1][0]) * (-0.5)) / det; coeffs[3][1] = - (-(mpt[2](0) - mpt[3](0)) * (-0.5) + (mpt[0](0) - mpt[1](0)) * (-0.5)) / + (-(mpt[2][0] - mpt[3][0]) * (-0.5) + (mpt[0][0] - mpt[1][0]) * (-0.5)) / det; - coeffs[0][2] = 0.25 - cpt(0) * coeffs[0][0] - cpt(1) * coeffs[0][1]; - coeffs[1][2] = 0.25 - cpt(0) * coeffs[1][0] - cpt(1) * coeffs[1][1]; - coeffs[2][2] = 0.25 - cpt(0) * coeffs[2][0] - cpt(1) * coeffs[2][1]; - coeffs[3][2] = 0.25 - cpt(0) * coeffs[3][0] - cpt(1) * coeffs[3][1]; + coeffs[0][2] = 0.25 - cpt[0] * coeffs[0][0] - cpt[1] * coeffs[0][1]; + coeffs[1][2] = 0.25 - cpt[0] * coeffs[1][0] - cpt[1] * coeffs[1][1]; + coeffs[2][2] = 0.25 - cpt[0] * coeffs[2][0] - cpt[1] * coeffs[2][1]; + coeffs[3][2] = 0.25 - cpt[0] * coeffs[3][0] - cpt[1] * coeffs[3][1]; return coeffs; } @@ -231,8 +231,8 @@ FE_P1NC::fill_fe_values( for (unsigned int i = 0; i < n_q_points; ++i) for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) output_data.shape_values[k][i] = - (coeffs[k][0] * mapping_data.quadrature_points[i](0) + - coeffs[k][1] * mapping_data.quadrature_points[i](1) + coeffs[k][2]); + (coeffs[k][0] * mapping_data.quadrature_points[i][0] + + coeffs[k][1] * mapping_data.quadrature_points[i][1] + coeffs[k][2]); if (flags & update_gradients) for (unsigned int i = 0; i < n_q_points; ++i) @@ -277,8 +277,8 @@ FE_P1NC::fill_fe_face_values( quadrature_on_face.point(i)); output_data.shape_values[k][i] = - (coeffs[k][0] * quadrature_point(0) + - coeffs[k][1] * quadrature_point(1) + coeffs[k][2]); + (coeffs[k][0] * quadrature_point[0] + + coeffs[k][1] * quadrature_point[1] + coeffs[k][2]); } if (flags & update_gradients) @@ -322,8 +322,8 @@ FE_P1NC::fill_fe_subface_values( cell, quadrature_on_subface.point(i)); output_data.shape_values[k][i] = - (coeffs[k][0] * quadrature_point(0) + - coeffs[k][1] * quadrature_point(1) + coeffs[k][2]); + (coeffs[k][0] * quadrature_point[0] + + coeffs[k][1] * quadrature_point[1] + coeffs[k][2]); } } diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index 7006d4cafc..1b2b49a3a6 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -188,7 +188,7 @@ FE_PolyTensor::FE_PolyTensor( , mapping_kind({MappingKind::mapping_none}) , poly_space(polynomials.clone()) { - cached_point(0) = -1; + cached_point[0] = -1; // Set up the table converting // components to base // components. Since we have only diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc index 198e85c6d1..33c6c9e5b2 100644 --- a/source/fe/fe_q.cc +++ b/source/fe/fe_q.cc @@ -134,7 +134,7 @@ FE_Q::get_name() const const QGaussLobatto<1> points_gl(this->degree + 1); bool gauss_lobatto = true; for (unsigned int j = 0; j <= this->degree; ++j) - if (points[j] != points_gl.point(j)(0)) + if (points[j] != points_gl.point(j)[0]) { gauss_lobatto = false; break; diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index afabd9881e..20dc0c5a2c 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -199,7 +199,7 @@ struct FE_Q_Base::Implementation for (unsigned int j = 0; j < q_deg + 1; ++j) { Point p; - p[0] = constraint_points[i](0); + p[0] = constraint_points[i][0]; fe.interface_constraints(i, face_index_map[j]) = fe.poly_space->compute_value(index_map_inverse[j], p); @@ -257,7 +257,7 @@ struct FE_Q_Base::Implementation const double step = 1. / q_deg; std::vector> line_support_points(n); for (unsigned int i = 0; i < n; ++i) - line_support_points[i](0) = (i + 1) * step; + line_support_points[i][0] = (i + 1) * step; Quadrature qline(line_support_points); // auxiliary points in 2d @@ -357,8 +357,8 @@ struct FE_Q_Base::Implementation for (unsigned int k = 0; k < dim - 1; ++k) { const int coord_int = - static_cast(constraint_points[i](k) * interval + 0.25); - constraint_point(k) = 1. * coord_int / interval; + static_cast(constraint_points[i][k] * interval + 0.25); + constraint_point[k] = 1. * coord_int / interval; // The following lines of code should eliminate the problems with // the constraints object which appeared for P>=4. The @@ -380,9 +380,9 @@ struct FE_Q_Base::Implementation // // For a different explanation of the problem, see the discussion // in the FiniteElement class for constraint matrices in 3d. - mirror[k] = (constraint_point(k) > 0.5); + mirror[k] = (constraint_point[k] > 0.5); if (mirror[k]) - constraint_point(k) = 1.0 - constraint_point(k); + constraint_point[k] = 1.0 - constraint_point[k]; } for (unsigned int j = 0; j < pnts; ++j) diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc index bb5e976090..3be49d1315 100644 --- a/source/fe/fe_q_bubbles.cc +++ b/source/fe/fe_q_bubbles.cc @@ -155,7 +155,7 @@ namespace internal fine.JxW(q); Point quad_tmp; for (unsigned int k = 0; k < dim; ++k) - quad_tmp(k) = fine.quadrature_point(q)(k); + quad_tmp[k] = fine.quadrature_point(q)[k]; coarse_rhs_matrix(gdi, j) += fine.shape_value(i, q) * fe.shape_value(j, quad_tmp) * fine.JxW(q); @@ -282,16 +282,16 @@ FE_Q_Bubbles::get_name() const // Decode the support points in one coordinate direction. for (unsigned int j = 0; j < dofs_per_cell; ++j) { - if ((dim > 1) ? (unit_support_points[j](1) == 0 && - ((dim > 2) ? unit_support_points[j](2) == 0 : true)) : + if ((dim > 1) ? (unit_support_points[j][1] == 0 && + ((dim > 2) ? unit_support_points[j][2] == 0 : true)) : true) { if (index == 0) - points[index] = unit_support_points[j](0); + points[index] = unit_support_points[j][0]; else if (index == 1) - points[n_points - 1] = unit_support_points[j](0); + points[n_points - 1] = unit_support_points[j][0]; else - points[index - 1] = unit_support_points[j](0); + points[index - 1] = unit_support_points[j][0]; ++index; } @@ -325,7 +325,7 @@ FE_Q_Bubbles::get_name() const const QGaussLobatto<1> points_gl(n_points); type = true; for (unsigned int j = 0; j < n_points; ++j) - if (points[j] != points_gl.point(j)(0)) + if (points[j] != points_gl.point(j)[0]) { type = false; break; diff --git a/source/fe/fe_q_dg0.cc b/source/fe/fe_q_dg0.cc index e5e898233b..355ebbfab0 100644 --- a/source/fe/fe_q_dg0.cc +++ b/source/fe/fe_q_dg0.cc @@ -108,16 +108,16 @@ FE_Q_DG0::get_name() const // Decode the support points in one coordinate direction. for (unsigned int j = 0; j < dofs_per_cell; ++j) { - if ((dim > 1) ? (unit_support_points[j](1) == 0 && - ((dim > 2) ? unit_support_points[j](2) == 0 : true)) : + if ((dim > 1) ? (unit_support_points[j][1] == 0 && + ((dim > 2) ? unit_support_points[j][2] == 0 : true)) : true) { if (index == 0) - points[index] = unit_support_points[j](0); + points[index] = unit_support_points[j][0]; else if (index == 1) - points[n_points - 1] = unit_support_points[j](0); + points[n_points - 1] = unit_support_points[j][0]; else - points[index - 1] = unit_support_points[j](0); + points[index - 1] = unit_support_points[j][0]; ++index; } @@ -150,7 +150,7 @@ FE_Q_DG0::get_name() const const QGaussLobatto<1> points_gl(n_points); type = true; for (unsigned int j = 0; j < n_points; ++j) - if (points[j] != points_gl.point(j)(0)) + if (points[j] != points_gl.point(j)[0]) { type = false; break; diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index 4c5e51ada7..fb44b9e2bd 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -862,28 +862,28 @@ FE_Q_Hierarchical::initialize_generalized_support_points() for (unsigned int ix = 0; ix <= this->degree; ++ix) { if (ix == 0) - p(0) = 0.; + p[0] = 0.; else if (ix == 1) - p(0) = 1.; + p[0] = 1.; else - p(0) = .5; + p[0] = .5; if (dim > 1) { if (iy == 0) - p(1) = 0.; + p[1] = 0.; else if (iy == 1) - p(1) = 1.; + p[1] = 1.; else - p(1) = .5; + p[1] = .5; } if (dim > 2) { if (iz == 0) - p(2) = 0.; + p[2] = 0.; else if (iz == 1) - p(2) = 1.; + p[2] = 1.; else - p(2) = .5; + p[2] = .5; } this->generalized_support_points[index_map_inverse[k++]] = p; } @@ -1933,28 +1933,28 @@ FE_Q_Hierarchical::initialize_generalized_face_support_points() for (unsigned int ix = 0; ix <= this->degree; ++ix) { if (ix == 0) - p(0) = 0.; + p[0] = 0.; else if (ix == 1) - p(0) = 1.; + p[0] = 1.; else - p(0) = .5; + p[0] = .5; if (codim > 1) { if (iy == 0) - p(1) = 0.; + p[1] = 0.; else if (iy == 1) - p(1) = 1.; + p[1] = 1.; else - p(1) = .5; + p[1] = .5; } if (codim > 2) { if (iz == 0) - p(2) = 0.; + p[2] = 0.; else if (iz == 1) - p(2) = 1.; + p[2] = 1.; else - p(2) = .5; + p[2] = .5; } this->generalized_face_support_points[face_no][face_renumber[k++]] = p; diff --git a/source/fe/fe_rannacher_turek.cc b/source/fe/fe_rannacher_turek.cc index 873b44c2c3..d6d3bf51be 100644 --- a/source/fe/fe_rannacher_turek.cc +++ b/source/fe/fe_rannacher_turek.cc @@ -96,13 +96,13 @@ FE_RannacherTurek::initialize_support_points() for (unsigned int q = 0; q < face_quadrature.size(); ++q) { this->generalized_support_points[0 * face_quadrature.size() + q] = - dealii::Point(0, 1 - face_quadrature.point(q)(0)); + dealii::Point(0, 1 - face_quadrature.point(q)[0]); this->generalized_support_points[1 * face_quadrature.size() + q] = - dealii::Point(1, 1 - face_quadrature.point(q)(0)); + dealii::Point(1, 1 - face_quadrature.point(q)[0]); this->generalized_support_points[2 * face_quadrature.size() + q] = - dealii::Point(face_quadrature.point(q)(0), 0); + dealii::Point(face_quadrature.point(q)[0], 0); this->generalized_support_points[3 * face_quadrature.size() + q] = - dealii::Point(face_quadrature.point(q)(0), 1); + dealii::Point(face_quadrature.point(q)[0], 1); } } diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index c4f9c5d0cc..e7bc4f7353 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -171,18 +171,18 @@ Mapping::project_real_point_to_unit_point_on_face( if (dim == 2) { if (unit_normal_direction == 0) - return Point{unit_cell_pt(1)}; + return Point{unit_cell_pt[1]}; else if (unit_normal_direction == 1) - return Point{unit_cell_pt(0)}; + return Point{unit_cell_pt[0]}; } else if (dim == 3) { if (unit_normal_direction == 0) - return Point{unit_cell_pt(1), unit_cell_pt(2)}; + return Point{unit_cell_pt[1], unit_cell_pt[2]}; else if (unit_normal_direction == 1) - return Point{unit_cell_pt(0), unit_cell_pt(2)}; + return Point{unit_cell_pt[0], unit_cell_pt[2]}; else if (unit_normal_direction == 2) - return Point{unit_cell_pt(0), unit_cell_pt(1)}; + return Point{unit_cell_pt[0], unit_cell_pt[1]}; } // We should never get here diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 2fbd42983b..62f3064c54 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -243,16 +243,16 @@ MappingCartesian::update_cell_extents( switch (dim) { case 1: - data.cell_extents[0] = cell->vertex(1)(0) - start(0); + data.cell_extents[0] = cell->vertex(1)[0] - start[0]; break; case 2: - data.cell_extents[0] = cell->vertex(1)(0) - start(0); - data.cell_extents[1] = cell->vertex(2)(1) - start(1); + data.cell_extents[0] = cell->vertex(1)[0] - start[0]; + data.cell_extents[1] = cell->vertex(2)[1] - start[1]; break; case 3: - data.cell_extents[0] = cell->vertex(1)(0) - start(0); - data.cell_extents[1] = cell->vertex(2)(1) - start(1); - data.cell_extents[2] = cell->vertex(4)(2) - start(2); + data.cell_extents[0] = cell->vertex(1)[0] - start[0]; + data.cell_extents[1] = cell->vertex(2)[1] - start[1]; + data.cell_extents[2] = cell->vertex(4)[2] - start[2]; break; default: Assert(false, ExcNotImplemented()); @@ -356,8 +356,8 @@ MappingCartesian::transform_quadrature_points( { quadrature_points[i] = start; for (unsigned int d = 0; d < dim; ++d) - quadrature_points[i](d) += - data.cell_extents[d] * data.quadrature_points[i + offset](d); + quadrature_points[i][d] += + data.cell_extents[d] * data.quadrature_points[i + offset][d]; } } @@ -1177,16 +1177,16 @@ MappingCartesian::transform_unit_to_real_cell( switch (dim) { case 1: - length[0] = cell->vertex(1)(0) - start(0); + length[0] = cell->vertex(1)[0] - start[0]; break; case 2: - length[0] = cell->vertex(1)(0) - start(0); - length[1] = cell->vertex(2)(1) - start(1); + length[0] = cell->vertex(1)[0] - start[0]; + length[1] = cell->vertex(2)[1] - start[1]; break; case 3: - length[0] = cell->vertex(1)(0) - start(0); - length[1] = cell->vertex(2)(1) - start(1); - length[2] = cell->vertex(4)(2) - start(2); + length[0] = cell->vertex(1)[0] - start[0]; + length[1] = cell->vertex(2)[1] - start[1]; + length[2] = cell->vertex(4)[2] - start[2]; break; default: Assert(false, ExcNotImplemented()); @@ -1194,7 +1194,7 @@ MappingCartesian::transform_unit_to_real_cell( Point p_real = cell->vertex(0); for (unsigned int d = 0; d < dim; ++d) - p_real(d) += length[d] * p(d); + p_real[d] += length[d] * p[d]; return p_real; } @@ -1218,16 +1218,16 @@ MappingCartesian::transform_real_to_unit_cell( switch (dim) { case 1: - real(0) /= cell->vertex(1)(0) - start(0); + real[0] /= cell->vertex(1)[0] - start[0]; break; case 2: - real(0) /= cell->vertex(1)(0) - start(0); - real(1) /= cell->vertex(2)(1) - start(1); + real[0] /= cell->vertex(1)[0] - start[0]; + real[1] /= cell->vertex(2)[1] - start[1]; break; case 3: - real(0) /= cell->vertex(1)(0) - start(0); - real(1) /= cell->vertex(2)(1) - start(1); - real(2) /= cell->vertex(4)(2) - start(2); + real[0] /= cell->vertex(1)[0] - start[0]; + real[1] /= cell->vertex(2)[1] - start[1]; + real[2] /= cell->vertex(4)[2] - start[2]; break; default: Assert(false, ExcNotImplemented()); diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 73444b89f9..943098bd80 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -531,8 +531,8 @@ MappingQ::transform_real_to_unit_cell( // the cell. only take the computed point if it lies // inside the reference cell const double eps = 1e-15; - if (-eps <= point(1) && point(1) <= 1 + eps && - -eps <= point(0) && point(0) <= 1 + eps) + if (-eps <= point[1] && point[1] <= 1 + eps && + -eps <= point[0] && point[0] <= 1 + eps) { return point; } diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc index c7ff940f1a..eea4b4d0af 100644 --- a/source/fe/mapping_q_eulerian.cc +++ b/source/fe/mapping_q_eulerian.cc @@ -182,7 +182,7 @@ MappingQEulerian::compute_mapping_support_points( { a[q] = fe_values.quadrature_point(q); for (unsigned int d = 0; d < spacedim; ++d) - a[q](d) += shift_vector[q](d); + a[q][d] += shift_vector[q][d]; } return a; diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index e6e3db8b81..2ece4e68e4 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -205,11 +205,11 @@ namespace GridGenerator Point<2>( 1)}}} /* dummy vector since we are asserting later*/), data.airfoil_length)) - , end_b0_x_u(airfoil_1D[0][n_cells_x_0](0)) - , end_b0_x_l(airfoil_1D[1][n_cells_x_0](0)) - , nose_x(airfoil_1D[0].front()(0)) - , tail_x(airfoil_1D[0].back()(0)) - , tail_y(airfoil_1D[0].back()(1)) + , end_b0_x_u(airfoil_1D[0][n_cells_x_0][0]) + , end_b0_x_l(airfoil_1D[1][n_cells_x_0][0]) + , nose_x(airfoil_1D[0].front()[0]) + , tail_x(airfoil_1D[0].back()[0]) + , tail_y(airfoil_1D[0].back()[1]) , center_mesh(0.5 * std::abs(end_b0_x_u + end_b0_x_l)) , length_b1_x(tail_x - center_mesh) , gamma(std::atan(height / @@ -227,8 +227,8 @@ namespace GridGenerator , H(tail_x, 0) , I(tail_x, -height) , J(tail_x + length_b2, 0) - , K(J(0), G(1)) - , L(J(0), I(1)) + , K(J[0], G[1]) + , L(J[0], I[1]) { Assert(data.airfoil_type == "Joukowski" || data.airfoil_type == "NACA", @@ -444,14 +444,14 @@ namespace GridGenerator // find index in vector to nose point (min) and tail point (max) for (unsigned int i = 0; i < jouk_points.size(); ++i) { - if (jouk_points[i](0) < nose_x_coordinate) + if (jouk_points[i][0] < nose_x_coordinate) { - nose_x_coordinate = jouk_points[i](0); + nose_x_coordinate = jouk_points[i][0]; nose_index = i; } - if (jouk_points[i](0) > tail_x_coordinate) + if (jouk_points[i][0] > tail_x_coordinate) { - tail_x_coordinate = jouk_points[i](0); + tail_x_coordinate = jouk_points[i][0]; tail_index = i; } } @@ -474,9 +474,9 @@ namespace GridGenerator // move nose to origin auto move_nose_to_origin = [](std::vector> &vector) { - const double nose_x_pos = vector.front()(0); + const double nose_x_pos = vector.front()[0]; for (auto &i : vector) - i(0) -= nose_x_pos; + i[0] -= nose_x_pos; }; move_nose_to_origin(airfoil_1D[1]); @@ -520,10 +520,10 @@ namespace GridGenerator // Calculate radius so that point (x=1|y=0) is enclosed - requirement // for Joukowski transform - const double radius = std::sqrt(center(1) * center(1) + - (1 - center(0)) * (1 - center(0))); + const double radius = std::sqrt(center[1] * center[1] + + (1 - center[0]) * (1 - center[0])); const double radius_test = std::sqrt( - center(1) * center(1) + (1 + center(0)) * (1 + center(0))); + center[1] * center[1] + (1 + center[0]) * (1 + center[0])); // Make sure point (x=-1|y=0) is enclosed by the circle (void)radius_test; AssertThrow( @@ -559,8 +559,8 @@ namespace GridGenerator // transform each point for (unsigned int i = 0; i < circle_points.size(); ++i) { - const double chi = circle_points[i](0); - const double eta = circle_points[i](1); + const double chi = circle_points[i][0]; + const double eta = circle_points[i][1]; const std::complex zeta(chi, eta); const std::complex z = zeta + 1. / zeta; @@ -849,8 +849,8 @@ namespace GridGenerator for (auto it = tria.begin_vertex(); it < tria.end_vertex(); ++it) { auto &point = it->vertex(); - const double xi = point(0); - const double eta = point(1); + const double xi = point[0]; + const double eta = point[1]; // bilinear mapping point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] + @@ -970,7 +970,7 @@ namespace GridGenerator // horizontal offset in order to place coarse-grid node A in the // origin - const Point<2, double> horizontal_offset(A(0), 0.0); + const Point<2, double> horizontal_offset(A[0], 0.0); // Move block 1 so that face BC coincides the x-axis const Point<2, double> trapeze_offset(0.0, @@ -994,7 +994,7 @@ namespace GridGenerator cell->material_id() == id_block_4) // block 1 and 4 { // step 1: rotate block 1 clockwise by gamma and move block - // 1 so that A(0) is on y-axis so that faces AD and BC are + // 1 so that A[0] is on y-axis so that faces AD and BC are // horizontal. This simplifies the computation of the // required indices for interpolation (all x-nodes are // positive) Move trapeze to be in first quadrant by adding @@ -1023,21 +1023,21 @@ namespace GridGenerator const double L = height / std::sin(gamma); const double l_a = std::cos(gamma) * edge_length; const double l_b = trapeze_height * std::tan(gamma); - const double x1 = std::abs(node_(1)) / std::tan(gamma); + const double x1 = std::abs(node_[1]) / std::tan(gamma); const double x2 = L - l_a - l_b; - const double x3 = std::abs(node_(1)) * std::tan(gamma); + const double x3 = std::abs(node_[1]) * std::tan(gamma); const double Dx = x1 + x2 + x3; const double deltax = - (trapeze_height - std::abs(node_(1))) / std::tan(gamma); + (trapeze_height - std::abs(node_[1])) / std::tan(gamma); const double dx = Dx / n_cells_x_0; const double dy = trapeze_height / n_cells_y; const int ix = - static_cast(std::round((node_(0) - deltax) / dx)); + static_cast(std::round((node_[0] - deltax) / dx)); const int iy = - static_cast(std::round(std::abs(node_(1)) / dy)); + static_cast(std::round(std::abs(node_[1]) / dy)); - node_(0) = numbers::PI / 2 * (1.0 * ix) / n_cells_x_0; - node_(1) = height * (1.0 * iy) / n_cells_y; + node_[0] = numbers::PI / 2 * (1.0 * ix) / n_cells_x_0; + node_[1] = height * (1.0 * iy) / n_cells_y; } // step 3: Interpolation between semicircle (of C-Mesh) and @@ -1046,12 +1046,12 @@ namespace GridGenerator const double dx = numbers::PI / 2 / n_cells_x_0; const double dy = height / n_cells_y; const int ix = - static_cast(std::round(node_(0) / dx)); + static_cast(std::round(node_[0] / dx)); const int iy = - static_cast(std::round(node_(1) / dy)); + static_cast(std::round(node_[1] / dy)); const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y); - const double theta = node_(0); + const double theta = node_[0]; const Point<2> p(-height * std::cos(theta) + center_mesh, ((cell->material_id() == id_block_1) ? (height) : @@ -1069,22 +1069,22 @@ namespace GridGenerator { // geometric parameters and indices for interpolation Assert( - (std::abs(D(1) - C(1)) == std::abs(F(1) - E(1))) && - (std::abs(C(1)) == std::abs(E(1))) && - (std::abs(G(1)) == std::abs(I(1))), + (std::abs(D[1] - C[1]) == std::abs(F[1] - E[1])) && + (std::abs(C[1]) == std::abs(E[1])) && + (std::abs(G[1]) == std::abs(I[1])), ExcMessage( "Points D,C,G and E,F,I are not defined symmetric to " "x-axis, which is required to interpolate block 2" " and 5 with same geometric computations.")); - const double l_y = D(1) - C(1); - const double l_h = D(1) - l_y; - const double by = -l_h / length_b1_x * (node(0) - H(0)); + const double l_y = D[1] - C[1]; + const double l_h = D[1] - l_y; + const double by = -l_h / length_b1_x * (node[0] - H[0]); const double dy = (height - by) / n_cells_y; const int iy = static_cast( - std::round((std::abs(node(1)) - by) / dy)); + std::round((std::abs(node[1]) - by) / dy)); const double dx = length_b1_x / n_cells_x_1; const int ix = static_cast( - std::round(std::abs(node(0) - center_mesh) / dx)); + std::round(std::abs(node[0] - center_mesh) / dx)); const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y); // define points on upper/lower horizontal far field side, @@ -1111,9 +1111,9 @@ namespace GridGenerator const double dx = length_b2 / n_cells_x_2; const double dy = height / n_cells_y; const int ix = static_cast( - std::round(std::abs(node(0) - H(0)) / dx)); + std::round(std::abs(node[0] - H[0]) / dx)); const int iy = - static_cast(std::round(std::abs(node(1)) / dy)); + static_cast(std::round(std::abs(node[1]) / dy)); const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y); const double alpha_x = @@ -1121,14 +1121,14 @@ namespace GridGenerator // define on upper/lower horizontal far field side at y = // +/- height, i.e. face GK or IL incline factor to move // points G and H to the right - const Point<2> p1(J(0) - (1 - incline_factor) * length_b2 * + const Point<2> p1(J[0] - (1 - incline_factor) * length_b2 * (alpha_x), ((cell->material_id() == id_block_3) ? (height) : (-height))); // define points on HJ but use tail_y as y-coordinate, in // case last airfoil point has y =/= 0 - const Point<2> p2(J(0) - alpha_x * length_b2, tail_y); + const Point<2> p2(J[0] - alpha_x * length_b2, tail_y); node = p1 * (1 - alpha_y) + p2 * alpha_y; } else @@ -1275,7 +1275,7 @@ namespace GridGenerator tria.begin(); cell != tria.end(); ++cell) - if (cell->center()(0) > 0) + if (cell->center()[0] > 0) cell->set_material_id(1); // boundary indicators are set to // 0 (left) and 1 (right) by default. @@ -1307,17 +1307,17 @@ namespace GridGenerator { const Point center(face->center()); - if (std::abs(center(0) - p1[0]) < epsilon) + if (std::abs(center[0] - p1[0]) < epsilon) face->set_boundary_id(0); - else if (std::abs(center(0) - p2[0]) < epsilon) + else if (std::abs(center[0] - p2[0]) < epsilon) face->set_boundary_id(1); - else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon) + else if (dim > 1 && std::abs(center[1] - p1[1]) < epsilon) face->set_boundary_id(2); - else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon) + else if (dim > 1 && std::abs(center[1] - p2[1]) < epsilon) face->set_boundary_id(3); - else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon) + else if (dim > 2 && std::abs(center[2] - p1[2]) < epsilon) face->set_boundary_id(4); - else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon) + else if (dim > 2 && std::abs(center[2] - p2[2]) < epsilon) face->set_boundary_id(5); else // triangulation says it @@ -1331,7 +1331,7 @@ namespace GridGenerator { types::material_id id = 0; for (unsigned int d = 0; d < dim; ++d) - if (cell->center()(d) > 0) + if (cell->center()[d] > 0) id += (1 << d); cell->set_material_id(id); } @@ -1464,7 +1464,7 @@ namespace GridGenerator continue; double radius = cell->face(f)->center().norm() - center.norm(); - if (std::fabs(cell->face(f)->center()(0)) < + if (std::fabs(cell->face(f)->center()[0]) < eps) // x = 0 set boundary 2 { cell->face(f)->set_boundary_id(2); @@ -1476,7 +1476,7 @@ namespace GridGenerator eps) cell->face(f)->line(j)->set_boundary_id(2); } - else if (std::fabs(cell->face(f)->center()(1)) < + else if (std::fabs(cell->face(f)->center()[1]) < eps) // y = 0 set boundary 3 { cell->face(f)->set_boundary_id(3); @@ -1488,7 +1488,7 @@ namespace GridGenerator eps) cell->face(f)->line(j)->set_boundary_id(3); } - else if (std::fabs(cell->face(f)->center()(2)) < + else if (std::fabs(cell->face(f)->center()[2]) < eps) // z = 0 set boundary 4 { cell->face(f)->set_boundary_id(4); @@ -1543,8 +1543,8 @@ namespace GridGenerator Point p1, p2; for (unsigned int i = 0; i < dim; ++i) { - p1(i) = std::min(p_1(i), p_2(i)); - p2(i) = std::max(p_1(i), p_2(i)); + p1[i] = std::min(p_1[i], p_2[i]); + p2[i] = std::max(p_1[i], p_2[i]); } std::vector> vertices(GeometryInfo::vertices_per_cell); @@ -1558,22 +1558,22 @@ namespace GridGenerator vertices[0] = vertices[1] = p1; vertices[2] = vertices[3] = p2; - vertices[1](0) = p2(0); - vertices[2](0) = p1(0); + vertices[1][0] = p2[0]; + vertices[2][0] = p1[0]; break; case 3: vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1; vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2; - vertices[1](0) = p2(0); - vertices[2](1) = p2(1); - vertices[3](0) = p2(0); - vertices[3](1) = p2(1); + vertices[1][0] = p2[0]; + vertices[2][1] = p2[1]; + vertices[3][0] = p2[0]; + vertices[3][1] = p2[1]; - vertices[4](0) = p1(0); - vertices[4](1) = p1(1); - vertices[5](1) = p1(1); - vertices[6](0) = p1(0); + vertices[4][0] = p1[0]; + vertices[4][1] = p1[1]; + vertices[5][1] = p1[1]; + vertices[6][0] = p1[0]; break; default: @@ -1608,8 +1608,8 @@ namespace GridGenerator Point p1, p2; for (unsigned int i = 0; i < dim; ++i) { - p1(i) = left; - p2(i) = right; + p1[i] = left; + p2[i] = right; } hyper_rectangle(tria, p1, p2, colorize); } @@ -1628,7 +1628,7 @@ namespace GridGenerator Tensor<2, dim> vector_matrix; for (unsigned int d = 0; d < dim; ++d) for (unsigned int c = 1; c <= dim; ++c) - vector_matrix[c - 1][d] = vertices[c](d) - vertices[0](d); + vector_matrix[c - 1][d] = vertices[c][d] - vertices[0][d]; Assert(determinant(vector_matrix) > 0., ExcMessage("Vertices of simplex must form a right handed system")); # endif @@ -2515,8 +2515,8 @@ namespace GridGenerator Point p1, p2; for (unsigned int i = 0; i < dim; ++i) { - p1(i) = std::min(p_1(i), p_2(i)); - p2(i) = std::max(p_1(i), p_2(i)); + p1[i] = std::min(p_1[i], p_2[i]); + p2[i] = std::max(p_1[i], p_2[i]); } // calculate deltas and validate input @@ -2671,9 +2671,9 @@ namespace GridGenerator for (unsigned int i = 0; i < dim; ++i) { - if (p1(i) > p2(i)) + if (p1[i] > p2[i]) { - std::swap(p1(i), p2(i)); + std::swap(p1[i], p2[i]); std::reverse(step_sizes[i].begin(), step_sizes[i].end()); } @@ -2681,7 +2681,7 @@ namespace GridGenerator double x = 0; for (unsigned int j = 0; j < step_sizes.at(i).size(); ++j) x += step_sizes[i][j]; - Assert(std::fabs(x - (p2(i) - p1(i))) <= 1e-12 * std::fabs(x), + Assert(std::fabs(x - (p2[i] - p1[i])) <= 1e-12 * std::fabs(x), ExcMessage( "The sequence of step sizes in coordinate direction " + Utilities::int_to_string(i) + @@ -3142,7 +3142,7 @@ namespace GridGenerator Point p1; Point p2; for (unsigned int d = 0; d < dim; ++d) - p2(d) = 1.; + p2[d] = 1.; // then check that all repetitions // are >= 1, and calculate deltas @@ -3852,7 +3852,7 @@ namespace GridGenerator { Point p; for (unsigned int d = 0; d < dim; ++d) - p(d) = 0.5 * dimensions[d] * + p[d] = 0.5 * dimensions[d] * GeometryInfo::unit_normal_orientation [GeometryInfo::vertex_to_face[i][d]]; points.push_back(p); @@ -3888,7 +3888,7 @@ namespace GridGenerator cells[cell_index].vertices[cellv] = points.size(); Point p = points[cells[cell_index].vertices[ocellv]]; - p(dir) += GeometryInfo::unit_normal_orientation[face] * + p[dir] += GeometryInfo::unit_normal_orientation[face] * dimensions[dir]; points.push_back(p); } @@ -4573,8 +4573,8 @@ namespace GridGenerator // If one the components is the same as the respective // component of the center, then this is part of the plane - if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius || - cell->face(i)->center()(1) < p(1) + 1.e-5 * radius) + if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius || + cell->face(i)->center()[1] < p[1] + 1.e-5 * radius) { cell->face(i)->set_boundary_id(1); cell->face(i)->set_manifold_id(numbers::flat_manifold_id); @@ -4639,7 +4639,7 @@ namespace GridGenerator continue; // If x is zero, then this is part of the plane - if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius) + if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius) { cell->face(i)->set_boundary_id(1); cell->face(i)->set_manifold_id(numbers::flat_manifold_id); @@ -5526,9 +5526,9 @@ namespace GridGenerator // Turn cylinder such that y->x for (auto &vertex : vertices) { - const double h = vertex(1); - vertex(1) = -vertex(0); - vertex(0) = h; + const double h = vertex[1]; + vertex[1] = -vertex[0]; + vertex[0] = h; } std::vector> cell_vertices; @@ -5585,7 +5585,7 @@ namespace GridGenerator for (const unsigned int i : GeometryInfo<3>::face_indices()) if (cell->at_boundary(i)) { - if (cell->face(i)->center()(0) > half_length - tolerance) + if (cell->face(i)->center()[0] > half_length - tolerance) { cell->face(i)->set_boundary_id(2); cell->face(i)->set_manifold_id(numbers::flat_manifold_id); @@ -5602,7 +5602,7 @@ namespace GridGenerator numbers::flat_manifold_id); } } - else if (cell->face(i)->center()(0) < -half_length + tolerance) + else if (cell->face(i)->center()[0] < -half_length + tolerance) { cell->face(i)->set_boundary_id(1); cell->face(i)->set_manifold_id(numbers::flat_manifold_id); @@ -5698,9 +5698,9 @@ namespace GridGenerator continue; // If x,y or z is zero, then this is part of the plane - if (cell->face(i)->center()(0) < center(0) + 1.e-5 * radius || - cell->face(i)->center()(1) < center(1) + 1.e-5 * radius || - cell->face(i)->center()(2) < center(2) + 1.e-5 * radius) + if (cell->face(i)->center()[0] < center[0] + 1.e-5 * radius || + cell->face(i)->center()[1] < center[1] + 1.e-5 * radius || + cell->face(i)->center()[2] < center[2] + 1.e-5 * radius) { cell->face(i)->set_boundary_id(1); cell->face(i)->set_manifold_id(numbers::flat_manifold_id); @@ -5809,7 +5809,7 @@ namespace GridGenerator // If the center is on the plane x=0, this is a planar element. set // its boundary indicator. also set the boundary indicators of the // bounding faces unless both vertices are on the perimeter - if (cell->face(i)->center()(0) < center(0) + 1.e-5 * radius) + if (cell->face(i)->center()[0] < center[0] + 1.e-5 * radius) { cell->face(i)->set_boundary_id(1); cell->face(i)->set_manifold_id(numbers::flat_manifold_id); @@ -6298,7 +6298,7 @@ namespace GridGenerator const Triangulation<3>::face_iterator face = cell->face(i); const Point<3> face_center(face->center()); - if (std::abs(face_center(0) - center(0)) > + if (std::abs(face_center[0] - center[0]) > 1.e-6 * face_center.norm()) { if (std::abs((face_center - center).norm() - inner_radius) < @@ -7306,8 +7306,8 @@ namespace GridGenerator for (const unsigned int f : cell->face_indices()) if (cell->face(f)->at_boundary()) { - const double dx = cell->face(f)->center()(0) - center(0); - const double dy = cell->face(f)->center()(1) - center(1); + const double dx = cell->face(f)->center()[0] - center[0]; + const double dy = cell->face(f)->center()[1] - center[1]; if (colorize) { if (std::abs(dx + outer_radius) < eps) @@ -7525,9 +7525,9 @@ namespace GridGenerator for (const unsigned int f : cell->face_indices()) if (cell->face(f)->at_boundary()) { - const double dx = cell->face(f)->center()(0); - const double dy = cell->face(f)->center()(1); - const double dz = cell->face(f)->center()(2); + const double dx = cell->face(f)->center()[0]; + const double dy = cell->face(f)->center()[1]; + const double dz = cell->face(f)->center()[2]; if (colorize) { @@ -7558,7 +7558,7 @@ namespace GridGenerator else { Point c = cell->face(f)->center(); - c(2) = 0; + c[2] = 0; const double d = c.norm(); if (d - inner_radius < 0) { diff --git a/source/grid/grid_generator_pipe_junction.cc b/source/grid/grid_generator_pipe_junction.cc index f5ba8d62c9..9ad5359c3e 100644 --- a/source/grid/grid_generator_pipe_junction.cc +++ b/source/grid/grid_generator_pipe_junction.cc @@ -242,15 +242,15 @@ namespace Manifold::push_forward(const Point<3> &chart_point) const { // Rotate the orthogonal direction by the given angle. - const double sine_r = chart_point(0) * std::sin(chart_point(1)); - const double cosine_r = chart_point(0) * std::cos(chart_point(1)); + const double sine_r = chart_point[0] * std::sin(chart_point[1]); + const double cosine_r = chart_point[0] * std::cos(chart_point[1]); const Tensor<1, spacedim> intermediate = normal_direction * cosine_r + dxn * sine_r; // Map the axial coordinate back to the pipe segment. const double lambda = - chart_point(2) * compute_z_expansion(cosine_r, sine_r, data); + chart_point[2] * compute_z_expansion(cosine_r, sine_r, data); // Finally, put everything together. return point_on_axis + direction * lambda + intermediate; diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index 7bad170c4f..130fbf8ede 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -212,11 +212,11 @@ GridIn::read_vtk(std::istream &in) { // VTK format always specifies vertex coordinates with 3 components Point<3> x; - in >> x(0) >> x(1) >> x(2); + in >> x[0] >> x[1] >> x[2]; vertices.emplace_back(); for (unsigned int d = 0; d < spacedim; ++d) - vertices.back()(d) = x(d); + vertices.back()[d] = x[d]; } } @@ -683,7 +683,7 @@ GridIn::read_unv(std::istream &in) vertices.emplace_back(); for (unsigned int d = 0; d < spacedim; ++d) - vertices.back()(d) = x[d]; + vertices.back()[d] = x[d]; vertex_indices[vertex_index] = n_vertices; @@ -950,7 +950,7 @@ GridIn::read_ucd(std::istream &in, // store vertex for (unsigned int d = 0; d < spacedim; ++d) - vertices[vertex](d) = x[d]; + vertices[vertex][d] = x[d]; // store mapping; note that // vertices_indices[i] is automatically // created upon first usage @@ -2270,7 +2270,7 @@ GridIn::read_msh(std::istream &in) in >> vertex_number >> x[0] >> x[1] >> x[2]; for (unsigned int d = 0; d < spacedim; ++d) - vertices[global_vertex](d) = x[d]; + vertices[global_vertex][d] = x[d]; // store mapping vertex_indices[vertex_number] = global_vertex; @@ -3239,13 +3239,13 @@ GridIn<2>::read_tecplot(std::istream &in) Utilities::break_text_into_lines(line, 1); char *endptr; for (unsigned int i = 1; i < first_var.size() + 1; ++i) - vertices[i](0) = std::strtod(first_var[i - 1].c_str(), &endptr); + vertices[i][0] = std::strtod(first_var[i - 1].c_str(), &endptr); // if there are many points, the data // for this var might continue in the // next line(s) for (unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j) - in >> vertices[j](next_index); + in >> vertices[j][next_index]; // now we got all values of the first // variable, so increase the counter next_index = 1; @@ -3268,7 +3268,7 @@ GridIn<2>::read_tecplot(std::istream &in) { // we need this line, read it in for (unsigned int j = 1; j < n_vertices + 1; ++j) - in >> vertices[j](next_index); + in >> vertices[j][next_index]; ++next_index; } else @@ -3296,7 +3296,7 @@ GridIn<2>::read_tecplot(std::istream &in) Utilities::break_text_into_lines(line, 1); char *endptr; for (unsigned int d = 0; d < dim; ++d) - vertices[1](d) = + vertices[1][d] = std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr); // read the remaining vertices from the @@ -3310,7 +3310,7 @@ GridIn<2>::read_tecplot(std::istream &in) // of coordinates in the list of // variables for (unsigned int i = 0; i < dim; ++i) - vertices[v](i) = vars[tecplot2deal[i]]; + vertices[v][i] = vars[tecplot2deal[i]]; } } @@ -4005,10 +4005,10 @@ GridIn<2>::debug_output_grid(const std::vector> &cells, const std::vector> &vertices, std::ostream &out) { - double min_x = vertices[cells[0].vertices[0]](0), - max_x = vertices[cells[0].vertices[0]](0), - min_y = vertices[cells[0].vertices[0]](1), - max_y = vertices[cells[0].vertices[0]](1); + double min_x = vertices[cells[0].vertices[0]][0], + max_x = vertices[cells[0].vertices[0]][0], + min_y = vertices[cells[0].vertices[0]][1], + max_y = vertices[cells[0].vertices[0]][1]; for (unsigned int i = 0; i < cells.size(); ++i) { @@ -4016,14 +4016,14 @@ GridIn<2>::debug_output_grid(const std::vector> &cells, { const Point<2> &p = vertices[vertex]; - if (p(0) < min_x) - min_x = p(0); - if (p(0) > max_x) - max_x = p(0); - if (p(1) < min_y) - min_y = p(1); - if (p(1) > max_y) - max_y = p(1); + if (p[0] < min_x) + min_x = p[0]; + if (p[0] > max_x) + max_x = p[0]; + if (p[1] < min_y) + min_y = p[1]; + if (p[1] > max_y) + max_y = p[1]; } out << "# cell " << i << std::endl; @@ -4032,21 +4032,21 @@ GridIn<2>::debug_output_grid(const std::vector> &cells, center += vertices[vertex]; center /= 4; - out << "set label \"" << i << "\" at " << center(0) << ',' << center(1) + out << "set label \"" << i << "\" at " << center[0] << ',' << center[1] << " center" << std::endl; // first two line right direction for (unsigned int f = 0; f < 2; ++f) - out << "set arrow from " << vertices[cells[i].vertices[f]](0) << ',' - << vertices[cells[i].vertices[f]](1) << " to " - << vertices[cells[i].vertices[(f + 1) % 4]](0) << ',' - << vertices[cells[i].vertices[(f + 1) % 4]](1) << std::endl; + out << "set arrow from " << vertices[cells[i].vertices[f]][0] << ',' + << vertices[cells[i].vertices[f]][1] << " to " + << vertices[cells[i].vertices[(f + 1) % 4]][0] << ',' + << vertices[cells[i].vertices[(f + 1) % 4]][1] << std::endl; // other two lines reverse direction for (unsigned int f = 2; f < 4; ++f) - out << "set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]](0) - << ',' << vertices[cells[i].vertices[(f + 1) % 4]](1) << " to " - << vertices[cells[i].vertices[f]](0) << ',' - << vertices[cells[i].vertices[f]](1) << std::endl; + out << "set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]][0] + << ',' << vertices[cells[i].vertices[(f + 1) % 4]][1] << " to " + << vertices[cells[i].vertices[f]][0] << ',' + << vertices[cells[i].vertices[f]][1] << std::endl; out << std::endl; } diff --git a/source/grid/grid_out.cc b/source/grid/grid_out.cc index ff628d015b..b8476649db 100644 --- a/source/grid/grid_out.cc +++ b/source/grid/grid_out.cc @@ -1377,8 +1377,8 @@ GridOut::write_xfig(const Triangulation<2> &tria, cell->vertex(GeometryInfo::ucd_to_deal[k % nv]); for (unsigned int d = 0; d < static_cast(dim); ++d) { - int val = static_cast(1200 * xfig_flags.scaling(d) * - (p(d) - xfig_flags.offset(d))); + int val = static_cast(1200 * xfig_flags.scaling[d] * + (p[d] - xfig_flags.offset[d])); out << '\t' << ((d == 0) ? val : -val); } out << std::endl; @@ -1423,8 +1423,8 @@ GridOut::write_xfig(const Triangulation<2> &tria, ++d) { int val = - static_cast(1200 * xfig_flags.scaling(d) * - (p(d) - xfig_flags.offset(d))); + static_cast(1200 * xfig_flags.scaling[d] * + (p[d] - xfig_flags.offset[d])); out << '\t' << ((d == 0) ? val : -val); } out << std::endl; @@ -4232,7 +4232,7 @@ namespace internal boundary_points[0][0] = 0; boundary_points[n_points - 1][0] = 1; for (unsigned int i = 1; i < n_points - 1; ++i) - boundary_points[i](0) = 1. * i / (n_points - 1); + boundary_points[i][0] = 1. * i / (n_points - 1); std::vector dummy_weights(n_points, 1. / n_points); Quadrature quadrature(boundary_points, dummy_weights); @@ -4349,7 +4349,7 @@ namespace internal boundary_points[0][0] = 0; boundary_points[n_points - 1][0] = 1; for (unsigned int i = 1; i < n_points - 1; ++i) - boundary_points[i](0) = 1. * i / (n_points - 1); + boundary_points[i][0] = 1. * i / (n_points - 1); std::vector dummy_weights(n_points, 1. / n_points); Quadrature<1> quadrature1d(boundary_points, dummy_weights); @@ -4755,8 +4755,8 @@ namespace internal // optimize away this // little kludge line_list.emplace_back( - Point<2>(line->vertex(0)(0), line->vertex(0)(1)), - Point<2>(line->vertex(1)(0), line->vertex(1)(1)), + Point<2>(line->vertex(0)[0], line->vertex(0)[1]), + Point<2>(line->vertex(1)[0], line->vertex(1)[1]), line->user_flag_set(), cell->level()); } @@ -4777,7 +4777,7 @@ namespace internal std::vector> boundary_points(n_points); for (unsigned int i = 0; i < n_points; ++i) - boundary_points[i](0) = 1. * (i + 1) / (n_points + 1); + boundary_points[i][0] = 1. * (i + 1) / (n_points + 1); Quadrature quadrature(boundary_points); Quadrature q_projector( @@ -4798,7 +4798,7 @@ namespace internal if (face->at_boundary()) { Point p0_dim(face->vertex(0)); - Point<2> p0(p0_dim(0), p0_dim(1)); + Point<2> p0(p0_dim[0], p0_dim[1]); // loop over // all pieces @@ -4811,7 +4811,7 @@ namespace internal const Point p1_dim( mapping->transform_unit_to_real_cell( cell, q_projector.point(offset + i))); - const Point<2> p1(p1_dim(0), p1_dim(1)); + const Point<2> p1(p1_dim[0], p1_dim[1]); line_list.emplace_back(p0, p1, @@ -4822,7 +4822,7 @@ namespace internal // generate last piece const Point p1_dim(face->vertex(1)); - const Point<2> p1(p1_dim(0), p1_dim(1)); + const Point<2> p1(p1_dim[0], p1_dim[1]); line_list.emplace_back(p0, p1, face->user_flag_set(), @@ -4914,9 +4914,9 @@ namespace internal // find out minimum and maximum x and // y coordinates to compute offsets // and scaling factors - double x_min = tria.begin_active()->vertex(0)(0); + double x_min = tria.begin_active()->vertex(0)[0]; double x_max = x_min; - double y_min = tria.begin_active()->vertex(0)(1); + double y_min = tria.begin_active()->vertex(0)[1]; double y_max = y_min; unsigned int max_level = line_list.begin()->level; @@ -4924,17 +4924,17 @@ namespace internal line != line_list.end(); ++line) { - x_min = std::min(x_min, line->first(0)); - x_min = std::min(x_min, line->second(0)); + x_min = std::min(x_min, line->first[0]); + x_min = std::min(x_min, line->second[0]); - x_max = std::max(x_max, line->first(0)); - x_max = std::max(x_max, line->second(0)); + x_max = std::max(x_max, line->first[0]); + x_max = std::max(x_max, line->second[0]); - y_min = std::min(y_min, line->first(1)); - y_min = std::min(y_min, line->second(1)); + y_min = std::min(y_min, line->first[1]); + y_min = std::min(y_min, line->second[1]); - y_max = std::max(y_max, line->first(1)); - y_max = std::max(y_max, line->second(1)); + y_max = std::max(y_max, line->first[1]); + y_max = std::max(y_max, line->second[1]); max_level = std::max(max_level, line->level); } @@ -5054,8 +5054,8 @@ namespace internal for (const auto &cell : tria.active_cell_iterators()) { - out << (cell->center()(0) - offset(0)) * scale << ' ' - << (cell->center()(1) - offset(1)) * scale << " m" << '\n' + out << (cell->center()[0] - offset[0]) * scale << ' ' + << (cell->center()[1] - offset[1]) * scale << " m" << '\n' << "[ [(Helvetica) 12.0 0.0 true true ("; if (eps_flags_2.write_cell_number_level) out << cell; @@ -5084,8 +5084,8 @@ namespace internal { treated_vertices.insert(cell->vertex_index(vertex_no)); - out << (cell->vertex(vertex_no)(0) - offset(0)) * scale << ' ' - << (cell->vertex(vertex_no)(1) - offset(1)) * scale + out << (cell->vertex(vertex_no)[0] - offset[0]) * scale << ' ' + << (cell->vertex(vertex_no)[1] - offset[1]) * scale << " m" << '\n' << "[ [(Helvetica) 10.0 0.0 true true (" << cell->vertex_index(vertex_no) << ")] " diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index dad68daa4e..b42d5a57d6 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -364,8 +364,8 @@ namespace GridTools cell->vertex_dof_index(vertex_no, 0), {}, (solve_for_absolute_positions ? - map_iter->second(i) : - map_iter->second(i) - vertex_point[i])); + map_iter->second[i] : + map_iter->second[i] - vertex_point[i])); } } } @@ -398,9 +398,9 @@ namespace GridTools cell->vertex_dof_index(vertex_no, 0); for (unsigned int i = 0; i < dim; ++i) if (solve_for_absolute_positions) - v(i) = us[i](dof_index); + v[i] = us[i](dof_index); else - v(i) += us[i](dof_index); + v[i] += us[i](dof_index); vertex_touched[cell->vertex_index(vertex_no)] = true; } @@ -525,7 +525,7 @@ namespace GridTools // first compute a random shift vector Point shift_vector; for (unsigned int d = 0; d < spacedim; ++d) - shift_vector(d) = uniform_distribution(rng); + shift_vector[d] = uniform_distribution(rng); shift_vector *= factor * minimal_length[global_vertex_no] / std::sqrt(shift_vector.square()); @@ -562,7 +562,7 @@ namespace GridTools // compute a random shift vector Point shift_vector; for (unsigned int d = 0; d < spacedim; ++d) - shift_vector(d) = uniform_distribution(rng); + shift_vector[d] = uniform_distribution(rng); shift_vector *= factor * minimal_length[vertex] / std::sqrt(shift_vector.square()); diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 4c6d30cc3d..46f47c93af 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -2387,7 +2387,7 @@ namespace GridTools if (matrix.m() == spacedim) for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < spacedim; ++j) - distance(i) += matrix(i, j) * point1(j); + distance[i] += matrix(i, j) * point1[j]; else distance = point1; @@ -2399,7 +2399,7 @@ namespace GridTools if (i == direction) continue; - if (std::abs(distance(i)) > 1.e-10) + if (std::abs(distance[i]) > 1.e-10) return false; } diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc index 6f3d4686a6..2e6cff31d4 100644 --- a/source/grid/grid_tools_nontemplates.cc +++ b/source/grid/grid_tools_nontemplates.cc @@ -49,13 +49,13 @@ namespace GridTools { if (vertex_indices.size() == 3) // triangle { - const double x[3] = {all_vertices[vertex_indices[0]](0), - all_vertices[vertex_indices[1]](0), - all_vertices[vertex_indices[2]](0)}; + const double x[3] = {all_vertices[vertex_indices[0]][0], + all_vertices[vertex_indices[1]][0], + all_vertices[vertex_indices[2]][0]}; - const double y[3] = {all_vertices[vertex_indices[0]](1), - all_vertices[vertex_indices[1]](1), - all_vertices[vertex_indices[2]](1)}; + const double y[3] = {all_vertices[vertex_indices[0]][1], + all_vertices[vertex_indices[1]][1], + all_vertices[vertex_indices[2]][1]}; return 0.5 * ((x[0] - x[2]) * (y[1] - y[0]) - (x[1] - x[0]) * (y[0] - y[2])); @@ -99,15 +99,15 @@ namespace GridTools additional optimization: divide by 2 only one time */ - const double x[4] = {all_vertices[vertex_indices[0]](0), - all_vertices[vertex_indices[1]](0), - all_vertices[vertex_indices[2]](0), - all_vertices[vertex_indices[3]](0)}; + const double x[4] = {all_vertices[vertex_indices[0]][0], + all_vertices[vertex_indices[1]][0], + all_vertices[vertex_indices[2]][0], + all_vertices[vertex_indices[3]][0]}; - const double y[4] = {all_vertices[vertex_indices[0]](1), - all_vertices[vertex_indices[1]](1), - all_vertices[vertex_indices[2]](1), - all_vertices[vertex_indices[3]](1)}; + const double y[4] = {all_vertices[vertex_indices[0]][1], + all_vertices[vertex_indices[1]][1], + all_vertices[vertex_indices[2]][1], + all_vertices[vertex_indices[3]][1]}; return (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] - x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]) / @@ -274,30 +274,30 @@ namespace GridTools AssertDimension(vertex_indices.size(), GeometryInfo<3>::vertices_per_cell); - const double x[8] = {all_vertices[vertex_indices[0]](0), - all_vertices[vertex_indices[1]](0), - all_vertices[vertex_indices[2]](0), - all_vertices[vertex_indices[3]](0), - all_vertices[vertex_indices[4]](0), - all_vertices[vertex_indices[5]](0), - all_vertices[vertex_indices[6]](0), - all_vertices[vertex_indices[7]](0)}; - const double y[8] = {all_vertices[vertex_indices[0]](1), - all_vertices[vertex_indices[1]](1), - all_vertices[vertex_indices[2]](1), - all_vertices[vertex_indices[3]](1), - all_vertices[vertex_indices[4]](1), - all_vertices[vertex_indices[5]](1), - all_vertices[vertex_indices[6]](1), - all_vertices[vertex_indices[7]](1)}; - const double z[8] = {all_vertices[vertex_indices[0]](2), - all_vertices[vertex_indices[1]](2), - all_vertices[vertex_indices[2]](2), - all_vertices[vertex_indices[3]](2), - all_vertices[vertex_indices[4]](2), - all_vertices[vertex_indices[5]](2), - all_vertices[vertex_indices[6]](2), - all_vertices[vertex_indices[7]](2)}; + const double x[8] = {all_vertices[vertex_indices[0]][0], + all_vertices[vertex_indices[1]][0], + all_vertices[vertex_indices[2]][0], + all_vertices[vertex_indices[3]][0], + all_vertices[vertex_indices[4]][0], + all_vertices[vertex_indices[5]][0], + all_vertices[vertex_indices[6]][0], + all_vertices[vertex_indices[7]][0]}; + const double y[8] = {all_vertices[vertex_indices[0]][1], + all_vertices[vertex_indices[1]][1], + all_vertices[vertex_indices[2]][1], + all_vertices[vertex_indices[3]][1], + all_vertices[vertex_indices[4]][1], + all_vertices[vertex_indices[5]][1], + all_vertices[vertex_indices[6]][1], + all_vertices[vertex_indices[7]][1]}; + const double z[8] = {all_vertices[vertex_indices[0]][2], + all_vertices[vertex_indices[1]][2], + all_vertices[vertex_indices[2]][2], + all_vertices[vertex_indices[3]][2], + all_vertices[vertex_indices[4]][2], + all_vertices[vertex_indices[5]][2], + all_vertices[vertex_indices[6]][2], + all_vertices[vertex_indices[7]][2]}; /* This is the same Maple script as in the barycenter method above diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 4a6a71bb1b..0a27044222 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -1113,13 +1113,13 @@ CylindricalManifold::push_forward( ExcMessage("CylindricalManifold can only be used for spacedim==3!")); // Rotate the orthogonal direction by the given angle - const double sine_r = std::sin(chart_point(1)) * chart_point(0); - const double cosine_r = std::cos(chart_point(1)) * chart_point(0); + const double sine_r = std::sin(chart_point[1]) * chart_point[0]; + const double cosine_r = std::cos(chart_point[1]) * chart_point[0]; const Tensor<1, spacedim> intermediate = normal_direction * cosine_r + dxn * sine_r; // Finally, put everything together. - return point_on_axis + direction * chart_point(2) + intermediate; + return point_on_axis + direction * chart_point[2] + intermediate; } @@ -1135,8 +1135,8 @@ CylindricalManifold::push_forward_gradient( Tensor<2, 3> derivatives; // Rotate the orthogonal direction by the given angle - const double sine = std::sin(chart_point(1)); - const double cosine = std::cos(chart_point(1)); + const double sine = std::sin(chart_point[1]); + const double cosine = std::cos(chart_point[1]); const Tensor<1, spacedim> intermediate = normal_direction * cosine + dxn * sine; @@ -1517,9 +1517,9 @@ template Point<3> TorusManifold::pull_back(const Point<3> &p) const { - double x = p(0); - double z = p(1); - double y = p(2); + double x = p[0]; + double z = p[1]; + double y = p[2]; double phi = std::atan2(y, x); double theta = std::atan2(z, std::sqrt(x * x + y * y) - R); double w = @@ -1535,9 +1535,9 @@ template Point<3> TorusManifold::push_forward(const Point<3> &chart_point) const { - double phi = chart_point(0); - double theta = chart_point(1); - double w = chart_point(2); + double phi = chart_point[0]; + double theta = chart_point[1]; + double w = chart_point[2]; return {std::cos(phi) * R + r * w * std::cos(theta) * std::cos(phi), r * w * std::sin(theta), @@ -1575,9 +1575,9 @@ TorusManifold::push_forward_gradient(const Point<3> &chart_point) const { DerivativeForm<1, spacedim, spacedim> DX; - double phi = chart_point(0); - double theta = chart_point(1); - double w = chart_point(2); + double phi = chart_point[0]; + double theta = chart_point[1]; + double w = chart_point[2]; DX[0][0] = -std::sin(phi) * R - r * w * std::cos(theta) * std::sin(phi); DX[0][1] = -r * w * std::sin(theta) * std::cos(phi); diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index b9a9c0cfd5..4e6d80b33c 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -186,14 +186,14 @@ namespace C(array(1..2, [xs, ys]), optimized); */ - const double x[4] = {accessor.vertex(0)(0), - accessor.vertex(1)(0), - accessor.vertex(2)(0), - accessor.vertex(3)(0)}; - const double y[4] = {accessor.vertex(0)(1), - accessor.vertex(1)(1), - accessor.vertex(2)(1), - accessor.vertex(3)(1)}; + const double x[4] = {accessor.vertex(0)[0], + accessor.vertex(1)[0], + accessor.vertex(2)[0], + accessor.vertex(3)[0]}; + const double y[4] = {accessor.vertex(0)[1], + accessor.vertex(1)[1], + accessor.vertex(2)[1], + accessor.vertex(3)[1]}; const double t1 = x[0] * x[1]; const double t3 = x[0] * x[0]; const double t5 = x[1] * x[1]; @@ -314,30 +314,30 @@ namespace z[] back to the standard ordering. */ - const double x[8] = {accessor.vertex(0)(0), - accessor.vertex(1)(0), - accessor.vertex(5)(0), - accessor.vertex(4)(0), - accessor.vertex(2)(0), - accessor.vertex(3)(0), - accessor.vertex(7)(0), - accessor.vertex(6)(0)}; - const double y[8] = {accessor.vertex(0)(1), - accessor.vertex(1)(1), - accessor.vertex(5)(1), - accessor.vertex(4)(1), - accessor.vertex(2)(1), - accessor.vertex(3)(1), - accessor.vertex(7)(1), - accessor.vertex(6)(1)}; - const double z[8] = {accessor.vertex(0)(2), - accessor.vertex(1)(2), - accessor.vertex(5)(2), - accessor.vertex(4)(2), - accessor.vertex(2)(2), - accessor.vertex(3)(2), - accessor.vertex(7)(2), - accessor.vertex(6)(2)}; + const double x[8] = {accessor.vertex(0)[0], + accessor.vertex(1)[0], + accessor.vertex(5)[0], + accessor.vertex(4)[0], + accessor.vertex(2)[0], + accessor.vertex(3)[0], + accessor.vertex(7)[0], + accessor.vertex(6)[0]}; + const double y[8] = {accessor.vertex(0)[1], + accessor.vertex(1)[1], + accessor.vertex(5)[1], + accessor.vertex(4)[1], + accessor.vertex(2)[1], + accessor.vertex(3)[1], + accessor.vertex(7)[1], + accessor.vertex(6)[1]}; + const double z[8] = {accessor.vertex(0)[2], + accessor.vertex(1)[2], + accessor.vertex(5)[2], + accessor.vertex(4)[2], + accessor.vertex(2)[2], + accessor.vertex(3)[2], + accessor.vertex(7)[2], + accessor.vertex(6)[2]}; double s1, s2, s3, s4, s5, s6, s7, s8; diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc index 7fd19bf129..e14dd7093d 100644 --- a/source/non_matching/quadrature_generator.cc +++ b/source/non_matching/quadrature_generator.cc @@ -1927,7 +1927,7 @@ namespace NonMatching const Point vertex0 = box.vertex(GeometryInfo::face_to_cell_vertices(face_index, 0)); - const double coordinate_value = vertex0(face_normal_direction); + const double coordinate_value = vertex0[face_normal_direction]; const Functions::CoordinateRestriction face_restriction( level_set, face_normal_direction, coordinate_value); diff --git a/source/numerics/data_out_rotation.cc b/source/numerics/data_out_rotation.cc index edd84d2112..05845bcedd 100644 --- a/source/numerics/data_out_rotation.cc +++ b/source/numerics/data_out_rotation.cc @@ -148,8 +148,8 @@ DataOutRotation::build_one_patch( { case 1: { - const double r1 = (*cell)->vertex(0)(0), - r2 = (*cell)->vertex(1)(0); + const double r1 = (*cell)->vertex(0)[0], + r2 = (*cell)->vertex(1)[0]; Assert(r1 >= 0, ExcRadialVariableHasNegativeValues(r1)); Assert(r2 >= 0, ExcRadialVariableHasNegativeValues(r2)); @@ -169,19 +169,19 @@ DataOutRotation::build_one_patch( const Point v = (*cell)->vertex(vertex); // make sure that the radial variable is nonnegative - Assert(v(0) >= 0, ExcRadialVariableHasNegativeValues(v(0))); + Assert(v[0] >= 0, ExcRadialVariableHasNegativeValues(v[0])); // now set the vertices of the patch my_patches[angle].vertices[vertex] = - v(0) * angle_directions[angle]; - my_patches[angle].vertices[vertex][0] = v(1); + v[0] * angle_directions[angle]; + my_patches[angle].vertices[vertex][0] = v[1]; my_patches[angle] .vertices[vertex + GeometryInfo::vertices_per_cell] = - v(0) * angle_directions[angle + 1]; + v[0] * angle_directions[angle + 1]; my_patches[angle] .vertices[vertex + GeometryInfo::vertices_per_cell] - [0] = v(1); + [0] = v[1]; } break; diff --git a/source/numerics/data_out_stack.cc b/source/numerics/data_out_stack.cc index 49fe7dd2c4..3bccab9038 100644 --- a/source/numerics/data_out_stack.cc +++ b/source/numerics/data_out_stack.cc @@ -339,34 +339,34 @@ DataOutStack::build_patches(const unsigned int nnnn_subdivisions) { case 1: patch->vertices[0] = - Point(cell->vertex(0)(0), parameter - parameter_step); + Point(cell->vertex(0)[0], parameter - parameter_step); patch->vertices[1] = - Point(cell->vertex(1)(0), parameter - parameter_step); - patch->vertices[2] = Point(cell->vertex(0)(0), parameter); - patch->vertices[3] = Point(cell->vertex(1)(0), parameter); + Point(cell->vertex(1)[0], parameter - parameter_step); + patch->vertices[2] = Point(cell->vertex(0)[0], parameter); + patch->vertices[3] = Point(cell->vertex(1)[0], parameter); break; case 2: - patch->vertices[0] = Point(cell->vertex(0)(0), - cell->vertex(0)(1), + patch->vertices[0] = Point(cell->vertex(0)[0], + cell->vertex(0)[1], parameter - parameter_step); - patch->vertices[1] = Point(cell->vertex(1)(0), - cell->vertex(1)(1), + patch->vertices[1] = Point(cell->vertex(1)[0], + cell->vertex(1)[1], parameter - parameter_step); - patch->vertices[2] = Point(cell->vertex(2)(0), - cell->vertex(2)(1), + patch->vertices[2] = Point(cell->vertex(2)[0], + cell->vertex(2)[1], parameter - parameter_step); - patch->vertices[3] = Point(cell->vertex(3)(0), - cell->vertex(3)(1), + patch->vertices[3] = Point(cell->vertex(3)[0], + cell->vertex(3)[1], parameter - parameter_step); patch->vertices[4] = - Point(cell->vertex(0)(0), cell->vertex(0)(1), parameter); + Point(cell->vertex(0)[0], cell->vertex(0)[1], parameter); patch->vertices[5] = - Point(cell->vertex(1)(0), cell->vertex(1)(1), parameter); + Point(cell->vertex(1)[0], cell->vertex(1)[1], parameter); patch->vertices[6] = - Point(cell->vertex(2)(0), cell->vertex(2)(1), parameter); + Point(cell->vertex(2)[0], cell->vertex(2)[1], parameter); patch->vertices[7] = - Point(cell->vertex(3)(0), cell->vertex(3)(1), parameter); + Point(cell->vertex(3)[0], cell->vertex(3)[1], parameter); break; default: diff --git a/source/particles/particle.cc b/source/particles/particle.cc index 26bb888eb6..e0685cbe77 100644 --- a/source/particles/particle.cc +++ b/source/particles/particle.cc @@ -87,12 +87,12 @@ namespace Particles Point location; for (unsigned int i = 0; i < spacedim; ++i) - location(i) = *pdata++; + location[i] = *pdata++; set_location(location); Point reference_location; for (unsigned int i = 0; i < dim; ++i) - reference_location(i) = *pdata++; + reference_location[i] = *pdata++; set_reference_location(reference_location); // See if there are properties to load @@ -259,12 +259,12 @@ namespace Particles Point location; for (unsigned int i = 0; i < spacedim; ++i) - location(i) = *pdata++; + location[i] = *pdata++; set_location(location); Point reference_location; for (unsigned int i = 0; i < dim; ++i) - reference_location(i) = *pdata++; + reference_location[i] = *pdata++; set_reference_location(reference_location); // See if there are properties to load -- 2.39.5