]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Testsuite: Clean up codim_one/bem; remove unused code fragments 1617/head
authorMatthias Maier <tamiko@43-1.org>
Tue, 15 Sep 2015 20:33:58 +0000 (15:33 -0500)
committerMatthias Maier <tamiko@43-1.org>
Tue, 15 Sep 2015 20:37:12 +0000 (15:37 -0500)
tests/codim_one/bem.cc

index 7b28ad6d92a2c3ee2250d800aa7929c71742596a..ea3dc570730d8646988fd8d9b09d8847cdb99ed5 100644 (file)
@@ -79,10 +79,8 @@ private:
   FE_DGQ<spacedim-1, spacedim>          fe;
   DoFHandler<spacedim-1, spacedim>      dof_handler;
 
-  // finite elements used to smoothen
-  // the solution (from piecewise
-  // constant to continuous piecewise
-  // quadratic)
+  // finite elements used to smoothen the solution (from piecewise constant
+  // to continuous piecewise quadratic)
   FE_Q<spacedim-1, spacedim>            fe_q;
   DoFHandler<spacedim-1, spacedim>      dof_handler_q;
 
@@ -153,64 +151,9 @@ BEM<spacedim>::run()
 
           tria.set_boundary(1);
 
-
-//  std::cout
-//      << "solution"
-//      << std::endl;
-//  for(unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-//      std::cout
-//    << i<< " - "
-//    <<solution(i)<< std::endl;
-
-//  std::cout
-//      << "smooth solution"
-//      << std::endl;
-//  for(unsigned int i=0; i<dof_handler_q.n_dofs(); ++i)
-//      std::cout
-//    << i<< " - "
-//    <<smooth_solution(i)<< std::endl;
-
-//  std::cout
-//      << "tangential velocity"
-//      << std::endl;
-//  for(unsigned int i=0; i<dof_handler_q.n_dofs(); ++i)
-//      std::cout
-//    << i<< " - "
-//    <<tangential_derivative(i)<< std::endl;
-
-//  std::cout<<tria.n_active_cells()<<std::endl;
-
-//  std::cout
-//      << "cell error" << std::endl
-//      << "solution // smooth_sol // tang_deriv // tang_vel // exact_sol // error"
-//      << std::endl;
-//  for(unsigned int i=0; i<dof_handler_q.n_dofs(); ++i)
-//      std::cout
-//    << i<< " - "
-//    << solution(i) << " // "
-//    << smooth_solution(i)<< " // "
-//    << tangential_derivative(i) << " // "
-//    << tangential_velocity(i) << " // "
-//    // exact solution
-//    << 2. * velocity[0]
-//       * sin (
-//           numbers::PI/2.
-//           - (1./2.+i) * 2.*numbers::PI/tria.n_active_cells()
-//             )
-//    << " // "
-//    << error(i)
-//    << std::endl;
-
-//      std::cout
-//    << "global error"
-//    <<std::endl;
-
           for (unsigned int i=0; i<dof_handler_q.n_dofs(); ++i)
             global_error +=
               pow(fabs(error(i)),2) * side_length;
-//      std::cout
-//    << global_error
-//    << std::endl;
 
           table.add_value("L^2 norm error", global_error);
 
@@ -256,9 +199,8 @@ BEM<spacedim>::assemble_system()
   const unsigned int        dofs_per_cell = fe.dofs_per_cell;
   const unsigned int        n_q_points    = quadrature_formula.size();
 
-  // The following two matrices will
-  // be subtracted to obtain the
-  // system matrix.
+  // The following two matrices will be subtracted to obtain the system
+  // matrix.
   FullMatrix<double>        DLP_matrix(dof_handler.n_dofs()); // DLP stands for double layer potential
   FullMatrix<double>        mass_matrix(dof_handler.n_dofs());
 
@@ -271,14 +213,6 @@ BEM<spacedim>::assemble_system()
   std::vector<types::global_dof_index> local_dof_indices_i (dofs_per_cell);
   std::vector<types::global_dof_index> local_dof_indices_j (dofs_per_cell);
 
