From: heltai Date: Fri, 3 Dec 2010 16:36:49 +0000 (+0000) Subject: Added tutorial step-38 on Laplace Beltrami Operator. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1765077c9d8b6f4c0ce44da160daa9ba4db0c0ec;p=dealii-svn.git Added tutorial step-38 on Laplace Beltrami Operator. git-svn-id: https://svn.dealii.org/trunk@22911 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-38/Makefile b/deal.II/examples/step-38/Makefile new file mode 100644 index 0000000000..e72451dab1 --- /dev/null +++ b/deal.II/examples/step-38/Makefile @@ -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========= $( $@" + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +./%.$(OBJEXT) : + @echo "==============optimized===== $( $@" + @$(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 index 0000000000..48a0f73876 --- /dev/null +++ b/deal.II/examples/step-38/doc/builds-on @@ -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 index 0000000000..eae18a9fea --- /dev/null +++ b/deal.II/examples/step-38/doc/intro.dox @@ -0,0 +1,15 @@ +
+ +This program was contributed by Andrea Bonito + + + + +

Introduction

+ +

The numerical approximation

+ +

Implementation

+ +

Testcase

+ diff --git a/deal.II/examples/step-38/doc/kind b/deal.II/examples/step-38/doc/kind new file mode 100644 index 0000000000..15a13db451 --- /dev/null +++ b/deal.II/examples/step-38/doc/kind @@ -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 index 0000000000..7ccb36ef7a --- /dev/null +++ b/deal.II/examples/step-38/doc/results.dox @@ -0,0 +1,25 @@ +

Results

+ +We ran the program using the following parameters.prm 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 index 0000000000..f1a5e8797a 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 index 0000000000..550a7eb4e4 --- /dev/null +++ b/deal.II/examples/step-38/doc/tooltip @@ -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 index 0000000000..33516965c6 --- /dev/null +++ b/deal.II/examples/step-38/step-38.cc @@ -0,0 +1,466 @@ +#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 + +#define deal_II_dimension 3 + +using std::cout; +using std::endl; +using namespace dealii; + + +template +class Solution : public Function +{ +public: + Solution () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + +}; + +template +double Solution::value (const Point &p, + const unsigned int) const +{ + return sin(numbers::PI * p(0))*cos(numbers::PI * p(1))*exp(p(2)); +} + +template +Tensor<1,dim> Solution::gradient (const Point &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 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 +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + +template +double RightHandSide::value (const Point &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 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 +class LaplaceBeltrami +{ + public: + LaplaceBeltrami (Triangulation *tria, Function &func_data, + unsigned int fe_degree = 1, unsigned int mapping_degree = 1, + Function *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 *pTria; + FE_Q fe; + DoFHandler dh; + MappingQ mapping; + + ConstraintMatrix matrix_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + Function &rhs_func; // function data + + Function *pExact; // exact solution if provided + +}; + + +template +LaplaceBeltrami::LaplaceBeltrami (Triangulation *tria, + Function &func_data, + unsigned int fe_degree, + unsigned int mapping_degree, + Function *pExact + ) + : + fe (fe_degree), + dh(*tria), + mapping(mapping_degree), + rhs_func(func_data) +{ + pTria = tria; + this->pExact = pExact; + +} + +template +LaplaceBeltrami::~LaplaceBeltrami () +{ + dh.clear (); +} + +template +void LaplaceBeltrami::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 +void LaplaceBeltrami::assemble_system () +{ + + system_matrix = 0; + system_rhs = 0; + + QGauss quadrature_formula(2); + + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector< double > rhs_values(n_q_points); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i 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 +void LaplaceBeltrami::solve () +{ + SolverControl solver_control (1000, 1e-7); + SolverCG<> cg (solver_control); + + cg.solve (system_matrix, solution, system_rhs, + PreconditionIdentity()); + + + matrix_constraints.distribute (solution); + +} + + + +template +void LaplaceBeltrami::output_results () const +{ + + std::string filename = "solution.vtk"; + + std::ofstream output (filename.c_str()); + + DataOut > data_out; + data_out.attach_dof_handler (dh); + + + + data_out.add_data_vector (solution, "solution",DataOut_DoFData,dim-1,dim>::type_dof_data); + data_out.build_patches (mapping, + mapping.get_degree()); + data_out.write_vtk (output); +} + + + +template +void LaplaceBeltrami::run () +{ + + setup_system(); + assemble_system(); + solve(); + output_results(); + +} + +//################################################################################// + +template +double LaplaceBeltrami::compute_error (VectorTools::NormType norm_type) const +{ + Assert(pExact != 0, ExcInternalError()); + + Vector 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< tria; + + // create a mesh consisting on the boundary of the half sphere... thx Seba + std::map< Triangulation::cell_iterator, + Triangulation::face_iterator> surface_to_volume_mapping; + + HyperBallBoundary boundary_description; + Triangulation 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 surface_description; + tria.set_boundary (1, surface_description); + tria.set_boundary (0, surface_description); + + std::set boundary_ids; + boundary_ids.insert(0); + + GridTools::extract_boundary_mesh (volume_mesh, tria, + surface_to_volume_mapping, + boundary_ids); + + + tria.refine_global(4); + + + RightHandSide rhs; + Solution exact; + + LaplaceBeltrami laplace_beltrami_2d(&tria,rhs,2,2,&exact); + + laplace_beltrami_2d.run(); + + std::cout<