From 9bf53cb8bb453987419125a1dcc90c7c5e96f062 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Tue, 15 Sep 2015 15:33:58 -0500 Subject: [PATCH] Testsuite: Clean up codim_one/bem; remove unused code fragments --- tests/codim_one/bem.cc | 355 ++++++----------------------------------- 1 file changed, 52 insertions(+), 303 deletions(-) diff --git a/tests/codim_one/bem.cc b/tests/codim_one/bem.cc index 7b28ad6d92..ea3dc57073 100644 --- a/tests/codim_one/bem.cc +++ b/tests/codim_one/bem.cc @@ -79,10 +79,8 @@ private: FE_DGQ fe; DoFHandler 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 fe_q; DoFHandler dof_handler_q; @@ -153,64 +151,9 @@ BEM::run() tria.set_boundary(1); - -// std::cout -// << "solution" -// << std::endl; -// for(unsigned int i=0; i::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 DLP_matrix(dof_handler.n_dofs()); // DLP stands for double layer potential FullMatrix mass_matrix(dof_handler.n_dofs()); @@ -271,14 +213,6 @@ BEM::assemble_system() std::vector local_dof_indices_i (dofs_per_cell); std::vector 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())<::active_cell_iterator cell_i = dof_handler.begin_active(), cell_j = dof_handler.begin_active(), @@ -292,20 +226,9 @@ BEM::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 a_unit(0.), @@ -321,35 +244,27 @@ BEM::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::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::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()<<","<index()<<"; " -// <::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::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_iindex(), 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 - // x-y<\tt> is always zero. This - // is because, if x<\tt> and - // 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 x-y<\tt> is always zero. This is because, + // if x<\tt> and 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)< 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 @@ -542,8 +374,7 @@ BEM::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::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::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::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; pntindex() -// << std::endl; - + gradient[pnt]+= + smooth_solution(local_dof_indices[fun]) + * + fe_values_q.shape_grad(fun, pnt); for (unsigned int pnt=0; pntindex())= - 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; pntindex()) -// << std::endl; + tangential_velocity(cell->index()) = velocity * cell_tangentials[0]; + } } - - -// DataOut > 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 @@ -707,7 +457,6 @@ BEM::output_results() dataout.add_data_vector(tangential_derivative, "tangential_velocity", DataOut >::type_cell_data); dataout.add_data_vector(error, "error", DataOut >::type_cell_data); dataout.build_patches(); -// dataout.build_patches(fe_q.degree); dataout.write_vtk(logfile); } -- 2.39.5