From c0c9f28d025c71ea9a3d40611c95d1f66516c91d Mon Sep 17 00:00:00 2001 From: Winnifried Wollner Date: Fri, 28 Feb 2020 17:29:49 +0100 Subject: [PATCH] indentation --- .../deal.II/numerics/vector_tools.templates.h | 149 ++--- tests/numerics/project_bv_curl_conf_04.cc | 569 ++++++++++-------- 2 files changed, 381 insertions(+), 337 deletions(-) diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 415952dee0..5d04209f72 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -5406,15 +5406,15 @@ namespace VectorTools fe.base_element(i).n_components(); } else - { - //Assert that the FE is in fact an FE_Nedelec, so that the default - //base_indices == (0,0) is correct. - Assert((dynamic_cast *>(&cell->get_fe()) != - nullptr) || - (dynamic_cast *>(&cell->get_fe()) != - nullptr), - ExcNotImplemented()); - } + { + // Assert that the FE is in fact an FE_Nedelec, so that the default + // base_indices == (0,0) is correct. + Assert((dynamic_cast *>(&cell->get_fe()) != + nullptr) || + (dynamic_cast *>(&cell->get_fe()) != + nullptr), + ExcNotImplemented()); + } // Store degree as fe.degree-1 // For nedelec elements FE_Nedelec (0) returns fe.degree = 1. // For FESystem get the degree from the base_element @@ -5438,7 +5438,8 @@ namespace VectorTools // element and the index within this ordering. // // We call the map associated_edge_dof_to_face_dof - std::vector associated_edge_dof_to_face_dof(degree + 1, numbers::invalid_dof_index); + std::vector associated_edge_dof_to_face_dof( + degree + 1, numbers::invalid_dof_index); // Lowest DoF in the base element allowed for this edge: const unsigned int lower_bound = @@ -5450,64 +5451,65 @@ namespace VectorTools .face_to_cell_index((line + 1) * (degree + 1) - 1, face); unsigned int associated_edge_dof_index = 0; - for (unsigned int line_dof_idx = 0; line_dof_idx < fe.dofs_per_line; ++line_dof_idx) + for (unsigned int line_dof_idx = 0; line_dof_idx < fe.dofs_per_line; + ++line_dof_idx) { // For each DoF associated with the (interior of) the line, we need - // to figure out which base element it belongs to and then if - // that's the correct base element. This is complicated by the - // fact that the FiniteElement class has functions that translate - // from face to cell, but not from edge to cell index systems. So - // we have to do that step by step. - // - // DoFs on a face in 3d are numbered in order by vertices then lines - // then faces. + // to figure out which base element it belongs to and then if + // that's the correct base element. This is complicated by the + // fact that the FiniteElement class has functions that translate + // from face to cell, but not from edge to cell index systems. So + // we have to do that step by step. + // + // DoFs on a face in 3d are numbered in order by vertices then lines + // then faces. // i.e. line 0 has degree+1 dofs numbered 0,..,degree // line 1 has degree+1 dofs numbered (degree+1),..,2*(degree+1) // and so on. - const unsigned int face_dof_idx = - GeometryInfo::vertices_per_face * fe.dofs_per_vertex + - line * fe.dofs_per_line + line_dof_idx; + const unsigned int face_dof_idx = + GeometryInfo::vertices_per_face * fe.dofs_per_vertex + + line * fe.dofs_per_line + line_dof_idx; // Note, assuming that the edge orientations are "standard" // i.e. cell->line_orientation(line) = true. - Assert(cell->line_orientation(line), - ExcMessage("Edge orientation doesnot meet expectation.")); - // Next, translate from face to cell. Note, this might be assuming + Assert(cell->line_orientation(line), + ExcMessage("Edge orientation doesnot meet expectation.")); + // Next, translate from face to cell. Note, this might be assuming // that the edge orientations are "standard" (not sure any more at // this time), i.e. // cell->line_orientation(line) = true. const unsigned int cell_dof_idx = fe.face_to_cell_index(face_dof_idx, face); - // Check this cell_dof_idx belongs to the correct base_element, component - // and line: We do this for each of the supported elements - // separately - bool dof_is_of_interest = false; + // Check this cell_dof_idx belongs to the correct base_element, + // component and line: We do this for each of the supported elements + // separately + bool dof_is_of_interest = false; if (dynamic_cast *>(&fe) != nullptr) - { - dof_is_of_interest = - (fe.system_to_base_index(cell_dof_idx).first == base_indices) && - (lower_bound <= fe.system_to_base_index(cell_dof_idx).second) && - (fe.system_to_base_index(cell_dof_idx).second <= upper_bound); - } - else if ((dynamic_cast *>(&fe) != nullptr) || - (dynamic_cast *>(&fe) != nullptr)) - { - Assert((line * (degree + 1) <= face_dof_idx) && - (face_dof_idx < (line + 1) * (degree + 1)), - ExcInternalError()); - dof_is_of_interest = true; - } - else - Assert(false, ExcNotImplemented()); - - if (dof_is_of_interest) - { - associated_edge_dof_to_face_dof[associated_edge_dof_index] = - face_dof_idx; - ++associated_edge_dof_index; - } + { + dof_is_of_interest = + (fe.system_to_base_index(cell_dof_idx).first == base_indices) && + (lower_bound <= fe.system_to_base_index(cell_dof_idx).second) && + (fe.system_to_base_index(cell_dof_idx).second <= upper_bound); + } + else if ((dynamic_cast *>(&fe) != nullptr) || + (dynamic_cast *>(&fe) != nullptr)) + { + Assert((line * (degree + 1) <= face_dof_idx) && + (face_dof_idx < (line + 1) * (degree + 1)), + ExcInternalError()); + dof_is_of_interest = true; + } + else + Assert(false, ExcNotImplemented()); + + if (dof_is_of_interest) + { + associated_edge_dof_to_face_dof[associated_edge_dof_index] = + face_dof_idx; + ++associated_edge_dof_index; + } } // Sanity check: const unsigned int n_associated_edge_dofs = associated_edge_dof_index; @@ -5706,15 +5708,15 @@ namespace VectorTools fe.base_element(i).n_components(); } else - { - //Assert that the FE is in fact an FE_Nedelec, so that the default - //base_indices == (0,0) is correct. - Assert((dynamic_cast *>(&cell->get_fe()) != - nullptr) || - (dynamic_cast *>(&cell->get_fe()) != - nullptr), - ExcNotImplemented()); - } + { + // Assert that the FE is in fact an FE_Nedelec, so that the default + // base_indices == (0,0) is correct. + Assert((dynamic_cast *>(&cell->get_fe()) != + nullptr) || + (dynamic_cast *>(&cell->get_fe()) != + nullptr), + ExcNotImplemented()); + } const unsigned int degree = fe.base_element(base_indices.first).degree - 1; @@ -5881,10 +5883,11 @@ namespace VectorTools fe.base_element(base_indices.first) .face_to_cell_index((line + 1) * (degree + 1) - 1, face); unsigned int associated_edge_dof_index = 0; - for (unsigned int line_dof_idx = 0; line_dof_idx < fe.dofs_per_line; + for (unsigned int line_dof_idx = 0; + line_dof_idx < fe.dofs_per_line; ++line_dof_idx) { - // For each DoF associated with the (interior of) the + // For each DoF associated with the (interior of) the // line, we need to figure out which base element it // belongs to and then if that's the correct base element. // This is complicated by the fact that the FiniteElement @@ -5902,9 +5905,9 @@ namespace VectorTools fe.dofs_per_vertex + line * fe.dofs_per_line + line_dof_idx; - const unsigned int cell_dof_idx = + const unsigned int cell_dof_idx = fe.face_to_cell_index(face_dof_idx, face); - // Check that this cell_idx belongs to the correct + // Check that this cell_idx belongs to the correct // base_element, component and line. We do this for each // of the supported elements separately bool dof_is_of_interest = false; @@ -5931,7 +5934,7 @@ namespace VectorTools else Assert(false, ExcNotImplemented()); - if (dof_is_of_interest) + if (dof_is_of_interest) { associated_edge_dof_to_face_dof [line][associated_edge_dof_index] = face_dof_idx; @@ -5965,12 +5968,12 @@ namespace VectorTools // Loop over the quad-interior dofs. unsigned int associated_face_dof_index = 0; for (unsigned int quad_dof_idx = 0; - quad_dof_idx < fe.dofs_per_quad; - ++quad_dof_idx) + quad_dof_idx < fe.dofs_per_quad; + ++quad_dof_idx) { - const unsigned int face_idx = - GeometryInfo::vertices_per_face * fe.dofs_per_vertex + - lines_per_face * fe.dofs_per_line + quad_dof_idx; + const unsigned int face_idx = + GeometryInfo::vertices_per_face * fe.dofs_per_vertex + + lines_per_face * fe.dofs_per_line + quad_dof_idx; const unsigned int cell_idx = fe.face_to_cell_index(face_idx, face); if (((dynamic_cast *>(&fe) != nullptr) && @@ -5979,8 +5982,8 @@ namespace VectorTools (dynamic_cast *>(&fe) != nullptr) || (dynamic_cast *>(&fe) != nullptr)) { - AssertIndexRange(associated_face_dof_index, - associated_face_dof_to_face_dof.size()); + AssertIndexRange(associated_face_dof_index, + associated_face_dof_to_face_dof.size()); associated_face_dof_to_face_dof [associated_face_dof_index] = face_idx; ++associated_face_dof_index; diff --git a/tests/numerics/project_bv_curl_conf_04.cc b/tests/numerics/project_bv_curl_conf_04.cc index d4e7c128f6..8f93314a55 100644 --- a/tests/numerics/project_bv_curl_conf_04.cc +++ b/tests/numerics/project_bv_curl_conf_04.cc @@ -1,378 +1,419 @@ -#include -#include -#include +#include + #include -#include +#include +#include #include +#include + +#include + #include -#include + #include using namespace dealii; -const unsigned int dim=3; +const unsigned int dim = 3; -void apply_boundary_values(DoFHandler& dof_handler, AffineConstraints & constraints, unsigned int n_comps, unsigned int start_comp, Mapping & mapping, Vector &dst) +void +apply_boundary_values(DoFHandler & dof_handler, + AffineConstraints &constraints, + unsigned int n_comps, + unsigned int start_comp, + Mapping & mapping, + Vector & dst) { constraints.clear(); - for(unsigned int i = 0; i < 6; i++) - VectorTools::project_boundary_values_curl_conforming_l2(dof_handler, start_comp, - Functions::ConstantFunction(1, n_comps), i, constraints, mapping); + for (unsigned int i = 0; i < 6; i++) + VectorTools::project_boundary_values_curl_conforming_l2( + dof_handler, + start_comp, + Functions::ConstantFunction(1, n_comps), + i, + constraints, + mapping); constraints.close(); - + constraints.distribute(dst); } -bool test_boundary_values(DoFHandler& dof_handler, Mapping & mapping, FiniteElement & fe, unsigned int n_comps, unsigned int start_comp, Vector &vec) +bool +test_boundary_values(DoFHandler & dof_handler, + Mapping & mapping, + FiniteElement &fe, + unsigned int n_comps, + unsigned int start_comp, + Vector & vec) { double ret = true; - //Initialize - QGaussLobatto quadrature(3); - FEFaceValues fe_values(mapping, fe, quadrature, - update_values | update_normal_vectors - | update_quadrature_points | update_JxW_values); + // Initialize + QGaussLobatto quadrature(3); + FEFaceValues fe_values(mapping, + fe, + quadrature, + update_values | update_normal_vectors | + update_quadrature_points | update_JxW_values); typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); - unsigned num_cells = 0; - std::vector > local_averages(6,Vector(dof_handler.n_dofs())); + dof_handler.begin_active(), + endc = dof_handler.end(); + unsigned num_cells = 0; + std::vector> local_averages( + 6, Vector(dof_handler.n_dofs())); const FEValuesExtractors::Vector nedelec(start_comp); - //For testing we take only one cell! + // For testing we take only one cell! assert(dof_handler.n_dofs() == cell->get_fe().dofs_per_cell); - - for (; cell != endc; ++cell) { - //For testing we take only one cell! Otherwise, local to global is missing - assert(num_cells == 0); - for (unsigned int face = 0; - face < dealii::GeometryInfo::faces_per_cell; ++face) { - if (cell->face(face)->at_boundary()) - { - unsigned int boundary_number = cell->face(face)->boundary_id(); - Tensor<1,dim> expected_value; - expected_value[0] = 1; - expected_value[1] = 1; - expected_value[2] = 1; - assert(boundary_number < 6); - fe_values.reinit(cell,face); - std::vector > local_values(fe_values.n_quadrature_points,Vector(n_comps)); - fe_values.get_function_values(vec,local_values); + for (; cell != endc; ++cell) + { + // For testing we take only one cell! Otherwise, local to global is + // missing + assert(num_cells == 0); + for (unsigned int face = 0; + face < dealii::GeometryInfo::faces_per_cell; + ++face) + { + if (cell->face(face)->at_boundary()) + { + unsigned int boundary_number = cell->face(face)->boundary_id(); + Tensor<1, dim> expected_value; + expected_value[0] = 1; + expected_value[1] = 1; + expected_value[2] = 1; + assert(boundary_number < 6); + fe_values.reinit(cell, face); - for(unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; q_point ++) - { - //Compare n x v values - Tensor<1,dim> nedelec_values; - nedelec_values[0] = local_values[q_point](start_comp); - nedelec_values[1] = local_values[q_point](start_comp+1); - nedelec_values[2] = local_values[q_point](start_comp+2); - Tensor<1,dim> normal = fe_values.normal_vector(q_point); - - for(unsigned int i = 0; i < cell->get_fe().dofs_per_cell; i++) - { - local_averages[boundary_number](i) += ((cross_product_3d(normal,nedelec_values)-cross_product_3d(normal,expected_value))*cross_product_3d(normal,fe_values[nedelec].value(i,q_point)))*fe_values.JxW(q_point); - } - } - } - } - num_cells++; - } - //Test if ok - for(unsigned int bc=0; bc < 6; bc++) - { - for(unsigned int basis=0; basis 1.e-10) - { - std::cout<<"Wrong values on boundary "<> local_values( + fe_values.n_quadrature_points, Vector(n_comps)); + fe_values.get_function_values(vec, local_values); + + for (unsigned int q_point = 0; + q_point < fe_values.n_quadrature_points; + q_point++) + { + // Compare n x v values + Tensor<1, dim> nedelec_values; + nedelec_values[0] = local_values[q_point](start_comp); + nedelec_values[1] = local_values[q_point](start_comp + 1); + nedelec_values[2] = local_values[q_point](start_comp + 2); + Tensor<1, dim> normal = fe_values.normal_vector(q_point); + + for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; + i++) + { + local_averages[boundary_number](i) += + ((cross_product_3d(normal, nedelec_values) - + cross_product_3d(normal, expected_value)) * + cross_product_3d( + normal, fe_values[nedelec].value(i, q_point))) * + fe_values.JxW(q_point); + } + } + } + } + num_cells++; + } + // Test if ok + for (unsigned int bc = 0; bc < 6; bc++) + { + for (unsigned int basis = 0; basis < dof_handler.n_dofs(); basis++) + { + if (fabs(local_averages[bc](basis)) > 1.e-10) + { + std::cout << "Wrong values on boundary " << bc + << " in basis function " << basis << " of size " + << fabs(local_averages[bc](basis)) << std::endl; + ret = false; + } + } } - } return ret; } -int main(int /*argc*/, char **/*argv*/) { - +int +main(int /*argc*/, char ** /*argv*/) +{ Triangulation triangulation; GridGenerator::hyper_cube(triangulation, 0, 1, true); - MappingQ mapping(1); - AffineConstraints< double > constraints; - Vector test_vec; - Vector test_vec2; - + MappingQ mapping(1); + AffineConstraints constraints; + Vector test_vec; + Vector test_vec2; + { - //Test 1: Only FE_Nedelec - std::cout<<"Testing FE_Nedelec(0): "; + // Test 1: Only FE_Nedelec + std::cout << "Testing FE_Nedelec(0): "; FE_Nedelec fe(0); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe); test_vec2.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,3,0,mapping,test_vec2); - if(test_boundary_values(dof_handler,mapping,fe,3,0,test_vec2)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(0),1); + // Test 2: FESystem(FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Nedelec(0)): "; + FESystem fe_system(FE_Nedelec(0), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,3,0,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,3,0,test_vec)) - { - std::cout<<"OK"< fe(1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe); test_vec2.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,3,0,mapping,test_vec2); - if(test_boundary_values(dof_handler,mapping,fe,3,0,test_vec2)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(1),1); + // Test 2b: FESystem(FE_Nedelec(1)) + std::cout << "Testing FESystem(FE_Nedelec(1)): "; + FESystem fe_system(FE_Nedelec(1), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,3,0,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,3,0,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(0),1,FE_Nedelec(0),1); + // Test 3: FESystem(FE_Nedelec(0), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Nedelec(0)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Nedelec(0), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,6,3,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,6,3,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(1),1,FE_Nedelec(0),1); + // Test 4: FESystem(FE_Nedelec(1), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Nedelec(1), FE_Nedelec(0)): "; + FESystem fe_system(FE_Nedelec(1), 1, FE_Nedelec(0), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,6,3,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,6,3,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(0),1,FE_Nedelec(1),1); + // Test 5: FESystem(FE_Nedelec(0), FE_Nedelec(1)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Nedelec(1)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Nedelec(1), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,6,3,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,6,3,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(0),1,FE_Q(1),1); + // Test 6: FESystem(FE_Nedelec(0), FE_Q(1)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Q(1)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Q(1), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,4,0,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,4,0,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Q(1),1,FE_Nedelec(0),1); + // Test 7: FESystem(FE_Q(1), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Q(1), FE_Nedelec(0)): "; + FESystem fe_system(FE_Q(1), 1, FE_Nedelec(0), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,4,1,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,4,1,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(0),1,FE_Q(2),1); + // Test 8: FESystem(FE_Nedelec(0), FE_Q(2)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Q(2)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Q(2), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,4,0,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,4,0,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Q(2),1,FE_Nedelec(0),1); + // Test 9: FESystem(FE_Q(2), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Q(2), FE_Nedelec(0)): "; + FESystem fe_system(FE_Q(2), 1, FE_Nedelec(0), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,4,1,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,4,1,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Nedelec(1),1,FE_Q(1),1); + { + // Test 10: FESystem(FE_Nedelec(1), FE_Q(1)) + std::cout << "Testing FESystem(FE_Nedelec(1), FE_Q(1)): "; + FESystem fe_system(FE_Nedelec(1), 1, FE_Q(1), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,4,0,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,4,0,test_vec)) - { - std::cout<<"OK"< fe_system(FE_Q(1),1,FE_Nedelec(1),1); + // Test 11: FESystem(FE_Q(1), FE_Nedelec(1)) + std::cout << "Testing FESystem(FE_Q(1), FE_Nedelec(1)): "; + FESystem fe_system(FE_Q(1), 1, FE_Nedelec(1), 1); DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(fe_system); test_vec.reinit(dof_handler.n_dofs()); - apply_boundary_values(dof_handler,constraints,4,1,mapping,test_vec); - if(test_boundary_values(dof_handler,mapping,fe_system,4,1,test_vec)) - { - std::cout<<"OK"<