-
-//     std::cout
-//  << "no. of cells: "<< tria.n_active_cells()<< std::endl
-//  << "no. of dofs: "<< dof_handler.n_dofs()<< std::endl
-//  << "dofs per cell: "<< dofs_per_cell<< std::endl
-//  << "side length (r=10): "<< 10 * 2*sin(numbers::PI/tria.n_cells())<<std::endl ;
-
-
   typename DoFHandler<spacedim-1, spacedim>::active_cell_iterator
   cell_i = dof_handler.begin_active(),
   cell_j = dof_handler.begin_active(),
@@ -292,20 +226,9 @@ BEM<spacedim>::assemble_system()
       cell_normals_i = fe_values_i.get_all_normal_vectors();
       cell_i->get_dof_indices (local_dof_indices_i);
 
-//  if (cell_i->index()%100==0)
-//    std::cout
-//        << cell_i->index()<< " || ";
-//      << cell_i->vertex(0)<< " || "
-//      << velocity*cell_normals_i[0] << "; "
-//      << local_dof_indices_i[0]
-//      << std::endl;
-
       cell_DLP_matrix = 0.;
       cell_mass_matrix = 0.;
 
-
-
-
       // assembling of the right hand side
       Point<spacedim-1>
       a_unit(0.),
@@ -321,35 +244,27 @@ BEM<spacedim>::assemble_system()
       cell_rhs = 0.;
 
 
-      // In order to obtain the Gram
-      // determinant it is used JxW and
-      // then it is divided by the weight.
-      // Both the jacobian and the normals
-      // are constant on the cell, so they
-      // are taken in the first quadrature
+      // In order to obtain the Gram determinant it is used JxW and then it
+      // is divided by the weight. Both the jacobian and the normals are
+      // constant on the cell, so they are taken in the first quadrature
       // point.
       double constant_factor =
         - pow(fe_values_i.JxW (0)/fe_values_i.get_quadrature().weight(0), 2)
         / numbers::PI
         * velocity*cell_normals_i[0];
 
-      // These are constant on the cell so
-      // there is no need to loop on
-      // quadrature points. Besides, the
-      // cell_rhs vector element index is
-      // chosen = 0 because it is actually
-      // a 1x1 matrix. There should be a
-      // loop over cell dofs, but this is
-      // necessary just for higher degree
+      // These are constant on the cell so there is no need to loop on
+      // quadrature points. Besides, the cell_rhs vector element index is
+      // chosen = 0 because it is actually a 1x1 matrix. There should be a
+      // loop over cell dofs, but this is necessary just for higher degree
       // elements.
       cell_rhs(0) +=
         constant_factor
         *(
           log( A.distance(B) ) + 8.*qlog.weight(0)
         );
