lex_normal.push_back(i);
lex_normal.push_back(dofs_per_face_normal);
- // 'direction' distingusishes between normal and tangential direction
+ // 'direction' distinguishes between normal and tangential direction
for (unsigned int direction = 0; direction < 2; ++direction)
{
UnivariateShapeData<Number> &univariate_shape_data =
return;
// grant write access to common univariate shape data
- auto &shape_values = univariate_shape_data.shape_values;
- auto &shape_gradients = univariate_shape_data.shape_gradients;
- auto &shape_hessians = univariate_shape_data.shape_hessians;
- auto &shape_gradients_collocation =
- univariate_shape_data.shape_gradients_collocation;
- auto &shape_hessians_collocation =
- univariate_shape_data.shape_hessians_collocation;
- auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
- auto &shape_data_on_face = univariate_shape_data.shape_data_on_face;
- auto &quadrature_data_on_face =
- univariate_shape_data.quadrature_data_on_face;
+ auto &shape_values = univariate_shape_data.shape_values;
+ auto &shape_gradients = univariate_shape_data.shape_gradients;
+ auto &shape_hessians = univariate_shape_data.shape_hessians;
+ auto &shape_data_on_face = univariate_shape_data.shape_data_on_face;
auto &values_within_subface = univariate_shape_data.values_within_subface;
auto &gradients_within_subface =
univariate_shape_data.gradients_within_subface;
fe.shape_grad_grad(my_i, q_point)[0][0];
}
- if (n_q_points_1d < 200)
- {
- quadrature_data_on_face[0].resize(quad.size() * 3);
- quadrature_data_on_face[1].resize(quad.size() * 3);
-
- const std::vector<Polynomials::Polynomial<double>> poly_coll =
- Polynomials::generate_complete_Lagrange_basis(quad.get_points());
-
- for (unsigned int i = 0; i < quad.size(); ++i)
- {
- std::array<double, 3> values;
- poly_coll[i].value(0.0, 2, values.data());
- for (unsigned int d = 0; d < 3; ++d)
- quadrature_data_on_face[0][i + d * quad.size()] = values[d];
- poly_coll[i].value(1.0, 2, values.data());
- for (unsigned int d = 0; d < 3; ++d)
- quadrature_data_on_face[1][i + d * quad.size()] = values[d];
- }
- }
-
if (dim > 1 && (dynamic_cast<const FE_Q<dim> *>(&fe) ||
dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe)))
{
}
}
- // get gradient and Hessian transformation matrix for the polynomial
- // space associated with the quadrature rule (collocation space). We
- // need to avoid the case with more than a few hundreds of quadrature
- // points when the Lagrange polynomials might underflow. Note that 200
- // is not an exact value, as different quadrature formulas behave
- // slightly differently, but 200 has been observed to be low enough for
- // all common quadrature formula types. For QGauss, the actual limit is
- // 517 points, for example.
- if (n_q_points_1d < 200)
- {
- shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d);
- shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d);
- const std::vector<Polynomials::Polynomial<double>> poly_coll =
- Polynomials::generate_complete_Lagrange_basis(quad.get_points());
- std::array<double, 3> values;
- for (unsigned int i = 0; i < n_q_points_1d; ++i)
- for (unsigned int q = 0; q < n_q_points_1d; ++q)
- {
- poly_coll[i].value(quad.get_points()[q][0], 2, values.data());
- shape_gradients_collocation[i * n_q_points_1d + q] = values[1];
- shape_hessians_collocation[i * n_q_points_1d + q] = values[2];
- }
-
- // compute the inverse shape functions in three steps: we first
- // change from the given quadrature formula and the associated
- // Lagrange polynomials to the Lagrange polynomials at quadrature
- // points. in this basis, we can then perform the second step, which
- // is the computation of a projection matrix from the potentially
- // higher polynomial degree associated to the quadrature points to a
- // polynomial space of degree equal to the degree of the given
- // elements. in the third step, we change from the Lagrange
- // polynomials in the Gauss quadrature points to the polynomial
- // space of the given element
-
- // step 1: change basis from the Lagrange polynomials at the given
- // quadrature points to the Lagrange basis at Gauss quadrature
- // points. this is often the identity operation as we often compute
- // with Gaussian quadrature, but not necessarily so
- QGauss<1> quad_gauss(n_q_points_1d);
- FullMatrix<double> transform_to_gauss(n_q_points_1d, n_q_points_1d);
- for (unsigned int i = 0; i < n_q_points_1d; ++i)
- for (unsigned int j = 0; j < n_q_points_1d; ++j)
- transform_to_gauss(i, j) =
- poly_coll[j].value(quad_gauss.point(i)[0]);
-
- // step 2: computation for the projection (in reference coordinates)
- // from higher to lower polynomial degree
- //
- // loop over quadrature points, multiply by q-weight on high degree
- // integrate loop going from high degree to low degree loop over new
- // points, multiply by inverse q-weight on low degree
- //
- // This projection step is for the special case of Lagrange
- // polynomials where most of the interpolation matrices are unit
- // matrices when applying the inverse mass matrix, so we do not need
- // to compute much.
- QGauss<1> quad_project(n_dofs_1d);
- const std::vector<Polynomials::Polynomial<double>> poly_project =
- Polynomials::generate_complete_Lagrange_basis(
- quad_project.get_points());
-
- FullMatrix<double> project_gauss(n_dofs_1d, n_q_points_1d);
-
- for (unsigned int i = 0; i < n_dofs_1d; ++i)
- for (unsigned int q = 0; q < n_q_points_1d; ++q)
- project_gauss(i, q) =
- poly_project[i].value(quad_gauss.get_points()[q][0]) *
- (quad_gauss.weight(q) / quad_project.weight(i));
- FullMatrix<double> project_to_dof_space(n_dofs_1d, n_q_points_1d);
- project_gauss.mmult(project_to_dof_space, transform_to_gauss);
-
- // step 3: change the basis back to the given finite element
- // space. we can use a shortcut for elements that define support
- // points, in which case we can evaluate the Lagrange polynomials of
- // the Gauss quadrature in those points. this will give more
- // accurate results than the inversion of a matrix. for more general
- // polynomial spaces, we must invert a matrix of a Vandermonde type,
- // which we do by a Householder transformation to keep roundoff
- // errors low.
- inverse_shape_values.resize_fast(array_size);
- FullMatrix<double> transform_from_gauss(n_dofs_1d, n_dofs_1d);
- if (fe.has_support_points())
- {
- for (unsigned int i = 0; i < n_dofs_1d; ++i)
- for (unsigned int j = 0; j < n_dofs_1d; ++j)
- transform_from_gauss(i, j) = poly_project[j].value(
- fe.get_unit_support_points()[scalar_lexicographic[i]][0]);
- FullMatrix<double> result(n_dofs_1d, n_q_points_1d);
- transform_from_gauss.mmult(result, project_to_dof_space);
-
- // set very small entries to zero - we are in reference space
- // with normalized numbers, so this is straight-forward to check
- // here
- for (unsigned int i = 0; i < n_dofs_1d; ++i)
- for (unsigned int q = 0; q < n_q_points_1d; ++q)
- inverse_shape_values[i * n_q_points_1d + q] =
- std::abs(result(i, q)) < 1e-15 ? 0 : result(i, q);
- }
- else
- {
- for (unsigned int i = 0; i < n_dofs_1d; ++i)
- for (unsigned int j = 0; j < n_dofs_1d; ++j)
- {
- Point<dim> q_point = unit_point;
- q_point[0] = quad_project.point(i)[0];
-
- transform_from_gauss(i, j) =
- fe.shape_value(scalar_lexicographic[j], q_point);
- }
- Householder<double> H(transform_from_gauss);
- Vector<double> in(n_dofs_1d), out(n_dofs_1d);
- for (unsigned int q = 0; q < n_q_points_1d; ++q)
- {
- for (unsigned int i = 0; i < n_dofs_1d; ++i)
- in(i) = project_to_dof_space(i, q);
- H.least_squares(out, in);
- for (unsigned int i = 0; i < n_dofs_1d; ++i)
- inverse_shape_values[i * n_q_points_1d + q] =
- std::abs(out(i)) < 1e-15 ? 0. : out(i);
- }
- }
- }
+ univariate_shape_data.evaluate_collocation_space(fe,
+ quad,
+ scalar_lexicographic);
if (element_type == tensor_general &&
- check_1d_shapes_symmetric(univariate_shape_data))
+ univariate_shape_data.check_and_set_shapes_symmetric())
{
if (dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe) &&
fe.tensor_degree() > 1)
element_type = tensor_symmetric_no_collocation;
- else if (check_1d_shapes_collocation(univariate_shape_data))
+ else if (univariate_shape_data.check_shapes_collocation())
element_type = tensor_symmetric_collocation;
else
element_type = tensor_symmetric;
}
}
else if (element_type == tensor_symmetric_plus_dg0)
- check_1d_shapes_symmetric(univariate_shape_data);
+ univariate_shape_data.check_and_set_shapes_symmetric();
nodal_at_cell_boundaries = true;
for (unsigned int i = 1; i < n_dofs_1d; ++i)
+ template <typename Number>
+ template <int dim, int spacedim>
+ void
+ UnivariateShapeData<Number>::evaluate_collocation_space(
+ const FiniteElement<dim, spacedim> &fe,
+ const Quadrature<1> & quad,
+ const std::vector<unsigned int> & lexicographic)
+ {
+ const unsigned int n_q_points_1d = quad.size();
+ const unsigned int n_dofs_1d =
+ std::min(fe.n_dofs_per_cell(), fe_degree + 1);
+
+ // get gradient and Hessian transformation matrix for the polynomial
+ // space associated with the quadrature rule (collocation space). We
+ // need to avoid the case with more than a few hundreds of quadrature
+ // points when the Lagrange polynomials might underflow. Note that 200
+ // is not an exact value, as different quadrature formulas behave
+ // slightly differently, but 200 has been observed to be low enough for
+ // all common quadrature formula types. For QGauss, the actual limit is
+ // 517 points, for example.
+ if (n_q_points_1d >= 200)
+ return;
+
+ shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d);
+ shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d);
+ const std::vector<Polynomials::Polynomial<double>> poly_coll =
+ Polynomials::generate_complete_Lagrange_basis(quad.get_points());
+ std::array<double, 3> values;
+ for (unsigned int i = 0; i < n_q_points_1d; ++i)
+ for (unsigned int q = 0; q < n_q_points_1d; ++q)
+ {
+ poly_coll[i].value(quad.get_points()[q][0], 2, values.data());
+ shape_gradients_collocation[i * n_q_points_1d + q] = values[1];
+ shape_hessians_collocation[i * n_q_points_1d + q] = values[2];
+ }
+
+ // compute the inverse shape functions in three steps: we first
+ // change from the given quadrature formula and the associated
+ // Lagrange polynomials to the Lagrange polynomials at quadrature
+ // points. in this basis, we can then perform the second step, which
+ // is the computation of a projection matrix from the potentially
+ // higher polynomial degree associated to the quadrature points to a
+ // polynomial space of degree equal to the degree of the given
+ // elements. in the third step, we change from the Lagrange
+ // polynomials in the Gauss quadrature points to the polynomial
+ // space of the given element
+
+ // step 1: change basis from the Lagrange polynomials at the given
+ // quadrature points to the Lagrange basis at Gauss quadrature
+ // points. this is often the identity operation as we often compute
+ // with Gaussian quadrature, but not necessarily so
+ QGauss<1> quad_gauss(n_q_points_1d);
+ FullMatrix<double> transform_to_gauss(n_q_points_1d, n_q_points_1d);
+ for (unsigned int i = 0; i < n_q_points_1d; ++i)
+ for (unsigned int j = 0; j < n_q_points_1d; ++j)
+ transform_to_gauss(i, j) = poly_coll[j].value(quad_gauss.point(i)[0]);
+
+ // step 2: computation for the projection (in reference coordinates)
+ // from higher to lower polynomial degree
+ //
+ // loop over quadrature points, multiply by q-weight on high degree
+ // integrate loop going from high degree to low degree loop over new
+ // points, multiply by inverse q-weight on low degree
+ //
+ // This projection step is for the special case of Lagrange
+ // polynomials where most of the interpolation matrices are unit
+ // matrices when applying the inverse mass matrix, so we do not need
+ // to compute much.
+ QGauss<1> quad_project(n_dofs_1d);
+ const std::vector<Polynomials::Polynomial<double>> poly_project =
+ Polynomials::generate_complete_Lagrange_basis(
+ quad_project.get_points());
+
+ FullMatrix<double> project_gauss(n_dofs_1d, n_q_points_1d);
+
+ for (unsigned int i = 0; i < n_dofs_1d; ++i)
+ for (unsigned int q = 0; q < n_q_points_1d; ++q)
+ project_gauss(i, q) =
+ poly_project[i].value(quad_gauss.get_points()[q][0]) *
+ (quad_gauss.weight(q) / quad_project.weight(i));
+ FullMatrix<double> project_to_dof_space(n_dofs_1d, n_q_points_1d);
+ project_gauss.mmult(project_to_dof_space, transform_to_gauss);
+
+ // step 3: change the basis back to the given finite element
+ // space. we can use a shortcut for elements that define support
+ // points, in which case we can evaluate the Lagrange polynomials of
+ // the Gauss quadrature in those points. this will give more
+ // accurate results than the inversion of a matrix. for more general
+ // polynomial spaces, we must invert a matrix of a Vandermonde type,
+ // which we do by a Householder transformation to keep roundoff
+ // errors low.
+ inverse_shape_values.resize_fast(shape_values.size());
+ FullMatrix<double> transform_from_gauss(n_dofs_1d, n_dofs_1d);
+ if (fe.has_support_points())
+ {
+ for (unsigned int i = 0; i < n_dofs_1d; ++i)
+ for (unsigned int j = 0; j < n_dofs_1d; ++j)
+ transform_from_gauss(i, j) = poly_project[j].value(
+ fe.get_unit_support_points()[lexicographic[i]][0]);
+ FullMatrix<double> result(n_dofs_1d, n_q_points_1d);
+ transform_from_gauss.mmult(result, project_to_dof_space);
+
+ // set very small entries to zero - we are in reference space
+ // with normalized numbers, so this is straight-forward to check
+ // here
+ for (unsigned int i = 0; i < n_dofs_1d; ++i)
+ for (unsigned int q = 0; q < n_q_points_1d; ++q)
+ inverse_shape_values[i * n_q_points_1d + q] =
+ std::abs(result(i, q)) < 1e-15 ? 0 : result(i, q);
+ }
+ else
+ {
+ for (unsigned int i = 0; i < n_dofs_1d; ++i)
+ for (unsigned int j = 0; j < n_dofs_1d; ++j)
+ {
+ Point<dim> q_point;
+ q_point[0] = quad_project.point(i)[0];
+
+ transform_from_gauss(i, j) =
+ fe.shape_value(lexicographic[j], q_point);
+ }
+ Householder<double> H(transform_from_gauss);
+ Vector<double> in(n_dofs_1d), out(n_dofs_1d);
+ for (unsigned int q = 0; q < n_q_points_1d; ++q)
+ {
+ for (unsigned int i = 0; i < n_dofs_1d; ++i)
+ in(i) = project_to_dof_space(i, q);
+ H.least_squares(out, in);
+ for (unsigned int i = 0; i < n_dofs_1d; ++i)
+ inverse_shape_values[i * n_q_points_1d + q] =
+ std::abs(out(i)) < 1e-15 ? 0. : out(i);
+ }
+ }
+ quadrature_data_on_face[0].resize(quad.size() * 3);
+ quadrature_data_on_face[1].resize(quad.size() * 3);
+
+ for (unsigned int i = 0; i < quad.size(); ++i)
+ {
+ std::array<double, 3> values;
+ poly_coll[i].value(0.0, 2, values.data());
+ for (unsigned int d = 0; d < 3; ++d)
+ quadrature_data_on_face[0][i + d * quad.size()] = values[d];
+ poly_coll[i].value(1.0, 2, values.data());
+ for (unsigned int d = 0; d < 3; ++d)
+ quadrature_data_on_face[1][i + d * quad.size()] = values[d];
+ }
+ }
+
+
+
template <typename Number>
bool
- ShapeInfo<Number>::check_1d_shapes_symmetric(
- UnivariateShapeData<Number> &univariate_shape_data)
+ UnivariateShapeData<Number>::check_and_set_shapes_symmetric()
{
- if (dofs_per_component_on_cell == 0)
+ if (fe_degree == 0)
return false;
- const auto n_q_points_1d = univariate_shape_data.n_q_points_1d;
- const auto fe_degree = univariate_shape_data.fe_degree;
- auto & shape_values = univariate_shape_data.shape_values;
- auto & shape_gradients = univariate_shape_data.shape_gradients;
- auto & shape_hessians = univariate_shape_data.shape_hessians;
- auto & shape_gradients_collocation =
- univariate_shape_data.shape_gradients_collocation;
- auto &shape_hessians_collocation =
- univariate_shape_data.shape_hessians_collocation;
- auto &shape_values_eo = univariate_shape_data.shape_values_eo;
- auto &shape_gradients_eo = univariate_shape_data.shape_gradients_eo;
- auto &shape_hessians_eo = univariate_shape_data.shape_hessians_eo;
- auto &shape_gradients_collocation_eo =
- univariate_shape_data.shape_gradients_collocation_eo;
- auto &shape_hessians_collocation_eo =
- univariate_shape_data.shape_hessians_collocation_eo;
- auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
- auto &inverse_shape_values_eo =
- univariate_shape_data.inverse_shape_values_eo;
-
const double zero_tol =
std::is_same_v<Number, double> == true ? 1e-12 : 1e-7;
// symmetry for values
template <typename Number>
bool
- ShapeInfo<Number>::check_1d_shapes_collocation(
- const UnivariateShapeData<Number> &univariate_shape_data) const
+ UnivariateShapeData<Number>::check_shapes_collocation() const
{
- if (dofs_per_component_on_cell != n_q_points)
+ if (fe_degree + 1 != n_q_points_1d)
return false;
- const auto fe_degree = univariate_shape_data.fe_degree;
- auto & shape_values = univariate_shape_data.shape_values;
-
const double zero_tol =
std::is_same_v<Number, double> == true ? 1e-12 : 1e-7;
// check: identity operation for shape values
projects / dealii.git / commitdiff