From: bangerth Date: Sun, 27 Jan 2008 04:39:35 +0000 (+0000) Subject: add X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=83227b92c94d68dc5d6a6e85bb5aee77bbe07f59;p=dealii-svn.git add git-svn-id: https://svn.dealii.org/trunk@15691 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/Makefile b/deal.II/examples/step-31/Makefile new file mode 100644 index 0000000000..26266ada33 --- /dev/null +++ b/deal.II/examples/step-31/Makefile @@ -0,0 +1,156 @@ +# $Id: Makefile,v 1.4 2006/02/10 17:53:05 wolf Exp $ + + +# 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 = off + + +# 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 something 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-2d.g) \ + $(lib-deal2-3d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-2d.o) \ + $(lib-deal2-3d.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-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox new file mode 100644 index 0000000000..c59cfdf609 --- /dev/null +++ b/deal.II/examples/step-31/doc/intro.dox @@ -0,0 +1,2 @@ + +

Introduction

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

Results

diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc new file mode 100644 index 0000000000..1fea326201 --- /dev/null +++ b/deal.II/examples/step-31/step-31.cc @@ -0,0 +1,792 @@ +/* $Id: step-22.cc 15679 2008-01-24 23:28:37Z bangerth $ */ +/* Author: Wolfgang Bangerth, Texas A&M University, 2008 */ + +/* $Id: step-22.cc 15679 2008-01-24 23:28:37Z bangerth $ */ +/* Version: $Name$ */ +/* */ +/* Copyright (C) 2008 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. */ + + + +#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 dealii; + + + +template +class StokesProblem +{ + public: + StokesProblem (const unsigned int degree); + void run (); + + private: + void setup_dofs (); + void assemble_system (); + void solve (); + void output_results (const unsigned int refinement_cycle) const; + void refine_mesh (); + + const unsigned int degree; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + ConstraintMatrix hanging_node_constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + + boost::shared_ptr A_preconditioner; +}; + + + + + +template +class PressureBoundaryValues : public Function +{ + public: + PressureBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + +template +double +PressureBoundaryValues::value (const Point &/*p*/, + const unsigned int /*component*/) const +{ + return 0; +} + + + +template +class BoundaryValues : public Function +{ + public: + BoundaryValues () : Function(dim+1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; +}; + + +template +double +BoundaryValues::value (const Point &p, + const unsigned int component) const +{ + if (component == 0) + return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0)); + return 0; +} + + +template +void +BoundaryValues::vector_value (const Point &p, + Vector &values) const +{ + for (unsigned int c=0; cn_components; ++c) + values(c) = BoundaryValues::value (p, c); +} + + + +template +class RightHandSide : public Function +{ + public: + RightHandSide () : Function(dim+1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int component) const +{ + return 0; +} + + +template +void +RightHandSide::vector_value (const Point &p, + Vector &values) const +{ + for (unsigned int c=0; cn_components; ++c) + values(c) = RightHandSide::value (p, c); +} + + + + + + +template +Tensor<1,dim> +extract_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<1,dim> tmp; + + const unsigned int component + = fe_values.get_fe().system_to_component_index(i).first; + + if (component < dim) + tmp[component] = fe_values.shape_value (i,q); + + return tmp; +} + + + +template +double +extract_div_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + const unsigned int component + = fe_values.get_fe().system_to_component_index(i).first; + + if (component < dim) + return fe_values.shape_grad (i,q)[component]; + else + return 0; +} + + + +template +Tensor<2,dim> +extract_grad_s_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<2,dim> tmp; + + const unsigned int component + = fe_values.get_fe().system_to_component_index(i).first; + + if (component < dim) + { + const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2; + + for (unsigned int e=0; e +double extract_p (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + const unsigned int component + = fe_values.get_fe().system_to_component_index(i).first; + + if (component == dim) + return fe_values.shape_value (i,q); + else + return 0; +} + + + + +template +class InverseMatrix : public Subscriptor +{ + public: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + const Preconditioner &preconditioner; + + mutable GrowingVectorMemory<> vector_memory; +}; + + +template +InverseMatrix::InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (preconditioner) +{} + + + +template +void InverseMatrix::vmult (Vector &dst, + const Vector &src) const +{ + SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); + SolverCG<> cg (solver_control, vector_memory); + + dst = 0; + + try + { + cg.solve (*matrix, dst, src, preconditioner); + } + catch (std::exception &e) + { + Assert (false, ExcMessage(e.what())); + } +} + + + +template +class SchurComplement : public Subscriptor +{ + public: + SchurComplement (const BlockSparseMatrix &A, + const InverseMatrix,Preconditioner> &Minv); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + const SmartPointer,Preconditioner> > m_inverse; + + mutable Vector tmp1, tmp2; +}; + + + +template +SchurComplement:: +SchurComplement (const BlockSparseMatrix &A, + const InverseMatrix,Preconditioner> &Minv) + : + system_matrix (&A), + m_inverse (&Minv), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) +{} + + +template +void SchurComplement::vmult (Vector &dst, + const Vector &src) const +{ + system_matrix->block(0,1).vmult (tmp1, src); + m_inverse->vmult (tmp2, tmp1); + system_matrix->block(1,0).vmult (dst, tmp2); +} + + + +template +StokesProblem::StokesProblem (const unsigned int degree) + : + degree (degree), + triangulation (Triangulation::maximum_smoothing), + fe (FE_Q(degree+1), dim, + FE_Q(degree), 1), + dof_handler (triangulation) +{} + + + + +template +void StokesProblem::setup_dofs () +{ + dof_handler.distribute_dofs (fe); + DoFRenumbering::component_wise (dof_handler); + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + + std::vector dofs_per_component (dim+1); + DoFTools::count_dofs_per_component (dof_handler, dofs_per_component); + const unsigned int n_u = dofs_per_component[0] * dim, + n_p = dofs_per_component[dim]; + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << " (" << n_u << '+' << n_p <<')' + << std::endl; + + const unsigned int + n_couplings = dof_handler.max_couplings_between_dofs(); + + system_matrix.clear (); + + sparsity_pattern.reinit (2,2); + sparsity_pattern.block(0,0).reinit (n_u, n_u, n_couplings); + sparsity_pattern.block(1,0).reinit (n_p, n_u, n_couplings); + sparsity_pattern.block(0,1).reinit (n_u, n_p, n_couplings); + sparsity_pattern.block(1,1).reinit (n_p, n_p, n_couplings); + + sparsity_pattern.collect_sizes(); + + + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + hanging_node_constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (2); + solution.block(0).reinit (n_u); + solution.block(1).reinit (n_p); + solution.collect_sizes (); + + system_rhs.reinit (2); + system_rhs.block(0).reinit (n_u); + system_rhs.block(1).reinit (n_p); + system_rhs.collect_sizes (); +} + + +template +double +scalar_product (const Tensor<2,dim> &a, + const Tensor<2,dim> &b) +{ + double tmp = 0; + for (unsigned int i=0; i +void StokesProblem::assemble_system () +{ + system_matrix=0; + system_rhs=0; + + QGauss quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); + + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); + FEFaceValues fe_face_values (fe, face_quadrature_formula, + update_values | update_normal_vectors | + 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(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const PressureBoundaryValues pressure_boundary_values; + + std::vector boundary_values (n_face_q_points); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + local_rhs = 0; + + for (unsigned int q=0; q phi_i_u = extract_u (fe_values, i, q); + + const Tensor<2,dim> phi_i_grads_u= extract_grad_s_u (fe_values, i, q); + const double div_phi_i_u = extract_div_u (fe_values, i, q); + const double phi_i_p = extract_p (fe_values, i, q); + + for (unsigned int j=0; j phi_j_grads_u = extract_grad_s_u (fe_values, j, q); + const double div_phi_j_u = extract_div_u (fe_values, j, q); + const double phi_j_p = extract_p (fe_values, j, q); + + local_matrix(i,j) += (scalar_product(phi_i_grads_u, phi_j_grads_u) + - div_phi_i_u * phi_j_p + - phi_i_p * div_phi_j_u + + phi_i_p * phi_j_p) + * fe_values.JxW(q); + } + } + } + + + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + if (cell->at_boundary(face_no)) + { + fe_face_values.reinit (cell, face_no); + + pressure_boundary_values + .value_list (fe_face_values.get_quadrature_points(), + boundary_values); + + for (unsigned int q=0; q + phi_i_u = extract_u (fe_face_values, i, q); + + local_rhs(i) += -(phi_i_u * + fe_face_values.normal_vector(q) * + boundary_values[q] * + fe_face_values.JxW(q)); + } + } + + cell->get_dof_indices (local_dof_indices); + + for (unsigned int i=0; i boundary_values; + std::vector component_mask (dim+1, true); + component_mask[dim] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 3, + BoundaryValues(), + boundary_values, + component_mask); + + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); + } + + std::cout << " Computing preconditioner..." << std::flush; + + A_preconditioner + = boost::shared_ptr(new SparseDirectUMFPACK()); + A_preconditioner->initialize (system_matrix.block(0,0)); + + std::cout << std::endl; +} + + + +template +void StokesProblem::solve () +{ + const InverseMatrix,SparseDirectUMFPACK> + A_inverse (system_matrix.block(0,0), *A_preconditioner); + Vector tmp (solution.block(0).size()); + Vector schur_rhs (solution.block(1).size()); + + { + A_inverse.vmult (tmp, system_rhs.block(0)); + system_matrix.block(1,0).vmult (schur_rhs, tmp); + schur_rhs -= system_rhs.block(1); + + SchurComplement + schur_complement (system_matrix, A_inverse); + + SolverControl solver_control (system_matrix.block(0,0).m(), + 1e-6*schur_rhs.l2_norm()); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize (system_matrix.block(1,1), 1.2); + + InverseMatrix,PreconditionSSOR<> > + m_inverse (system_matrix.block(1,1), preconditioner); + + try + { + cg.solve (schur_complement, solution.block(1), schur_rhs, + m_inverse); + } + catch (...) + { + abort (); + } + + // produce a consistent flow field + hanging_node_constraints.distribute (solution); + + std::cout << " " + << solver_control.last_step() + << " CG Schur complement iterations for pressure." + << std::endl; + } + + { + system_matrix.block(0,1).vmult (tmp, solution.block(1)); + tmp *= -1; + tmp += system_rhs.block(0); + + A_inverse.vmult (solution.block(0), tmp); + + // produce a consistent pressure field + hanging_node_constraints.distribute (solution); + } +} + + + +template +void +StokesProblem::output_results (const unsigned int refinement_cycle) const +{ + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("p"); + + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + + std::vector + data_component_interpretation + (dim+1, DataComponentInterpretation::component_is_scalar); + for (unsigned int i=0; i::type_dof_data, + data_component_interpretation); + + data_out.build_patches (); + + std::ostringstream filename; + filename << "solution-" + << Utilities::int_to_string (refinement_cycle, 2) + << ".vtk"; + + std::ofstream output (filename.str().c_str()); + data_out.write_vtk (output); +} + + + +template +void +StokesProblem::refine_mesh () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + std::vector component_mask (dim+1, false); + component_mask[dim] = true; + KellyErrorEstimator::estimate (dof_handler, + QGauss(degree+1), + typename FunctionMap::type(), + solution, + estimated_error_per_cell, + component_mask); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.0); + triangulation.execute_coarsening_and_refinement (); +} + + +template +void StokesProblem::run () +{ + switch (dim) + { + case 2: + { + std::vector subdivisions (dim, 1); + subdivisions[0] = 4; + + GridGenerator::subdivided_hyper_rectangle (triangulation, + subdivisions, + Point(-2,-1), + Point(2,0), + true); + + triangulation.refine_global (2); + + break; + } + + case 3: + { + GridGenerator::hyper_shell (triangulation, + Point(), 0.5, 1.0); + + static HyperShellBoundary boundary; + triangulation.set_boundary (0, boundary); + + triangulation.refine_global (2); + + break; + } + + default: + Assert (false, ExcNotImplemented()); + } + + + for (unsigned int refinement_cycle = 0; refinement_cycle<9; + ++refinement_cycle) + { + std::cout << "Refinement cycle " << refinement_cycle << std::endl; + + if (refinement_cycle > 0) + refine_mesh (); + + setup_dofs (); + + std::cout << " Assembling..." << std::endl; + assemble_system (); + + std::cout << " Solving..." << std::endl; + solve (); + + output_results (refinement_cycle); + + std::cout << std::endl; + } +} + + + +int main () +{ + try + { + deallog.depth_console (0); + + StokesProblem<2> flow_problem(1); + flow_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; +}