From aafabbfc233ebcdaa7db019b218d06317665dd30 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 18 Apr 2002 16:08:34 +0000 Subject: [PATCH] Check-in to allow nightly buildsto succeed. git-svn-id: https://svn.dealii.org/trunk@5678 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-14/Makefile | 159 ++++++ deal.II/examples/step-14/step-14.cc | 851 ++++++++++++++++++++++++++++ 2 files changed, 1010 insertions(+) create mode 100644 deal.II/examples/step-14/Makefile create mode 100644 deal.II/examples/step-14/step-14.cc diff --git a/deal.II/examples/step-14/Makefile b/deal.II/examples/step-14/Makefile new file mode 100644 index 0000000000..c02ef0f148 --- /dev/null +++ b/deal.II/examples/step-14/Makefile @@ -0,0 +1,159 @@ +# $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 + + +# 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. +%.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-14/step-14.cc b/deal.II/examples/step-14/step-14.cc new file mode 100644 index 0000000000..cc79c0902a --- /dev/null +++ b/deal.II/examples/step-14/step-14.cc @@ -0,0 +1,851 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 2001, 2002 */ + +/* $Id$ */ +/* Version: $Name$ */ +/* */ +/* Copyright (C) 2001, 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 + +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + + + +namespace Evaluation +{ + + template + class EvaluationBase + { + public: + virtual ~EvaluationBase (); + + void set_refinement_cycle (const unsigned int refinement_cycle); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const = 0; + protected: + unsigned int refinement_cycle; + }; + + + template + EvaluationBase::~EvaluationBase () + {}; + + + + template + void + EvaluationBase::set_refinement_cycle (const unsigned int step) + { + refinement_cycle = step; + }; + + + + template + class PointValueEvaluation : public EvaluationBase + { + public: + PointValueEvaluation (const Point &evaluation_point, + TableHandler &results_table); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + private: + const Point evaluation_point; + TableHandler &results_table; + }; + + + template + PointValueEvaluation:: + PointValueEvaluation (const Point &evaluation_point, + TableHandler &results_table) + : + evaluation_point (evaluation_point), + results_table (results_table) + {}; + + + + template + void + PointValueEvaluation:: + operator () (const DoFHandler &dof_handler, + const Vector &solution) const + { + double point_value = 1e20; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + bool evaluation_point_found = false; + for (; (cell!=endc) && !evaluation_point_found; ++cell) + for (unsigned int vertex=0; + vertex::vertices_per_cell; + ++vertex) + if (cell->vertex(vertex) == evaluation_point) + { + point_value = solution(cell->vertex_dof_index(vertex,0)); + + evaluation_point_found = true; + break; + }; + + AssertThrow (evaluation_point_found, + ExcEvaluationPointNotFound(evaluation_point)); + + results_table.add_value ("DoFs", dof_handler.n_dofs()); + results_table.add_value ("u(x_0)", point_value); + }; + + + + + + template + class SolutionOutput : public EvaluationBase + { + public: + SolutionOutput (const std::string &output_name_base, + const typename DataOut::OutputFormat output_format); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const; + private: + const std::string output_name_base; + const typename DataOut::OutputFormat output_format; + }; + + + template + SolutionOutput:: + SolutionOutput (const std::string &output_name_base, + const typename DataOut::OutputFormat output_format) + : + output_name_base (output_name_base), + output_format (output_format) + {}; + + + template + void + SolutionOutput::operator () (const DoFHandler &dof_handler, + const Vector &solution) const + { + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream filename; +#else + std::ostrstream filename; +#endif + filename << output_name_base << "-" + << refinement_cycle + << data_out.default_suffix (output_format) + << std::ends; +#ifdef HAVE_STD_STRINGSTREAM + std::ofstream out (filename.str().c_str()); +#else + std::ofstream out (filename.str()); +#endif + + data_out.write (out, output_format); + }; + + + + +}; + + + +namespace LaplaceSolver +{ + + template + class Base + { + public: + Base (Triangulation &coarse_grid); + virtual ~Base (); + + virtual void solve_problem () = 0; + virtual void postprocess (const Evaluation::EvaluationBase &postprocessor) const = 0; + virtual void refine_grid () = 0; + virtual unsigned int n_dofs () const = 0; + + protected: + const SmartPointer > triangulation; + }; + + + template + Base::Base (Triangulation &coarse_grid) + : + triangulation (&coarse_grid) + {}; + + + template + Base::~Base () + {}; + + + + template + class Solver : public virtual Base + { + public: + Solver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &boundary_values); + virtual + ~Solver (); + + virtual + void + solve_problem (); + + virtual + void + postprocess (const Evaluation::EvaluationBase &postprocessor) const; + + virtual + unsigned int + n_dofs () const; + + protected: + const SmartPointer > fe; + const SmartPointer > quadrature; + DoFHandler dof_handler; + Vector solution; + const SmartPointer > boundary_values; + + virtual void assemble_rhs (Vector &rhs) const = 0; + + private: + struct LinearSystem + { + LinearSystem (const DoFHandler &dof_handler); + + void solve (Vector &solution) const; + + ConstraintMatrix hanging_node_constraints; + SparsityPattern sparsity_pattern; + SparseMatrix matrix; + Vector rhs; + }; + + void + assemble_linear_system (LinearSystem &linear_system); + + void + assemble_matrix (LinearSystem &linear_system, + const typename DoFHandler::active_cell_iterator &begin_cell, + const typename DoFHandler::active_cell_iterator &end_cell, + Threads::ThreadMutex &mutex) const ; + }; + + + + template + Solver::Solver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &boundary_values) + : + Base (triangulation), + fe (&fe), + quadrature (&quadrature), + dof_handler (triangulation), + boundary_values (&boundary_values) + {}; + + + template + Solver::~Solver () + { + dof_handler.clear (); + }; + + + template + void + Solver::solve_problem () + { + dof_handler.distribute_dofs (*fe); + solution.reinit (dof_handler.n_dofs()); + + LinearSystem linear_system (dof_handler); + assemble_linear_system (linear_system); + linear_system.solve (solution); + }; + + + template + void + Solver:: + postprocess (const Evaluation::EvaluationBase &postprocessor) const + { + postprocessor (dof_handler, solution); + }; + + + template + unsigned int + Solver::n_dofs () const + { + return dof_handler.n_dofs(); + }; + + + template + void + Solver::assemble_linear_system (LinearSystem &linear_system) + { + typedef + typename DoFHandler::active_cell_iterator + active_cell_iterator; + + const unsigned int n_threads = multithread_info.n_default_threads; + std::vector > + thread_ranges + = Threads::split_range (dof_handler.begin_active (), + dof_handler.end (), + n_threads); + + Threads::ThreadMutex mutex; + Threads::ThreadManager thread_manager; + for (unsigned int thread=0; thread::assemble_matrix) + .collect_args (this, + linear_system, + thread_ranges[thread].first, + thread_ranges[thread].second, + mutex)); + + assemble_rhs (linear_system.rhs); + linear_system.hanging_node_constraints.condense (linear_system.rhs); + + std::map boundary_value_map; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + *boundary_values, + boundary_value_map); + + + thread_manager.wait (); + linear_system.hanging_node_constraints.condense (linear_system.matrix); + + MatrixTools::apply_boundary_values (boundary_value_map, + linear_system.matrix, + solution, + linear_system.rhs); + + }; + + + template + void + Solver::assemble_matrix (LinearSystem &linear_system, + const typename DoFHandler::active_cell_iterator &begin_cell, + const typename DoFHandler::active_cell_iterator &end_cell, + Threads::ThreadMutex &mutex) const + { + FEValues fe_values (*fe, *quadrature, + UpdateFlags(update_gradients | + update_JxW_values)); + + const unsigned int dofs_per_cell = fe->dofs_per_cell; + const unsigned int n_q_points = quadrature->n_quadrature_points; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + for (typename DoFHandler::active_cell_iterator cell=begin_cell; + cell!=end_cell; ++cell) + { + cell_matrix.clear (); + + fe_values.reinit (cell); + const std::vector > > + & shape_grads = fe_values.get_shape_grads(); + const std::vector + & JxW_values = fe_values.get_JxW_values(); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + mutex.acquire (); + for (unsigned int i=0; i + Solver::LinearSystem:: + LinearSystem (const DoFHandler &dof_handler) + { + hanging_node_constraints.clear (); + + void (*mhnc_p) (const DoFHandler &, + ConstraintMatrix &) + = &DoFTools::make_hanging_node_constraints; + + Threads::ThreadManager thread_manager; + Threads::spawn (thread_manager, + Threads::encapsulate (mhnc_p) + .collect_args (dof_handler, + hanging_node_constraints)); + + 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); + + thread_manager.wait (); + hanging_node_constraints.close (); + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + matrix.reinit (sparsity_pattern); + rhs.reinit (dof_handler.n_dofs()); + }; + + + + template + void + Solver::LinearSystem::solve (Vector &solution) const + { + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(matrix, 1.2); + + cg.solve (matrix, solution, rhs, preconditioner); + + hanging_node_constraints.distribute (solution); + }; + + + + + + template + class PrimalSolver : public Solver + { + public: + PrimalSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values); + protected: + const SmartPointer > rhs_function; + virtual void assemble_rhs (Vector &rhs) const; + }; + + + template + PrimalSolver:: + PrimalSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (triangulation), + Solver (triangulation, fe, + quadrature, boundary_values), + rhs_function (&rhs_function) + {}; + + + + template + void + PrimalSolver:: + assemble_rhs (Vector &rhs) const + { + FEValues fe_values (*fe, *quadrature, + UpdateFlags(update_values | + update_q_points | + update_JxW_values)); + + const unsigned int dofs_per_cell = fe->dofs_per_cell; + const unsigned int n_q_points = quadrature->n_quadrature_points; + + Vector cell_rhs (dofs_per_cell); + std::vector rhs_values (n_q_points); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_rhs.clear (); + + fe_values.reinit (cell); + const FullMatrix + & shape_values = fe_values.get_shape_values(); + const std::vector + & JxW_values = fe_values.get_JxW_values(); + const std::vector > + & q_points = fe_values.get_quadrature_points(); + + rhs_function->value_list (q_points, rhs_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i + class RefinementGlobal : public PrimalSolver + { + public: + RefinementGlobal (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values); + + virtual void refine_grid (); + }; + + + + template + RefinementGlobal:: + RefinementGlobal (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (coarse_grid), + PrimalSolver (coarse_grid, fe, quadrature, + rhs_function, boundary_values) + {}; + + + + template + void + RefinementGlobal::refine_grid () + { + triangulation->refine_global (1); + }; + + + + template + class RefinementKelly : public PrimalSolver + { + public: + RefinementKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values); + + virtual void refine_grid (); + }; + + + + template + RefinementKelly:: + RefinementKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (coarse_grid), + PrimalSolver (coarse_grid, fe, quadrature, + rhs_function, boundary_values) + {}; + + + + template + void + RefinementKelly::refine_grid () + { + Vector estimated_error_per_cell (triangulation->n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + QGauss3(), + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + GridRefinement::refine_and_coarsen_fixed_number (*triangulation, + estimated_error_per_cell, + 0.3, 0.03); + triangulation->execute_coarsening_and_refinement (); + }; + +}; + + + + + +template +class Solution : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +Solution::value (const Point &p, + const unsigned int /*component*/) const +{ + double q = p(0); + for (unsigned int i=1; i +class RightHandSide : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double q = p(0); + for (unsigned int i=1; i +void +run_simulation (LaplaceSolver::Base &solver, + const std::list *> &postprocessor_list) +{ + std::cout << "Refinement cycle: "; + + for (unsigned int step=0; true; ++step) + { + std::cout << step << " " << std::flush; + + solver.solve_problem (); + + for (typename std::list *>::const_iterator + i = postprocessor_list.begin(); + i != postprocessor_list.end(); ++i) + { + (*i)->set_refinement_cycle (step); + solver.postprocess (**i); + }; + + + if (solver.n_dofs() < 20000) + solver.refine_grid (); + else + break; + }; + + std::cout << std::endl; +}; + + + +template +void solve_problem (const std::string &solver_name) +{ + const std::string header = "Running tests with \"" + solver_name + + "\" refinement criterion:"; + std::cout << header << std::endl + << std::string (header.size(), '-') << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (2); + const FE_Q fe(1); + const QGauss4 quadrature; + const RightHandSide rhs_function; + const Solution boundary_values; + + LaplaceSolver::Base * solver = 0; + if (solver_name == "global") + solver = new LaplaceSolver::RefinementGlobal (triangulation, fe, + quadrature, + rhs_function, + boundary_values); + else if (solver_name == "kelly") + solver = new LaplaceSolver::RefinementKelly (triangulation, fe, + quadrature, + rhs_function, + boundary_values); + else + AssertThrow (false, ExcNotImplemented()); + + TableHandler results_table; + Evaluation::PointValueEvaluation + postprocessor1 (Point(0.5,0.5), results_table); + + Evaluation::SolutionOutput + postprocessor2 (std::string("solution-")+solver_name, + DataOut::gnuplot); + + std::list *> postprocessor_list; + postprocessor_list.push_back (&postprocessor1); + postprocessor_list.push_back (&postprocessor2); + + run_simulation (*solver, postprocessor_list); + + results_table.write_text (std::cout); + delete solver; + + std::cout << std::endl; +}; + + + +int main () +{ + try + { + deallog.depth_console (0); + + solve_problem<2> ("global"); + solve_problem<2> ("kelly"); + } + 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