From: bangerth Date: Tue, 14 Apr 2009 21:50:52 +0000 (+0000) Subject: Reindent. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1b8ac5d54848773ecb7d2d53793a89459808659e;p=dealii-svn.git Reindent. git-svn-id: https://svn.dealii.org/trunk@18612 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index 07509c97c9..1fcc39fac7 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -66,122 +66,175 @@ using namespace dealii; // @sect3{Single and double layer operator kernels} - // First, let us define a bit of the boundary integral equation machinery. - - // The following two functions are the actual calculations of the - // single and double layer potential kernels, that is $G$ and $\nabla - // G$. They are well defined only if the vector $R = \mathbf{y}-\mathbf{x}$ is - // different from zero. - // - // Whenever the integration is performed with the singularity - // inside the given cell, then a special quadrature formula is - // used that allows one to integrate arbitrary functions against a - // singular weight on the reference cell. - // There are two options when the integral is singular. One could - // take into account the singularity inside the quadrature formula - // as a weigthing function, or one could use a quadrature formula - // that is taylored to integrate singular objects, but where the - // actual weighting function is one. The use of the first method - // requires the user to provide "desingularized" single and - // double layer potentials which can then be integrated on the - // given cell. When the @p factor_out_singularity parameter is set - // to true, then the computed kernels do not conatain the singular - // factor, which is included in the quadrature formulas as a - // weighting function. This works best in two dimension, where the - // singular integrals are integrals along a segment of a - // logarithmic singularity. - // + // First, let us define a bit of the + // boundary integral equation + // machinery. + + // The following two functions are + // the actual calculations of the + // single and double layer potential + // kernels, that is $G$ and $\nabla + // G$. They are well defined only if + // the vector $R = + // \mathbf{y}-\mathbf{x}$ is + // different from zero. + // + // Whenever the integration is + // performed with the singularity + // inside the given cell, then a + // special quadrature formula is used + // that allows one to integrate + // arbitrary functions against a + // singular weight on the reference + // cell. There are two options when + // the integral is singular. One + // could take into account the + // singularity inside the quadrature + // formula as a weigthing function, + // or one could use a quadrature + // formula that is taylored to + // integrate singular objects, but + // where the actual weighting + // function is one. The use of the + // first method requires the user to + // provide "desingularized" single + // and double layer potentials which + // can then be integrated on the + // given cell. When the @p + // factor_out_singularity parameter + // is set to true, then the computed + // kernels do not conatain the + // singular factor, which is included + // in the quadrature formulas as a + // weighting function. This works + // best in two dimension, where the + // singular integrals are integrals + // along a segment of a logarithmic + // singularity. + // //TODO: Can you elaborate in formulas? - // These integrals are somewhat delicate, because inserting a - // factor Jx in the variable of integration does not result only - // in a factor J appearing as a constant factor on the entire - // integral, but also on an additional integral to be added, that - // contains the logarithm of J. For this reason in two dimensions - // we opt for the desingularized kernel, and use the QGaussLogR - // quadrature formula, that takes care of integrating the correct - // weight for us. - // - // In the three dimensional case the singular integral is taken - // care of using the QGaussOneOverR quadrature formula. We could - // use the desingularized kernel here as well, but this would - // require us to be careful about the different scaling of $r$ in - // the reference cell and in real space. The quadrature formula - // uses as weight $1/r$ in local coordinates, while we need to - // integrate $1/R$ in real coordinates. A factor of $r/R$ has to be - // introduced in the quadrature formula. This can be done - // manually, or we simply calculate the standard kernels and then - // use a desingularized quadrature formula, i.e., one which is - // taylored for singular integrals, but whose weight is 1 instead - // of the singularity. - // - // Notice that the QGaussLog quadrature formula is made to integrate - // $f(x)\ln |\mathbf{x}-\mathbf{x}_0|$, but the kernel for two dimensional - // problems has the opposite sign. This is taken care of by switching the - // sign of the two dimensional desingularized kernel. - // - // The last argument to both functions is simply ignored in three - // dimensions. + // These integrals are somewhat + // delicate, because inserting a + // factor Jx in the variable of + // integration does not result only + // in a factor J appearing as a + // constant factor on the entire + // integral, but also on an + // additional integral to be added, + // that contains the logarithm of + // J. For this reason in two + // dimensions we opt for the + // desingularized kernel, and use the + // QGaussLogR quadrature formula, + // that takes care of integrating the + // correct weight for us. + // + // In the three dimensional case the + // singular integral is taken care of + // using the QGaussOneOverR + // quadrature formula. We could use + // the desingularized kernel here as + // well, but this would require us to + // be careful about the different + // scaling of $r$ in the reference + // cell and in real space. The + // quadrature formula uses as weight + // $1/r$ in local coordinates, while + // we need to integrate $1/R$ in real + // coordinates. A factor of $r/R$ has + // to be introduced in the quadrature + // formula. This can be done + // manually, or we simply calculate + // the standard kernels and then use + // a desingularized quadrature + // formula, i.e., one which is + // taylored for singular integrals, + // but whose weight is 1 instead of + // the singularity. + // + // Notice that the QGaussLog + // quadrature formula is made to + // integrate $f(x)\ln + // |\mathbf{x}-\mathbf{x}_0|$, but + // the kernel for two dimensional + // problems has the opposite + // sign. This is taken care of by + // switching the sign of the two + // dimensional desingularized kernel. + // + // The last argument to both + // functions is simply ignored in + // three dimensions. namespace LaplaceKernel { -template -double single_layer(const Point &R, - const bool factor_out_2d_singularity = false) { + template + double single_layer(const Point &R, + const bool factor_out_2d_singularity = false) { switch(dim) { - case 2: - if (factor_out_2d_singularity == true) - return -1./(2*numbers::PI); - else - return (-std::log(R.norm()) / (2*numbers::PI) ); - - case 3: - return (1./( R.norm()*4*numbers::PI ) ); - - default: - Assert(false, ExcInternalError()); - return 0.; + case 2: + if (factor_out_2d_singularity == true) + return -1./(2*numbers::PI); + else + return (-std::log(R.norm()) / (2*numbers::PI) ); + + case 3: + return (1./( R.norm()*4*numbers::PI ) ); + + default: + Assert(false, ExcInternalError()); + return 0.; } -} + } -template -Point double_layer(const Point &R, - const bool factor_out_2d_singularity = false) { - switch(dim) { - case 2: - if (factor_out_2d_singularity) - return Point(); - else - return R / (-2*numbers::PI * R.square()); - case 3: - return R / ( -4*numbers::PI * R.square()*R.norm() ); + template + Point double_layer(const Point &R, + const bool factor_out_2d_singularity = false) { + switch(dim) { + case 2: + if (factor_out_2d_singularity) + return Point(); + else + return R / (-2*numbers::PI * R.square()); + case 3: + return R / ( -4*numbers::PI * R.square()*R.norm() ); default: - Assert(false, ExcInternalError()); - return Point(); + Assert(false, ExcInternalError()); + return Point(); } -} + } } // @sect3{The BEMProblem class} - // The structure of a boundary element method code is very similar to the - // structure of a finite element code, and so the member functions of this - // class are like those of most of the other tutorial programs. In - // particular, by now you should be familiar with reading parameters from - // an external file, and with the splitting of the different tasks into - // different modules. The same applies to boundary element methods, and we - // won't comment too much on them, except on the differences. + // The structure of a boundary + // element method code is very + // similar to the structure of a + // finite element code, and so the + // member functions of this class are + // like those of most of the other + // tutorial programs. In particular, + // by now you should be familiar with + // reading parameters from an + // external file, and with the + // splitting of the different tasks + // into different modules. The same + // applies to boundary element + // methods, and we won't comment too + // much on them, except on the + // differences. template class BEMProblem { -public: + public: BEMProblem(); void run(); -private: + private: void read_parameters (const std::string &filename); @@ -189,154 +242,274 @@ private: void refine_and_resize(); - // The only really different function that we find here is the - // assembly routine. We wrote this function in the most possible - // general way, in order to allow for easy generalization to - // higher order methods and to different fundamental solutions - // (e.g., Stokes or Maxwell). - // - // The most noticeable difference is the fact that the final - // matrix is full, and that we have a nested loop inside the usual - // loop on cells that visits all support points of the degrees of - // freedom. Moreover, when the support point lies inside the cell - // which we are visiting, then the integral we perform becomes - // singular. - // - // The practical consequence is that we have two sets of - // quadrature formulas, finite element values and temporary - // storage, one for standard integration and one for the singular - // integration, which are used where necessary. + // The only really different + // function that we find here is + // the assembly routine. We wrote + // this function in the most + // possible general way, in order + // to allow for easy + // generalization to higher order + // methods and to different + // fundamental solutions (e.g., + // Stokes or Maxwell). + // + // The most noticeable difference + // is the fact that the final + // matrix is full, and that we + // have a nested loop inside the + // usual loop on cells that + // visits all support points of + // the degrees of freedom. + // Moreover, when the support + // point lies inside the cell + // which we are visiting, then + // the integral we perform + // becomes singular. + // + // The practical consequence is + // that we have two sets of + // quadrature formulas, finite + // element values and temporary + // storage, one for standard + // integration and one for the + // singular integration, which + // are used where necessary. void assemble_system(); - // Notwithstanding the fact that the matrix is full, we use a SparseMatrix - // object and the SparseDirectUMFPACK solver, since in our experience it - // works better than using, for example, the LapackFullMatrix class. Of - // course, using a SparseMatrix object to store the matrix is wasteful, - // but at least for the moment that is all the SparseDirectUMFPACK class - // can deal with. - // - // An alternative approach would be the - // use of the GMRES method; however the construction of an efficient - // preconditioner for boundary element methods is not a trivial - // issue, and we won't treat this problem here. - // - // Moreover, we should notice that the solution we will obtain will only - // be unique up to an additive constant. This is taken care of in the - // solve_system() method, which filters out the mean value of - // the solution at the end of the computation. + // Notwithstanding the fact that + // the matrix is full, we use a + // SparseMatrix object and the + // SparseDirectUMFPACK solver, + // since in our experience it + // works better than using, for + // example, the LapackFullMatrix + // class. Of course, using a + // SparseMatrix object to store + // the matrix is wasteful, but at + // least for the moment that is + // all the SparseDirectUMFPACK + // class can deal with. + // + // An alternative approach would + // be the use of the GMRES + // method; however the + // construction of an efficient + // preconditioner for boundary + // element methods is not a + // trivial issue, and we won't + // treat this problem here. + // + // Moreover, we should notice + // that the solution we will + // obtain will only be unique up + // to an additive constant. This + // is taken care of in the + // solve_system() + // method, which filters out the + // mean value of the solution at + // the end of the computation. void solve_system(); - // Once we obtained the solution, we compute the $L^2$ error of - // the computed potential as well as the $L^\infty$ error of the - // approximation of the solid angle. The mesh we are using is an - // approximation of a smooth curve, therefore the computed - // diagonal matrix of fraction of angles or solid angles - // $\alpha(\mathbf{x})$ should be constantly equal to $\frac - // 12$. In this routine we output the error on the potential and - // the error in the approximation of the computed angle. Notice - // that the latter error is actually not the error in the - // computation of the angle, but a measure of how well we are - // approximating the sphere and the circle. - // - // Experimenting a little with the computation of the angles gives very - // accurate results for simpler geometries. To verify this you can comment - // out, in the read_domain() method, the tria.set_boundary(1, boundary) - // line, and check the alpha that is generated by the program. By removing - // this call, whenever the mesh is refined new nodes will be placed along - // the straight lines that made up the coarse mesh, rather than be pulled - // onto the surface that we really want to approximate. In the three - // dimensional case, the coarse grid of the sphere is obtained starting - // from a cube, and the obtained values of alphas are exactly $\frac 12$ - // on the nodes of the faces, $\frac 14$ on the nodes of the edges and - // $\frac 18$ on the 8 nodes of the vertices. + // Once we obtained the solution, + // we compute the $L^2$ error of + // the computed potential as well + // as the $L^\infty$ error of the + // approximation of the solid + // angle. The mesh we are using + // is an approximation of a + // smooth curve, therefore the + // computed diagonal matrix of + // fraction of angles or solid + // angles $\alpha(\mathbf{x})$ + // should be constantly equal to + // $\frac 12$. In this routine we + // output the error on the + // potential and the error in the + // approximation of the computed + // angle. Notice that the latter + // error is actually not the + // error in the computation of + // the angle, but a measure of + // how well we are approximating + // the sphere and the circle. + // + // Experimenting a little with + // the computation of the angles + // gives very accurate results + // for simpler geometries. To + // verify this you can comment + // out, in the read_domain() + // method, the + // tria.set_boundary(1, boundary) + // line, and check the alpha that + // is generated by the + // program. By removing this + // call, whenever the mesh is + // refined new nodes will be + // placed along the straight + // lines that made up the coarse + // mesh, rather than be pulled + // onto the surface that we + // really want to approximate. In + // the three dimensional case, + // the coarse grid of the sphere + // is obtained starting from a + // cube, and the obtained values + // of alphas are exactly $\frac + // 12$ on the nodes of the faces, + // $\frac 14$ on the nodes of the + // edges and $\frac 18$ on the 8 + // nodes of the vertices. void compute_errors(const unsigned int cycle); - // Once we obtained a solution on the codimension one domain, we - // want to interpolate it to the rest of the - // space. This is done by performing again the convolution of the - // solution with the kernel in the compute_exterior_solution() function. - // - // We would like to plot the velocity variable which is the - // gradient of the potential solution. The potential solution is - // only known on the boundary, but we use the convolution with the - // fundamental solution to interpolate it on a standard dim - // dimensional continuous finite element space. The plot of the - // gradient of the extrapolated solution will give us the velocity - // we want. - // - // In addition to the solution on the exterior domain, we also output the - // solution on the domain's boundary in the output_results() function, of - // course. + // Once we obtained a solution on + // the codimension one domain, we + // want to interpolate it to the + // rest of the space. This is + // done by performing again the + // convolution of the solution + // with the kernel in the + // compute_exterior_solution() + // function. + // + // We would like to plot the + // velocity variable which is the + // gradient of the potential + // solution. The potential + // solution is only known on the + // boundary, but we use the + // convolution with the + // fundamental solution to + // interpolate it on a standard + // dim dimensional continuous + // finite element space. The plot + // of the gradient of the + // extrapolated solution will + // give us the velocity we want. + // + // In addition to the solution on + // the exterior domain, we also + // output the solution on the + // domain's boundary in the + // output_results() function, of + // course. void compute_exterior_solution(); void output_results(const unsigned int cycle); - // The usual deal.II classes can be used for boundary element - // methods by specifying the "codimension" of the problem. This is - // done by setting the optional second template arguments to - // Triangulation, FiniteElement and DoFHandler to the dimension of - // the embedding space. In our case we generate either 1 or 2 - // dimensional meshes embedded in 2 or 3 dimensional spaces. - // - // The optional argument by default is equal to the first - // argument, and produces the usual finite element classes that we - // saw in all previous examples. + // The usual deal.II classes can + // be used for boundary element + // methods by specifying the + // "codimension" of the + // problem. This is done by + // setting the optional second + // template arguments to + // Triangulation, FiniteElement + // and DoFHandler to the + // dimension of the embedding + // space. In our case we generate + // either 1 or 2 dimensional + // meshes embedded in 2 or 3 + // dimensional spaces. + // + // The optional argument by + // default is equal to the first + // argument, and produces the + // usual finite element classes + // that we saw in all previous + // examples. Triangulation tria; FE_Q fe; DoFHandler dh; - // In BEM methods, the matrix that is generated is - // dense. Depending on the size of the problem, the final system - // might be solved by direct LU decomposition, or by iterative - // methods. In this example we use the SparseDirectUMFPACK solver, - // applied to a "fake" sparse matrix (a sparse matrix will all - // entries different from zero). We found that this method is - // faster than using a LapackFullMatrix object. + // In BEM methods, the matrix + // that is generated is + // dense. Depending on the size + // of the problem, the final + // system might be solved by + // direct LU decomposition, or by + // iterative methods. In this + // example we use the + // SparseDirectUMFPACK solver, + // applied to a "fake" sparse + // matrix (a sparse matrix will + // all entries different from + // zero). We found that this + // method is faster than using a + // LapackFullMatrix object. SparsityPattern sparsity; SparseMatrix system_matrix; Vector system_rhs; - // The next two variables will denote the - // solution $\phi$ as well as a vector - // that will hold the values of - // $\alpha(\mathbf x)$ (the fraction of - // space visible from a point $\mathbf - // x$) at the support points of our shape + // The next two variables will + // denote the solution $\phi$ as + // well as a vector that will + // hold the values of + // $\alpha(\mathbf x)$ (the + // fraction of space visible from + // a point $\mathbf x$) at the + // support points of our shape // functions. Vector phi; Vector alpha; - // The convergence table is used to output errors in the exact - // solution and in the computed alphas. + // The convergence table is used + // to output errors in the exact + // solution and in the computed + // alphas. ConvergenceTable convergence_table; - // The following variables are the ones that we fill through a - // parameter file. The new objects that we use in this example - // are the Functions::ParsedFunction object and the QuadratureSelector - // object. - // - // The Functions::ParsedFunction class allows us to easily and quickly - // define new function objects via parameter files, with custom - // definitions which can be very complex (see the documentation of that - // class for all the available options). - // - // We will allocate the quadrature object using the QuadratureSelector - // class that allows us to generate quadrature formulas based on an - // identifying string and on the possible degree of the formula itself. We - // used this to allow custom selection of the quadrature formulas for the - // standard integration, and to define the order of the singular - // quadrature rule. - // - // Notice that the pointer given below for the quadrature rule is - // only used for non singular integrals. Whenever the integral is - // singular, then only the degree of the quadrature pointer is - // used, and the integration is a special one (see the - // assemble_matrix() function below for further details). - // - // We also define a couple of parameters which are used in case we - // wanted to extend the solution to the entire domain. + // The following variables are + // the ones that we fill through + // a parameter file. The new + // objects that we use in this + // example are the + // Functions::ParsedFunction + // object and the + // QuadratureSelector object. + // + // The Functions::ParsedFunction + // class allows us to easily and + // quickly define new function + // objects via parameter files, + // with custom definitions which + // can be very complex (see the + // documentation of that class + // for all the available + // options). + // + // We will allocate the + // quadrature object using the + // QuadratureSelector class that + // allows us to generate + // quadrature formulas based on + // an identifying string and on + // the possible degree of the + // formula itself. We used this + // to allow custom selection of + // the quadrature formulas for + // the standard integration, and + // to define the order of the + // singular quadrature rule. + // + // Notice that the pointer given + // below for the quadrature rule + // is only used for non singular + // integrals. Whenever the + // integral is singular, then + // only the degree of the + // quadrature pointer is used, + // and the integration is a + // special one (see the + // assemble_matrix() function + // below for further details). + // + // We also define a couple of + // parameters which are used in + // case we wanted to extend the + // solution to the entire domain. Functions::ParsedFunction wind; Functions::ParsedFunction exact_solution; @@ -353,534 +526,696 @@ private: // @sect4{BEMProblem::BEMProblem and BEMProblem::read_parameters} -// The constructor initializes the variuous object in much the same way as -// done in the finite element programs such as step-4 or step-6. The only new -// ingredient here is the ParsedFunction object, which needs, at construction -// time, the specification of the number of components. -// -// For the exact solution the number of vector components is one, and no -// action is required since one is the default value for a ParsedFunction -// object. The wind, however, requires dim components to be specified. Notice -// that when declaring entries in a parameter file for the expression of the -// Functions::ParsedFunction, we need to specify the number of components -// explicitly, since the function -// Functions::ParsedFunction::declare_parameters is static, and has no -// knowledge of the number of components. + // The constructor initializes the + // variuous object in much the same + // way as done in the finite element + // programs such as step-4 or + // step-6. The only new ingredient + // here is the ParsedFunction object, + // which needs, at construction time, + // the specification of the number of + // components. + // + // For the exact solution the number + // of vector components is one, and + // no action is required since one is + // the default value for a + // ParsedFunction object. The wind, + // however, requires dim components + // to be specified. Notice that when + // declaring entries in a parameter + // file for the expression of the + // Functions::ParsedFunction, we need + // to specify the number of + // components explicitly, since the + // function + // Functions::ParsedFunction::declare_parameters + // is static, and has no knowledge of + // the number of components. template BEMProblem::BEMProblem() : - fe(1), - dh(tria), - wind(dim) + fe(1), + dh(tria), + wind(dim) {} template void BEMProblem::read_parameters (const std::string &filename) { - deallog << std::endl << "Parsing parameter file " << filename << std::endl - << "for a " << dim << " dimensional simulation. " << std::endl; + deallog << std::endl << "Parsing parameter file " << filename << std::endl + << "for a " << dim << " dimensional simulation. " << std::endl; - ParameterHandler prm; + ParameterHandler prm; - prm.declare_entry("Number of cycles", "4", - Patterns::Integer()); - prm.declare_entry("External refinement", "5", - Patterns::Integer()); - prm.declare_entry("Extend solution on the -2,2 box", "true", - Patterns::Bool()); - prm.declare_entry("Run 2d simulation", "true", - Patterns::Bool()); - prm.declare_entry("Run 3d simulation", "true", - Patterns::Bool()); + prm.declare_entry("Number of cycles", "4", + Patterns::Integer()); + prm.declare_entry("External refinement", "5", + Patterns::Integer()); + prm.declare_entry("Extend solution on the -2,2 box", "true", + Patterns::Bool()); + prm.declare_entry("Run 2d simulation", "true", + Patterns::Bool()); + prm.declare_entry("Run 3d simulation", "true", + Patterns::Bool()); - prm.enter_subsection("Quadrature rules"); - { - prm.declare_entry("Quadrature type", "gauss", - Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names())); - prm.declare_entry("Quadrature order", "4", Patterns::Integer()); - prm.declare_entry("Singular quadrature order", "5", Patterns::Integer()); - } - prm.leave_subsection(); + prm.enter_subsection("Quadrature rules"); + { + prm.declare_entry("Quadrature type", "gauss", + Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names())); + prm.declare_entry("Quadrature order", "4", Patterns::Integer()); + prm.declare_entry("Singular quadrature order", "5", Patterns::Integer()); + } + prm.leave_subsection(); - // For both two and three dimensions, we set the default input data to be - // such that the solution is $x+y$ or $x+y+z$. The actually computed - // solution will differ from this by a constant (remember that for the - // velocity $\mathbf{\tilde v}$ we only need the gradient of the potential - // $\phi$, so an additive constant is of no concern to us) but we will - // remove it after solving for $\phi$ to make the solution function have a - // mean value of zero. - // - // The use of the Functions::ParsedFunction object is pretty straight - // forward. The Functions::ParsedFunction::declare_parameters function - // takes an additional integer argument that specifies the number of - // components of the given function. Its default value is one. When the - // corresponding Functions::ParsedFunction::parse_parameters method is - // called, the calling object has to have the same number of components - // defined here, otherwise an exception is thrown. - // - // When declaring entries, we declare both 2 and three dimensional - // functions. However only the dim-dimensional one is ultimately - // parsed. This allows us to have only one parameter file for both 2 and 3 - // dimensional problems. - prm.enter_subsection("Wind function 2d"); - { - Functions::ParsedFunction<2>::declare_parameters(prm, 2); - prm.set("Function expression", "1; 1"); - } - prm.leave_subsection(); + // For both two and three + // dimensions, we set the default + // input data to be such that the + // solution is $x+y$ or + // $x+y+z$. The actually computed + // solution will differ from this + // by a constant (remember that for + // the velocity $\mathbf{\tilde v}$ + // we only need the gradient of the + // potential $\phi$, so an additive + // constant is of no concern to us) + // but we will remove it after + // solving for $\phi$ to make the + // solution function have a mean + // value of zero. + // + // The use of the + // Functions::ParsedFunction object + // is pretty straight forward. The + // Functions::ParsedFunction::declare_parameters + // function takes an additional + // integer argument that specifies + // the number of components of the + // given function. Its default + // value is one. When the + // corresponding + // Functions::ParsedFunction::parse_parameters + // method is called, the calling + // object has to have the same + // number of components defined + // here, otherwise an exception is + // thrown. + // + // When declaring entries, we + // declare both 2 and three + // dimensional functions. However + // only the dim-dimensional one is + // ultimately parsed. This allows + // us to have only one parameter + // file for both 2 and 3 + // dimensional problems. + prm.enter_subsection("Wind function 2d"); + { + Functions::ParsedFunction<2>::declare_parameters(prm, 2); + prm.set("Function expression", "1; 1"); + } + prm.leave_subsection(); - prm.enter_subsection("Wind function 3d"); - { - Functions::ParsedFunction<3>::declare_parameters(prm, 3); - prm.set("Function expression", "1; 1; 1"); - } - prm.leave_subsection(); + prm.enter_subsection("Wind function 3d"); + { + Functions::ParsedFunction<3>::declare_parameters(prm, 3); + prm.set("Function expression", "1; 1; 1"); + } + prm.leave_subsection(); - prm.enter_subsection("Exact solution 2d"); - { - Functions::ParsedFunction<2>::declare_parameters(prm); - prm.set("Function expression", "x+y"); - } - prm.leave_subsection(); + prm.enter_subsection("Exact solution 2d"); + { + Functions::ParsedFunction<2>::declare_parameters(prm); + prm.set("Function expression", "x+y"); + } + prm.leave_subsection(); - prm.enter_subsection("Exact solution 3d"); - { - Functions::ParsedFunction<3>::declare_parameters(prm); - prm.set("Function expression", "x+y+z"); - } - prm.leave_subsection(); - - // After declaring all these parameters - // to the ParameterHandler object, let's - // read an input file that will give the - // parameters their values. We then - // proceed to extract these values from - // the ParameterHandler object: - prm.read_input(filename); - - n_cycles = prm.get_integer("Number of cycles"); - external_refinement = prm.get_integer("External refinement"); - extend_solution = prm.get_bool("Extend solution on the -2,2 box"); + prm.enter_subsection("Exact solution 3d"); + { + Functions::ParsedFunction<3>::declare_parameters(prm); + prm.set("Function expression", "x+y+z"); + } + prm.leave_subsection(); + + // After declaring all these + // parameters to the + // ParameterHandler object, let's + // read an input file that will + // give the parameters their + // values. We then proceed to + // extract these values from the + // ParameterHandler object: + prm.read_input(filename); + + n_cycles = prm.get_integer("Number of cycles"); + external_refinement = prm.get_integer("External refinement"); + extend_solution = prm.get_bool("Extend solution on the -2,2 box"); - prm.enter_subsection("Quadrature rules"); - { - quadrature = - std_cxx0x::shared_ptr > - (new QuadratureSelector (prm.get("Quadrature type"), - prm.get_integer("Quadrature order"))); - singular_quadrature_order = prm.get_integer("Singular quadrature order"); - } - prm.leave_subsection(); + prm.enter_subsection("Quadrature rules"); + { + quadrature = + std_cxx0x::shared_ptr > + (new QuadratureSelector (prm.get("Quadrature type"), + prm.get_integer("Quadrature order"))); + singular_quadrature_order = prm.get_integer("Singular quadrature order"); + } + prm.leave_subsection(); - prm.enter_subsection(std::string("Wind function ")+ - Utilities::int_to_string(dim)+std::string("d")); - { - wind.parse_parameters(prm); - } - prm.leave_subsection(); - - prm.enter_subsection(std::string("Exact solution ")+ - Utilities::int_to_string(dim)+std::string("d")); - { - exact_solution.parse_parameters(prm); - } - prm.leave_subsection(); + prm.enter_subsection(std::string("Wind function ")+ + Utilities::int_to_string(dim)+std::string("d")); + { + wind.parse_parameters(prm); + } + prm.leave_subsection(); - // Finally, here's another example of how to use parameter files in - // dimension independent programming. If we wanted to switch off one of - // the two simulations, we could do this by setting the corresponding "Run - // 2d simulation" or "Run 3d simulation" flag to false: - run_in_this_dimension = prm.get_bool("Run " + - Utilities::int_to_string(dim) + - "d simulation"); + prm.enter_subsection(std::string("Exact solution ")+ + Utilities::int_to_string(dim)+std::string("d")); + { + exact_solution.parse_parameters(prm); + } + prm.leave_subsection(); + + // Finally, here's another example + // of how to use parameter files in + // dimension independent + // programming. If we wanted to + // switch off one of the two + // simulations, we could do this by + // setting the corresponding "Run + // 2d simulation" or "Run 3d + // simulation" flag to false: + run_in_this_dimension = prm.get_bool("Run " + + Utilities::int_to_string(dim) + + "d simulation"); } // @sect4{BEMProblem::read_domain} - // A boundary element method triangulation is basically the same - // as a (dim-1) dimensional triangulation, with the difference that the - // vertices belong to a (dim) dimensional space. - // - // Some of the mesh formats supported in deal.II use by default three - // dimensional points to describe meshes. These are the formats which are - // compatible with the boundary element method capabilities of deal.II. In - // particular we can use either UCD or GMSH formats. In both cases, we - // have to be particularly careful with the orientation of the mesh, - // because, unlike in the standard finite element case, no reordering or - // compatibility check is performed here. All meshes are considered as - // oriented, because they are embedded in a higher dimensional space. (See - // the documentation of the GridIn and of the Triangulation for further - // details on orientation of cells in a triangulation.) In our case, the - // normals to the mesh are external to both the circle in 2d or the sphere - // in 3d. - // - // The other detail that is required for appropriate refinement of the - // boundary element mesh, is an accurate description of the manifold that - // the mesh is approximating. We already saw this several times for the - // boundary of standard finite element meshes (for example in step-5 and - // step-6), and here the principle and usage is the same, except that the - // HyperBallBoundary class takes an additional template parameter that - // specifies the embedding space dimension. The function object still has - // to be static to live at least as long as the triangulation object to - // which it is attached. + // A boundary element method + // triangulation is basically the + // same as a (dim-1) dimensional + // triangulation, with the difference + // that the vertices belong to a + // (dim) dimensional space. + // + // Some of the mesh formats supported + // in deal.II use by default three + // dimensional points to describe + // meshes. These are the formats + // which are compatible with the + // boundary element method + // capabilities of deal.II. In + // particular we can use either UCD + // or GMSH formats. In both cases, we + // have to be particularly careful + // with the orientation of the mesh, + // because, unlike in the standard + // finite element case, no reordering + // or compatibility check is + // performed here. All meshes are + // considered as oriented, because + // they are embedded in a higher + // dimensional space. (See the + // documentation of the GridIn and of + // the Triangulation for further + // details on orientation of cells in + // a triangulation.) In our case, the + // normals to the mesh are external + // to both the circle in 2d or the + // sphere in 3d. + // + // The other detail that is required + // for appropriate refinement of the + // boundary element mesh, is an + // accurate description of the + // manifold that the mesh is + // approximating. We already saw this + // several times for the boundary of + // standard finite element meshes + // (for example in step-5 and + // step-6), and here the principle + // and usage is the same, except that + // the HyperBallBoundary class takes + // an additional template parameter + // that specifies the embedding space + // dimension. The function object + // still has to be static to live at + // least as long as the triangulation + // object to which it is attached. template void BEMProblem::read_domain() { - static HyperBallBoundary boundary(Point(),1.); - - std::ifstream in; - switch (dim) - { - case 2: - in.open ("coarse_circle.inp"); - break; + static HyperBallBoundary boundary(Point(),1.); + + std::ifstream in; + switch (dim) + { + case 2: + in.open ("coarse_circle.inp"); + break; - case 3: - in.open ("coarse_sphere.inp"); - break; + case 3: + in.open ("coarse_sphere.inp"); + break; - default: - Assert (false, ExcNotImplemented()); - } + default: + Assert (false, ExcNotImplemented()); + } - GridIn gi; - gi.attach_triangulation (tria); - gi.read_ucd (in); - tria.set_boundary(1, boundary); + GridIn gi; + gi.attach_triangulation (tria); + gi.read_ucd (in); + tria.set_boundary(1, boundary); } // @sect4{BEMProblem::refine_and_resize} - // This function globally refines the mesh, - // distributes degrees of freedom, and - // resizes matrices and vectors. + // This function globally refines the + // mesh, distributes degrees of + // freedom, and resizes matrices and + // vectors. // - // Note that the matrix is a full matrix and that consequently we have to - // build a sparsity pattern that contains every single - // entry. Notwithstanding this fact, the SparseMatrix class coupled with - // the SparseDirectUMFPACK solver are still faster than Lapack solvers for - // full matrices. The drawback is that we need to assemble a full - // SparsityPattern, which is not the most efficient way to store a full - // matrix. + // Note that the matrix is a full + // matrix and that consequently we + // have to build a sparsity pattern + // that contains every single + // entry. Notwithstanding this fact, + // the SparseMatrix class coupled + // with the SparseDirectUMFPACK + // solver are still faster than + // Lapack solvers for full + // matrices. The drawback is that we + // need to assemble a full + // SparsityPattern, which is not the + // most efficient way to store a full + // matrix. template void BEMProblem::refine_and_resize() { - tria.refine_global(1); + tria.refine_global(1); - dh.distribute_dofs(fe); + dh.distribute_dofs(fe); - const unsigned int n_dofs = dh.n_dofs(); + const unsigned int n_dofs = dh.n_dofs(); - system_matrix.clear(); - sparsity.reinit(n_dofs, n_dofs, n_dofs); - for(unsigned int i=0; i void BEMProblem::assemble_system() { - std::vector > sing_quadratures_3d; - for(unsigned int i=0; i<4; ++i) { - sing_quadratures_3d.push_back - (QGaussOneOverR<2>(singular_quadrature_order, i, true)); - } + std::vector > sing_quadratures_3d; + for(unsigned int i=0; i<4; ++i) { + sing_quadratures_3d.push_back + (QGaussOneOverR<2>(singular_quadrature_order, i, true)); + } - // Next, we initialize an FEValues object with the quadrature formula for - // the integration of the kernel in non singular cells. This quadrature is - // selected with the parameter file, and needs to be quite precise, since - // the functions we are integrating are not polynomial functions. - FEValues fe_v(fe, *quadrature, - update_values | - update_cell_normal_vectors | - update_quadrature_points | - update_JxW_values); + // Next, we initialize an FEValues + // object with the quadrature + // formula for the integration of + // the kernel in non singular + // cells. This quadrature is + // selected with the parameter + // file, and needs to be quite + // precise, since the functions we + // are integrating are not + // polynomial functions. + FEValues fe_v(fe, *quadrature, + update_values | + update_cell_normal_vectors | + update_quadrature_points | + update_JxW_values); - const unsigned int n_q_points = fe_v.n_quadrature_points; + const unsigned int n_q_points = fe_v.n_quadrature_points; - std::vector local_dof_indices(fe.dofs_per_cell); + std::vector local_dof_indices(fe.dofs_per_cell); - std::vector > cell_wind(n_q_points, Vector(dim) ); - double normal_wind; + std::vector > cell_wind(n_q_points, Vector(dim) ); + double normal_wind; - // Unlike in finite element methods, if we use a collocation boundary - // element method, then in each assembly loop we only assemble the - // information that refers to the coupling between one degree of freedom - // (the degree associated with support point $i$) and the current - // cell. This is done using a vector of fe.dofs_per_cell elements, which - // will then be distributed to the matrix in the global row $i$. The - // following object will hold this information: - Vector local_matrix_row_i(fe.dofs_per_cell); + // Unlike in finite element + // methods, if we use a collocation + // boundary element method, then in + // each assembly loop we only + // assemble the information that + // refers to the coupling between + // one degree of freedom (the + // degree associated with support + // point $i$) and the current + // cell. This is done using a + // vector of fe.dofs_per_cell + // elements, which will then be + // distributed to the matrix in the + // global row $i$. The following + // object will hold this + // information: + Vector local_matrix_row_i(fe.dofs_per_cell); - // The index $i$ runs on the collocation points, which are the support - // points of the $i$th basis function, while $j$ runs on inner integration - // points. We perform the following check to ensure that we are not trying - // to use this code for high order elements. It will only work with Q1 - // elements, that is, for fe.dofs_per_cell == - // GeometryInfo::vertices_per_cell. - AssertThrow(fe.dofs_per_cell == GeometryInfo::vertices_per_cell, - ExcMessage("The code in this function can only be used for " - "the usual Q1 elements.")); + // The index $i$ runs on the + // collocation points, which are + // the support points of the $i$th + // basis function, while $j$ runs + // on inner integration points. We + // perform the following check to + // ensure that we are not trying to + // use this code for high order + // elements. It will only work with + // Q1 elements, that is, for + // fe.dofs_per_cell == + // GeometryInfo::vertices_per_cell. + AssertThrow(fe.dofs_per_cell == GeometryInfo::vertices_per_cell, + ExcMessage("The code in this function can only be used for " + "the usual Q1 elements.")); - // Now that we have checked that the number of vertices is equal to the - // number of degrees of freedom, we construct a vector of support - // points which will be used in the local integrations: - std::vector > support_points(dh.n_dofs()); - DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, - dh, support_points); - - // After doing so, we can start the - // integration loop over all cells, where - // we first initialize the FEValues - // object and get the values of - // $\mathbf{\tilde v}$ at the quadrature - // points (this vector field should be - // constant, but it doesn't hurt to be - // more general): - typename DoFHandler::active_cell_iterator - cell = dh.begin_active(), - endc = dh.end(); + // Now that we have checked that + // the number of vertices is equal + // to the number of degrees of + // freedom, we construct a vector + // of support points which will be + // used in the local integrations: + std::vector > support_points(dh.n_dofs()); + DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, + dh, support_points); + + // After doing so, we can start the + // integration loop over all cells, + // where we first initialize the + // FEValues object and get the + // values of $\mathbf{\tilde v}$ at + // the quadrature points (this + // vector field should be constant, + // but it doesn't hurt to be more + // general): + typename DoFHandler::active_cell_iterator + cell = dh.begin_active(), + endc = dh.end(); - for(cell = dh.begin_active(); cell != endc; ++cell) { + for(cell = dh.begin_active(); cell != endc; ++cell) { - fe_v.reinit(cell); - cell->get_dof_indices(local_dof_indices); + fe_v.reinit(cell); + cell->get_dof_indices(local_dof_indices); - const std::vector > &q_points = fe_v.get_quadrature_points(); - const std::vector > &normals = fe_v.get_cell_normal_vectors(); - wind.vector_value_list(q_points, cell_wind); + const std::vector > &q_points = fe_v.get_quadrature_points(); + const std::vector > &normals = fe_v.get_cell_normal_vectors(); + wind.vector_value_list(q_points, cell_wind); - // We then form the integral over the current cell for all degrees - // of freedom (note that this includes degrees of freedom not - // located on the current cell, a deviation from the usual finite - // element integrals). The integral that we need to perform is - // singular if one of the local degrees of freedom is the same as - // the support point $i$. A the beginning of the loop we therefore - // check wether this is the case, and we store which one is the - // singular index: - for(unsigned int i=0; i R = q_points[q] - support_points[i]; + const Point R = q_points[q] - support_points[i]; - system_rhs(i) += ( LaplaceKernel::single_layer(R) * - normal_wind * - fe_v.JxW(q) ); + system_rhs(i) += ( LaplaceKernel::single_layer(R) * + normal_wind * + fe_v.JxW(q) ); - for(unsigned int j=0; j or a - // QGaussOneOverR<2> to a - // Quadrature. C++ - // doesn't allow this right - // away, and neither is a - // static_cast - // possible. However, we can - // attempt a dynamic_cast: the - // implementation will then - // look up at run time whether - // the conversion is possible - // (which we know it - // is) and if that isn't the - // case simply return a null - // pointer. To be sure we can - // then add a safety check at - // the end: - Assert(singular_index != numbers::invalid_unsigned_int, - ExcInternalError()); - - const Quadrature * - singular_quadrature - = (dim == 2 - ? - dynamic_cast*>( - new QGaussLogR<1>(singular_quadrature_order, - Point<1>((double)singular_index), - 1./cell->measure())) - : - (dim == 3 - ? - dynamic_cast*>( - &sing_quadratures_3d[singular_index]) - : - 0)); - Assert(singular_quadrature, ExcInternalError()); + local_matrix_row_i(j) += ( ( LaplaceKernel::double_layer(R) * + normals[q] ) * + fe_v.shape_value(j,q) * + fe_v.JxW(q) ); + } + } + } else { + // Now we treat the more + // delicate case. If we are + // here, this means that the + // cell that runs on the $j$ + // index contains + // support_point[i]. In this + // case both the single and + // the double layer potential + // are singular, and they + // require special treatment, + // as explained in the + // introduction. + // + // In the two dimensional + // case we perform the + // integration using a + // QGaussLogR quadrature + // formula, which is + // specifically designed to + // integrate logarithmic + // singularities on the unit + // interval, while in three + // dimensions we use the + // QGaussOneOverR class, + // which allows us to + // integrate 1/R + // singularities on the + // vertices of the reference + // element. Since we don't + // want to rebuild the two + // dimensional quadrature + // formula at each singular + // integration, we have built + // them outside the loop on + // the cells, and we only use + // a pointer to that + // quadrature here. + // + // Notice that in one + // dimensional integration + // this is not possible, + // since we need to know the + // scaling parameter for the + // quadrature, which is not + // known a priori. Here, the + // singular quadrature rule + // depends also on the size + // of the current cell. For + // this reason, it is + // necessary to create a new + // quadrature for each + // singular + // integration. Since we + // create it using the new + // operator of C++, we also + // need to destroy it using + // the dual of new: + // delete. This is done at + // the end, and only if dim + // == 2. + // + // Putting all this into a + // dimension independent + // framework requires a little + // trick. The problem is that, + // depending on dimension, we'd + // like to either assign a + // QGaussLogR<1> or a + // QGaussOneOverR<2> to a + // Quadrature. C++ + // doesn't allow this right + // away, and neither is a + // static_cast + // possible. However, we can + // attempt a dynamic_cast: the + // implementation will then + // look up at run time whether + // the conversion is possible + // (which we know it + // is) and if that isn't the + // case simply return a null + // pointer. To be sure we can + // then add a safety check at + // the end: + Assert(singular_index != numbers::invalid_unsigned_int, + ExcInternalError()); + + const Quadrature * + singular_quadrature + = (dim == 2 + ? + dynamic_cast*>( + new QGaussLogR<1>(singular_quadrature_order, + Point<1>((double)singular_index), + 1./cell->measure())) + : + (dim == 3 + ? + dynamic_cast*>( + &sing_quadratures_3d[singular_index]) + : + 0)); + Assert(singular_quadrature, ExcInternalError()); - FEValues fe_v_singular (fe, *singular_quadrature, - update_jacobians | - update_values | - update_cell_normal_vectors | - update_quadrature_points ); + FEValues fe_v_singular (fe, *singular_quadrature, + update_jacobians | + update_values | + update_cell_normal_vectors | + update_quadrature_points ); - fe_v_singular.reinit(cell); + fe_v_singular.reinit(cell); - std::vector > singular_cell_wind( (*singular_quadrature).size(), - Vector(dim) ); + std::vector > singular_cell_wind( (*singular_quadrature).size(), + Vector(dim) ); - const std::vector > &singular_normals = fe_v_singular.get_cell_normal_vectors(); - const std::vector > &singular_q_points = fe_v_singular.get_quadrature_points(); + const std::vector > &singular_normals = fe_v_singular.get_cell_normal_vectors(); + const std::vector > &singular_q_points = fe_v_singular.get_quadrature_points(); - wind.vector_value_list(singular_q_points, singular_cell_wind); + wind.vector_value_list(singular_q_points, singular_cell_wind); - for(unsigned int q=0; qsize(); ++q) { - const Point R = singular_q_points[q]- support_points[i]; - double normal_wind = 0; - for(unsigned int d=0; dsize(); ++q) { + const Point R = singular_q_points[q]- support_points[i]; + double normal_wind = 0; + for(unsigned int d=0; d ones(dh.n_dofs()); - ones.add(-1.); - - system_matrix.vmult(alpha, ones); - for(unsigned int i = 0; i ones(dh.n_dofs()); + ones.add(-1.); + + system_matrix.vmult(alpha, ones); + for(unsigned int i = 0; i::assemble_system() { // vector entry to normalize it. template void BEMProblem::solve_system() { - SparseDirectUMFPACK inverse_matrix; - inverse_matrix.initialize (system_matrix); - inverse_matrix.vmult (phi, system_rhs); + SparseDirectUMFPACK inverse_matrix; + inverse_matrix.initialize (system_matrix); + inverse_matrix.vmult (phi, system_rhs); - phi.add(-phi.mean_value()); + phi.add(-phi.mean_value()); } // @sect4{BEMProblem::solve_system} - // The computation of the errors is exactly the same in all other - // example programs, and we won't comment too much. Notice how the - // same methods that are used in the finite element methods can be - // used here. + // The computation of the errors is + // exactly the same in all other + // example programs, and we won't + // comment too much. Notice how the + // same methods that are used in the + // finite element methods can be used + // here. template void BEMProblem::compute_errors(const unsigned int cycle) { - Vector difference_per_cell (tria.n_active_cells()); - VectorTools::integrate_difference (dh, phi, - exact_solution, - difference_per_cell, - QGauss<(dim-1)>(3), - VectorTools::L2_norm); - const double L2_error = difference_per_cell.l2_norm(); + Vector difference_per_cell (tria.n_active_cells()); + VectorTools::integrate_difference (dh, phi, + exact_solution, + difference_per_cell, + QGauss<(dim-1)>(3), + VectorTools::L2_norm); + const double L2_error = difference_per_cell.l2_norm(); - // The error in the alpha vector can be computed directly using the - // Vector::linfty_norm() function, since on each node, the value should be - // $\frac 12$. All errors are then output and appended to our - // ConvergenceTable object for later computation of convergence rates: - Vector difference_per_node(alpha); - difference_per_node.add(-.5); + // The error in the alpha vector + // can be computed directly using + // the Vector::linfty_norm() + // function, since on each node, + // the value should be $\frac + // 12$. All errors are then output + // and appended to our + // ConvergenceTable object for + // later computation of convergence + // rates: + Vector difference_per_node(alpha); + difference_per_node.add(-.5); - const double alpha_error = difference_per_node.linfty_norm(); - const unsigned int n_active_cells=tria.n_active_cells(); - const unsigned int n_dofs=dh.n_dofs(); + const double alpha_error = difference_per_node.linfty_norm(); + const unsigned int n_active_cells=tria.n_active_cells(); + const unsigned int n_dofs=dh.n_dofs(); - deallog << "Cycle " << cycle << ':' - << std::endl - << " Number of active cells: " - << n_active_cells - << std::endl - << " Number of degrees of freedom: " - << n_dofs - << std::endl; + deallog << "Cycle " << cycle << ':' + << std::endl + << " Number of active cells: " + << n_active_cells + << std::endl + << " Number of degrees of freedom: " + << n_dofs + << std::endl; - convergence_table.add_value("cycle", cycle); - convergence_table.add_value("cells", n_active_cells); - convergence_table.add_value("dofs", n_dofs); - convergence_table.add_value("L2(phi)", L2_error); - convergence_table.add_value("Linfty(alpha)", alpha_error); + convergence_table.add_value("cycle", cycle); + convergence_table.add_value("cells", n_active_cells); + convergence_table.add_value("dofs", n_dofs); + convergence_table.add_value("L2(phi)", L2_error); + convergence_table.add_value("Linfty(alpha)", alpha_error); } // @sect4{BEMProblem::compute_exterior_solution} - // We'd like to also know something about the - // value of the potential $\phi$ in the - // exterior domain: after all our motivation - // to consider the boundary integral problem - // was that we wanted to know the velocity in - // the exterior domain! + // We'd like to also know something + // about the value of the potential + // $\phi$ in the exterior domain: + // after all our motivation to + // consider the boundary integral + // problem was that we wanted to know + // the velocity in the exterior + // domain! // - // To this end, let us assume here that the boundary element domain is - // contained in the box $[-2,2]^{\text{dim}}$, and we extrapolate the actual - // solution inside this box using the convolution with the fundamental - // solution. The formula for this is given in the introduction. - // - // The reconstruction of the solution in the entire space is done on a - // continuous finite element grid of dimension dim. These are the usual - // ones, and we don't comment any further on them. At the end of the - // function, we output this exterior solution in, again, much the usual - // way. + // To this end, let us assume here + // that the boundary element domain + // is contained in the box + // $[-2,2]^{\text{dim}}$, and we + // extrapolate the actual solution + // inside this box using the + // convolution with the fundamental + // solution. The formula for this is + // given in the introduction. + // + // The reconstruction of the solution + // in the entire space is done on a + // continuous finite element grid of + // dimension dim. These are the usual + // ones, and we don't comment any + // further on them. At the end of the + // function, we output this exterior + // solution in, again, much the usual + // way. template void BEMProblem::compute_exterior_solution() { - Triangulation external_tria; - GridGenerator::hyper_cube(external_tria, -2, 2); + Triangulation external_tria; + GridGenerator::hyper_cube(external_tria, -2, 2); - FE_Q external_fe(1); - DoFHandler external_dh (external_tria); - Vector external_phi; + FE_Q external_fe(1); + DoFHandler external_dh (external_tria); + Vector external_phi; - external_tria.refine_global(external_refinement); - external_dh.distribute_dofs(external_fe); - external_phi.reinit(external_dh.n_dofs()); + external_tria.refine_global(external_refinement); + external_dh.distribute_dofs(external_fe); + external_phi.reinit(external_dh.n_dofs()); - typename DoFHandler::active_cell_iterator - cell = dh.begin_active(), - endc = dh.end(); + typename DoFHandler::active_cell_iterator + cell = dh.begin_active(), + endc = dh.end(); - FEValues fe_v(fe, *quadrature, - update_values | - update_cell_normal_vectors | - update_quadrature_points | - update_JxW_values); + FEValues fe_v(fe, *quadrature, + update_values | + update_cell_normal_vectors | + update_quadrature_points | + update_JxW_values); - const unsigned int n_q_points = fe_v.n_quadrature_points; + const unsigned int n_q_points = fe_v.n_quadrature_points; - std::vector dofs(fe.dofs_per_cell); + std::vector dofs(fe.dofs_per_cell); - std::vector local_phi(n_q_points); - std::vector normal_wind(n_q_points); - std::vector > local_wind(n_q_points, Vector(dim) ); + std::vector local_phi(n_q_points); + std::vector normal_wind(n_q_points); + std::vector > local_wind(n_q_points, Vector(dim) ); - typename DoFHandler::active_cell_iterator - external_cell = external_dh.begin_active(), - external_endc = external_dh.end(); + typename DoFHandler::active_cell_iterator + external_cell = external_dh.begin_active(), + external_endc = external_dh.end(); - std::vector > external_support_points(external_dh.n_dofs()); - DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, - external_dh, external_support_points); + std::vector > external_support_points(external_dh.n_dofs()); + DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, + external_dh, external_support_points); - for(cell = dh.begin_active(); cell != endc; ++cell) { - fe_v.reinit(cell); + for(cell = dh.begin_active(); cell != endc; ++cell) { + fe_v.reinit(cell); - const std::vector > &q_points = fe_v.get_quadrature_points(); - const std::vector > &normals = fe_v.get_cell_normal_vectors(); + const std::vector > &q_points = fe_v.get_quadrature_points(); + const std::vector > &normals = fe_v.get_cell_normal_vectors(); - cell->get_dof_indices(dofs); - fe_v.get_function_values(phi, local_phi); + cell->get_dof_indices(dofs); + fe_v.get_function_values(phi, local_phi); - wind.vector_value_list(q_points, local_wind); + wind.vector_value_list(q_points, local_wind); - for(unsigned int q=0; q R = q_points[q] - external_support_points[i]; + const Point R = q_points[q] - external_support_points[i]; - external_phi(i) += ( ( LaplaceKernel::single_layer(R) * - normal_wind[q] - + - (LaplaceKernel::double_layer(R) * - normals[q] ) * - local_phi[q] ) * - fe_v.JxW(q) ); - } - } + external_phi(i) += ( ( LaplaceKernel::single_layer(R) * + normal_wind[q] + + + (LaplaceKernel::double_layer(R) * + normals[q] ) * + local_phi[q] ) * + fe_v.JxW(q) ); + } } + } - DataOut data_out; + DataOut data_out; - data_out.attach_dof_handler(external_dh); - data_out.add_data_vector(external_phi, "external_phi"); - data_out.build_patches(); + data_out.attach_dof_handler(external_dh); + data_out.add_data_vector(external_phi, "external_phi"); + data_out.build_patches(); - const std::string - filename = Utilities::int_to_string(dim) + "d_external.vtk"; - std::ofstream file(filename.c_str()); + const std::string + filename = Utilities::int_to_string(dim) + "d_external.vtk"; + std::ofstream file(filename.c_str()); - data_out.write_vtk(file); + data_out.write_vtk(file); } // @sect4{BEMProblem::output_results} - // Outputting the results of our computations - // is a rather mechanical tasks. All the - // components of this function have been - // discussed before. + // Outputting the results of our + // computations is a rather + // mechanical tasks. All the + // components of this function have + // been discussed before. template void BEMProblem::output_results(const unsigned int cycle) { - DataOut > dataout; + DataOut > dataout; - dataout.attach_dof_handler(dh); - dataout.add_data_vector(phi, "phi"); - dataout.add_data_vector(alpha, "alpha"); - dataout.build_patches(); + dataout.attach_dof_handler(dh); + dataout.add_data_vector(phi, "phi"); + dataout.add_data_vector(alpha, "alpha"); + dataout.build_patches(); - std::string filename = ( Utilities::int_to_string(dim) + - "d_boundary_solution_" + - Utilities::int_to_string(cycle) + - ".vtk" ); - std::ofstream file(filename.c_str()); + std::string filename = ( Utilities::int_to_string(dim) + + "d_boundary_solution_" + + Utilities::int_to_string(cycle) + + ".vtk" ); + std::ofstream file(filename.c_str()); - dataout.write_vtk(file); + dataout.write_vtk(file); - if(cycle == n_cycles-1) { - convergence_table.set_precision("L2(phi)", 3); - convergence_table.set_precision("Linfty(alpha)", 3); + if(cycle == n_cycles-1) { + convergence_table.set_precision("L2(phi)", 3); + convergence_table.set_precision("Linfty(alpha)", 3); - convergence_table.set_scientific("L2(phi)", true); - convergence_table.set_scientific("Linfty(alpha)", true); + convergence_table.set_scientific("L2(phi)", true); + convergence_table.set_scientific("Linfty(alpha)", true); - convergence_table - .evaluate_convergence_rates("L2(phi)", ConvergenceTable::reduction_rate_log2); - convergence_table - .evaluate_convergence_rates("Linfty(alpha)", ConvergenceTable::reduction_rate_log2); - deallog << std::endl; - convergence_table.write_text(std::cout); - } + convergence_table + .evaluate_convergence_rates("L2(phi)", ConvergenceTable::reduction_rate_log2); + convergence_table + .evaluate_convergence_rates("Linfty(alpha)", ConvergenceTable::reduction_rate_log2); + deallog << std::endl; + convergence_table.write_text(std::cout); + } } // @sect4{BEMProblem::run} - // This is the main function. It should be - // self explanatory in its briefness: + // This is the main function. It + // should be self explanatory in its + // briefness: template void BEMProblem::run() { - read_parameters("parameters.prm"); + read_parameters("parameters.prm"); - if(run_in_this_dimension == false) - { - deallog << "Run in dimension " << dim - << " explicitly disabled in parameter file. " - << std::endl; - return; - } + if(run_in_this_dimension == false) + { + deallog << "Run in dimension " << dim + << " explicitly disabled in parameter file. " + << std::endl; + return; + } - read_domain(); + read_domain(); - for(unsigned int cycle=0; cycle::run() { int main () { try - { + { deallog.depth_console (3); BEMProblem<2> laplace_problem_2d; laplace_problem_2d.run(); BEMProblem<3> laplace_problem_3d; laplace_problem_3d.run(); - } + } catch (std::exception &exc) { std::cerr << std::endl << std::endl