From: wolf Date: Wed, 26 Jan 2000 23:15:27 +0000 (+0000) Subject: Add step-7. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ffb1cd645640fa1954201e971faac64072943e52;p=dealii-svn.git Add step-7. git-svn-id: https://svn.dealii.org/trunk@2274 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 index aa281775ed..2f33358c23 100644 --- 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 @@ -811,6 +811,7 @@ void LaplaceProblem::output_results (const unsigned int cycle) const }; + template void LaplaceProblem::run () { diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-7/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-7/Makefile new file mode 100644 index 0000000000..be82c25ebd --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-7/Makefile @@ -0,0 +1,98 @@ +# $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.g = $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs = $(lib-deal2-2d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# 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-7/step-7.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc new file mode 100644 index 0000000000..c14b80d9ff --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc @@ -0,0 +1,486 @@ +/* $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 LaplaceProblem +{ + public: + enum RefinementMode { + global_refinement, adaptive_refinement + }; + + LaplaceProblem (const FiniteElement &fe, + const RefinementMode refinement_mode); + ~LaplaceProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void process_solution (const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + //... + SmartPointer > fe; + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + RefinementMode refinement_mode; +}; + + + +template +class SolutionBase +{ + protected: + static const unsigned int n_source_centers = 3; + static const Point source_centers[n_source_centers]; + static const double width; +}; + + +template +class Solution : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class RightHandSide : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template <> +const Point<1> +SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] += { Point<1>(-1.0 / 3.0), + Point<1>(0.0), + Point<1>(+1.0 / 3.0) }; + +template <> +const Point<2> +SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] += { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) }; + +template +const double SolutionBase::width = 0.15; + + + +template +double Solution::value (const Point &p, + const unsigned int) const +{ + double return_value = 0; + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + + return_value += exp(-shifted_point.square() / (width*width)); + }; + + return return_value; +}; + + + +template +Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const +{ + Tensor<1,dim> return_value; + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + + return_value += (-2 / (width*width) * + exp(-shifted_point.square() / (width*width)) * + shifted_point); + }; + + return return_value; +}; + + + +template +double RightHandSide::value (const Point &p, + const unsigned int) const +{ + double return_value = 0; + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + + return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / (width*width) * + exp(-shifted_point.square() / (width*width))); + }; + + return return_value; +}; + + + +template +LaplaceProblem::LaplaceProblem (const FiniteElement &fe, + const RefinementMode refinement_mode) : + dof_handler (triangulation), + fe (&fe), + refinement_mode (refinement_mode) +{}; + + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +}; + + + +template +void LaplaceProblem::setup_system () +{ + dof_handler.distribute_dofs (*fe); + // Renumber the degrees of freedom... + DoFRenumbering::Cuthill_McKee (dof_handler); + + 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 () +{ + 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; + + RightHandSide right_hand_side; + vector rhs_values (n_q_points); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + vector local_dof_indices (dofs_per_cell); + + 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(); + + right_hand_side.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 boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + Solution(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + 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 () +{ + switch (refinement_mode) + { + case global_refinement: + { + triangulation.refine_global (1); + break; + }; + + case adaptive_refinement: + { + 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 (); + + break; + }; + }; +}; + + +#include + +template +void LaplaceProblem::process_solution (const unsigned int cycle) const +{ + Vector difference_per_cell (triangulation.n_active_cells()); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss3(), + L2_norm); + const double L2_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss3(), + H1_seminorm); + const double H1_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss3(), + Linfty_norm); + const double Linfty_error = difference_per_cell.linfty_norm(); + + cout << "Cycle " << cycle << ':' + << endl + << " Number of active cells: " + << triangulation.n_active_cells() + << endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << endl; + + cout << " L2 error: " << L2_error << endl + << " H1 error: " << H1_error << endl + << " Linfty error: " << Linfty_error << endl; +}; + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<12; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (1); + } + else + refine_grid (); + + setup_system (); + + assemble_system (); + solve (); + process_solution (cycle); + }; + + string filename; + switch (refinement_mode) + { + case global_refinement: + filename = "solution-global"; + break; + case adaptive_refinement: + filename = "solution-adaptive"; + break; + default: + Assert (false, ExcInternalError()); + }; + filename += ".gmv"; + + ofstream output (filename.c_str()); + + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + data_out.write_gmv (output); +}; + + + +int main () +{ + try + { + deallog.depth_console (0); + + FEQ1<2> fe; + LaplaceProblem<2> laplace_problem_2d (fe, LaplaceProblem<2>::adaptive_refinement); + laplace_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; +}; diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/intro.html b/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/intro.html new file mode 100644 index 0000000000..e69de29bb2 diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/results.html b/deal.II/doc/tutorial/chapter-2.step-by-step/step-7.data/results.html new file mode 100644 index 0000000000..e69de29bb2 diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc index aa281775ed..2f33358c23 100644 --- a/deal.II/examples/step-6/step-6.cc +++ b/deal.II/examples/step-6/step-6.cc @@ -811,6 +811,7 @@ void LaplaceProblem::output_results (const unsigned int cycle) const }; + template void LaplaceProblem::run () { diff --git a/deal.II/examples/step-7/Makefile b/deal.II/examples/step-7/Makefile new file mode 100644 index 0000000000..be82c25ebd --- /dev/null +++ b/deal.II/examples/step-7/Makefile @@ -0,0 +1,98 @@ +# $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.g = $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs = $(lib-deal2-2d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# 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-7/step-7.cc b/deal.II/examples/step-7/step-7.cc new file mode 100644 index 0000000000..c14b80d9ff --- /dev/null +++ b/deal.II/examples/step-7/step-7.cc @@ -0,0 +1,486 @@ +/* $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 LaplaceProblem +{ + public: + enum RefinementMode { + global_refinement, adaptive_refinement + }; + + LaplaceProblem (const FiniteElement &fe, + const RefinementMode refinement_mode); + ~LaplaceProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void process_solution (const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + //... + SmartPointer > fe; + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + RefinementMode refinement_mode; +}; + + + +template +class SolutionBase +{ + protected: + static const unsigned int n_source_centers = 3; + static const Point source_centers[n_source_centers]; + static const double width; +}; + + +template +class Solution : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class RightHandSide : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template <> +const Point<1> +SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] += { Point<1>(-1.0 / 3.0), + Point<1>(0.0), + Point<1>(+1.0 / 3.0) }; + +template <> +const Point<2> +SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] += { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) }; + +template +const double SolutionBase::width = 0.15; + + + +template +double Solution::value (const Point &p, + const unsigned int) const +{ + double return_value = 0; + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + + return_value += exp(-shifted_point.square() / (width*width)); + }; + + return return_value; +}; + + + +template +Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const +{ + Tensor<1,dim> return_value; + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + + return_value += (-2 / (width*width) * + exp(-shifted_point.square() / (width*width)) * + shifted_point); + }; + + return return_value; +}; + + + +template +double RightHandSide::value (const Point &p, + const unsigned int) const +{ + double return_value = 0; + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + + return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / (width*width) * + exp(-shifted_point.square() / (width*width))); + }; + + return return_value; +}; + + + +template +LaplaceProblem::LaplaceProblem (const FiniteElement &fe, + const RefinementMode refinement_mode) : + dof_handler (triangulation), + fe (&fe), + refinement_mode (refinement_mode) +{}; + + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +}; + + + +template +void LaplaceProblem::setup_system () +{ + dof_handler.distribute_dofs (*fe); + // Renumber the degrees of freedom... + DoFRenumbering::Cuthill_McKee (dof_handler); + + 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 () +{ + 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; + + RightHandSide right_hand_side; + vector rhs_values (n_q_points); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + vector local_dof_indices (dofs_per_cell); + + 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(); + + right_hand_side.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 boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + Solution(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + 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 () +{ + switch (refinement_mode) + { + case global_refinement: + { + triangulation.refine_global (1); + break; + }; + + case adaptive_refinement: + { + 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 (); + + break; + }; + }; +}; + + +#include + +template +void LaplaceProblem::process_solution (const unsigned int cycle) const +{ + Vector difference_per_cell (triangulation.n_active_cells()); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss3(), + L2_norm); + const double L2_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss3(), + H1_seminorm); + const double H1_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss3(), + Linfty_norm); + const double Linfty_error = difference_per_cell.linfty_norm(); + + cout << "Cycle " << cycle << ':' + << endl + << " Number of active cells: " + << triangulation.n_active_cells() + << endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << endl; + + cout << " L2 error: " << L2_error << endl + << " H1 error: " << H1_error << endl + << " Linfty error: " << Linfty_error << endl; +}; + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<12; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (1); + } + else + refine_grid (); + + setup_system (); + + assemble_system (); + solve (); + process_solution (cycle); + }; + + string filename; + switch (refinement_mode) + { + case global_refinement: + filename = "solution-global"; + break; + case adaptive_refinement: + filename = "solution-adaptive"; + break; + default: + Assert (false, ExcInternalError()); + }; + filename += ".gmv"; + + ofstream output (filename.c_str()); + + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + data_out.write_gmv (output); +}; + + + +int main () +{ + try + { + deallog.depth_console (0); + + FEQ1<2> fe; + LaplaceProblem<2> laplace_problem_2d (fe, LaplaceProblem<2>::adaptive_refinement); + laplace_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; +};