-      // A Gauss integration is performed to compute the
-      // integrals on K_i X K_j in the case when i!=j.
-      // cycle on j index
+      // A Gauss integration is performed to compute the integrals on K_i X
+      // K_j in the case when i!=j. cycle on j index
       for (cell_j=dof_handler.begin_active(); cell_j!=endc; ++cell_j)
         {
           fe_values_j.reinit (cell_j);
@@ -358,35 +273,29 @@ BEM<spacedim>::assemble_system()
 
           if (cell_j != cell_i)
             {
-              // The mass matrix has only diagonal
-              // elements.
+              // The mass matrix has only diagonal elements.
               mass_matrix(cell_i->index(), cell_j->index())= 0.;
 
-              // with constant elements there is
-              // only 1 dof per cell, so there is
-              // no real cycle over cell dofs
+              // with constant elements there is only 1 dof per cell, so
+              // there is no real cycle over cell dofs
               for (unsigned int a=0; a<dofs_per_cell; ++a)
                 for (unsigned int q_point_i=0; q_point_i<n_q_points; q_point_i++)
                   for (unsigned int q_point_j=0; q_point_j<n_q_points; q_point_j++)
                     {
-                      // If the integration is performed
-                      // on two different elements there
-                      // are no singularities in the
-                      // domain of integration, so usual
-                      // gauss formulas can be used.
+                      // If the integration is performed on two different
+                      // elements there are no singularities in the domain
+                      // of integration, so usual gauss formulas can be
+                      // used.
 
                       for (unsigned int b=0; b<dofs_per_cell; ++b)
                         {
                           // assembling of double layer potential
 
                           cell_DLP_matrix(a,b) +=
-                            1./numbers::PI
-                            * contract(
-                              cell_normals_i[q_point_i],
-                              (fe_values_j.quadrature_point(q_point_j)
-                               -
-                               fe_values_i.quadrature_point(q_point_i))
-                            )
+                            1./numbers::PI * cell_normals_i[q_point_i] *
+                            (fe_values_j.quadrature_point(q_point_j)
+                             -
+                             fe_values_i.quadrature_point(q_point_i))
                             / pow(
                               (fe_values_j.quadrature_point(q_point_j)
                                -
@@ -396,22 +305,6 @@ BEM<spacedim>::assemble_system()
                             )
                             * fe_values_i.JxW(q_point_i)
                             * fe_values_j.JxW(q_point_j);
-
-//            if ( (cell_i->index()==0) && (cell_j->index()==5) )
-//          std::cout
-//              << "("<< cell_i->index()<<","<<cell_j->index()<<"; "
-//              <<q_point_i<<","<<q_point_j<<") "
-//              << cell_DLP_matrix(a,b)<< " | "
-//              << fe_values_i.quadrature_point(q_point_i)<< "||"
-//              << fe_values_j.quadrature_point(q_point_j)<< "|<"
-//              << contract(
-//            cell_normals_i[q_point_i],
-//            (fe_values_j.quadrature_point(q_point_j)
-//             -
-//             fe_values_i.quadrature_point(q_point_i))
-//                   )<< "<<"
-//              ;
-
                         }
 
 
@@ -428,9 +321,6 @@ BEM<spacedim>::assemble_system()
                         )
                         * fe_values_i.JxW(q_point_i)
                         * fe_values_j.JxW(q_point_j);
-//    std::cout
-//        << cell_j->index()<< " -- "
-//        << cell_rhs(i);
                     }
               for (unsigned int a=0; a<dofs_per_cell; ++a)
                 for (unsigned int b=0; b<dofs_per_cell; ++b)
@@ -440,93 +330,35 @@ BEM<spacedim>::assemble_system()
             }
           else      // case when cell_i=cell_j
             {
-              // The mass matrix is simply a
-              // diagonal matrix with the area of
-              // each element as entries.
+              // The mass matrix is simply a diagonal matrix with the area
+              // of each element as entries.
               for (unsigned q_point_i=0; q_point_i<n_q_points; ++q_point_i)
                 mass_matrix(cell_i->index(), cell_i->index())
                 +=
                   fe_values_i.JxW(q_point_i);
 
-              // The double layer potential matrix
-              // has no diagonal terms since the
-              // scalar product between the cell
-              // normal and the difference
-              // <tt>x-y<\tt> is always zero. This
-              // is because, if <tt>x<\tt> and
-              // <tt>y<\tt> are on the same
-              // element and the element is a
-              // straight segment, their
-              // difference is always orthogonal
-              // to the cell normal.
+              // The double layer potential matrix has no diagonal terms
+              // since the scalar product between the cell normal and the
+              // difference <tt>x-y<\tt> is always zero. This is because,
+              // if <tt>x<\tt> and <tt>y<\tt> are on the same element and
+              // the element is a straight segment, their difference is
+              // always orthogonal to the cell normal.
               DLP_matrix(cell_i->index(), cell_i->index()) = 0.;
             }
 
           cell_DLP_matrix = 0.;
-
-//    std::cout
-//        << cell_DLP_matrix(0,0)<< " - ";
-
         }
 
