From: Wolfgang Bangerth Date: Sat, 14 Aug 2010 02:44:11 +0000 (+0000) Subject: Simple style changes: indent according to our usual style, break lines X-Git-Tag: v8.0.0~5732 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=771330ceef5ad5675fd1ca6af6083f9fa0ece4c0;p=dealii.git Simple style changes: indent according to our usual style, break lines at 80 characters if possible. git-svn-id: https://svn.dealii.org/trunk@21654 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 0260ff79de..6b1504ac55 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -2922,816 +2922,1142 @@ namespace internal } -namespace internals { - namespace VectorTools { - - // This function computes the projection of the - // boundary function on edges for 3D. +namespace internals +{ + namespace VectorTools + { + + // This function computes the + // projection of the boundary + // function on edges for 3D. template void - compute_edge_projection (const cell_iterator& cell, const unsigned int face, const - unsigned int line, FEValues& fe_values, const Quadrature& quadrature, - const Function& boundary_function, const unsigned int first_vector_component, - std::vector& dof_values) { - fe_values.reinit (cell); - - // Initialize the required objects. - std::vector > jacobians = fe_values.get_jacobians (); - std::vector > tangentials (fe_values.n_quadrature_points); - std::vector > quadrature_points = fe_values.get_quadrature_points (); - std::vector > values (fe_values.n_quadrature_points, - Vector (dim)); - - // Get boundary function values at quadrature points. - boundary_function.vector_value_list (quadrature_points, values); - quadrature_points = quadrature.get_points (); - - const unsigned int superdegree = cell->get_fe ().degree; - const unsigned int degree = superdegree - 1; - Point shifted_reference_point_1; - Point shifted_reference_point_2; - unsigned int edge_coordinate_direction[4]; - - // Get coordinate directions of the edges of the face. - switch (face) { - case 0: case 1: { - edge_coordinate_direction[0] = 2; - edge_coordinate_direction[1] = 2; - edge_coordinate_direction[2] = 1; - edge_coordinate_direction[3] = 1; - break; - } - - case 2: case 3: { - edge_coordinate_direction[0] = 0; - edge_coordinate_direction[1] = 0; - edge_coordinate_direction[2] = 2; - edge_coordinate_direction[3] = 2; - break; - } - - default: { - edge_coordinate_direction[0] = 1; - edge_coordinate_direction[1] = 1; - edge_coordinate_direction[2] = 0; - edge_coordinate_direction[3] = 0; - } - } - - // The interpolation for the lowest order edge shape - // functions is just the mean value of the tangential - // components of the boundary function on the edge. - for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; - ++q_point) { - // Therefore compute the tangential of the edge at the - // quadrature point. - for (unsigned int d = 0; d < dim; ++d) { - shifted_reference_point_1 (d) = quadrature_points[q_point] (d); - shifted_reference_point_2 (d) = quadrature_points[q_point] (d); - } - - shifted_reference_point_1 (edge_coordinate_direction[line]) += 1e-13; - shifted_reference_point_2 (edge_coordinate_direction[line]) -= 1e-13; - tangentials[q_point] = 2e13 * - (fe_values.get_mapping ().transform_unit_to_real_cell (cell, - shifted_reference_point_1) - - fe_values.get_mapping ().transform_unit_to_real_cell (cell, - shifted_reference_point_2)); - tangentials[q_point] /= std::sqrt (tangentials[q_point].square ()); - // Compute the mean value. - dof_values[line * superdegree] += fe_values.JxW (q_point) - * (values[q_point] (0) * tangentials[q_point] (0) + values[q_point] (1) - * tangentials[q_point] (1) + values[q_point] (2) * tangentials[q_point] (2)) - / (jacobians[q_point][0][edge_coordinate_direction[line]] - * jacobians[q_point][0][edge_coordinate_direction[line]] - + jacobians[q_point][1][edge_coordinate_direction[line]] - * jacobians[q_point][1][edge_coordinate_direction[line]] - + jacobians[q_point][2][edge_coordinate_direction[line]] - * jacobians[q_point][2][edge_coordinate_direction[line]]); - } - - // If there are also higher order shape functions we have - // still some work left. - if (degree > 0) { - const FEValuesExtractors::Vector vec (first_vector_component); - FullMatrix assembling_matrix (degree, fe_values.n_quadrature_points); - Tensor<1, dim> shape_value; - Tensor<1, dim> tmp; - Vector assembling_vector (fe_values.n_quadrature_points); - - // We set up a linear system of equations to get the values - // for the remaining degrees of freedom associated with - // the edge. - for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; - ++q_point) { - // The right hand side of the corresponding problem is the - // tangential components of the residual of the boundary - // function and the interpolated part above. - tmp = std::sqrt (fe_values.JxW (q_point) - / (jacobians[q_point][0][edge_coordinate_direction[line]] - * jacobians[q_point][0][edge_coordinate_direction[line]] - + jacobians[q_point][1][edge_coordinate_direction[line]] - * jacobians[q_point][1][edge_coordinate_direction[line]] - + jacobians[q_point][2][edge_coordinate_direction[line]] - * jacobians[q_point][2][edge_coordinate_direction[line]])) - * tangentials[q_point]; - shape_value - = fe_values[vec].value (cell->get_fe ().face_to_cell_index (line * superdegree, face), - q_point); - // In the weak form the right hand side function is multiplicated - // by the higher order shape functions. - assembling_vector (q_point) = (values[q_point] (0) - - dof_values[line * superdegree] * shape_value[0]) * tmp[0] + (values[q_point] (1) - - dof_values[line * superdegree] * shape_value[1]) * tmp[1] + (values[q_point] (2) - - dof_values[line * superdegree] * shape_value[2]) * tmp[2]; - - for (unsigned int i = 0; i < degree; ++i) - assembling_matrix (i, q_point) - = fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + line * superdegree + 1, face), - q_point) * tmp; - } - - FullMatrix cell_matrix (degree, degree); - - // Create the system matrix by multiplying the assembling - // matrix with its transposed. - assembling_matrix.mTmult (cell_matrix, assembling_matrix); - - Vector cell_rhs (degree); - - // Create the system right hand side vector by multiplying - // the assembling matrix with the assembling vector. - assembling_matrix.vmult (cell_rhs, assembling_vector); - - PreconditionJacobi > precondition; - - // Use Jacobi preconditioner with the PCG method to solve the - // problem. - precondition.initialize (cell_matrix); - - SolverControl solver_control (degree, 1e-15, false, false); - SolverCG<> cg (solver_control); - Vector solution (degree); - - cg.solve (cell_matrix, solution, cell_rhs, precondition); - - // Store the computed values. - for (unsigned int i = 0; i < degree; ++i) - dof_values[i + line * superdegree + 1] = solution (i); - } + compute_edge_projection (const cell_iterator& cell, + const unsigned int face, + const unsigned int line, + FEValues& fe_values, + const Quadrature& quadrature, + const Function& boundary_function, + const unsigned int first_vector_component, + std::vector& dof_values) + { + fe_values.reinit (cell); + + // Initialize the required + // objects. + std::vector > jacobians = fe_values.get_jacobians (); + std::vector > tangentials (fe_values.n_quadrature_points); + std::vector > quadrature_points = fe_values.get_quadrature_points (); + std::vector > values (fe_values.n_quadrature_points, + Vector (dim)); + + // Get boundary function values + // at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + quadrature_points = quadrature.get_points (); + + const unsigned int superdegree = cell->get_fe ().degree; + const unsigned int degree = superdegree - 1; + Point shifted_reference_point_1; + Point shifted_reference_point_2; + unsigned int edge_coordinate_direction[4]; + + // Get coordinate directions of + // the edges of the face. + switch (face) + { + case 0: + case 1: + { + edge_coordinate_direction[0] = 2; + edge_coordinate_direction[1] = 2; + edge_coordinate_direction[2] = 1; + edge_coordinate_direction[3] = 1; + break; + } + + case 2: + case 3: + { + edge_coordinate_direction[0] = 0; + edge_coordinate_direction[1] = 0; + edge_coordinate_direction[2] = 2; + edge_coordinate_direction[3] = 2; + break; + } + + default: + { + edge_coordinate_direction[0] = 1; + edge_coordinate_direction[1] = 1; + edge_coordinate_direction[2] = 0; + edge_coordinate_direction[3] = 0; + } + } + + // The interpolation for the + // lowest order edge shape + // functions is just the mean + // value of the tangential + // components of the boundary + // function on the edge. + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + { + // Therefore compute the + // tangential of the edge at + // the quadrature point. + for (unsigned int d = 0; d < dim; ++d) + { + shifted_reference_point_1 (d) = quadrature_points[q_point] (d); + shifted_reference_point_2 (d) = quadrature_points[q_point] (d); + } + + shifted_reference_point_1 (edge_coordinate_direction[line]) += 1e-13; + shifted_reference_point_2 (edge_coordinate_direction[line]) -= 1e-13; + tangentials[q_point] + = (2e13 * + (fe_values.get_mapping () + .transform_unit_to_real_cell (cell, + shifted_reference_point_1) + - + fe_values.get_mapping () + .transform_unit_to_real_cell (cell, + shifted_reference_point_2))); + tangentials[q_point] + /= std::sqrt (tangentials[q_point].square ()); + + // Compute the mean value. + dof_values[line * superdegree] + += (fe_values.JxW (q_point) + * (values[q_point] (0) * tangentials[q_point] (0) + + values[q_point] (1) * tangentials[q_point] (1) + + values[q_point] (2) * tangentials[q_point] (2)) + / (jacobians[q_point][0][edge_coordinate_direction[line]] + * jacobians[q_point][0][edge_coordinate_direction[line]] + + + jacobians[q_point][1][edge_coordinate_direction[line]] + * jacobians[q_point][1][edge_coordinate_direction[line]] + + + jacobians[q_point][2][edge_coordinate_direction[line]] + * jacobians[q_point][2][edge_coordinate_direction[line]])); + } + + // If there are also higher + // order shape functions we + // have still some work left. + if (degree > 0) + { + const FEValuesExtractors::Vector vec (first_vector_component); + FullMatrix assembling_matrix (degree, fe_values.n_quadrature_points); + Tensor<1, dim> shape_value; + Tensor<1, dim> tmp; + Vector assembling_vector (fe_values.n_quadrature_points); + + // We set up a linear system + // of equations to get the + // values for the remaining + // degrees of freedom + // associated with the edge. + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + { + // The right hand side of + // the corresponding + // problem is the + // tangential components of + // the residual of the + // boundary function and + // the interpolated part + // above. + tmp = + std::sqrt (fe_values.JxW (q_point) + / (jacobians[q_point][0][edge_coordinate_direction[line]] + * jacobians[q_point][0][edge_coordinate_direction[line]] + + + jacobians[q_point][1][edge_coordinate_direction[line]] + * jacobians[q_point][1][edge_coordinate_direction[line]] + + + jacobians[q_point][2][edge_coordinate_direction[line]] + * jacobians[q_point][2][edge_coordinate_direction[line]])) + * tangentials[q_point]; + shape_value + = fe_values[vec].value (cell->get_fe () + .face_to_cell_index (line * superdegree, face), + q_point); + // In the weak form the + // right hand side function + // is multiplicated by the + // higher order shape + // functions. + assembling_vector (q_point) + = ((values[q_point] (0) + - + dof_values[line * superdegree] * shape_value[0]) * tmp[0] + + + (values[q_point] (1) + - + dof_values[line * superdegree] * shape_value[1]) * tmp[1] + + + (values[q_point] (2) + - + dof_values[line * superdegree] * shape_value[2]) * tmp[2]); + + for (unsigned int i = 0; i < degree; ++i) + assembling_matrix (i, q_point) + = fe_values[vec].value (cell->get_fe () + .face_to_cell_index (i + line * superdegree + 1, + face), + q_point) * tmp; + } + + FullMatrix cell_matrix (degree, degree); + + // Create the system matrix + // by multiplying the + // assembling matrix with its + // transposed. + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + + Vector cell_rhs (degree); + + // Create the system right + // hand side vector by + // multiplying the assembling + // matrix with the assembling + // vector. + assembling_matrix.vmult (cell_rhs, assembling_vector); + + PreconditionJacobi > precondition; + + // Use Jacobi preconditioner + // with the PCG method to + // solve the problem. + precondition.initialize (cell_matrix); + + SolverControl solver_control (degree, 1e-15, false, false); + SolverCG<> cg (solver_control); + Vector solution (degree); + + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + // Store the computed values. + for (unsigned int i = 0; i < degree; ++i) + dof_values[i + line * superdegree + 1] = solution (i); + } } - - // This function computes the projection of the - // boundary function on the interior of faces in - // 3D. + + + + // This function computes the + // projection of the boundary + // function on the interior of + // faces in 3D. template void - compute_face_projection (const cell_iterator& cell, const unsigned int face, - FEValues& fe_values, const Function& boundary_function, const unsigned - int first_vector_component, - std::vector& dof_values) { - fe_values.reinit (cell); - - // Initialize the required objects. - std::vector > jacobians = fe_values.get_jacobians (); - std::vector > quadrature_points = fe_values.get_quadrature_points (); - std::vector > values (fe_values.n_quadrature_points, - Vector (dim)); - - // Get boundary function values at quadrature points. - boundary_function.vector_value_list (quadrature_points, values); - - const FEValuesExtractors::Vector vec (first_vector_component); - const unsigned int superdegree = cell->get_fe ().degree; - const unsigned int degree = superdegree - 1; - double JxW; - FullMatrix assembling_matrix (degree * superdegree, - dim * fe_values.n_quadrature_points); - Vector assembling_vector (assembling_matrix.n ()); - Vector cell_rhs (assembling_matrix.m ()); - FullMatrix cell_matrix (assembling_matrix.m (), assembling_matrix.m ()); - Vector solution (cell_matrix.m ()); - SolverControl solver_control (cell_matrix.m (), 1e-15, false, false); - SolverCG<> cg (solver_control); - PreconditionJacobi > precondition; - Tensor<1, dim> tmp; - Tensor<1, dim> shape_value; - unsigned int global_face_coordinate_directions[2]; - unsigned int local_face_coordinate_directions[2]; - - // Get coordinate directions of the face. - switch (face) { - case 0: case 1: { - global_face_coordinate_directions[0] = 1; - global_face_coordinate_directions[1] = 2; - local_face_coordinate_directions[0] = 1; - local_face_coordinate_directions[1] = 0; - break; - } - - case 2: case 3: { - global_face_coordinate_directions[0] = 0; - global_face_coordinate_directions[1] = 2; - local_face_coordinate_directions[0] = 0; - local_face_coordinate_directions[1] = 1; - break; - } - - default: { - global_face_coordinate_directions[0] = 0; - global_face_coordinate_directions[1] = 1; - local_face_coordinate_directions[0] = 1; - local_face_coordinate_directions[1] = 0; - } - } - - // The projection is divided into two steps. In the first step we - // project the boundary function on the horizontal shape functions. - // Then the bounary function is projected on the vertical shape - // functions. - // We begin with the horizontal shape functions and set up a linear - // system of equations to get the values for degrees of freedom - // associated with the interior of the face. - for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; - ++q_point) { - // The right hand side of the corresponding problem is the - // residual of the boundary function and the already - // interpolated part on the edges. - for (unsigned int d = 0; d < dim; ++d) - tmp[d] = values[q_point] (d); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j <= degree; ++j) - tmp -= dof_values[(i + 2 * local_face_coordinate_directions[0]) * superdegree + j] - * fe_values[vec].value (cell->get_fe ().face_to_cell_index - ((i + 2 * local_face_coordinate_directions[0]) * superdegree + j, - face), q_point); - - JxW = std::sqrt (fe_values.JxW (q_point) - / ((jacobians[q_point][0][global_face_coordinate_directions[0]] - * jacobians[q_point][0][global_face_coordinate_directions[0]] - + jacobians[q_point][1][global_face_coordinate_directions[0]] - * jacobians[q_point][1][global_face_coordinate_directions[0]] - + jacobians[q_point][2][global_face_coordinate_directions[0]] - * jacobians[q_point][2][global_face_coordinate_directions[0]]) - * (jacobians[q_point][0][global_face_coordinate_directions[1]] - * jacobians[q_point][0][global_face_coordinate_directions[1]] - + jacobians[q_point][1][global_face_coordinate_directions[1]] - * jacobians[q_point][1][global_face_coordinate_directions[1]] - + jacobians[q_point][2][global_face_coordinate_directions[1]] - * jacobians[q_point][2][global_face_coordinate_directions[1]]))); - - // In the weak form the right hand side function is multiplicated - // by the horizontal shape functions defined in the interior of the - // face. - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = JxW * tmp[d]; - - for (unsigned int i = 0; i <= degree; ++i) - for (unsigned int j = 0; j < degree; ++j) { - shape_value = JxW - * fe_values[vec].value (cell->get_fe ().face_to_cell_index - ((i + GeometryInfo::lines_per_face) * degree + j - + GeometryInfo::lines_per_face, face), q_point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i * degree + j, dim * q_point + d) = shape_value[d]; - } - } - - // Create the system matrix by multiplying the assembling - // matrix with its transposed and the right hand side vector - // by mutliplying the assembling matrix with the assembling - // vector. The problem is solved by the PCG method. - assembling_matrix.mTmult (cell_matrix, assembling_matrix); - assembling_matrix.vmult (cell_rhs, assembling_vector); - precondition.initialize (cell_matrix); - cg.solve (cell_matrix, solution, cell_rhs, precondition); - - // Store the computed values. - for (unsigned int i = 0; i <= degree; ++i) - for (unsigned int j = 0; j < degree; ++j) - dof_values[(i + GeometryInfo::lines_per_face) * degree + j - + GeometryInfo::lines_per_face] = solution (i * degree + j); - - // Now we do the same as above with the vertical shape functions - // instead of the horizontal ones. - for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; - ++q_point) { - for (unsigned int d = 0; d < dim; ++d) - tmp[d] = values[q_point] (d); - - for (unsigned int i = 0; i < 2; ++i) - for (unsigned int j = 0; j <= degree; ++j) - tmp - -= dof_values[(i + 2 * local_face_coordinate_directions[1]) * superdegree + j] - * fe_values[vec].value (cell->get_fe ().face_to_cell_index - ((i + 2 * local_face_coordinate_directions[1]) * superdegree + j, - face), q_point); - - JxW = std::sqrt (fe_values.JxW (q_point) - / ((jacobians[q_point][0][global_face_coordinate_directions[0]] - * jacobians[q_point][0][global_face_coordinate_directions[0]] - + jacobians[q_point][1][global_face_coordinate_directions[0]] - * jacobians[q_point][1][global_face_coordinate_directions[0]] - + jacobians[q_point][2][global_face_coordinate_directions[0]] - * jacobians[q_point][2][global_face_coordinate_directions[0]]) - * (jacobians[q_point][0][global_face_coordinate_directions[1]] - * jacobians[q_point][0][global_face_coordinate_directions[1]] - + jacobians[q_point][1][global_face_coordinate_directions[1]] - * jacobians[q_point][1][global_face_coordinate_directions[1]] - + jacobians[q_point][2][global_face_coordinate_directions[1]] - * jacobians[q_point][2][global_face_coordinate_directions[1]]))); - - for (unsigned int d = 0; d < dim; ++d) - assembling_vector (dim * q_point + d) = JxW * tmp[d]; - - for (unsigned int i = 0; i < degree; ++i) - for (unsigned int j = 0; j <= degree; ++j) { - shape_value = JxW - * fe_values[vec].value (cell->get_fe ().face_to_cell_index - ((i + degree + GeometryInfo::lines_per_face) * superdegree + j, - face), q_point); - - for (unsigned int d = 0; d < dim; ++d) - assembling_matrix (i * superdegree + j, dim * q_point + d) = shape_value[d]; - } - } - - assembling_matrix.mTmult (cell_matrix, assembling_matrix); - assembling_matrix.vmult (cell_rhs, assembling_vector); - precondition.initialize (cell_matrix); - cg.solve (cell_matrix, solution, cell_rhs, precondition); - - for (unsigned int i = 0; i < degree; ++i) - for (unsigned int j = 0; j <= degree; ++j) - dof_values[(i + degree + GeometryInfo::lines_per_face) * superdegree + j] - = solution (i * superdegree + j); + compute_face_projection (const cell_iterator& cell, + const unsigned int face, + FEValues& fe_values, + const Function& boundary_function, + const unsigned int first_vector_component, + std::vector& dof_values) + { + fe_values.reinit (cell); + + // Initialize the required objects. + std::vector > jacobians = fe_values.get_jacobians (); + std::vector > quadrature_points = fe_values.get_quadrature_points (); + std::vector > values (fe_values.n_quadrature_points, + Vector (dim)); + + // Get boundary function values + // at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + + const FEValuesExtractors::Vector vec (first_vector_component); + const unsigned int superdegree = cell->get_fe ().degree; + const unsigned int degree = superdegree - 1; + double JxW; + FullMatrix assembling_matrix (degree * superdegree, + dim * fe_values.n_quadrature_points); + Vector assembling_vector (assembling_matrix.n ()); + Vector cell_rhs (assembling_matrix.m ()); + FullMatrix cell_matrix (assembling_matrix.m (), + assembling_matrix.m ()); + Vector solution (cell_matrix.m ()); + SolverControl solver_control (cell_matrix.m (), 1e-15, false, false); + SolverCG<> cg (solver_control); + PreconditionJacobi > precondition; + Tensor<1, dim> tmp; + Tensor<1, dim> shape_value; + unsigned int global_face_coordinate_directions[2]; + unsigned int local_face_coordinate_directions[2]; + + // Get coordinate directions of + // the face. + switch (face) + { + case 0: + case 1: + { + global_face_coordinate_directions[0] = 1; + global_face_coordinate_directions[1] = 2; + local_face_coordinate_directions[0] = 1; + local_face_coordinate_directions[1] = 0; + break; + } + + case 2: + case 3: + { + global_face_coordinate_directions[0] = 0; + global_face_coordinate_directions[1] = 2; + local_face_coordinate_directions[0] = 0; + local_face_coordinate_directions[1] = 1; + break; + } + + default: + { + global_face_coordinate_directions[0] = 0; + global_face_coordinate_directions[1] = 1; + local_face_coordinate_directions[0] = 1; + local_face_coordinate_directions[1] = 0; + } + } + + // The projection is divided + // into two steps. In the first + // step we project the boundary + // function on the horizontal + // shape functions. Then the + // bounary function is + // projected on the vertical + // shape functions. We begin + // with the horizontal shape + // functions and set up a + // linear system of equations + // to get the values for + // degrees of freedom + // associated with the interior + // of the face. + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + { + // The right hand side of the + // corresponding problem is + // the residual of the + // boundary function and the + // already interpolated part + // on the edges. + for (unsigned int d = 0; d < dim; ++d) + tmp[d] = values[q_point] (d); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= degree; ++j) + tmp -= dof_values[(i + 2 * local_face_coordinate_directions[0]) * superdegree + j] + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + 2 * local_face_coordinate_directions[0]) * superdegree + j, + face), q_point); + + JxW = std::sqrt (fe_values.JxW (q_point) + / ((jacobians[q_point][0][global_face_coordinate_directions[0]] + * jacobians[q_point][0][global_face_coordinate_directions[0]] + + + jacobians[q_point][1][global_face_coordinate_directions[0]] + * jacobians[q_point][1][global_face_coordinate_directions[0]] + + + jacobians[q_point][2][global_face_coordinate_directions[0]] + * jacobians[q_point][2][global_face_coordinate_directions[0]]) + * + (jacobians[q_point][0][global_face_coordinate_directions[1]] + * jacobians[q_point][0][global_face_coordinate_directions[1]] + + + jacobians[q_point][1][global_face_coordinate_directions[1]] + * jacobians[q_point][1][global_face_coordinate_directions[1]] + + + jacobians[q_point][2][global_face_coordinate_directions[1]] + * jacobians[q_point][2][global_face_coordinate_directions[1]]))); + + // In the weak form the right + // hand side function is + // multiplicated by the + // horizontal shape functions + // defined in the interior of + // the face. + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = JxW * tmp[d]; + + for (unsigned int i = 0; i <= degree; ++i) + for (unsigned int j = 0; j < degree; ++j) + { + shape_value + = (JxW + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + GeometryInfo::lines_per_face) + * degree + + j + + GeometryInfo::lines_per_face, + face), + q_point)); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i * degree + j, + dim * q_point + d) + = shape_value[d]; + } + } + + // Create the system matrix by + // multiplying the assembling + // matrix with its transposed + // and the right hand side + // vector by mutliplying the + // assembling matrix with the + // assembling vector. The + // problem is solved by the PCG + // method. + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + assembling_matrix.vmult (cell_rhs, assembling_vector); + precondition.initialize (cell_matrix); + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + // Store the computed values. + for (unsigned int i = 0; i <= degree; ++i) + for (unsigned int j = 0; j < degree; ++j) + dof_values[(i + GeometryInfo::lines_per_face) * degree + j + + GeometryInfo::lines_per_face] + = solution (i * degree + j); + + // Now we do the same as above + // with the vertical shape + // functions instead of the + // horizontal ones. + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + { + for (unsigned int d = 0; d < dim; ++d) + tmp[d] = values[q_point] (d); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= degree; ++j) + tmp + -= dof_values[(i + 2 * local_face_coordinate_directions[1]) * superdegree + j] + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + 2 * local_face_coordinate_directions[1]) * superdegree + j, + face), q_point); + + JxW = std::sqrt (fe_values.JxW (q_point) + / ((jacobians[q_point][0][global_face_coordinate_directions[0]] + * jacobians[q_point][0][global_face_coordinate_directions[0]] + + + jacobians[q_point][1][global_face_coordinate_directions[0]] + * jacobians[q_point][1][global_face_coordinate_directions[0]] + + + jacobians[q_point][2][global_face_coordinate_directions[0]] + * jacobians[q_point][2][global_face_coordinate_directions[0]]) + * + (jacobians[q_point][0][global_face_coordinate_directions[1]] + * jacobians[q_point][0][global_face_coordinate_directions[1]] + + + jacobians[q_point][1][global_face_coordinate_directions[1]] + * jacobians[q_point][1][global_face_coordinate_directions[1]] + + + jacobians[q_point][2][global_face_coordinate_directions[1]] + * jacobians[q_point][2][global_face_coordinate_directions[1]]))); + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = JxW * tmp[d]; + + for (unsigned int i = 0; i < degree; ++i) + for (unsigned int j = 0; j <= degree; ++j) + { + shape_value = JxW + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + degree + GeometryInfo::lines_per_face) * superdegree + j, + face), q_point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i * superdegree + j, dim * q_point + d) + = shape_value[d]; + } + } + + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + assembling_matrix.vmult (cell_rhs, assembling_vector); + precondition.initialize (cell_matrix); + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + for (unsigned int i = 0; i < degree; ++i) + for (unsigned int j = 0; j <= degree; ++j) + dof_values[(i + degree + GeometryInfo::lines_per_face) * superdegree + j] + = solution (i * superdegree + j); } - - - // This function computes the projection of the - // boundary function on the faces in 2D. + + + + + // This function computes the + // projection of the boundary + // function on the faces in 2D. template void - compute_face_projection (const cell_iterator& cell, const unsigned int face, - FEValues& fe_values, const Quadrature& quadrature, const Function& - boundary_function, const unsigned int first_vector_component, std::vector& - dof_values) { - fe_values.reinit (cell); - - // Initialize the required objects. - std::vector > tangentials (fe_values.n_quadrature_points); - std::vector > jacobians = fe_values.get_jacobians (); - std::vector > quadrature_points = fe_values.get_quadrature_points (); - std::vector > values (fe_values.n_quadrature_points, - Vector (dim)); - - // Get boundary function values at quadrature points. - boundary_function.vector_value_list (quadrature_points, values); - quadrature_points = quadrature.get_points (); - - const unsigned int degree = cell->get_fe ().degree - 1; - Point shifted_reference_point_1; - Point shifted_reference_point_2; - Tensor<1, dim> tmp; - Tensor<1, dim> shape_value; - unsigned int face_coordinate_direction; - - // Get coordinate directions of the face. - switch (face) { - case 0: case 1: { - face_coordinate_direction = 1; - break; - } - - default: - face_coordinate_direction = 0; - } - - // The interpolation for the lowest order face shape - // functions is just the mean value of the tangential - // components of the boundary function on the edge. - for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; - ++q_point) { - // Therefore compute the tangential of the face at the - // quadrature point. - for (unsigned int d = 0; d < dim; ++d) { - shifted_reference_point_1 (d) = quadrature_points[q_point] (d); - shifted_reference_point_2 (d) = quadrature_points[q_point] (d); - } - - shifted_reference_point_1 (face_coordinate_direction) += 1e-13; - shifted_reference_point_2 (face_coordinate_direction) -= 1e-13; - tangentials[q_point] = 2e13 - * (fe_values.get_mapping ().transform_unit_to_real_cell (cell, shifted_reference_point_1) - - fe_values.get_mapping ().transform_unit_to_real_cell (cell, shifted_reference_point_2)); - tangentials[q_point] /= std::sqrt (tangentials[q_point].square ()); - // Compute the mean value. - dof_values[0] += fe_values.JxW (q_point) * (values[q_point] (0) - * tangentials[q_point] (0) + values[q_point] (1) * tangentials[q_point] (1)) - / (jacobians[q_point][0][face_coordinate_direction] - * jacobians[q_point][0][face_coordinate_direction] - + jacobians[q_point][1][face_coordinate_direction] - * jacobians[q_point][1][face_coordinate_direction]); - } - - // If there are also higher order shape functions we have - // still some work left. - if (degree > 0) { - const FEValuesExtractors::Vector vec (first_vector_component); - FullMatrix assembling_matrix (degree, fe_values.n_quadrature_points); - Vector assembling_vector (fe_values.n_quadrature_points); - - // We set up a linear system of equations to get the values - // for the remaining degrees of freedom associated with - // the face. - for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; - ++q_point) { - // The right hand side of the corresponding problem is the - // tangential components of the residual of the boundary - // function and the interpolated part above. - tmp = std::sqrt (fe_values.JxW (q_point) - / std::sqrt (jacobians[q_point][0][face_coordinate_direction] - * jacobians[q_point][0][face_coordinate_direction] - + jacobians[q_point][1][face_coordinate_direction] - * jacobians[q_point][1][face_coordinate_direction])) * tangentials[q_point]; - shape_value - = fe_values[vec].value (cell->get_fe ().face_to_cell_index (0, face), q_point); - assembling_vector (q_point) = (values[q_point] (0) - - dof_values[0] * shape_value[0]) * tmp[0] + (values[q_point] (1) - - dof_values[1] * shape_value[1]) * tmp[1]; - - // In the weak form the right hand side function is multiplicated - // by the higher order shape functions. - for (unsigned int i = 0; i < degree; ++i) - assembling_matrix (i, q_point) - = fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + 1, face), - q_point) * tmp; - } - - FullMatrix cell_matrix (degree, degree); - - // Create the system matrix by multiplying the assembling - // matrix with its transposed. - assembling_matrix.mTmult (cell_matrix, assembling_matrix); - - Vector cell_rhs (degree); - - // Create the system right hand side vector by multiplying - // the assembling matrix with the assembling vector. - assembling_matrix.vmult (cell_rhs, assembling_vector); - - PreconditionJacobi > precondition; - - // Use Jacobi preconditioner with the PCG method to solve the - // problem. - precondition.initialize (cell_matrix); - - SolverControl solver_control (degree, 1e-15, false, false); - SolverCG<> cg (solver_control); - Vector solution (degree); - - cg.solve (cell_matrix, solution, cell_rhs, precondition); - - // Store the computed values. - for (unsigned int i = 0; i < degree; ++i) - dof_values[i + 1] = solution (i); - } + compute_face_projection (const cell_iterator& cell, + const unsigned int face, + FEValues& fe_values, + const Quadrature& quadrature, + const Function& boundary_function, + const unsigned int first_vector_component, + std::vector& dof_values) + { + fe_values.reinit (cell); + + // Initialize the required objects. + std::vector > tangentials (fe_values.n_quadrature_points); + std::vector > jacobians = fe_values.get_jacobians (); + std::vector > quadrature_points = fe_values.get_quadrature_points (); + std::vector > values (fe_values.n_quadrature_points, + Vector (dim)); + + // Get boundary function values + // at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + quadrature_points = quadrature.get_points (); + + const unsigned int degree = cell->get_fe ().degree - 1; + Point shifted_reference_point_1; + Point shifted_reference_point_2; + Tensor<1, dim> tmp; + Tensor<1, dim> shape_value; + unsigned int face_coordinate_direction; + + // Get coordinate directions of the face. + switch (face) + { + case 0: + case 1: + { + face_coordinate_direction = 1; + break; + } + + default: + face_coordinate_direction = 0; + } + + // The interpolation for the + // lowest order face shape + // functions is just the mean + // value of the tangential + // components of the boundary + // function on the edge. + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + { + // Therefore compute the + // tangential of the face at + // the quadrature point. + for (unsigned int d = 0; d < dim; ++d) + { + shifted_reference_point_1 (d) = quadrature_points[q_point] (d); + shifted_reference_point_2 (d) = quadrature_points[q_point] (d); + } + + shifted_reference_point_1 (face_coordinate_direction) += 1e-13; + shifted_reference_point_2 (face_coordinate_direction) -= 1e-13; + tangentials[q_point] = 2e13 + * (fe_values.get_mapping () + .transform_unit_to_real_cell (cell, + shifted_reference_point_1) + - + fe_values.get_mapping () + .transform_unit_to_real_cell (cell, + shifted_reference_point_2)); + tangentials[q_point] /= std::sqrt (tangentials[q_point].square ()); + // Compute the mean value. + dof_values[0] += fe_values.JxW (q_point) + * (values[q_point] (0) + * tangentials[q_point] (0) + + + values[q_point] (1) * tangentials[q_point] (1)) + / (jacobians[q_point][0][face_coordinate_direction] + * jacobians[q_point][0][face_coordinate_direction] + + jacobians[q_point][1][face_coordinate_direction] + * jacobians[q_point][1][face_coordinate_direction]); + } + + // If there are also higher + // order shape functions we + // have still some work left. + if (degree > 0) + { + const FEValuesExtractors::Vector vec (first_vector_component); + FullMatrix assembling_matrix (degree, fe_values.n_quadrature_points); + Vector assembling_vector (fe_values.n_quadrature_points); + + // We set up a linear system + // of equations to get the + // values for the remaining + // degrees of freedom + // associated with the face. + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + { + // The right hand side of + // the corresponding + // problem is the + // tangential components of + // the residual of the + // boundary function and + // the interpolated part + // above. + tmp = std::sqrt (fe_values.JxW (q_point) + / std::sqrt (jacobians[q_point][0][face_coordinate_direction] + * jacobians[q_point][0][face_coordinate_direction] + + jacobians[q_point][1][face_coordinate_direction] + * jacobians[q_point][1][face_coordinate_direction])) * tangentials[q_point]; + shape_value + = fe_values[vec].value (cell->get_fe ().face_to_cell_index (0, face), q_point); + assembling_vector (q_point) = (values[q_point] (0) + - dof_values[0] * shape_value[0]) * tmp[0] + (values[q_point] (1) + - dof_values[1] * shape_value[1]) * tmp[1]; + + // In the weak form the + // right hand side function + // is multiplicated by the + // higher order shape + // functions. + for (unsigned int i = 0; i < degree; ++i) + assembling_matrix (i, q_point) + = fe_values[vec].value (cell->get_fe () + .face_to_cell_index (i + 1, face), + q_point) * tmp; + } + + FullMatrix cell_matrix (degree, degree); + + // Create the system matrix + // by multiplying the + // assembling matrix with its + // transposed. + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + + Vector cell_rhs (degree); + + // Create the system right + // hand side vector by + // multiplying the assembling + // matrix with the assembling + // vector. + assembling_matrix.vmult (cell_rhs, assembling_vector); + + PreconditionJacobi > precondition; + + // Use Jacobi preconditioner + // with the PCG method to + // solve the problem. + precondition.initialize (cell_matrix); + + SolverControl solver_control (degree, 1e-15, false, false); + SolverCG<> cg (solver_control); + Vector solution (degree); + + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + // Store the computed values. + for (unsigned int i = 0; i < degree; ++i) + dof_values[i + 1] = solution (i); + } } } } - - - // Projection-based interpolation is performed in two (in 2D) - // respectively three (in 3D) steps. First the tangential - // component of the function is interpolated on each edge. - // This gives the values for the degrees of freedom corresponding - // to the lowest order edge shape functions. Then the interpolated - // part of the function is subtracted and we project the tangential - // component of the residual onto the space of the remaining - // (higher order) edge shape functions. This is done by building - // a linear system of equations of dimension degree. The - // solution gives us the values for the degrees of freedom - // corresponding to the remaining edge shape functions. Now we are - // done for 2D, but in 3D we possibly have also degrees of freedom, - // which are located in the interior of the faces. Therefore we - // compute the residual of the function describing the boundary - // values and the interpolated part, which we have computed in the - // last two steps. On the faces there are two kinds of shape - // functions, the horizontal and the vertical ones. Thus we have - // two solve two linear systems of equations of size - // degree * (degree + 1) to obtain the values for the - // corresponding degrees of freedom. + + template -void VectorTools::project_boundary_values_curl_conforming (const DoFHandler& dof_handler, - const unsigned int first_vector_component, - const Function& boundary_function, - const unsigned char boundary_component, - ConstraintMatrix& constraints, - const Mapping& mapping) +void +VectorTools:: +project_boundary_values_curl_conforming (const DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const Mapping& mapping) { - std::vector dof_values; - std::vector face_dof_indices; - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); - unsigned int dofs_per_face; - unsigned int superdegree; - - switch (dim) { - case 2: { - for (; cell != dof_handler.end (); ++cell) - if (cell->at_boundary ()) - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) - if (cell->face (face)->boundary_indicator () == boundary_component) { - // this is only implemented, if the FE is a Nedelec element - typedef FiniteElement FEL; - - AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), - typename FEL::ExcInterpolationNotImplemented ()); - - dofs_per_face = cell->get_fe ().dofs_per_face; - dof_values.resize (dofs_per_face); - - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - dof_values[dof] = 0.0; - - superdegree = cell->get_fe ().degree; - - QGauss reference_face_quadrature (2 * superdegree); - Quadrature face_quadrature - = QProjector::project_to_face (reference_face_quadrature, face); - FEValues fe_face_values (mapping, cell->get_fe (), face_quadrature, - update_jacobians | update_JxW_values | update_quadrature_points | update_values); - - // Compute the projection of the boundary function on the edge. - internals::VectorTools::compute_face_projection (cell, face, fe_face_values, - face_quadrature, boundary_function, first_vector_component, dof_values); - face_dof_indices.resize (dofs_per_face); - cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); - - // Add the computed constraints to the constraint matrix. - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) { - constraints.add_line (face_dof_indices[dof]); - - if (std::abs (dof_values[dof]) > 1e-14) - constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); - } - } - - break; + // Projection-based interpolation + // is performed in two (in 2D) + // respectively three (in 3D) + // steps. First the tangential + // component of the function is + // interpolated on each edge. This + // gives the values for the degrees + // of freedom corresponding to the + // lowest order edge shape + // functions. Then the interpolated + // part of the function is + // subtracted and we project the + // tangential component of the + // residual onto the space of the + // remaining (higher order) edge + // shape functions. This is done by + // building a linear system of + // equations of dimension + // degree. The solution + // gives us the values for the + // degrees of freedom corresponding + // to the remaining edge shape + // functions. Now we are done for + // 2D, but in 3D we possibly have + // also degrees of freedom, which + // are located in the interior of + // the faces. Therefore we compute + // the residual of the function + // describing the boundary values + // and the interpolated part, which + // we have computed in the last two + // steps. On the faces there are + // two kinds of shape functions, + // the horizontal and the vertical + // ones. Thus we have two solve two + // linear systems of equations of + // size degree * (degree + + // 1) to obtain the values for + // the corresponding degrees of + // freedom. + std::vector dof_values; + std::vector face_dof_indices; + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + unsigned int dofs_per_face; + unsigned int superdegree; + + switch (dim) + { + case 2: + { + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) + { + // this is only + // implemented, if the + // FE is a Nedelec + // element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + superdegree = cell->get_fe ().degree; + + QGauss reference_face_quadrature (2 * superdegree); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, face); + FEValues fe_face_values (mapping, cell->get_fe (), + face_quadrature, + update_jacobians | + update_JxW_values | + update_quadrature_points | + update_values); + + // Compute the + // projection of the + // boundary function on + // the edge. + internals::VectorTools + ::compute_face_projection (cell, face, fe_face_values, + face_quadrature, boundary_function, + first_vector_component, dof_values); + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, + cell->active_fe_index ()); + + // Add the computed + // constraints to the + // constraint matrix. + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + { + constraints.add_line (face_dof_indices[dof]); + + if (std::abs (dof_values[dof]) > 1e-14) + constraints.set_inhomogeneity (face_dof_indices[dof], + dof_values[dof]); + } + } + + break; } - - case 3: { - const unsigned int n_dofs = dof_handler.n_dofs (); - std::vector computed_constraints (n_dofs); - std::vector projected_dofs (n_dofs); - unsigned int degree; - unsigned int superdegree; - - for (unsigned int dof = 0; dof < n_dofs; ++dof) - projected_dofs[dof] = -1; - - for (; cell != dof_handler.end (); ++cell) - if (cell->at_boundary ()) - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) - if (cell->face (face)->boundary_indicator () == boundary_component) { - // this is only implemented, if the FE is a Nedelec element - typedef FiniteElement FEL; - - AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), - typename FEL::ExcInterpolationNotImplemented ()); - - superdegree = cell->get_fe ().degree; - degree = superdegree - 1; - - QGauss reference_edge_quadrature (2 * superdegree); - - dofs_per_face = cell->get_fe ().dofs_per_face; - dof_values.resize (dofs_per_face); - - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - dof_values[dof] = 0.0; - - face_dof_indices.resize (dofs_per_face); - cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); - - // First we compute the projection on the edges. - for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) { - // If we have reached this edge through another cell before, we do - // not do here anything unless we have a good reason, i.e. a higher - // polynomial degree. - if (projected_dofs[face_dof_indices[line * superdegree]] < (int) degree) { - Quadrature edge_quadrature - = QProjector::project_to_face (QProjector::project_to_face - (reference_edge_quadrature, line), face); - FEValues fe_edge_values (mapping, cell->get_fe (), edge_quadrature, - update_JxW_values | update_jacobians | update_quadrature_points | update_values); - // Compute the projection of the boundary function on the edge. - internals::VectorTools::compute_edge_projection (cell, face, line, - fe_edge_values, edge_quadrature, boundary_function, first_vector_component, - dof_values); - // Mark the projected degrees of freedom. - for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; - ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } - - // If we have computed the values in a previous step of the loop, - // we just copy the values in the local vector. - else - for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; - ++dof) - dof_values[dof] = computed_constraints[face_dof_indices[dof]]; - } - - // If there are higher order shape functions, there is still some - // work left. - if (degree > 0) { - QGauss reference_face_quadrature (2 * superdegree); - Quadrature face_quadrature - = QProjector::project_to_face (reference_face_quadrature, face); - FEValues fe_face_values (mapping, cell->get_fe (), face_quadrature, - update_JxW_values | update_jacobians | update_quadrature_points | update_values); - - // Compute the projection of the boundary function on the interior - // of the face. - internals::VectorTools::compute_face_projection (cell, face, fe_face_values, - boundary_function, first_vector_component, dof_values); - - // Mark the projected degrees of freedom. - for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; - dof < dofs_per_face; ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } - - // Store the computed values in the global vector. - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - if (std::abs (dof_values[dof]) > 1e-14) - computed_constraints[face_dof_indices[dof]] = dof_values[dof]; - } - - // Add the computed constraints to the constraint matrix. - for (unsigned int dof = 0; dof < n_dofs; ++dof) - if (projected_dofs[dof] != -1) { - constraints.add_line (dof); - constraints.set_inhomogeneity (dof, computed_constraints[dof]); - } + + case 3: + { + const unsigned int n_dofs = dof_handler.n_dofs (); + std::vector computed_constraints (n_dofs); + std::vector projected_dofs (n_dofs); + unsigned int degree; + unsigned int superdegree; + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + projected_dofs[dof] = -1; + + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) + { + // this is only + // implemented, if the + // FE is a Nedelec + // element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + superdegree = cell->get_fe ().degree; + degree = superdegree - 1; + + QGauss reference_edge_quadrature (2 * superdegree); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, + cell->active_fe_index ()); + + // First we compute the + // projection on the + // edges. + for (unsigned int line = 0; + line < GeometryInfo<3>::lines_per_face; ++line) + { + // If we have reached + // this edge through + // another cell + // before, we do not + // do here anything + // unless we have a + // good reason, + // i.e. a higher + // polynomial degree. + if (projected_dofs[face_dof_indices[line * superdegree]] + < + (int) degree) + { + Quadrature edge_quadrature + = QProjector::project_to_face (QProjector::project_to_face + (reference_edge_quadrature, line), face); + FEValues fe_edge_values (mapping, cell->get_fe (), + edge_quadrature, + update_JxW_values | + update_jacobians | + update_quadrature_points | + update_values); + // Compute the + // projection of + // the boundary + // function on the + // edge. + internals::VectorTools + ::compute_edge_projection (cell, face, line, + fe_edge_values, + edge_quadrature, + boundary_function, + first_vector_component, + dof_values); + // Mark the + // projected + // degrees of + // freedom. + for (unsigned int dof = line * superdegree; + dof < (line + 1) * superdegree; ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + // If we have + // computed the + // values in a + // previous step of + // the loop, we just + // copy the values in + // the local vector. + else + for (unsigned int dof = line * superdegree; + dof < (line + 1) * superdegree; + ++dof) + dof_values[dof] = computed_constraints[face_dof_indices[dof]]; + } + + // If there are higher + // order shape + // functions, there is + // still some work + // left. + if (degree > 0) + { + QGauss reference_face_quadrature (2 * superdegree); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, + face); + FEValues fe_face_values (mapping, cell->get_fe (), + face_quadrature, + update_JxW_values | + update_jacobians | + update_quadrature_points | + update_values); + + // Compute the + // projection of the + // boundary function + // on the interior of + // the face. + internals::VectorTools + ::compute_face_projection (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values); + + // Mark the projected + // degrees of + // freedom. + for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; + dof < dofs_per_face; ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + // Store the computed + // values in the global + // vector. + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + if (std::abs (dof_values[dof]) > 1e-14) + computed_constraints[face_dof_indices[dof]] = dof_values[dof]; + } + + // Add the computed constraints + // to the constraint matrix. + for (unsigned int dof = 0; dof < n_dofs; ++dof) + if (projected_dofs[dof] != -1) + { + constraints.add_line (dof); + constraints.set_inhomogeneity (dof, computed_constraints[dof]); + } } - } + } } + template -void VectorTools::project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, - const unsigned int first_vector_component, - const Function& boundary_function, - const unsigned char boundary_component, - ConstraintMatrix& constraints, - const hp::MappingCollection& mapping_collection) +void +VectorTools:: +project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const hp::MappingCollection& mapping_collection) { - std::vector dof_values; - std::vector face_dof_indices; - typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); - unsigned int dofs_per_face; - - switch (dim) { - case 2: { - for (; cell != dof_handler.end (); ++cell) - if (cell->at_boundary ()) - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) - if (cell->face (face)->boundary_indicator () == boundary_component) { - // this is only implemented, if the FE is a Nédélec element - typedef FiniteElement FEL; - - AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), - typename FEL::ExcInterpolationNotImplemented ()); - - dofs_per_face = cell->get_fe ().dofs_per_face; - dof_values.resize (dofs_per_face); - - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - dof_values[dof] = 0.0; - - QGauss reference_face_quadrature (2 * (cell->get_fe ().degree)); - Quadrature face_quadrature - = QProjector::project_to_face (reference_face_quadrature, face); - FEValues fe_face_values (mapping_collection[cell->active_fe_index ()], - cell->get_fe (), face_quadrature, update_jacobians | update_JxW_values - | update_quadrature_points | update_values); - - internals::VectorTools::compute_face_projection (cell, face, fe_face_values, - face_quadrature, boundary_function, first_vector_component, dof_values); - face_dof_indices.resize (dofs_per_face); - cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); - - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) { - constraints.add_line (face_dof_indices[dof]); - - if (std::abs (dof_values[dof]) > 1e-14) - constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); - } - } - - break; + std::vector dof_values; + std::vector face_dof_indices; + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + unsigned int dofs_per_face; + + switch (dim) + { + case 2: + { + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) + { + // this is only + // implemented, if the + // FE is a Nedelec + // element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + QGauss reference_face_quadrature (2 * (cell->get_fe ().degree)); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, face); + FEValues fe_face_values (mapping_collection[cell->active_fe_index ()], + cell->get_fe (), face_quadrature, + update_jacobians | + update_JxW_values | + update_quadrature_points | + update_values); + + internals::VectorTools + ::compute_face_projection (cell, face, fe_face_values, + face_quadrature, + boundary_function, + first_vector_component, + dof_values); + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, + cell->active_fe_index ()); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + { + constraints.add_line (face_dof_indices[dof]); + + if (std::abs (dof_values[dof]) > 1e-14) + constraints.set_inhomogeneity (face_dof_indices[dof], + dof_values[dof]); + } + } + + break; } - - case 3: { - const unsigned int n_dofs = dof_handler.n_dofs (); - std::vector computed_constraints (n_dofs); - std::vector projected_dofs (n_dofs); - unsigned int degree; - unsigned int superdegree; - - for (unsigned int dof = 0; dof < n_dofs; ++dof) - projected_dofs[dof] = -1; - - for (; cell != dof_handler.end (); ++cell) - if (cell->at_boundary ()) - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) - if (cell->face (face)->boundary_indicator () == boundary_component) { - // this is only implemented, if the FE is a Nédélec element - typedef FiniteElement FEL; - - AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), - typename FEL::ExcInterpolationNotImplemented ()); - - superdegree = cell->get_fe ().degree; - degree = superdegree - 1; - - QGauss reference_edge_quadrature (2 * superdegree); - - dofs_per_face = cell->get_fe ().dofs_per_face; - dof_values.resize (dofs_per_face); - - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - dof_values[dof] = 0.0; - - face_dof_indices.resize (dofs_per_face); - cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); - - for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) { - if (projected_dofs[face_dof_indices[line * superdegree]] < (int) degree) { - Quadrature edge_quadrature = QProjector::project_to_face - (QProjector::project_to_face (reference_edge_quadrature, line), - face); - FEValues fe_edge_values (mapping_collection[cell->active_fe_index ()], - cell->get_fe (), edge_quadrature, update_JxW_values | update_jacobians - | update_quadrature_points | update_values); - - internals::VectorTools::compute_edge_projection (cell, face, line, - fe_edge_values, edge_quadrature, boundary_function, first_vector_component, - dof_values); - - for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; - ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } - - else - for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; - ++dof) - dof_values[dof] = computed_constraints[face_dof_indices[dof]]; - } - - if (degree > 0) { - QGauss reference_face_quadrature (2 * superdegree); - Quadrature face_quadrature - = QProjector::project_to_face (reference_face_quadrature, face); - FEValues fe_face_values (mapping_collection[cell->active_fe_index ()], - cell->get_fe (), face_quadrature, update_JxW_values | update_jacobians - | update_quadrature_points | update_values); - - internals::VectorTools::compute_face_projection (cell, face, fe_face_values, - boundary_function, first_vector_component, dof_values); - - for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; - dof < dofs_per_face; ++dof) - projected_dofs[face_dof_indices[dof]] = degree; - } - - for (unsigned int dof = 0; dof < dofs_per_face; ++dof) - if (std::abs (dof_values[dof]) > 1e-14) - computed_constraints[face_dof_indices[dof]] = dof_values[dof]; - } - - for (unsigned int dof = 0; dof < n_dofs; ++dof) - if (projected_dofs[dof] != -1) { - constraints.add_line (dof); - constraints.set_inhomogeneity (dof, computed_constraints[dof]); - } + + case 3: + { + const unsigned int n_dofs = dof_handler.n_dofs (); + std::vector computed_constraints (n_dofs); + std::vector projected_dofs (n_dofs); + unsigned int degree; + unsigned int superdegree; + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + projected_dofs[dof] = -1; + + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) + { + // this is only + // implemented, if the + // FE is a Nedelec + // element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + superdegree = cell->get_fe ().degree; + degree = superdegree - 1; + + QGauss reference_edge_quadrature (2 * superdegree); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, + cell->active_fe_index ()); + + for (unsigned int line = 0; + line < GeometryInfo::lines_per_face; ++line) + { + if (projected_dofs[face_dof_indices[line * superdegree]] + < + (int) degree) + { + Quadrature edge_quadrature = QProjector::project_to_face + (QProjector::project_to_face (reference_edge_quadrature, line), + face); + FEValues fe_edge_values (mapping_collection[cell->active_fe_index ()], + cell->get_fe (), + edge_quadrature, + update_JxW_values | + update_jacobians | + update_quadrature_points | + update_values); + + internals::VectorTools + ::compute_edge_projection (cell, face, line, + fe_edge_values, + edge_quadrature, + boundary_function, + first_vector_component, + dof_values); + + for (unsigned int dof = line * superdegree; + dof < (line + 1) * superdegree; + ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + else + for (unsigned int dof = line * superdegree; + dof < (line + 1) * superdegree; + ++dof) + dof_values[dof] = computed_constraints[face_dof_indices[dof]]; + } + + if (degree > 0) + { + QGauss reference_face_quadrature (2 * superdegree); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, + face); + FEValues fe_face_values (mapping_collection[cell->active_fe_index ()], + cell->get_fe (), + face_quadrature, + update_JxW_values | + update_jacobians | + update_quadrature_points | + update_values); + + internals::VectorTools + ::compute_face_projection (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values); + + for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; + dof < dofs_per_face; ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + if (std::abs (dof_values[dof]) > 1e-14) + computed_constraints[face_dof_indices[dof]] = dof_values[dof]; + } + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + if (projected_dofs[dof] != -1) + { + constraints.add_line (dof); + constraints.set_inhomogeneity (dof, + computed_constraints[dof]); + } } - } + } } + template class DH, int spacedim> void -VectorTools::compute_no_normal_flux_constraints (const DH &dof_handler, - const unsigned int first_vector_component, - const std::set &boundary_ids, - ConstraintMatrix &constraints, - const Mapping &mapping) +VectorTools:: +compute_no_normal_flux_constraints (const DH &dof_handler, + const unsigned int first_vector_component, + const std::set &boundary_ids, + ConstraintMatrix &constraints, + const Mapping &mapping) { Assert (dim > 1, ExcMessage ("This function is not useful in 1d because it amounts "