From 58816da7d80932418aac96397843b2264cd2842f Mon Sep 17 00:00:00 2001 From: heltai Date: Mon, 23 Feb 2009 08:32:15 +0000 Subject: [PATCH] Added preliminary version of step-34. BEM. git-svn-id: https://svn.dealii.org/trunk@18413 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-34/Makefile | 156 +++++++ deal.II/examples/step-34/coarse_sphere.inp | 15 + deal.II/examples/step-34/doc/intro.dox | 9 + deal.II/examples/step-34/doc/results.dox | 3 + deal.II/examples/step-34/step-34.cc | 450 +++++++++++++++++++++ 5 files changed, 633 insertions(+) create mode 100644 deal.II/examples/step-34/Makefile create mode 100644 deal.II/examples/step-34/coarse_sphere.inp create mode 100644 deal.II/examples/step-34/doc/intro.dox create mode 100644 deal.II/examples/step-34/doc/results.dox create mode 100644 deal.II/examples/step-34/step-34.cc diff --git a/deal.II/examples/step-34/Makefile b/deal.II/examples/step-34/Makefile new file mode 100644 index 0000000000..96fe808f62 --- /dev/null +++ b/deal.II/examples/step-34/Makefile @@ -0,0 +1,156 @@ +# $Id: Makefile 11909 2005-12-21 13:30:36Z guido $ + + +# For the small projects Makefile, you basically need to fill in only +# four fields. +# +# The first is the name of the application. It is assumed that the +# application name is the same as the base file name of the single C++ +# file from which the application is generated. +target = $(basename $(shell echo step-*.cc)) + +# The second field determines whether you want to run your program in +# debug or optimized mode. The latter is significantly faster, but no +# run-time checking of parameters and internal states is performed, so +# you should set this value to `on' while you develop your program, +# and to `off' when running production computations. +debug-mode = on + + +# As third field, we need to give the path to the top-level deal.II +# directory. You need to adjust this to your needs. Since this path is +# probably the most often needed one in the Makefile internals, it is +# designated by a single-character variable, since that can be +# reference using $D only, i.e. without the parentheses that are +# required for most other parameters, as e.g. in $(target). +D = ../../ + + +# The last field specifies the names of data and other files that +# shall be deleted when calling `make clean'. Object and backup files, +# executables and the like are removed anyway. Here, we give a list of +# files in the various output formats that deal.II supports. +clean-up-files = *gmv *gnuplot *gpl *eps *pov + + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# The next statement tell the `make' program where to find the +# deal.II top level directory and to include the file with the global +# settings +include $D/common/Make.global_options + + +################################################################ +# Since the whole project consists of only one file, we need not +# consider difficult dependencies. We only have to declare the +# libraries which we want to link to the object file, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2-3d.g) \ + $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-3d.o) \ + $(lib-deal2-2d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# We now use the variable defined above which switch between debug and +# optimized mode to select the set of libraries to link with. Included +# in the list of libraries is the name of the object file which we +# will produce from the single C++ file. Note that by default we use +# the extension .g.o for object files compiled in debug mode and .o for +# object files in optimized mode (or whatever the local default on your +# system is instead of .o). +ifeq ($(debug-mode),on) + libraries = $(target).g.$(OBJEXT) $(libs.g) +else + libraries = $(target).$(OBJEXT) $(libs.o) +endif + + +# Now comes the first production rule: how to link the single object +# file produced from the single C++ file into the executable. Since +# this is the first rule in the Makefile, it is the one `make' selects +# if you call it without arguments. +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) + + +# To make running the application somewhat independent of the actual +# program name, we usually declare a rule `run' which simply runs the +# program. You can then run it by typing `make run'. This is also +# useful if you want to call the executable with arguments which do +# not change frequently. You may then want to add them to the +# following rule: +run: $(target) + @echo ============================ Running $< + @./$(target)$(EXEEXT) + + +# As a last rule to the `make' program, we define what to do when +# cleaning up a directory. This usually involves deleting object files +# and other automatically created files such as the executable itself, +# backup files, and data files. Since the latter are not usually quite +# diverse, you needed to declare them at the top of this file. +clean: + -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) + + +# Since we have not yet stated how to make an object file from a C++ +# file, we should do so now. Since the many flags passed to the +# compiler are usually not of much interest, we suppress the actual +# command line using the `at' sign in the first column of the rules +# and write the string indicating what we do instead. +./%.g.$(OBJEXT) : + @echo ==============debug========= $( $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then : else rm $@ ; fi + + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + diff --git a/deal.II/examples/step-34/coarse_sphere.inp b/deal.II/examples/step-34/coarse_sphere.inp new file mode 100644 index 0000000000..5ea249c4dc --- /dev/null +++ b/deal.II/examples/step-34/coarse_sphere.inp @@ -0,0 +1,15 @@ +8 6 0 0 0 +1 -0.577350269 -0.577350269 -0.577350269 +2 0.577350269 -0.577350269 -0.577350269 +3 -0.577350269 0.577350269 -0.577350269 +4 0.577350269 0.577350269 -0.577350269 +5 -0.577350269 -0.577350269 0.577350269 +6 0.577350269 -0.577350269 0.577350269 +7 -0.577350269 0.577350269 0.577350269 +8 0.577350269 0.577350269 0.577350269 +1 1 quad 3 4 2 1 +2 1 quad 5 6 8 7 +3 1 quad 1 2 6 5 +4 1 quad 3 7 8 4 +5 1 quad 5 7 3 1 +6 1 quad 2 4 8 6 diff --git a/deal.II/examples/step-34/doc/intro.dox b/deal.II/examples/step-34/doc/intro.dox new file mode 100644 index 0000000000..0d85ffb09c --- /dev/null +++ b/deal.II/examples/step-34/doc/intro.dox @@ -0,0 +1,9 @@ + +

