From fd551996b10495da8405343d42d9b7aa128072a0 Mon Sep 17 00:00:00 2001 From: heltai Date: Tue, 3 Mar 2009 18:07:13 +0000 Subject: [PATCH] Ameliorated parameter selection, and assembly of the system matrix. Now it is slightly more readable. git-svn-id: https://svn.dealii.org/trunk@18444 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-34/parameters.prm | 12 +- deal.II/examples/step-34/step-34.cc | 235 +++++++++++++++++------- 2 files changed, 177 insertions(+), 70 deletions(-) diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index b4995437f6..d91513a956 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -1,8 +1,14 @@ # Listing of Parameters # --------------------- - set Number of cycles = 3 + +subsection Inner quadrature rule + set Quadrature order = 2 + set Quadrature type = gauss +end + + subsection Outer quadrature rule set Quadrature order = 0 set Quadrature type = midpoint @@ -11,11 +17,11 @@ end subsection Wind function # Any constant used inside the function which is not a variable name. - set Function constants = angle=0 + set Function constants = # Separate vector valued expressions by ';' as ',' is used internally by the # function parser. - set Function expression = cos(angle); sin(angle); 0 + set Function expression = 1; 0; 0 # The name of the variables as they will be used in the function, separated # by ','. diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index c39ccb761d..eb87d30313 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -82,8 +82,16 @@ public: // \f] void compute_SD_integral_on_cell(vector > > &dst, typename DoFHandler::active_cell_iterator &cell, - const vector > &q); + const vector > &q, + const Function &rhs); + // The following two functions are the actual calculations of the + // single and double layer potential kernels, with a minus sign in + // front of them. They are well defined only if the vector $r = + // x-y$ is different from zero. + double nS(const Point &R); + Point nD(const Point &R); + private: // The following two helper functions should only be called when // dim=3. If this is not the case, the default implementation is @@ -170,7 +178,8 @@ private: // quadrature. We use a parsed function, for its ease of // definition, and the quadrature formula Functions::ParsedFunction wind; - SmartPointer > quadrature_pointer; + SmartPointer > outer_quadrature_pointer; + SmartPointer > inner_quadrature_pointer; unsigned int n_cycles; }; @@ -204,6 +213,13 @@ void BEMProblem::read_parameters(std::string filename) { Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names())); prm.declare_entry("Quadrature order", "0", Patterns::Integer()); prm.leave_subsection(); + + + prm.enter_subsection("Inner quadrature rule"); + prm.declare_entry("Quadrature type", "midpoint", + Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names())); + prm.declare_entry("Quadrature order", "0", Patterns::Integer()); + prm.leave_subsection(); prm.enter_subsection("Wind function"); Functions::ParsedFunction::declare_parameters(prm, dim); @@ -214,7 +230,14 @@ void BEMProblem::read_parameters(std::string filename) { n_cycles = prm.get_integer("Number of cycles"); prm.enter_subsection("Outer quadrature rule"); - static QuadratureSelector quadrature + static QuadratureSelector outer_quadrature + (prm.get("Quadrature type"), + prm.get_integer("Quadrature order")); + prm.leave_subsection(); + + + prm.enter_subsection("Inner quadrature rule"); + static QuadratureSelector inner_quadrature (prm.get("Quadrature type"), prm.get_integer("Quadrature order")); prm.leave_subsection(); @@ -223,9 +246,41 @@ void BEMProblem::read_parameters(std::string filename) { wind.parse_parameters(prm); prm.leave_subsection(); - quadrature_pointer = &quadrature; + outer_quadrature_pointer = &outer_quadrature; + inner_quadrature_pointer = &inner_quadrature; } + + +template +double LaplaceKernelIntegration::nS(const Point &R) { + if(dim == 2) + return (-std::log(R.norm()) / numbers::PI); + else if(dim == 3) + return (1./(R.norm()*numbers::PI) ); + else { + Assert(false, ExcInternalError()); + } + return 0.; +} + + + +template +Point LaplaceKernelIntegration::nD(const Point &R) { + Point D(R); + if(dim == 2) + D /= -numbers::PI * R.square(); + else if(dim == 3) + D /= -2*numbers::PI * R.square() * R.norm(); + else { + Assert(false, ExcInternalError()); + } + return D; +} + + + template <> LaplaceKernelIntegration<3>::LaplaceKernelIntegration(FiniteElement<2,3> &fe) { @@ -247,6 +302,7 @@ LaplaceKernelIntegration<3>::LaplaceKernelIntegration(FiniteElement<2,3> &fe) fe_values = new FEValues<2,3>(fe,quadrature, update_values | update_jacobians | + update_cell_normal_vectors | update_quadrature_points ); } @@ -312,12 +368,25 @@ template <> void LaplaceKernelIntegration<3>::compute_SD_integral_on_cell(vector > > &dstvv, DoFHandler<2,3>::active_cell_iterator &cell, - const vector > &q_points) + const vector > &q_points, + const Function<3> &rhs) { fe_values->reinit(cell); - vector > jacobians = fe_values->get_jacobians(); - vector > quad_points = fe_values->get_quadrature_points(); - + const vector > &jacobians = fe_values->get_jacobians(); + const vector > &quad_points = fe_values->get_quadrature_points(); + const vector > &normals = fe_values->get_cell_normal_vectors(); + + static vector > cell_wind + ( (*fe_values).n_quadrature_points, Vector(3) ); + static vector normal_wind(quad_points.size()); + + rhs.vector_value_list(quad_points, cell_wind); + + for(unsigned int q=0; q r,a1,a2,n,r_c,n_c; Assert(dstvv.size() == fe_values->dofs_per_cell, @@ -348,7 +417,7 @@ LaplaceKernelIntegration<3>::compute_SD_integral_on_cell(vectorshape_value(i,inner_q_point); + i_S[inner_q_point]= term_S(r,a1,a2,n,rn_c) * normal_wind[inner_q_point]; i_D[inner_q_point]= term_D(r,a1,a2) * fe_values->shape_value(i,inner_q_point); } dst[0] = (i_S[3]-i_S[1]-i_S[2]+i_S[0]); @@ -410,39 +479,34 @@ void BEMProblem::assemble_system() { typename DoFHandler::active_cell_iterator celli = dh.begin_active(), cellj = dh.begin_active(), - cellvi = dhv.begin_active(), endc = dh.end(); // Outer quadrature rule. If we choose midpoint quadrature rule, // then this is a collocation method. If we choose any other // Quadrature rule, then this is Galerkin method. - Quadrature &quadrature_outer = *quadrature_pointer; - QMidpoint quadrature_mid; - + Quadrature &outer_quadrature = *outer_quadrature_pointer; + Quadrature &inner_quadrature = *inner_quadrature_pointer; - FEValues fe_outer(fe, quadrature_outer, + FEValues fe_outer(fe, outer_quadrature, update_values | update_cell_normal_vectors | - update_quadrature_points); + update_quadrature_points | + update_JxW_values); - FEValues fe_inner(fe, quadrature_mid, + FEValues fe_inner(fe, inner_quadrature, update_values | update_cell_normal_vectors | - update_quadrature_points); + update_quadrature_points | + update_JxW_values); const unsigned int n_q_points_outer = fe_outer.n_quadrature_points; const unsigned int n_q_points_inner = fe_inner.n_quadrature_points; vector dofs_i(fe.dofs_per_cell); vector dofs_j(fe.dofs_per_cell); - vector dofs_v_i(fev.dofs_per_cell); - - vector > > single_double_layer_potentials - (fe.dofs_per_cell, vector > - (n_q_points_outer, vector (2, 0.) ) ); - vector > cell_wind(n_q_points_inner, Vector(dim) ); - vector normal_wind(n_q_points_inner); + vector > inner_cell_wind(n_q_points_inner, Vector(dim) ); + double inner_normal_wind; Vector local_rhs(fe.dofs_per_cell); FullMatrix local_matrix(fe.dofs_per_cell, fe.dofs_per_cell); @@ -450,28 +514,32 @@ void BEMProblem::assemble_system() { // The kernel. LaplaceKernelIntegration kernel(fe); - // i runs on outer integration, j runs on inner integration. - for(; celli != endc; ++celli, ++cellvi) { + vector > > single_double_layer_potentials + (fe.dofs_per_cell, vector > + (n_q_points_outer, vector (2, 0.) ) ); + + Point R; + + + // The index i runs on outer integration, while j runs on inner integration. + for(; celli != endc; ++celli) { fe_outer.reinit(celli); const vector > &q_points_outer = fe_outer.get_quadrature_points(); const vector > &normals_outer = fe_outer.get_cell_normal_vectors(); celli->get_dof_indices(dofs_i); - cellvi->get_dof_indices(dofs_v_i); - // Now the mass matrix and the single and double layer - // potentials. Notice that instead of integrating the single - // layer potential against the normal velocity, we integrate - // it against the average value of the velocity in the given - // cell. - // - // The reason for proceeding in this way is that in the - // current three-dimensional formulation we can only integrate - // the constants against the single layer potential. While - // this is certainly a rough approximation, it suffices for - // the purpose of this example. for(cellj = dh.begin_active(); cellj != endc; ++cellj) { + + // If we are on the same cell, then the integrals we are + // performing are singular, and they require a special + // treatment, as explained in the introduction. + // + // In all other cases, standard Gauss quadrature rules can + // be used. + bool is_singular = (cellj->index() == celli->index()); + local_rhs = 0; local_matrix = 0; @@ -480,40 +548,73 @@ void BEMProblem::assemble_system() { const vector > &q_points_inner = fe_inner.get_quadrature_points(); const vector > &normals_inner = fe_inner.get_cell_normal_vectors(); - wind.vector_value_list(q_points_inner, cell_wind); - - - for(unsigned int q=0; qindex() != cellj->index()) - local_matrix(i,j) += ( -single_double_layer_potentials[j][qo][1] * - fe_outer.shape_value(i,qo) * - fe_outer.JxW(qo) ); - // When the indices are the same, we assemble - // also the mass matrix. - if(celli->index() == cellj->index()) - local_matrix(i,j) += ( fe_outer.shape_value(i,qo) * - fe_outer.shape_value(j,qo) * - fe_outer.JxW(qo) ); + for(unsigned int j=0; j