From: heltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 23 Feb 2009 08:32:15 +0000 (+0000)
Subject: Added preliminary version of step-34. BEM.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=58816da7d80932418aac96397843b2264cd2842f;p=dealii-svn.git

Added preliminary version of step-34. BEM.

git-svn-id: https://svn.dealii.org/trunk@18413 0785d39b-7218-0410-832d-ea1e28bc413d
---

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========= $(<F)
+	@$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+	@echo ==============optimized===== $(<F)
+	@$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule to created a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# If the creation of Makefile.dep fails, blow it away and fail
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $D/*/include/*/*.h)
+	@echo ============================ Remaking $@
+	@$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
+		> $@ \
+	  || (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 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+<h3>The problem</h3>
+
+<h3>Galerkin BEM</h3>
+
+<h3>What the program does</h3>
+
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 @@
+<h1>Results</h1>
+
+<h3>Simple Example</h3>
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 <fstream>
+#include <base/logstream.h>
+
+#include <base/convergence_table.h>
+#include <base/quadrature_lib.h>
+#include <base/table.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_renumbering.h>
+#include <fe/fe_dgp.h>
+#include <fe/fe_system.h>
+#include <fe/fe_tools.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_in.h>
+#include <grid/grid_out.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/tria_iterator.h>
+#include <lac/full_matrix.h>
+#include <lac/precondition.h>
+#include <lac/solver_cg.h>
+#include <lac/vector.h>
+#include <lac/sparse_direct.h>
+#include <lac/lapack_full_matrix.h>
+#include <numerics/data_out.h>
+#include <base/smartpointer.h>
+
+#include <cmath>
+#include <iostream>
+#include <math.h>
+#include <string>
+
+using namespace std;
+using namespace dealii;
+
+template <int dim>
+class LaplaceKernelIntegration
+{
+public:
+
+    LaplaceKernelIntegration();
+    ~LaplaceKernelIntegration();
+  
+    void run();
+    
+    void compute_SD_integral_on_cell(vector<double> &dst,
+				     typename DoFHandler<dim,dim+1>::active_cell_iterator &cell,
+				     const Point<dim+1> &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<FEValues<dim,dim+1> > fe_values;
+};
+
+template <int dim> 
+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<dim-1, dim> 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<dim-1,dim> dh;
+    
+    /** The velocity degrees of freedom. */
+    DoFHandler<dim-1,dim> 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<LAPACKFullMatrix<double> > system_matrix;
+    
+    /** Single layer potential matrix. */
+    FullMatrix<double> single_layer_matrix;
+    
+    /** Normal component of the wind field. */
+    Vector<double> Vn;
+    
+    /** System rhs. */
+    Vector<double> rhs;
+    
+    /** Potential. */
+    Vector<double> phi;
+    
+    /** Something else. */
+    Vector<double> vs;
+};
+
+template <int dim>
+BEMProblem<dim>::BEMProblem() :
+    fe(0),
+    fev(FE_DGP<dim-1,dim>(0), dim),
+    dh(tria),
+    dhv(tria)
+{}
+
+
+template <int dim>
+BEMProblem<dim>::~BEMProblem() {
+    LAPACKFullMatrix<double> * p = system_matrix;
+    system_matrix = 0;
+    delete p;
+}
+
+
+template <>
+LaplaceKernelIntegration<2>::LaplaceKernelIntegration()
+{
+    static FE_DGP<2,3> fe(0);
+    vector<Point<2> > 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<double> 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 <int dim>
+LaplaceKernelIntegration<dim>::~LaplaceKernelIntegration() {
+    FEValues<dim,dim+1> *fp = fe_values;
+    fe_values = 0;
+    delete fp;
+}
+
+
+template <>
+void
+LaplaceKernelIntegration<2>::compute_SD_integral_on_cell(vector<double> &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<Tensor<2,3> > jacobians = fe_values->get_jacobians();
+    vector<Point<3> > 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<double> i_S(4);
+    vector<double> 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 <int dim>
+double
+LaplaceKernelIntegration<dim>::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 <int dim>
+double
+LaplaceKernelIntegration<dim>::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 <int dim>
+void BEMProblem<dim>::read_domain() {
+    // Center of the ball. It is the origin by default.
+    Point<dim> p;
+    static HyperBallBoundary<dim-1, dim> boundary(p,1.);    
+
+    // Read the sphere from
+    GridIn<dim-1, dim> 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 <int dim>
+void BEMProblem<dim>::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<double> * p = system_matrix;
+	system_matrix = 0;
+	delete p;
+    }
+    
+    system_matrix = new LAPACKFullMatrix<double>(ndofs, ndofs);
+
+    Vn.reinit(ndofs);
+    rhs.reinit(ndofs);
+    phi.reinit(ndofs);
+    vs.reinit(nvdofs);
+}    
+
+template <int dim>
+void BEMProblem<dim>::assemble_system() {
+    
+    Point<dim> wind;
+    wind[0] = 1.;
+    
+    typename DoFHandler<dim-1,dim>::active_cell_iterator
+	celli = dh.begin_active(),
+	cellj = dh.begin_active(),
+	cellv = dhv.begin_active(),
+	endc = dh.end();
+    
+    QMidpoint<dim-1> midpoint;
+    FEValues<dim-1,dim> fe_mid(fe, midpoint,
+			       update_values |
+			       update_cell_normal_vectors |
+			       update_quadrature_points);
+    
+    vector<unsigned int> dofsi(fe_mid.n_quadrature_points);
+    vector<unsigned int> dofsj(fe_mid.n_quadrature_points);
+    vector<unsigned int> dofsv(dim);
+    
+    // Temporary matrix 
+    FullMatrix<double> _B(dh.n_dofs(), dh.n_dofs());
+    
+    // The kernel.
+    LaplaceKernelIntegration<dim-1> kernel;
+    
+    // i runs on points, j runs on cells.
+    for(; celli != endc; ++celli, ++cellv) {
+	    fe_mid.reinit(celli);
+	    Point<dim> ci = celli->center();
+	    Point<dim> 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<double> 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 <int dim>
+void BEMProblem<dim>::solve_system() {
+    phi.swap(rhs);
+    system_matrix->compute_lu_factorization();
+    system_matrix->apply_lu_factorization(phi, false);
+}
+
+
+template <int dim>
+void BEMProblem<dim>::output_results(unsigned int cycle) {
+    
+    DataOut<dim-1, DoFHandler<dim-1, dim> > 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 <int dim>
+void BEMProblem<dim>::run() {
+    read_domain();
+    
+    const unsigned int number_of_cycles = 4;
+    
+    for(unsigned int cycle=0; cycle<number_of_cycles; ++cycle) {
+	refine_and_resize();
+	assemble_system();
+	solve_system();
+	output_results(cycle);
+    }
+}
+
+
+int main () 
+{
+  try
+  {
+      deallog.depth_console (3);
+      
+      BEMProblem<3> 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;
+}