From 6aa2a46df89b65d9d6a3f03ff323ea141bd60fbe Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 8 Nov 2002 23:49:28 +0000 Subject: [PATCH] Add step-15 in its initial stages. git-svn-id: https://svn.dealii.org/trunk@6745 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-15/Makefile | 150 ++++++++ deal.II/examples/step-15/step-15.cc | 513 ++++++++++++++++++++++++++++ 2 files changed, 663 insertions(+) create mode 100644 deal.II/examples/step-15/Makefile create mode 100644 deal.II/examples/step-15/step-15.cc diff --git a/deal.II/examples/step-15/Makefile b/deal.II/examples/step-15/Makefile new file mode 100644 index 0000000000..76278fb380 --- /dev/null +++ b/deal.II/examples/step-15/Makefile @@ -0,0 +1,150 @@ +# $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 = 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-1d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-1d.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 .go for object files compiled in debug mode and .o for +# object files in optimized mode. +ifeq ($(debug-mode),on) + libraries = $(target).g.o $(libs.g) +else + libraries = $(target).o $(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 $@ $^ $(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) + + +# 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. +./%.g.o : + @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-15/step-15.cc b/deal.II/examples/step-15/step-15.cc new file mode 100644 index 0000000000..00743ae483 --- /dev/null +++ b/deal.II/examples/step-15/step-15.cc @@ -0,0 +1,513 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 2002 */ + +/* $Id$ */ +/* Version: $Name$ */ +/* */ +/* Copyright (C) 2002 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 + +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + + + +template +class MinimizationProblem +{ + public: + MinimizationProblem (); + ~MinimizationProblem (); + void run (); + + private: + void setup_system (); + void assemble_newton_step (bool p); + void solve (); + void initialize (); + void refine_grid (); + double energy () const; + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + + FE_Q fe; + DoFHandler dof_handler; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix newton_matrix; + + Vector present_solution; + Vector newton_residual; + Vector newton_update; +}; + + + +class InitializationValues : public Function<1> +{ + public: + InitializationValues () : Function<1>() {}; + + virtual double value (const Point<1> &p, + const unsigned int component = 0) const; +}; + + + +double InitializationValues::value (const Point<1> &p, + const unsigned int) const +{ + return std::pow(p(0), 1./3.); +}; + + + +template +MinimizationProblem::MinimizationProblem () : + fe (1), + dof_handler (triangulation) +{}; + + +template +MinimizationProblem::~MinimizationProblem () +{ + dof_handler.clear (); +}; + + + +template +void MinimizationProblem::setup_system () +{ + 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(); +}; + + +template +double gradient_power (const Tensor<1,dim> &v, + const unsigned int n) +{ + Assert ((n/2)*2 == n, ExcMessage ("Value of 'n' must be even")); + double p = 1; + for (unsigned int k=0; k +void MinimizationProblem::assemble_newton_step (const bool p) +{ + if (p) + newton_matrix.reinit (sparsity_pattern); + newton_update.reinit (dof_handler.n_dofs()); + newton_residual.reinit (dof_handler.n_dofs()); + + QGauss3 quadrature_formula; + + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + std::vector local_solution_values (n_q_points); + std::vector > local_solution_grads (n_q_points); + + typename 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); + + fe_values.get_function_values (present_solution, + local_solution_values); + fe_values.get_function_grads (present_solution, + local_solution_grads); + + for (unsigned int q_point=0; q_point u_prime = local_solution_grads[q_point]; + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + if (dim == 1) + VectorTools::interpolate_boundary_values (dof_handler, + 1, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + newton_matrix, + newton_update, + newton_residual); +// std::cout << " res=" << newton_residual.l2_norm() << std::endl; +}; + + + +template +void MinimizationProblem::solve () +{ + SolverControl solver_control (1000, + 1e-3*newton_residual.l2_norm()); + PrimitiveVectorMemory<> vector_memory; + SolverGMRES<> gmres (solver_control, vector_memory); + + PreconditionJacobi<> preconditioner; + preconditioner.initialize(newton_matrix, 1.2); + +// std::cout << "------------------------------" << std::endl; +// newton_matrix.print_formatted(cout); +// std::cout << std::endl; +// for (unsigned int i=0; i +void MinimizationProblem::initialize () +{ + dof_handler.distribute_dofs (fe); + present_solution.reinit (dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, + InitializationValues(), + present_solution); +}; + + + +template +void MinimizationProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + typename FunctionMap::type neumann_boundary; + + KellyErrorEstimator::estimate (dof_handler, + Quadrature(1), + neumann_boundary, + present_solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + + SolutionTransfer solution_transfer(dof_handler); + triangulation.prepare_coarsening_and_refinement(); + solution_transfer.prepare_for_coarsening_and_refinement (present_solution); + triangulation.execute_coarsening_and_refinement (); + dof_handler.distribute_dofs (fe); + + Vector tmp (dof_handler.n_dofs()); + solution_transfer.interpolate (present_solution, tmp); + present_solution = tmp; +}; + + + +template +double MinimizationProblem::energy () const +{ + double energy = 0.; + + QGauss3 quadrature_formula; + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + std::vector local_solution_values (n_q_points); + std::vector > local_solution_grads (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (present_solution, + local_solution_values); + fe_values.get_function_grads (present_solution, + local_solution_grads); + + for (unsigned int q_point=0; q_point +void MinimizationProblem::output_results (const unsigned int cycle) const +{ + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (present_solution, "solution"); + data_out.build_patches (); + +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream filename; +#else + std::ostrstream filename; +#endif + filename << "solution-" + << cycle + << ".gnuplot" + << std::ends; +#ifdef HAVE_STD_STRINGSTREAM + std::ofstream out (filename.str().c_str()); +#else + std::ofstream out (filename.str()); +#endif + + data_out.write_gnuplot (out); +}; + + + +template +void MinimizationProblem::run () +{ + GridGenerator::hyper_cube (triangulation, 0., 1.); + triangulation.refine_global (1); + initialize (); + + for (unsigned int refinement_cycle=0; refinement_cycle<50; + ++refinement_cycle) + { + std::cout << "Cycle " << refinement_cycle << ':' << std::endl; + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + unsigned int iteration=0; + for (; iteration<1000; ++iteration) + { + assemble_newton_step (true); + solve (); + + if (newton_residual.l2_norm() < 1.e-12) + break; + }; + output_results (refinement_cycle); + + std::cout << " Iterations : " + << iteration + << std::endl; + std::cout << " Energy : " + << energy () + << std::endl; + + refine_grid (); + }; +}; + + +int main () +{ + try + { + deallog.depth_console (0); + + MinimizationProblem<1> minimization_problem_1d; + minimization_problem_1d.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; +}; -- 2.39.5