From e4051c01a86a2fb184572573152f135df8a03949 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 16 Mar 2004 21:10:03 +0000 Subject: [PATCH] Add new example program. git-svn-id: https://svn.dealii.org/trunk@8794 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-17/.cvsignore | 4 + deal.II/examples/step-17/Makefile | 151 +++++++++ deal.II/examples/step-17/step-17.cc | 474 ++++++++++++++++++++++++++++ 3 files changed, 629 insertions(+) create mode 100644 deal.II/examples/step-17/.cvsignore create mode 100644 deal.II/examples/step-17/Makefile create mode 100644 deal.II/examples/step-17/step-17.cc diff --git a/deal.II/examples/step-17/.cvsignore b/deal.II/examples/step-17/.cvsignore new file mode 100644 index 0000000000..4847beaca4 --- /dev/null +++ b/deal.II/examples/step-17/.cvsignore @@ -0,0 +1,4 @@ +*.o *.go Makefile.dep *.gnuplot *.gmv *.eps +step-8 +*.ii +*.ti diff --git a/deal.II/examples/step-17/Makefile b/deal.II/examples/step-17/Makefile new file mode 100644 index 0000000000..a92d87c241 --- /dev/null +++ b/deal.II/examples/step-17/Makefile @@ -0,0 +1,151 @@ +# $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 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========= $( 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-17/step-17.cc b/deal.II/examples/step-17/step-17.cc new file mode 100644 index 0000000000..7e827bed64 --- /dev/null +++ b/deal.II/examples/step-17/step-17.cc @@ -0,0 +1,474 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ + +/* $Id$ */ +/* Version: $Name$ */ +/* */ +/* Copyright (C) 2004 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 + + // xxx +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + // xxx +#include + +#include +#include + + +template +class ElasticProblem +{ + public: + ElasticProblem (); + ~ElasticProblem (); + 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; + + FESystem fe; + + ConstraintMatrix hanging_node_constraints; + + // xxx no sparsity + PETScWrappers::SparseMatrix system_matrix; + + PETScWrappers::Vector solution; + PETScWrappers::Vector system_rhs; + + // xxx + const unsigned int n_partitions; + const unsigned int this_partition; +}; + + +template +class RightHandSide : public Function +{ + public: + RightHandSide (); + + virtual void vector_value (const Point &p, + Vector &values) const; + + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; +}; + + +template +RightHandSide::RightHandSide () : + Function (dim) +{} + + +template +inline +void RightHandSide::vector_value (const Point &p, + Vector &values) const +{ + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + Assert (dim >= 2, ExcInternalError()); + + Point point_1, point_2; + point_1(0) = 0.5; + point_2(0) = -0.5; + + if (((p-point_1).square() < 0.2*0.2) || + ((p-point_2).square() < 0.2*0.2)) + values(0) = 1; + else + values(0) = 0; + + if (p.square() < 0.2*0.2) + values(1) = 1; + else + values(1) = 0; +} + + + +template +void RightHandSide::vector_value_list (const std::vector > &points, + std::vector > &value_list) const +{ + const unsigned int n_points = points.size(); + + Assert (value_list.size() == n_points, + ExcDimensionMismatch (value_list.size(), n_points)); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); +} + + + + +template +ElasticProblem::ElasticProblem () + : + dof_handler (triangulation), + fe (FE_Q(1), dim), + // xxx + n_partitions (8), + this_partition (0) +{} + + + +template +ElasticProblem::~ElasticProblem () +{ + dof_handler.clear (); +} + + +template +void ElasticProblem::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 (); + + // no sparsity pattern + system_matrix.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + +template +void ElasticProblem::assemble_system () +{ + // move to front + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(dim), + boundary_values); + + QGauss2 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 lambda_values (n_q_points); + std::vector mu_values (n_q_points); + + ConstantFunction lambda(1.), mu(1.); + + RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim)); + + + 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); + + lambda.value_list (fe_values.get_quadrature_points(), lambda_values); + mu.value_list (fe_values.get_quadrature_points(), mu_values); + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + //xxx + MatrixTools::local_apply_boundary_values (boundary_values, + local_dof_indices, + cell_matrix, + cell_rhs, + false); + + // xxx + hanging_node_constraints + .distribute_local_to_global (cell_matrix, + local_dof_indices, + system_matrix); + + hanging_node_constraints + .distribute_local_to_global (cell_rhs, + local_dof_indices, + system_rhs); + } + + //xxx no condense necessary, no apply_b_v + //either + + + // xxx + system_matrix.compress (); + system_rhs.compress (); +} + + + +template +void ElasticProblem::solve () +{ + // xxx + SolverControl solver_control (1000, 1e-10); + PETScWrappers::SolverCG cg (solver_control); + + PETScWrappers::PreconditionSSOR preconditioner(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +} + + + +template +void ElasticProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + typename FunctionMap::type neumann_boundary; + KellyErrorEstimator::estimate (dof_handler, + QGauss2(), + neumann_boundary, + 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 (); + + // xxx + GridTools::partition_triangulation (n_partitions, triangulation); +} + + +template +void ElasticProblem::output_results (const unsigned int cycle) const +{ + std::string filename = "solution-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".gmv"; + std::ofstream output (filename.c_str()); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + + + + std::vector solution_names; + switch (dim) + { + case 1: + solution_names.push_back ("displacement"); + break; + case 2: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + break; + case 3: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + solution_names.push_back ("z_displacement"); + break; + default: + Assert (false, ExcInternalError()); + }; + + data_out.add_data_vector (solution, solution_names); + data_out.build_patches (); + data_out.write_gmv (output); +} + + + +template +void ElasticProblem::run () +{ + for (unsigned int cycle=0; cycle<12; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (2); + + // xxx + GridTools::partition_triangulation (n_partitions, triangulation); + } + else + refine_grid (); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + assemble_system (); + solve (); + output_results (cycle); + }; +} + + +int main (int argc, char **argv) +{ + try + { + // xxx + PetscInitialize(&argc,&argv,0,0); + deallog.depth_console (0); + + // xxx localize scope + { + ElasticProblem<2> elastic_problem_2d; + elastic_problem_2d.run (); + } + + // xxx + PetscFinalize(); + } + 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