From: heltai Date: Fri, 13 Mar 2009 16:10:25 +0000 (+0000) Subject: Added better comments to step-34, created interpolate function. It seems to work... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0839b55e8e2186c66cfedfc805236a6c42accb1c;p=dealii-svn.git Added better comments to step-34, created interpolate function. It seems to work in 2d. Not completely tested in 3d yet. It does not yet work for higher order approximations. Only for dgp(0). git-svn-id: https://svn.dealii.org/trunk@18484 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-34/Makefile b/deal.II/examples/step-34/Makefile index 3dc5c73fbf..1089cac0b3 100644 --- a/deal.II/examples/step-34/Makefile +++ b/deal.II/examples/step-34/Makefile @@ -60,10 +60,12 @@ include $D/common/Make.global_options # program for other dimensions, or when using third party libraries libs.g = $(lib-deal2-3d.g) \ $(lib-deal2-2d.g) \ + $(lib-deal2-1d.g) \ $(lib-lac.g) \ $(lib-base.g) libs.o = $(lib-deal2-3d.o) \ $(lib-deal2-2d.o) \ + $(lib-deal2-1d.o) \ $(lib-lac.o) \ $(lib-base.o) diff --git a/deal.II/examples/step-34/coarse_circle.inp b/deal.II/examples/step-34/coarse_circle.inp new file mode 100644 index 0000000000..dda97cde9e --- /dev/null +++ b/deal.II/examples/step-34/coarse_circle.inp @@ -0,0 +1,21 @@ +10 10 0 0 0 +1 0 1.0 0 +2 0.587785252292 0.809016994375 0 +3 0.951056516295 0.309016994375 0 +4 0.951056516295 -0.309016994375 0 +5 0.587785252292 -0.809016994375 0 +6 0 -1.0 0 +7 -0.587785252292 -0.809016994375 0 +8 -0.951056516295 -0.309016994375 0 +9 -0.951056516295 0.309016994375 0 +10 -0.587785252292 0.809016994375 0 +1 1 line 1 2 +2 1 line 2 3 +3 1 line 3 4 +4 1 line 4 5 +5 1 line 5 6 +6 1 line 6 7 +7 1 line 7 8 +8 1 line 8 9 +9 1 line 9 10 +10 1 line 10 1 diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index d91513a956..9bd40dc837 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -1,6 +1,7 @@ # Listing of Parameters # --------------------- -set Number of cycles = 3 +set Number of cycles = 3 +set External refinement = 5 subsection Inner quadrature rule @@ -14,8 +15,20 @@ subsection Outer quadrature rule set Quadrature type = midpoint end +subsection Wind function 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 = 1; 0 + + # The name of the variables as they will be used in the function, separated + # by ','. + set Variable names = x,y,t +end -subsection Wind function +subsection Wind function 3d # Any constant used inside the function which is not a variable name. set Function constants = diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index eb87d30313..5fdf13e53f 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -13,72 +13,213 @@ // //---------------------------- step-34.cc --------------------------- - -// - -#include #include - +#include #include #include #include #include -#include -#include -#include #include +#include + +#include +#include +#include +#include +#include +#include +#include + #include #include #include #include #include -#include -#include -#include -#include -#include -#include -#include + +#include +#include + #include #include #include #include #include #include + #include -#include #include #include +#include #include #include using namespace std; using namespace dealii; + template -class LaplaceKernelIntegration +class LaplaceKernelIntegration; + + +template +class BEMProblem { public: + BEMProblem(const unsigned int degree = 0); + ~BEMProblem(); + + // The structure of a boundary element method code is very similar + // to the structure of a finite element code. By now you should be + // familiar with reading paramaters from an external file, and + // with the splitting of the different tasks into different + // modules. The same applyes to boundary element methods, and we + // won't comment too much on them, except on the differences. - LaplaceKernelIntegration(FiniteElement &fe); - ~LaplaceKernelIntegration(); - + void read_parameters(std::string filename); + void run(); + void read_domain(); + + 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 two nested loops on cells + // instead of the usual one we have in finite element method. + // + // The reason for this is that while the basis functions have a + // compact support, their convolution with the fundamental + // solution of the laplace equation is global, and needs to be + // integrated against all other basis functions. + // + // The practical consequence is that we have two sets of + // quadrature formulas, finite element values and temporary + // elements, one for the inner integration and one for the outer + // integration. We allow for different quadrature rules to be used + // in the two integrations to preserve generality and to allow, + // for example, the use of collocation method (by specifying midpoint + // quadrature formula on the outer integration). + void assemble_system(); + + // The only difference in the solution of the system is that the + // matrix is a LAPACKFullMatrix, which requires a different + // treatment with respect to what we saw in most of the other + // examples. Besides from this detail, things proceeds pretty much + // in the same way as usual. + void solve_system(); + + // 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 interpolate() 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. + void interpolate(); + + void output_results(unsigned int cycle); + +private: + // 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 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_DGP 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. Just for the purpose of illustrating the use of the + // LAPACK classes, we opt for LU decomposition of the final + // system. Note that this will be very inefficient when the number + // of dofs grows, since it is of order $n^3$. + + SmartPointer > system_matrix; + Vector system_rhs; + Vector phi; + + // 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. + + Triangulation external_tria; + FE_Q external_fe; + DoFHandler external_dh; + Vector external_phi; + + // 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 object. + // + // The 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). + // + // The QuadratureSelector class 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 quadrature formulas for the inner as well as the + // outer integration in the calculation of the boundary element + // matrix. + // + // Notice that selecting the midpoint rule as the outer + // integration formula on uniformly refined meshes is equivalent + // (up to a scaling factor) to solving the boundary element method + // via collocation instead of Galerkin technique. + Functions::ParsedFunction wind; + SmartPointer > outer_quadrature_pointer; + SmartPointer > inner_quadrature_pointer; + unsigned int n_cycles; + unsigned int external_refinement; +}; + + + +template +class LaplaceKernelIntegration +{ +public: + + LaplaceKernelIntegration(const FiniteElement &fe); + ~LaplaceKernelIntegration(); + // This functions computes the integral of the single and double // layer potentials on the cell given as a parameter, at the // quadrature points @p q. In practice this function produces the objects // // \f[ - // \text{dst}_{ik0} := \int_{\text{cell}} G(y - \text[q]_k) \phi_i dy + // \text{dst}_{ik0} := \int_{\text{cell}} G(y - \text[q]_k) rhs(y) dy // \f] // // and // // \f[ // \text{dst}_{ik1} := \int_{\text{cell}} \frac{\partial - // G}{\partial \textbf{n}} (y - \text[q]_k) \phi_i dy + // G}{\partial \textbf{n}} (y - \text[q]_k) \phi_i(y) dy // \f] void compute_SD_integral_on_cell(vector > > &dst, typename DoFHandler::active_cell_iterator &cell, @@ -87,7 +228,7 @@ public: // The following two functions are the actual calculations of the // single and double layer potential kernels, with a minus sign in - // front of them. They are well defined only if the vector $r = + // front of them. They are well defined only if the vector $R = // x-y$ is different from zero. double nS(const Point &R); Point nD(const Point &R); @@ -114,83 +255,19 @@ private: AssertThrow(false, ExcImpossibleInDim()); return 0; }; - + + SmartPointer > fe; SmartPointer > fe_values; }; -template -class BEMProblem -{ -public: - BEMProblem(); - ~BEMProblem(); - - // Read parameters. - void read_parameters(std::string filename); - - // Starts the Boundary Element Method Computation. - void run(); - - // Initialize mesh and vector space. - void read_domain(); - - // Refine and resize all vectors for the active step. - void refine_and_resize(); - - // Assemble the two system matrices as well as the system right - // hands side. - void assemble_system(); - // Solve the system. - void solve_system(); - - // Output results for the given cycle. - void output_results(unsigned int cycle); - -private: - // The boundary element method triangulation. - Triangulation tria; - - // The finite element spaces for the potential and the velocity. - FE_DGP fe; - FESystem fev; - - // Finite element space used to smoothen the potential solution - // (from piecewise constant to continuous piecewise quadratic) - FE_Q fe_q; - - // And the relevant degrees of freedom. - DoFHandler dh; - DoFHandler dhv; - DoFHandler dhq; - - // The system matrix. This is I-C. Since the LAPACKFullMatrix does - // not have a reinit method, we need to work around this a little. - SmartPointer > system_matrix; - - // The right hand side, the potential and its smoothed version - Vector system_rhs; - Vector phi; - Vector smooth_phi; - - // These are the parameters that we read in from a parameter file. - // In particular we define the wind function and the outer - // quadrature. We use a parsed function, for its ease of - // definition, and the quadrature formula - Functions::ParsedFunction wind; - SmartPointer > outer_quadrature_pointer; - SmartPointer > inner_quadrature_pointer; - unsigned int n_cycles; -}; template -BEMProblem::BEMProblem() : - fe(0), - fev(FE_DGP(0), dim), - fe_q(FE_Q(2)), +BEMProblem::BEMProblem(const unsigned int degree) : + fe(degree), dh(tria), - dhv(tria), - dhq(tria), + external_fe(1), + external_dh(external_tria), wind(dim) {} @@ -207,6 +284,7 @@ void BEMProblem::read_parameters(std::string filename) { ParameterHandler prm; prm.declare_entry("Number of cycles", "4", Patterns::Integer()); + prm.declare_entry("External refinement", "5", Patterns::Integer()); prm.enter_subsection("Outer quadrature rule"); prm.declare_entry("Quadrature type", "midpoint", @@ -216,18 +294,23 @@ void BEMProblem::read_parameters(std::string filename) { prm.enter_subsection("Inner quadrature rule"); - prm.declare_entry("Quadrature type", "midpoint", + prm.declare_entry("Quadrature type", "gauss", Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names())); - prm.declare_entry("Quadrature order", "0", Patterns::Integer()); + prm.declare_entry("Quadrature order", "2", Patterns::Integer()); prm.leave_subsection(); - prm.enter_subsection("Wind function"); - Functions::ParsedFunction::declare_parameters(prm, dim); + prm.enter_subsection("Wind function 2d"); + Functions::ParsedFunction<2>::declare_parameters(prm, 2); + prm.leave_subsection(); + + prm.enter_subsection("Wind function 3d"); + Functions::ParsedFunction<3>::declare_parameters(prm, 3); prm.leave_subsection(); prm.read_input(filename); n_cycles = prm.get_integer("Number of cycles"); + external_refinement = prm.get_integer("External refinement"); prm.enter_subsection("Outer quadrature rule"); static QuadratureSelector outer_quadrature @@ -242,7 +325,8 @@ void BEMProblem::read_parameters(std::string filename) { prm.get_integer("Quadrature order")); prm.leave_subsection(); - prm.enter_subsection("Wind function"); + prm.enter_subsection(std::string("Wind function ")+ + Utilities::int_to_string(dim)+std::string("d")); wind.parse_parameters(prm); prm.leave_subsection(); @@ -282,7 +366,8 @@ Point LaplaceKernelIntegration::nD(const Point &R) { template <> -LaplaceKernelIntegration<3>::LaplaceKernelIntegration(FiniteElement<2,3> &fe) +LaplaceKernelIntegration<3>::LaplaceKernelIntegration(const FiniteElement<2,3> &fe) : + fe(&fe) { // In order to perform the two dimensional singular integration on // the given cell, we use standard formulas derived by Morino and @@ -306,6 +391,15 @@ LaplaceKernelIntegration<3>::LaplaceKernelIntegration(FiniteElement<2,3> &fe) update_quadrature_points ); } + +// The one dimensional singular integration can be calculated +// exploiting QGaussLogR quadrature formula. The quadrature formula +// is constructed in each step, so the constructor is empty. +template <> +LaplaceKernelIntegration<2>::LaplaceKernelIntegration(const FiniteElement<1,2> &fe) : + fe(&fe) +{} + template LaplaceKernelIntegration::~LaplaceKernelIntegration() { // We delete the pointer. Since this was created via the new @@ -313,9 +407,11 @@ LaplaceKernelIntegration::~LaplaceKernelIntegration() { // not take smart pointers, which implies we need to first remove // detach the smart pointer from the fe_values object, and then // delete it. - FEValues *fp = fe_values; - fe_values = 0; - delete fp; + if(fe_values) { + FEValues *fp = fe_values; + fe_values = 0; + delete fp; + } } @@ -426,13 +522,42 @@ LaplaceKernelIntegration<3>::compute_SD_integral_on_cell(vector void BEMProblem::read_domain() { - // Center of the ball. It is the origin by default. + + // A boundary element method triangulation is basically the same + // as a (dim-1) 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 + // the orientation. + // + // 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, and here the principle and usage is the same, except + // that the Boundary description class takes an additional + // template parameter that specifies the embedding space + // dimension. + Point p; static HyperBallBoundary boundary(p,1.); - // Read the sphere from GridIn gi; gi.attach_triangulation (tria); if(dim == 3) { @@ -452,14 +577,11 @@ void BEMProblem::refine_and_resize() { tria.refine_global(1); dh.distribute_dofs(fe); - dhv.distribute_dofs(fev); const unsigned int ndofs = dh.n_dofs(); - const unsigned int nvdofs = dhv.n_dofs(); deallog << "Levels: " << tria.n_levels() - << ", potential dofs: " << ndofs - << ", velocity dofs: " << nvdofs << endl; + << ", potential dofs: " << ndofs << endl; if(system_matrix) { LAPACKFullMatrix * p = system_matrix; @@ -589,27 +711,74 @@ void BEMProblem::assemble_system() { // layer potential are singular, and they require a // special treatment, as explained in the // introduction. - - kernel.compute_SD_integral_on_cell(single_double_layer_potentials, - cellj, q_points_outer, wind); - - for(unsigned int i=0; i singular_quad(inner_quadrature.size(), + outer_quadrature.point(q_outer), + 1./cellj->measure()); + FEValues<1,2> fe_v_singular(fe, singular_quad, + update_jacobians | + update_cell_normal_vectors | + update_quadrature_points ); + fe_v_singular.reinit(cellj); - for(unsigned int j=0; j > singular_cell_wind(singular_quad.size(), + Vector(dim) ); + + const vector > &singular_normals = fe_v_singular.get_cell_normal_vectors(); + const vector > &singular_q_points = fe_v_singular.get_quadrature_points(); + + wind.vector_value_list(singular_q_points, singular_cell_wind); + + for(unsigned int i=0; i::solve_system() { } +// 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. +template +void BEMProblem::interpolate() { + // Generate the mesh, refine it and distribute dofs on it. + GridGenerator::hyper_cube(external_tria, -2, 2); + 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(); + + + Quadrature &quadrature = *inner_quadrature_pointer; + + 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; + + vector dofs(fe.dofs_per_cell); + + vector local_phi(n_q_points); + vector > local_wind(n_q_points, Vector(dim) ); + double normal_wind; + + LaplaceKernelIntegration kernel(fe); + Point R; + + + typename DoFHandler::active_cell_iterator + external_cell = external_dh.begin_active(), + external_endc = external_dh.end(); + + vector external_dofs(external_fe.dofs_per_cell); + vector dof_is_treated(external_dh.n_dofs(), false); + + + for(; external_cell != external_endc; ++external_cell) { + external_cell->get_dof_indices(external_dofs); + + for(unsigned int i=0; i > &q_points = fe_v.get_quadrature_points(); + const vector > &normals = fe_v.get_cell_normal_vectors(); + + cell->get_dof_indices(dofs); + fe_v.get_function_values(phi, local_phi); + + wind.vector_value_list(q_points, local_wind); + + for(unsigned int q=0; qvertex(i) - q_points[q]; + + external_phi(external_dofs[i]) += ( ( - kernel.nS(R) * + normal_wind - + // + ( kernel.nD(R) * + normals[q] ) * + local_phi[q] ) * + fe_v.JxW(q) ); + } + } + } + } + DataOut > dataout; + + dataout.attach_dof_handler(external_dh); + dataout.add_data_vector(external_phi, "external_phi"); + dataout.build_patches(); + + std::string filename = Utilities::int_to_string(dim) + "d_external.vtk"; + std::ofstream file(filename.c_str()); + dataout.write_vtk(file); +} + + template void BEMProblem::output_results(unsigned int cycle) { @@ -641,9 +905,11 @@ void BEMProblem::output_results(unsigned int cycle) { dataout.add_data_vector(phi, "phi"); dataout.build_patches(); - char fname[100]; - sprintf(fname, "test_%02d.vtk", cycle); - std::ofstream file(fname); + 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); } @@ -660,6 +926,8 @@ void BEMProblem::run() { solve_system(); output_results(cycle); } + + interpolate(); } @@ -668,9 +936,11 @@ int main () try { deallog.depth_console (3); - - BEMProblem<3> laplace_problem; - laplace_problem.run(); + BEMProblem<2> laplace_problem_2d; + // BEMProblem<3> laplace_problem_3d; + + laplace_problem_2d.run(); + // laplace_problem_3d.run(); } catch (std::exception &exc) {