From: Wolfgang Bangerth Date: Wed, 10 May 2000 13:45:15 +0000 (+0000) Subject: Add step-9, which is by no means complete, however. X-Git-Tag: v8.0.0~20573 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4c30f2ca643f1e3d5f8ea662eb9fe0e20d9b73cb;p=dealii.git Add step-9, which is by no means complete, however. git-svn-id: https://svn.dealii.org/trunk@2832 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-9/Makefile b/deal.II/examples/step-9/Makefile new file mode 100644 index 0000000000..e24b9541e7 --- /dev/null +++ b/deal.II/examples/step-9/Makefile @@ -0,0 +1,167 @@ +# $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 + + + + +# +# +# 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-lac.g) \ + $(lib-base.g) +libs.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 correct compiler flags and 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 .go for object files +# compiled in debug mode and .o for object files in optimized mode. +ifeq ($(debug-mode),on) + libraries = $(target).go $(libs.g) + flags = $(CXXFLAGS.g) +else + libraries = $(target).o $(libs.o) + flags = $(CXXFLAGS.o) +endif + + +# If in multithread mode, add the ACE library to the libraries which +# we need to link with: +ifneq ($(with-multithreading),no) + libraries += $(lib-ACE) +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) $(flags) -o $@ $^ + + +# 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) + + +# 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 *.o *.go *~ Makefile.dep $(target) $(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. +%.go : %.cc + @echo ==============debug========= $( Makefile.dep + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-9/step-9.cc b/deal.II/examples/step-9/step-9.cc new file mode 100644 index 0000000000..7de8c2edf7 --- /dev/null +++ b/deal.II/examples/step-9/step-9.cc @@ -0,0 +1,533 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ + +#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 + + +template +class AdvectionProblem +{ + public: + AdvectionProblem (); + ~AdvectionProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + + FEQ1 fe; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + +template +class AdvectionField +{ + public: + Point value (const Point &p) const; + + void value_list (const vector > &points, + vector > &values) const; + + /** + * Exception + */ + DeclException2 (ExcVectorHasWrongSize, + int, int, + << "The vector has size " << arg1 << " but should have " + << arg2 << " elements."); +}; + + + +template +Point +AdvectionField::value (const Point &p) const +{ + Point value; + value[0] = 2; + for (unsigned int i=1; i +void +AdvectionField::value_list (const vector > &points, + vector > &values) const +{ + Assert (values.size() == points.size(), + ExcVectorHasWrongSize (values.size(), points.size())); + + for (unsigned int i=0; i::value (points[i]); +}; + + + + +template +class RightHandSide : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const vector > &points, + vector &values, + const unsigned int component = 0) const; + + private: + static const Point center_point; +}; + + +template <> +const Point<1> RightHandSide<1>::center_point = Point<1> (-0.75); + +template <> +const Point<2> RightHandSide<2>::center_point = Point<2> (-0.75, -0.75); + +template <> +const Point<3> RightHandSide<3>::center_point = Point<3> (-0.75, -0.75, -0.75); + + + +template +double +RightHandSide::value (const Point &p, + const unsigned int component) const +{ + Assert (component == 0, ExcWrongComponent (component, 1)); + const double diameter = 0.1; + return ( (p-center_point).square() < diameter*diameter ? + .1/pow(diameter,dim) : + 0); +}; + + + +template +void +RightHandSide::value_list (const vector > &points, + vector &values, + const unsigned int component) const +{ + Assert (values.size() == points.size(), + ExcVectorHasWrongSize (values.size(), points.size())); + + for (unsigned int i=0; i::value (points[i], component); +}; + + + +template +class BoundaryValues : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const vector > &points, + vector &values, + const unsigned int component = 0) const; +}; + + + +template +double +BoundaryValues::value (const Point &p, + const unsigned int component) const +{ + Assert (component == 0, ExcWrongComponent (component, 1)); + + const double sine_term = sin(16*M_PI*sqrt(p.square())); + const double weight = exp(-5*p.square()) / exp(-5.); + return sine_term * weight; +}; + + + +template +void +BoundaryValues::value_list (const vector > &points, + vector &values, + const unsigned int component) const +{ + Assert (values.size() == points.size(), + ExcVectorHasWrongSize (values.size(), points.size())); + + for (unsigned int i=0; i::value (points[i], component); +}; + + + +template +AdvectionProblem::AdvectionProblem () : + dof_handler (triangulation) +{}; + + + +template +AdvectionProblem::~AdvectionProblem () +{ + dof_handler.clear (); +}; + + + +template +void AdvectionProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +}; + + + +template +void AdvectionProblem::assemble_system () +{ + AdvectionField advection_field; + RightHandSide right_hand_side; + BoundaryValues boundary_values; + + QGauss3 quadrature_formula; + QGauss3 face_quadrature_formula; + + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + FEFaceValues fe_face_values (fe, face_quadrature_formula, + UpdateFlags (update_values | + update_q_points | + update_JxW_values | + update_normal_vectors)); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + const unsigned int n_face_q_points = face_quadrature_formula.n_quadrature_points; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + vector local_dof_indices (dofs_per_cell); + + vector rhs_values (n_q_points); + vector > advection_directions (n_q_points); + vector face_boundary_values (n_face_q_points); + vector > face_advection_directions (n_face_q_points); + + DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix.clear (); + cell_rhs.clear (); + + fe_values.reinit (cell); + const FullMatrix + & shape_values = fe_values.get_shape_values(); + const vector > > + & shape_grads = fe_values.get_shape_grads(); + const vector + & JxW_values = fe_values.get_JxW_values(); + const vector > + & q_points = fe_values.get_quadrature_points(); + + advection_field.value_list (q_points, advection_directions); + right_hand_side.value_list (q_points, rhs_values); + + const double delta = 0.1 * cell->diameter (); + + for (unsigned int q_point=0; q_point::faces_per_cell; ++face) + { + fe_face_values.reinit (cell, face); + + const FullMatrix + & face_shape_values = fe_face_values.get_shape_values(); + const vector + & face_JxW_values = fe_face_values.get_JxW_values(); + const vector > + & face_q_points = fe_face_values.get_quadrature_points(); + const vector > + & normal_vectors = fe_face_values.get_normal_vectors(); + + boundary_values.value_list (face_q_points, face_boundary_values); + advection_field.value_list (face_q_points, face_advection_directions); + + for (unsigned int q_point=0; q_pointface(face)->at_boundary () && + (normal_vectors[q_point] * face_advection_directions[q_point] < 0)) + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i +void AdvectionProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverBicgstab<> bicgstab (solver_control, vector_memory); + + PreconditionRelaxation<> + preconditioner(system_matrix, + &SparseMatrix::template precondition_Jacobi, + 1.0); + + bicgstab.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +}; + + +template +void AdvectionProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::FunctionMap neumann_boundary; + + KellyErrorEstimator::estimate (dof_handler, + QGauss3(), + neumann_boundary, + solution, + estimated_error_per_cell); + + triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell, + 0.3, 0.03); + + triangulation.execute_coarsening_and_refinement (); +}; + + + +template +void AdvectionProblem::output_results (const unsigned int cycle) const +{ + string filename = "grid-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".eps"; + ofstream output (filename.c_str()); + + GridOut grid_out; + grid_out.write_eps (triangulation, output); +}; + + + +template +void AdvectionProblem::run () +{ + for (unsigned int cycle=0; cycle<6; ++cycle) + { + cout << "Cycle " << cycle << ':' << endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (4); + } + else + { + refine_grid (); + }; + + + cout << " Number of active cells: " + << triangulation.n_active_cells() + << endl; + + setup_system (); + + cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << endl; + + assemble_system (); + solve (); + output_results (cycle); + }; + + DataOut::EpsFlags eps_flags; + eps_flags.z_scaling = 4; + + DataOut data_out; + data_out.set_flags (eps_flags); + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + ofstream output ("final-solution.gmv"); + data_out.write_gmv (output); +}; + + + + // The ``main'' function is exactly + // like in previous examples, with + // the only difference in the name of + // the main class that actually does + // the computation. +int main () +{ + try + { + deallog.depth_console (0); + + AdvectionProblem<2> advection_problem_2d; + advection_problem_2d.run (); + } + catch (exception &exc) + { + cerr << endl << endl + << "----------------------------------------------------" + << endl; + cerr << "Exception on processing: " << endl + << exc.what() << endl + << "Aborting!" << endl + << "----------------------------------------------------" + << endl; + return 1; + } + catch (...) + { + cerr << endl << endl + << "----------------------------------------------------" + << endl; + cerr << "Unknown exception!" << endl + << "Aborting!" << endl + << "----------------------------------------------------" + << endl; + return 1; + }; + + return 0; +};