]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Reindent.. 1666/head
authorAlexander Grayver <agrayver@erdw.ethz.ch>
Thu, 24 Sep 2015 14:51:23 +0000 (16:51 +0200)
committerAlexander Grayver <agrayver@erdw.ethz.ch>
Thu, 24 Sep 2015 14:51:23 +0000 (16:51 +0200)
tests/fe/fe_project_3d.cc

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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.