-//      std::cout
-//    << cell_rhs(0)<<std::endl;
-
       for (unsigned int a=0; a<dofs_per_cell; ++a)
         system_rhs(local_dof_indices_i[a]) += cell_rhs(a);
 
       // end of assembling of the right hand side
-
     }
 
-//     TableIndices<2> table_indices(dof_handler.n_dofs(), dof_handler.n_dofs());
-
   system_matrix.add(
     1, mass_matrix,
     -1, DLP_matrix
   );
-
-//     std::cout
-//  << "rhs"<< std::endl;
-//     for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-//  std::cout
-//      << system_rhs(i)<< " | "
-//      << std::endl;
-
-//     std::cout
-//  << "DLP matrix"<< std::endl;
-//     for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-//     {
-//  for (unsigned int j=0; j<dof_handler.n_dofs(); ++j)
-//      std::cout
-//    << DLP_matrix(i,j)<< " | ";
-//  std::cout
-//      << std::endl;
-//     }
-
-//     std::cout
-//  << "mass matrix"<< std::endl;
-//     for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-//     {
-//  for (unsigned int j=0; j<dof_handler.n_dofs(); ++j)
-//      std::cout
-//    << mass_matrix(i,j)<< " | ";
-//  std::cout
-//      << std::endl;
-//     }
-
-//     std::cout
-//  << "system matrix"<< std::endl;
-//     for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-//     {
-//  for (unsigned int j=0; j<dof_handler.n_dofs(); ++j)
-//      std::cout
-//    << system_matrix(i,j)<< " | ";
-//  std::cout
-//      << std::endl;
-//     }
-
-
 }
 
 template <int spacedim>
@@ -542,8 +374,7 @@ BEM<spacedim>::solve()
             system_rhs,
             PreconditionIdentity() );
 
-  // Smoothen the piecewise constant
-  // solution to a continuous piecewise
+  // Smoothen the piecewise constant solution to a continuous piecewise
   // quadratic function.
   smooth_solution.reinit(dof_handler_q.n_dofs());
 
@@ -552,10 +383,8 @@ BEM<spacedim>::solve()
                         dof_handler_q,
                         smooth_solution);
 
-  // Calculate the tangential derivative
-  // of the smoothened potential, that
-  // is the tangential component of the
-  // perturbation induced in the
+  // Calculate the tangential derivative of the smoothened potential, that
+  // is the tangential component of the perturbation induced in the
   // velocity.
   tangential_velocity.reinit(tria.n_active_cells());
   tangential_derivative.reinit(tria.n_active_cells());
@@ -582,20 +411,6 @@ BEM<spacedim>::solve()
       fe_values_q.reinit(cell);
       cell->get_dof_indices (local_dof_indices);
 
-//QUADRATIC INTERPOLATION OF THE POTENTIAL
-//  std::cout
-//      << "indices - ";
-//  for(unsigned int i=0; i<fe_q.dofs_per_cell; ++i)
-//      std::cout
-//    << local_dof_indices[(i==0? 0
-//              :
-//              i==1 ?  2
-//              :
-//              1
-//             )]<< ", ";
-//  std::cout<<std::endl;
-
-
       cell_normals = fe_values_q.get_all_normal_vectors();
       for (unsigned int i=0; i<q_iterated.size(); ++i)
         {
@@ -603,97 +418,32 @@ BEM<spacedim>::solve()
           cell_tangentials[i][1] = -cell_normals[i][0];
         }
 
-      // Create a vector where gradients at
-      // quadrature points are
-      // stored. Notice that the first
-      // factor (smooth_solution..) is taken
-      // so that it is the coefficient of
-      // the fun_th shape function of the cell.
+      // Create a vector where gradients at quadrature points are stored.
+      // Notice that the first factor (smooth_solution..) is taken so that
+      // it is the coefficient of the fun_th shape function of the cell.
       std::vector< Tensor<1,spacedim> > gradient(q_iterated.size());
       for (unsigned int pnt=0; pnt<q_iterated.size(); ++pnt)
         for (unsigned int fun=0; fun<fe_q.dofs_per_cell; ++fun)
