From ff4aa557ed4bcd49d873987f3b880ff22f39ac02 Mon Sep 17 00:00:00 2001 From: Alexander Grayver Date: Thu, 24 Sep 2015 16:51:23 +0200 Subject: [PATCH] Reindent.. --- tests/fe/fe_project_3d.cc | 278 ++++++++++++++++++++------------------ 1 file changed, 144 insertions(+), 134 deletions(-) diff --git a/tests/fe/fe_project_3d.cc b/tests/fe/fe_project_3d.cc index f058bfa3db..a96d88a14f 100644 --- a/tests/fe/fe_project_3d.cc +++ b/tests/fe/fe_project_3d.cc @@ -32,18 +32,19 @@ /* * This program projects a function into FE spaces defined on meshes * consisting of: rectangular cells, affine cells, non-affine cells - * + * * The error, curl and divergence are then numerically calculated * on a series of globally refined meshes and output. - * + * * Among FE spaces tested are FE_Nedelec and FE_RaviartThomas - * + * * Alexander Grayver, Maien Hamed */ using namespace dealii; -static const Point<3> vertices_affine[] = { +static const Point<3> vertices_affine[] = +{ Point<3> (-1., -1., -1.), Point<3> (0., -1., -1.), Point<3> (1., -1., -1.), @@ -81,7 +82,8 @@ static const Point<3> vertices_affine[] = { Point<3> (0.6, 1., 1) }; -static const Point<3> vertices_nonaffine[] = { +static const Point<3> vertices_nonaffine[] = +{ Point<3> (-1., -1., -1.), Point<3> (0., -1., -1.), Point<3> (1., -1., -1.), @@ -119,7 +121,8 @@ static const Point<3> vertices_nonaffine[] = { Point<3> (1., 1., 1.) }; -static const Point<3> vertices_rectangular[] = { +static const Point<3> vertices_rectangular[] = +{ Point<3> (-1., -1., -1.), Point<3> (0., -1., -1.), Point<3> (1., -1., -1.), @@ -167,7 +170,7 @@ public: VectorFunction() : Function(dim) {} virtual double value (const Point &p, const unsigned int component) const; virtual void vector_value(const Point &p, Vector &values) const; - virtual Tensor< 1, dim > gradient (const Point< dim > &p, const unsigned int component=0) const; + virtual Tensor< 1, dim > gradient (const Point< dim > &p, const unsigned int component=0) const; }; template<> @@ -177,19 +180,25 @@ double VectorFunction<3>::value(const Point<3> &p, const unsigned int component) const double PI = numbers::PI; double val = 0.0; - switch(component) - { - case 0: val = -sin(PI*p(0))*cos(PI*p(1))*cos(PI*p(2)); break; - case 1: val = -cos(PI*p(0))*sin(PI*p(1))*cos(PI*p(2)); break; - case 2: val = 2*cos(PI*p(0))*cos(PI*p(1))*sin(PI*p(2)); break; - } + switch (component) + { + case 0: + val = -sin(PI*p(0))*cos(PI*p(1))*cos(PI*p(2)); + break; + case 1: + val = -cos(PI*p(0))*sin(PI*p(1))*cos(PI*p(2)); + break; + case 2: + val = 2*cos(PI*p(0))*cos(PI*p(1))*sin(PI*p(2)); + break; + } return val; } template void VectorFunction::vector_value(const Point &p, Vector &values) const { - for(int i = 0; i < dim; ++i) + for (int i = 0; i < dim; ++i) values(i) = value(p, i); } @@ -200,34 +209,35 @@ Tensor<1, 3> VectorFunction<3>::gradient(const Point<3> &p, const unsigned int c Tensor<1, 3> val; double x = p(0), y = p(1), z = p(2); - switch(component) - { - case 0: - val[0] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); - val[1] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); - val[2] = -2*PI*cos(PI*y)*sin(PI*x)*sin(PI*z); - break; - case 1: - val[0] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); - val[1] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); - val[2] = -2*PI*cos(PI*x)*sin(PI*y)*sin(PI*z); - break; - case 2: - val[0] = PI*cos(PI*y)*sin(PI*x)*sin(PI*z); - val[1] = PI*cos(PI*x)*sin(PI*y)*sin(PI*z); - val[2] = 2*PI*cos(PI*x)*cos(PI*y)*cos(PI*z); - break; - } + switch (component) + { + case 0: + val[0] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + val[1] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); + val[2] = -2*PI*cos(PI*y)*sin(PI*x)*sin(PI*z); + break; + case 1: + val[0] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); + val[1] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + val[2] = -2*PI*cos(PI*x)*sin(PI*y)*sin(PI*z); + break; + case 2: + val[0] = PI*cos(PI*y)*sin(PI*x)*sin(PI*z); + val[1] = PI*cos(PI*x)*sin(PI*y)*sin(PI*z); + val[2] = 2*PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + break; + } return val; } -void create_tria(Triangulation<3>& triangulation, const Point<3> *vertices_parallelograms) +void create_tria(Triangulation<3> &triangulation, const Point<3> *vertices_parallelograms) { const std::vector > vertices (&vertices_parallelograms[0], - &vertices_parallelograms[n_vertices]); + &vertices_parallelograms[n_vertices]); // create grid with all possible combintations of face_flip, face_orientation and face_rotation flags - static const int cell_vertices[][GeometryInfo<3>::vertices_per_cell] = { + static const int cell_vertices[][GeometryInfo<3>::vertices_per_cell] = + { {0, 1, 9, 10, 3, 4, 12, 13}, // cell 1 standard {1, 2, 10, 11, 4, 5, 13, 14}, // cell 2 standard //{10, 11, 13, 14, 1, 2, 4, 5}, // cell 2 rotated by 270 deg @@ -243,17 +253,17 @@ void create_tria(Triangulation<3>& triangulation, const Point<3> *vertices_paral std::vector > cells (n_cells, CellData<3>()); for (unsigned i = 0; i::vertices_per_cell; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - } + { + for (unsigned int j=0; j::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } triangulation.create_triangulation (vertices, cells, SubCellData()); } template -void test(const FiniteElement& fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) +void test(const FiniteElement &fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) { deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; deallog << "DoFs\t\t||u-u_h||_1\tcurl(u_h)\ttangentials\tcurl(curl(u_h))\tcurl_curl_traces\tdiv(u_h)\tboundary_flux" << std::endl; @@ -278,112 +288,112 @@ void test(const FiniteElement& fe, unsigned n_cycles, bool global, const Po std::vector > face_values (n_face_q_points); std::vector::curl_type> face_curls (n_face_q_points); - for(unsigned cycle = 0; cycle < n_cycles; ++cycle) - { - dof_handler.distribute_dofs(fe); + for (unsigned cycle = 0; cycle < n_cycles; ++cycle) + { + dof_handler.distribute_dofs(fe); - ConstraintMatrix constraints; - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - constraints.close(); + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); - Vector v(dof_handler.n_dofs()); - VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); + Vector v(dof_handler.n_dofs()); + VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); - Vector diff(triangulation.n_active_cells()); - VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, - QGauss(fe.degree + 2), VectorTools::L1_norm); + Vector diff(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, + QGauss(fe.degree + 2), VectorTools::L1_norm); - typename FEValuesViews::Vector::curl_type total_curl, boundary_tangentials; - Tensor<1, dim> total_curl_curl, boundary_curl_curl_traces; - double total_div = 0; - double boundary_flux = 0; - total_curl *= 0.; - boundary_tangentials *= 0.; + typename FEValuesViews::Vector::curl_type total_curl, boundary_tangentials; + Tensor<1, dim> total_curl_curl, boundary_curl_curl_traces; + double total_div = 0; + double boundary_flux = 0; + total_curl *= 0.; + boundary_tangentials *= 0.; - FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_hessians); - FEFaceValues fe_face_values(mapping, fe, face_quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_normal_vectors ); - unsigned int cell_index = 0; + FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_hessians); + FEFaceValues fe_face_values(mapping, fe, face_quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_normal_vectors ); + unsigned int cell_index = 0; - for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); - cell != dof_handler.end (); ++cell, ++cell_index) - { - fe_values.reinit (cell); - const std::vector& JxW_values = fe_values.get_JxW_values (); - fe_values[vec].get_function_divergences (v, div_v); - fe_values[vec].get_function_curls (v, curl_v); - fe_values[vec].get_function_hessians (v, hessians); - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + cell != dof_handler.end (); ++cell, ++cell_index) { - total_div += JxW_values[q_point] * div_v[q_point]; - total_curl += JxW_values[q_point] * curl_v[q_point]; - if(dim == 3) - { - total_curl_curl[0] += JxW_values[q_point] * (hessians[q_point][1][0][1] + hessians[q_point][2][0][2] - hessians[q_point][0][1][1] - hessians[q_point][0][2][2]); - total_curl_curl[1] += JxW_values[q_point] * (hessians[q_point][2][1][2] + hessians[q_point][0][0][1] - hessians[q_point][1][2][2] - hessians[q_point][1][0][0]); - total_curl_curl[2] += JxW_values[q_point] * (hessians[q_point][0][0][2] + hessians[q_point][1][1][2] - hessians[q_point][2][0][0] - hessians[q_point][2][1][1]); - } - } + fe_values.reinit (cell); + const std::vector &JxW_values = fe_values.get_JxW_values (); + fe_values[vec].get_function_divergences (v, div_v); + fe_values[vec].get_function_curls (v, curl_v); + fe_values[vec].get_function_hessians (v, hessians); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + total_div += JxW_values[q_point] * div_v[q_point]; + total_curl += JxW_values[q_point] * curl_v[q_point]; + if (dim == 3) + { + total_curl_curl[0] += JxW_values[q_point] * (hessians[q_point][1][0][1] + hessians[q_point][2][0][2] - hessians[q_point][0][1][1] - hessians[q_point][0][2][2]); + total_curl_curl[1] += JxW_values[q_point] * (hessians[q_point][2][1][2] + hessians[q_point][0][0][1] - hessians[q_point][1][2][2] - hessians[q_point][1][0][0]); + total_curl_curl[2] += JxW_values[q_point] * (hessians[q_point][0][0][2] + hessians[q_point][1][1][2] - hessians[q_point][2][0][0] - hessians[q_point][2][1][1]); + } + } - for (unsigned int face=0; face::faces_per_cell; ++face) - { - fe_face_values.reinit(cell,face); - const std::vector& face_JxW_values = fe_face_values.get_JxW_values (); - fe_face_values[vec].get_function_values (v, face_values); - if(dim==3) fe_face_values[vec].get_function_curls (v, face_curls); - for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) - { - const Tensor<1,dim> &normal = fe_face_values.normal_vector(q_point); - - // boundary flux - if(cell->at_boundary(face)) - boundary_flux += face_JxW_values[q_point] * (face_values[q_point] * normal); - else - total_div -= face_JxW_values[q_point] * (face_values[q_point] * normal); - - // boundary tangentials (curl traces) - typename FEValuesViews::Vector::curl_type n_x_v; - if(dim==2) - n_x_v[0] = (-normal[1]*face_values[q_point][0] + normal[0]*face_values[q_point][1]); - else if(dim==3) - cross_product(*reinterpret_cast*>(&n_x_v), normal, face_values[q_point]); - - if(cell->at_boundary(face)) - boundary_tangentials += face_JxW_values[q_point] * n_x_v; - else - total_curl -= face_JxW_values[q_point] * n_x_v; - - // boundary curl curl traces - if(dim==3) + for (unsigned int face=0; face::faces_per_cell; ++face) { - Tensor<1,dim> n_x_curl_u; - cross_product(n_x_curl_u, normal, *reinterpret_cast*>(&face_curls[q_point])); - if(cell->at_boundary(face)) - boundary_curl_curl_traces += face_JxW_values[q_point] * n_x_curl_u; - else - total_curl_curl -= face_JxW_values[q_point] * n_x_curl_u; + fe_face_values.reinit(cell,face); + const std::vector &face_JxW_values = fe_face_values.get_JxW_values (); + fe_face_values[vec].get_function_values (v, face_values); + if (dim==3) fe_face_values[vec].get_function_curls (v, face_curls); + for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) + { + const Tensor<1,dim> &normal = fe_face_values.normal_vector(q_point); + + // boundary flux + if (cell->at_boundary(face)) + boundary_flux += face_JxW_values[q_point] * (face_values[q_point] * normal); + else + total_div -= face_JxW_values[q_point] * (face_values[q_point] * normal); + + // boundary tangentials (curl traces) + typename FEValuesViews::Vector::curl_type n_x_v; + if (dim==2) + n_x_v[0] = (-normal[1]*face_values[q_point][0] + normal[0]*face_values[q_point][1]); + else if (dim==3) + cross_product(*reinterpret_cast*>(&n_x_v), normal, face_values[q_point]); + + if (cell->at_boundary(face)) + boundary_tangentials += face_JxW_values[q_point] * n_x_v; + else + total_curl -= face_JxW_values[q_point] * n_x_v; + + // boundary curl curl traces + if (dim==3) + { + Tensor<1,dim> n_x_curl_u; + cross_product(n_x_curl_u, normal, *reinterpret_cast*>(&face_curls[q_point])); + if (cell->at_boundary(face)) + boundary_curl_curl_traces += face_JxW_values[q_point] * n_x_curl_u; + else + total_curl_curl -= face_JxW_values[q_point] * n_x_curl_u; + } + } } - } } - } - deallog << dof_handler.n_dofs() << "\t\t" - << diff.l1_norm() << "\t" - << total_curl.norm() << "\t" - << boundary_tangentials.norm() << "\t" - << total_curl_curl.norm() << "\t" - << boundary_curl_curl_traces.norm() << "\t" - << total_div << "\t" - << boundary_flux << std::endl; - - if(global) - triangulation.refine_global(); - else - { - GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); - triangulation.prepare_coarsening_and_refinement(); - triangulation.execute_coarsening_and_refinement(); + deallog << dof_handler.n_dofs() << "\t\t" + << diff.l1_norm() << "\t" + << total_curl.norm() << "\t" + << boundary_tangentials.norm() << "\t" + << total_curl_curl.norm() << "\t" + << boundary_curl_curl_traces.norm() << "\t" + << total_div << "\t" + << boundary_flux << std::endl; + + if (global) + triangulation.refine_global(); + else + { + GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement(); + } } - } } int main () -- 2.39.5