From 604ba67c7edff1f5696c127364510015ea90253d Mon Sep 17 00:00:00 2001 From: Winnifried Wollner Date: Fri, 28 Feb 2020 10:26:44 +0100 Subject: [PATCH] New Test-Case for Curl-Conforming BV and FESystem --- tests/numerics/project_bv_curl_conf_04.cc | 250 ++++++++++++++++++ tests/numerics/project_bv_curl_conf_04.output | 0 2 files changed, 250 insertions(+) create mode 100644 tests/numerics/project_bv_curl_conf_04.cc create mode 100644 tests/numerics/project_bv_curl_conf_04.output diff --git a/tests/numerics/project_bv_curl_conf_04.cc b/tests/numerics/project_bv_curl_conf_04.cc new file mode 100644 index 0000000000..0e19fbc650 --- /dev/null +++ b/tests/numerics/project_bv_curl_conf_04.cc @@ -0,0 +1,250 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace dealii; + +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) +{ + 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); + 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) +{ + double ret = true; + //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())); + const FEValuesExtractors::Vector nedelec(start_comp); + + //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(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 "< triangulation; + GridGenerator::hyper_cube(triangulation, 0, 1, true); + MappingQ mapping(1); + AffineConstraints< double > constraints; + Vector test_vec; + Vector test_vec2; + + { + //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); + 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); + 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); + 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); + 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); + 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); + 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"<