case 1:
return p;
case 2:
- return Point<dim>(p(0) * (1 + p(1)), p(1) * (1 + p(0)));
+ return Point<dim>(p[0] * (1 + p[1]), p[1] * (1 + p[0]));
case 3:
- return Point<dim>(p(0) * (1 + p(1)) * (1 + p(2)),
- p(1) * (1 + p(0)) * (1 + p(2)),
- p(2) * (1 + p(0)) * (1 + p(1)));
+ return Point<dim>(p[0] * (1 + p[1]) * (1 + p[2]),
+ p[1] * (1 + p[0]) * (1 + p[2]),
+ p[2] * (1 + p[0]) * (1 + p[1]));
default:
AssertThrow(false, ExcNotImplemented());
return Point<dim>();
double
AutoSinExp<dim>::value(const Point<dim> &p, const unsigned int) const
{
- return std::sin(2 * p(0)) * std::exp(3 * p(1));
+ return std::sin(2 * p[0]) * std::exp(3 * p[1]);
}
ExactSinExp<dim>::gradient(const Point<dim> &p, const unsigned int) const
{
Tensor<1, dim> grad;
- grad[0] = 2 * std::cos(2 * p(0)) * std::exp(3 * p(1));
- grad[1] = 3 * std::sin(2 * p(0)) * std::exp(3 * p(1));
+ grad[0] = 2 * std::cos(2 * p[0]) * std::exp(3 * p[1]);
+ grad[1] = 3 * std::sin(2 * p[0]) * std::exp(3 * p[1]);
return grad;
}
Point<dim> wave_vector;
for (unsigned int d = 0; d < dim; ++d)
- wave_vector(d) = d + 2.;
+ wave_vector[d] = d + 2.;
Functions::FourierSineFunction<dim> f(wave_vector);
{
deallog << "Direction " << d << std::endl;
Point<dim> dir;
- dir(d) = 1.;
+ dir[d] = 1.;
deallog.push("Euler");
FunctionDerivative<dim> df(f, dir, 1.e-4);
check_derivative_order(gradients, df, quadrature, d, 2);
for (unsigned int ix = 0; ix < 2; ++ix)
{
if (dim > 0)
- patches[0].vertices[vertex_number](0) = -1. + 2. * ix;
+ patches[0].vertices[vertex_number][0] = -1. + 2. * ix;
if (dim > 1)
- patches[0].vertices[vertex_number](1) = -1. + 2. * iy;
+ patches[0].vertices[vertex_number][1] = -1. + 2. * iy;
if (dim > 2)
- patches[0].vertices[vertex_number](2) = -1. + 2. * iz;
+ patches[0].vertices[vertex_number][2] = -1. + 2. * iz;
++vertex_number;
}
for (const unsigned int i : GeometryInfo<dim>::face_indices())
for (unsigned int ix = 0; ix <= sub; ++ix)
{
if (dim > 0)
- points[vertex_number](0) = -1. + ix * h;
+ points[vertex_number][0] = -1. + ix * h;
if (dim > 1)
- points[vertex_number](1) = -1. + iy * h;
+ points[vertex_number][1] = -1. + iy * h;
if (dim > 2)
- points[vertex_number](2) = -1. + iz * h;
+ points[vertex_number][2] = -1. + iz * h;
++vertex_number;
}
Point<dim> point;
for (int i = 0; i < dim; ++i)
- point(i) = i;
+ point[i] = i;
deallog << "->value:" << std::endl;
PrintTensor<rank, dim>::print_tensor(foo->value(point));
points.push_back(point);
for (int i = 0; i < dim; ++i)
- point(i) = dim - i;
+ point[i] = dim - i;
points.push_back(point);
std::vector<Tensor<rank, dim>> tensors;
auto p_3 = Point<dim, Number>::unit_vector(0);
check[3] = p_3.norm_square();
- auto entry_1 = p_1(0);
+ auto entry_1 = p_1[0];
check[4] = entry_1;
- p_1(0) = Number{1.};
+ p_1[0] = Number{1.};
check[5] = p_1.norm_square();
auto p_4 = p_1 + Tensor<1, dim, Number>{};
check[6] = p_4.norm_square();
p[i] = i;
for (unsigned int i = 0; i < dim; ++i)
- deallog << p(i) << ' ';
+ deallog << p[i] << ' ';
deallog << std::endl;
}
Point<3> p(point);
for (unsigned int i = 0; i < 3; ++i)
- deallog << p(i) << ' ';
+ deallog << p[i] << ' ';
deallog << std::endl;
}
Point<2> p(point);
for (unsigned int i = 0; i < 2; ++i)
- deallog << p(i) << ' ';
+ deallog << p[i] << ' ';
deallog << std::endl;
}
Point<1> p(point);
- deallog << p(0) << ' ';
+ deallog << p[0] << ' ';
deallog << std::endl;
}
for (unsigned int i = 0; i < dim; ++i)
- deallog << p_float(i) << ' ';
+ deallog << p_float[i] << ' ';
deallog << std::endl;
for (unsigned int i = 0; i < dim; ++i)
- deallog << p_double(i) << ' ';
+ deallog << p_double[i] << ' ';
deallog << std::endl;
for (unsigned int i = 0; i < dim; ++i)
- deallog << p_complex_float(i) << ' ';
+ deallog << p_complex_float[i] << ' ';
deallog << std::endl;
for (unsigned int i = 0; i < dim; ++i)
- deallog << p_complex_double(i) << ' ';
+ deallog << p_complex_double[i] << ' ';
deallog << std::endl;
}
p_complex_double = p_double;
for (unsigned int i = 0; i < dim; ++i)
- deallog << p_double(i) << ' ';
+ deallog << p_double[i] << ' ';
deallog << std::endl;
for (unsigned int i = 0; i < dim; ++i)
- deallog << p_complex_float(i) << ' ';
+ deallog << p_complex_float[i] << ' ';
deallog << std::endl;
for (unsigned int i = 0; i < dim; ++i)
- deallog << p_complex_double(i) << ' ';
+ deallog << p_complex_double[i] << ' ';
deallog << std::endl;
}
}
template <int dim>
constexpr bool val =
Point<dim>((Point<dim>::unit_vector(0) + Point<dim>::unit_vector(0)) -
- 2. * Point<dim>::unit_vector(0))(0) == 0.;
+ 2. * Point<dim>::unit_vector(0))[0] == 0.;
int
double sum = 0.;
for (unsigned int i = 0; i < gauss.size(); ++i)
{
- double x = gauss.point(i)(0);
+ double x = gauss.point(i)[0];
double P1 = p1.value(x);
double P2 = p2.value(x);
sum += gauss.weight(i) * P1 * P2;
for (unsigned int k = 0; k < x.size(); ++k)
{
deallog << '.';
- const double y = p[i].value(x[k](0));
+ const double y = p[i].value(x[k][0]);
if (i == k)
{
if (std::fabs(y - 1.) > 2.e-10)
std::vector<Point<1>> x(n + 1);
const double h = 1. / n;
for (unsigned int i = 0; i <= n; ++i)
- x[i](0) = h * i;
+ x[i][0] = h * i;
check_interpolation(p, x);
}
for (unsigned int k = 0; k < x.size(); ++k)
{
deallog << '.';
- const double y = p[i].value(x[k](0));
+ const double y = p[i].value(x[k][0]);
if (i == k)
{
if (std::fabs(y - 1.) > 1e-13)
std::vector<Point<1>> x(n + 1);
const double h = 1. / n;
for (unsigned int i = 0; i <= n; ++i)
- x[i](0) = h * i;
+ x[i][0] = h * i;
check_interpolation(p, x);
check_constant(p);
deallog << std::endl;
PolynomialsBDM<dim> p3(3);
PolynomialsBDM<dim> p4(4);
- x(0) = 2.;
- x(1) = 3.;
+ x[0] = 2.;
+ x[1] = 3.;
if (dim > 2)
- x(2) = 4;
+ x[2] = 4;
check_point(x, p1);
check_point(x, p2);
PolynomialsRaviartThomas<dim> p2(2);
PolynomialsRaviartThomas<dim> p3(3);
- x(0) = 2.;
- x(1) = 3.;
+ x[0] = 2.;
+ x[1] = 3.;
if (dim > 2)
- x(2) = 4;
+ x[2] = 4;
check_point(x, p0);
check_point(x, p1);
{
Point<dim> p1;
Point<dim> p2;
- p1(0) = 1.;
- p2(0) = 7.;
+ p1[0] = 1.;
+ p2[0] = 7.;
if (dim > 1)
{
- p1(1) = 3;
- p2(1) = -5.;
+ p1[1] = 3;
+ p2[1] = -5.;
}
if (dim > 2)
{
- p1(2) = 0;
- p2(2) = 10.;
+ p1[2] = 0;
+ p2[2] = 10.;
}
Quadrature<dim> q = QProjector<dim>::project_to_line(
ReferenceCells::get_hypercube<dim>(), quadrature, p1, p2);
long double f = 1.;
for (unsigned int x = 0; x < quadrature.size(); ++x)
{
- f = std::pow(static_cast<long double>(points[x](0)), i * 1.0L);
+ f = std::pow(static_cast<long double>(points[x][0]), i * 1.0L);
quadrature_int += f * static_cast<long double>(weights[x]);
}
err = std::fabs(quadrature_int - exact_monomials[i]);
long double f = 1.;
for (unsigned int x = 0; x < quadrature.size(); ++x)
{
- f = std::pow(static_cast<long double>(points[x](0)), i * 1.0L);
+ f = std::pow(static_cast<long double>(points[x][0]), i * 1.0L);
quadrature_int += f * static_cast<long double>(weights[x]);
}
err = std::fabs(quadrature_int - exact_monomials[i]);
switch (dim)
{
case 3:
- f *= std::pow(static_cast<double>(points[x](2)), i * 1.0);
+ f *= std::pow(static_cast<double>(points[x][2]), i * 1.0);
case 2:
- f *= std::pow(static_cast<double>(points[x](1)), i * 1.0);
+ f *= std::pow(static_cast<double>(points[x][1]), i * 1.0);
case 1:
- f *= std::pow(static_cast<double>(points[x](0)), i * 1.0);
+ f *= std::pow(static_cast<double>(points[x][0]), i * 1.0);
}
quadrature_int += f * weights[x];
}
switch (dim)
{
case 3:
- f *= std::pow((long double)points[x](2), i * 1.0L);
+ f *= std::pow((long double)points[x][2], i * 1.0L);
case 2:
- f *= std::pow((long double)points[x](1), i * 1.0L);
+ f *= std::pow((long double)points[x][1], i * 1.0L);
case 1:
- f *= std::pow((long double)points[x](0), i * 1.0L);
+ f *= std::pow((long double)points[x][0], i * 1.0L);
}
quadrature_int += f * weights[x];
}
switch (dim)
{
case 3:
- f *= std::pow(static_cast<double>(points[x](2)), i * 1.0);
+ f *= std::pow(static_cast<double>(points[x][2]), i * 1.0);
case 2:
- f *= std::pow(static_cast<double>(points[x](1)), i * 1.0);
+ f *= std::pow(static_cast<double>(points[x][1]), i * 1.0);
case 1:
- f *= std::pow(static_cast<double>(points[x](0)), i * 1.0);
+ f *= std::pow(static_cast<double>(points[x][0]), i * 1.0);
}
quadrature_int += f * weights[x];
}
bool in_order = true;
for (unsigned int x = 1; x < quadrature.size(); ++x)
{
- if (points[x](0) <= points[x - 1](0))
+ if (points[x][0] <= points[x - 1][0])
in_order = false;
}
if (!in_order)
switch (dim)
{
case 3:
- f *= std::pow((long double)points[x](2), i * 1.0L);
+ f *= std::pow((long double)points[x][2], i * 1.0L);
case 2:
- f *= std::pow((long double)points[x](1), i * 1.0L);
+ f *= std::pow((long double)points[x][1], i * 1.0L);
case 1:
- f *= std::pow((long double)points[x](0), i * 1.0L);
+ f *= std::pow((long double)points[x][0], i * 1.0L);
}
quadrature_int += f * weights[x];
}
Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2));
for (unsigned int i = 0; i < 2; ++i)
- values(i) = p(i) * p(i);
+ values(i) = p[i] * p[i];
}
virtual double
value(const Point<2> &p, const unsigned int) const
{
- return std::sin(3.14159 * p(0)) * std::sin(3.14159 * p(1));
+ return std::sin(3.14159 * p[0]) * std::sin(3.14159 * p[1]);
}
};
++v)
{
const Point<dim> vertex = face->vertex(v);
- const Tensor<1, dim> tangent_1({-vertex(2), 0., vertex(0)});
+ const Tensor<1, dim> tangent_1({-vertex[2], 0., vertex[0]});
const Tensor<1, dim> tangent_2 = vertex - Point<dim>(0, 3, 0);
// get the normal vector and test it
{
double delta = 0.05;
double x, y, r;
- x = p(0);
- y = p(1);
+ x = p[0];
+ y = p[1];
r = std::sqrt(x * x + y * y);
return 0.5 * (1 - std::tanh((r - 0.5) / (2 * M_SQRT2 * delta)));
}
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
{
unsigned int material = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
material |= (1 << d);
AssertThrow(material < (1 << dim), ExcInternalError());
TestFunction<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
- double val = base[0].value(p(0));
+ double val = base[0].value(p[0]);
for (unsigned int i = 1; i < dim; ++i)
- val *= base[i].value(p(i));
+ val *= base[i].value(p[i]);
return val;
}
grid_in1.read_msh(input_file1);
Point<2> ePos;
- ePos(0) = 0.0653630060373507487669897386695;
- ePos(1) = 1125.59175030825804242340382189;
+ ePos[0] = 0.0653630060373507487669897386695;
+ ePos[1] = 1125.59175030825804242340382189;
MappingQ<2> mapping(1);
MappingQ<2> &mapping2 = StaticMappingQ1<2>::mapping;
deallog.precision(16);
Point<2> ePos;
- ePos(0) = 0.0653630060373507487669897386695;
- ePos(1) = 1125.59175030825804242340382189;
+ ePos[0] = 0.0653630060373507487669897386695;
+ ePos[1] = 1125.59175030825804242340382189;
MappingQ<2> mapping(1);
MappingQ<2> &mapping2 = StaticMappingQ1<2>::mapping;
GridGenerator::hyper_cube(tr);
Point<2> p;
- p(0) = -0.1;
- p(1) = 0.5;
+ p[0] = -0.1;
+ p[1] = 0.5;
MappingQ<2> mapping(1);
for (; cell != endc; ++cell)
{
Point<2> cell_center = cell->center();
- if (std::abs(cell_center(0) - 1500) < 550)
+ if (std::abs(cell_center[0] - 1500) < 550)
{
cell->set_refine_flag();
}
value(const Point<dim> &p, const unsigned int component = 0) const override
{
if (dim == 3)
- return std::sin(p(0) + 1.) * std::sin(p(1) + 2.) * std::sin(p(2) + 3.);
- return std::sin(p(0) + 1.) * std::sin(p(1) + 2.);
+ return std::sin(p[0] + 1.) * std::sin(p[1] + 2.) * std::sin(p[2] + 3.);
+ return std::sin(p[0] + 1.) * std::sin(p[1] + 2.);
}
};
Vector<double> value2(1);
Point<2> point1;
- point1(0) = -numbers::PI + 2. * i / n_points + eps;
- point1(1) = -numbers::PI;
+ point1[0] = -numbers::PI + 2. * i / n_points + eps;
+ point1[1] = -numbers::PI;
Point<2> point2;
- point2(0) = -numbers::PI + 2. * i / n_points + eps;
- point2(1) = numbers::PI;
+ point2[0] = -numbers::PI + 2. * i / n_points + eps;
+ point2[1] = numbers::PI;
VectorTools::point_value(dof_handler, solution, point1, value1);
VectorTools::point_value(dof_handler, solution, point2, value2);
Vector<double> value2(1);
Point<2> point1;
- point1(1) = -numbers::PI + 2. * i / n_points + eps;
- point1(0) = -numbers::PI;
+ point1[1] = -numbers::PI + 2. * i / n_points + eps;
+ point1[0] = -numbers::PI;
Point<2> point2;
- point2(1) = -numbers::PI + 2. * i / n_points + eps;
- point2(0) = numbers::PI;
+ point2[1] = -numbers::PI + 2. * i / n_points + eps;
+ point2[0] = numbers::PI;
VectorTools::point_value(dof_handler, solution, point1, value1);
VectorTools::point_value(dof_handler, solution, point2, value2);
value(const Point<dim> &p, const unsigned int component = 0) const override
{
if (dim == 3)
- return std::sin(p(0) + 1.) * std::sin(p(1) + 2.) * std::sin(p(2) + 3.);
- return std::sin(p(0) + 1.) * std::sin(p(1) + 2.);
+ return std::sin(p[0] + 1.) * std::sin(p[1] + 2.) * std::sin(p[2] + 3.);
+ return std::sin(p[0] + 1.) * std::sin(p[1] + 2.);
}
};
Vector<double> value2(1);
Point<3> point1;
- point1(0) = -numbers::PI + 2. * i / n_points + eps;
- point1(1) = -numbers::PI;
- point1(2) = -numbers::PI + 2. * j / n_points + eps;
+ point1[0] = -numbers::PI + 2. * i / n_points + eps;
+ point1[1] = -numbers::PI;
+ point1[2] = -numbers::PI + 2. * j / n_points + eps;
Point<3> point2;
- point2(0) = -numbers::PI + 2. * i / n_points + eps;
- point2(1) = numbers::PI;
- point2(2) = -numbers::PI + 2. * j / n_points + eps;
+ point2[0] = -numbers::PI + 2. * i / n_points + eps;
+ point2[1] = numbers::PI;
+ point2[2] = -numbers::PI + 2. * j / n_points + eps;
VectorTools::point_value(dof_handler, solution, point1, value1);
VectorTools::point_value(dof_handler, solution, point2, value2);
Vector<double> value2(1);
Point<3> point1;
- point1(2) = -numbers::PI + 2. * j / n_points + eps;
- point1(1) = -numbers::PI + 2. * i / n_points + eps;
- point1(0) = -numbers::PI;
+ point1[2] = -numbers::PI + 2. * j / n_points + eps;
+ point1[1] = -numbers::PI + 2. * i / n_points + eps;
+ point1[0] = -numbers::PI;
Point<3> point2;
- point2(2) = -numbers::PI + 2. * j / n_points + eps;
- point2(1) = -numbers::PI + 2. * i / n_points + eps;
- point2(0) = numbers::PI;
+ point2[2] = -numbers::PI + 2. * j / n_points + eps;
+ point2[1] = -numbers::PI + 2. * i / n_points + eps;
+ point2[0] = numbers::PI;
VectorTools::point_value(dof_handler, solution, point1, value1);
VectorTools::point_value(dof_handler, solution, point2, value2);
Vector<double> value2(1);
Point<3> point1;
- point1(0) = -numbers::PI + 2. * j / n_points + eps;
- point1(1) = -numbers::PI + 2. * i / n_points + eps;
- point1(2) = -numbers::PI;
+ point1[0] = -numbers::PI + 2. * j / n_points + eps;
+ point1[1] = -numbers::PI + 2. * i / n_points + eps;
+ point1[2] = -numbers::PI;
Point<3> point2;
- point2(0) = -numbers::PI + 2. * j / n_points + eps;
- point2(1) = -numbers::PI + 2. * i / n_points + eps;
- point2(2) = numbers::PI;
+ point2[0] = -numbers::PI + 2. * j / n_points + eps;
+ point2[1] = -numbers::PI + 2. * i / n_points + eps;
+ point2[2] = numbers::PI;
VectorTools::point_value(dof_handler, solution, point1, value1);
VectorTools::point_value(dof_handler, solution, point2, value2);
{
Tensor<1, dim> return_value;
for (unsigned int i = 0; i < dim; ++i)
- return_value[i] = 2 * (component + 1) * p(component);
+ return_value[i] = 2 * (component + 1) * p[component];
return return_value;
}
{
Tensor<1, dim> return_value;
for (unsigned int i = 0; i < dim; ++i)
- return_value[i] = 2 * (component + 1) * p(component);
+ return_value[i] = 2 * (component + 1) * p[component];
return return_value;
}
{
Tensor<1, dim> return_value;
for (unsigned int i = 0; i < dim; ++i)
- return_value[i] = 2 * (component + 1) * p(component);
+ return_value[i] = 2 * (component + 1) * p[component];
return return_value;
}
{
Tensor<1, dim> return_value;
for (unsigned int i = 0; i < dim; ++i)
- return_value[i] = 2 * (component + 1) * p(component);
+ return_value[i] = 2 * (component + 1) * p[component];
return return_value;
}
for (int i = 0; i < 11; ++i)
{
Point<dim> testpoint;
- testpoint(0) = testcoord[i][0];
- testpoint(1) = testcoord[i][1];
+ testpoint[0] = testcoord[i][0];
+ testpoint[1] = testcoord[i][1];
if (dim == 3)
- testpoint(2) = testcoord[i][2];
+ testpoint[2] = testcoord[i][2];
bool res = cell->point_inside(testpoint);
deallog << testpoint << " inside " << res << std::endl;
// Now get the cell
const typename Triangulation<dim>::cell_iterator cell = triangulation.begin();
- cell->vertex(0)(0) = -1.;
+ cell->vertex(0)[0] = -1.;
// and test it.
double testcoord[14][3] = {{0.5, 0.5, 0.5},
for (int i = 0; i < 14; ++i)
{
Point<dim> testpoint;
- testpoint(0) = testcoord[i][0];
- testpoint(1) = testcoord[i][1];
+ testpoint[0] = testcoord[i][0];
+ testpoint[1] = testcoord[i][1];
if (dim == 3)
- testpoint(2) = testcoord[i][2];
+ testpoint[2] = testcoord[i][2];
bool res = cell->point_inside(testpoint);
deallog << testpoint << " \t inside " << res << " expected "
endc = triangulation.end();
for (; cell != endc; ++cell)
for (const unsigned int face : GeometryInfo<dim>::face_indices())
- if ((cell->face(face)->center()(0) == -1) ||
- (cell->face(face)->center()(1) == -1))
+ if ((cell->face(face)->center()[0] == -1) ||
+ (cell->face(face)->center()[1] == -1))
cell->face(face)->set_boundary_id(1);
}
else
const Point<dim> &p = points[i];
Point<dim> &beta = values[i];
- beta(0) = -p(1);
- beta(1) = p(0);
+ beta[0] = -p[1];
+ beta[1] = p[0];
beta /= std::sqrt(beta.square());
}
}
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.;
Solution<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
- 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;
}
RightHandSide<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
- 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;
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 4 * std::pow(p(i), 4);
+ return_value += 4 * std::pow(p[i], 4);
return return_value;
}
Vector<double> value2(1);
Point<2> point1;
- point1(0) = 2 * (1. * i / n_points + eps) - 1;
- point1(1) = -1.;
+ point1[0] = 2 * (1. * i / n_points + eps) - 1;
+ point1[1] = -1.;
Point<2> point2;
- point2(0) = 2 * (1. * i / n_points + eps) - 1;
- point2(1) = 1.;
+ point2[0] = 2 * (1. * i / n_points + eps) - 1;
+ point2[1] = 1.;
VectorTools::point_value(dof_handler, solution, point1, value1);
VectorTools::point_value(dof_handler, solution, point2, value2);
Vector<double> value2(1);
Point<2> point1;
- point1(0) = 2 * (1. * i / n_points + eps) - 1;
- point1(1) = -1.;
+ point1[0] = 2 * (1. * i / n_points + eps) - 1;
+ point1[1] = -1.;
Point<2> point2;
- point2(0) = 2 * (1. * i / n_points + eps) - 1;
- point2(1) = 1.;
+ point2[0] = 2 * (1. * i / n_points + eps) - 1;
+ point2[1] = 1.;
VectorTools::point_value(dof_handler, solution, point1, value1);
VectorTools::point_value(dof_handler, solution, point2, value2);
for (const unsigned int face : GeometryInfo<dim>::face_indices())
if (cell->face(face)->at_boundary())
for (unsigned int d = 0; d < dim; ++d)
- if ((std::fabs(cell->face(face)->center()(d) - (1)) < 1e-12))
+ if ((std::fabs(cell->face(face)->center()[d] - (1)) < 1e-12))
cell->face(face)->set_boundary_id(1);
}
for (const unsigned int face : GeometryInfo<dim>::face_indices())
if (cell->face(face)->at_boundary())
for (unsigned int d = 0; d < dim; ++d)
- if ((std::fabs(cell->face(face)->center()(d) - (1)) < 1e-12))
+ if ((std::fabs(cell->face(face)->center()[d] - (1)) < 1e-12))
cell->face(face)->set_boundary_id(1);
}
endc = triangulation.end();
for (; cell != endc; ++cell)
for (const unsigned int face : GeometryInfo<dim>::face_indices())
- if ((cell->face(face)->center()(0) == -1) ||
- (cell->face(face)->center()(1) == -1))
+ if ((cell->face(face)->center()[0] == -1) ||
+ (cell->face(face)->center()[1] == -1))
cell->face(face)->set_boundary_id(1);
}
else
Assert(dim >= 2, ExcNotImplemented());
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))
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
virtual double
value(const Point<dim> &p, const unsigned int component) const
{
- return p(component);
+ return p[component];
}
virtual void
vector_value(const Point<dim> &p, Vector<double> &values) const
{
for (unsigned int i = 0; i < dim; ++i)
- values(i) = p(i);
+ values(i) = p[i];
}
};
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
virtual double
value(const Point<dim> &p, const unsigned int component) const
{
- return p(0) * p(0) + 2.0 * p(0) * p(1);
+ return p[0] * p[0] + 2.0 * p[0] * p[1];
}
virtual void
value(const Point<dim> &p, const unsigned int component) const
{
double val = 0.0;
- if (std::abs(p(1) - 1.0) < 1e-5)
+ if (std::abs(p[1] - 1.0) < 1e-5)
val = 2.0;
deallog << "evaluate normal derivative at " << p << " with value " << val
// Turn cylinder such that y->x
for (unsigned int i = 0; i < 16; ++i)
{
- const double h = vertices[i](1);
- vertices[i](1) = -vertices[i](0);
- vertices[i](0) = h;
+ const double h = vertices[i][1];
+ vertices[i][1] = -vertices[i][0];
+ vertices[i][0] = h;
}
int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
const std::vector<Point<3>> &vertices = tria.get_vertices();
for (unsigned int i = 0; i < vertices.size(); ++i)
- if (vertices[i](2) > 1e-7)
+ if (vertices[i][2] > 1e-7)
std::cout << "Error!" << std::endl;
for (unsigned int i = 0; i < ncells; ++i)
for (unsigned int j = 0; j < spacedim; ++j)
- p.vertices[i](j) =
+ p.vertices[i][j] =
PatchInfo<dim>::vertices[i][j] + PatchInfo<dim>::offsets[c][j];
for (const unsigned int i : GeometryInfo<dim>::face_indices())
m[i].set_inhomogeneity(face->dof_index(0),
(face->center() *
(r_i / face->center().norm() -
- 1))(i));
+ 1))[i]);
}
else if (std::fabs(face->vertex(1).norm() - r_a) < eps)
for (unsigned int i = 0; i < 2; ++i)
Assert(dim >= 2, ExcNotImplemented());
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)
{
Assert(dim >= 2, ExcNotImplemented());
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)
{
for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
for (unsigned int d = 0; d < spacedim; ++d)
- patch.vertices[v](d) =
+ patch.vertices[v][d] =
p + cell_coordinates[d][v] + ((d >= dim) ? v : 0);
unsigned int n1 = (dim > 0) ? nsubp : 1;
Point<dim> midpoint;
for (unsigned int d = 0; d < dim; ++d)
- midpoint(d) = n_cells / 2.;
+ midpoint[d] = n_cells / 2.;
Functions::CutOffFunctionCinfty<dim> function(2., midpoint);
{
Point<dim> p = trapez.point(k);
if (dim >= 1)
- p(0) += i1;
+ p[0] += i1;
if (dim >= 2)
- p(1) += i2;
+ p[1] += i2;
if (dim >= 3)
- p(2) += i3;
+ p[2] += i3;
patch.vertices[k] = p;
}
std::vector<Point<dim>> points = trapezsub.get_points();
for (unsigned int k = 0; k < points.size(); ++k)
{
if (dim >= 1)
- points[k](0) += i1;
+ points[k][0] += i1;
if (dim >= 2)
- points[k](1) += i2;
+ points[k][1] += i2;
if (dim >= 3)
- points[k](2) += i3;
+ points[k][2] += i3;
}
std::vector<double> values(points.size());
function.value_list(points, values);
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:
vertices[GeometryInfo<dim>::vertices_per_cell + 2] =
vertices[GeometryInfo<dim>::vertices_per_cell + 3] = p2;
- vertices[GeometryInfo<dim>::vertices_per_cell + 1](0) = p2(0);
- vertices[GeometryInfo<dim>::vertices_per_cell + 2](0) = p1(0);
+ vertices[GeometryInfo<dim>::vertices_per_cell + 1][0] = p2[0];
+ vertices[GeometryInfo<dim>::vertices_per_cell + 2][0] = p1[0];
break;
case 3:
vertices[GeometryInfo<dim>::vertices_per_cell + 0] =
vertices[GeometryInfo<dim>::vertices_per_cell + 6] =
vertices[GeometryInfo<dim>::vertices_per_cell + 7] = p2;
- vertices[GeometryInfo<dim>::vertices_per_cell + 1](0) = p2(0);
- vertices[GeometryInfo<dim>::vertices_per_cell + 2](1) = p2(1);
- vertices[GeometryInfo<dim>::vertices_per_cell + 3](0) = p2(0);
- vertices[GeometryInfo<dim>::vertices_per_cell + 3](1) = p2(1);
+ vertices[GeometryInfo<dim>::vertices_per_cell + 1][0] = p2[0];
+ vertices[GeometryInfo<dim>::vertices_per_cell + 2][1] = p2[1];
+ vertices[GeometryInfo<dim>::vertices_per_cell + 3][0] = p2[0];
+ vertices[GeometryInfo<dim>::vertices_per_cell + 3][1] = p2[1];
- vertices[GeometryInfo<dim>::vertices_per_cell + 4](0) = p1(0);
- vertices[GeometryInfo<dim>::vertices_per_cell + 4](1) = p1(1);
- vertices[GeometryInfo<dim>::vertices_per_cell + 5](1) = p1(1);
- vertices[GeometryInfo<dim>::vertices_per_cell + 6](0) = p1(0);
+ vertices[GeometryInfo<dim>::vertices_per_cell + 4][0] = p1[0];
+ vertices[GeometryInfo<dim>::vertices_per_cell + 4][1] = p1[1];
+ vertices[GeometryInfo<dim>::vertices_per_cell + 5][1] = p1[1];
+ vertices[GeometryInfo<dim>::vertices_per_cell + 6][0] = p1[0];
break;
default:
if (cell->is_locally_owned())
for (const unsigned int face : GeometryInfo<dim>::face_indices())
{
- if (std::fabs(cell->face(face)->center()(0) - 0.0) < 1e-12)
+ if (std::fabs(cell->face(face)->center()[0] - 0.0) < 1e-12)
cell->face(face)->set_all_boundary_ids(1);
- if (std::fabs(cell->face(face)->center()(0) - 1.0) < 1e-12)
+ if (std::fabs(cell->face(face)->center()[0] - 1.0) < 1e-12)
cell->face(face)->set_all_boundary_ids(2);
- if (std::fabs(cell->face(face)->center()(1) - 0.0) < 1e-12)
+ if (std::fabs(cell->face(face)->center()[1] - 0.0) < 1e-12)
cell->face(face)->set_all_boundary_ids(3);
- if (std::fabs(cell->face(face)->center()(1) - 1.0) < 1e-12)
+ if (std::fabs(cell->face(face)->center()[1] - 1.0) < 1e-12)
cell->face(face)->set_all_boundary_ids(4);
}
Point<dim> middle = FlatManifold<dim>::get_new_point_on_line(line);
for (int i = 0; i < dim; ++i)
- middle(i) -= .5;
+ middle[i] -= .5;
middle *=
std::sqrt(static_cast<double>(dim)) / (std::sqrt(middle.square()) * 2);
for (int i = 0; i < dim; ++i)
- middle(i) += .5;
+ middle[i] += .5;
return middle;
}
Point<dim> middle = FlatManifold<dim>::get_new_point_on_quad(quad);
for (int i = 0; i < dim; ++i)
- middle(i) -= .5;
+ middle[i] -= .5;
middle *=
std::sqrt(static_cast<double>(dim)) / (std::sqrt(middle.square()) * 2);
for (int i = 0; i < dim; ++i)
- middle(i) += .5;
+ middle[i] += .5;
return middle;
}
// 0 or 1, then the z-values of all
// vertices of the line is like that
if (dim >= 3)
- if (((middle(2) == 0) || (middle(2) == 1))
+ if (((middle[2] == 0) || (middle[2] == 1))
// find out, if the line is in the
// interior of the top or bottom face
// of the domain, or at the edge.
return middle;
- double x = middle(0), y = middle(1);
+ double x = middle[0], y = middle[1];
if (y < x)
if (y < 1 - x)
- middle(1) = 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 0.04 * std::sin(6 * numbers::PI * middle[0]);
else
- middle(0) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[1]);
else if (y < 1 - x)
- middle(0) = 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 0.04 * std::sin(6 * numbers::PI * middle[1]);
else
- middle(1) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[0]);
return middle;
}
// z-value of the midpoint is either
// 0 or 1, then the z-values of all
// vertices of the quad is like that
- if (dim == 3 && ((middle(dim - 1) == 0) || (middle(dim - 1) == 1)))
+ if (dim == 3 && ((middle[dim - 1] == 0) || (middle[dim - 1] == 1)))
return middle;
- double x = middle(0), y = middle(1);
+ double x = middle[0], y = middle[1];
if (y < x)
if (y < 1 - x)
- middle(1) = 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 0.04 * std::sin(6 * numbers::PI * middle[0]);
else
- middle(0) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[1]);
else if (y < 1 - x)
- middle(0) = 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 0.04 * std::sin(6 * numbers::PI * middle[1]);
else
- middle(1) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[0]);
return middle;
}
for (unsigned int ref = 0; ref < 2; ++ref)
{
for (auto &cell : tria.active_cell_iterators())
- if (cell->is_locally_owned() && cell->center()(0) < .5 &&
- cell->center()(1) < .5)
+ if (cell->is_locally_owned() && cell->center()[0] < .5 &&
+ cell->center()[1] < .5)
cell->set_refine_flag();
tria.execute_coarsening_and_refinement();
}
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<3>::face_indices())
{
- if (cell_1->face(j)->center()(2) > 2.9)
+ if (cell_1->face(j)->center()[2] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(2) < -2.9)
+ if (cell_2->face(j)->center()[2] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
{
for (const unsigned int j : GeometryInfo<dim>::face_indices())
{
- if (cell->face(j)->center()(dim == 2 ? 1 : 2) > 2.9)
+ if (cell->face(j)->center()[dim == 2 ? 1 : 2] > 2.9)
face_1 = cell->face(j);
- if (cell->face(j)->center()(dim == 2 ? 1 : 2) < -2.9)
+ if (cell->face(j)->center()[dim == 2 ? 1 : 2] < -2.9)
face_2 = cell->face(j);
}
}
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<3>::face_indices())
{
- if (cell_1->face(j)->center()(2) > 2.9)
+ if (cell_1->face(j)->center()[2] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(2) < -2.9)
+ if (cell_2->face(j)->center()[2] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
{
for (const unsigned int j : GeometryInfo<dim>::face_indices())
{
- if (cell->face(j)->center()(dim == 2 ? 1 : 2) > 2.9)
+ if (cell->face(j)->center()[dim == 2 ? 1 : 2] > 2.9)
face_1 = cell->face(j);
- if (cell->face(j)->center()(dim == 2 ? 1 : 2) < -2.9)
+ if (cell->face(j)->center()[dim == 2 ? 1 : 2] < -2.9)
face_2 = cell->face(j);
}
}
triangulation.begin_active();
cell != triangulation.end();
++cell)
- if (cell->center()(0) > 0.49)
+ if (cell->center()[0] > 0.49)
cell->set_refine_flag();
triangulation.prepare_coarsening_and_refinement();
for (unsigned int ref = 0; ref < 2; ++ref)
{
for (auto &cell : tria.active_cell_iterators())
- if (cell->is_locally_owned() && cell->center()(0) < .5 &&
- cell->center()(1) < .5)
+ if (cell->is_locally_owned() && cell->center()[0] < .5 &&
+ cell->center()[1] < .5)
cell->set_refine_flag();
tria.execute_coarsening_and_refinement();
}
bool
pred_d(const typename DoFHandler<dim>::active_cell_iterator &cell)
{
- return (cell->center()(0) < 0.49 && cell->center()(1) < 0.49);
+ return (cell->center()[0] < 0.49 && cell->center()[1] < 0.49);
}
template <int dim>
bool
pred_d(const typename DoFHandler<dim>::active_cell_iterator &cell)
{
- return (cell->center()(0) < 0.5 && cell->center()(1) < 0.5);
+ return (cell->center()[0] < 0.5 && cell->center()[1] < 0.5);
}
bool
pred_left(const typename DoFHandler<dim>::active_cell_iterator &cell)
{
- return (cell->center()(0) < 0.49);
+ return (cell->center()[0] < 0.49);
}
template <int dim>
bool
pred_right(const typename DoFHandler<dim>::active_cell_iterator &cell)
{
- return (cell->center()(0) > 0.51);
+ return (cell->center()[0] > 0.51);
}
bool
pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
{
- return (cell->center()(0) < 0.49 && cell->center()(1) < 0.49) ||
- (cell->center()(0) > 0.49 && cell->center()(1) > 0.49);
+ return (cell->center()[0] < 0.49 && cell->center()[1] < 0.49) ||
+ (cell->center()[0] > 0.49 && cell->center()[1] > 0.49);
}
bool
pred_d(const typename DoFHandler<dim>::active_cell_iterator &cell)
{
- return (cell->center()(0) < 0.49);
+ return (cell->center()[0] < 0.49);
}
template <int dim>
bool
pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
{
- return (cell->center()(0) < 0.49 && cell->center()(1) < 0.49) ||
- (cell->center()(0) > 0.49 && cell->center()(1) > 0.49);
+ return (cell->center()[0] < 0.49 && cell->center()[1] < 0.49) ||
+ (cell->center()[0] > 0.49 && cell->center()[1] > 0.49);
}
value(const Point<2> &p, unsigned int c = 0) const override
{
(void)c;
- return p(0) * p(1);
+ return p[0] * p[1];
}
double left = -1.0, right = 1.0;
Point<2> left_point, right_point;
for (unsigned int i = 0; i < 2; ++i)
- left_point(i) = left, right_point(i) = right;
+ left_point[i] = left, right_point[i] = right;
GridGenerator::subdivided_hyper_cube(tr, 4, left, right);
FE_Hermite<2> herm(2 * regularity + 1);
Point<2>
stretch_coordinates(const Point<2> p)
{
- return Point<2>(2 * p(0), p(1));
+ return Point<2>(2 * p[0], p[1]);
}
Point<2>
tilt_coordinates(const Point<2> p)
{
- return Point<2>(p(0) + p(1), p(1));
+ return Point<2>(p[0] + p[1], p[1]);
}
TestDef1<dim>::value(const Point<dim> &p, const unsigned int component) const
{
Point<2> center;
- center(0) = 0.5;
- center(1) = 0.5;
+ center[0] = 0.5;
+ center[1] = 0.5;
double rad = p.distance(center),
- phi_p = atan2(p(0) - center(0), p(1) - center(1));
+ phi_p = atan2(p[0] - center[0], p[1] - center[1]);
if (component == 0)
return rad * (sin(phi + phi_p) - sin(phi_p));
double
TestDef2<dim>::value(const Point<dim> &p, const unsigned int component) const
{
- double x = p(0), y = p(1);
+ double x = p[0], y = p[1];
if (component == 0)
return scale * x;
double
TestDef3<dim>::value(const Point<dim> &p, const unsigned int component) const
{
- double y = p(1);
+ double y = p[1];
if (component == 0)
return scale * y;
double
TestPoly<dim>::value(const Point<dim> &p, const unsigned int component) const
{
- double x = p(0), y = p(1);
+ double x = p[0], y = p[1];
if (component == 0)
return polys[0].value(x) + polys[1].value(y);
switch (component)
{
case 0:
- val = cos(numbers::PI * p(0)) * sin(numbers::PI * p(1)) + bc_constant;
+ val = cos(numbers::PI * p[0]) * sin(numbers::PI * p[1]) + bc_constant;
case 1:
- val = -sin(numbers::PI * p(0)) * cos(numbers::PI * p(1)) + bc_constant;
+ val = -sin(numbers::PI * p[0]) * cos(numbers::PI * p[1]) + bc_constant;
}
return val;
}
Vector<double> &result) const
{
Assert(dim >= 2, ExcNotImplemented());
- result(0) = cos(numbers::PI * p(0)) * sin(numbers::PI * p(1)) + bc_constant;
- result(1) = -sin(numbers::PI * p(0)) * cos(numbers::PI * p(1)) + bc_constant;
+ result(0) = cos(numbers::PI * p[0]) * sin(numbers::PI * p[1]) + bc_constant;
+ result(1) = -sin(numbers::PI * p[0]) * cos(numbers::PI * p[1]) + bc_constant;
}
template <int dim>
void
{
case 0:
values[i] =
- cos(numbers::PI * p(0)) * sin(numbers::PI * p(1)) + bc_constant;
+ cos(numbers::PI * p[0]) * sin(numbers::PI * p[1]) + bc_constant;
case 1:
values[i] =
- -sin(numbers::PI * p(0)) * cos(numbers::PI * p(1)) + bc_constant;
+ -sin(numbers::PI * p[0]) * cos(numbers::PI * p[1]) + bc_constant;
}
}
}
{
const Point<dim> &p = points[i];
values[i](0) =
- cos(numbers::PI * p(0)) * sin(numbers::PI * p(1)) + bc_constant;
+ cos(numbers::PI * p[0]) * sin(numbers::PI * p[1]) + bc_constant;
values[i](1) =
- -sin(numbers::PI * p(0)) * cos(numbers::PI * p(1)) + bc_constant;
+ -sin(numbers::PI * p[0]) * cos(numbers::PI * p[1]) + bc_constant;
}
}
// END EXACT SOLUTION MEMBERS
Assert(dim >= 2, ExcNotImplemented());
// 2D solution
- values(0) = (2 * numbers::PI * numbers::PI + 1) * cos(numbers::PI * p(0)) *
- sin(numbers::PI * p(1)) +
+ values(0) = (2 * numbers::PI * numbers::PI + 1) * cos(numbers::PI * p[0]) *
+ sin(numbers::PI * p[1]) +
bc_constant;
- values(1) = -(2 * numbers::PI * numbers::PI + 1) * sin(numbers::PI * p(0)) *
- cos(numbers::PI * p(1)) +
+ values(1) = -(2 * numbers::PI * numbers::PI + 1) * sin(numbers::PI * p[0]) *
+ cos(numbers::PI * p[1]) +
bc_constant;
}
template <int dim>
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
{
Point<dim> &point = cell->vertex(i);
- if (std::abs(point(dim - 1) - 1.0) < 1e-5)
- point(dim - 1) += 0.15;
+ if (std::abs(point[dim - 1] - 1.0) < 1e-5)
+ point[dim - 1] += 0.15;
}
FE_ABF<dim> fe(1);
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
{
Point<dim> &point = cell->vertex(i);
- if (std::abs(point(dim - 1) - 1.0) < 1e-5)
- point(dim - 1) += 0.15;
+ if (std::abs(point[dim - 1] - 1.0) < 1e-5)
+ point[dim - 1] += 0.15;
}
FE_BernardiRaugel<dim> fe(1);
const Point<dim> &p = points[i];
/* quadratic: */
- value_list[i](0) = p(0) * p(0);
- value_list[i](1) = p(1) * p(1);
- value_list[i](2) = p(2) * p(2);
+ value_list[i](0) = p[0] * p[0];
+ value_list[i](1) = p[1] * p[1];
+ value_list[i](2) = p[2] * p[2];
}
}
// Additional functions to create Neumann conditions, zero in this case.
switch (component)
{
case 0:
- val = cos(PI * p(0)) * sin(PI * p(1));
+ val = cos(PI * p[0]) * sin(PI * p[1]);
break;
case 1:
- val = -sin(PI * p(0)) * cos(PI * p(1));
+ val = -sin(PI * p[0]) * cos(PI * p[1]);
break;
}
return val;
switch (component)
{
case 0:
- val = -sin(PI * p(0)) * cos(PI * p(1)) * cos(PI * p(2));
+ val = -sin(PI * p[0]) * cos(PI * p[1]) * cos(PI * p[2]);
break;
case 1:
- val = -cos(PI * p(0)) * sin(PI * p(1)) * cos(PI * p(2));
+ val = -cos(PI * p[0]) * sin(PI * p[1]) * cos(PI * p[2]);
break;
case 2:
- val = 2 * cos(PI * p(0)) * cos(PI * p(1)) * sin(PI * p(2));
+ val = 2 * cos(PI * p[0]) * cos(PI * p[1]) * sin(PI * p[2]);
break;
}
return val;
{
const double PI = numbers::PI;
Tensor<1, 3> val;
- double x = p(0), y = p(1), z = p(2);
+ double x = p[0], y = p[1], z = p[2];
switch (component)
{
{
double return_value = 1.;
for (unsigned int i = 0; i < dim; ++i)
- return_value *= (1 - p(i) * p(i));
- return_value *= std::pow(p(m_direction), m_degree - 1);
+ return_value *= (1 - p[i] * p[i]);
+ return_value *= std::pow(p[m_direction], m_degree - 1);
return return_value;
}
grad[d] = 1.;
// compute grad(\prod_{i=1}^d (1-x_i^2))(p)
for (unsigned j = 0; j < dim; ++j)
- grad[d] *= (d == j ? -2 * p(j) : (1 - p(j) * p(j)));
+ grad[d] *= (d == j ? -2 * p[j] : (1 - p[j] * p[j]));
// and multiply with x_i^{r-1}
- grad[d] *= std::pow(p(m_direction), m_degree - 1);
+ grad[d] *= std::pow(p[m_direction], m_degree - 1);
}
if (m_degree >= 2)
// add \prod_{i=1}^d (1-x_i^2))(p)
double value = 1.;
for (unsigned int j = 0; j < dim; ++j)
- value *= (1 - p(j) * p(j));
+ value *= (1 - p[j] * p[j]);
// and multiply with grad(x_i^{r-1})
grad[m_direction] +=
- value * (m_degree - 1) * std::pow(p(m_direction), m_degree - 2);
+ value * (m_degree - 1) * std::pow(p[m_direction], m_degree - 2);
}
return grad;
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
{
Point<dim> &point = cell->vertex(i);
- if (std::abs(point(dim - 1) - 1.0) < 1e-5)
- point(dim - 1) += 0.15;
+ if (std::abs(point[dim - 1] - 1.0) < 1e-5)
+ point[dim - 1] += 0.15;
}
FE_RaviartThomas<dim> fe(2);
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
{
Point<dim> &point = cell->vertex(i);
- if (std::abs(point(dim - 1) - 1.0) < 1e-5)
- point(dim - 1) += 0.15;
+ if (std::abs(point[dim - 1] - 1.0) < 1e-5)
+ point[dim - 1] += 0.15;
}
FE_RaviartThomas<dim> fe(2);
TestFunction::value(const Point<2> &p, const unsigned int component) const
{
Assert(component == 0, ExcInternalError());
- return std::sin(pi * p(0)) * std::cos(pi * p(1));
+ return std::sin(pi * p[0]) * std::cos(pi * p[1]);
}
// Jacobian contains many nonzero entries
Point<dim> quad_p;
for (int d = 0; d < dim; ++d)
- quad_p(d) = 0.42 + 0.11 * d;
+ quad_p[d] = 0.42 + 0.11 * d;
Quadrature<dim> quad(quad_p);
Point<dim - 1> f_quad_p;
for (int d = 0; d < dim - 1; ++d)
- f_quad_p(d) = 0.42 + 0.11 * d;
+ f_quad_p[d] = 0.42 + 0.11 * d;
Quadrature<dim - 1> f_quad(f_quad_p);
switch (component)
{
case 0:
- val = pow(p(0), 3);
+ val = pow(p[0], 3);
break;
case 1:
- val = pow(p(1), 2) * p(0);
+ val = pow(p[1], 2) * p[0];
break;
}
return val;
switch (component)
{
case 0:
- val = pow(p(0), 3);
+ val = pow(p[0], 3);
break;
case 1:
- val = pow(p(1), 2) * p(0);
+ val = pow(p[1], 2) * p[0];
break;
case 2:
- val = p(2) * p(1) * p(0);
+ val = p[2] * p[1] * p[0];
break;
}
return val;
double result = 1.;
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int k = 0; k < degree; ++k)
- result *= p(d) + c;
+ result *= p[d] + c;
return result;
}
double result = 1.;
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int k = 0; k < degree; ++k)
- result *= p(d) + c;
+ result *= p[d] + c;
values[i] = result;
}
}
double result = 1.;
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int k = 0; k < degree; ++k)
- result *= p(d);
+ result *= p[d];
values[i](c) = result;
}
}
// contains many nonzero entries
Point<dim> quad_p;
for (int d = 0; d < dim; ++d)
- quad_p(d) = 0.42 + 0.11 * d;
+ quad_p[d] = 0.42 + 0.11 * d;
Quadrature<dim> quad(quad_p);
{
// Jacobian contains many nonzero entries
Point<dim - 1> quad_p;
for (int d = 0; d < dim - 1; ++d)
- quad_p(d) = 0.42 + 0.11 * d;
+ quad_p[d] = 0.42 + 0.11 * d;
Quadrature<dim - 1> quad(quad_p);
{
// Jacobian contains many nonzero entries
Point<dim - 1> quad_p;
for (int d = 0; d < dim - 1; ++d)
- quad_p(d) = 0.42 + 0.11 * d;
+ quad_p[d] = 0.42 + 0.11 * d;
Quadrature<dim - 1> quad(quad_p);
{
// Jacobian contains many nonzero entries
Point<dim - 1> quad_p;
for (int d = 0; d < dim - 1; ++d)
- quad_p(d) = 0.42 + 0.11 * d;
+ quad_p[d] = 0.42 + 0.11 * d;
Quadrature<dim - 1> quad(quad_p);
{
Point<2>
stretch_coordinates(const Point<2> p)
{
- return Point<2>(2 * p(0), p(1));
+ return Point<2>(2 * p[0], p[1]);
}
Point<2>
tilt_coordinates(const Point<2> p)
{
- return Point<2>(p(0) + p(1), p(1));
+ return Point<2>(p[0] + p[1], p[1]);
}
const Point<dim> &p = points[i];
/* quadratic: */
- value_list[i](0) = p(0) * p(0);
- value_list[i](1) = p(1) * p(1);
- value_list[i](2) = p(2) * p(2);
+ value_list[i](0) = p[0] * p[0];
+ value_list[i](1) = p[1] * p[1];
+ value_list[i](2) = p[2] * p[2];
}
}
// Additional functions to create Neumann conditions, zero in this case.
stretch(const Point<dim> &p)
{
Point<dim> q = p;
- q(dim - 1) *= 2.;
+ q[dim - 1] *= 2.;
return q;
}
Point<dim> q = p;
if (dim >= 2)
{
- q(0) = 2. * p(0) + p(1);
- q(1) = p(0) + 3. * p(1);
+ q[0] = 2. * p[0] + p[1];
+ q[1] = p[0] + 3. * p[1];
}
return q;
Point<dim> q = p;
if (dim >= 2)
{
- q(0) = 16. * p(0) + 4. * p(1) - 10. * p(0) * p(1);
- q(1) = 6. * p(1) + 4. * p(0) * p(1);
+ q[0] = 16. * p[0] + 4. * p[1] - 10. * p[0] * p[1];
+ q[1] = 6. * p[1] + 4. * p[0] * p[1];
}
return q;
stretch(const Point<dim> &p)
{
Point<dim> q = p;
- q(dim - 1) *= 2.;
+ q[dim - 1] *= 2.;
return q;
}
Point<dim> q = p;
if (dim >= 2)
{
- q(0) = 2. * p(0) + p(1);
- q(1) = p(0) + 3. * p(1);
+ q[0] = 2. * p[0] + p[1];
+ q[1] = p[0] + 3. * p[1];
}
return q;
Point<dim> q = p;
if (dim >= 2)
{
- q(0) = 16. * p(0) + 4. * p(1) - 10. * p(0) * p(1);
- q(1) = 6. * p(1) + 4. * p(0) * p(1);
+ q[0] = 16. * p[0] + 4. * p[1] - 10. * p[0] * p[1];
+ q[1] = 6. * p[1] + 4. * p[0] * p[1];
}
return q;
Point<2>
stretch_coordinates(const Point<2> p)
{
- return Point<2>(2 * p(0), p(1));
+ return Point<2>(2 * p[0], p[1]);
}
Point<2>
tilt_coordinates(const Point<2> p)
{
- return Point<2>(p(0) + p(1), p(1));
+ return Point<2>(p[0] + p[1], p[1]);
}
TestDef1<dim>::value(const Point<dim> &p, const unsigned int component) const
{
Point<2> center;
- center(0) = 0.5;
- center(1) = 0.5;
+ center[0] = 0.5;
+ center[1] = 0.5;
double rad = p.distance(center),
- phi_p = atan2(p(0) - center(0), p(1) - center(1));
+ phi_p = atan2(p[0] - center[0], p[1] - center[1]);
if (component == 0)
return rad * (sin(phi + phi_p) - sin(phi_p));
double
TestDef2<dim>::value(const Point<dim> &p, const unsigned int component) const
{
- double x = p(0), y = p(1);
+ double x = p[0], y = p[1];
if (component == 0)
return scale * x;
double
TestDef3<dim>::value(const Point<dim> &p, const unsigned int component) const
{
- double y = p(1);
+ double y = p[1];
if (component == 0)
return scale * y;
double
TestPoly<dim>::value(const Point<dim> &p, const unsigned int component) const
{
- double x = p(0), y = p(1);
+ double x = p[0], y = p[1];
// Ugly hack, but should do the job ...
if (component == 0)
Point<2>
stretch_coordinates(const Point<2> p)
{
- return Point<2>(2 * p(0), p(1));
+ return Point<2>(2 * p[0], p[1]);
}
Point<2>
tilt_coordinates(const Point<2> p)
{
- return Point<2>(p(0) + p(1), p(1));
+ return Point<2>(p[0] + p[1], p[1]);
}
case 1:
return p;
case 2:
- return Point<dim>(p(0) * (1 + p(1)), p(1) * (1 + p(0)));
+ return Point<dim>(p[0] * (1 + p[1]), p[1] * (1 + p[0]));
case 3:
- return Point<dim>(p(0) * (1 + p(1)) * (1 + p(2)),
- p(1) * (1 + p(0)) * (1 + p(2)),
- p(2) * (1 + p(0)) * (1 + p(1)));
+ return Point<dim>(p[0] * (1 + p[1]) * (1 + p[2]),
+ p[1] * (1 + p[0]) * (1 + p[2]),
+ p[2] * (1 + p[0]) * (1 + p[1]));
default:
Assert(false, ExcNotImplemented());
return Point<dim>();
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.;
Assert(dim >= 2, ExcNotImplemented());
Point<dim> wind_field;
- wind_field(0) = -p(1);
- wind_field(1) = p(0);
+ wind_field[0] = -p[1];
+ wind_field[1] = p[0];
wind_field /= wind_field.norm();
return wind_field;
Solution<2>::value(const Point<2> &p, const unsigned int component) const
{
using numbers::PI;
- const double x = p(0);
- const double y = p(1);
+ const double x = p[0];
+ const double y = p[1];
// zero on BD's
if (component == 0)
return PI * sin(PI * x) * sin(PI * x) * sin(2.0 * PI * y);
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. * PI * sin(PI * x) * sin(PI * x) * sin(2.0 * PI * y) *
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);
+ double x = p[0];
+ double y = p[1];
double nu = 1.0;
// RHS for 0 BD's
Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1));
using numbers::PI;
- double x = p(0);
- double y = p(1);
- double z = p(2);
+ double x = p[0];
+ double y = p[1];
+ double z = p[2];
if (component == 0)
return 4. * PI * PI * PI *
auto endc = tria.end();
for (; cell != endc; ++cell)
for (const unsigned int face_number : GeometryInfo<dim>::face_indices())
- if (std::fabs(cell->face(face_number)->center()(0) - left) < 1e-12)
+ if (std::fabs(cell->face(face_number)->center()[0] - left) < 1e-12)
cell->face(face_number)->set_all_boundary_ids(1);
- else if (std::fabs(cell->face(face_number)->center()(0) - right) <
+ else if (std::fabs(cell->face(face_number)->center()[0] - right) <
1e-12)
cell->face(face_number)->set_all_boundary_ids(2);
else if (dim >= 2 &&
- std::fabs(cell->face(face_number)->center()(1) - left) < 1e-12)
+ std::fabs(cell->face(face_number)->center()[1] - left) < 1e-12)
cell->face(face_number)->set_all_boundary_ids(3);
- else if (dim >= 2 && std::fabs(cell->face(face_number)->center()(1) -
+ else if (dim >= 2 && std::fabs(cell->face(face_number)->center()[1] -
right) < 1e-12)
cell->face(face_number)->set_all_boundary_ids(4);
else if (dim >= 3 &&
- std::fabs(cell->face(face_number)->center()(2) - left) < 1e-12)
+ std::fabs(cell->face(face_number)->center()[2] - left) < 1e-12)
cell->face(face_number)->set_all_boundary_ids(5);
- else if (dim >= 3 && std::fabs(cell->face(face_number)->center()(2) -
+ else if (dim >= 3 && std::fabs(cell->face(face_number)->center()[2] -
right) < 1e-12)
cell->face(face_number)->set_all_boundary_ids(6);
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<2> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
break;
}
case 2:
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<2> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
Point<2> &v3 = tria.begin_active()->vertex(3);
- v3(0) = 4.;
+ v3[0] = 4.;
break;
}
case 3:
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<3> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
break;
}
case 3:
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<2> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
Point<2> &v2 = tria.begin_active()->vertex(3);
- v2(0) = 5.;
- v2(1) = 4.;
+ v2[0] = 5.;
+ v2[1] = 4.;
break;
}
default:
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<3> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
break;
}
default:
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < 2; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << 2), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < 2; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << 2), ExcInternalError());
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < 2; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << 2), ExcInternalError());
for (unsigned int d = 0; d < dim; ++d)
{
- left(d) = -1.;
- right(d) = d + 2;
+ left[d] = -1.;
+ right[d] = d + 2;
subdivisions[d] = n * (d + 3);
}
Triangulation<dim> tria;
// Turn cylinder such that y->x
for (unsigned int i = 0; i < 16; ++i)
{
- const double h = vertices[i](1);
- vertices[i](1) = -vertices[i](0);
- vertices[i](0) = h;
+ const double h = vertices[i][1];
+ vertices[i][1] = -vertices[i][0];
+ vertices[i][0] = h;
}
int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
dim_2(std::ostream &os)
{
std::vector<Point<2>> vertices(4);
- vertices[0](0) = -1.;
- vertices[0](1) = -1.;
- vertices[1](0) = 1.;
- vertices[1](1) = -1.5;
- vertices[2](0) = 1.5;
- vertices[2](1) = 1.5;
- vertices[3](0) = 2.;
- vertices[3](1) = 0.5;
+ vertices[0][0] = -1.;
+ vertices[0][1] = -1.;
+ vertices[1][0] = 1.;
+ vertices[1][1] = -1.5;
+ vertices[2][0] = 1.5;
+ vertices[2][1] = 1.5;
+ vertices[3][0] = 2.;
+ vertices[3][1] = 0.5;
Triangulation<2> tria;
GridGenerator::general_cell<2>(tria, vertices);
dim_3(std::ostream &os)
{
std::vector<Point<3>> vertices(8);
- vertices[0](0) = -1.;
- vertices[0](1) = -1.;
- vertices[0](2) = -1.;
- vertices[1](0) = 1.;
- vertices[1](1) = -1.5;
- vertices[1](2) = -1.5;
- vertices[2](0) = 2.;
- vertices[2](1) = 1.5;
- vertices[2](2) = -2.;
- vertices[3](0) = 2.5;
- vertices[3](1) = 0.5;
- vertices[3](2) = -3.;
- vertices[4](0) = -1.;
- vertices[4](1) = -1.;
- vertices[4](2) = 1.;
- vertices[5](0) = 1.;
- vertices[5](1) = -1.5;
- vertices[5](2) = 1.5;
- vertices[6](0) = 2.;
- vertices[6](1) = 1.5;
- vertices[6](2) = 2.;
- vertices[7](0) = 2.;
- vertices[7](1) = 0.5;
- vertices[7](2) = 3.;
+ vertices[0][0] = -1.;
+ vertices[0][1] = -1.;
+ vertices[0][2] = -1.;
+ vertices[1][0] = 1.;
+ vertices[1][1] = -1.5;
+ vertices[1][2] = -1.5;
+ vertices[2][0] = 2.;
+ vertices[2][1] = 1.5;
+ vertices[2][2] = -2.;
+ vertices[3][0] = 2.5;
+ vertices[3][1] = 0.5;
+ vertices[3][2] = -3.;
+ vertices[4][0] = -1.;
+ vertices[4][1] = -1.;
+ vertices[4][2] = 1.;
+ vertices[5][0] = 1.;
+ vertices[5][1] = -1.5;
+ vertices[5][2] = 1.5;
+ vertices[6][0] = 2.;
+ vertices[6][1] = 1.5;
+ vertices[6][2] = 2.;
+ vertices[7][0] = 2.;
+ vertices[7][1] = 0.5;
+ vertices[7][2] = 3.;
Triangulation<3> tria;
GridGenerator::general_cell<3>(tria, vertices);
dim_2_3(std::ostream &os)
{
std::vector<Point<3>> vertices(4);
- vertices[0](0) = -1.;
- vertices[0](1) = -1.;
- vertices[0](0) = 0.;
+ vertices[0][0] = -1.;
+ vertices[0][1] = -1.;
+ vertices[0][0] = 0.;
- vertices[1](0) = 1.;
- vertices[1](1) = -1.5;
- vertices[1](2) = 0.;
+ vertices[1][0] = 1.;
+ vertices[1][1] = -1.5;
+ vertices[1][2] = 0.;
- vertices[2](0) = 1.5;
- vertices[2](1) = 1.5;
- vertices[2](2) = 0.;
+ vertices[2][0] = 1.5;
+ vertices[2][1] = 1.5;
+ vertices[2][2] = 0.;
- vertices[3](0) = 2.;
- vertices[3](1) = 0.5;
- vertices[3](2) = 0.;
+ vertices[3][0] = 2.;
+ vertices[3][1] = 0.5;
+ vertices[3][2] = 0.;
Triangulation<2, 3> tria;
GridGenerator::general_cell<2, 3>(tria, vertices);
dim_1_3(std::ostream &os)
{
std::vector<Point<3>> vertices(2);
- vertices[0](0) = -1.;
- vertices[0](1) = -1.;
- vertices[0](2) = -1.;
- vertices[1](0) = 1.;
- vertices[1](1) = -1.5;
- vertices[1](2) = -1.5;
+ vertices[0][0] = -1.;
+ vertices[0][1] = -1.;
+ vertices[0][2] = -1.;
+ vertices[1][0] = 1.;
+ vertices[1][1] = -1.5;
+ vertices[1][2] = -1.5;
Triangulation<1, 3> tria;
GridGenerator::general_cell<1, 3>(tria, vertices);
dim_1_2(std::ostream &os)
{
std::vector<Point<2>> vertices(2);
- vertices[0](0) = -1.;
- vertices[0](1) = -1.;
+ vertices[0][0] = -1.;
+ vertices[0][1] = -1.;
- vertices[1](0) = 1.;
- vertices[1](1) = -1.5;
+ vertices[1][0] = 1.;
+ vertices[1][1] = -1.5;
Triangulation<1, 2> tria;
GridGenerator::general_cell<1, 2>(tria, vertices);
Triangulation<d> tr;
std::vector<Point<d>> vertices(d + 1);
- vertices[1](0) = 0.5;
- vertices[1](1) = .85;
- vertices[2](0) = -0.5;
- vertices[2](1) = .85;
+ vertices[1][0] = 0.5;
+ vertices[1][1] = .85;
+ vertices[2][0] = -0.5;
+ vertices[2][1] = .85;
GridGenerator::simplex(tr, vertices);
GridOut gout;
Triangulation<d> tr;
std::vector<Point<d>> vertices(d + 1);
- vertices[0](0) = 1.;
- vertices[0](1) = 0.;
- vertices[0](2) = .7;
- vertices[1](0) = -1.;
- vertices[1](1) = 0.;
- vertices[1](2) = .7;
- vertices[2](0) = 0.;
- vertices[2](1) = 1.;
- vertices[2](2) = -.7;
- vertices[3](0) = 0.;
- vertices[3](1) = -1.;
- vertices[3](2) = -.7;
+ vertices[0][0] = 1.;
+ vertices[0][1] = 0.;
+ vertices[0][2] = .7;
+ vertices[1][0] = -1.;
+ vertices[1][1] = 0.;
+ vertices[1][2] = .7;
+ vertices[2][0] = 0.;
+ vertices[2][1] = 1.;
+ vertices[2][2] = -.7;
+ vertices[3][0] = 0.;
+ vertices[3][1] = -1.;
+ vertices[3][2] = -.7;
GridGenerator::simplex(tr, vertices);
GridOut gout;
test(bool second_case = false)
{
std::vector<Point<dim>> vertices(GeometryInfo<dim>::vertices_per_cell);
- vertices[1](1) = 1;
- vertices[2](0) = 1;
- vertices[2](1) = 1;
- vertices[3](0) = 1;
+ vertices[1][1] = 1;
+ vertices[2][0] = 1;
+ vertices[2][1] = 1;
+ vertices[3][0] = 1;
if (dim == 3)
{
for (unsigned int i = 4; i < GeometryInfo<dim>::vertices_per_cell; ++i)
- vertices[i](2) = -1;
- vertices[5](1) = 1;
- vertices[6](0) = 1;
- vertices[6](1) = 1;
- vertices[7](0) = 1;
+ vertices[i][2] = -1;
+ vertices[5][1] = 1;
+ vertices[6][0] = 1;
+ vertices[6][1] = 1;
+ vertices[7][0] = 1;
}
std::vector<CellData<dim>> cells(1);
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
std::swap(cells[0].vertices[1], cells[0].vertices[3]);
std::swap(cells[0].vertices[5], cells[0].vertices[7]);
for (unsigned int i = 4; i < GeometryInfo<dim>::vertices_per_cell; ++i)
- vertices[i](2) = 1;
+ vertices[i][2] = 1;
}
SubCellData subcelldata;
{
const static int dim = 2;
std::vector<Point<dim>> vertices(6);
- vertices[1](1) = 1;
- vertices[2](0) = 1;
- vertices[3](0) = 1;
- vertices[3](1) = 1;
- vertices[4](0) = 2;
- vertices[5](0) = 2;
- vertices[5](1) = 1;
+ vertices[1][1] = 1;
+ vertices[2][0] = 1;
+ vertices[3][0] = 1;
+ vertices[3][1] = 1;
+ vertices[4][0] = 2;
+ vertices[5][0] = 2;
+ vertices[5][1] = 1;
std::vector<CellData<dim>> cells(2);
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
FlatManifold<dim>::get_new_point(surrounding_points, weights);
for (int i = 0; i < dim; ++i)
- middle(i) -= .5;
+ middle[i] -= .5;
middle *=
std::sqrt(static_cast<double>(dim)) / (std::sqrt(middle.square()) * 2);
for (int i = 0; i < dim; ++i)
- middle(i) += .5;
+ middle[i] += .5;
return middle;
}
// 0 or 1, then the z-values of all
// vertices of the line is like that
if (dim >= 3)
- if (((middle(2) == 0) || (middle(2) == 1))
+ if (((middle[2] == 0) || (middle[2] == 1))
// find out, if the line is in the
// interior of the top or bottom face
// of the domain, or at the edge.
return middle;
- double x = middle(0), y = middle(1);
+ double x = middle[0], y = middle[1];
if (y < x)
if (y < 1 - x)
- middle(1) = 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 0.04 * std::sin(6 * numbers::PI * middle[0]);
else
- middle(0) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[1]);
else if (y < 1 - x)
- middle(0) = 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 0.04 * std::sin(6 * numbers::PI * middle[1]);
else
- middle(1) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[0]);
return middle;
}
// z-value of the midpoint is either
// 0 or 1, then the z-values of all
// vertices of the quad is like that
- if (dim == 3 && ((middle(dim - 1) == 0) || (middle(dim - 1) == 1)))
+ if (dim == 3 && ((middle[dim - 1] == 0) || (middle[dim - 1] == 1)))
return middle;
- double x = middle(0), y = middle(1);
+ double x = middle[0], y = middle[1];
if (y < x)
if (y < 1 - x)
- middle(1) = 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 0.04 * std::sin(6 * numbers::PI * middle[0]);
else
- middle(0) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[1]);
else if (y < 1 - x)
- middle(0) = 0.04 * std::sin(6 * numbers::PI * middle(1));
+ middle[0] = 0.04 * std::sin(6 * numbers::PI * middle[1]);
else
- middle(1) = 1 + 0.04 * std::sin(6 * numbers::PI * middle(0));
+ middle[1] = 1 + 0.04 * std::sin(6 * numbers::PI * middle[0]);
return middle;
}
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<3>::face_indices())
{
- if (cell_1->face(j)->center()(2) > 2.9)
+ if (cell_1->face(j)->center()[2] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(2) < -2.9)
+ if (cell_2->face(j)->center()[2] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<2>::face_indices())
{
- if (cell_1->face(j)->center()(1) > 2.9)
+ if (cell_1->face(j)->center()[1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(1) < -2.9)
+ if (cell_2->face(j)->center()[1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<3>::face_indices())
{
- if (cell_1->face(j)->center()(2) > 2.9)
+ if (cell_1->face(j)->center()[2] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(2) < -2.9)
+ if (cell_2->face(j)->center()[2] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
Point<2>
trans_func(const Point<2> &p)
{
- Point<2> r(p(0) + p(1) * p(1), p(1));
+ Point<2> r(p[0] + p[1] * p[1], p[1]);
return r;
}
template <>
Point<3>
trans_func(const Point<3> &p)
{
- Point<3> r(p(0) + p(1) * p(1), p(1), p(2));
+ Point<3> r(p[0] + p[1] * p[1], p[1], p[2]);
return r;
}
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<2> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
Point<2> &v2 = tria.begin_active()->vertex(3);
- v2(0) = 5.;
- v2(1) = 4.;
+ v2[0] = 5.;
+ v2[1] = 4.;
// exact_areas.push_back(7.);
break;
}
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<3> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
break;
}
default:
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<2> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
break;
}
case 2:
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<2> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
Point<2> &v3 = tria.begin_active()->vertex(3);
- v3(0) = 4.;
+ v3[0] = 4.;
break;
}
default:
{
GridGenerator::hyper_cube(tria, 1., 3.);
Point<3> &v0 = tria.begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
break;
}
default:
Point<3>
distort_planar(Point<3> p)
{
- if (p(1) > 0.5 && p(2) > 0.5)
+ if (p[1] > 0.5 && p[2] > 0.5)
{
- p(1) += 1;
+ p[1] += 1;
}
return p;
}
Point<3>
distort_twisted(Point<3> p)
{
- if (p(2) > 0.5 && ((p(0) > 0.5) ^ (p(1) > 0.5)))
+ if (p[2] > 0.5 && ((p[0] > 0.5) ^ (p[1] > 0.5)))
{
- p(2) += 1;
+ p[2] += 1;
}
return p;
}
Assert(dim == 2 || dim == 3, ExcNotImplemented());
std::vector<Point<dim>> vertices(GeometryInfo<dim>::vertices_per_cell);
- vertices[0](0) = 0;
- vertices[0](1) = 0;
- vertices[1](0) = 2;
- vertices[1](1) = 1;
- vertices[2](0) = 3;
- vertices[2](1) = 3;
- vertices[3](0) = 0;
- vertices[3](1) = 1;
+ vertices[0][0] = 0;
+ vertices[0][1] = 0;
+ vertices[1][0] = 2;
+ vertices[1][1] = 1;
+ vertices[2][0] = 3;
+ vertices[2][1] = 3;
+ vertices[3][0] = 0;
+ vertices[3][1] = 1;
if (dim == 3)
{
// for the new numbering
// for the old numbering
for (unsigned int i = 0; i < 4; ++i)
{
- std::swap(vertices[i](1), vertices[i](2));
+ std::swap(vertices[i][1], vertices[i][2]);
vertices[i + 4] = vertices[i];
- vertices[i + 4](1) = 1;
+ vertices[i + 4][1] = 1;
}
}
{
unsigned int subdomain = 0;
for (unsigned int d = 0; d < dim; ++d)
- if (cell->center()(d) > 0)
+ if (cell->center()[d] > 0)
subdomain |= (1 << d);
AssertThrow(subdomain < (1 << dim), ExcInternalError());
double
ExactSolution<dim>::value(const Point<dim> &p, const unsigned int) const
{
- return p(0) * p(0);
+ return p[0] * p[0];
}
template <int dim>
virtual double
value(const Point<dim> &p, const unsigned int) const
{
- return p(0);
+ return p[0];
}
};
/*Make a square*/
Point<dim> point_1, point_2;
- point_1(0) = 0;
- point_1(1) = 0;
- point_2(0) = 1;
- point_2(1) = 1;
+ point_1[0] = 0;
+ point_1[1] = 0;
+ point_2[0] = 1;
+ point_2[1] = 1;
GridGenerator::hyper_rectangle(triangulation, point_1, point_2);
Triangulation<dim> triangulation_temp;
- point_1(0) = 1;
- point_2(0) = 2;
+ point_1[0] = 1;
+ point_2[0] = 2;
GridGenerator::hyper_rectangle(triangulation_temp, point_1, point_2);
/*glue squares together*/
GridGenerator::merge_triangulations(triangulation_temp,
// Define a rectangular shape
for (unsigned int d = 0; d < dim; ++d)
{
- p1(d) = 0;
- p2(d) = (d == 0) ? 2.0 : 1.0;
+ p1[d] = 0;
+ p2[d] = (d == 0) ? 2.0 : 1.0;
sub_div.push_back((d == 0) ? 2 : 1);
}
GridGenerator::subdivided_hyper_rectangle(tria, sub_div, p1, p2, true);
const Point<dim> &p = points[i];
Point<dim> &beta = values[i];
- beta(0) = -p(1);
- beta(1) = p(0);
+ beta[0] = -p[1];
+ beta[1] = p[0];
beta /= std::sqrt(beta.square());
}
}
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.;
Solution<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
- 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;
}
RightHandSide<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
- 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;
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 4 * std::pow(p(i), 4);
+ return_value += 4 * std::pow(p[i], 4);
return return_value;
}
endc = triangulation.end();
for (; cell != endc; ++cell)
for (const unsigned int face : GeometryInfo<dim>::face_indices())
- if ((cell->face(face)->center()(0) == -1) ||
- (cell->face(face)->center()(1) == -1))
+ if ((cell->face(face)->center()[0] == -1) ||
+ (cell->face(face)->center()[1] == -1))
cell->face(face)->set_boundary_id(1);
}
else
Assert(dim >= 2, ExcNotImplemented());
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))
++cell)
{
real_cell_center = cell->center();
- if (real_cell_center(0) < 0.5 && real_cell_center(1) < 0.5)
+ if (real_cell_center[0] < 0.5 && real_cell_center[1] < 0.5)
cell->set_material_id(0); // solid
else
cell->set_material_id(1); // fluid
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 4 * std::pow(p(i), 4);
+ return_value += 4 * std::pow(p[i], 4);
return return_value;
}
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 2 * std::pow(p(i), 2);
+ return_value += 2 * std::pow(p[i], 2);
return return_value;
}
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 2 * std::pow(p(i), 2);
+ return_value += 2 * std::pow(p[i], 2);
return return_value;
}
Point<dim> p;
- p(0) = 0.45;
- p(1) = 0.16;
+ p[0] = 0.45;
+ p[1] = 0.16;
for (unsigned int step = 0; step < 10; ++step)
{
std::cout << "step " << step << std::endl;
tria_ptr.push_back(tria);
GridGenerator::hyper_cube(*tria, 1., 3.);
Point<2> &v0 = tria->begin_active()->vertex(0);
- v0(0) = 0.;
+ v0[0] = 0.;
Point<2> &v3 = tria->begin_active()->vertex(3);
- v3(0) = 5.;
- v3(1) = 4.;
+ v3[0] = 5.;
+ v3[1] = 4.;
exact_areas.push_back(7.);
show[1][0] = 1;
}
GridGenerator::hyper_cube(*tria, 1., 5.);
Point<2> &v2 = tria->begin_active()->vertex(2);
Point<2> &v3 = tria->begin_active()->vertex(3);
- v2(0) = 1.;
- v2(1) = 3.;
- v3(0) = 3.;
- v3(1) = 3.;
+ v2[0] = 1.;
+ v2[1] = 3.;
+ v3[0] = 3.;
+ v3[1] = 3.;
tria->set_manifold(1, *boundary1);
tria->set_manifold(2, *boundary2);
tria->begin_active()->face(1)->set_manifold_id(1);
tria_ptr.push_back(tria);
Point<2> p0;
Point<2> p1;
- p0(0) = 1.;
- p0(1) = 2.5;
- p1(0) = 2.;
- p1(1) = 4.;
+ p0[0] = 1.;
+ p0[1] = 2.5;
+ p1[0] = 2.;
+ p1[1] = 4.;
GridGenerator::hyper_rectangle(*tria, p0, p1);
exact_areas.push_back(1.5);
show[3][4] = 1;
GridGenerator::hyper_cube(*tria, 0., 1.);
Point<2> &v2 = tria->begin_active()->vertex(2);
Point<2> &v3 = tria->begin_active()->vertex(3);
- v2(0) = 2.;
- v2(1) = 1.;
- v3(0) = 0.5;
- v3(1) = 1.5;
+ v2[0] = 2.;
+ v2[1] = 1.;
+ v3[0] = 0.5;
+ v3[1] = 1.5;
tria->set_manifold(1, *boundary1);
tria->begin_active()->face(1)->set_manifold_id(1);
exact_areas.push_back(0.);
tria_ptr.push_back(tria);
GridGenerator::hyper_cube(*tria, 1., 3.);
Point<3> &v = tria->begin()->vertex(7);
- v(0) = 5.;
- v(1) = 4.;
- v(2) = 4.5;
+ v[0] = 5.;
+ v[1] = 4.;
+ v[2] = 4.5;
exact_areas.push_back(12.5);
}
tria_ptr.push_back(tria);
Point<3> p0;
Point<3> p1;
- p0(0) = 1.;
- p0(1) = 2.5;
- p0(2) = 3.;
- p1(0) = 2.;
- p1(1) = 4.;
- p1(2) = 6.;
+ p0[0] = 1.;
+ p0[1] = 2.5;
+ p0[2] = 3.;
+ p1[0] = 2.;
+ p1[1] = 4.;
+ p1[2] = 6.;
GridGenerator::hyper_rectangle(*tria, p0, p1);
exact_areas.push_back(4.5);
show[3][4] = 1;
// shift one point of the cell
// somehow
if (dim > 1)
- tr.begin_active()->vertex(dim == 2 ? 3 : 5)(dim - 1) += 1. / std::sqrt(2.);
+ tr.begin_active()->vertex(dim == 2 ? 3 : 5)[dim - 1] += 1. / std::sqrt(2.);
DoFHandler<dim> dof(tr);
dof.distribute_dofs(fe);
ImposedDisplacement<2>::vector_value(const Point<2> &p,
Vector<double> &value) const
{
- double radius = 1 + (sqrt(5) - 1) * p(0);
- double angle = 0.5 * numbers::PI * (1 - p(1));
- value(0) = radius * sin(angle) - p(0);
- value(1) = radius * cos(angle) - p(1);
+ double radius = 1 + (sqrt(5) - 1) * p[0];
+ double angle = 0.5 * numbers::PI * (1 - p[1]);
+ value(0) = radius * sin(angle) - p[0];
+ value(1) = radius * cos(angle) - p[1];
}
double
VectorFunction<3>::value(const Point<3> &p, const unsigned int component) const
{
- 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]);
}
template <>
double
VectorFunction<2>::value(const Point<2> &p, const unsigned int component) const
{
- return (1 - p(0) * p(0)) * (1 - p(1) * p(1));
+ return (1 - p[0] * p[0]) * (1 - p[1] * p[1]);
}
template <int dim>
for (unsigned int j = 0; j < 19; ++j)
{
/// all points are inside
- points[19 * i + j](0) = -0.7 + (i + 1) * .07;
- points[19 * i + j](1) = -0.7 + (j + 1) * .07;
+ points[19 * i + j][0] = -0.7 + (i + 1) * .07;
+ points[19 * i + j][1] = -0.7 + (j + 1) * .07;
}
points[95] = p;
fe_function.value_list(points, m); // <<<< this fails at point[95] but only if
double sinval = deformation;
for (unsigned int d = 0; d < dim; ++d)
sinval *= std::sin(frequency * dealii::numbers::PI *
- (chart_point(d) - left) / (right - left));
+ (chart_point[d] - left) / (right - left));
dealii::Point<dim> space_point;
for (unsigned int d = 0; d < dim; ++d)
- space_point(d) = chart_point(d) + sinval;
+ space_point[d] = chart_point[d] + sinval;
return space_point;
}
dealii::Point<dim> x = space_point;
dealii::Point<dim> one;
for (unsigned int d = 0; d < dim; ++d)
- one(d) = 1.;
+ one[d] = 1.;
// Newton iteration to solve the nonlinear equation given by the point
dealii::Tensor<1, dim> sinvals;
for (unsigned int d = 0; d < dim; ++d)
- sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x(d) - left) /
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x[d] - left) /
(right - left));
double sinval = deformation;
double sinval_der = deformation * frequency / (right - left) *
dealii::numbers::PI *
std::cos(frequency * dealii::numbers::PI *
- (x(d) - left) / (right - left));
+ (x[d] - left) / (right - left));
for (unsigned int e = 0; e < dim; ++e)
if (e != d)
sinval_der *= sinvals[e];
for (unsigned int d = 0; d < dim; ++d)
sinvals[d] = std::sin(frequency * dealii::numbers::PI *
- (x(d) - left) / (right - left));
+ (x[d] - left) / (right - left));
sinval = deformation;
for (unsigned int d = 0; d < dim; ++d)
double sinval = deformation;
for (unsigned int d = 0; d < dim; ++d)
sinval *= std::sin(frequency * dealii::numbers::PI *
- (chart_point(d) - left) / (right - left));
+ (chart_point[d] - left) / (right - left));
dealii::Point<dim> space_point;
for (unsigned int d = 0; d < dim; ++d)
- space_point(d) = chart_point(d) + sinval;
+ space_point[d] = chart_point[d] + sinval;
return space_point;
}
dealii::Point<dim> x = space_point;
dealii::Point<dim> one;
for (unsigned int d = 0; d < dim; ++d)
- one(d) = 1.;
+ one[d] = 1.;
// Newton iteration to solve the nonlinear equation given by the point
dealii::Tensor<1, dim> sinvals;
for (unsigned int d = 0; d < dim; ++d)
- sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x(d) - left) /
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x[d] - left) /
(right - left));
double sinval = deformation;
double sinval_der = deformation * frequency / (right - left) *
dealii::numbers::PI *
std::cos(frequency * dealii::numbers::PI *
- (x(d) - left) / (right - left));
+ (x[d] - left) / (right - left));
for (unsigned int e = 0; e < dim; ++e)
if (e != d)
sinval_der *= sinvals[e];
for (unsigned int d = 0; d < dim; ++d)
sinvals[d] = std::sin(frequency * dealii::numbers::PI *
- (x(d) - left) / (right - left));
+ (x[d] - left) / (right - left));
sinval = deformation;
for (unsigned int d = 0; d < dim; ++d)
double sinval = deformation;
for (unsigned int d = 0; d < dim; ++d)
sinval *= std::sin(frequency * dealii::numbers::PI *
- (chart_point(d) - left) / (right - left));
+ (chart_point[d] - left) / (right - left));
dealii::Point<dim> space_point;
for (unsigned int d = 0; d < dim; ++d)
- space_point(d) = chart_point(d) + sinval;
+ space_point[d] = chart_point[d] + sinval;
return space_point;
}
dealii::Point<dim> x = space_point;
dealii::Point<dim> one;
for (unsigned int d = 0; d < dim; ++d)
- one(d) = 1.;
+ one[d] = 1.;
// Newton iteration to solve the nonlinear equation given by the point
dealii::Tensor<1, dim> sinvals;
for (unsigned int d = 0; d < dim; ++d)
- sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x(d) - left) /
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x[d] - left) /
(right - left));
double sinval = deformation;
double sinval_der = deformation * frequency / (right - left) *
dealii::numbers::PI *
std::cos(frequency * dealii::numbers::PI *
- (x(d) - left) / (right - left));
+ (x[d] - left) / (right - left));
for (unsigned int e = 0; e < dim; ++e)
if (e != d)
sinval_der *= sinvals[e];
for (unsigned int d = 0; d < dim; ++d)
sinvals[d] = std::sin(frequency * dealii::numbers::PI *
- (x(d) - left) / (right - left));
+ (x[d] - left) / (right - left));
sinval = deformation;
for (unsigned int d = 0; d < dim; ++d)
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
{
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
}
return_value[i] = 2.0 * viscosity_function.value(p);
}
{
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
- p(d) = pressure.quadrature_point(q)(d)[i];
+ p[d] = pressure.quadrature_point(q)[d][i];
return_value[i] = 1.0 / viscosity_function.value(p);
}
one_over_viscosity(cell, q) = return_value;
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
{
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
}
return_value[i] = 2.0 * viscosity_function.value(p);
}
{
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
Vector<double> rhs_temp(dim + 1);
right_hand_side.vector_value(p, rhs_temp);
auto point = gpu_data->get_quadrature_point(cell, q_point);
dst[pos] =
- dim == 2 ? point(0) + point(1) : point(0) + point(1) + point(2);
+ dim == 2 ? point[0] + point[1] : point[0] + point[1] + point[2];
});
}
const unsigned int pos =
gpu_data_host.local_q_point_id(cell_id, n_q_points_per_cell, i);
auto p = gpu_data_host.get_quadrature_point(cell_id, i);
- const double p_val = dim == 2 ? p(0) + p(1) : p(0) + p(1) + p(2);
+ const double p_val = dim == 2 ? p[0] + p[1] : p[0] + p[1] + p[2];
AssertThrow(std::abs(coef[pos] - p_val) < 1e-12,
ExcInternalError());
}
if (!vertex_moved[global_vertex_no] &&
locally_owned_vertices[global_vertex_no])
{
- cell->vertex(vertex_no)(0) += 1.e-1;
+ cell->vertex(vertex_no)[0] += 1.e-1;
vertex_moved[global_vertex_no] = true;
}
}
for (const unsigned int vertex_no : GeometryInfo<dim>::vertex_indices())
{
Point<dim> point = cell->vertex(vertex_no);
- point(0) -= 1.e-1;
+ point[0] -= 1.e-1;
non_artificial_vertices_new[cell->vertex_index(vertex_no)] = point;
}
const Point<dim> &p = cell->center();
bool all_positive = true;
for (unsigned int d = 0; d < dim; ++d)
- if (p(d) <= 0.)
+ if (p[d] <= 0.)
all_positive = false;
if (all_positive)
cell->set_refine_flag();
{
double return_value = 0.;
for (unsigned int d = 0; d < dim; ++d)
- return_value += std::pow(std::abs(.5 - p(d)), degree);
+ return_value += std::pow(std::abs(.5 - p[d]), degree);
return return_value;
}
for (unsigned int n_ref = 0; n_ref < 2; ++n_ref)
{
for (const auto &cell : triangulation.active_cell_iterators())
- if (cell->is_locally_owned() && cell->center()(1) < 0.5)
+ if (cell->is_locally_owned() && cell->center()[1] < 0.5)
cell->set_refine_flag();
triangulation.prepare_coarsening_and_refinement();
for (const auto &cell : dof_handler.active_cell_iterators() |
IteratorFilters::LocallyOwnedCell())
{
- if (cell->center()(0) < 0.5)
+ if (cell->center()[0] < 0.5)
cell->set_active_fe_index(0);
else
cell->set_active_fe_index(1);
double
value(const Point<dim> &p, const unsigned int) const
{
- return 1 + p(0) * p(0);
+ return 1 + p[0] * p[0];
}
};
for (const auto &cell : tria.cell_iterators())
for (const unsigned int face_index : GeometryInfo<dim>::face_indices())
{
- if (std::abs(cell->face(face_index)->center()(face_index / 2)) < 1e-12)
+ if (std::abs(cell->face(face_index)->center()[face_index / 2]) < 1e-12)
cell->face(face_index)->set_all_boundary_ids(face_index);
- if (std::abs(cell->face(face_index)->center()(face_index / 2) - 1.) <
+ if (std::abs(cell->face(face_index)->center()[face_index / 2] - 1.) <
1e-12)
cell->face(face_index)->set_all_boundary_ids(face_index);
}
TemperatureInitialValues<dim>::value(const Point<dim> &p,
const unsigned int) const
{
- return p(0) * T1 + p(1) * (T0 - T1); // simple
+ return p[0] * T1 + p[1] * (T0 - T1); // simple
}
TemperatureInitialValues<dim>::value(const Point<dim> &p,
const unsigned int) const
{
- return p(0) * T1 + p(1) * (T0 - T1); // simple
+ return p[0] * T1 + p[1] * (T0 - T1); // simple
}
Vector<PetscScalar> value2(1);
Point<2> point1;
- point1(0) = 1. * i / n_points + eps;
- point1(1) = 0.;
+ point1[0] = 1. * i / n_points + eps;
+ point1[1] = 0.;
Point<2> point2;
- point2(0) = 1. * i / n_points + eps;
- point2(1) = 1.;
+ point2[0] = 1. * i / n_points + eps;
+ point2[1] = 1.;
get_point_value(point1, 0, value1);
get_point_value(point2, 0, value2);
Vector<PetscScalar> value4(1);
Point<3> point1;
- point1(0) = 1. * i / n_points + eps;
- point1(1) = 1. * j / n_points + eps;
- point1(2) = 0;
+ point1[0] = 1. * i / n_points + eps;
+ point1[1] = 1. * j / n_points + eps;
+ point1[2] = 0;
Point<3> point2;
- point2(0) = 1. * i / n_points + eps;
- point2(1) = 1. * j / n_points + eps;
- point2(2) = 1.;
+ point2[0] = 1. * i / n_points + eps;
+ point2[1] = 1. * j / n_points + eps;
+ point2[2] = 1.;
Point<3> point3;
- point3(0) = 1. * i / n_points + eps;
- point3(1) = 0.;
- point3(2) = 1. * j / n_points + eps;
+ point3[0] = 1. * i / n_points + eps;
+ point3[1] = 0.;
+ point3[2] = 1. * j / n_points + eps;
Point<3> point4;
- point4(0) = 1. * i / n_points + eps;
- point4(1) = 1.;
- point4(2) = 1. * j / n_points + eps;
+ point4[0] = 1. * i / n_points + eps;
+ point4[1] = 1.;
+ point4[2] = 1. * j / n_points + eps;
get_point_value(point1, 0, value1);
get_point_value(point2, 0, value2);
Point<dim> p1;
Point<dim> p2;
for (unsigned int i = 0; i < dim; ++i)
- p2(i) = 1.0;
+ p2[i] = 1.0;
GridGenerator::subdivided_hyper_rectangle(
triangulation, reps, p1, p2, true);
Vector<double> value2(3);
Point<2> point1;
- point1(0) = 0;
- point1(1) = .5 * (1. + 1. * i / n_points + eps);
+ point1[0] = 0;
+ point1[1] = .5 * (1. + 1. * i / n_points + eps);
Point<2> point2;
- point2(0) = .5 * (1. + 1. * i / n_points + eps);
- point2(1) = 0.;
+ point2[0] = .5 * (1. + 1. * i / n_points + eps);
+ point2[1] = 0.;
get_point_value(point1, 0, value1);
get_point_value(point2, 0, value2);
fe_face_values.get_quadrature_points();
for (unsigned int i = 0; i < tmp_points.size(); ++i)
for (unsigned int c = 0; c < dim; ++c)
- local_quad_points_first.push_back(tmp_points[i](c));
+ local_quad_points_first.push_back(tmp_points[i][c]);
}
else if (face->boundary_id() == 4)
{
fe_face_values.get_quadrature_points();
for (unsigned int i = 0; i < tmp_points.size(); ++i)
for (unsigned int c = 0; c < dim; ++c)
- local_quad_points_second.push_back(tmp_points[i](c));
+ local_quad_points_second.push_back(tmp_points[i][c]);
}
}
}
for (unsigned int c = 0; c < dim; ++c)
{
vector_point_1(c) = global_quad_points_first[i + c];
- point_1(c) = vector_point_1(c);
+ point_1[c] = vector_point_1(c);
}
Vector<double> vector_point_2(dim);
rot_matrix.Tvmult(vector_point_2, vector_point_1);
for (unsigned int c = 0; c < dim; ++c)
- point_2(c) = vector_point_2(c) + offset[c];
+ point_2[c] = vector_point_2(c) + offset[c];
get_point_value(point_1, 0, value_1);
get_point_value(point_2, 0, value_2);
for (unsigned int c = 0; c < dim; ++c)
{
vector_point_1(c) = global_quad_points_second[i + c];
- point_1(c) = vector_point_1(c);
+ point_1[c] = vector_point_1(c);
}
Vector<double> vector_point_2(dim);
for (unsigned int c = 0; c < dim; ++c)
vector_point_1(c) -= offset[c];
rot_matrix.vmult(vector_point_2, vector_point_1);
for (unsigned int c = 0; c < dim; ++c)
- point_2(c) = vector_point_2(c);
+ point_2[c] = vector_point_2(c);
get_point_value(point_1, 0, value_1);
get_point_value(point_2, 0, value_2);
// Look for the two outermost faces:
for (const unsigned int j : GeometryInfo<dim>::face_indices())
{
- if (cell_1->face(j)->center()(dim - 1) > 2.9)
+ if (cell_1->face(j)->center()[dim - 1] > 2.9)
face_1 = cell_1->face(j);
- if (cell_2->face(j)->center()(dim - 1) < -2.9)
+ if (cell_2->face(j)->center()[dim - 1] < -2.9)
face_2 = cell_2->face(j);
}
face_1->set_boundary_id(42);
// now refine and check if the neighboring faces are correctly found
typename Triangulation<dim>::active_cell_iterator cell;
for (cell = triangulation.begin_active(); cell != triangulation.end(); ++cell)
- if (cell->is_locally_owned() && cell->center()(dim - 1) > 0)
+ if (cell->is_locally_owned() && cell->center()[dim - 1] > 0)
cell->set_refine_flag();
triangulation.execute_coarsening_and_refinement();
const unsigned int face_no_2 = it->second.first.second;
const Point<dim> face_center_1 = cell_1->face(face_no_1)->center();
const Point<dim> face_center_2 = cell_2->face(face_no_2)->center();
- Assert(std::min(std::abs(face_center_1(dim - 1) - 3.),
- std::abs(face_center_1(dim - 1) + 3.)) < 1.e-8,
+ Assert(std::min(std::abs(face_center_1[dim - 1] - 3.),
+ std::abs(face_center_1[dim - 1] + 3.)) < 1.e-8,
ExcInternalError());
- Assert(std::min(std::abs(face_center_2(dim - 1) - 3.),
- std::abs(face_center_2(dim - 1) + 3.)) < 1.e-8,
+ Assert(std::min(std::abs(face_center_2[dim - 1] - 3.),
+ std::abs(face_center_2[dim - 1] + 3.)) < 1.e-8,
ExcInternalError());
if (cell_1->level() == cell_2->level())
for (unsigned int c = 0; c < dim - 1; ++c)
- if (std::abs(face_center_1(c) - face_center_2(c)) > 1.e-8)
+ if (std::abs(face_center_1[c] - face_center_2[c]) > 1.e-8)
{
std::cout << "face_center_1: " << face_center_1 << std::endl;
std::cout << "face_center_2: " << face_center_2 << std::endl;
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
{
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
}
return_value[i] = 2.0 * viscosity_function.value(p);
}
{
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
- p(d) = pressure.quadrature_point(q)(d)[i];
+ p[d] = pressure.quadrature_point(q)[d][i];
return_value[i] = 1.0 / viscosity_function.value(p);
}
one_over_viscosity(cell, q) = return_value;
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
{
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
}
return_value[i] = 2.0 * viscosity_function.value(p);
}
{
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
Vector<double> rhs_temp(dim + 1);
right_hand_side.vector_value(p, rhs_temp);
mg_dof_handler.distribute_dofs(fe);
mg_dof_handler.distribute_mg_dofs();
Point<dim> a;
- a(0) = 1;
+ a[0] = 1;
for (unsigned int level = 0; level < tria.n_levels(); ++level)
DoFRenumbering::downstream(mg_dof_handler, level, a);
}
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
++cell)
{
const Point<dim> p = cell->center();
- bool positive = p(0) > 0;
+ bool positive = p[0] > 0;
if (positive)
{
cell->set_refine_flag();
Point<dim> center;
for (int d = 0; d < dim; ++d)
- center(d) = .5;
+ center[d] = .5;
// For each coordinate direction set up a plane through the center.
for (int plane_direction = 0; plane_direction < dim; ++plane_direction)
Tensor<1, dim> plane_normal;
for (unsigned int i = 1; i < dim; i++)
{
- point_on_zero_contour(i) = 0.25;
+ point_on_zero_contour[i] = 0.25;
plane_normal[i] = 1;
}
Point<dim> center;
for (int d = 0; d < dim; ++d)
{
- center(d) = .5;
+ center[d] = .5;
}
for (int direction = 0; direction < dim; ++direction)
{
for (int i = 0; i < dim; ++i)
{
normal[i] = 1;
- center(i) += epsilon;
+ center[i] += epsilon;
}
const Functions::SignedDistance::Plane<dim> level_set(center, normal);
: Functions::ConstantFunction<dim>(1)
{
for (int d = 0; d < dim; ++d)
- unit_box_center(d) = .5;
+ unit_box_center[d] = .5;
}
SymmetricTensor<2, dim>
double
value(const Point<1> &point, const unsigned int component = 0) const override
{
- return C * std::pow(point(0) - x_0, 2) + y_0;
+ return C * std::pow(point[0] - x_0, 2) + y_0;
};
Tensor<1, 1>
const unsigned int component = 0) const override
{
Tensor<1, 1> grad;
- grad[0] = 2 * C * (point(0) - x_0);
+ grad[0] = 2 * C * (point[0] - x_0);
return grad;
};
Point<dim - 1> point;
for (int i = 0; i < dim - 1; ++i)
{
- point(i) = 10 * (i + 1);
+ point[i] = 10 * (i + 1);
}
// Both points in the 1D-quadrature have weight 1/2 so
// this should also be the weight of the points in the final
Point<dim> p(tria.begin_active()->center());
for (unsigned int i = 0; i < dim; ++i)
- orientation(i) = i;
+ orientation[i] = i;
Vector<double> vector(dof.n_dofs());
Point<dim> p(tria.begin_active()->center());
for (unsigned int i = 0; i < dim; ++i)
- orientation(i) = i;
+ orientation[i] = i;
Vector<double> vector(dof.n_dofs());
virtual double
value(const Point<spacedim> &p, const unsigned int component = 0) const
{
- return p(0) * p(0);
+ return p[0] * p[0];
}
};
virtual double
value(const Point<spacedim> &p, const unsigned int component = 0) const
{
- return p(0) * p(0);
+ return p[0] * p[0];
}
};
continue;
double radius = cell->face(f)->center().norm() - center.norm();
- if (std::fabs(cell->face(f)->center()(2) -
- sqrt(3.) * cell->face(f)->center()(0)) <
+ if (std::fabs(cell->face(f)->center()[2] -
+ sqrt(3.) * cell->face(f)->center()[0]) <
eps) // z = sqrt(3)x set boundary 2
{
cell->face(f)->set_boundary_id(2);
cell->face(f)->line(j)->vertex(1).norm()) > eps)
cell->face(f)->line(j)->set_boundary_id(2);
}
- else if (std::fabs(cell->face(f)->center()(2) +
- sqrt(3.) * cell->face(f)->center()(0)) <
+ else if (std::fabs(cell->face(f)->center()[2] +
+ sqrt(3.) * cell->face(f)->center()[0]) <
eps) // z = -sqrt(3)x set boundary 3
{
cell->face(f)->set_boundary_id(3);
cell->face(f)->line(j)->vertex(1).norm()) > eps)
cell->face(f)->line(j)->set_boundary_id(3);
}
- else if (std::fabs(cell->face(f)->center()(2) -
- sqrt(3.) * cell->face(f)->center()(1)) <
+ else if (std::fabs(cell->face(f)->center()[2] -
+ sqrt(3.) * cell->face(f)->center()[1]) <
eps) // z = sqrt(3)y set boundary 4
{
cell->face(f)->set_boundary_id(4);
cell->face(f)->line(j)->vertex(1).norm()) > eps)
cell->face(f)->line(j)->set_boundary_id(4);
}
- else if (std::fabs(cell->face(f)->center()(2) +
- sqrt(3.) * cell->face(f)->center()(1)) <
+ else if (std::fabs(cell->face(f)->center()[2] +
+ sqrt(3.) * cell->face(f)->center()[1]) <
eps) // z = -sqrt(3)y set boundary 5
{
cell->face(f)->set_boundary_id(5);
{
if (cell->face(face)->at_boundary())
{
- if ((std::fabs(cell->face(face)->center()(0)) < 0.1) &&
- (std::fabs(cell->face(face)->center()(dim - 1)) <
+ if ((std::fabs(cell->face(face)->center()[0]) < 0.1) &&
+ (std::fabs(cell->face(face)->center()[dim - 1]) <
1e-12))
{
cell->face(face)->set_boundary_id(1);
}
- if ((std::fabs(cell->face(face)->center()(0)) < 1e-12) &&
- (std::fabs(cell->face(face)->center()(dim - 1)) < 0.1))
+ if ((std::fabs(cell->face(face)->center()[0]) < 1e-12) &&
+ (std::fabs(cell->face(face)->center()[dim - 1]) < 0.1))
{
cell->face(face)->set_boundary_id(1);
}
- if ((std::fabs(1.0 - cell->face(face)->center()(0)) <
+ if ((std::fabs(1.0 - cell->face(face)->center()[0]) <
0.1) &&
- (std::fabs(1.0 - cell->face(face)->center()(dim - 1)) <
+ (std::fabs(1.0 - cell->face(face)->center()[dim - 1]) <
1e-12))
{
cell->face(face)->set_boundary_id(2);
}
- if ((std::fabs(1.0 - cell->face(face)->center()(0)) <
+ if ((std::fabs(1.0 - cell->face(face)->center()[0]) <
1e-12) &&
- (std::fabs(1.0 - cell->face(face)->center()(dim - 1)) <
+ (std::fabs(1.0 - cell->face(face)->center()[dim - 1]) <
0.1))
{
cell->face(face)->set_boundary_id(2);
// no normal flux boundary
- if ((std::fabs(cell->face(face)->center()(0)) >= 0.1 &&
- std::fabs(cell->face(face)->center()(0)) <= 1.0) &&
- (std::fabs(cell->face(face)->center()(dim - 1)) <
+ if ((std::fabs(cell->face(face)->center()[0]) >= 0.1 &&
+ std::fabs(cell->face(face)->center()[0]) <= 1.0) &&
+ (std::fabs(cell->face(face)->center()[dim - 1]) <
1e-12))
{
cell->face(face)->set_boundary_id(3);
}
- if ((std::fabs(cell->face(face)->center()(0)) >= 0.0 &&
- std::fabs(cell->face(face)->center()(0)) <= 0.9) &&
- (std::fabs(1.0 - cell->face(face)->center()(dim - 1)) <
+ if ((std::fabs(cell->face(face)->center()[0]) >= 0.0 &&
+ std::fabs(cell->face(face)->center()[0]) <= 0.9) &&
+ (std::fabs(1.0 - cell->face(face)->center()[dim - 1]) <
1e-12))
{
cell->face(face)->set_boundary_id(5);
}
- if ((std::fabs(1.0 - cell->face(face)->center()(0)) <
+ if ((std::fabs(1.0 - cell->face(face)->center()[0]) <
1e-12) &&
- (std::fabs(cell->face(face)->center()(dim - 1)) >=
+ (std::fabs(cell->face(face)->center()[dim - 1]) >=
0.0 &&
- std::fabs(cell->face(face)->center()(dim - 1)) <= 0.9))
+ std::fabs(cell->face(face)->center()[dim - 1]) <= 0.9))
{
cell->face(face)->set_boundary_id(4);
}
- if ((std::fabs(cell->face(face)->center()(0)) < 1e-12) &&
- (std::fabs(cell->face(face)->center()(dim - 1)) >=
+ if ((std::fabs(cell->face(face)->center()[0]) < 1e-12) &&
+ (std::fabs(cell->face(face)->center()[dim - 1]) >=
0.1 &&
- std::fabs(cell->face(face)->center()(dim - 1)) <= 1.0))
+ std::fabs(cell->face(face)->center()[dim - 1]) <= 1.0))
{
cell->face(face)->set_boundary_id(6);
}
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
{
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
}
return_value[i] = 2.0 * viscosity_function.value(p);
}
{
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
- p(d) = pressure.quadrature_point(q)(d)[i];
+ p[d] = pressure.quadrature_point(q)[d][i];
return_value[i] = 1.0 / viscosity_function.value(p);
}
one_over_viscosity(cell, q) = return_value;
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
{
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
}
return_value[i] = 2.0 * viscosity_function.value(p);
}
{
Point<dim> p;
for (unsigned int d = 0; d < dim; ++d)
- p(d) = velocity.quadrature_point(q)(d)[i];
+ p[d] = velocity.quadrature_point(q)[d][i];
Vector<double> rhs_temp(dim + 1);
right_hand_side.vector_value(p, rhs_temp);
++q_point)
{
cell_pole(dof) += (fe_values.shape_value(dof, q_point) *
- dof_locations[q_point](dof_component % dim));
+ dof_locations[q_point][dof_component % dim]);
}
solution(local_dof_indices[dof]) = 1; // start all solutions at 1
poles(local_dof_indices[dof]) -= cell_pole(dof);
finite_element.system_to_component_index(dof).first;
poles(local_dof_indices[dof]) =
- -dof_locations[dof](dof_component % dim);
+ -dof_locations[dof][dof_component % dim];
if (dof_component == dim) // components start numbering at 0
poles(local_dof_indices[dof]) =
// simplifying the remaining lines
poles(local_dof_indices[dof]) =
- -dof_locations[dof](dof_component % dim);
+ -dof_locations[dof][dof_component % dim];
if (dof_component == dim) // components start numbering at 0
poles(local_dof_indices[dof]) =
if (degree < 2)
{
for (unsigned int d = 0; d < dim; ++d)
- values[k](d) = points[k](d) - d;
+ values[k](d) = points[k][d] - d;
}
else
{
{
Point<dim> p = points[k];
for (unsigned int dd = 0; dd < dim; ++dd)
- p(dd) -= d;
+ p[dd] -= d;
const double r2 = p.square();
values[k](d) = std::pow(r2, (int)degree / 2);
}
ComplexNumber
value(const Position &x, const unsigned int) const
{
- return x(0) + x(1) - ComplexNumber(0.2, 0.7) * x(2);
+ return x[0] + x[1] - ComplexNumber(0.2, 0.7) * x[2];
}
};
operator()(const Point<spacedim> p) const
{
Point<spacedim> q;
- q[0] = std::cos(angle) * p(0) - std::sin(angle) * p(1);
- q[1] = std::sin(angle) * p(0) + std::cos(angle) * p(1);
+ q[0] = std::cos(angle) * p[0] - std::sin(angle) * p[1];
+ q[1] = std::sin(angle) * p[0] + std::cos(angle) * p[1];
for (unsigned d = 2; d < spacedim; ++d)
q[d] = p[d];
return q;
operator()(const Point<spacedim> p) const
{
Point<spacedim> q;
- q[0] = std::cos(angle) * p(0) - std::sin(angle) * p(1);
- q[1] = std::sin(angle) * p(0) + std::cos(angle) * p(1);
+ q[0] = std::cos(angle) * p[0] - std::sin(angle) * p[1];
+ q[1] = std::sin(angle) * p[0] + std::cos(angle) * p[1];
for (unsigned d = 2; d < spacedim; ++d)
q[d] = p[d];
return q;
for (unsigned int i = 0; i < dim; ++i)
{
- position[0](i) = 0.125;
- position[1](i) = 0.525;
+ position[0][i] = 0.125;
+ position[1][i] = 0.525;
}
Particles::Particle<dim, spacedim> particle1(position[0],
for (unsigned int i = 0; i < dim; ++i)
{
- position[0](i) = 0.125;
+ position[0][i] = 0.125;
}
for (unsigned int i = 0; i < spacedim; ++i)
if (Utilities::MPI::this_mpi_process(tr.get_communicator()) == 0)
for (unsigned int i = 0; i < dim; ++i)
- position(i) = 0.475;
+ position[i] = 0.475;
else
for (unsigned int i = 0; i < dim; ++i)
- position(i) = 0.525;
+ position[i] = 0.525;
Particles::Particle<dim, spacedim> particle(
position,
std::vector<Point<dim>> particle_reference_locations(1, Point<dim>());
for (unsigned int i = 0; i < dim; ++i)
- particle_reference_locations[0](i) = 0.5;
+ particle_reference_locations[0][i] = 0.5;
Particles::Generators::regular_reference_locations(
tr, particle_reference_locations, particle_handler);
Particles::Particle<dim> particle;
Point<dim> position;
- position(0) = 1.0;
+ position[0] = 1.0;
particle.set_location(position);
deallog << "Particle location: " << particle.get_location() << std::endl;
{
{
Point<2> position;
- position(0) = 0.3;
- position(1) = 0.5;
+ position[0] = 0.3;
+ position[1] = 0.5;
Point<2> reference_position;
- reference_position(0) = 0.2;
- reference_position(1) = 0.4;
+ reference_position[0] = 0.2;
+ reference_position[1] = 0.4;
const types::particle_index index(7);
Particles::PropertyPool<dim> pool(n_properties_per_particle);
Point<2> position;
- position(0) = 0.3;
- position(1) = 0.5;
+ position[0] = 0.3;
+ position[1] = 0.5;
Point<2> reference_position;
- reference_position(0) = 0.2;
- reference_position(1) = 0.4;
+ reference_position[0] = 0.2;
+ reference_position[1] = 0.4;
const types::particle_index index(7);
{
Point<spacedim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (spacedim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
if (spacedim > 2)
- position(2) = 0.7;
+ position[2] = 0.7;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
if (dim > 2)
- reference_position(2) = 0.6;
+ reference_position[2] = 0.6;
const types::particle_index index(7);
{
Point<spacedim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (spacedim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
if (spacedim > 2)
- position(2) = 0.7;
+ position[2] = 0.7;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
if (dim > 2)
- reference_position(2) = 0.6;
+ reference_position[2] = 0.6;
const types::particle_index index(7);
Point<spacedim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (spacedim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
if (spacedim > 2)
- position(2) = 0.7;
+ position[2] = 0.7;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
if (dim > 2)
- reference_position(2) = 0.6;
+ reference_position[2] = 0.6;
const types::particle_index index(7);
Particles::PropertyPool<dim> pool(n_properties_per_particle);
Point<2> position;
- position(0) = 0.3;
- position(1) = 0.5;
+ position[0] = 0.3;
+ position[1] = 0.5;
Point<2> reference_position;
- reference_position(0) = 0.2;
- reference_position(1) = 0.4;
+ reference_position[0] = 0.2;
+ reference_position[1] = 0.4;
const types::particle_index index(7);
{
{
Point<2> position;
- position(0) = 0.3;
- position(1) = 0.5;
+ position[0] = 0.3;
+ position[1] = 0.5;
Point<2> reference_position;
- reference_position(0) = 0.2;
- reference_position(1) = 0.4;
+ reference_position[0] = 0.2;
+ reference_position[1] = 0.4;
const types::particle_index index(7);
Point<spacedim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (spacedim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
if (spacedim > 2)
- position(2) = 0.7;
+ position[2] = 0.7;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
if (dim > 2)
- reference_position(2) = 0.6;
+ reference_position[2] = 0.6;
Particles::Particle<dim, spacedim> particle(position,
reference_position,
Point<spacedim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (spacedim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
if (spacedim > 2)
- position(2) = 0.7;
+ position[2] = 0.7;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
if (dim > 2)
- reference_position(2) = 0.6;
+ reference_position[2] = 0.6;
Particles::Particle<dim, spacedim> particle(position,
reference_position,
particle_handler.insert_particle(particle, cell_position.first);
particle_handler.insert_particle(particle, cell_position.first);
- position(0) = 0.7;
+ position[0] = 0.7;
Particles::Particle<dim, spacedim> particle2(position,
reference_position,
9);
for (unsigned int i = 0; i < dim; ++i)
{
- position[0](i) = 0.25;
- position[1](i) = 0.75;
+ position[0][i] = 0.25;
+ position[1][i] = 0.75;
}
Particles::Particle<dim, spacedim> particle1(position[0],
// move particle 2 out of the domain. Note that we need to change the
// coordinate dim-1 despite having a spacedim point.
Point<spacedim> shift;
- shift(dim - 1) = 0.5;
+ shift[dim - 1] = 0.5;
for (auto &particle : particle_handler)
particle.set_location(particle.get_location() + shift);
for (unsigned int i = 0; i < dim; ++i)
{
- position[0](i) = 0.125;
- position[1](i) = 0.525;
+ position[0][i] = 0.125;
+ position[1][i] = 0.525;
}
Particles::Particle<dim, spacedim> particle1(position[0],
// move particle 2 out of the domain. Note that we need to change the
// coordinate dim-1 despite having a spacedim point.
Point<spacedim> shift;
- shift(dim - 1) = 0.5;
+ shift[dim - 1] = 0.5;
for (auto &particle : particle_handler)
particle.set_location(particle.get_location() + shift);
if (Utilities::MPI::this_mpi_process(tr.get_communicator()) == 0)
for (unsigned int i = 0; i < dim; ++i)
- position(i) = 0.475;
+ position[i] = 0.475;
else
for (unsigned int i = 0; i < dim; ++i)
- position(i) = 0.525;
+ position[i] = 0.525;
Particles::Particle<dim, spacedim> particle(
position,
if (Utilities::MPI::this_mpi_process(tr.get_communicator()) == 0)
for (unsigned int i = 0; i < dim; ++i)
- position(i) = 0.475;
+ position[i] = 0.475;
else
for (unsigned int i = 0; i < dim; ++i)
- position(i) = 0.525;
+ position[i] = 0.525;
Particles::Particle<dim, spacedim> particle(
position,
if (Utilities::MPI::this_mpi_process(tr.get_communicator()) == 0)
{
for (unsigned int i = 0; i < dim; ++i)
- position(i) = 0.410 + 0.01 * p;
+ position[i] = 0.410 + 0.01 * p;
Particles::Particle<dim, spacedim> particle(
position,
for (unsigned int i = 0; i < dim; ++i)
{
- position[0](i) = 0.125;
- position[1](i) = 0.525;
+ position[0][i] = 0.125;
+ position[1][i] = 0.525;
}
Particles::Particle<dim, spacedim> particle1(position[0],
// move particle 2 out of the domain. Note that we need to change the
// coordinate dim-1 despite having a spacedim point.
Point<spacedim> shift;
- shift(dim - 1) = 0.5;
+ shift[dim - 1] = 0.5;
for (auto &particle : particle_handler)
particle.set_location(particle.get_location() + shift);
Point<spacedim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (spacedim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
if (spacedim > 2)
- position(2) = 0.7;
+ position[2] = 0.7;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
if (dim > 2)
- reference_position(2) = 0.6;
+ reference_position[2] = 0.6;
Particles::Particle<dim, spacedim> particle(position, reference_position, 7);
deallog << "Particle location: " << particle.get_location() << std::endl;
Point<spacedim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (spacedim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
if (spacedim > 2)
- position(2) = 0.7;
+ position[2] = 0.7;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
if (dim > 2)
- reference_position(2) = 0.6;
+ reference_position[2] = 0.6;
Particles::Particle<dim, spacedim> particle(position,
reference_position,
particle_handler.insert_particle(particle, cell_position.first);
particle_handler.insert_particle(particle, cell_position.first);
- position(0) = 0.7;
+ position[0] = 0.7;
Particles::Particle<dim, spacedim> particle2(position,
reference_position,
9);
for (unsigned int i = 0; i < dim; ++i)
{
- position[0](i) = 0.25;
- position[1](i) = 0.75;
+ position[0][i] = 0.25;
+ position[1][i] = 0.75;
}
Particles::Particle<dim, spacedim> particle1(position[0],
// move particle 2 out of the domain. Note that we need to change the
// coordinate dim-1 despite having a spacedim point.
Point<spacedim> shift;
- shift(dim - 1) = 0.5;
+ shift[dim - 1] = 0.5;
for (auto &particle : particle_handler)
particle.set_location(particle.get_location() + shift);
for (unsigned int i = 0; i < dim; ++i)
{
- position[0](i) = 0.125;
- position[1](i) = 0.525;
+ position[0][i] = 0.125;
+ position[1][i] = 0.525;
}
Particles::Particle<dim, spacedim> particle1(position[0],
// move particle 2 out of the domain. Note that we need to change the
// coordinate dim-1 despite having a spacedim point.
Point<spacedim> shift;
- shift(dim - 1) = 0.5;
+ shift[dim - 1] = 0.5;
for (auto &particle : particle_handler)
particle.set_location(particle.get_location() + shift);
n_properties_per_particle);
Point<dim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (dim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
const types::particle_index index(7);
Particles::PropertyPool<dim> pool(n_properties_per_particle);
Point<dim> position;
- position(0) = 0.3;
+ position[0] = 0.3;
if (dim > 1)
- position(1) = 0.5;
+ position[1] = 0.5;
Point<dim> reference_position;
- reference_position(0) = 0.2;
+ reference_position[0] = 0.2;
if (dim > 1)
- reference_position(1) = 0.4;
+ reference_position[1] = 0.4;
const types::particle_index index(7);
for (unsigned int i = 0; i < dim; ++i)
{
- particle_reference_locations[0](i) = 0.25;
- particle_reference_locations[1](i) = 0.5;
- particle_reference_locations[2](i) = 0.75;
+ particle_reference_locations[0][i] = 0.25;
+ particle_reference_locations[1][i] = 0.5;
+ particle_reference_locations[2][i] = 0.75;
}
Particles::Generators::regular_reference_locations(
const double T = 4;
const double t = this->get_time();
- const double px = numbers::PI * point(0);
- const double py = numbers::PI * point(1);
+ const double px = numbers::PI * point[0];
+ const double py = numbers::PI * point[1];
const double pt = numbers::PI / T * t;
values[0] = -2 * cos(pt) * pow(sin(px), 2) * sin(py) * cos(py);
value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
{
if (dim == 2)
- return -2. * M_PI * M_PI * std::sin(M_PI * p(0)) * std::sin(M_PI * p(1));
+ return -2. * M_PI * M_PI * std::sin(M_PI * p[0]) * std::sin(M_PI * p[1]);
else /* if(dim == 3)*/
- return -3. * M_PI * M_PI * std::sin(M_PI * p(0)) * std::sin(M_PI * p(1)) *
- std::sin(M_PI * p(2));
+ return -3. * M_PI * M_PI * std::sin(M_PI * p[0]) * std::sin(M_PI * p[1]) *
+ std::sin(M_PI * p[2]);
}
};
#ifdef DEBUG
// all vertices should be in a plane
for (const auto &vertex : vertices)
- Assert(midpoint(coordinate) == vertex(coordinate), ExcInternalError());
+ Assert(midpoint[coordinate] == vertex[coordinate], ExcInternalError());
#endif
// add another vertex as tip of triangle/pyramid
Point<spacedim> tip = midpoint;
- tip(coordinate) += (face_no % 2 == 1) ? 1. : -1.;
+ tip[coordinate] += (face_no % 2 == 1) ? 1. : -1.;
vertices.push_back(tip);
CellData<dim> simplex(vertices.size());
{
double u = 1.0;
for (int d = 0; d < dim; ++d)
- u *= std::sin(numbers::PI * p(d));
+ u *= std::sin(numbers::PI * p[d]);
return u;
}
{
double u = 1.0;
for (int d = 0; d < dim; ++d)
- u *= std::sin(numbers::PI * p(d));
+ u *= std::sin(numbers::PI * p[d]);
return (1 + 0.0 * dim * numbers::PI * numbers::PI) * u;
}
{
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;
}
{
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;
}
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);
}
}
Assert(dim >= 2, ExcNotImplemented());
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)
{
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.;
Assert(dim >= 2, ExcNotImplemented());
Point<dim> wind_field;
- wind_field(0) = -p(1);
- wind_field(1) = p(0);
+ wind_field[0] = -p[1];
+ wind_field[1] = p[0];
if (wind_field.norm() > 1e-10)
wind_field /= wind_field.norm();
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.;
Assert(dim >= 2, ExcNotImplemented());
Point<dim> wind_field;
- wind_field(0) = -p(1);
- wind_field(1) = p(0);
+ wind_field[0] = -p[1];
+ wind_field[1] = p[0];
if (wind_field.norm() > 1e-6)
wind_field /= wind_field.norm();
const Point<dim> &p = points[i];
Point<dim> &beta = values[i];
- beta(0) = -p(1);
- beta(1) = p(0);
+ beta[0] = -p[1];
+ beta[1] = p[0];
if (beta.norm() > 1e-10)
beta /= std::sqrt(beta.square());
for (unsigned int i = 0; i < std::min(values.size(), points.size()); ++i)
{
- if (points[i](0) < 0.5)
+ if (points[i][0] < 0.5)
values[i] = 1.;
else
values[i] = 0.;
const Point<dim> &p = points[i];
Point<dim> &beta = values[i];
- beta(0) = -p(1);
- beta(1) = p(0);
+ beta[0] = -p[1];
+ beta[1] = p[0];
if (beta.norm() > 1e-10)
beta /= std::sqrt(beta.square());
for (unsigned int i = 0; i < std::min(values.size(), points.size()); ++i)
{
- if (points[i](0) < 0.5)
+ if (points[i][0] < 0.5)
values[i] = 1.;
else
values[i] = 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))
double
Solution<2>::value(const Point<2> &p, const unsigned int) const
{
- return (-2. * p(0) * p(1));
+ return (-2. * p[0] * p[1]);
}
template <>
Solution<2>::gradient(const Point<2> &p, 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;
}
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]));
}
template <>
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;
}
RightHandSide<2>::value(const Point<2> &p,
const unsigned int /*component*/) const
{
- return (-8. * p(0) * p(1));
+ return (-8. * p[0] * p[1]);
}
template <>
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();
const double T = 4;
const double t = this->get_time();
- const double px = numbers::PI * point(0);
- const double py = numbers::PI * point(1);
+ const double px = numbers::PI * point[0];
+ const double py = numbers::PI * point[1];
const double pt = numbers::PI / T * t;
values[0] = -2 * cos(pt) * pow(sin(px), 2) * sin(py) * cos(py);
deallog
<< "This sentence should appear once when the corresponding face is visited only once on cell "
<< cell->index() << std::endl;
- return (std::abs(cell->face(face_index)->center()(0)) < 0.01);
+ return (std::abs(cell->face(face_index)->center()[0]) < 0.01);
}
template <int dim>
const Point<dim> &p = cell->center();
bool negative = true;
for (unsigned int d = 0; d < dim; ++d)
- if (p(d) >= 0.)
+ if (p[d] >= 0.)
negative = false;
if (negative)
cell->set_refine_flag();
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 4 * std::pow(p(i), 4);
+ return_value += 4 * std::pow(p[i], 4);
return return_value;
}
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 2 * std::pow(p(i), 2);
+ return_value += 2 * std::pow(p[i], 2);
return return_value;
}
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 4 * std::pow(p(i), 4);
+ return_value += 4 * std::pow(p[i], 4);
return return_value;
}
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 4 * std::pow(p(i), 4);
+ return_value += 4 * std::pow(p[i], 4);
return return_value;
}
{
double return_value = 0;
for (unsigned int i = 0; i < dim; ++i)
- return_value += 4 * std::pow(p(i), 4);
+ return_value += 4 * std::pow(p[i], 4);
return return_value;
}
projects / dealii.git / commitdiff