template <int dim, int spacedim>
void
FE_NedelecSZ<dim, spacedim>::evaluate(
- const std::vector<Point<dim>> /*p_list*/,
- const UpdateFlags /*update_flags*/,
+ const std::vector<Point<dim>> &p_list,
+ const UpdateFlags update_flags,
std::unique_ptr<
- typename dealii::FiniteElement<dim, spacedim>::InternalDataBase>
- & /*data_ptr*/) const
-{
- DEAL_II_NOT_IMPLEMENTED();
-}
-
-
-
-template <>
-void
-FE_NedelecSZ<2, 2>::evaluate(
- const std::vector<Point<2>> &p_list,
- const UpdateFlags update_flags,
- std::unique_ptr<typename dealii::FiniteElement<2, 2>::InternalDataBase>
- &data_ptr) const
+ typename dealii::FiniteElement<dim, spacedim>::InternalDataBase> &data_ptr)
+ const
{
auto &data = dynamic_cast<InternalData &>(*data_ptr);
data.update_each = requires_update_flags(update_flags);
// Useful quantities:
- const unsigned int dim = 2;
const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1
- const unsigned int vertices_per_cell =
- GeometryInfo<2 /*dim*/>::vertices_per_cell;
- const unsigned int lines_per_cell = GeometryInfo<2 /*dim*/>::lines_per_cell;
+ const unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
+ const unsigned int lines_per_cell = GeometryInfo<dim>::lines_per_cell;
+ const unsigned int faces_per_cell = GeometryInfo<dim>::faces_per_cell;
const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
+ // we assume that all quads have the same number of dofs
+ const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
const unsigned int n_q_points = p_list.size();
std::vector<std::vector<double>> lambda(n_q_points,
std::vector<double>(lines_per_cell));
- for (unsigned int q = 0; q < n_q_points; ++q)
+ switch (dim)
{
- sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
- sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
- sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
- sigma[q][3] = p_list[q][0] + p_list[q][1];
-
- lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
- lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
- lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
- lambda[q][3] = p_list[q][0] * p_list[q][1];
- for (unsigned int i = 0; i < vertices_per_cell; ++i)
- for (unsigned int j = 0; j < vertices_per_cell; ++j)
- data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
- }
+ case 2:
+ {
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
+ sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
+ sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
+ sigma[q][3] = p_list[q][0] + p_list[q][1];
+
+ lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
+ lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
+ lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
+ lambda[q][3] = p_list[q][0] * p_list[q][1];
+ for (unsigned int i = 0; i < vertices_per_cell; ++i)
+ for (unsigned int j = 0; j < vertices_per_cell; ++j)
+ data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
+ }
- // Calculate the gradient of sigma_imj_values[q][i][j] =
- // sigma[q][i]-sigma[q][j]
- // - this depends on the component and the direction of the
- // corresponding edge.
- // - the direction of the edge is determined by
- // sigma_imj_sign[i][j].
- // Helper arrays:
- const int sigma_comp_signs[vertices_per_cell][dim] = {{-1, -1},
- {1, -1},
- {-1, 1},
- {1, 1}};
- int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
- unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell];
-
- for (unsigned int i = 0; i < vertices_per_cell; ++i)
- for (unsigned int j = 0; j < vertices_per_cell; ++j)
- {
- // sigma_imj_sign is the sign (+/-) of the coefficient of
- // x/y/z in sigma_imj_values Due to the numbering of vertices
- // on the reference element it is easy to find edges in the
- // positive direction are from smaller to higher local vertex
- // numbering.
- sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
- sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
-
- // Now store the component which the sigma_i - sigma_j
- // corresponds to:
- sigma_imj_component[i][j] = 0;
- for (unsigned int d = 0; d < dim; ++d)
- {
- int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
- // Only interested in the first non-zero
- // as if there is a second, it can not be a valid edge.
- if (temp_imj != 0)
+ // Calculate the gradient of sigma_imj_values[q][i][j] =
+ // sigma[q][i]-sigma[q][j]
+ // - this depends on the component and the direction of the
+ // corresponding edge.
+ // - the direction of the edge is determined by
+ // sigma_imj_sign[i][j].
+ // Helper arrays:
+ const int sigma_comp_signs[GeometryInfo<2>::vertices_per_cell][2] = {
+ {-1, -1}, {1, -1}, {-1, 1}, {1, 1}};
+ int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
+ unsigned int sigma_imj_component[vertices_per_cell]
+ [vertices_per_cell];
+
+ for (unsigned int i = 0; i < vertices_per_cell; ++i)
+ for (unsigned int j = 0; j < vertices_per_cell; ++j)
{
- sigma_imj_component[i][j] = d;
- break;
+ // sigma_imj_sign is the sign (+/-) of the coefficient of
+ // x/y/z in sigma_imj_values Due to the numbering of vertices
+ // on the reference element it is easy to find edges in the
+ // positive direction are from smaller to higher local vertex
+ // numbering.
+ sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
+ sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
+
+ // Now store the component which the sigma_i - sigma_j
+ // corresponds to:
+ sigma_imj_component[i][j] = 0;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ int temp_imj =
+ sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
+ // Only interested in the first non-zero
+ // as if there is a second, it can not be a valid edge.
+ if (temp_imj != 0)
+ {
+ sigma_imj_component[i][j] = d;
+ break;
+ }
+ }
+ // Can now calculate the gradient, only non-zero in the
+ // component given: Note some i,j combinations will be
+ // incorrect, but only on invalid edges.
+ data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
+ 2.0 * sigma_imj_sign[i][j];
}
- }
- // Can now calculate the gradient, only non-zero in the
- // component given: Note some i,j combinations will be
- // incorrect, but only on invalid edges.
- data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
- 2.0 * sigma_imj_sign[i][j];
- }
- // Now compute the edge parameterisations for a single element
- // with global numbering matching that of the reference element:
+ // Now compute the edge parameterisations for a single element
+ // with global numbering matching that of the reference element:
- // Resize the edge parameterisations
- data.edge_sigma_values.resize(lines_per_cell,
- std::vector<double>(n_q_points));
- data.edge_sigma_grads.resize(lines_per_cell, std::vector<double>(dim));
+ // Resize the edge parameterisations
+ data.edge_sigma_values.resize(lines_per_cell,
+ std::vector<double>(n_q_points));
+ data.edge_sigma_grads.resize(lines_per_cell,
+ std::vector<double>(dim));
- // Fill the values for edge lambda and edge sigma:
- const unsigned int edge_sigma_direction[lines_per_cell] = {1, 1, 0, 0};
+ // Fill the values for edge lambda and edge sigma:
+ const unsigned int
+ edge_sigma_direction[GeometryInfo<2>::lines_per_cell] = {1,
+ 1,
+ 0,
+ 0};
- data.edge_lambda_values.resize(lines_per_cell,
- std::vector<double>(n_q_points));
+ data.edge_lambda_values.resize(lines_per_cell,
+ std::vector<double>(n_q_points));
- data.edge_lambda_grads_2d.resize(lines_per_cell, std::vector<double>(dim));
+ data.edge_lambda_grads_2d.resize(lines_per_cell,
+ std::vector<double>(dim));
- for (unsigned int m = 0; m < lines_per_cell; ++m)
- {
- // e1=max(reference vertex numbering on this edge)
- // e2=min(reference vertex numbering on this edge)
- // Which is guaranteed to be:
- const unsigned int e1(GeometryInfo<dim>::line_to_cell_vertices(m, 1));
- const unsigned int e2(GeometryInfo<dim>::line_to_cell_vertices(m, 0));
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- data.edge_sigma_values[m][q] = data.sigma_imj_values[q][e2][e1];
- data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
- }
-
- data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
- }
-
- data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
- data.edge_lambda_grads_2d[1] = {1.0, 0.0};
- data.edge_lambda_grads_2d[2] = {0.0, -1.0};
- data.edge_lambda_grads_2d[3] = {0.0, 1.0};
+ for (unsigned int m = 0; m < lines_per_cell; ++m)
+ {
+ // e1=max(reference vertex numbering on this edge)
+ // e2=min(reference vertex numbering on this edge)
+ // Which is guaranteed to be:
+ const unsigned int e1(
+ GeometryInfo<dim>::line_to_cell_vertices(m, 1));
+ const unsigned int e2(
+ GeometryInfo<dim>::line_to_cell_vertices(m, 0));
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ data.edge_sigma_values[m][q] =
+ data.sigma_imj_values[q][e2][e1];
+ data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
+ }
- // If the polynomial order is 0, then no more work to do:
- if (degree < 1)
- return;
+ data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
+ }
- // Otherwise, we can compute the non-cell dependent shape functions.
- //
- // Note: the local dof numberings follow the usual order of lines ->
- // faces -> cells
- // (we have no vertex-based DoFs in this element).
- // For a given cell we have:
- // n_line_dofs = dofs_per_line*lines_per_cell.
- // n_face_dofs = dofs_per_face*faces_per_cell.
- // n_cell_dofs = dofs_per_quad (2d)
- // = dofs_per_hex (3d)
- //
- // i.e. For the local dof numbering:
- // the first line dof is 0,
- // the first face dof is n_line_dofs,
- // the first cell dof is n_line_dofs + n_face_dofs.
- //
- // On a line, DoFs are ordered first by line_dof and then line_index:
- // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
- //
- // and similarly for faces:
- // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
- //
- // HOWEVER, we have different types of DoFs on a line/face/cell.
- // On a line we have two types, lowest order and higher order
- // gradients.
- // - The numbering is such the lowest order is first, then higher
- // order.
- // This is simple enough as there is only 1 lowest order and
- // degree higher orders DoFs per line.
- //
- // On a 2d cell, we have 3 types: Type 1/2/3:
- // - The ordering done by type:
- // - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
- // Numbered: ij1 = i1 + j1*(degree). i.e. cell_dof_index
- // = ij1.
- // - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
- // Numbered: ij2 = i2 + j2*(degree). i.e. cell_dof_index
- // = degree^2 + ij2
- // - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
- // Numbered: ij3 = i3. i.e. cell_dof_index
- // = 2*(degree^2) + ij3.
- //
- // These then fit into the local dof numbering described above:
- // - local dof numberings are:
- // line_dofs: local_dof = line_dof_index. 0 <= local_dof <
- // dofs_per_line*lines_per_cell face_dofs: local_dof =
- // n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof
- // = n_lines_dof + n_face_dofs + cell_dof_index.
- //
- // The cell-based shape functions are:
- //
- // Type 1 (gradients):
- // \phi^{C_{1}}_{ij} = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
- //
- // 0 <= i,j < degree.
- //
- // NOTE: The derivative produced by IntegratedLegendrePolynomials does
- // not account for the
- // (2*x-1) or (2*y-1) so we must take this into account when
- // taking derivatives.
- const unsigned int cell_type1_offset = n_line_dofs;
-
- // Type 2:
- // \phi^{C_{2}}_{ij} = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
- // - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
- //
- // 0 <= i,j < degree.
- const unsigned int cell_type2_offset = cell_type1_offset + degree * degree;
+ data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
+ data.edge_lambda_grads_2d[1] = {1.0, 0.0};
+ data.edge_lambda_grads_2d[2] = {0.0, -1.0};
+ data.edge_lambda_grads_2d[3] = {0.0, 1.0};
- // Type 3 (two subtypes):
- // \phi^{C_{3}}_{j} = L_{j+2}(2y-1) \mathbf{e}_{x}
- //
- // \phi^{C_{3}}_{j+degree} = L_{j+2}(2x-1) \mathbf{e}_{y}
- //
- // 0 <= j < degree
- const unsigned int cell_type3_offset1 = cell_type2_offset + degree * degree;
- const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
+ // If the polynomial order is 0, then no more work to do:
+ if (degree < 1)
+ break;
- if (flags & (update_values | update_gradients | update_hessians))
- {
- // compute all points we must evaluate the 1d polynomials at:
- std::vector<Point<dim>> cell_points(n_q_points);
- for (unsigned int q = 0; q < n_q_points; ++q)
- for (unsigned int d = 0; d < dim; ++d)
- cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
-
- // Loop through quad points:
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- // pre-compute values & required derivatives at this quad
- // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1),
+ // Otherwise, we can compute the non-cell dependent shape functions.
+ //
+ // Note: the local dof numberings follow the usual order of lines ->
+ // faces -> cells
+ // (we have no vertex-based DoFs in this element).
+ // For a given cell we have:
+ // n_line_dofs = dofs_per_line*lines_per_cell.
+ // n_face_dofs = dofs_per_face*faces_per_cell.
+ // n_cell_dofs = dofs_per_quad (2d)
+ // = dofs_per_hex (3d)
+ //
+ // i.e. For the local dof numbering:
+ // the first line dof is 0,
+ // the first face dof is n_line_dofs,
+ // the first cell dof is n_line_dofs + n_face_dofs.
+ //
+ // On a line, DoFs are ordered first by line_dof and then line_index:
+ // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
+ //
+ // and similarly for faces:
+ // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
+ //
+ // HOWEVER, we have different types of DoFs on a line/face/cell.
+ // On a line we have two types, lowest order and higher order
+ // gradients.
+ // - The numbering is such the lowest order is first, then higher
+ // order.
+ // This is simple enough as there is only 1 lowest order and
+ // degree higher orders DoFs per line.
+ //
+ // On a 2d cell, we have 3 types: Type 1/2/3:
+ // - The ordering done by type:
+ // - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
+ // Numbered: ij1 = i1 + j1*(degree). i.e. cell_dof_index
+ // = ij1.
+ // - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
+ // Numbered: ij2 = i2 + j2*(degree). i.e. cell_dof_index
+ // = degree^2 + ij2
+ // - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
+ // Numbered: ij3 = i3. i.e. cell_dof_index
+ // = 2*(degree^2) + ij3.
+ //
+ // These then fit into the local dof numbering described above:
+ // - local dof numberings are:
+ // line_dofs: local_dof = line_dof_index. 0 <= local_dof <
+ // dofs_per_line*lines_per_cell face_dofs: local_dof =
+ // n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof
+ // = n_lines_dof + n_face_dofs + cell_dof_index.
//
- // for each polyc[d], c=x,y, contains the d-th derivative with
- // respect to the coordinate c.
+ // The cell-based shape functions are:
+ //
+ // Type 1 (gradients):
+ // \phi^{C_{1}}_{ij} = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
+ //
+ // 0 <= i,j < degree.
+ //
+ // NOTE: The derivative produced by IntegratedLegendrePolynomials does
+ // not account for the
+ // (2*x-1) or (2*y-1) so we must take this into account when
+ // taking derivatives.
+ const unsigned int cell_type1_offset = n_line_dofs;
+
+ // Type 2:
+ // \phi^{C_{2}}_{ij} = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
+ // - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
+ //
+ // 0 <= i,j < degree.
+ const unsigned int cell_type2_offset =
+ cell_type1_offset + degree * degree;
- // We only need poly values and 1st derivative for
- // update_values, but need the 2nd derivative too for
- // update_gradients. For update_hessians we also need the 3rd
- // derivatives.
- const unsigned int poly_length =
- (flags & update_hessians) ? 4 :
- ((flags & update_gradients) ? 3 : 2);
+ // Type 3 (two subtypes):
+ // \phi^{C_{3}}_{j} = L_{j+2}(2y-1) \mathbf{e}_{x}
+ //
+ // \phi^{C_{3}}_{j+degree} = L_{j+2}(2x-1) \mathbf{e}_{y}
+ //
+ // 0 <= j < degree
+ const unsigned int cell_type3_offset1 =
+ cell_type2_offset + degree * degree;
+ const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
- std::vector<std::vector<double>> polyx(
- degree, std::vector<double>(poly_length));
- std::vector<std::vector<double>> polyy(
- degree, std::vector<double>(poly_length));
- for (unsigned int i = 0; i < degree; ++i)
- {
- // Compute all required 1d polynomials and their
- // derivatives, starting at degree 2. e.g. to access
- // L'_{3}(2x-1) use polyx[1][1].
- IntegratedLegendrePolynomials[i + 2].value(cell_points[q][0],
- polyx[i]);
- IntegratedLegendrePolynomials[i + 2].value(cell_points[q][1],
- polyy[i]);
- }
- // Now use these to compute the shape functions:
- if (flags & update_values)
- {
- for (unsigned int j = 0; j < degree; ++j)
- {
- const unsigned int shift_j(j * degree);
- for (unsigned int i = 0; i < degree; ++i)
- {
- const unsigned int shift_ij(i + shift_j);
-
- // Type 1:
- const unsigned int dof_index1(cell_type1_offset +
- shift_ij);
- data.shape_values[dof_index1][q][0] =
- 2.0 * polyx[i][1] * polyy[j][0];
- data.shape_values[dof_index1][q][1] =
- 2.0 * polyx[i][0] * polyy[j][1];
-
- // Type 2:
- const unsigned int dof_index2(cell_type2_offset +
- shift_ij);
- data.shape_values[dof_index2][q][0] =
- data.shape_values[dof_index1][q][0];
- data.shape_values[dof_index2][q][1] =
- -1.0 * data.shape_values[dof_index1][q][1];
- }
- // Type 3:
- const unsigned int dof_index3_1(cell_type3_offset1 + j);
- data.shape_values[dof_index3_1][q][0] = polyy[j][0];
- data.shape_values[dof_index3_1][q][1] = 0.0;
-
- const unsigned int dof_index3_2(cell_type3_offset2 + j);
- data.shape_values[dof_index3_2][q][0] = 0.0;
- data.shape_values[dof_index3_2][q][1] = polyx[j][0];
- }
- }
- if (flags & update_gradients)
+ if (flags & (update_values | update_gradients | update_hessians))
{
- for (unsigned int j = 0; j < degree; ++j)
+ // compute all points we must evaluate the 1d polynomials at:
+ std::vector<Point<dim>> cell_points(n_q_points);
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ for (unsigned int d = 0; d < dim; ++d)
+ cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
+
+ // Loop through quad points:
+ for (unsigned int q = 0; q < n_q_points; ++q)
{
- const unsigned int shift_j(j * degree);
+ // pre-compute values & required derivatives at this quad
+ // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1),
+ //
+ // for each polyc[d], c=x,y, contains the d-th derivative with
+ // respect to the coordinate c.
+
+ // We only need poly values and 1st derivative for
+ // update_values, but need the 2nd derivative too for
+ // update_gradients. For update_hessians we also need the 3rd
+ // derivatives.
+ const unsigned int poly_length =
+ (flags & update_hessians) ?
+ 4 :
+ ((flags & update_gradients) ? 3 : 2);
+
+ std::vector<std::vector<double>> polyx(
+ degree, std::vector<double>(poly_length));
+ std::vector<std::vector<double>> polyy(
+ degree, std::vector<double>(poly_length));
for (unsigned int i = 0; i < degree; ++i)
{
- const unsigned int shift_ij(i + shift_j);
-
- // Type 1:
- const unsigned int dof_index1(cell_type1_offset +
- shift_ij);
- data.shape_grads[dof_index1][q][0][0] =
- 4.0 * polyx[i][2] * polyy[j][0];
- data.shape_grads[dof_index1][q][0][1] =
- 4.0 * polyx[i][1] * polyy[j][1];
- data.shape_grads[dof_index1][q][1][0] =
- data.shape_grads[dof_index1][q][0][1];
- data.shape_grads[dof_index1][q][1][1] =
- 4.0 * polyx[i][0] * polyy[j][2];
-
- // Type 2:
- const unsigned int dof_index2(cell_type2_offset +
- shift_ij);
- data.shape_grads[dof_index2][q][0][0] =
- data.shape_grads[dof_index1][q][0][0];
- data.shape_grads[dof_index2][q][0][1] =
- data.shape_grads[dof_index1][q][0][1];
- data.shape_grads[dof_index2][q][1][0] =
- -1.0 * data.shape_grads[dof_index1][q][1][0];
- data.shape_grads[dof_index2][q][1][1] =
- -1.0 * data.shape_grads[dof_index1][q][1][1];
+ // Compute all required 1d polynomials and their
+ // derivatives, starting at degree 2. e.g. to access
+ // L'_{3}(2x-1) use polyx[1][1].
+ IntegratedLegendrePolynomials[i + 2].value(
+ cell_points[q][0], polyx[i]);
+ IntegratedLegendrePolynomials[i + 2].value(
+ cell_points[q][1], polyy[i]);
}
- // Type 3:
- const unsigned int dof_index3_1(cell_type3_offset1 + j);
- data.shape_grads[dof_index3_1][q][0][0] = 0.0;
- data.shape_grads[dof_index3_1][q][0][1] = 2.0 * polyy[j][1];
- data.shape_grads[dof_index3_1][q][1][0] = 0.0;
- data.shape_grads[dof_index3_1][q][1][1] = 0.0;
-
- const unsigned int dof_index3_2(cell_type3_offset2 + j);
- data.shape_grads[dof_index3_2][q][0][0] = 0.0;
- data.shape_grads[dof_index3_2][q][0][1] = 0.0;
- data.shape_grads[dof_index3_2][q][1][0] = 2.0 * polyx[j][1];
- data.shape_grads[dof_index3_2][q][1][1] = 0.0;
- }
- }
- if (flags & update_hessians)
- {
- for (unsigned int j = 0; j < degree; ++j)
- {
- const unsigned int shift_j(j * degree);
- for (unsigned int i = 0; i < degree; ++i)
+ // Now use these to compute the shape functions:
+ if (flags & update_values)
{
- const unsigned int shift_ij(i + shift_j);
-
- // Type 1:
- const unsigned int dof_index1(cell_type1_offset +
- shift_ij);
- data.shape_hessians[dof_index1][q][0][0][0] =
- 8.0 * polyx[i][3] * polyy[j][0];
- data.shape_hessians[dof_index1][q][1][0][0] =
- 8.0 * polyx[i][2] * polyy[j][1];
-
- data.shape_hessians[dof_index1][q][0][1][0] =
- data.shape_hessians[dof_index1][q][1][0][0];
- data.shape_hessians[dof_index1][q][1][1][0] =
- 8.0 * polyx[i][1] * polyy[j][2];
-
- data.shape_hessians[dof_index1][q][0][0][1] =
- data.shape_hessians[dof_index1][q][1][0][0];
- data.shape_hessians[dof_index1][q][1][0][1] =
- data.shape_hessians[dof_index1][q][1][1][0];
-
- data.shape_hessians[dof_index1][q][0][1][1] =
- data.shape_hessians[dof_index1][q][1][1][0];
- data.shape_hessians[dof_index1][q][1][1][1] =
- 8.0 * polyx[i][0] * polyy[j][3];
+ for (unsigned int j = 0; j < degree; ++j)
+ {
+ const unsigned int shift_j(j * degree);
+ for (unsigned int i = 0; i < degree; ++i)
+ {
+ const unsigned int shift_ij(i + shift_j);
+ // Type 1:
+ const unsigned int dof_index1(cell_type1_offset +
+ shift_ij);
+ data.shape_values[dof_index1][q][0] =
+ 2.0 * polyx[i][1] * polyy[j][0];
+ data.shape_values[dof_index1][q][1] =
+ 2.0 * polyx[i][0] * polyy[j][1];
+ // Type 2:
+ const unsigned int dof_index2(cell_type2_offset +
+ shift_ij);
+ data.shape_values[dof_index2][q][0] =
+ data.shape_values[dof_index1][q][0];
+ data.shape_values[dof_index2][q][1] =
+ -1.0 * data.shape_values[dof_index1][q][1];
+ }
+ // Type 3:
+ const unsigned int dof_index3_1(cell_type3_offset1 +
+ j);
+ data.shape_values[dof_index3_1][q][0] = polyy[j][0];
+ data.shape_values[dof_index3_1][q][1] = 0.0;
- // Type 2:
- const unsigned int dof_index2(cell_type2_offset +
- shift_ij);
- for (unsigned int d = 0; d < dim; ++d)
- {
- data.shape_hessians[dof_index2][q][0][0][d] =
- data.shape_hessians[dof_index1][q][0][0][d];
- data.shape_hessians[dof_index2][q][0][1][d] =
- data.shape_hessians[dof_index1][q][0][1][d];
- data.shape_hessians[dof_index2][q][1][0][d] =
- -1.0 * data.shape_hessians[dof_index1][q][1][0][d];
- data.shape_hessians[dof_index2][q][1][1][d] =
- -1.0 * data.shape_hessians[dof_index1][q][1][1][d];
+ const unsigned int dof_index3_2(cell_type3_offset2 +
+ j);
+ data.shape_values[dof_index3_2][q][0] = 0.0;
+ data.shape_values[dof_index3_2][q][1] = polyx[j][0];
}
}
- // Type 3:
- const unsigned int dof_index3_1(cell_type3_offset1 + j);
- data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0;
- data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0;
- data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0;
- data.shape_hessians[dof_index3_1][q][0][1][1] =
- 4.0 * polyy[j][2];
- data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0;
- data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0;
- data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0;
- data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0;
-
- const unsigned int dof_index3_2(cell_type3_offset2 + j);
- data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0;
- data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0;
- data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0;
- data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0;
- data.shape_hessians[dof_index3_2][q][1][0][0] =
- 4.0 * polyx[j][2];
- data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0;
- data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0;
- data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0;
- }
- }
- }
- }
-}
-
-
-
-template <>
-void
-FE_NedelecSZ<3, 3>::evaluate(
- const std::vector<Point<3>> p_list,
- const UpdateFlags update_flags,
- std::unique_ptr<typename dealii::FiniteElement<3, 3>::InternalDataBase>
- &data_ptr) const
-{
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
- data.update_each = requires_update_flags(update_flags);
+ if (flags & update_gradients)
+ {
+ for (unsigned int j = 0; j < degree; ++j)
+ {
+ const unsigned int shift_j(j * degree);
+ for (unsigned int i = 0; i < degree; ++i)
+ {
+ const unsigned int shift_ij(i + shift_j);
- // Useful quantities:
- const unsigned int dim = 3;
- const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1
+ // Type 1:
+ const unsigned int dof_index1(cell_type1_offset +
+ shift_ij);
+ data.shape_grads[dof_index1][q][0][0] =
+ 4.0 * polyx[i][2] * polyy[j][0];
+ data.shape_grads[dof_index1][q][0][1] =
+ 4.0 * polyx[i][1] * polyy[j][1];
+ data.shape_grads[dof_index1][q][1][0] =
+ data.shape_grads[dof_index1][q][0][1];
+ data.shape_grads[dof_index1][q][1][1] =
+ 4.0 * polyx[i][0] * polyy[j][2];
- const unsigned int vertices_per_cell =
- GeometryInfo<3 /*dim*/>::vertices_per_cell;
- const unsigned int lines_per_cell = GeometryInfo<3 /*dim*/>::lines_per_cell;
- const unsigned int faces_per_cell = GeometryInfo<3 /*dim*/>::faces_per_cell;
+ // Type 2:
+ const unsigned int dof_index2(cell_type2_offset +
+ shift_ij);
+ data.shape_grads[dof_index2][q][0][0] =
+ data.shape_grads[dof_index1][q][0][0];
+ data.shape_grads[dof_index2][q][0][1] =
+ data.shape_grads[dof_index1][q][0][1];
+ data.shape_grads[dof_index2][q][1][0] =
+ -1.0 * data.shape_grads[dof_index1][q][1][0];
+ data.shape_grads[dof_index2][q][1][1] =
+ -1.0 * data.shape_grads[dof_index1][q][1][1];
+ }
+ // Type 3:
+ const unsigned int dof_index3_1(cell_type3_offset1 +
+ j);
+ data.shape_grads[dof_index3_1][q][0][0] = 0.0;
+ data.shape_grads[dof_index3_1][q][0][1] =
+ 2.0 * polyy[j][1];
+ data.shape_grads[dof_index3_1][q][1][0] = 0.0;
+ data.shape_grads[dof_index3_1][q][1][1] = 0.0;
- const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
+ const unsigned int dof_index3_2(cell_type3_offset2 +
+ j);
+ data.shape_grads[dof_index3_2][q][0][0] = 0.0;
+ data.shape_grads[dof_index3_2][q][0][1] = 0.0;
+ data.shape_grads[dof_index3_2][q][1][0] =
+ 2.0 * polyx[j][1];
+ data.shape_grads[dof_index3_2][q][1][1] = 0.0;
+ }
+ }
+ if (flags & update_hessians)
+ {
+ for (unsigned int j = 0; j < degree; ++j)
+ {
+ const unsigned int shift_j(j * degree);
+ for (unsigned int i = 0; i < degree; ++i)
+ {
+ const unsigned int shift_ij(i + shift_j);
- // we assume that all quads have the same number of dofs
- const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
- const unsigned int n_q_points = p_list.size();
+ // Type 1:
+ const unsigned int dof_index1(cell_type1_offset +
+ shift_ij);
+ data.shape_hessians[dof_index1][q][0][0][0] =
+ 8.0 * polyx[i][3] * polyy[j][0];
+ data.shape_hessians[dof_index1][q][1][0][0] =
+ 8.0 * polyx[i][2] * polyy[j][1];
- const UpdateFlags flags(data.update_each);
+ data.shape_hessians[dof_index1][q][0][1][0] =
+ data.shape_hessians[dof_index1][q][1][0][0];
+ data.shape_hessians[dof_index1][q][1][1][0] =
+ 8.0 * polyx[i][1] * polyy[j][2];
- // Resize the internal data storage:
- data.sigma_imj_values.resize(
- n_q_points,
- std::vector<std::vector<double>>(vertices_per_cell,
- std::vector<double>(vertices_per_cell)));
+ data.shape_hessians[dof_index1][q][0][0][1] =
+ data.shape_hessians[dof_index1][q][1][0][0];
+ data.shape_hessians[dof_index1][q][1][0][1] =
+ data.shape_hessians[dof_index1][q][1][1][0];
- data.sigma_imj_grads.resize(vertices_per_cell,
- std::vector<std::vector<double>>(
- vertices_per_cell, std::vector<double>(dim)));
+ data.shape_hessians[dof_index1][q][0][1][1] =
+ data.shape_hessians[dof_index1][q][1][1][0];
+ data.shape_hessians[dof_index1][q][1][1][1] =
+ 8.0 * polyx[i][0] * polyy[j][3];
- // Resize shape function arrays according to update flags:
- if (flags & update_values)
- data.shape_values.resize(this->n_dofs_per_cell(),
- std::vector<Tensor<1, dim>>(n_q_points));
- if (flags & update_gradients)
- data.shape_grads.resize(this->n_dofs_per_cell(),
- std::vector<DerivativeForm<1, dim, dim>>(
- n_q_points));
- if (flags & update_hessians)
- data.shape_hessians.resize(this->n_dofs_per_cell(),
- std::vector<DerivativeForm<2, dim, dim>>(
- n_q_points));
-
- // Compute values of sigma & lambda and the sigma differences and
- // lambda additions.
- std::vector<std::vector<double>> sigma(n_q_points,
- std::vector<double>(lines_per_cell));
- std::vector<std::vector<double>> lambda(n_q_points,
- std::vector<double>(lines_per_cell));
+ // Type 2:
+ const unsigned int dof_index2(cell_type2_offset +
+ shift_ij);
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ data.shape_hessians[dof_index2][q][0][0][d] =
+ data.shape_hessians[dof_index1][q][0][0][d];
+ data.shape_hessians[dof_index2][q][0][1][d] =
+ data.shape_hessians[dof_index1][q][0][1][d];
+ data.shape_hessians[dof_index2][q][1][0][d] =
+ -1.0 *
+ data.shape_hessians[dof_index1][q][1][0][d];
+ data.shape_hessians[dof_index2][q][1][1][d] =
+ -1.0 *
+ data.shape_hessians[dof_index1][q][1][1][d];
+ }
+ }
+ // Type 3:
+ const unsigned int dof_index3_1(cell_type3_offset1 +
+ j);
+ data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0;
+ data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0;
+ data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0;
+ data.shape_hessians[dof_index3_1][q][0][1][1] =
+ 4.0 * polyy[j][2];
+ data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0;
+ data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0;
+ data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0;
+ data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0;
+
+ const unsigned int dof_index3_2(cell_type3_offset2 +
+ j);
+ data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0;
+ data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0;
+ data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0;
+ data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0;
+ data.shape_hessians[dof_index3_2][q][1][0][0] =
+ 4.0 * polyx[j][2];
+ data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0;
+ data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0;
+ data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0;
+ }
+ }
+ }
+ }
+ break;
+ }
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- sigma[q][0] =
- (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
- sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
- sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
- sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
- sigma[q][4] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
- sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
- sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
- sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
-
- lambda[q][0] =
- (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
- lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
- lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
- lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
- lambda[q][4] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
- lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
- lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
- lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
-
- // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
- // and lambda_ipj = \lambda_{i} + \lambda_{j}.
- for (unsigned int i = 0; i < vertices_per_cell; ++i)
- for (unsigned int j = 0; j < vertices_per_cell; ++j)
- data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
- }
+ case 3:
+ {
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) +
+ (1 - p_list[q][2]);
+ sigma[q][1] =
+ p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
+ sigma[q][2] =
+ (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
+ sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
+ sigma[q][4] =
+ (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
+ sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
+ sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
+ sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
+
+ lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) *
+ (1.0 - p_list[q][2]);
+ lambda[q][1] =
+ p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
+ lambda[q][2] =
+ (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
+ lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
+ lambda[q][4] =
+ (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
+ lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
+ lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
+ lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
+
+ // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
+ // and lambda_ipj = \lambda_{i} + \lambda_{j}.
+ for (unsigned int i = 0; i < vertices_per_cell; ++i)
+ for (unsigned int j = 0; j < vertices_per_cell; ++j)
+ data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
+ }
- // We now want some additional information about
- // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to
- // calculate values & derivatives of the shape functions we need to
- // know:
- // - The component the sigma_imj value corresponds to - this varies
- // with i & j.
- // - The gradient of the sigma_imj value
- // - this depends on the component and the direction of the
- // corresponding edge.
- // - the direction of the edge is determined by
- // sigma_imj_sign[i][j].
- //
- // Note that not every i,j combination is a valid edge (there are only
- // 12 valid edges in 3d), but we compute them all as it simplifies
- // things.
-
- // store the sign of each component x, y, z in the sigma list.
- // can use this to fill in the sigma_imj_component data.
- const int sigma_comp_signs[vertices_per_cell][dim] = {{-1, -1, -1},
- {1, -1, -1},
- {-1, 1, -1},
- {1, 1, -1},
- {-1, -1, 1},
- {1, -1, 1},
- {-1, 1, 1},
- {1, 1, 1}};
-
- int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
- unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell];
-
- for (unsigned int i = 0; i < vertices_per_cell; ++i)
- for (unsigned int j = 0; j < vertices_per_cell; ++j)
- {
- // sigma_imj_sign is the sign (+/-) of the coefficient of
- // x/y/z in sigma_imj. Due to the numbering of vertices on the
- // reference element this is easy to work out because edges in
- // the positive direction go from smaller to higher local
- // vertex numbering.
- sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
- sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
-
- // Now store the component which the sigma_i - sigma_j
- // corresponds to:
- sigma_imj_component[i][j] = 0;
- for (unsigned int d = 0; d < dim; ++d)
- {
- int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
- // Only interested in the first non-zero
- // as if there is a second, it will not be a valid edge.
- if (temp_imj != 0)
+ // We now want some additional information about
+ // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to
+ // calculate values & derivatives of the shape functions we need to
+ // know:
+ // - The component the sigma_imj value corresponds to - this varies
+ // with i & j.
+ // - The gradient of the sigma_imj value
+ // - this depends on the component and the direction of the
+ // corresponding edge.
+ // - the direction of the edge is determined by
+ // sigma_imj_sign[i][j].
+ //
+ // Note that not every i,j combination is a valid edge (there are only
+ // 12 valid edges in 3d), but we compute them all as it simplifies
+ // things.
+
+ // store the sign of each component x, y, z in the sigma list.
+ // can use this to fill in the sigma_imj_component data.
+ const int sigma_comp_signs[GeometryInfo<3>::vertices_per_cell][3] = {
+ {-1, -1, -1},
+ {1, -1, -1},
+ {-1, 1, -1},
+ {1, 1, -1},
+ {-1, -1, 1},
+ {1, -1, 1},
+ {-1, 1, 1},
+ {1, 1, 1}};
+
+ int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
+ unsigned int sigma_imj_component[vertices_per_cell]
+ [vertices_per_cell];
+
+ for (unsigned int i = 0; i < vertices_per_cell; ++i)
+ for (unsigned int j = 0; j < vertices_per_cell; ++j)
{
- sigma_imj_component[i][j] = d;
- break;
+ // sigma_imj_sign is the sign (+/-) of the coefficient of
+ // x/y/z in sigma_imj. Due to the numbering of vertices on the
+ // reference element this is easy to work out because edges in
+ // the positive direction go from smaller to higher local
+ // vertex numbering.
+ sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
+ sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
+
+ // Now store the component which the sigma_i - sigma_j
+ // corresponds to:
+ sigma_imj_component[i][j] = 0;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ int temp_imj =
+ sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
+ // Only interested in the first non-zero
+ // as if there is a second, it will not be a valid edge.
+ if (temp_imj != 0)
+ {
+ sigma_imj_component[i][j] = d;
+ break;
+ }
+ }
+ // Can now calculate the gradient, only non-zero in the
+ // component given: Note some i,j combinations will be
+ // incorrect, but only on invalid edges.
+ data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
+ 2.0 * sigma_imj_sign[i][j];
}
- }
- // Can now calculate the gradient, only non-zero in the
- // component given: Note some i,j combinations will be
- // incorrect, but only on invalid edges.
- data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
- 2.0 * sigma_imj_sign[i][j];
- }
-
- // Now compute the edge parameterisations for a single element
- // with global numbering matching that of the reference element:
-
- // resize the edge parameterisations
- data.edge_sigma_values.resize(lines_per_cell,
- std::vector<double>(n_q_points));
- data.edge_lambda_values.resize(lines_per_cell,
- std::vector<double>(n_q_points));
- data.edge_sigma_grads.resize(lines_per_cell, std::vector<double>(dim));
- data.edge_lambda_grads_3d.resize(
- lines_per_cell,
- std::vector<std::vector<double>>(n_q_points, std::vector<double>(dim)));
- data.edge_lambda_gradgrads_3d.resize(
- lines_per_cell,
- std::vector<std::vector<double>>(dim, std::vector<double>(dim)));
-
- // Fill the values:
- const unsigned int edge_sigma_direction[lines_per_cell] = {
- 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
-
- for (unsigned int m = 0; m < lines_per_cell; ++m)
- {
- // e1=max(reference vertex numbering on this edge)
- // e2=min(reference vertex numbering on this edge)
- // Which is guaranteed to be:
- const unsigned int e1(GeometryInfo<dim>::line_to_cell_vertices(m, 1));
- const unsigned int e2(GeometryInfo<dim>::line_to_cell_vertices(m, 0));
-
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- data.edge_sigma_values[m][q] = data.sigma_imj_values[q][e2][e1];
- data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
- }
- data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
- }
+ // Now compute the edge parameterisations for a single element
+ // with global numbering matching that of the reference element:
+
+ // resize the edge parameterisations
+ data.edge_sigma_values.resize(lines_per_cell,
+ std::vector<double>(n_q_points));
+ data.edge_lambda_values.resize(lines_per_cell,
+ std::vector<double>(n_q_points));
+ data.edge_sigma_grads.resize(lines_per_cell,
+ std::vector<double>(dim));
+ data.edge_lambda_grads_3d.resize(
+ lines_per_cell,
+ std::vector<std::vector<double>>(n_q_points,
+ std::vector<double>(dim)));
+ data.edge_lambda_gradgrads_3d.resize(
+ lines_per_cell,
+ std::vector<std::vector<double>>(dim, std::vector<double>(dim)));
+
+ // Fill the values:
+ const unsigned int
+ edge_sigma_direction[GeometryInfo<3>::lines_per_cell] = {
+ 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
- // edge_lambda_grads
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- double x(p_list[q][0]);
- double y(p_list[q][1]);
- double z(p_list[q][2]);
- data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0};
- data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x};
- data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0};
- data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y};
- data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x};
- data.edge_lambda_grads_3d[5][q] = {z, 0.0, x};
- data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y};
- data.edge_lambda_grads_3d[7][q] = {0.0, z, y};
- data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0};
- data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0};
- data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
- data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
- }
+ for (unsigned int m = 0; m < lines_per_cell; ++m)
+ {
+ // e1=max(reference vertex numbering on this edge)
+ // e2=min(reference vertex numbering on this edge)
+ // Which is guaranteed to be:
+ const unsigned int e1(
+ GeometryInfo<dim>::line_to_cell_vertices(m, 1));
+ const unsigned int e2(
+ GeometryInfo<dim>::line_to_cell_vertices(m, 0));
- // edge_lambda gradgrads:
- const int edge_lambda_sign[lines_per_cell] = {
- 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1}; // sign of the 2nd derivative for
- // each edge.
-
- const unsigned int edge_lambda_directions[lines_per_cell][2] = {
- {0, 2},
- {0, 2},
- {1, 2},
- {1, 2},
- {0, 2},
- {0, 2},
- {1, 2},
- {1, 2},
- {0, 1},
- {0, 1},
- {0, 1},
- {0, 1}}; // component which edge_lambda[m] depends on.
-
- for (unsigned int m = 0; m < lines_per_cell; ++m)
- {
- data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
- [edge_lambda_directions[m][1]] =
- edge_lambda_sign[m];
- data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
- [edge_lambda_directions[m][0]] =
- edge_lambda_sign[m];
- }
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ data.edge_sigma_values[m][q] =
+ data.sigma_imj_values[q][e2][e1];
+ data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
+ }
- // If the polynomial order is 0, then no more work to do:
- if (degree < 1)
- return;
+ data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
+ }
- // resize required data:
- data.face_lambda_values.resize(faces_per_cell,
- std::vector<double>(n_q_points));
- data.face_lambda_grads.resize(faces_per_cell, std::vector<double>(dim));
+ // edge_lambda_grads
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ double x(p_list[q][0]);
+ double y(p_list[q][1]);
+ double z(p_list[q][2]);
+ data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0};
+ data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x};
+ data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0};
+ data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y};
+ data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x};
+ data.edge_lambda_grads_3d[5][q] = {z, 0.0, x};
+ data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y};
+ data.edge_lambda_grads_3d[7][q] = {0.0, z, y};
+ data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0};
+ data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0};
+ data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
+ data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
+ }
- // Fill in the values (these don't change between cells).
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- double x(p_list[q][0]);
- double y(p_list[q][1]);
- double z(p_list[q][2]);
- data.face_lambda_values[0][q] = 1.0 - x;
- data.face_lambda_values[1][q] = x;
- data.face_lambda_values[2][q] = 1.0 - y;
- data.face_lambda_values[3][q] = y;
- data.face_lambda_values[4][q] = 1.0 - z;
- data.face_lambda_values[5][q] = z;
- }
+ // edge_lambda gradgrads:
+ const int edge_lambda_sign[GeometryInfo<3>::lines_per_cell] = {
+ 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1}; // sign of the 2nd
+ // derivative for each
+ // edge.
- // gradients are constant:
- data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
- data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
- data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
- data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
- data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
- data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
-
- // for cell-based shape functions:
- // these don't depend on the cell, so can precompute all here:
- if (flags & (update_values | update_gradients | update_hessians))
- {
- // Cell-based shape functions:
- //
- // Type-1 (gradients):
- // \phi^{C_{1}}_{ijk} = grad(
- // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
- //
- // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
- const unsigned int cell_type1_offset(n_line_dofs + n_face_dofs);
- // Type-2:
- //
- // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
- // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1,
- // -1)\phi^{C_{1}}_{ijk}
- //
- // 0 <= i,j,k < degree. (subtypes in groups of
- // degree*degree*degree)
- //
- // here we order so that all of subtype 1 comes first, then
- // subtype 2.
- const unsigned int cell_type2_offset1(cell_type1_offset +
- degree * degree * degree);
- const unsigned int cell_type2_offset2(cell_type2_offset1 +
- degree * degree * degree);
- // Type-3
- // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
- // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
- // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
- //
- // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
- //
- // again we order so we compute all of subtype 1 first, then
- // subtype 2, etc.
- const unsigned int cell_type3_offset1(cell_type2_offset2 +
- degree * degree * degree);
- const unsigned int cell_type3_offset2(cell_type3_offset1 +
- degree * degree);
- const unsigned int cell_type3_offset3(cell_type3_offset2 +
- degree * degree);
-
- // compute all points we must evaluate the 1d polynomials at:
- std::vector<Point<dim>> cell_points(n_q_points);
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- for (unsigned int d = 0; d < dim; ++d)
+ const unsigned int
+ edge_lambda_directions[GeometryInfo<3>::lines_per_cell][2] = {
+ {0, 2},
+ {0, 2},
+ {1, 2},
+ {1, 2},
+ {0, 2},
+ {0, 2},
+ {1, 2},
+ {1, 2},
+ {0, 1},
+ {0, 1},
+ {0, 1},
+ {0, 1}}; // component which edge_lambda[m] depends on.
+
+ for (unsigned int m = 0; m < lines_per_cell; ++m)
{
- cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
+ data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
+ [edge_lambda_directions[m][1]] =
+ edge_lambda_sign[m];
+ data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
+ [edge_lambda_directions[m][0]] =
+ edge_lambda_sign[m];
}
- }
- // We only need poly values and 1st derivative for
- // update_values, but need the 2nd derivative too for
- // update_gradients. For update_hessians we also need 3rd
- // derivative.
- const unsigned int poly_length =
- (flags & update_hessians) ? 4 : ((flags & update_gradients) ? 3 : 2);
+ // If the polynomial order is 0, then no more work to do:
+ if (degree < 1)
+ break;
- // Loop through quad points:
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- // pre-compute values & required derivatives at this quad
- // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy =
- // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1).
- //
- // for each polyc[d], c=x,y,z, contains the d-th
- // derivative with respect to the coordinate c.
- std::vector<std::vector<double>> polyx(
- degree, std::vector<double>(poly_length));
- std::vector<std::vector<double>> polyy(
- degree, std::vector<double>(poly_length));
- std::vector<std::vector<double>> polyz(
- degree, std::vector<double>(poly_length));
- for (unsigned int i = 0; i < degree; ++i)
+ // resize required data:
+ data.face_lambda_values.resize(faces_per_cell,
+ std::vector<double>(n_q_points));
+ data.face_lambda_grads.resize(faces_per_cell,
+ std::vector<double>(dim));
+
+ // Fill in the values (these don't change between cells).
+ for (unsigned int q = 0; q < n_q_points; ++q)
{
- // compute all required 1d polynomials for i
- IntegratedLegendrePolynomials[i + 2].value(cell_points[q][0],
- polyx[i]);
- IntegratedLegendrePolynomials[i + 2].value(cell_points[q][1],
- polyy[i]);
- IntegratedLegendrePolynomials[i + 2].value(cell_points[q][2],
- polyz[i]);
+ double x(p_list[q][0]);
+ double y(p_list[q][1]);
+ double z(p_list[q][2]);
+ data.face_lambda_values[0][q] = 1.0 - x;
+ data.face_lambda_values[1][q] = x;
+ data.face_lambda_values[2][q] = 1.0 - y;
+ data.face_lambda_values[3][q] = y;
+ data.face_lambda_values[4][q] = 1.0 - z;
+ data.face_lambda_values[5][q] = z;
}
- // Now use these to compute the shape functions:
- if (flags & update_values)
+
+ // gradients are constant:
+ data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
+ data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
+ data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
+ data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
+ data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
+ data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
+
+ // for cell-based shape functions:
+ // these don't depend on the cell, so can precompute all here:
+ if (flags & (update_values | update_gradients | update_hessians))
{
- for (unsigned int k = 0; k < degree; ++k)
+ // Cell-based shape functions:
+ //
+ // Type-1 (gradients):
+ // \phi^{C_{1}}_{ijk} = grad(
+ // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
+ //
+ // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
+ const unsigned int cell_type1_offset(n_line_dofs + n_face_dofs);
+ // Type-2:
+ //
+ // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
+ // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1,
+ // -1)\phi^{C_{1}}_{ijk}
+ //
+ // 0 <= i,j,k < degree. (subtypes in groups of
+ // degree*degree*degree)
+ //
+ // here we order so that all of subtype 1 comes first, then
+ // subtype 2.
+ const unsigned int cell_type2_offset1(cell_type1_offset +
+ degree * degree * degree);
+ const unsigned int cell_type2_offset2(cell_type2_offset1 +
+ degree * degree * degree);
+ // Type-3
+ // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
+ // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
+ // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
+ //
+ // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
+ //
+ // again we order so we compute all of subtype 1 first, then
+ // subtype 2, etc.
+ const unsigned int cell_type3_offset1(cell_type2_offset2 +
+ degree * degree * degree);
+ const unsigned int cell_type3_offset2(cell_type3_offset1 +
+ degree * degree);
+ const unsigned int cell_type3_offset3(cell_type3_offset2 +
+ degree * degree);
+
+ // compute all points we must evaluate the 1d polynomials at:
+ std::vector<Point<dim>> cell_points(n_q_points);
+ for (unsigned int q = 0; q < n_q_points; ++q)
{
- const unsigned int shift_k(k * degree * degree);
- const unsigned int shift_j(k *
- degree); // Used below when subbing
- // k for j (type 3)
- for (unsigned int j = 0; j < degree; ++j)
+ for (unsigned int d = 0; d < dim; ++d)
{
- const unsigned int shift_jk(j * degree + shift_k);
- for (unsigned int i = 0; i < degree; ++i)
- {
- const unsigned int shift_ijk(shift_jk + i);
-
- // Type 1:
- const unsigned int dof_index1(cell_type1_offset +
- shift_ijk);
-
- data.shape_values[dof_index1][q][0] =
- 2.0 * polyx[i][1] * polyy[j][0] * polyz[k][0];
- data.shape_values[dof_index1][q][1] =
- 2.0 * polyx[i][0] * polyy[j][1] * polyz[k][0];
- data.shape_values[dof_index1][q][2] =
- 2.0 * polyx[i][0] * polyy[j][0] * polyz[k][1];
-
- // Type 2:
- const unsigned int dof_index2_1(cell_type2_offset1 +
- shift_ijk);
- const unsigned int dof_index2_2(cell_type2_offset2 +
- shift_ijk);
-
- data.shape_values[dof_index2_1][q][0] =
- data.shape_values[dof_index1][q][0];
- data.shape_values[dof_index2_1][q][1] =
- -1.0 * data.shape_values[dof_index1][q][1];
- data.shape_values[dof_index2_1][q][2] =
- data.shape_values[dof_index1][q][2];
-
- data.shape_values[dof_index2_2][q][0] =
- data.shape_values[dof_index1][q][0];
- data.shape_values[dof_index2_2][q][1] =
- -1.0 * data.shape_values[dof_index1][q][1];
- data.shape_values[dof_index2_2][q][2] =
- -1.0 * data.shape_values[dof_index1][q][2];
- }
- // Type 3: (note we re-use k and j for
- // convenience):
- const unsigned int shift_ij(j +
- shift_j); // here we've subbed
- // j for i, k for j.
- const unsigned int dof_index3_1(cell_type3_offset1 +
- shift_ij);
- const unsigned int dof_index3_2(cell_type3_offset2 +
- shift_ij);
- const unsigned int dof_index3_3(cell_type3_offset3 +
- shift_ij);
-
- data.shape_values[dof_index3_1][q][0] =
- polyy[j][0] * polyz[k][0];
- data.shape_values[dof_index3_1][q][1] = 0.0;
- data.shape_values[dof_index3_1][q][2] = 0.0;
-
- data.shape_values[dof_index3_2][q][0] = 0.0;
- data.shape_values[dof_index3_2][q][1] =
- polyx[j][0] * polyz[k][0];
- data.shape_values[dof_index3_2][q][2] = 0.0;
-
- data.shape_values[dof_index3_3][q][0] = 0.0;
- data.shape_values[dof_index3_3][q][1] = 0.0;
- data.shape_values[dof_index3_3][q][2] =
- polyx[j][0] * polyy[k][0];
+ cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
}
}
- }
- if (flags & update_gradients)
- {
- for (unsigned int k = 0; k < degree; ++k)
+
+ // We only need poly values and 1st derivative for
+ // update_values, but need the 2nd derivative too for
+ // update_gradients. For update_hessians we also need 3rd
+ // derivative.
+ const unsigned int poly_length =
+ (flags & update_hessians) ?
+ 4 :
+ ((flags & update_gradients) ? 3 : 2);
+
+ // Loop through quad points:
+ for (unsigned int q = 0; q < n_q_points; ++q)
{
- const unsigned int shift_k(k * degree * degree);
- const unsigned int shift_j(k *
- degree); // Used below when subbing
- // k for j (type 3)
- for (unsigned int j = 0; j < degree; ++j)
+ // pre-compute values & required derivatives at this quad
+ // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy =
+ // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1).
+ //
+ // for each polyc[d], c=x,y,z, contains the d-th
+ // derivative with respect to the coordinate c.
+ std::vector<std::vector<double>> polyx(
+ degree, std::vector<double>(poly_length));
+ std::vector<std::vector<double>> polyy(
+ degree, std::vector<double>(poly_length));
+ std::vector<std::vector<double>> polyz(
+ degree, std::vector<double>(poly_length));
+ for (unsigned int i = 0; i < degree; ++i)
{
- const unsigned int shift_jk(j * degree + shift_k);
- for (unsigned int i = 0; i < degree; ++i)
+ // compute all required 1d polynomials for i
+ IntegratedLegendrePolynomials[i + 2].value(
+ cell_points[q][0], polyx[i]);
+ IntegratedLegendrePolynomials[i + 2].value(
+ cell_points[q][1], polyy[i]);
+ IntegratedLegendrePolynomials[i + 2].value(
+ cell_points[q][2], polyz[i]);
+ }
+ // Now use these to compute the shape functions:
+ if (flags & update_values)
+ {
+ for (unsigned int k = 0; k < degree; ++k)
{
- const unsigned int shift_ijk(shift_jk + i);
-
- // Type 1:
- const unsigned int dof_index1(cell_type1_offset +
- shift_ijk);
-
- data.shape_grads[dof_index1][q][0][0] =
- 4.0 * polyx[i][2] * polyy[j][0] * polyz[k][0];
- data.shape_grads[dof_index1][q][0][1] =
- 4.0 * polyx[i][1] * polyy[j][1] * polyz[k][0];
- data.shape_grads[dof_index1][q][0][2] =
- 4.0 * polyx[i][1] * polyy[j][0] * polyz[k][1];
-
- data.shape_grads[dof_index1][q][1][0] =
- data.shape_grads[dof_index1][q][0][1];
- data.shape_grads[dof_index1][q][1][1] =
- 4.0 * polyx[i][0] * polyy[j][2] * polyz[k][0];
- data.shape_grads[dof_index1][q][1][2] =
- 4.0 * polyx[i][0] * polyy[j][1] * polyz[k][1];
-
- data.shape_grads[dof_index1][q][2][0] =
- data.shape_grads[dof_index1][q][0][2];
- data.shape_grads[dof_index1][q][2][1] =
- data.shape_grads[dof_index1][q][1][2];
- data.shape_grads[dof_index1][q][2][2] =
- 4.0 * polyx[i][0] * polyy[j][0] * polyz[k][2];
-
- // Type 2:
- const unsigned int dof_index2_1(cell_type2_offset1 +
- shift_ijk);
- const unsigned int dof_index2_2(cell_type2_offset2 +
- shift_ijk);
-
- for (unsigned int d = 0; d < dim; ++d)
+ const unsigned int shift_k(k * degree * degree);
+ const unsigned int shift_j(
+ k * degree); // Used below when subbing
+ // k for j (type 3)
+ for (unsigned int j = 0; j < degree; ++j)
{
- data.shape_grads[dof_index2_1][q][0][d] =
- data.shape_grads[dof_index1][q][0][d];
- data.shape_grads[dof_index2_1][q][1][d] =
- -1.0 * data.shape_grads[dof_index1][q][1][d];
- data.shape_grads[dof_index2_1][q][2][d] =
- data.shape_grads[dof_index1][q][2][d];
-
- data.shape_grads[dof_index2_2][q][0][d] =
- data.shape_grads[dof_index1][q][0][d];
- data.shape_grads[dof_index2_2][q][1][d] =
- -1.0 * data.shape_grads[dof_index1][q][1][d];
- data.shape_grads[dof_index2_2][q][2][d] =
- -1.0 * data.shape_grads[dof_index1][q][2][d];
+ const unsigned int shift_jk(j * degree + shift_k);
+ for (unsigned int i = 0; i < degree; ++i)
+ {
+ const unsigned int shift_ijk(shift_jk + i);
+
+ // Type 1:
+ const unsigned int dof_index1(
+ cell_type1_offset + shift_ijk);
+
+ data.shape_values[dof_index1][q][0] =
+ 2.0 * polyx[i][1] * polyy[j][0] *
+ polyz[k][0];
+ data.shape_values[dof_index1][q][1] =
+ 2.0 * polyx[i][0] * polyy[j][1] *
+ polyz[k][0];
+ data.shape_values[dof_index1][q][2] =
+ 2.0 * polyx[i][0] * polyy[j][0] *
+ polyz[k][1];
+
+ // Type 2:
+ const unsigned int dof_index2_1(
+ cell_type2_offset1 + shift_ijk);
+ const unsigned int dof_index2_2(
+ cell_type2_offset2 + shift_ijk);
+
+ data.shape_values[dof_index2_1][q][0] =
+ data.shape_values[dof_index1][q][0];
+ data.shape_values[dof_index2_1][q][1] =
+ -1.0 * data.shape_values[dof_index1][q][1];
+ data.shape_values[dof_index2_1][q][2] =
+ data.shape_values[dof_index1][q][2];
+
+ data.shape_values[dof_index2_2][q][0] =
+ data.shape_values[dof_index1][q][0];
+ data.shape_values[dof_index2_2][q][1] =
+ -1.0 * data.shape_values[dof_index1][q][1];
+ data.shape_values[dof_index2_2][q][2] =
+ -1.0 * data.shape_values[dof_index1][q][2];
+ }
+ // Type 3: (note we re-use k and j for
+ // convenience):
+ const unsigned int shift_ij(
+ j + shift_j); // here we've subbed
+ // j for i, k for j.
+ const unsigned int dof_index3_1(
+ cell_type3_offset1 + shift_ij);
+ const unsigned int dof_index3_2(
+ cell_type3_offset2 + shift_ij);
+ const unsigned int dof_index3_3(
+ cell_type3_offset3 + shift_ij);
+
+ data.shape_values[dof_index3_1][q][0] =
+ polyy[j][0] * polyz[k][0];
+ data.shape_values[dof_index3_1][q][1] = 0.0;
+ data.shape_values[dof_index3_1][q][2] = 0.0;
+
+ data.shape_values[dof_index3_2][q][0] = 0.0;
+ data.shape_values[dof_index3_2][q][1] =
+ polyx[j][0] * polyz[k][0];
+ data.shape_values[dof_index3_2][q][2] = 0.0;
+
+ data.shape_values[dof_index3_3][q][0] = 0.0;
+ data.shape_values[dof_index3_3][q][1] = 0.0;
+ data.shape_values[dof_index3_3][q][2] =
+ polyx[j][0] * polyy[k][0];
}
}
- // Type 3: (note we re-use k and j for
- // convenience):
- const unsigned int shift_ij(j +
- shift_j); // here we've subbed
- // j for i, k for j.
- const unsigned int dof_index3_1(cell_type3_offset1 +
- shift_ij);
- const unsigned int dof_index3_2(cell_type3_offset2 +
- shift_ij);
- const unsigned int dof_index3_3(cell_type3_offset3 +
- shift_ij);
- for (unsigned int d1 = 0; d1 < dim; ++d1)
+ }
+ if (flags & update_gradients)
+ {
+ for (unsigned int k = 0; k < degree; ++k)
{
- for (unsigned int d2 = 0; d2 < dim; ++d2)
+ const unsigned int shift_k(k * degree * degree);
+ const unsigned int shift_j(
+ k * degree); // Used below when subbing
+ // k for j (type 3)
+ for (unsigned int j = 0; j < degree; ++j)
{
- data.shape_grads[dof_index3_1][q][d1][d2] = 0.0;
- data.shape_grads[dof_index3_2][q][d1][d2] = 0.0;
- data.shape_grads[dof_index3_3][q][d1][d2] = 0.0;
+ const unsigned int shift_jk(j * degree + shift_k);
+ for (unsigned int i = 0; i < degree; ++i)
+ {
+ const unsigned int shift_ijk(shift_jk + i);
+
+ // Type 1:
+ const unsigned int dof_index1(
+ cell_type1_offset + shift_ijk);
+
+ data.shape_grads[dof_index1][q][0][0] =
+ 4.0 * polyx[i][2] * polyy[j][0] *
+ polyz[k][0];
+ data.shape_grads[dof_index1][q][0][1] =
+ 4.0 * polyx[i][1] * polyy[j][1] *
+ polyz[k][0];
+ data.shape_grads[dof_index1][q][0][2] =
+ 4.0 * polyx[i][1] * polyy[j][0] *
+ polyz[k][1];
+
+ data.shape_grads[dof_index1][q][1][0] =
+ data.shape_grads[dof_index1][q][0][1];
+ data.shape_grads[dof_index1][q][1][1] =
+ 4.0 * polyx[i][0] * polyy[j][2] *
+ polyz[k][0];
+ data.shape_grads[dof_index1][q][1][2] =
+ 4.0 * polyx[i][0] * polyy[j][1] *
+ polyz[k][1];
+
+ data.shape_grads[dof_index1][q][2][0] =
+ data.shape_grads[dof_index1][q][0][2];
+ data.shape_grads[dof_index1][q][2][1] =
+ data.shape_grads[dof_index1][q][1][2];
+ data.shape_grads[dof_index1][q][2][2] =
+ 4.0 * polyx[i][0] * polyy[j][0] *
+ polyz[k][2];
+
+ // Type 2:
+ const unsigned int dof_index2_1(
+ cell_type2_offset1 + shift_ijk);
+ const unsigned int dof_index2_2(
+ cell_type2_offset2 + shift_ijk);
+
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ data.shape_grads[dof_index2_1][q][0][d] =
+ data.shape_grads[dof_index1][q][0][d];
+ data.shape_grads[dof_index2_1][q][1][d] =
+ -1.0 *
+ data.shape_grads[dof_index1][q][1][d];
+ data.shape_grads[dof_index2_1][q][2][d] =
+ data.shape_grads[dof_index1][q][2][d];
+
+ data.shape_grads[dof_index2_2][q][0][d] =
+ data.shape_grads[dof_index1][q][0][d];
+ data.shape_grads[dof_index2_2][q][1][d] =
+ -1.0 *
+ data.shape_grads[dof_index1][q][1][d];
+ data.shape_grads[dof_index2_2][q][2][d] =
+ -1.0 *
+ data.shape_grads[dof_index1][q][2][d];
+ }
+ }
+ // Type 3: (note we re-use k and j for
+ // convenience):
+ const unsigned int shift_ij(
+ j + shift_j); // here we've subbed
+ // j for i, k for j.
+ const unsigned int dof_index3_1(
+ cell_type3_offset1 + shift_ij);
+ const unsigned int dof_index3_2(
+ cell_type3_offset2 + shift_ij);
+ const unsigned int dof_index3_3(
+ cell_type3_offset3 + shift_ij);
+ for (unsigned int d1 = 0; d1 < dim; ++d1)
+ {
+ for (unsigned int d2 = 0; d2 < dim; ++d2)
+ {
+ data
+ .shape_grads[dof_index3_1][q][d1][d2] =
+ 0.0;
+ data
+ .shape_grads[dof_index3_2][q][d1][d2] =
+ 0.0;
+ data
+ .shape_grads[dof_index3_3][q][d1][d2] =
+ 0.0;
+ }
+ }
+ data.shape_grads[dof_index3_1][q][0][1] =
+ 2.0 * polyy[j][1] * polyz[k][0];
+ data.shape_grads[dof_index3_1][q][0][2] =
+ 2.0 * polyy[j][0] * polyz[k][1];
+
+ data.shape_grads[dof_index3_2][q][1][0] =
+ 2.0 * polyx[j][1] * polyz[k][0];
+ data.shape_grads[dof_index3_2][q][1][2] =
+ 2.0 * polyx[j][0] * polyz[k][1];
+
+ data.shape_grads[dof_index3_3][q][2][0] =
+ 2.0 * polyx[j][1] * polyy[k][0];
+ data.shape_grads[dof_index3_3][q][2][1] =
+ 2.0 * polyx[j][0] * polyy[k][1];
}
}
- data.shape_grads[dof_index3_1][q][0][1] =
- 2.0 * polyy[j][1] * polyz[k][0];
- data.shape_grads[dof_index3_1][q][0][2] =
- 2.0 * polyy[j][0] * polyz[k][1];
-
- data.shape_grads[dof_index3_2][q][1][0] =
- 2.0 * polyx[j][1] * polyz[k][0];
- data.shape_grads[dof_index3_2][q][1][2] =
- 2.0 * polyx[j][0] * polyz[k][1];
-
- data.shape_grads[dof_index3_3][q][2][0] =
- 2.0 * polyx[j][1] * polyy[k][0];
- data.shape_grads[dof_index3_3][q][2][1] =
- 2.0 * polyx[j][0] * polyy[k][1];
}
- }
- }
- if (flags & update_hessians)
- {
- for (unsigned int k = 0; k < degree; ++k)
- {
- const unsigned int shift_k(k * degree * degree);
- const unsigned int shift_j(k *
- degree); // Used below when subbing
- // k for j type 3
-
- for (unsigned int j = 0; j < degree; ++j)
+ if (flags & update_hessians)
{
- const unsigned int shift_jk(j * degree + shift_k);
- for (unsigned int i = 0; i < degree; ++i)
+ for (unsigned int k = 0; k < degree; ++k)
{
- const unsigned int shift_ijk(shift_jk + i);
-
- // Type 1:
- const unsigned int dof_index1(cell_type1_offset +
- shift_ijk);
-
- data.shape_hessians[dof_index1][q][0][0][0] =
- 8.0 * polyx[i][3] * polyy[j][0] * polyz[k][0];
- data.shape_hessians[dof_index1][q][1][0][0] =
- 8.0 * polyx[i][2] * polyy[j][1] * polyz[k][0];
- data.shape_hessians[dof_index1][q][2][0][0] =
- 8.0 * polyx[i][2] * polyy[j][0] * polyz[k][1];
-
- data.shape_hessians[dof_index1][q][0][1][0] =
- data.shape_hessians[dof_index1][q][1][0][0];
- data.shape_hessians[dof_index1][q][1][1][0] =
- 8.0 * polyx[i][1] * polyy[j][2] * polyz[k][0];
- data.shape_hessians[dof_index1][q][2][1][0] =
- 8.0 * polyx[i][1] * polyy[j][1] * polyz[k][1];
-
- data.shape_hessians[dof_index1][q][0][2][0] =
- data.shape_hessians[dof_index1][q][2][0][0];
- data.shape_hessians[dof_index1][q][1][2][0] =
- data.shape_hessians[dof_index1][q][2][1][0];
- data.shape_hessians[dof_index1][q][2][2][0] =
- 8.0 * polyx[i][1] * polyy[j][0] * polyz[k][2];
-
-
- data.shape_hessians[dof_index1][q][0][0][1] =
- data.shape_hessians[dof_index1][q][1][0][0];
- data.shape_hessians[dof_index1][q][1][0][1] =
- data.shape_hessians[dof_index1][q][1][1][0];
- data.shape_hessians[dof_index1][q][2][0][1] =
- data.shape_hessians[dof_index1][q][2][1][0];
-
- data.shape_hessians[dof_index1][q][0][1][1] =
- data.shape_hessians[dof_index1][q][1][1][0];
- data.shape_hessians[dof_index1][q][1][1][1] =
- 8.0 * polyx[i][0] * polyy[j][3] * polyz[k][0];
- data.shape_hessians[dof_index1][q][2][1][1] =
- 8.0 * polyx[i][0] * polyy[j][2] * polyz[k][1];
-
- data.shape_hessians[dof_index1][q][0][2][1] =
- data.shape_hessians[dof_index1][q][2][1][0];
- data.shape_hessians[dof_index1][q][1][2][1] =
- data.shape_hessians[dof_index1][q][2][1][1];
- data.shape_hessians[dof_index1][q][2][2][1] =
- 8.0 * polyx[i][0] * polyy[j][1] * polyz[k][2];
-
-
- data.shape_hessians[dof_index1][q][0][0][2] =
- data.shape_hessians[dof_index1][q][2][0][0];
- data.shape_hessians[dof_index1][q][1][0][2] =
- data.shape_hessians[dof_index1][q][2][1][0];
- data.shape_hessians[dof_index1][q][2][0][2] =
- data.shape_hessians[dof_index1][q][2][2][0];
-
- data.shape_hessians[dof_index1][q][0][1][2] =
- data.shape_hessians[dof_index1][q][2][1][0];
- data.shape_hessians[dof_index1][q][1][1][2] =
- data.shape_hessians[dof_index1][q][2][1][1];
- data.shape_hessians[dof_index1][q][2][1][2] =
- data.shape_hessians[dof_index1][q][2][2][1];
-
- data.shape_hessians[dof_index1][q][0][2][2] =
- data.shape_hessians[dof_index1][q][2][2][0];
- data.shape_hessians[dof_index1][q][1][2][2] =
- data.shape_hessians[dof_index1][q][2][2][1];
- data.shape_hessians[dof_index1][q][2][2][2] =
- 8.0 * polyx[i][0] * polyy[j][0] * polyz[k][3];
-
-
- // Type 2:
- const unsigned int dof_index2_1(cell_type2_offset1 +
- shift_ijk);
- const unsigned int dof_index2_2(cell_type2_offset2 +
- shift_ijk);
+ const unsigned int shift_k(k * degree * degree);
+ const unsigned int shift_j(
+ k * degree); // Used below when subbing
+ // k for j type 3
- for (unsigned int d1 = 0; d1 < dim; ++d1)
+ for (unsigned int j = 0; j < degree; ++j)
{
- for (unsigned int d2 = 0; d2 < dim; ++d2)
+ const unsigned int shift_jk(j * degree + shift_k);
+ for (unsigned int i = 0; i < degree; ++i)
{
- data.shape_hessians[dof_index2_1][q][0][d1]
- [d2] =
- data
- .shape_hessians[dof_index1][q][0][d1][d2];
- data.shape_hessians[dof_index2_1][q][1][d1]
- [d2] =
- -1.0 *
- data
- .shape_hessians[dof_index1][q][1][d1][d2];
- data.shape_hessians[dof_index2_1][q][2][d1]
- [d2] =
- data
- .shape_hessians[dof_index1][q][2][d1][d2];
-
- data.shape_hessians[dof_index2_2][q][0][d1]
- [d2] =
- data
- .shape_hessians[dof_index1][q][0][d1][d2];
- data.shape_hessians[dof_index2_2][q][1][d1]
- [d2] =
- -1.0 *
- data
- .shape_hessians[dof_index1][q][1][d1][d2];
- data.shape_hessians[dof_index2_2][q][2][d1]
- [d2] =
- -1.0 *
- data
- .shape_hessians[dof_index1][q][2][d1][d2];
+ const unsigned int shift_ijk(shift_jk + i);
+
+ // Type 1:
+ const unsigned int dof_index1(
+ cell_type1_offset + shift_ijk);
+
+ data.shape_hessians[dof_index1][q][0][0][0] =
+ 8.0 * polyx[i][3] * polyy[j][0] *
+ polyz[k][0];
+ data.shape_hessians[dof_index1][q][1][0][0] =
+ 8.0 * polyx[i][2] * polyy[j][1] *
+ polyz[k][0];
+ data.shape_hessians[dof_index1][q][2][0][0] =
+ 8.0 * polyx[i][2] * polyy[j][0] *
+ polyz[k][1];
+
+ data.shape_hessians[dof_index1][q][0][1][0] =
+ data.shape_hessians[dof_index1][q][1][0][0];
+ data.shape_hessians[dof_index1][q][1][1][0] =
+ 8.0 * polyx[i][1] * polyy[j][2] *
+ polyz[k][0];
+ data.shape_hessians[dof_index1][q][2][1][0] =
+ 8.0 * polyx[i][1] * polyy[j][1] *
+ polyz[k][1];
+
+ data.shape_hessians[dof_index1][q][0][2][0] =
+ data.shape_hessians[dof_index1][q][2][0][0];
+ data.shape_hessians[dof_index1][q][1][2][0] =
+ data.shape_hessians[dof_index1][q][2][1][0];
+ data.shape_hessians[dof_index1][q][2][2][0] =
+ 8.0 * polyx[i][1] * polyy[j][0] *
+ polyz[k][2];
+
+
+ data.shape_hessians[dof_index1][q][0][0][1] =
+ data.shape_hessians[dof_index1][q][1][0][0];
+ data.shape_hessians[dof_index1][q][1][0][1] =
+ data.shape_hessians[dof_index1][q][1][1][0];
+ data.shape_hessians[dof_index1][q][2][0][1] =
+ data.shape_hessians[dof_index1][q][2][1][0];
+
+ data.shape_hessians[dof_index1][q][0][1][1] =
+ data.shape_hessians[dof_index1][q][1][1][0];
+ data.shape_hessians[dof_index1][q][1][1][1] =
+ 8.0 * polyx[i][0] * polyy[j][3] *
+ polyz[k][0];
+ data.shape_hessians[dof_index1][q][2][1][1] =
+ 8.0 * polyx[i][0] * polyy[j][2] *
+ polyz[k][1];
+
+ data.shape_hessians[dof_index1][q][0][2][1] =
+ data.shape_hessians[dof_index1][q][2][1][0];
+ data.shape_hessians[dof_index1][q][1][2][1] =
+ data.shape_hessians[dof_index1][q][2][1][1];
+ data.shape_hessians[dof_index1][q][2][2][1] =
+ 8.0 * polyx[i][0] * polyy[j][1] *
+ polyz[k][2];
+
+
+ data.shape_hessians[dof_index1][q][0][0][2] =
+ data.shape_hessians[dof_index1][q][2][0][0];
+ data.shape_hessians[dof_index1][q][1][0][2] =
+ data.shape_hessians[dof_index1][q][2][1][0];
+ data.shape_hessians[dof_index1][q][2][0][2] =
+ data.shape_hessians[dof_index1][q][2][2][0];
+
+ data.shape_hessians[dof_index1][q][0][1][2] =
+ data.shape_hessians[dof_index1][q][2][1][0];
+ data.shape_hessians[dof_index1][q][1][1][2] =
+ data.shape_hessians[dof_index1][q][2][1][1];
+ data.shape_hessians[dof_index1][q][2][1][2] =
+ data.shape_hessians[dof_index1][q][2][2][1];
+
+ data.shape_hessians[dof_index1][q][0][2][2] =
+ data.shape_hessians[dof_index1][q][2][2][0];
+ data.shape_hessians[dof_index1][q][1][2][2] =
+ data.shape_hessians[dof_index1][q][2][2][1];
+ data.shape_hessians[dof_index1][q][2][2][2] =
+ 8.0 * polyx[i][0] * polyy[j][0] *
+ polyz[k][3];
+
+
+ // Type 2:
+ const unsigned int dof_index2_1(
+ cell_type2_offset1 + shift_ijk);
+ const unsigned int dof_index2_2(
+ cell_type2_offset2 + shift_ijk);
+
+ for (unsigned int d1 = 0; d1 < dim; ++d1)
+ {
+ for (unsigned int d2 = 0; d2 < dim; ++d2)
+ {
+ data.shape_hessians[dof_index2_1][q]
+ [0][d1][d2] =
+ data.shape_hessians[dof_index1][q]
+ [0][d1][d2];
+ data.shape_hessians[dof_index2_1][q]
+ [1][d1][d2] =
+ -1.0 *
+ data.shape_hessians[dof_index1][q]
+ [1][d1][d2];
+ data.shape_hessians[dof_index2_1][q]
+ [2][d1][d2] =
+ data.shape_hessians[dof_index1][q]
+ [2][d1][d2];
+
+ data.shape_hessians[dof_index2_2][q]
+ [0][d1][d2] =
+ data.shape_hessians[dof_index1][q]
+ [0][d1][d2];
+ data.shape_hessians[dof_index2_2][q]
+ [1][d1][d2] =
+ -1.0 *
+ data.shape_hessians[dof_index1][q]
+ [1][d1][d2];
+ data.shape_hessians[dof_index2_2][q]
+ [2][d1][d2] =
+ -1.0 *
+ data.shape_hessians[dof_index1][q]
+ [2][d1][d2];
+ }
+ }
}
- }
- }
- // Type 3: (note we re-use k and j for
- // convenience):
- const unsigned int shift_ij(j +
- shift_j); // here we've subbed
- // j for i, k for j.
- const unsigned int dof_index3_1(cell_type3_offset1 +
- shift_ij);
- const unsigned int dof_index3_2(cell_type3_offset2 +
- shift_ij);
- const unsigned int dof_index3_3(cell_type3_offset3 +
- shift_ij);
- for (unsigned int d1 = 0; d1 < dim; ++d1)
- {
- for (unsigned int d2 = 0; d2 < dim; ++d2)
- {
- for (unsigned int d3 = 0; d3 < dim; ++d3)
+ // Type 3: (note we re-use k and j for
+ // convenience):
+ const unsigned int shift_ij(
+ j + shift_j); // here we've subbed
+ // j for i, k for j.
+ const unsigned int dof_index3_1(
+ cell_type3_offset1 + shift_ij);
+ const unsigned int dof_index3_2(
+ cell_type3_offset2 + shift_ij);
+ const unsigned int dof_index3_3(
+ cell_type3_offset3 + shift_ij);
+ for (unsigned int d1 = 0; d1 < dim; ++d1)
{
- data.shape_hessians[dof_index3_1][q][d1][d2]
- [d3] = 0.0;
- data.shape_hessians[dof_index3_2][q][d1][d2]
- [d3] = 0.0;
- data.shape_hessians[dof_index3_3][q][d1][d2]
- [d3] = 0.0;
+ for (unsigned int d2 = 0; d2 < dim; ++d2)
+ {
+ for (unsigned int d3 = 0; d3 < dim; ++d3)
+ {
+ data.shape_hessians[dof_index3_1][q]
+ [d1][d2][d3] = 0.0;
+ data.shape_hessians[dof_index3_2][q]
+ [d1][d2][d3] = 0.0;
+ data.shape_hessians[dof_index3_3][q]
+ [d1][d2][d3] = 0.0;
+ }
+ }
}
- }
- }
- data.shape_hessians[dof_index3_1][q][0][1][1] =
- 4.0 * polyy[j][2] * polyz[k][0];
- data.shape_hessians[dof_index3_1][q][0][1][2] =
- 4.0 * polyy[j][1] * polyz[k][1];
+ data.shape_hessians[dof_index3_1][q][0][1][1] =
+ 4.0 * polyy[j][2] * polyz[k][0];
+ data.shape_hessians[dof_index3_1][q][0][1][2] =
+ 4.0 * polyy[j][1] * polyz[k][1];
- data.shape_hessians[dof_index3_1][q][0][2][1] =
- data.shape_hessians[dof_index3_1][q][0][1][2];
- data.shape_hessians[dof_index3_1][q][0][2][2] =
- 4.0 * polyy[j][0] * polyz[k][2];
+ data.shape_hessians[dof_index3_1][q][0][2][1] =
+ data.shape_hessians[dof_index3_1][q][0][1][2];
+ data.shape_hessians[dof_index3_1][q][0][2][2] =
+ 4.0 * polyy[j][0] * polyz[k][2];
- data.shape_hessians[dof_index3_2][q][1][0][0] =
- 4.0 * polyx[j][2] * polyz[k][0];
- data.shape_hessians[dof_index3_2][q][1][0][2] =
- 4.0 * polyx[j][1] * polyz[k][1];
+ data.shape_hessians[dof_index3_2][q][1][0][0] =
+ 4.0 * polyx[j][2] * polyz[k][0];
+ data.shape_hessians[dof_index3_2][q][1][0][2] =
+ 4.0 * polyx[j][1] * polyz[k][1];
- data.shape_hessians[dof_index3_2][q][1][2][0] =
- data.shape_hessians[dof_index3_2][q][1][0][2];
- data.shape_hessians[dof_index3_2][q][1][2][2] =
- 4.0 * polyx[j][0] * polyz[k][2];
+ data.shape_hessians[dof_index3_2][q][1][2][0] =
+ data.shape_hessians[dof_index3_2][q][1][0][2];
+ data.shape_hessians[dof_index3_2][q][1][2][2] =
+ 4.0 * polyx[j][0] * polyz[k][2];
- data.shape_hessians[dof_index3_3][q][2][0][0] =
- 4.0 * polyx[j][2] * polyy[k][0];
- data.shape_hessians[dof_index3_3][q][2][0][1] =
- 4.0 * polyx[j][1] * polyy[k][1];
+ data.shape_hessians[dof_index3_3][q][2][0][0] =
+ 4.0 * polyx[j][2] * polyy[k][0];
+ data.shape_hessians[dof_index3_3][q][2][0][1] =
+ 4.0 * polyx[j][1] * polyy[k][1];
- data.shape_hessians[dof_index3_3][q][2][1][0] =
- data.shape_hessians[dof_index3_3][q][2][0][1];
- data.shape_hessians[dof_index3_3][q][2][1][1] =
- 4.0 * polyx[j][0] * polyy[k][2];
+ data.shape_hessians[dof_index3_3][q][2][1][0] =
+ data.shape_hessians[dof_index3_3][q][2][0][1];
+ data.shape_hessians[dof_index3_3][q][2][1][1] =
+ 4.0 * polyx[j][0] * polyy[k][2];
+ }
+ }
}
}
}
+ break;
+ }
+
+ default:
+ {
+ DEAL_II_NOT_IMPLEMENTED();
+ break;
}
}
}
typename dealii::FiniteElement<dim, spacedim>::InternalDataBase>
data_ptr = std::make_unique<InternalData>();
- const unsigned int n_q_points = quadrature.size();
-
+ const unsigned int n_q_points = quadrature.size();
std::vector<Point<dim>> p_list(n_q_points);
p_list = quadrature.get_points();
IntegratedLegendreSZ::generate_complete_basis(degree + 1);
}
+
+
// explicit instantiations
#include "fe_nedelec_sz.inst"
projects / dealii.git / commitdiff