From: Wolfgang Bangerth Date: Sun, 2 Jan 2000 20:34:42 +0000 (+0000) Subject: Add step-6. X-Git-Tag: v8.0.0~21219 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=262d62bb3c677de126d5deef479227e655285fd5;p=dealii.git Add step-6. git-svn-id: https://svn.dealii.org/trunk@2141 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile new file mode 100644 index 0000000000..dc2798f839 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile @@ -0,0 +1,94 @@ +# $Id$ +# Copyright W. Bangerth, University of Heidelberg, 1999 + +# Template for makefiles for the examples subdirectory. In principle, +# everything should be done automatically if you set the target file +# here correctly. We get deduce it from the files in the present +# directory: +target = $(basename $(shell echo step-*.cc)) + +# All dependencies between files should be updated by the included +# file Makefile.dep if necessary. Object files are compiled into +# the archives ./Obj.a and ./Obj.g.a. By default, the debug version +# is used to link. It you don't like that, change the following +# variable to "off" +debug-mode = on + + + +############################################################################### +# Internals + +#deal include base path +D = ../../../.. + +include $D/common/Make.global_options + +# get lists of files we need + + +# list of libraries needed to link with +libs = -ldeal_II_2d -llac -lbase +libs.g = -ldeal_II_2d.g -llac.g -lbase.g + + +# check whether we use debug mode or not +ifeq ($(debug-mode),on) + libraries = $(target).go $(libs.g) + flags = $(CXXFLAGS.g) +else + libraries = $(target).go $(libs) + flags = $(CXXFLAGS.o) +endif + + + +# make rule for the target. $^ is the object file $(target).g?o +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) $(flags) -o $@ $^ + +# rule how to run the program +run: $(target) + @echo ============================ Running $@ + @./$(target) + + +# rule to make object files +%.go : %.cc + @echo ============================ Compiling with debugging information: $< + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +%.o : %.cc + @echo ============================ Compiling with optimization: $< + @$(CXX) $(CXXFLAGS) -c $< -o $@ + + +clean: + -rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps + + + +.PHONY: clean + + +# Rule to generate the dependency file. 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. +# +# Since the script prefixes the output names by lib/g?o, we have to +# strip that again (the script was written for the main libraries and +# large projects where object files are put into subdirs) +Makefile.dep: $(target).cc Makefile \ + $(shell echo $D/base/include/base/*.h \ + $D/lac/include/lac/*.h \ + $D/deal.II/include/*/*.h) + @echo ============================ Remaking Makefile + @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \ + | perl -pi -e 's!lib/g?o/!!g;' \ + > Makefile.dep + + +include Makefile.dep + diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc new file mode 100644 index 0000000000..b87329d1dc --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc @@ -0,0 +1,349 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ + + + // still unfinished +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + // out statt in +#include + +#include + + //... +#include +#include + +#include + + + +template +class LaplaceProblem +{ + public: + LaplaceProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + FEQ1 fe; + DoFHandler dof_handler; + + ConstraintMatrix hanging_node_constraints; + + SparseMatrixStruct sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + +template +class Coefficient : 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 Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 20; + else + return 1; +}; + + + +template +void Coefficient::value_list (const vector > &points, + vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcVectorHasWrongSize (values.size(), n_points)); + + Assert (component == 0, + ExcWrongComponent (component, 1)); + + for (unsigned int i=0; i +LaplaceProblem::LaplaceProblem () : + dof_handler (triangulation) +{}; + + + +template +void LaplaceProblem::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 LaplaceProblem::assemble_system () +{ + const Coefficient coefficient; + + 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); + + vector local_dof_indices (dofs_per_cell); + + vector coefficient_values (n_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(); + + coefficient.value_list (q_points, coefficient_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); + + hanging_node_constraints.condense (system_matrix); + hanging_node_constraints.condense (system_rhs); +}; + + + +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionRelaxation<> + preconditioner(system_matrix, + &SparseMatrix::template precondition_SSOR, + 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +}; + + + //... +template +void LaplaceProblem::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 LaplaceProblem::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 LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) + { + cout << "Cycle " << cycle << ':' << endl; + + if (cycle == 0) + { + //... + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + triangulation.refine_global (1); + } + else + // ... + refine_grid (); + + cout << " Number of active cells: " + << triangulation.n_active_cells() + << endl; + + setup_system (); + 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.eps"); + data_out.write_eps (output); +}; + + + +int main () +{ + deallog.depth_console (0); + + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + + return 0; +}; diff --git a/deal.II/examples/step-6/Makefile b/deal.II/examples/step-6/Makefile new file mode 100644 index 0000000000..dc2798f839 --- /dev/null +++ b/deal.II/examples/step-6/Makefile @@ -0,0 +1,94 @@ +# $Id$ +# Copyright W. Bangerth, University of Heidelberg, 1999 + +# Template for makefiles for the examples subdirectory. In principle, +# everything should be done automatically if you set the target file +# here correctly. We get deduce it from the files in the present +# directory: +target = $(basename $(shell echo step-*.cc)) + +# All dependencies between files should be updated by the included +# file Makefile.dep if necessary. Object files are compiled into +# the archives ./Obj.a and ./Obj.g.a. By default, the debug version +# is used to link. It you don't like that, change the following +# variable to "off" +debug-mode = on + + + +############################################################################### +# Internals + +#deal include base path +D = ../../../.. + +include $D/common/Make.global_options + +# get lists of files we need + + +# list of libraries needed to link with +libs = -ldeal_II_2d -llac -lbase +libs.g = -ldeal_II_2d.g -llac.g -lbase.g + + +# check whether we use debug mode or not +ifeq ($(debug-mode),on) + libraries = $(target).go $(libs.g) + flags = $(CXXFLAGS.g) +else + libraries = $(target).go $(libs) + flags = $(CXXFLAGS.o) +endif + + + +# make rule for the target. $^ is the object file $(target).g?o +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) $(flags) -o $@ $^ + +# rule how to run the program +run: $(target) + @echo ============================ Running $@ + @./$(target) + + +# rule to make object files +%.go : %.cc + @echo ============================ Compiling with debugging information: $< + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +%.o : %.cc + @echo ============================ Compiling with optimization: $< + @$(CXX) $(CXXFLAGS) -c $< -o $@ + + +clean: + -rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps + + + +.PHONY: clean + + +# Rule to generate the dependency file. 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. +# +# Since the script prefixes the output names by lib/g?o, we have to +# strip that again (the script was written for the main libraries and +# large projects where object files are put into subdirs) +Makefile.dep: $(target).cc Makefile \ + $(shell echo $D/base/include/base/*.h \ + $D/lac/include/lac/*.h \ + $D/deal.II/include/*/*.h) + @echo ============================ Remaking Makefile + @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \ + | perl -pi -e 's!lib/g?o/!!g;' \ + > Makefile.dep + + +include Makefile.dep + diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc new file mode 100644 index 0000000000..b87329d1dc --- /dev/null +++ b/deal.II/examples/step-6/step-6.cc @@ -0,0 +1,349 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ + + + // still unfinished +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + // out statt in +#include + +#include + + //... +#include +#include + +#include + + + +template +class LaplaceProblem +{ + public: + LaplaceProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + FEQ1 fe; + DoFHandler dof_handler; + + ConstraintMatrix hanging_node_constraints; + + SparseMatrixStruct sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + +template +class Coefficient : 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 Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 20; + else + return 1; +}; + + + +template +void Coefficient::value_list (const vector > &points, + vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcVectorHasWrongSize (values.size(), n_points)); + + Assert (component == 0, + ExcWrongComponent (component, 1)); + + for (unsigned int i=0; i +LaplaceProblem::LaplaceProblem () : + dof_handler (triangulation) +{}; + + + +template +void LaplaceProblem::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 LaplaceProblem::assemble_system () +{ + const Coefficient coefficient; + + 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); + + vector local_dof_indices (dofs_per_cell); + + vector coefficient_values (n_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(); + + coefficient.value_list (q_points, coefficient_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); + + hanging_node_constraints.condense (system_matrix); + hanging_node_constraints.condense (system_rhs); +}; + + + +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionRelaxation<> + preconditioner(system_matrix, + &SparseMatrix::template precondition_SSOR, + 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +}; + + + //... +template +void LaplaceProblem::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 LaplaceProblem::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 LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) + { + cout << "Cycle " << cycle << ':' << endl; + + if (cycle == 0) + { + //... + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + triangulation.refine_global (1); + } + else + // ... + refine_grid (); + + cout << " Number of active cells: " + << triangulation.n_active_cells() + << endl; + + setup_system (); + 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.eps"); + data_out.write_eps (output); +}; + + + +int main () +{ + deallog.depth_console (0); + + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + + return 0; +};