set Extend solution on the -2,2 box = true
set External refinement = 5
set Number of cycles = 4
+set Run 2d simulation = true
+set Run 3d simulation = true
+
+
+subsection Exact solution 2d
+ # Any constant used inside the function which is not a variable name.
+ set Function constants =
+
+ # Separate vector valued expressions by ';' as ',' is used internally by the
+ # function parser.
+ set Function expression = x+y # default: 0
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,t
+end
+
+
+subsection Exact solution 3d
+ # Any constant used inside the function which is not a variable name.
+ set Function constants =
+
+ # Separate vector valued expressions by ';' as ',' is used internally by the
+ # function parser.
+ set Function expression = x+y+z # default: 0
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,z,t
+end
subsection Quadrature rules
# Separate vector valued expressions by ';' as ',' is used internally by the
# function parser.
- set Function expression = 1; 1
+ set Function expression = 1; 1 # default: 0; 0
# The name of the variables as they will be used in the function, separated
# by ','.
# Separate vector valued expressions by ';' as ',' is used internally by the
# function parser.
- set Function expression = 1; 1; 1
+ set Function expression = 1; 1; 1 # default: 0; 0; 0
# The name of the variables as they will be used in the function, separated
# by ','.
When we run the program, the following is printed on screen:
@verbatim
-
+DEAL::
DEAL::Parsing parameter file parameters.prm
DEAL::for a 2 dimensional simulation.
-DEAL::Levels: 2, potential dofs: 20
-DEAL::Levels: 3, potential dofs: 40
-DEAL::Levels: 4, potential dofs: 80
-DEAL::Levels: 5, potential dofs: 160
-
+DEAL::Cycle 0:
+DEAL:: Number of active cells: 20
+DEAL:: Number of degrees of freedom: 20
+DEAL::Cycle 1:
+DEAL:: Number of active cells: 40
+DEAL:: Number of degrees of freedom: 40
+DEAL::Cycle 2:
+DEAL:: Number of active cells: 80
+DEAL:: Number of degrees of freedom: 80
+DEAL::Cycle 3:
+DEAL:: Number of active cells: 160
+DEAL:: Number of degrees of freedom: 160
+DEAL::
+cycle cells dofs L2(phi) Linfty(alpha)
+ 0 20 20 1.149e-07 - 5.000e-02 -
+ 1 40 40 4.330e-08 1.41 2.500e-02 1.00
+ 2 80 80 1.848e-08 1.23 1.250e-02 1.00
+ 3 160 160 8.798e-09 1.07 6.250e-03 1.00
+DEAL::
DEAL::Parsing parameter file parameters.prm
DEAL::for a 3 dimensional simulation.
-DEAL::Levels: 2, potential dofs: 26
-DEAL::Levels: 3, potential dofs: 98
-DEAL::Levels: 4, potential dofs: 386
-DEAL::Levels: 5, potential dofs: 1538
+DEAL::Cycle 0:
+DEAL:: Number of active cells: 24
+DEAL:: Number of degrees of freedom: 26
+DEAL::Cycle 1:
+DEAL:: Number of active cells: 96
+DEAL:: Number of degrees of freedom: 98
+DEAL::Cycle 2:
+DEAL:: Number of active cells: 384
+DEAL:: Number of degrees of freedom: 386
+DEAL::Cycle 3:
+DEAL:: Number of active cells: 1536
+DEAL:: Number of degrees of freedom: 1538
+DEAL::
+cycle cells dofs L2(phi) Linfty(alpha)
+ 0 24 26 3.415e-06 - 2.327e-01 -
+ 1 96 98 7.248e-07 2.24 1.239e-01 0.91
+ 2 384 386 1.512e-07 2.26 6.319e-02 0.97
+ 3 1536 1538 6.576e-08 1.20 3.176e-02 0.99
+
@endverbatim
+As we can see from the convergence table in 2d, if we choose
+quadrature formulas which are accurate enough, then the error we
+obtain for $\alpha(\mathbf{x})$ should be exactly the inverse of the
+number of elements. The approximation of the circle with N segments of
+equal size generates a regular polygon with N faces, whose angles are
+exactly $\pi-\frac {2\pi}{N}$, therefore the error we commit should be
+exactly $\frac 12 - (\frac 12 -\frac 1N) = \frac 1N$. In fact this is
+a very good indicator that we are performing the singular integrals in
+an appropiate manner.
+
+The error in the approximation of the potential $\phi$ is largely due
+to approximation of the domain. A much better approximation could be
+obtained by using higher order mappings.
+
The result of running these computations is a bunch of output files that we
can pass to our visualization program of choice.
together with a countur plot:
@image html step-34_3d.png
-
-The exact solution for this model problem should be $x+y$ or
-$x+y+z$ on the surface, which is precisely what we obtain.
#include <lac/matrix_lib.h>
#include <numerics/data_out.h>
+#include <numerics/vectors.h>
#include <cmath>
#include <iostream>
// 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. 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
DoFHandler<dim> external_dh;
Vector<double> external_phi;
+ // 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 one that we fill through a
// parameter file. The new objects that we use in this example
// are the ParsedFunction object and the QuadratureSelector
// 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;
+
SmartPointer<Quadrature<dim-1> > quadrature_pointer;
unsigned int singular_quadrature_order;
+
unsigned int n_cycles;
unsigned int external_refinement;
+
bool run_in_this_dimension;
bool extend_solution;
};
};
-
+// The constructor initializes the variuous object in the same way of
+// finite element problems. 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 this 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
+// ParsedFunction, we need to specify the number of components
+// explicitly, since the function
+// ParsedFunction<dim>::declare_parameters is static, and has no
+// knowledge of the number of components.
template <int dim>
BEMProblem<dim>::BEMProblem() :
fe(1),
prm.declare_entry("Singular quadrature order", "5", Patterns::Integer());
prm.leave_subsection();
- // For both two and three dimensions, we set the input data to be
- // such that the solution is $x+y+c$ or $x+y+z+c$.
+ // For both two and three dimensions, we set the default input
+ // data to be such that the solution is $x+y+c$ or $x+y+z+c$.
+ //
+ // The use of the ParsedFunction object is pretty straight
+ // forward. The 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
+ // correspending 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 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");
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 3d");
+ Functions::ParsedFunction<3>::declare_parameters(prm);
+ prm.set("Function expression","x+y+z");
+ prm.leave_subsection();
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");
+
+ // 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.
+ //
+ // This is another example of how to use parameter files in
+ // dimensio independent programming.
run_in_this_dimension = prm.get_bool("Run " +
Utilities::int_to_string(dim) +
"d simulation");
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();
+
quadrature_pointer = &quadrature;
}
// 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
- // the orientation.
+ // the orientation. In our case, the normals to the mesh are
+ // external to both the circle and the sphere.
//
// The other detail that is required for appropriate refinement of
// the boundary element mesh, is an accurate description of the
dh.distribute_dofs(fe);
- const unsigned int ndofs = dh.n_dofs();
-
- deallog << "Levels: " << tria.n_levels()
- << ", potential dofs: " << ndofs << endl;
+ const unsigned int n_dofs = dh.n_dofs();
// The matrix is a full matrix. Notwithstanding this fact, the
// SparseMatrix class coupled with the SparseDirectUMFPACK solver
// are still faster than Lapack solvers. The drawback is that we
// need to assemble a full SparsityPattern.
system_matrix.clear();
- sparsity.reinit(ndofs, ndofs, ndofs);
- for(unsigned int i=0; i<ndofs;++i)
- for(unsigned int j=0; j<ndofs; ++j)
+ 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(ndofs);
- phi.reinit(ndofs);
- alpha.reinit(ndofs);
+ system_rhs.reinit(n_dofs);
+ phi.reinit(n_dofs);
+ alpha.reinit(n_dofs);
}
template <int dim>
for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
local_matrix_row_i(j) += (( kernel.double_layer(R, is_singular) *
- singular_normals[q]) *
- fe_v_singular.shape_value(j,q) *
- fe_v_singular.JxW(q) );
+ singular_normals[q]) *
+ fe_v_singular.shape_value(j,q) *
+ fe_v_singular.JxW(q) );
}
}
if(dim==2) {
// only known up to a constant potential. We solve this issue by
// subtracting the mean value of the vector from each vector
// entry.
- double mean_value = phi.mean_value();
- for(unsigned int i=0; i<phi.size(); ++i)
- phi(i) -= mean_value;
+ phi.add(-phi.mean_value());
}
+
+template <int dim>
+void BEMProblem<dim>::compute_errors(const unsigned int cycle) {
+ // 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.
+
+ 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 linfty_norm method of Vector<double>, since on each node,
+ // the value should be $\frac 12$.
+ 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();
+
+ 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);
+}
+
// We 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.
external_phi(i) += ( ( kernel.single_layer(R) *
normal_wind[q] +
//
- (- kernel.double_layer(R) *
+ (kernel.double_layer(R) *
normals[q] ) *
local_phi[q] ) *
fe_v.JxW(q) );
std::ofstream file(filename.c_str());
dataout.write_vtk(file);
+
+ 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);
+
+ if(cycle == n_cycles-1) {
+ 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);
+ }
}
template <int dim>
refine_and_resize();
assemble_system();
solve_system();
+ compute_errors(cycle);
output_results(cycle);
}
if(extend_solution == true)
interpolate();
} else {
- deallog << "Run in dimension " << dim
+ deallog << "Run in dimension " << dim
<< " explicitly disabled in parameter file. "
<< std::endl;
}
projects / dealii-svn.git / commitdiff