]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added tutorial step-38 on Laplace Beltrami Operator.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 3 Dec 2010 16:36:49 +0000 (16:36 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 3 Dec 2010 16:36:49 +0000 (16:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@22911 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-38/Makefile [new file with mode: 0644]
deal.II/examples/step-38/doc/builds-on [new file with mode: 0644]
deal.II/examples/step-38/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-38/doc/kind [new file with mode: 0644]
deal.II/examples/step-38/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-38/doc/step-XX_0.png [new file with mode: 0644]
deal.II/examples/step-38/doc/tooltip [new file with mode: 0644]
deal.II/examples/step-38/step-38.cc [new file with mode: 0755]

diff --git a/deal.II/examples/step-38/Makefile b/deal.II/examples/step-38/Makefile
new file mode 100644 (file)
index 0000000..e72451d
--- /dev/null
@@ -0,0 +1,142 @@
+# $Id: Makefile 22601 2010-11-04 03:26:19Z bangerth $
+
+
+# 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 *vtk *ucd *.d2
+
+
+
+
+#
+#
+# Usually, you will not need to change anything beyond this point.
+#
+#
+# The next statement tells 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. deal.II has two
+# libraries: one for the debug mode version of the
+# application and one for optimized mode.
+libs.g   := $(lib-deal2.g)
+libs.o   := $(lib-deal2.o)
+
+
+# We now use the variable defined above to 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 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)$(EXEEXT) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@ $^ $(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)$(EXEEXT)
+       @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 creates 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/deal.II/*/*.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-38/doc/builds-on b/deal.II/examples/step-38/doc/builds-on
new file mode 100644 (file)
index 0000000..48a0f73
--- /dev/null
@@ -0,0 +1 @@
+step-4
diff --git a/deal.II/examples/step-38/doc/intro.dox b/deal.II/examples/step-38/doc/intro.dox
new file mode 100644 (file)
index 0000000..eae18a9
--- /dev/null
@@ -0,0 +1,15 @@
+<br>
+
+<i>This program was contributed by Andrea Bonito
+</i>
+
+<a name="Intro"></a>
+
+<h1>Introduction</h1>
+
+<h3>The numerical approximation</h3>
+
+<h3>Implementation</h3>
+
+<h3>Testcase</h3>
+
diff --git a/deal.II/examples/step-38/doc/kind b/deal.II/examples/step-38/doc/kind
new file mode 100644 (file)
index 0000000..15a13db
--- /dev/null
@@ -0,0 +1 @@
+basic
diff --git a/deal.II/examples/step-38/doc/results.dox b/deal.II/examples/step-38/doc/results.dox
new file mode 100644 (file)
index 0000000..7ccb36e
--- /dev/null
@@ -0,0 +1,25 @@
+<h1>Results</h1>
+
+We ran the program using the following <code>parameters.prm</code> file (which
+can also be found in the directory in which all the other source files are):
+@verbatim
+# Listing of Parameters
+# ---------------------
+@endverbatim
+
+When we run the program, the following is printed on screen:
+
+@verbatim
+================================
+          LB PROBLEM 
+================================
+
+DEAL:cg::Starting value 0.677406
+DEAL:cg::Convergence step 580 value 9.71224e-08
+0.00552653
+
+@endverbatim
+
+Here is a plot of the result:
+
+@image html step-XX_0.png
diff --git a/deal.II/examples/step-38/doc/step-XX_0.png b/deal.II/examples/step-38/doc/step-XX_0.png
new file mode 100644 (file)
index 0000000..f1a5e87
Binary files /dev/null and b/deal.II/examples/step-38/doc/step-XX_0.png differ
diff --git a/deal.II/examples/step-38/doc/tooltip b/deal.II/examples/step-38/doc/tooltip
new file mode 100644 (file)
index 0000000..550a7eb
--- /dev/null
@@ -0,0 +1 @@
+Solve the Laplace Beltrami operator on a Half Sphere.
diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc
new file mode 100755 (executable)
index 0000000..3351696
--- /dev/null
@@ -0,0 +1,466 @@
+#include <base/polynomial_space.h>
+#include <base/parsed_function.h>
+#include <base/smartpointer.h>
+#include <base/convergence_table.h>
+#include <base/quadrature_lib.h>
+#include <base/quadrature_selector.h>
+#include <base/utilities.h>
+
+#include <lac/full_matrix.h>
+#include <lac/vector.h>
+#include <lac/solver_control.h>
+#include <lac/solver_gmres.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <lac/constraint_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/compressed_sparsity_pattern.h>
+
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_tools.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_values.h>
+
+#include <numerics/data_out.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+#include <fstream>
+#include <iostream>
+
+#include <sys/types.h>
+#include <sys/stat.h>
+
+#define deal_II_dimension 3
+
+using std::cout;
+using std::endl;
+using namespace dealii;
+
+
+template <int dim>
+class Solution  : public Function<dim>
+{
+public:
+  Solution () : Function<dim>() {}
+  
+  virtual double value (const Point<dim>   &p,
+                       const unsigned int  component = 0) const;
+  
+  virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                 const unsigned int  component = 0) const;
+
+};
+template <int dim>
+double Solution<dim>::value (const Point<dim> &p,
+                         const unsigned int) const 
+{
+  return sin(numbers::PI * p(0))*cos(numbers::PI * p(1))*exp(p(2)); 
+}
+
+template <int dim>
+Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
+                                      const unsigned int) const
+{
+  double dPi = numbers::PI;
+
+  Tensor<1,dim> return_value;
+
+  return_value[0] = dPi *cos(dPi * p(0))*cos(dPi * p(1))*exp(p(2));
+  return_value[1] = -dPi *sin(dPi * p(0))*sin(dPi * p(1))*exp(p(2));
+  return_value[2] = sin(dPi * p(0))*cos(dPi * p(1))*exp(p(2));
+  
+  // tangential gradient: nabla u - (nabla u nu)nu
+  Point<dim> normal;
+  double dLength;  
+
+  dLength = sqrt(p(0)*p(0)+p(1)*p(1)+p(2)*p(2));
+  
+  normal[0] = p(0)/dLength;
+  normal[1] = p(1)/dLength;
+  normal[2] = p(2)/dLength;
+  
+  return return_value - (return_value*normal)*normal;
+}
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim>() {}
+  
+  virtual double value (const Point<dim>   &p,
+                       const unsigned int  component = 0) const;
+};
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+                                 const unsigned int comp) const 
+{
+  Assert(dim == 3,  ExcInternalError());
+  
+  double dPi = numbers::PI;
+
+  // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
+
+  Tensor<2,dim> hessian;
+
+  hessian[0][0] = -dPi*dPi*sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2));
+  hessian[1][1] = -dPi*dPi*sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2));
+  hessian[2][2] = sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2));
+
+  hessian[0][1] = -dPi*dPi*cos(dPi*p(0))*sin(dPi*p(1))*exp(p(2));
+  hessian[1][0] = -dPi*dPi*cos(dPi*p(0))*sin(dPi*p(1))*exp(p(2));
+
+  hessian[0][2] = dPi*cos(dPi*p(0))*cos(dPi*p(1))*exp(p(2));
+  hessian[2][0] = dPi*cos(dPi*p(0))*cos(dPi*p(1))*exp(p(2));
+
+  hessian[1][2] = -dPi*sin(dPi*p(0))*sin(dPi*p(1))*exp(p(2));
+  hessian[2][1] = -dPi*sin(dPi*p(0))*sin(dPi*p(1))*exp(p(2));
+
+  Tensor<1,dim> gradient;
+  gradient[0] = dPi * cos(dPi*p(0))*cos(dPi*p(1))*exp(p(2));
+  gradient[1] = - dPi * sin(dPi*p(0))*sin(dPi*p(1))*exp(p(2));
+  gradient[2] = sin(dPi*p(0))*cos(dPi*p(1))*exp(p(2));
+
+  double curvature;
+  Point<dim> normal;
+  double dLength;  
+
+  curvature = dim-1;
+  dLength = sqrt(p(0)*p(0)+p(1)*p(1)+p(2)*p(2));
+  
+  normal[0] = p(0)/dLength;
+  normal[1] = p(1)/dLength;
+  normal[2] = p(2)/dLength;
+  
+  return -trace(hessian) + (hessian * normal) * normal + (gradient * normal)*curvature;
+}
+
+
+
+
+
+
+template <int dim>
+class LaplaceBeltrami 
+{
+  public:
+  LaplaceBeltrami (Triangulation<dim-1,dim> *tria, Function<dim> &func_data, 
+                  unsigned int fe_degree = 1, unsigned int mapping_degree = 1,
+                  Function<dim> *pExact = 0);
+  // arguments are:
+  // triangulation
+  // right-hand side
+  // fe_degree for solution
+  // fe_degree for mapping
+  // exact solution is known
+  
+  ~LaplaceBeltrami ();
+  void run ();
+  double compute_error (VectorTools::NormType norm_type) const;
+  
+  private:
+  
+  void setup_system ();
+  void assemble_system ();
+  void solve ();
+  void output_results () const;
+  
+  
+  
+  Triangulation<dim-1,dim>   *pTria;
+  FE_Q<dim-1,dim>            fe;
+  DoFHandler<dim-1,dim>      dh;
+  MappingQ<dim-1, dim>       mapping;
+
+  ConstraintMatrix     matrix_constraints;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+  
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
+
+  Function<dim> &rhs_func; // function data
+
+  Function<dim> *pExact; // exact solution if provided
+
+};
+
+
+template <int dim>
+LaplaceBeltrami<dim>::LaplaceBeltrami (Triangulation<dim-1,dim> *tria,
+                                      Function<dim> &func_data,
+                                      unsigned int fe_degree,
+                                      unsigned int mapping_degree,
+                                      Function<dim> *pExact
+                                      )
+               :
+  fe (fe_degree),
+  dh(*tria),
+  mapping(mapping_degree),
+  rhs_func(func_data)
+{
+  pTria = tria;
+  this->pExact = pExact;
+
+}
+
+template <int dim>
+LaplaceBeltrami<dim>::~LaplaceBeltrami () 
+{
+  dh.clear ();
+}
+
+template <int dim>
+void LaplaceBeltrami<dim>::setup_system ()
+{
+  dh.distribute_dofs (fe);
+
+  matrix_constraints.clear ();
+  matrix_constraints.close ();
+
+  CompressedSparsityPattern csp (dh.n_dofs(),
+                                 dh.n_dofs());
+
+  DoFTools::make_sparsity_pattern (dh, csp);
+  matrix_constraints.condense (csp);
+
+  sparsity_pattern.copy_from (csp);
+
+  system_matrix.reinit (sparsity_pattern);
+  solution.reinit (dh.n_dofs());
+  system_rhs.reinit (dh.n_dofs());
+}
+
+template <int dim>
+void LaplaceBeltrami<dim>::assemble_system () 
+{  
+
+  system_matrix = 0;
+  system_rhs = 0;
+  
+  QGauss<dim-1>  quadrature_formula(2);
+
+  FEValues<dim-1,dim> fe_values (mapping,fe, quadrature_formula, 
+                                update_values   | update_cell_normal_vectors |
+                                update_gradients |
+                                update_quadrature_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector< double > rhs_values(n_q_points);
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  typename DoFHandler<dim-1,dim>::active_cell_iterator
+    cell = dh.begin_active(),
+    endc = dh.end();
+
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      fe_values.reinit (cell);
+
+      rhs_func.value_list (fe_values.get_quadrature_points(), rhs_values); 
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+      {          
+       for (unsigned int j=0; j<dofs_per_cell; ++j) 
+       {             
+         for (unsigned int q_point=0; q_point<n_q_points;
+              ++q_point)
+         {
+           cell_matrix(i,j) 
+             += fe_values.shape_grad(i,q_point) *
+             fe_values.shape_grad(j,q_point) *
+             fe_values.JxW(q_point);
+         }
+       }
+      }
+
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+      {
+         
+       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+         cell_rhs(i) += fe_values.shape_value(i,q_point) *
+           rhs_values[q_point]*
+           fe_values.JxW(q_point);
+      }
+
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+      {
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         system_matrix.add (local_dof_indices[i],
+                            local_dof_indices[j],
+                            cell_matrix(i,j));
+       
+       system_rhs(local_dof_indices[i]) += cell_rhs(i);
+      }
+    }
+
+
+   std::map<unsigned int,double> bdy_values; 
+   VectorTools::interpolate_boundary_values (mapping,dh,
+                                            0,
+                                            *pExact,
+                                            bdy_values
+                                            );
+
+   MatrixTools::apply_boundary_values (bdy_values,
+                                      system_matrix,
+                                      solution,
+                                      system_rhs,false);
+
+  // condense matrices
+  matrix_constraints.condense (system_matrix);
+  matrix_constraints.condense (system_rhs);
+}
+
+
+template <int dim>
+void LaplaceBeltrami<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-7);
+  SolverCG<>              cg (solver_control);
+
+  cg.solve (system_matrix, solution, system_rhs,
+           PreconditionIdentity());
+
+
+  matrix_constraints.distribute (solution);
+
+}
+
+
+
+template <int dim>
+void LaplaceBeltrami<dim>::output_results () const
+{
+  
+  std::string filename = "solution.vtk";
+
+  std::ofstream output (filename.c_str());
+
+  DataOut<dim-1,DoFHandler<dim-1,dim> > data_out;
+  data_out.attach_dof_handler (dh);
+
+
+  data_out.add_data_vector (solution, "solution",DataOut_DoFData<DoFHandler<dim-1,dim>,dim-1,dim>::type_dof_data);
+  data_out.build_patches (mapping,
+                         mapping.get_degree());
+  data_out.write_vtk (output);
+}
+
+
+
+template <int dim>
+void LaplaceBeltrami<dim>::run () 
+{
+
+  setup_system();
+  assemble_system();
+  solve();
+  output_results();
+
+}
+
+//################################################################################//
+
+template <int dim>
+double LaplaceBeltrami<dim>::compute_error (VectorTools::NormType norm_type) const
+{  
+  Assert(pExact != 0, ExcInternalError());
+  
+  Vector<float> difference_per_cell (pTria->n_active_cells());
+  VectorTools::integrate_difference (mapping, dh, solution,
+                                    *pExact, difference_per_cell,
+                                    QGauss<(dim-1)>(2*fe.degree+1),
+                                    norm_type);
+
+  return difference_per_cell.l2_norm();
+}
+
+
+
+
+int main ( int argc, char **argv )
+{
+  std::cout<<std::endl;
+  std::cout<<"================================"<<std::endl;
+  std::cout<<"          LB PROBLEM "<<std::endl;
+  std::cout<<"================================"<<std::endl;
+  std::cout<<std::endl;
+  
+  
+  Triangulation<deal_II_dimension-1,deal_II_dimension> tria;
+  
+  // create a mesh consisting on the boundary of the half sphere... thx Seba
+  std::map< Triangulation<deal_II_dimension-1,deal_II_dimension>::cell_iterator,
+    Triangulation<deal_II_dimension,deal_II_dimension>::face_iterator> surface_to_volume_mapping;
+
+  HyperBallBoundary<deal_II_dimension> boundary_description;
+  Triangulation<deal_II_dimension> volume_mesh;
+  GridGenerator::half_hyper_ball(volume_mesh);
+  
+  volume_mesh.set_boundary (1, boundary_description);
+  volume_mesh.set_boundary (0, boundary_description);
+  volume_mesh.refine_global (1);
+  
+  static HyperBallBoundary<deal_II_dimension-1,deal_II_dimension> surface_description;
+  tria.set_boundary (1, surface_description);
+  tria.set_boundary (0, surface_description);
+  
+  std::set<unsigned char> boundary_ids;
+  boundary_ids.insert(0);
+  
+  GridTools::extract_boundary_mesh (volume_mesh, tria,
+                                   surface_to_volume_mapping,
+                                   boundary_ids);
+
+  tria.refine_global(4);
+      
+
+  RightHandSide<deal_II_dimension> rhs;
+  Solution<deal_II_dimension> exact;
+  
+  LaplaceBeltrami<deal_II_dimension> laplace_beltrami_2d(&tria,rhs,2,2,&exact);
+  
+  laplace_beltrami_2d.run();
+  
+  std::cout<<laplace_beltrami_2d.compute_error(VectorTools::H1_norm)<<std::endl;
+
+  tria.set_boundary(0);
+  return 0;
+}
+
+
+
+
+
+
+
+
+
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.