// @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 <int dim>
-double single_layer(const Point<dim> &R,
- const bool factor_out_2d_singularity = false) {
+ template <int dim>
+ double single_layer(const Point<dim> &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 <int dim>
-Point<dim> double_layer(const Point<dim> &R,
- const bool factor_out_2d_singularity = false) {
- switch(dim) {
- case 2:
- if (factor_out_2d_singularity)
- return Point<dim>();
- else
- return R / (-2*numbers::PI * R.square());
- case 3:
- return R / ( -4*numbers::PI * R.square()*R.norm() );
+ template <int dim>
+ Point<dim> double_layer(const Point<dim> &R,
+ const bool factor_out_2d_singularity = false) {
+ switch(dim) {
+ case 2:
+ if (factor_out_2d_singularity)
+ return Point<dim>();
+ 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<dim>();
+ Assert(false, ExcInternalError());
+ return Point<dim>();
}
-}
+ }
}
// @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 <int dim>
class BEMProblem
{
-public:
+ public:
BEMProblem();
void run();
-private:
+ private:
void read_parameters (const std::string &filename);
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
- // <code>solve_system()</code> 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
+ // <code>solve_system()</code>
+ // 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<dim-1, dim> tria;
FE_Q<dim-1,dim> fe;
DoFHandler<dim-1,dim> 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<double> system_matrix;
Vector<double> 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<double> phi;
Vector<double> 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<dim> wind;
Functions::ParsedFunction<dim> exact_solution;
// @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 <int dim>
BEMProblem<dim>::BEMProblem() :
- fe(1),
- dh(tria),
- wind(dim)
+ fe(1),
+ dh(tria),
+ wind(dim)
{}
template <int dim>
void BEMProblem<dim>::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<Quadrature<dim-1> >
- (new QuadratureSelector<dim-1> (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<Quadrature<dim-1> >
+ (new QuadratureSelector<dim-1> (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 <int dim>
void BEMProblem<dim>::read_domain() {
- static HyperBallBoundary<dim-1, dim> boundary(Point<dim>(),1.);
-
- std::ifstream in;
- switch (dim)
- {
- case 2:
- in.open ("coarse_circle.inp");
- break;
+ static HyperBallBoundary<dim-1, dim> boundary(Point<dim>(),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<dim-1, dim> gi;
- gi.attach_triangulation (tria);
- gi.read_ucd (in);
- tria.set_boundary(1, boundary);
+ GridIn<dim-1, dim> 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 <int dim>
void BEMProblem<dim>::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<n_dofs;++i)
- for(unsigned int j=0; j<n_dofs; ++j)
- sparsity.add(i,j);
- sparsity.compress();
- system_matrix.reinit(sparsity);
+ system_matrix.clear();
+ sparsity.reinit(n_dofs, n_dofs, n_dofs);
+ for(unsigned int i=0; i<n_dofs;++i)
+ for(unsigned int j=0; j<n_dofs; ++j)
+ sparsity.add(i,j);
+ sparsity.compress();
+ system_matrix.reinit(sparsity);
- system_rhs.reinit(n_dofs);
- phi.reinit(n_dofs);
- alpha.reinit(n_dofs);
+ system_rhs.reinit(n_dofs);
+ phi.reinit(n_dofs);
+ alpha.reinit(n_dofs);
}
// @sect4{BEMProblem::assemble_system}
- // The following is the main function of this
- // program, assembling the matrix that
- // corresponds to the boundary integral
- // equation.
+ // The following is the main function
+ // of this program, assembling the
+ // matrix that corresponds to the
+ // boundary integral equation.
//
- // At the beginning, we create the singular quadratures for the three
- // dimensional problem (note that a 3d boundary integral problem requires
- // a 2d quadrature formula!), since in this case they only depend on the
- // reference element. This quadrature is a standard Gauss quadrature
- // formula reparametrized in such a way that allows one to integrate
- // singularities of the kind $1/R$ centered at one of the vertices. Here
- // we define a vector of four such quadratures (one per vertex of the two
- // dimensional cells for a surface in 3d) that will be used later on;
- // note, however, that these objects will only be used in the three
- // dimensional case.
+ // At the beginning, we create the
+ // singular quadratures for the three
+ // dimensional problem (note that a
+ // 3d boundary integral problem
+ // requires a 2d quadrature
+ // formula!), since in this case they
+ // only depend on the reference
+ // element. This quadrature is a
+ // standard Gauss quadrature formula
+ // reparametrized in such a way that
+ // allows one to integrate
+ // singularities of the kind $1/R$
+ // centered at one of the
+ // vertices. Here we define a vector
+ // of four such quadratures (one per
+ // vertex of the two dimensional
+ // cells for a surface in 3d) that
+ // will be used later on; note,
+ // however, that these objects will
+ // only be used in the three
+ // dimensional case.
template <int dim>
void BEMProblem<dim>::assemble_system() {
- std::vector<QGaussOneOverR<2> > 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<QGaussOneOverR<2> > 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<dim-1,dim> 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<dim-1,dim> 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<unsigned int> local_dof_indices(fe.dofs_per_cell);
+ std::vector<unsigned int> local_dof_indices(fe.dofs_per_cell);
- std::vector<Vector<double> > cell_wind(n_q_points, Vector<double>(dim) );
- double normal_wind;
+ std::vector<Vector<double> > cell_wind(n_q_points, Vector<double>(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<double> 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<double> 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<dim>::vertices_per_cell.
- AssertThrow(fe.dofs_per_cell == GeometryInfo<dim-1>::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<dim>::vertices_per_cell.
+ AssertThrow(fe.dofs_per_cell == GeometryInfo<dim-1>::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<Point<dim> > support_points(dh.n_dofs());
- DoFTools::map_dofs_to_support_points<dim-1, dim>( StaticMappingQ1<dim-1, dim>::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<dim-1,dim>::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<Point<dim> > support_points(dh.n_dofs());
+ DoFTools::map_dofs_to_support_points<dim-1, dim>( StaticMappingQ1<dim-1, dim>::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<dim-1,dim>::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<Point<dim> > &q_points = fe_v.get_quadrature_points();
-
const std::vector<Point<dim> > &normals = fe_v.get_cell_normal_vectors();
- wind.vector_value_list(q_points, cell_wind);
+ const std::vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
+ const std::vector<Point<dim> > &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<dh.n_dofs() ; ++i) {
+ // 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<dh.n_dofs() ; ++i) {
- local_matrix_row_i = 0;
+ local_matrix_row_i = 0;
- bool is_singular = false;
- unsigned int singular_index = numbers::invalid_unsigned_int;
+ bool is_singular = false;
+ unsigned int singular_index = numbers::invalid_unsigned_int;
- for(unsigned int j=0; j<fe.dofs_per_cell; ++j)
- if(local_dof_indices[j] == i) {
- singular_index = j;
- is_singular = true;
- break;
- }
-
- // We then perform the
- // integral. If the index $i$ is
- // not one of the local degrees
- // of freedom, we simply have to
- // add the single layer terms to
- // the right hand side, and the
- // double layer terms to the
- // matrix:
- if(is_singular == false) {
- for(unsigned int q=0; q<n_q_points; ++q) {
- normal_wind = 0;
- for(unsigned int d=0; d<dim; ++d)
- normal_wind += normals[q][d]*cell_wind[q](d);
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ if(local_dof_indices[j] == i) {
+ singular_index = j;
+ is_singular = true;
+ break;
+ }
+
+ // We then perform the
+ // integral. If the index $i$
+ // is not one of the local
+ // degrees of freedom, we
+ // simply have to add the
+ // single layer terms to the
+ // right hand side, and the
+ // double layer terms to the
+ // matrix:
+ if(is_singular == false) {
+ for(unsigned int q=0; q<n_q_points; ++q) {
+ normal_wind = 0;
+ for(unsigned int d=0; d<dim; ++d)
+ normal_wind += normals[q][d]*cell_wind[q](d);
- const Point<dim> R = q_points[q] - support_points[i];
+ const Point<dim> 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<fe.dofs_per_cell; ++j) {
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
- 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<dim-1>. 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 <em>know</em> 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<dim-1> *
- singular_quadrature
- = (dim == 2
- ?
- dynamic_cast<Quadrature<dim-1>*>(
- new QGaussLogR<1>(singular_quadrature_order,
- Point<1>((double)singular_index),
- 1./cell->measure()))
- :
- (dim == 3
- ?
- dynamic_cast<Quadrature<dim-1>*>(
- &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<dim-1>. 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 <em>know</em> 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<dim-1> *
+ singular_quadrature
+ = (dim == 2
+ ?
+ dynamic_cast<Quadrature<dim-1>*>(
+ new QGaussLogR<1>(singular_quadrature_order,
+ Point<1>((double)singular_index),
+ 1./cell->measure()))
+ :
+ (dim == 3
+ ?
+ dynamic_cast<Quadrature<dim-1>*>(
+ &sing_quadratures_3d[singular_index])
+ :
+ 0));
+ Assert(singular_quadrature, ExcInternalError());
- FEValues<dim-1,dim> fe_v_singular (fe, *singular_quadrature,
- update_jacobians |
- update_values |
- update_cell_normal_vectors |
- update_quadrature_points );
+ FEValues<dim-1,dim> 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<Vector<double> > singular_cell_wind( (*singular_quadrature).size(),
- Vector<double>(dim) );
+ std::vector<Vector<double> > singular_cell_wind( (*singular_quadrature).size(),
+ Vector<double>(dim) );
-
const std::vector<Point<dim> > &singular_normals = fe_v_singular.get_cell_normal_vectors();
-
const std::vector<Point<dim> > &singular_q_points = fe_v_singular.get_quadrature_points();
+ const std::vector<Point<dim> > &singular_normals = fe_v_singular.get_cell_normal_vectors();
+ const std::vector<Point<dim> > &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; q<singular_quadrature->size(); ++q) {
- const Point<dim> R = singular_q_points[q]- support_points[i];
- double normal_wind = 0;
- for(unsigned int d=0; d<dim; ++d)
- normal_wind += (singular_cell_wind[q](d)*
- singular_normals[q][d]);
+ for(unsigned int q=0; q<singular_quadrature->size(); ++q) {
+ const Point<dim> R = singular_q_points[q]- support_points[i];
+ double normal_wind = 0;
+ for(unsigned int d=0; d<dim; ++d)
+ normal_wind += (singular_cell_wind[q](d)*
+ singular_normals[q][d]);
- system_rhs(i) += ( LaplaceKernel::single_layer(R, is_singular) *
- normal_wind *
- fe_v_singular.JxW(q) );
+ system_rhs(i) += ( LaplaceKernel::single_layer(R, is_singular) *
+ normal_wind *
+ fe_v_singular.JxW(q) );
- for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
- local_matrix_row_i(j) += (( LaplaceKernel::double_layer(R, is_singular) *
- singular_normals[q]) *
- fe_v_singular.shape_value(j,q) *
- fe_v_singular.JxW(q) );
- }
- }
- if(dim==2) {
- delete singular_quadrature;
- }
- }
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
+ local_matrix_row_i(j) += (( LaplaceKernel::double_layer(R, is_singular) *
+ singular_normals[q]) *
+ fe_v_singular.shape_value(j,q) *
+ fe_v_singular.JxW(q) );
+ }
+ }
+ if(dim==2) {
+ delete singular_quadrature;
+ }
+ }
- // Finally, we need to add the contributions of the current cell
- // to the global matrix:
- for(unsigned int j=0; j<fe.dofs_per_cell; ++j)
- system_matrix.add(i,
- local_dof_indices[j],
- local_matrix_row_i(j));
- }
+ // Finally, we need to add the
+ // contributions of the current
+ // cell to the global matrix:
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ system_matrix.add(i,
+ local_dof_indices[j],
+ local_matrix_row_i(j));
}
+ }
- // The second part of the integral
- // operator is the term
- // $\alpha(\mathbf{x}_i)
- // \phi_j(\mathbf{x}_i)$. Since we use a
- // collocation scheme,
- // $\phi_j(\mathbf{x}_i)=\delta_{ij}$ and
- // the corresponding matrix is a diagonal
- // one with entries equal to
- // $\alpha(\mathbf{x}_i)$.
+ // The second part of the integral
+ // operator is the term
+ // $\alpha(\mathbf{x}_i)
+ // \phi_j(\mathbf{x}_i)$. Since we
+ // use a collocation scheme,
+ // $\phi_j(\mathbf{x}_i)=\delta_{ij}$
+ // and the corresponding matrix is
+ // a diagonal one with entries
+ // equal to $\alpha(\mathbf{x}_i)$.
- // One quick way to compute this diagonal matrix of the solid angles, is
- // to use the Neumann matrix itself. It is enough to multiply the matrix
- // with a vector of elements all equal to -1, to get the diagonal matrix
- // of the alpha angles, or solid angles (see the formula in the
- // introduction for this). The result is then added back onto the system
- // matrix object to yield the final form of the matrix:
- Vector<double> ones(dh.n_dofs());
- ones.add(-1.);
-
- system_matrix.vmult(alpha, ones);
- for(unsigned int i = 0; i<dh.n_dofs(); ++i)
- system_matrix.add(i,i,alpha(i));
+ // One quick way to compute this
+ // diagonal matrix of the solid
+ // angles, is to use the Neumann
+ // matrix itself. It is enough to
+ // multiply the matrix with a
+ // vector of elements all equal to
+ // -1, to get the diagonal matrix
+ // of the alpha angles, or solid
+ // angles (see the formula in the
+ // introduction for this). The
+ // result is then added back onto
+ // the system matrix object to
+ // yield the final form of the
+ // matrix:
+ Vector<double> ones(dh.n_dofs());
+ ones.add(-1.);
+
+ system_matrix.vmult(alpha, ones);
+ for(unsigned int i = 0; i<dh.n_dofs(); ++i)
+ system_matrix.add(i,i,alpha(i));
}
// @sect4{BEMProblem::solve_system}
- // The next function simply solves the linear
- // system. As described, we use the
- // SparseDirectUMFPACK direct solver to
- // compute the inverse of the matrix (in
- // reality it only produces an LU
- // decomposition) and then apply this inverse
- // to the right hand side to yield the
- // solution.
+ // The next function simply solves
+ // the linear system. As described,
+ // we use the SparseDirectUMFPACK
+ // direct solver to compute the
+ // inverse of the matrix (in reality
+ // it only produces an LU
+ // decomposition) and then apply this
+ // inverse to the right hand side to
+ // yield the solution.
//
// As mentioned in the introduction,
// the solution is only known up to a
// vector entry to normalize it.
template <int dim>
void BEMProblem<dim>::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 <int dim>
void BEMProblem<dim>::compute_errors(const unsigned int cycle) {
- Vector<float> 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<float> 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<double> 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<double> 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 <int dim>
void BEMProblem<dim>::compute_exterior_solution() {
- Triangulation<dim> external_tria;
- GridGenerator::hyper_cube(external_tria, -2, 2);
+ Triangulation<dim> external_tria;
+ GridGenerator::hyper_cube(external_tria, -2, 2);
- FE_Q<dim> external_fe(1);
- DoFHandler<dim> external_dh (external_tria);
- Vector<double> external_phi;
+ FE_Q<dim> external_fe(1);
+ DoFHandler<dim> external_dh (external_tria);
+ Vector<double> 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<dim-1,dim>::active_cell_iterator
- cell = dh.begin_active(),
- endc = dh.end();
+ typename DoFHandler<dim-1,dim>::active_cell_iterator
+ cell = dh.begin_active(),
+ endc = dh.end();
- FEValues<dim-1,dim> fe_v(fe, *quadrature,
- update_values |
- update_cell_normal_vectors |
- update_quadrature_points |
- update_JxW_values);
+ FEValues<dim-1,dim> 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<unsigned int> dofs(fe.dofs_per_cell);
+ std::vector<unsigned int> dofs(fe.dofs_per_cell);
- std::vector<double> local_phi(n_q_points);
- std::vector<double> normal_wind(n_q_points);
- std::vector<Vector<double> > local_wind(n_q_points, Vector<double>(dim) );
+ std::vector<double> local_phi(n_q_points);
+ std::vector<double> normal_wind(n_q_points);
+ std::vector<Vector<double> > local_wind(n_q_points, Vector<double>(dim) );
- typename DoFHandler<dim>::active_cell_iterator
- external_cell = external_dh.begin_active(),
- external_endc = external_dh.end();
+ typename DoFHandler<dim>::active_cell_iterator
+ external_cell = external_dh.begin_active(),
+ external_endc = external_dh.end();
-
std::vector<Point<dim> > external_support_points(external_dh.n_dofs());
- DoFTools::map_dofs_to_support_points<dim>( StaticMappingQ1<dim>::mapping,
- external_dh, external_support_points);
+ std::vector<Point<dim> > external_support_points(external_dh.n_dofs());
+ DoFTools::map_dofs_to_support_points<dim>( StaticMappingQ1<dim>::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<Point<dim> > &q_points = fe_v.get_quadrature_points();
-
const std::vector<Point<dim> > &normals = fe_v.get_cell_normal_vectors();
+ const std::vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
+ const std::vector<Point<dim> > &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<n_q_points; ++q){
- normal_wind[q] = 0;
- for(unsigned int d=0; d<dim; ++d)
- normal_wind[q] += normals[q][d]*local_wind[q](d);
- }
+ for(unsigned int q=0; q<n_q_points; ++q){
+ normal_wind[q] = 0;
+ for(unsigned int d=0; d<dim; ++d)
+ normal_wind[q] += normals[q][d]*local_wind[q](d);
+ }
- for(unsigned int i=0; i<external_dh.n_dofs(); ++i) {
+ for(unsigned int i=0; i<external_dh.n_dofs(); ++i) {
- for(unsigned int q=0; q<n_q_points; ++q) {
+ for(unsigned int q=0; q<n_q_points; ++q) {
- const Point<dim> R = q_points[q] - external_support_points[i];
+ const Point<dim> 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<dim> data_out;
+ DataOut<dim> 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 <int dim>
void BEMProblem<dim>::output_results(const unsigned int cycle) {
- DataOut<dim-1, DoFHandler<dim-1, dim> > dataout;
+ DataOut<dim-1, DoFHandler<dim-1, dim> > 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 <int dim>
void BEMProblem<dim>::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<n_cycles; ++cycle) {
- refine_and_resize();
- assemble_system();
- solve_system();
- compute_errors(cycle);
- output_results(cycle);
- }
+ for(unsigned int cycle=0; cycle<n_cycles; ++cycle) {
+ refine_and_resize();
+ assemble_system();
+ solve_system();
+ compute_errors(cycle);
+ output_results(cycle);
+ }
- if(extend_solution == true)
- compute_exterior_solution();
+ if(extend_solution == true)
+ compute_exterior_solution();
}
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
projects / dealii-svn.git / commitdiff