Introduction

+ +

The problem

+ +

Galerkin BEM

+ +

What the program does

+ diff --git a/deal.II/examples/step-34/doc/results.dox b/deal.II/examples/step-34/doc/results.dox new file mode 100644 index 0000000000..a1e1e35249 --- /dev/null +++ b/deal.II/examples/step-34/doc/results.dox @@ -0,0 +1,3 @@ +

Results

+ +

Simple Example

diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc new file mode 100644 index 0000000000..abb217697d --- /dev/null +++ b/deal.II/examples/step-34/step-34.cc @@ -0,0 +1,450 @@ +//---------------------------- bem_integration.cc --------------------------- +// $Id: testsuite.html 13373 2006-07-13 13:12:08Z manigrasso $ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +// Authors: Luca Heltai, Cataldo Manigrasso +// +//---------------------------- bem_integration.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 + +using namespace std; +using namespace dealii; + +template +class LaplaceKernelIntegration +{ +public: + + LaplaceKernelIntegration(); + ~LaplaceKernelIntegration(); + + void run(); + + void compute_SD_integral_on_cell(vector &dst, + typename DoFHandler::active_cell_iterator &cell, + const Point &point); + +private: + double term_S(const Point<3> &r, + const Point<3> &a1, + const Point<3> &a2, + const Point<3> &n, + const double &rn_c); + + double term_D(const Point<3> &r, + const Point<3> &a1, + const Point<3> &a2); + + SmartPointer > fe_values; +}; + +template +class BEMProblem +{ +public: + BEMProblem(); + ~BEMProblem(); + + /** 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 space for the potential. */ + FE_DGP<2,3> fe; + + /** The finite element space for the velocity. */ + FESystem<2,3> fev; + + /** The potential degrees of freedom. */ + DoFHandler dh; + + /** The velocity degrees of freedom. */ + DoFHandler dhv; + + /** 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; + + /** Single layer potential matrix. */ + FullMatrix single_layer_matrix; + + /** Normal component of the wind field. */ + Vector Vn; + + /** System rhs. */ + Vector rhs; + + /** Potential. */ + Vector phi; + + /** Something else. */ + Vector vs; +}; + +template +BEMProblem::BEMProblem() : + fe(0), + fev(FE_DGP(0), dim), + dh(tria), + dhv(tria) +{} + + +template +BEMProblem::~BEMProblem() { + LAPACKFullMatrix * p = system_matrix; + system_matrix = 0; + delete p; +} + + +template <> +LaplaceKernelIntegration<2>::LaplaceKernelIntegration() +{ + static FE_DGP<2,3> fe(0); + vector > qps(5); + qps[0] = Point<2>(0,0); + qps[1] = Point<2>(0,1); + qps[2] = Point<2>(1,0); + qps[3] = Point<2>(1,1); + qps[4] = Point<2>(.5,.5); + vector ws(5,1.); + static Quadrature<2> quadrature(qps, ws); + fe_values = new FEValues<2,3>(fe,quadrature, + update_values | + update_jacobians | + update_quadrature_points ); +} + +template +LaplaceKernelIntegration::~LaplaceKernelIntegration() { + FEValues *fp = fe_values; + fe_values = 0; + delete fp; +} + + +template <> +void +LaplaceKernelIntegration<2>::compute_SD_integral_on_cell(vector &dst, + DoFHandler<2,3>::active_cell_iterator &cell, + const Point<3> &point) +{ + Assert(dst.size() == 2, + ExcDimensionMismatch(dst.size(), 2)); + fe_values->reinit(cell); + vector > jacobians = fe_values->get_jacobians(); + vector > quad_points = fe_values->get_quadrature_points(); + Point<3> r,a1,a2,n,r_c,n_c; + r_c = point-cell->center(); + n_c = jacobians[4][2]; + double rn_c = r_c*n_c; + vector i_S(4); + vector i_D(4); + for (unsigned int q_point=0; q_point < 4; ++q_point) + { + r = point-quad_points[q_point]; + a1 = jacobians[q_point][0]; + a2 = jacobians[q_point][1]; + n = jacobians[q_point][2]; + i_S[q_point]=term_S(r,a1,a2,n,rn_c); + i_D[q_point]=term_D(r,a1,a2); + } + dst[0] = (i_S[3]-i_S[1]-i_S[2]+i_S[0]); + dst[1] = (i_D[3]-i_D[1]-i_D[2]+i_D[0]); + +} + +template +double +LaplaceKernelIntegration::term_S (const Point<3> &r, + const Point<3> &a1, + const Point<3> &a2, + const Point<3> &n, + const double &rn_c) +{ + Point<3> ra1, ra2, a12; + + cross_product(ra1,r,a1); + cross_product(ra2,r,a2); + cross_product(a12,a1,a2); + + double integral = + -1./2./numbers::PI + *( + - ra1*n/a1.norm() * asinh( r*a1/ra1.norm() ) + + ra2*n/a2.norm() * asinh( r*a2/ra2.norm() ) + + rn_c * atan( ra1*ra2 / (r.norm()* (r*(a12)))) + ); + + return integral; + +} + +template +double +LaplaceKernelIntegration::term_D (const Point<3> &r, + const Point<3> &a1, + const Point<3> &a2) +{ + Point<3> ra1, ra2, a12; + + cross_product(ra1,r,a1); + cross_product(ra2,r,a2); + cross_product(a12,a1,a2); + + double integral = -1./2./numbers::PI + *atan( ra1*ra2 / (r.norm()* (r*(a12)))); + + return integral; + +} + +template +void BEMProblem::read_domain() { + // Center of the ball. It is the origin by default. + Point p; + static HyperBallBoundary boundary(p,1.); + + // Read the sphere from + GridIn gi; + gi.attach_triangulation (tria); + if(dim == 3) { + std::ifstream in ("coarse_sphere.inp"); + gi.read_ucd (in); + } else if(dim == 2) { + std::ifstream in ("coarse_circle.inp"); + gi.read_ucd (in); + } + tria.set_boundary(1, boundary); +} + + + +template +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; + + if(system_matrix) { + LAPACKFullMatrix * p = system_matrix; + system_matrix = 0; + delete p; + } + + system_matrix = new LAPACKFullMatrix(ndofs, ndofs); + + Vn.reinit(ndofs); + rhs.reinit(ndofs); + phi.reinit(ndofs); + vs.reinit(nvdofs); +} + +template +void BEMProblem::assemble_system() { + + Point wind; + wind[0] = 1.; + + typename DoFHandler::active_cell_iterator + celli = dh.begin_active(), + cellj = dh.begin_active(), + cellv = dhv.begin_active(), + endc = dh.end(); + + QMidpoint midpoint; + FEValues fe_mid(fe, midpoint, + update_values | + update_cell_normal_vectors | + update_quadrature_points); + + vector dofsi(fe_mid.n_quadrature_points); + vector dofsj(fe_mid.n_quadrature_points); + vector dofsv(dim); + + // Temporary matrix + FullMatrix _B(dh.n_dofs(), dh.n_dofs()); + + // The kernel. + LaplaceKernelIntegration kernel; + + // i runs on points, j runs on cells. + for(; celli != endc; ++celli, ++cellv) { + fe_mid.reinit(celli); + Point ci = celli->center(); + Point ni = fe_mid.cell_normal_vector(0); + + celli->get_dof_indices(dofsi); + cellv->get_dof_indices(dofsv); + + // Vn vector: + Vn(dofsi[0]) = wind*ni; + vs(dofsv[0]) = ni[0]; + vs(dofsv[1]) = ni[1]; + vs(dofsv[2]) = ni[2]; + + // Now the two matrices. + for(cellj = dh.begin_active(); cellj != endc; ++cellj) { + vector SD(2,0.); // Single and Double layers. + cellj->get_dof_indices(dofsj); + kernel.compute_SD_integral_on_cell(SD, + cellj, + ci); + _B (dofsi[0], dofsj[0]) += -SD[0]; + if(dofsi[0] != dofsj[0]) + (*system_matrix)(dofsi[0], dofsj[0]) += -SD[1]; + if(dofsi[0] == dofsj[0]) + (*system_matrix)(dofsi[0], dofsj[0]) += 1.; + } + } + _B.vmult(rhs, Vn); +} + +template +void BEMProblem::solve_system() { + phi.swap(rhs); + system_matrix->compute_lu_factorization(); + system_matrix->apply_lu_factorization(phi, false); +} + + +template +void BEMProblem::output_results(unsigned int cycle) { + + DataOut > dataout; + + dataout.attach_dof_handler(dh); + dataout.add_data_vector(phi, "phi"); + dataout.build_patches(); + + char fname[100]; + sprintf(fname, "test_%02d.vtk", cycle); + std::ofstream file(fname); + + dataout.write_vtk(file); +} + +template +void BEMProblem::run() { + read_domain(); + + const unsigned int number_of_cycles = 4; + + for(unsigned int cycle=0; cycle laplace_problem; + laplace_problem.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} -- 2.39.5