-          {
-            gradient[pnt]+=
-              smooth_solution(local_dof_indices[fun])
-//QUADRATIC INTERPOLATION OF THE POTENTIAL
-//    smooth_solution(local_dof_indices[(fun==0? 0
-//               :
-//               fun==1 ? 2
-//               :
-//               1
-//             )]
-//                   )
-              *
-              fe_values_q.shape_grad(fun, pnt);
-
-//      std::cout
-//    << smooth_solution(local_dof_indices[(fun==0? 0
-//                  :
-//                  fun==1 ? 2
-//                  :
-//                  1
-//                )])<< " | "
-//    << fe_values_q.shape_grad(fun, pnt)
-//    << std::endl;
-          }
-
-//  std::cout
-//      << "gradienti "<<std::endl;
-//  for(unsigned int pnt=0; pnt<q_iterated.size(); ++pnt)
-//      std::cout
-//    << pnt << " - "
-//    << gradient[pnt]
-//    << std::endl;
-
-//    std::cout
-//      << "derivata tangenziale "
-//      << cell->index()
-//      << std::endl;
-
+          gradient[pnt]+=
+            smooth_solution(local_dof_indices[fun])
+            *
+            fe_values_q.shape_grad(fun, pnt);
 
       for (unsigned int pnt=0; pnt<q_iterated.size(); ++pnt)
         {
-          tangential_derivative(cell->index())=
-            contract(
-              gradient[pnt],
-              cell_tangentials[pnt]
-            )
-            +
-            contract(velocity,
-                     cell_tangentials[0]
-                    );
-          error(cell->index())=
-            tangential_derivative(cell->index())
-            -
-            2. * velocity[0]
-            * sin (numbers::PI/2. - (1./2.+cell->index()) * 2.*numbers::PI/tria.n_active_cells() );
-          tangential_velocity(cell->index())=
-            contract(velocity,
-                     cell_tangentials[0]
-                    );
+          tangential_derivative(cell->index()) =
+            gradient[pnt] * cell_tangentials[pnt] +
+            velocity * cell_tangentials[0];
 
-        }
+          error(cell->index()) =
+            tangential_derivative(cell->index()) -
+            2. * velocity[0] * sin(numbers::PI / 2. -
+                                   (1. / 2. + cell->index()) * 2. *
+                                   numbers::PI / tria.n_active_cells());
 
-// //   for(unsigned int pnt=0; pnt<q_iterated.size(); ++pnt)
-//      std::cout
-// //     << pnt <<" - "
-//    << tangential_derivative(cell->index())
-//    << std::endl;
+          tangential_velocity(cell->index()) = velocity * cell_tangentials[0];
+        }
     }
-
-
-//     DataOut<spacedim-1, DoFHandler<spacedim-1,spacedim> > dataout;
-//     dataout.attach_dof_handler(dof_handler_q);
-//     dataout.add_data_vector(smooth_solution, "quadratic_potential");
-//     dataout.build_patches(fe_q.degree);
-//     char outname[50];
-//     sprintf(outname, "bem_gradient.vtk");
-//     std::ofstream file(outname);
-//     dataout.write_vtk(file);
-//     dataout.write_vtk(logfile);
-
-
-
 }
 
 template <int spacedim>
@@ -707,7 +457,6 @@ BEM<spacedim>::output_results()
   dataout.add_data_vector(tangential_derivative, "tangential_velocity", DataOut<spacedim-1, DoFHandler<spacedim-1,spacedim> >::type_cell_data);
   dataout.add_data_vector(error, "error", DataOut<spacedim-1, DoFHandler<spacedim-1,spacedim> >::type_cell_data);
   dataout.build_patches();
-//    dataout.build_patches(fe_q.degree);
   dataout.write_vtk(logfile);
 
 }

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.