Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Replace a ReferenceCell function call by a non-deprecated one. #14770

Merged
merged 1 commit into from
Feb 13, 2023
Merged

Replace a ReferenceCell function call by a non-deprecated one. #14770

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
192 changes: 94 additions & 98 deletions source/base/qprojector.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -644,133 +644,129 @@ QProjector<3>::project_to_all_faces(const ReferenceCell & reference_cell,
const hp::QCollection<2> &quadrature)
{
const auto support_points_tri =
[](const auto & face,
[](const std::pair<std::vector<Point<3>>, double> &face,
const unsigned char orientation) -> std::vector<Point<3>> {
std::array<Point<3>, 3> vertices;
std::copy_n(face.first.begin(), face.first.size(), vertices.begin());
const auto temp =
ReferenceCells::Triangle.permute_according_orientation(vertices,
orientation);
return std::vector<Point<3>>(temp.begin(),
temp.begin() + face.first.size());
const boost::container::small_vector<Point<3>, 8> temp =
ReferenceCells::Triangle.permute_by_combined_orientation<Point<3>>(
face.first, orientation);
return std::vector<Point<3>>(temp.begin(), temp.end());
};

const auto support_points_quad =
[](const auto & face,
[](const std::pair<std::vector<Point<3>>, double> &face,
const unsigned char orientation) -> std::vector<Point<3>> {
std::array<Point<3>, 4> vertices;
std::copy_n(face.first.begin(), face.first.size(), vertices.begin());
const auto temp =
ReferenceCells::Quadrilateral.permute_according_orientation(vertices,
orientation);
return std::vector<Point<3>>(temp.begin(),
temp.begin() + face.first.size());
const boost::container::small_vector<Point<3>, 8> temp =
ReferenceCells::Quadrilateral.permute_by_combined_orientation<Point<3>>(
face.first, orientation);
return std::vector<Point<3>>(temp.begin(), temp.end());
};

const auto process = [&](const auto &faces) {
// new (projected) quadrature points and weights
std::vector<Point<3>> points;
std::vector<double> weights;

const auto poly_tri = BarycentricPolynomials<2>::get_fe_p_basis(1);
const TensorProductPolynomials<2> poly_quad(
Polynomials::generate_complete_Lagrange_basis(
{Point<1>(0.0), Point<1>(1.0)}));

// loop over all faces (triangles) ...
for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
{
// linear polynomial to map the reference quadrature points correctly
// on faces
const unsigned int n_shape_functions = faces[face_no].first.size();
const auto process =
[&](const std::vector<std::pair<std::vector<Point<3>>, double>> &faces) {
// new (projected) quadrature points and weights
std::vector<Point<3>> points;
std::vector<double> weights;

const auto &poly =
n_shape_functions == 3 ?
static_cast<const ScalarPolynomialsBase<2> &>(poly_tri) :
static_cast<const ScalarPolynomialsBase<2> &>(poly_quad);
const auto poly_tri = BarycentricPolynomials<2>::get_fe_p_basis(1);
const TensorProductPolynomials<2> poly_quad(
Polynomials::generate_complete_Lagrange_basis(
{Point<1>(0.0), Point<1>(1.0)}));

// ... and over all possible orientations
for (unsigned char orientation = 0;
orientation < reference_cell.n_face_orientations(face_no);
++orientation)
{
const auto &face = faces[face_no];
// loop over all faces (triangles) ...
for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
{
// linear polynomial to map the reference quadrature points correctly
// on faces
const unsigned int n_shape_functions = faces[face_no].first.size();

const auto &poly =
n_shape_functions == 3 ?
static_cast<const ScalarPolynomialsBase<2> &>(poly_tri) :
static_cast<const ScalarPolynomialsBase<2> &>(poly_quad);

// ... and over all possible orientations
for (unsigned char orientation = 0;
orientation < reference_cell.n_face_orientations(face_no);
++orientation)
{
const auto &face = faces[face_no];

const auto support_points =
n_shape_functions == 3 ? support_points_tri(face, orientation) :
support_points_quad(face, orientation);
const auto support_points =
n_shape_functions == 3 ? support_points_tri(face, orientation) :
support_points_quad(face, orientation);

// the quadrature rule to be projected ...
const auto &sub_quadrature_points =
quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
const auto &sub_quadrature_weights =
quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();
// the quadrature rule to be projected ...
const auto &sub_quadrature_points =
quadrature[quadrature.size() == 1 ? 0 : face_no].get_points();
const auto &sub_quadrature_weights =
quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights();

// loop over all quadrature points
for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
{
Point<3> mapped_point;
// loop over all quadrature points
for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j)
{
Point<3> mapped_point;

// map reference quadrature point
for (unsigned int i = 0; i < n_shape_functions; ++i)
mapped_point +=
support_points[i] *
poly.compute_value(i, sub_quadrature_points[j]);
// map reference quadrature point
for (unsigned int i = 0; i < n_shape_functions; ++i)
mapped_point +=
support_points[i] *
poly.compute_value(i, sub_quadrature_points[j]);

points.push_back(mapped_point);
points.push_back(mapped_point);

// scale quadrature weight
const double scaling = [&]() {
const auto & supp_pts = support_points;
const unsigned int dim_ = 2;
const unsigned int spacedim = 3;
// scale quadrature weight
const double scaling = [&]() {
const auto & supp_pts = support_points;
const unsigned int dim_ = 2;
const unsigned int spacedim = 3;

double result[spacedim][dim_];
double result[spacedim][dim_];

std::vector<Tensor<1, dim_>> shape_derivatives(
n_shape_functions);
std::vector<Tensor<1, dim_>> shape_derivatives(
n_shape_functions);

for (unsigned int i = 0; i < n_shape_functions; ++i)
shape_derivatives[i] =
poly.compute_1st_derivative(i, sub_quadrature_points[j]);
for (unsigned int i = 0; i < n_shape_functions; ++i)
shape_derivatives[i] =
poly.compute_1st_derivative(i,
sub_quadrature_points[j]);

for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim_; ++j)
result[i][j] = shape_derivatives[0][j] * supp_pts[0][i];
for (unsigned int k = 1; k < n_shape_functions; ++k)
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim_; ++j)
result[i][j] +=
shape_derivatives[k][j] * supp_pts[k][i];
result[i][j] = shape_derivatives[0][j] * supp_pts[0][i];
for (unsigned int k = 1; k < n_shape_functions; ++k)
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim_; ++j)
result[i][j] +=
shape_derivatives[k][j] * supp_pts[k][i];

DerivativeForm<1, dim_, spacedim> contravariant;
DerivativeForm<1, dim_, spacedim> contravariant;

for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim_; ++j)
contravariant[i][j] = result[i][j];
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim_; ++j)
contravariant[i][j] = result[i][j];


Tensor<1, spacedim> DX_t[dim_];
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim_; ++j)
DX_t[j][i] = contravariant[i][j];
Tensor<1, spacedim> DX_t[dim_];
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim_; ++j)
DX_t[j][i] = contravariant[i][j];

Tensor<2, dim_> G;
for (unsigned int i = 0; i < dim_; ++i)
for (unsigned int j = 0; j < dim_; ++j)
G[i][j] = DX_t[i] * DX_t[j];
Tensor<2, dim_> G;
for (unsigned int i = 0; i < dim_; ++i)
for (unsigned int j = 0; j < dim_; ++j)
G[i][j] = DX_t[i] * DX_t[j];

return std::sqrt(determinant(G));
}();
return std::sqrt(determinant(G));
}();

weights.push_back(sub_quadrature_weights[j] * scaling);
}
}
}
weights.push_back(sub_quadrature_weights[j] * scaling);
}
}
}

// construct new quadrature rule
return Quadrature<3>(points, weights);
};
// construct new quadrature rule
return Quadrature<3>(points, weights);
};

if (reference_cell == ReferenceCells::Tetrahedron)
{
Expand Down