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.;
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.;
{
(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;
}
{
(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;
MeshWorker::IntegrationInfo<dim> &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),
Assert(dim >= 2, ExcInternalError());
Point<dim> 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))
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.;
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<dim>();
- values[i](0) = -points[i](1);
+ values[i][0] = -points[i][1];
}
}
}
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];
}
{
// Create the rectangular domain.
Point<dim> 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<unsigned int> repetitions(dim, 1);
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;
{
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.;
{
(void)component;
AssertIndexRange(component, 1);
- return 25. - p(0);
+ return 25. - p[0];
}
template <int dim>
{
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)
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]);
}
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;
}
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]));
}
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;
}
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]);
}
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();
{
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;
}
(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;
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;
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());
// 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 <int dim>
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) /
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;
}
}
Point<2> operator()(const Point<2> &in) const
{
- return {in(0), trans(in(1))};
+ return {in[0], trans(in[1])};
}
};
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);
}
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) +
{
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 =
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);
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);
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)
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)
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)
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);
f << "# y u_x u_y" << std::endl;
Point<dim> 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<double> 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);
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);
}
}
double Solution<dim>::value(const Point<dim> &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;
}
Tensor<1, dim> Solution<dim>::gradient(const Point<dim> &p,
const unsigned int component) const
{
- return Point<dim>(-2 * p(component));
+ return Point<dim>(-2 * p[component]);
}
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
}
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
}
{
#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;
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);
// 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);
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);
// these areas. Note that upon construction of the Point objects, all
// components are set to zero.
Point<dim> 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)
{
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());
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());
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());
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());
* 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]))));
}
};
Number sum = internal::NumberType<Number>::value(0.0);
for (unsigned int i = 0; i < dim; ++i)
{
- const Number diff = static_cast<Number>(this->values[i]) - p(i);
+ const Number diff = static_cast<Number>(this->values[i]) - p[i];
sum += numbers::NumberTraits<Number>::abs_square(diff);
}
std::vector<double> 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];
}
*reinterpret_cast<Tensor<1, dim> *>(&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
{
std::vector<double> 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];
std::vector<std::vector<double>> v(dim, std::vector<double>(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<order, dim> derivative;
switch (order)
// 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)
// 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;
}
{
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
{
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
std::vector<Point<dim>> 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<dim> q_coarse(q_points_coarse,
fine.get_JxW_values());
FEValues<dim, spacedim> coarse(mapping,
std::vector<Point<dim>> 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<dim> q_coarse(q_points_coarse, fine.get_JxW_values());
FEValues<dim, spacedim> coarse(
coarse_cell->reference_cell()
const Point<spacedim> &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]));
}
// 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 -
Point<dim> 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 *
Point<spacedim> location;
for (unsigned int i = 0; i < spacedim; ++i)
- location(i) = *pdata++;
+ location[i] = *pdata++;
set_location(location);
Point<dim> 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
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;
{
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];
}
}
{
float data[dim];
for (unsigned int d = 0; d < dim; ++d)
- data[d] = p(d);
+ data[d] = p[d];
stream.write(reinterpret_cast<const char *>(data), dim * sizeof(*data));
}
else
{
for (unsigned int d = 0; d < dim; ++d)
- stream << p(d) << '\t';
+ stream << p[d] << '\t';
stream << '\n';
}
}
{
Assert(selected_component != numbers::invalid_unsigned_int,
ExcNotInitialized());
- stream << p(selected_component) << ' ';
+ stream << p[selected_component] << ' ';
}
{
Assert(selected_component != numbers::invalid_unsigned_int,
ExcNotInitialized());
- stream << p(selected_component) << '\n';
+ stream << p[selected_component] << '\n';
}
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 ";
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)
// / /\ / /\ /
// ____/ __/ \/ / \__/
- // { 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)
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;
}
}
// 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';
}
}
<< " 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';
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());
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 )
-(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)
{
// 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
// 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
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;
}
}
// 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.;
for (unsigned int k = 0; k < n; ++k)
{
const Point<dim> &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);
}
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);
for (unsigned int k = 0; k < n; ++k)
{
const Point<dim> &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);
}
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
for (unsigned int k = 0; k < n; ++k)
{
const Point<dim> &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);
}
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);
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))
{
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))
{
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);
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);
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);
Q1WedgeFunction<dim>::value(const Point<dim> &p, const unsigned int) const
{
Assert(dim >= 2, ExcInternalError());
- return p(0) * p(1);
+ return p[0] * p[1];
}
for (unsigned int i = 0; i < points.size(); ++i)
{
const Point<dim> &p = points[i];
- values[i] = p(0) * p(1);
+ values[i] = p[0] * p[1];
}
}
for (unsigned int i = 0; i < points.size(); ++i)
{
const Point<dim> &p = points[i];
- values[i](0) = p(0) * p(1);
+ values[i](0) = p[0] * p[1];
}
}
{
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;
}
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];
}
}
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];
}
}
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());
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:
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());
}
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());
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());
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());
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());
}
{
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());
}
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());
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());
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:
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;
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:
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;
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());
}
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());
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());
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());
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());
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());
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;
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;
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());
}
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());
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());
}
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());
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;
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;
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.;
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.;
{
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.;
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.);
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.);
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.);
{
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;
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;
{
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.);
SlitSingularityFunction<dim>::value(const Point<dim> &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;
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;
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;
SlitSingularityFunction<dim>::gradient(const Point<dim> &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.);
for (unsigned int i = 0; i < points.size(); ++i)
{
const Point<dim> &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.);
ExcDimensionMismatch(gradients[i].size(), 1));
const Point<dim> &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.);
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;
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;
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;
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.);
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.);
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.);
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());
double
JumpFunction<dim>::value(const Point<dim> &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);
}
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);
}
}
double
JumpFunction<dim>::laplacian(const Point<dim> &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);
}
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);
}
Tensor<1, dim>
JumpFunction<dim>::gradient(const Point<dim> &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;
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;
{
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) ?
{
const Point<dim> &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) ?
RayleighKotheVortex<dim>::vector_value(const Point<dim> &point,
Vector<double> &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) *
AssertIndexRange(component_in_dim_plus_1, dim + 1);
Point<dim + 1> 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<dim>(
component_in_dim_plus_1, d);
- output(component_to_write_to) = point(d);
+ output[component_to_write_to] = point[d];
}
return output;
const unsigned int component) const
{
const Point<dim + 1> 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);
}
const unsigned int component) const
{
const Point<dim + 1> 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);
const unsigned int component) const
{
const Point<dim + 1> 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);
{
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;
{
// acos returns the angle in the range [0,\pi]
if (scoord[0] > std::numeric_limits<double>::min())
- scoord[2] = std::acos(position(2) / scoord[0]);
+ scoord[2] = std::acos(position[2] / scoord[0]);
else
scoord[2] = 0.0;
}
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;
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;
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
std::vector<Point<1>> points(n + 1);
const double one_over_n = 1. / n;
for (unsigned int k = 0; k <= n; ++k)
- points[k](0) = static_cast<double>(k) * one_over_n;
+ points[k][0] = static_cast<double>(k) * one_over_n;
return points;
}
} // namespace LagrangeEquidistantImplementation
// 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;
}
std::vector<double> 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)
std::vector<double> 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)
{
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)
// directions.
Point<dim> 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,
PolynomialsAdini<dim>::compute_value(const unsigned int i,
const Point<dim> &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 +
PolynomialsAdini<dim>::compute_grad(const unsigned int i,
const Point<dim> &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 +
PolynomialsAdini<dim>::compute_grad_grad(const unsigned int i,
const Point<dim> &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 +
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];
}
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.;
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]);
}
{
// 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.;
{
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];
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];
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];
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.;
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];
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];
// shifted one step in positive direction
Point<dim> p;
- p(0) = unit_point(1);
- p(1) = unit_point(0);
+ p[0] = unit_point[1];
+ p[1] = unit_point[0];
std::vector<double> p_values((values.empty()) ? 0 : n_basis);
std::vector<Tensor<1, dim>> p_grads((grads.empty()) ? 0 : n_basis);
std::vector<Tensor<2, dim>> 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,
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());
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
{
// through all directions.
Point<dim> 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,
// 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)
{
// 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);
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]);
}
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);
}
switch (face_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
+ q_points[p] = Point<dim>(quadrature.point(p)[0], 0);
break;
case 1:
q_points[p] =
- Point<dim>(1 - quadrature.point(p)(0), quadrature.point(p)(0));
+ Point<dim>(1 - quadrature.point(p)[0], quadrature.point(p)[0]);
break;
case 2:
- q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0));
+ q_points[p] = Point<dim>(0, 1 - quadrature.point(p)[0]);
break;
default:
Assert(false, ExcInternalError());
switch (face_no)
{
case 0:
- q_points[p] = Point<dim>(0, quadrature.point(p)(0));
+ q_points[p] = Point<dim>(0, quadrature.point(p)[0]);
break;
case 1:
- q_points[p] = Point<dim>(1, quadrature.point(p)(0));
+ q_points[p] = Point<dim>(1, quadrature.point(p)[0]);
break;
case 2:
- q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
+ q_points[p] = Point<dim>(quadrature.point(p)[0], 0);
break;
case 3:
- q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
+ q_points[p] = Point<dim>(quadrature.point(p)[0], 1);
break;
default:
Assert(false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
+ q_points[p] = Point<dim>(quadrature.point(p)[0] / 2, 0);
break;
case 1:
q_points[p] =
- Point<dim>(0.5 + quadrature.point(p)(0) / 2, 0);
+ Point<dim>(0.5 + quadrature.point(p)[0] / 2, 0);
break;
default:
Assert(false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(1 - quadrature.point(p)(0) / 2,
- quadrature.point(p)(0) / 2);
+ q_points[p] = Point<dim>(1 - quadrature.point(p)[0] / 2,
+ quadrature.point(p)[0] / 2);
break;
case 1:
- q_points[p] = Point<dim>(0.5 - quadrature.point(p)(0) / 2,
- 0.5 + quadrature.point(p)(0) / 2);
+ q_points[p] = Point<dim>(0.5 - quadrature.point(p)[0] / 2,
+ 0.5 + quadrature.point(p)[0] / 2);
break;
default:
Assert(false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(0, 1 - quadrature.point(p)(0) / 2);
+ q_points[p] = Point<dim>(0, 1 - quadrature.point(p)[0] / 2);
break;
case 1:
q_points[p] =
- Point<dim>(0, 0.5 - quadrature.point(p)(0) / 2);
+ Point<dim>(0, 0.5 - quadrature.point(p)[0] / 2);
break;
default:
Assert(false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
+ q_points[p] = Point<dim>(0, quadrature.point(p)[0] / 2);
break;
case 1:
q_points[p] =
- Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
+ Point<dim>(0, quadrature.point(p)[0] / 2 + 0.5);
break;
default:
Assert(false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
+ q_points[p] = Point<dim>(1, quadrature.point(p)[0] / 2);
break;
case 1:
q_points[p] =
- Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
+ Point<dim>(1, quadrature.point(p)[0] / 2 + 0.5);
break;
default:
Assert(false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
+ q_points[p] = Point<dim>(quadrature.point(p)[0] / 2, 0);
break;
case 1:
q_points[p] =
- Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
+ Point<dim>(quadrature.point(p)[0] / 2 + 0.5, 0);
break;
default:
Assert(false, ExcInternalError());
switch (subface_no)
{
case 0:
- q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
+ q_points[p] = Point<dim>(quadrature.point(p)[0] / 2, 1);
break;
case 1:
q_points[p] =
- Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
+ Point<dim>(quadrature.point(p)[0] / 2 + 0.5, 1);
break;
default:
Assert(false, ExcInternalError());
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);
// 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);
{
// 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);
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);
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());
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());
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());
++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] =
{
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());
}
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]);
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;
}
std::vector<double> 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];
}
std::vector<double> 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];
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;
}
// coordinate direction
ndarray<double, dim, 2> 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)
ndarray<double, dim, 3> 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)
// 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;
}
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)
// 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;
}
{
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
{
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
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 <>
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];
}
}
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
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;
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;
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;
{
Point<codim> 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;
}
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)
{
this->generalized_support_points[q_point +
(i + 4 * j) * n_edge_points] =
- Point<dim>(i, reference_edge_quadrature.point(q_point)(0), j);
+ Point<dim>(i, reference_edge_quadrature.point(q_point)[0], j);
this->generalized_support_points[q_point + (i + 4 * j + 2) *
n_edge_points] =
- Point<dim>(reference_edge_quadrature.point(q_point)(0), i, j);
+ Point<dim>(reference_edge_quadrature.point(q_point)[0], i, j);
this->generalized_support_points[q_point + (i + 2 * (j + 4)) *
n_edge_points] =
- Point<dim>(i, j, reference_edge_quadrature.point(q_point)(0));
+ Point<dim>(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.
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]);
}
}
{
this->generalized_support_points[q_point +
(i + 4 * j) * n_edge_points] =
- Point<dim>(i, reference_edge_quadrature.point(q_point)(0), j);
+ Point<dim>(i, reference_edge_quadrature.point(q_point)[0], j);
this->generalized_support_points[q_point + (i + 4 * j + 2) *
n_edge_points] =
- Point<dim>(reference_edge_quadrature.point(q_point)(0), i, j);
+ Point<dim>(reference_edge_quadrature.point(q_point)[0], i, j);
this->generalized_support_points[q_point + (i + 2 * (j + 4)) *
n_edge_points] =
- Point<dim>(i, j, reference_edge_quadrature.point(q_point)(0));
+ Point<dim>(i, j, reference_edge_quadrature.point(q_point)[0]);
}
}
}
{
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<dim> 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);
else
{
Point<dim> 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);
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<double> system_matrix(deg, deg);
{
const double weight = edge_quadrature.weight(q_point);
const Point<dim> quadrature_point_0(
- i, edge_quadrature_points[q_point](0));
+ i, edge_quadrature_points[q_point][0]);
const Point<dim> 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<dim> quadrature_point_2(
- i, 2.0 * edge_quadrature_points[q_point](0));
+ i, 2.0 * edge_quadrature_points[q_point][0]);
tmp(0) =
weight *
quadrature_point_0,
1);
quadrature_point_2 = Point<dim>(
- 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(
Point<dim> 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 *
this->shape_value_component(
(i + 2) * this->degree, quadrature_point_1, 0);
quadrature_point_2 = Point<dim>(
- 2.0 * edge_quadrature_points[q_point](0) - 1.0,
+ 2.0 * edge_quadrature_points[q_point][0] - 1.0,
i);
tmp(3) =
weight *
{
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;
{
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]);
}
}
{
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<dim> 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);
else
{
const Point<dim> 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);
}
}
- else if (quadrature_points[q_point](1) < 0.5)
+ else if (quadrature_points[q_point][1] < 0.5)
{
const Point<dim> 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,
else
{
const Point<dim> 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,
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)
{
{
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<dim> 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,
weight *
this->shape_value_component(dof, quadrature_point, 1);
quadrature_point =
- Point<dim>(2.0 * edge_quadrature_points[q_point](0),
+ Point<dim>(2.0 * edge_quadrature_points[q_point][0],
i,
j);
this->restriction[index][2 * (i + 2 * j)](
quadrature_point =
Point<dim>(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) +=
for (unsigned int j = 0; j < 2; ++j)
{
Point<dim> 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,
weight *
this->shape_value_component(dof, quadrature_point, 1);
quadrature_point = Point<dim>(
- 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<dim>(
- 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 *
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<double> system_matrix(deg, deg);
{
const double weight = edge_quadrature.weight(q_point);
const Point<dim> quadrature_point_0(
- i, edge_quadrature_points[q_point](0), j);
+ i, edge_quadrature_points[q_point][0], j);
const Point<dim> quadrature_point_1(
- edge_quadrature_points[q_point](0), i, j);
+ edge_quadrature_points[q_point][0], i, j);
const Point<dim> 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<dim> 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(
quadrature_point_0,
1);
quadrature_point_3 = Point<dim>(
- 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(
quadrature_point_1,
0);
quadrature_point_3 = Point<dim>(
- 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(
Point<dim> 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 *
quadrature_point_1,
0);
quadrature_point_3 = Point<dim>(
- 2.0 * edge_quadrature_points[q_point](0) - 1.0,
+ 2.0 * edge_quadrature_points[q_point][0] - 1.0,
i,
j);
tmp(3) =
quadrature_point_3 = Point<dim>(
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(
{
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;
{
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]);
}
}
{
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<dim> 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<dim>(
- 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<dim>(
- 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);
{
Point<dim> 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(
tmp(3) += 2.0 * this->shape_value_component(
dof, quadrature_point_0, 2);
quadrature_point_0 = Point<dim>(
- 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<dim>(
- 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(
}
}
- else if (face_quadrature_points[q_point](1) < 0.5)
+ else if (face_quadrature_points[q_point][1] < 0.5)
{
Point<dim> 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<dim>(
- 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<dim>(
- 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);
{
Point<dim> 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<dim>(
- 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<dim>(
- 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);
const Point<dim> 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<dim> 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<dim> 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)
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)
{
{
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]);
}
}
}
{
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<dim> 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);
else
{
const Point<dim> 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);
}
}
- else if (quadrature_points[q_point](2) < 0.5)
+ else if (quadrature_points[q_point][2] < 0.5)
{
const Point<dim> 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);
else
{
const Point<dim> 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);
}
}
- 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<dim> 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);
else
{
const Point<dim> 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);
}
}
- else if (quadrature_points[q_point](2) < 0.5)
+ else if (quadrature_points[q_point][2] < 0.5)
{
const Point<dim> 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,
else
{
const Point<dim> 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,
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)
{
++q_point)
{
const Point<dim> 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) *
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<double> system_matrix(source_fe.degree - 1,
{
const Point<dim> 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<dim> 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,
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);
{
Point<dim> 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);
this->shape_value_component(
this->face_to_cell_index(dof, 4), quadrature_point, 1);
quadrature_point =
- Point<dim>(0.5 * (edge_quadrature_points[q_point](0) +
+ Point<dim>(0.5 * (edge_quadrature_points[q_point][0] +
shifts[subface][0]),
0.5 * (i + shifts[subface][1]),
0.0);
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<double> system_matrix(source_fe.degree - 1,
{
Point<dim> 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<dim> quadrature_point_1(
- i, edge_quadrature_points[q_point](0), 0.0);
+ i, edge_quadrature_points[q_point][0], 0.0);
tmp(i) =
weight *
i * source_fe.degree, quadrature_point_1, 1));
quadrature_point_0 =
Point<dim>(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<dim>(edge_quadrature_points[q_point](0),
+ Point<dim>(edge_quadrature_points[q_point][0],
i,
0.0);
tmp(i + 2) =
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<dim>::lines_per_face;
{
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]);
}
}
{
Point<dim> 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),
quadrature_point,
1);
quadrature_point =
- Point<dim>(quadrature_points[q_point](0),
- quadrature_points[q_point](1),
+ Point<dim>(quadrature_points[q_point][0],
+ quadrature_points[q_point][1],
0.0);
for (unsigned int i = 0; i < 2; ++i)
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)
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]);
}
}
}
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<double> system_matrix_inv(this->degree - 1,
this->degree - 1);
legendre_polynomials[i].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::lines_per_cell *
- n_edge_points](0)) *
+ n_edge_points][0]) *
lobatto_polynomials[j + 2].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::lines_per_cell *
- n_edge_points](1)) *
+ n_edge_points][1]) *
lobatto_polynomials_grad[k].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::lines_per_cell *
- n_edge_points](0)) *
+ n_edge_points][0]) *
lobatto_polynomials[l + 2].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::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);
lobatto_polynomials_grad[i].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::lines_per_cell *
- n_edge_points](0)) *
+ n_edge_points][0]) *
lobatto_polynomials[j + 2].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::lines_per_cell *
- n_edge_points](1));
+ n_edge_points][1]);
}
solution.reinit(system_matrix.m());
lobatto_polynomials_grad[i].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::lines_per_cell *
- n_edge_points](1)) *
+ n_edge_points][1]) *
lobatto_polynomials[j + 2].value(
this->generalized_support_points
[q_point + GeometryInfo<dim>::lines_per_cell *
- n_edge_points](0));
+ n_edge_points][0]);
}
system_matrix_inv.vmult(solution, system_rhs);
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<double> system_matrix_inv(this->degree - 1,
this->degree - 1);
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);
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](0)) *
+ n_face_points][0]) *
lobatto_polynomials[j + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](1)) *
+ n_face_points][1]) *
lobatto_polynomials[k + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](2)) *
+ n_face_points][2]) *
lobatto_polynomials_grad[l].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](0)) *
+ n_face_points][0]) *
lobatto_polynomials[m + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](1)) *
+ n_face_points][1]) *
lobatto_polynomials[n + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::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);
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](0)) *
+ n_face_points][0]) *
lobatto_polynomials[j + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](1)) *
+ n_face_points][1]) *
lobatto_polynomials[k + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](2));
+ n_face_points][2]);
}
solution.reinit(system_rhs.size());
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](1)) *
+ n_face_points][1]) *
lobatto_polynomials[j + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](0)) *
+ n_face_points][0]) *
lobatto_polynomials[k + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](2));
+ n_face_points][2]);
}
system_matrix_inv.vmult(solution, system_rhs);
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](2)) *
+ n_face_points][2]) *
lobatto_polynomials[j + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](0)) *
+ n_face_points][0]) *
lobatto_polynomials[k + 2].value(
this->generalized_support_points
[q_point +
GeometryInfo<dim>::lines_per_cell *
n_edge_points +
GeometryInfo<dim>::faces_per_cell *
- n_face_points](1));
+ n_face_points][1]);
}
system_matrix_inv.vmult(solution, system_rhs);
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<double, 4, 3> 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;
}
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)
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)
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]);
}
}
, 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
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;
for (unsigned int j = 0; j < q_deg + 1; ++j)
{
Point<dim> 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);
const double step = 1. / q_deg;
std::vector<Point<dim - 2>> 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<dim - 2> qline(line_support_points);
// auxiliary points in 2d
for (unsigned int k = 0; k < dim - 1; ++k)
{
const int coord_int =
- static_cast<int>(constraint_points[i](k) * interval + 0.25);
- constraint_point(k) = 1. * coord_int / interval;
+ static_cast<int>(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
//
// 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)
fine.JxW(q);
Point<dim> 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);
// 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;
}
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;
// 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;
}
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;
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;
}
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;
for (unsigned int q = 0; q < face_quadrature.size(); ++q)
{
this->generalized_support_points[0 * face_quadrature.size() + q] =
- dealii::Point<dim>(0, 1 - face_quadrature.point(q)(0));
+ dealii::Point<dim>(0, 1 - face_quadrature.point(q)[0]);
this->generalized_support_points[1 * face_quadrature.size() + q] =
- dealii::Point<dim>(1, 1 - face_quadrature.point(q)(0));
+ dealii::Point<dim>(1, 1 - face_quadrature.point(q)[0]);
this->generalized_support_points[2 * face_quadrature.size() + q] =
- dealii::Point<dim>(face_quadrature.point(q)(0), 0);
+ dealii::Point<dim>(face_quadrature.point(q)[0], 0);
this->generalized_support_points[3 * face_quadrature.size() + q] =
- dealii::Point<dim>(face_quadrature.point(q)(0), 1);
+ dealii::Point<dim>(face_quadrature.point(q)[0], 1);
}
}
if (dim == 2)
{
if (unit_normal_direction == 0)
- return Point<dim - 1>{unit_cell_pt(1)};
+ return Point<dim - 1>{unit_cell_pt[1]};
else if (unit_normal_direction == 1)
- return Point<dim - 1>{unit_cell_pt(0)};
+ return Point<dim - 1>{unit_cell_pt[0]};
}
else if (dim == 3)
{
if (unit_normal_direction == 0)
- return Point<dim - 1>{unit_cell_pt(1), unit_cell_pt(2)};
+ return Point<dim - 1>{unit_cell_pt[1], unit_cell_pt[2]};
else if (unit_normal_direction == 1)
- return Point<dim - 1>{unit_cell_pt(0), unit_cell_pt(2)};
+ return Point<dim - 1>{unit_cell_pt[0], unit_cell_pt[2]};
else if (unit_normal_direction == 2)
- return Point<dim - 1>{unit_cell_pt(0), unit_cell_pt(1)};
+ return Point<dim - 1>{unit_cell_pt[0], unit_cell_pt[1]};
}
// We should never get here
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());
{
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];
}
}
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());
Point<dim> 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;
}
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());
// 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;
}
{
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;
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 /
, 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",
// 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;
}
}
// move nose to origin
auto move_nose_to_origin = [](std::vector<Point<2>> &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]);
// 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(
// 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<double> zeta(chi, eta);
const std::complex<double> z = zeta + 1. / zeta;
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] +
// 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,
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
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<int>(std::round((node_(0) - deltax) / dx));
+ static_cast<int>(std::round((node_[0] - deltax) / dx));
const int iy =
- static_cast<int>(std::round(std::abs(node_(1)) / dy));
+ static_cast<int>(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
const double dx = numbers::PI / 2 / n_cells_x_0;
const double dy = height / n_cells_y;
const int ix =
- static_cast<int>(std::round(node_(0) / dx));
+ static_cast<int>(std::round(node_[0] / dx));
const int iy =
- static_cast<int>(std::round(node_(1) / dy));
+ static_cast<int>(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) :
{
// 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<int>(
- 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<int>(
- 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,
const double dx = length_b2 / n_cells_x_2;
const double dy = height / n_cells_y;
const int ix = static_cast<int>(
- std::round(std::abs(node(0) - H(0)) / dx));
+ std::round(std::abs(node[0] - H[0]) / dx));
const int iy =
- static_cast<int>(std::round(std::abs(node(1)) / dy));
+ static_cast<int>(std::round(std::abs(node[1]) / dy));
const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
const double alpha_x =
// 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
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.
{
const Point<spacedim> 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
{
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);
}
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);
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);
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);
Point<spacedim> 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<Point<spacedim>> vertices(GeometryInfo<dim>::vertices_per_cell);
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:
Point<dim> 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);
}
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
Point<spacedim> 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
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());
}
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) +
Point<spacedim> p1;
Point<spacedim> 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
{
Point<spacedim> p;
for (unsigned int d = 0; d < dim; ++d)
- p(d) = 0.5 * dimensions[d] *
+ p[d] = 0.5 * dimensions[d] *
GeometryInfo<dim>::unit_normal_orientation
[GeometryInfo<dim>::vertex_to_face[i][d]];
points.push_back(p);
cells[cell_index].vertices[cellv] = points.size();
Point<spacedim> p = points[cells[cell_index].vertices[ocellv]];
- p(dir) += GeometryInfo<dim>::unit_normal_orientation[face] *
+ p[dir] += GeometryInfo<dim>::unit_normal_orientation[face] *
dimensions[dir];
points.push_back(p);
}
// 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);
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);
// 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<std::vector<int>> cell_vertices;
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);
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);
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);
// 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);
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) <
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)
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)
{
else
{
Point<dim> c = cell->face(f)->center();
- c(2) = 0;
+ c[2] = 0;
const double d = c.norm();
if (d - inner_radius < 0)
{
Manifold<dim, spacedim>::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;
{
// 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];
}
}
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;
// 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
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;
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;
{
// 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
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
// 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]];
}
}
const std::vector<Point<2>> &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)
{
{
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;
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;
}
cell->vertex(GeometryInfo<dim>::ucd_to_deal[k % nv]);
for (unsigned int d = 0; d < static_cast<unsigned int>(dim); ++d)
{
- int val = static_cast<int>(1200 * xfig_flags.scaling(d) *
- (p(d) - xfig_flags.offset(d)));
+ int val = static_cast<int>(1200 * xfig_flags.scaling[d] *
+ (p[d] - xfig_flags.offset[d]));
out << '\t' << ((d == 0) ? val : -val);
}
out << std::endl;
++d)
{
int val =
- static_cast<int>(1200 * xfig_flags.scaling(d) *
- (p(d) - xfig_flags.offset(d)));
+ static_cast<int>(1200 * xfig_flags.scaling[d] *
+ (p[d] - xfig_flags.offset[d]));
out << '\t' << ((d == 0) ? val : -val);
}
out << std::endl;
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<double> dummy_weights(n_points, 1. / n_points);
Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
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<double> dummy_weights(n_points, 1. / n_points);
Quadrature<1> quadrature1d(boundary_points, dummy_weights);
// 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());
}
std::vector<Point<dim - 1>> 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<dim - 1> quadrature(boundary_points);
Quadrature<dim> q_projector(
if (face->at_boundary())
{
Point<dim> 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
const Point<dim> 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,
// generate last piece
const Point<dim> 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(),
// 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;
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);
}
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;
{
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) << ")] "
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]));
}
}
}
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;
}
// first compute a random shift vector
Point<spacedim> 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());
// compute a random shift vector
Point<spacedim> 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());
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;
if (i == direction)
continue;
- if (std::abs(distance(i)) > 1.e-10)
+ if (std::abs(distance[i]) > 1.e-10)
return false;
}
{
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]));
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]) /
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
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;
}
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;
Point<3>
TorusManifold<dim>::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 =
Point<3>
TorusManifold<dim>::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),
{
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);
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];
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;
const Point<dim> vertex0 =
box.vertex(GeometryInfo<dim>::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<dim - 1> face_restriction(
level_set, face_normal_direction, coordinate_value);
{
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));
const Point<dim> 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<dim>::vertices_per_cell] =
- v(0) * angle_directions[angle + 1];
+ v[0] * angle_directions[angle + 1];
my_patches[angle]
.vertices[vertex + GeometryInfo<dim>::vertices_per_cell]
- [0] = v(1);
+ [0] = v[1];
}
break;
{
case 1:
patch->vertices[0] =
- Point<dim + 1>(cell->vertex(0)(0), parameter - parameter_step);
+ Point<dim + 1>(cell->vertex(0)[0], parameter - parameter_step);
patch->vertices[1] =
- Point<dim + 1>(cell->vertex(1)(0), parameter - parameter_step);
- patch->vertices[2] = Point<dim + 1>(cell->vertex(0)(0), parameter);
- patch->vertices[3] = Point<dim + 1>(cell->vertex(1)(0), parameter);
+ Point<dim + 1>(cell->vertex(1)[0], parameter - parameter_step);
+ patch->vertices[2] = Point<dim + 1>(cell->vertex(0)[0], parameter);
+ patch->vertices[3] = Point<dim + 1>(cell->vertex(1)[0], parameter);
break;
case 2:
- patch->vertices[0] = Point<dim + 1>(cell->vertex(0)(0),
- cell->vertex(0)(1),
+ patch->vertices[0] = Point<dim + 1>(cell->vertex(0)[0],
+ cell->vertex(0)[1],
parameter - parameter_step);
- patch->vertices[1] = Point<dim + 1>(cell->vertex(1)(0),
- cell->vertex(1)(1),
+ patch->vertices[1] = Point<dim + 1>(cell->vertex(1)[0],
+ cell->vertex(1)[1],
parameter - parameter_step);
- patch->vertices[2] = Point<dim + 1>(cell->vertex(2)(0),
- cell->vertex(2)(1),
+ patch->vertices[2] = Point<dim + 1>(cell->vertex(2)[0],
+ cell->vertex(2)[1],
parameter - parameter_step);
- patch->vertices[3] = Point<dim + 1>(cell->vertex(3)(0),
- cell->vertex(3)(1),
+ patch->vertices[3] = Point<dim + 1>(cell->vertex(3)[0],
+ cell->vertex(3)[1],
parameter - parameter_step);
patch->vertices[4] =
- Point<dim + 1>(cell->vertex(0)(0), cell->vertex(0)(1), parameter);
+ Point<dim + 1>(cell->vertex(0)[0], cell->vertex(0)[1], parameter);
patch->vertices[5] =
- Point<dim + 1>(cell->vertex(1)(0), cell->vertex(1)(1), parameter);
+ Point<dim + 1>(cell->vertex(1)[0], cell->vertex(1)[1], parameter);
patch->vertices[6] =
- Point<dim + 1>(cell->vertex(2)(0), cell->vertex(2)(1), parameter);
+ Point<dim + 1>(cell->vertex(2)[0], cell->vertex(2)[1], parameter);
patch->vertices[7] =
- Point<dim + 1>(cell->vertex(3)(0), cell->vertex(3)(1), parameter);
+ Point<dim + 1>(cell->vertex(3)[0], cell->vertex(3)[1], parameter);
break;
default:
Point<spacedim> location;
for (unsigned int i = 0; i < spacedim; ++i)
- location(i) = *pdata++;
+ location[i] = *pdata++;
set_location(location);
Point<dim> 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
Point<spacedim> location;
for (unsigned int i = 0; i < spacedim; ++i)
- location(i) = *pdata++;
+ location[i] = *pdata++;
set_location(location);
Point<dim> 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
projects / dealii.git / commitdiff