From 2c4da4bc1609495168add2cd5fec058010dfdc0e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 15 Aug 2023 04:21:41 -0600 Subject: [PATCH] Modernize a test. --- tests/fe/rt_normal_02.cc | 154 ++++++++------------------------------- 1 file changed, 31 insertions(+), 123 deletions(-) diff --git a/tests/fe/rt_normal_02.cc b/tests/fe/rt_normal_02.cc index 00cd9a72ff..a3337afd15 100644 --- a/tests/fe/rt_normal_02.cc +++ b/tests/fe/rt_normal_02.cc @@ -58,7 +58,8 @@ */ void -EvaluateNormal2(DoFHandler<2> *dof_handler, Vector &solution) +evaluate_normal_component(const DoFHandler<2> &dof_handler, + Vector & solution) { // This quadrature rule determines the points, where the // continuity will be tested. @@ -67,41 +68,38 @@ EvaluateNormal2(DoFHandler<2> *dof_handler, Vector &solution) QProjector<2>::project_to_all_faces(ReferenceCells::Quadrilateral, quad); FEFaceValues<2> fe_v_face( - dof_handler->get_fe(), + dof_handler.get_fe(), quad, UpdateFlags(update_values | update_quadrature_points | update_gradients | update_normal_vectors | update_JxW_values)); - FEValues<2> fe_v(dof_handler->get_fe(), + FEValues<2> fe_v(dof_handler.get_fe(), qproject, UpdateFlags(update_values | update_quadrature_points | update_gradients | update_JxW_values)); - FEValues<2> fe_v_n(dof_handler->get_fe(), + FEValues<2> fe_v_n(dof_handler.get_fe(), qproject, UpdateFlags(update_values | update_quadrature_points | update_gradients | update_JxW_values)); const unsigned int n_q_face = quad.size(); const unsigned int n_q_proj = qproject.size(); - const unsigned int n_components = dof_handler->get_fe().n_components(); - const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell; + const unsigned int n_components = dof_handler.get_fe().n_components(); + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; deallog << "Quad Points Face " << n_q_face << ", Quad Points Proj. " << n_q_proj << std::endl; // Cell iterators - DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(), - endc = dof_handler->end(); - std::vector local_dof_indices(dofs_per_cell); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { cell->get_dof_indices(local_dof_indices); fe_v.reinit(cell); - for (const unsigned int f : GeometryInfo<2>::face_indices()) + for (const unsigned int f : cell->face_indices()) { if (!cell->face(f)->at_boundary()) { @@ -121,7 +119,7 @@ EvaluateNormal2(DoFHandler<2> *dof_handler, Vector &solution) const QProjector<2>::DataSetDescriptor offset_n = (QProjector<2>::DataSetDescriptor::face( - ReferenceCells::Quadrilateral, + cell_n->reference_cell(), neighbor, cell_n->face_orientation(neighbor), cell_n->face_flip(neighbor), @@ -140,107 +138,24 @@ EvaluateNormal2(DoFHandler<2> *dof_handler, Vector &solution) for (const auto q_point : fe_v_face.quadrature_point_indices()) { - Tensor<1, 2> vn = fe_v_face.normal_vector(q_point); - double nx = vn[0]; - double ny = vn[1]; - - double u = this_value[q_point + offset](0); - double v = this_value[q_point + offset](1); - - double u_n = this_value_n[q_point + offset_n](0); - double v_n = this_value_n[q_point + offset_n](1); - double un1 = u * nx + v * ny, un2 = u_n * nx + v_n * ny; - - deallog << "QP " << q_point - << ", Error: " << (u - u_n) * nx + (v - v_n) * ny - << ", u " << un1 << ", un " << un2 << ", Rat " - << un2 / un1 << std::endl; - - Assert(std::fabs((u - u_n) * nx + (v - v_n) * ny) < 1e-12, - ExcInternalError()); - } - } - } - } -} - - + const Tensor<1, 2> n = fe_v_face.normal_vector(q_point); + const double nx = n[0]; + const double ny = n[1]; -/* - * Check if the normal component is continuous over element edges. - */ + const Tensor<1, 2> u({this_value[q_point + offset](0), + this_value[q_point + offset](1)}); -void -EvaluateNormal(DoFHandler<2> *dof_handler, Vector &solution) -{ - // This quadrature rule determines the points, where the - // continuity will be tested. - QGauss<1> quad(6); - FEFaceValues<2> fe_v_face( - dof_handler->get_fe(), - quad, - UpdateFlags(update_values | update_quadrature_points | update_gradients | - update_normal_vectors | update_JxW_values)); + const Tensor<1, 2> u_n({this_value_n[q_point + offset_n](0), + this_value_n[q_point + offset_n](1)}); - FEFaceValues<2> fe_v_face_n( - dof_handler->get_fe(), - quad, - UpdateFlags(update_values | update_quadrature_points | update_gradients | - update_normal_vectors | update_JxW_values)); + const double un1 = u * n; + const double un2 = u_n * n; - const unsigned int n_q_face = quad.size(); - const unsigned int n_components = dof_handler->get_fe().n_components(); - const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell; - - // Cell iterators - DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(), - endc = dof_handler->end(); - - std::vector local_dof_indices(dofs_per_cell); - - for (; cell != endc; ++cell) - { - cell->get_dof_indices(local_dof_indices); - - for (const unsigned int f : GeometryInfo<2>::face_indices()) - { - if (!cell->face(f)->at_boundary()) - { - fe_v_face.reinit(cell, f); - - const unsigned int neighbor = cell->neighbor_of_neighbor(f); - fe_v_face_n.reinit(cell->neighbor(f), neighbor); - - // Get values from solution vector (For Trap.Rule) - std::vector> this_value( - n_q_face, Vector(n_components)); - fe_v_face.get_function_values(solution, this_value); - - // Same for neighbor cell - std::vector> this_value_n( - n_q_face, Vector(n_components)); - fe_v_face_n.get_function_values(solution, this_value_n); - - for (const auto q_point : fe_v_face.quadrature_point_indices()) - { - Tensor<1, 2> vn = fe_v_face.normal_vector(q_point); - double nx = vn[0]; - double ny = vn[1]; - - double u = this_value[q_point](0); - double v = this_value[q_point](1); - - double u_n = this_value_n[q_point](0); - double v_n = this_value_n[q_point](1); - double un1 = u * nx + v * ny, un2 = u_n * nx + v_n * ny; - - deallog << "QP " << q_point - << ", Error: " << (u - u_n) * nx + (v - v_n) * ny + deallog << "QP " << q_point << ", Error: " << (u - u_n) * n << ", u " << un1 << ", un " << un2 << ", Rat " << un2 / un1 << std::endl; - Assert(std::fabs((u - u_n) * nx + (v - v_n) * ny) < 1e-12, - ExcInternalError()); + Assert(std::fabs((u - u_n) * n) < 1e-12, ExcInternalError()); } } } @@ -248,20 +163,17 @@ EvaluateNormal(DoFHandler<2> *dof_handler, Vector &solution) } + int -main(int /*argc*/, char ** /*argv*/) +main() { initlog(); deallog << std::setprecision(PRECISION); deallog << std::fixed; - Triangulation<2> tria_test; - DoFHandler<2> * dof_handler; - Point<2> p1(0, 0), p2(1, 1); - std::vector sub_div; - - sub_div.push_back(1); - sub_div.push_back(4); + Triangulation<2> tria_test; + const Point<2> p1(0, 0), p2(1, 1); + const std::vector sub_div = {1, 4}; GridGenerator::subdivided_hyper_rectangle(tria_test, sub_div, p1, p2); tria_test.refine_global(2); @@ -269,21 +181,17 @@ main(int /*argc*/, char ** /*argv*/) // Create a DoFHandler FE_RaviartThomas<2> fe(1); - dof_handler = new DoFHandler<2>(tria_test); - dof_handler->distribute_dofs(fe); + DoFHandler<2> dof_handler(tria_test); + dof_handler.distribute_dofs(fe); // Alloc some DoFs - Vector solution; - solution.reinit(dof_handler->n_dofs()); + Vector solution(dof_handler.n_dofs()); // Fill solution vector with random values between 0 and 1. - for (unsigned int i = 0; i < dof_handler->n_dofs(); ++i) + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) solution(i) = random_value(); // Now check if the function is continuous in normal // direction. - EvaluateNormal2(dof_handler, solution); - - delete (dof_handler); - return (0); + evaluate_normal_component(dof_handler, solution); } -- 2.39.5