From: bangerth Date: Fri, 18 Feb 2011 01:31:39 +0000 (+0000) Subject: Add a stripped down version of step-22 as a first step of a tutorial program that... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a780120b83a57e537a01581a32ba23229d00baff;p=dealii-svn.git Add a stripped down version of step-22 as a first step of a tutorial program that tests the FE_Nothing element. git-svn-id: https://svn.dealii.org/trunk@23385 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-46/Makefile b/deal.II/examples/step-46/Makefile new file mode 100644 index 0000000000..575312d3ad --- /dev/null +++ b/deal.II/examples/step-46/Makefile @@ -0,0 +1,143 @@ +# $Id$ + + +# 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-46/doc/builds-on b/deal.II/examples/step-46/doc/builds-on new file mode 100644 index 0000000000..3c3e873724 --- /dev/null +++ b/deal.II/examples/step-46/doc/builds-on @@ -0,0 +1 @@ +step-8 step-22 step-27 diff --git a/deal.II/examples/step-46/doc/intro.dox b/deal.II/examples/step-46/doc/intro.dox new file mode 100644 index 0000000000..1518606e43 --- /dev/null +++ b/deal.II/examples/step-46/doc/intro.dox @@ -0,0 +1,3 @@ + +

Introduction

+ diff --git a/deal.II/examples/step-46/doc/kind b/deal.II/examples/step-46/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/deal.II/examples/step-46/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/deal.II/examples/step-46/doc/results.dox b/deal.II/examples/step-46/doc/results.dox new file mode 100644 index 0000000000..b5eaba9377 --- /dev/null +++ b/deal.II/examples/step-46/doc/results.dox @@ -0,0 +1,2 @@ +

Results

+ diff --git a/deal.II/examples/step-46/doc/tooltip b/deal.II/examples/step-46/doc/tooltip new file mode 100644 index 0000000000..ecd9c06522 --- /dev/null +++ b/deal.II/examples/step-46/doc/tooltip @@ -0,0 +1 @@ +Coupling different models in different parts of the domain diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc new file mode 100644 index 0000000000..5ad37baf86 --- /dev/null +++ b/deal.II/examples/step-46/step-46.cc @@ -0,0 +1,477 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, Texas A&M University, 2011 */ + +/* $Id$ */ +/* */ +/* Copyright (C) 2011 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 + +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 constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + +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 +{ + Assert (component < this->n_components, + ExcIndexRange (component, 0, this->n_components)); + + 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 +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 () +{ + system_matrix.clear (); + + dof_handler.distribute_dofs (fe); + DoFRenumbering::Cuthill_McKee (dof_handler); + + std::vector block_component (dim+1,0); + block_component[dim] = 1; + DoFRenumbering::component_wise (dof_handler, block_component); + + { + constraints.clear (); + std::vector component_mask (dim+1, true); + component_mask[dim] = false; + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 1, + BoundaryValues(), + constraints, + component_mask); + } + + constraints.close (); + + std::vector dofs_per_block (2); + DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); + const unsigned int n_u = dofs_per_block[0], + n_p = dofs_per_block[1]; + + 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; + + { + CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(), + dof_handler.n_dofs()); + + DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); + sparsity_pattern.copy_from (csp); + } + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + + +template +void StokesProblem::assemble_system () +{ + system_matrix=0; + system_rhs=0; + + QGauss quadrature_formula(degree+2); + + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + const unsigned int n_q_points = 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 RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + std::vector > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + 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; + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, local_rhs, + local_dof_indices, + system_matrix, system_rhs); + } +} + + + + +template +void StokesProblem::solve () +{ + SparseDirectUMFPACK direct_solver; + direct_solver.initialize (system_matrix); + direct_solver.vmult (solution, system_rhs); +} + + + + +template +void +StokesProblem::output_results (const unsigned int refinement_cycle) const +{ + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, solution_names, + DataOut::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 () +{ + { + std::vector subdivisions (dim, 1); + subdivisions[0] = 4; + + const Point bottom_left = (dim == 2 ? + Point(-2,-1) : + Point(-2,0,-1)); + const Point top_right = (dim == 2 ? + Point(2,0) : + Point(2,1,0)); + + GridGenerator::subdivided_hyper_rectangle (triangulation, + subdivisions, + bottom_left, + top_right); + } + + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->center()[dim-1] == 0) + cell->face(f)->set_all_boundary_indicators(1); + + + triangulation.refine_global (4-dim); + + for (unsigned int refinement_cycle = 0; refinement_cycle<6; + ++refinement_cycle) + { + std::cout << "Refinement cycle " << refinement_cycle << std::endl; + + if (refinement_cycle > 0) + refine_mesh (); + + setup_dofs (); + + std::cout << " Assembling..." << std::endl << std::flush; + assemble_system (); + + std::cout << " Solving..." << std::flush; + 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; +}