From 14f562fd0ed939edf7fe8d696d950b2c5c3c9188 Mon Sep 17 00:00:00 2001 From: guido Date: Wed, 26 Sep 2001 13:35:53 +0000 Subject: [PATCH] Big-tests removed git-svn-id: https://svn.dealii.org/trunk@5100 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/big-tests/Makefile | 44 - tests/big-tests/convergence/.cvsignore | 6 - tests/big-tests/convergence/Makefile | 165 ---- tests/big-tests/convergence/convergence.cc | 535 ----------- tests/big-tests/convergence/make_ps | 52 -- tests/big-tests/dof/.cvsignore | 5 - tests/big-tests/dof/Makefile | 168 ---- tests/big-tests/dof/dof_test.cc | 451 --------- tests/big-tests/dof/dof_test.prm | 4 - tests/big-tests/dof/make_ps | 122 --- tests/big-tests/dof/results/.cvsignore | 1 - tests/big-tests/error-estimation/.cvsignore | 4 - tests/big-tests/error-estimation/Makefile | 168 ---- tests/big-tests/error-estimation/ee.gauss.prm | 8 - tests/big-tests/error-estimation/ee.kink.prm | 8 - .../error-estimation/ee.singular.prm | 8 - .../error-estimation/error-estimation.cc | 759 ---------------- tests/big-tests/error-estimation/make_ps | 94 -- .../big-tests/error-estimation/strip_comments | 1 - tests/big-tests/grid/.cvsignore | 4 - tests/big-tests/grid/Makefile | 175 ---- tests/big-tests/grid/grid_test.cc | 329 ------- tests/big-tests/grid/make_ps | 43 - tests/big-tests/grid/results/.cvsignore | 1 - tests/big-tests/multigrid/.cvsignore | 4 - tests/big-tests/multigrid/Makefile | 165 ---- tests/big-tests/multigrid/make_ps | 52 -- tests/big-tests/multigrid/multigrid.cc | 515 ----------- .../fixed-point-iteration/.cvsignore | 3 - .../nonlinear/fixed-point-iteration/Makefile | 138 --- .../fixed-point-iteration/nonlinear.cc | 258 ------ tests/big-tests/poisson/.cvsignore | 6 - tests/big-tests/poisson/Makefile | 136 --- tests/big-tests/poisson/equation.cc | 86 -- tests/big-tests/poisson/poisson.cc | 31 - tests/big-tests/poisson/poisson.h | 100 -- tests/big-tests/poisson/poisson.prm | 5 - tests/big-tests/poisson/problem.cc | 622 ------------- tests/big-tests/poisson/results/.cvsignore | 2 - tests/big-tests/problem_base.h | 857 ------------------ 40 files changed, 6135 deletions(-) delete mode 100644 tests/big-tests/Makefile delete mode 100644 tests/big-tests/convergence/.cvsignore delete mode 100644 tests/big-tests/convergence/Makefile delete mode 100644 tests/big-tests/convergence/convergence.cc delete mode 100644 tests/big-tests/convergence/make_ps delete mode 100644 tests/big-tests/dof/.cvsignore delete mode 100644 tests/big-tests/dof/Makefile delete mode 100644 tests/big-tests/dof/dof_test.cc delete mode 100644 tests/big-tests/dof/dof_test.prm delete mode 100644 tests/big-tests/dof/make_ps delete mode 100644 tests/big-tests/dof/results/.cvsignore delete mode 100644 tests/big-tests/error-estimation/.cvsignore delete mode 100644 tests/big-tests/error-estimation/Makefile delete mode 100644 tests/big-tests/error-estimation/ee.gauss.prm delete mode 100644 tests/big-tests/error-estimation/ee.kink.prm delete mode 100644 tests/big-tests/error-estimation/ee.singular.prm delete mode 100644 tests/big-tests/error-estimation/error-estimation.cc delete mode 100644 tests/big-tests/error-estimation/make_ps delete mode 100755 tests/big-tests/error-estimation/strip_comments delete mode 100644 tests/big-tests/grid/.cvsignore delete mode 100644 tests/big-tests/grid/Makefile delete mode 100644 tests/big-tests/grid/grid_test.cc delete mode 100644 tests/big-tests/grid/make_ps delete mode 100644 tests/big-tests/grid/results/.cvsignore delete mode 100644 tests/big-tests/multigrid/.cvsignore delete mode 100644 tests/big-tests/multigrid/Makefile delete mode 100644 tests/big-tests/multigrid/make_ps delete mode 100644 tests/big-tests/multigrid/multigrid.cc delete mode 100644 tests/big-tests/nonlinear/fixed-point-iteration/.cvsignore delete mode 100644 tests/big-tests/nonlinear/fixed-point-iteration/Makefile delete mode 100644 tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc delete mode 100644 tests/big-tests/poisson/.cvsignore delete mode 100644 tests/big-tests/poisson/Makefile delete mode 100644 tests/big-tests/poisson/equation.cc delete mode 100644 tests/big-tests/poisson/poisson.cc delete mode 100644 tests/big-tests/poisson/poisson.h delete mode 100644 tests/big-tests/poisson/poisson.prm delete mode 100644 tests/big-tests/poisson/problem.cc delete mode 100644 tests/big-tests/poisson/results/.cvsignore delete mode 100644 tests/big-tests/problem_base.h diff --git a/tests/big-tests/Makefile b/tests/big-tests/Makefile deleted file mode 100644 index b3bcb7f5ee..0000000000 --- a/tests/big-tests/Makefile +++ /dev/null @@ -1,44 +0,0 @@ -# Generated automatically from Makefile.in by configure. -############################################################ -# $Id$ -# Copyright (C) 2000 by the deal.II authors -############################################################ - -############################################################ -# Include general settings for including DEAL libraries -############################################################ - -D = /home/atlas2/wolf/program/testdeal/deal.II/tests/.. - -include $D/common/Make.global_options - - -# list the directories we want to visit -subdirs = grid/ dof/ poisson/ convergence/ error-estimation/ multigrid/ nonlinear/fixed-point-iteration/ - -# define lists of targets: for each directory we produce a target name -# for compilation, running and cleaning by appending the action to -# the directory name (replacing the slash by ".action") -compile = $(subdirs:/=.compile) -run = $(subdirs:/=.run) -clean = $(subdirs:/=.clean) - -# define global targets which are to be excuted in every subdirectory -compile: $(compile) -run : $(run) -# for cleaning up: do this also for the present directory -clean : $(clean) - -rm -f *~ - - -# define the action of the targets for the specific subdirectories -$(compile) : - cd $(@:.compile=) ; $(MAKE) D=$D - -$(run) : - cd $(@:.run=) ; $(MAKE) run D=$D - -$(clean) : - -cd $(@:.clean=) ; $(MAKE) clean D=$D - - diff --git a/tests/big-tests/convergence/.cvsignore b/tests/big-tests/convergence/.cvsignore deleted file mode 100644 index 0cecd7ae85..0000000000 --- a/tests/big-tests/convergence/.cvsignore +++ /dev/null @@ -1,6 +0,0 @@ -convergence -Makefile.dep -*.go -*.o -*.ucd -*.history diff --git a/tests/big-tests/convergence/Makefile b/tests/big-tests/convergence/Makefile deleted file mode 100644 index 90a506661a..0000000000 --- a/tests/big-tests/convergence/Makefile +++ /dev/null @@ -1,165 +0,0 @@ -# $Id$ -# Copyright by the deal.II authors 1998, 1999, 2000 - -############################################################ -# $D contains the root of the deal distribution and should -# be supplied as a command line argument D=../../... -############################################################ - - -# 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 = convergence - -# 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 - - -# 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 *ucd *history *~ - - - - -# -# -# 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).go $(libs.o) - flags = $(CXXFLAGS.o) -endif - - -# If in multithread mode, add the ACE library to the libraries which -# we need to link with: -ifneq ($(with-multithreading),no) - libraries += $(lib-ACE) -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) $(LDFLAGS) -o $@ $^ $(LIBS) - - -# 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) - gnuplot make_ps - - -# 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/tests/big-tests/convergence/convergence.cc b/tests/big-tests/convergence/convergence.cc deleted file mode 100644 index 67b9a3dab8..0000000000 --- a/tests/big-tests/convergence/convergence.cc +++ /dev/null @@ -1,535 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "../problem_base.h" -#include -#include - -#include -#include -#include -#include -#include - - - - - -template -class PoissonEquation : public Equation { - public: - PoissonEquation (const Function &rhs) : - Equation(1), - right_hand_side (rhs) {}; - - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - protected: - const Function &right_hand_side; -}; - - - - - - -template -class PoissonProblem : public ProblemBase { - public: - PoissonProblem (unsigned int order); - ~PoissonProblem (); - - void clear (); - void create_new (); - int run (unsigned int level); - void print_history (string filename) const; - - protected: - Triangulation *tria; - DoFHandler *dof; - - Function *rhs; - Function *boundary_values; - - vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; - vector n_dofs; - - unsigned int order; -}; - - - - - -/** - Right hand side constructed such that the exact solution is - $sin(2 pi x) + sin(2 pi y)$ - */ -template -class RHSPoly : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const; -}; - - - -template -class Solution : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const; - /** - * Return the gradient of the function - * at the given point. - */ - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component) const; -}; - - - - -template <> -double RHSPoly<2>::value (const Point<2> &p, - const unsigned int) const { - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y)); -}; - - - -template <> -double Solution<2>::value (const Point<2> &p, - const unsigned int) const { - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return sin(2*pi*x)+sin(2*pi*y); -}; - - -template <> -Tensor<1,2> Solution<2>::gradient (const Point<2> &p, - const unsigned int) const { - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return Point<2> (2*pi*cos(2*pi*x), - 2*pi*cos(2*pi*y)); -}; - - - - - -template <> -void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues<2> &fe_values, - const DoFHandler<2>::cell_iterator &) const { - for (unsigned int point=0; point -void PoissonEquation::assemble (FullMatrix &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, ExcPureFunctionCalled()); -}; - - - -template -void PoissonEquation::assemble (Vector &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, ExcPureFunctionCalled()); -}; - - - - - - - - - -template -PoissonProblem::PoissonProblem (unsigned int order) : - tria(0), dof(0), rhs(0), - boundary_values(0), order(order) {}; - - -template -PoissonProblem::~PoissonProblem () -{ - clear (); -}; - - - -template -void PoissonProblem::clear () { - if (dof != 0) { - delete dof; - dof = 0; - }; - - if (tria != 0) { - delete tria; - tria = 0; - }; - - - // make it known to the underlying - // ProblemBase that tria and dof - // are already deleted - set_tria_and_dof (tria, dof); - - if (rhs != 0) - { - delete rhs; - rhs = 0; - }; - - if (boundary_values != 0) - { - delete boundary_values; - boundary_values = 0; - }; - - ProblemBase::clear (); -}; - - - - -template -void PoissonProblem::create_new () { - clear (); - - tria = new Triangulation(); - dof = new DoFHandler (*tria); - set_tria_and_dof (tria, dof); -}; - - - - -template -int PoissonProblem::run (const unsigned int level) { - create_new (); - - cout << "Refinement level = " << level - << ", using elements of type <"; - switch (order) - { - default: - cout << "Lagrange-" << order; - break; - }; - cout << ">" << endl; - - cout << " Making grid... "; - GridGenerator::hyper_ball (*tria); - HyperBallBoundary boundary_description; - tria->set_boundary (0, boundary_description); - tria->begin_active()->set_refine_flag(); - (++(++(tria->begin_active())))->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - tria->refine_global (level); - cout << tria->n_active_cells() << " active cells." << endl; - - rhs = new RHSPoly(); - boundary_values = new Solution (); - - - FiniteElement *fe; - PoissonEquation equation (*rhs); - Quadrature *quadrature; - Quadrature *boundary_quadrature; - switch (order) { - case 1: - fe = new FEQ1(); - quadrature = new QGauss3(); - boundary_quadrature = new QGauss2(); - break; - case 2: - fe = new FEQ2(); - quadrature = new QGauss4(); - boundary_quadrature = new QGauss3(); - break; - case 3: - fe = new FEQ3(); - quadrature = new QGauss5(); - boundary_quadrature = new QGauss4(); - break; - case 4: - fe = new FEQ4(); - quadrature = new QGauss6(); - boundary_quadrature = new QGauss5(); - break; - default: - return 100000; - }; - - cout << " Distributing dofs... "; - dof->distribute_dofs (*fe); - cout << dof->n_dofs() << " degrees of freedom." << endl; - n_dofs.push_back (dof->n_dofs()); - - cout << " Assembling matrices..." << endl; - UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | - update_gradients | update_JxW_values); - - ProblemBase::FunctionMap dirichlet_bc; - dirichlet_bc[0] = boundary_values; - assemble (equation, *quadrature, update_flags, dirichlet_bc); - - cout << " Solving..." << endl; - solve (); - - Solution sol; - Vector l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell; - Vector h1_seminorm_error_per_cell, h1_error_per_cell; - - cout << " Calculating L1 error... "; - VectorTools::integrate_difference (*dof_handler, - solution, sol, - l1_error_per_cell, - *quadrature, L1_norm); - cout << l1_error_per_cell.l1_norm() << endl; - l1_error.push_back (l1_error_per_cell.l1_norm()); - - cout << " Calculating L2 error... "; - VectorTools::integrate_difference (*dof_handler, - solution, sol, - l2_error_per_cell, - *quadrature, L2_norm); - cout << l2_error_per_cell.l2_norm() << endl; - l2_error.push_back (l2_error_per_cell.l2_norm()); - - cout << " Calculating L-infinity error... "; - VectorTools::integrate_difference (*dof_handler, - solution, sol, - linfty_error_per_cell, - *quadrature, Linfty_norm); - cout << linfty_error_per_cell.linfty_norm() << endl; - linfty_error.push_back (linfty_error_per_cell.linfty_norm()); - - cout << " Calculating H1-seminorm error... "; - VectorTools::integrate_difference (*dof_handler, - solution, sol, - h1_seminorm_error_per_cell, - *quadrature, H1_seminorm); - cout << h1_seminorm_error_per_cell.l2_norm() << endl; - h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); - - cout << " Calculating H1 error... "; - VectorTools::integrate_difference (*dof_handler, - solution, sol, - h1_error_per_cell, - *quadrature, H1_norm); - cout << h1_error_per_cell.l2_norm() << endl; - h1_error.push_back (h1_error_per_cell.l2_norm()); - - if (dof->n_dofs()<=5000) - { - Vector l1_error_per_dof(dof->n_dofs()); - Vector l2_error_per_dof(dof->n_dofs()); - Vector linfty_error_per_dof(dof->n_dofs()); - Vector h1_seminorm_error_per_dof(dof->n_dofs()); - Vector h1_error_per_dof(dof->n_dofs()); - DoFTools::distribute_cell_to_dof_vector (*dof, l1_error_per_cell, l1_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, l2_error_per_cell, l2_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, linfty_error_per_cell, - linfty_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, h1_seminorm_error_per_cell, - h1_seminorm_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, h1_error_per_cell, h1_error_per_dof); - -// Vector projected_solution; -// ConstraintMatrix constraints; -// constraints.close (); -// VectorTools::project (*dof, constraints, *fe, -// StraightBoundary(), *quadrature, -// sol, projected_solution, false, -// *boundary_quadrature); -// cout << " Calculating L2 error of projected solution... "; -// VectorTools::integrate_difference (*dof_handler, -// projected_solution, sol, -// l2_error_per_cell, -// *quadrature, *fe, L2_norm); -// cout << l2_error_per_cell.l2_norm() << endl; - - - string filename; - filename = ('0'+order); - filename += "."; - filename += ('0'+level); - filename += ".ucd"; - cout << " Writing error plots to <" << filename << ">..." << endl; - - DataOut out; - ofstream o(filename.c_str()); - fill_data (out); - out.add_data_vector (l1_error_per_dof, "L1_Error"); - out.add_data_vector (l2_error_per_dof, "L2_Error"); - out.add_data_vector (linfty_error_per_dof, "Linfty_Error"); - out.add_data_vector (h1_seminorm_error_per_dof, "H1_seminorm_Error"); - out.add_data_vector (h1_error_per_dof, "H1_Error"); - out.build_patches (); - out.write_ucd (o); - o.close (); - } - else - cout << " Not writing error as grid." << endl; - - cout << endl; - - const unsigned int n_dofs = dof->n_dofs(); - // release the lock that the dof object - // has to the finite element object - dof->clear (); - tria->set_boundary (0); - - delete fe; - delete quadrature; - delete boundary_quadrature; - - return n_dofs; -}; - - -template -void PoissonProblem::print_history (string filename) const { - ofstream out(filename.c_str()); - out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error" - << endl; - for (unsigned int i=0; ih/2:" << endl; - cout << " L1 error : " << 1./average_l1 << endl - << " L2 error : " << 1./average_l2 << endl - << " Linfty error : " << 1./average_linfty << endl - << " H1 seminorm error: " << 1./average_h1_semi << endl - << " H1 error : " << 1./average_h1 << endl; - cout << "==========================================================\n"; - cout << "==========================================================\n"; -}; - - - - -int main () { - deallog.depth_console (0); - for (unsigned int order=1; order<5; ++order) - { - PoissonProblem<2> problem (order); - - unsigned int level=0; - unsigned int n_dofs; - do - n_dofs = problem.run (level++); - while (n_dofs<25000); - - string filename; - switch (order) - { - case 1: - filename = "linear"; - break; - case 2: - filename = "quadratic"; - break; - case 3: - filename = "cubic"; - break; - case 4: - filename = "quartic"; - break; - }; - filename += ".history"; - - cout << endl << "Printing convergence history to <" - << filename << ">..." << endl; - problem.print_history (filename); - cout << endl << endl << endl; - }; - - return 0; -}; diff --git a/tests/big-tests/convergence/make_ps b/tests/big-tests/convergence/make_ps deleted file mode 100644 index 76c13a8624..0000000000 --- a/tests/big-tests/convergence/make_ps +++ /dev/null @@ -1,52 +0,0 @@ -set term postscript eps -set xlabel "Number of degrees of freedom" -set data style linespoints -set logscale xy - - - -set ylabel "Error" - -set output "criss-cross.eps" - -plot "criss_cross.history" using 1:2 title "L1 error","criss_cross.history" using 1:3 title "L2 error","criss_cross.history" using 1:4 title "Linfty error","criss_cross.history" using 1:5 title "H1 seminorm error","criss_cross.history" using 1:6 title "H1 error" - - - -set output "linear.eps" - -plot "linear.history" using 1:2 title "L1 error","linear.history" using 1:3 title "L2 error","linear.history" using 1:4 title "Linfty error","linear.history" using 1:5 title "H1 seminorm error","linear.history" using 1:6 title "H1 error" - - - -set output "quadratic.eps" - -plot "quadratic.history" using 1:2 title "L1 error","quadratic.history" using 1:3 title "L2 error","quadratic.history" using 1:4 title "Linfty error","quadratic.history" using 1:5 title "H1 seminorm error","quadratic.history" using 1:6 title "H1 error" - - - -set output "cubic.eps" - -plot "cubic.history" using 1:2 title "L1 error","cubic.history" using 1:3 title "L2 error","cubic.history" using 1:4 title "Linfty error","cubic.history" using 1:5 title "H1 seminorm error","cubic.history" using 1:6 title "H1 error" - - - -set output "quartic.eps" - -plot "quartic.history" using 1:2 title "L1 error","quartic.history" using 1:3 title "L2 error","quartic.history" using 1:4 title "Linfty error","quartic.history" using 1:5 title "H1 seminorm error","quartic.history" using 1:6 title "H1 error" - - - -set output "l2error.eps" -set ylabel "L2-error" - -plot "criss_cross.history" using 1:3 title "Criss-cross elements", "linear.history" using 1:3 title "Linear elements", "quadratic.history" using 1:3 title "Quadratic elements", "cubic.history" using 1:3 title "Cubic elements", "quartic.history" using 1:3 title "Quartic elements" - - - -set output "h1error.eps" -set ylabel "H1-error" - -plot "criss_cross.history" using 1:6 title "Criss-cross elements", "linear.history" using 1:6 title "Linear elements", "quadratic.history" using 1:6 title "Quadratic elements", "cubic.history" using 1:6 title "Cubic elements", "quartic.history" using 1:6 title "Quartic elements" - - diff --git a/tests/big-tests/dof/.cvsignore b/tests/big-tests/dof/.cvsignore deleted file mode 100644 index 6105d89a5b..0000000000 --- a/tests/big-tests/dof/.cvsignore +++ /dev/null @@ -1,5 +0,0 @@ -dof_test -Makefile.dep -Makefile.dep -*.go -*.o diff --git a/tests/big-tests/dof/Makefile b/tests/big-tests/dof/Makefile deleted file mode 100644 index 1ddc3b6d21..0000000000 --- a/tests/big-tests/dof/Makefile +++ /dev/null @@ -1,168 +0,0 @@ -# $Id$ -# Copyright by the deal.II authors 1998, 1999, 2000 - -############################################################ -# $D contains the root of the deal distribution and should -# be supplied as a command line argument D=../../... -############################################################ - - -# 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 = dof_test - -# 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 - - -# 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-deal2-3d.g) \ - $(lib-lac.g) \ - $(lib-base.g) -libs.o = $(lib-deal2-2d.o) \ - $(lib-deal2-3d.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).go $(libs.o) - flags = $(CXXFLAGS.o) -endif - - -# If in multithread mode, add the ACE library to the libraries which -# we need to link with: -ifneq ($(with-multithreading),no) - libraries += $(lib-ACE) -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) $(LDFLAGS) -o $@ $^ $(LIBS) - - -# 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) 2 $(target).prm - ./$(target) 3 $(target).prm - gnuplot make_ps - - -# 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/tests/big-tests/dof/dof_test.cc b/tests/big-tests/dof/dof_test.cc deleted file mode 100644 index a9513d4b60..0000000000 --- a/tests/big-tests/dof/dof_test.cc +++ /dev/null @@ -1,451 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - - - -// 1: continuous refinement of the unit square always in the middle -// 2: refinement of the circle at the boundary -// 2: refinement of a wiggled area at the boundary -// 4: random refinement - - - - - - -template -class Ball : - public StraightBoundary { - public: - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const { - Point middle = StraightBoundary::get_new_point_on_line(line); - - for (int i=0; i - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { - Point middle = StraightBoundary::get_new_point_on_quad(quad); - - for (int i=0; i -class CurvedLine : - public StraightBoundary { - public: - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - virtual Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; -}; - - - -template -Point -CurvedLine::get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - Point middle = StraightBoundary::get_new_point_on_line (line); - - // if the line is at the top of bottom - // face: do a special treatment on - // this line. Note that if the - // z-value of the midpoint is either - // 0 or 1, then the z-values of all - // vertices of the line is like that - if (dim>=3) - if (((middle(2) == 0) || (middle(2) == 1)) - // find out, if the line is in the - // interior of the top or bottom face - // of the domain, or at the edge. - // lines at the edge need to undergo - // the usual treatment, while for - // interior lines taking the midpoint - // is sufficient - // - // note: the trick with the boundary - // id was invented after the above was - // written, so we are not very strict - // here with using these flags - && (line->boundary_indicator() == 1)) - return middle; - - - double x=middle(0), - y=middle(1); - - if (y -Point -CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const -{ - Point middle = StraightBoundary::get_new_point_on_quad (quad); - - // if the face is at the top of bottom - // face: do not move the midpoint in - // x/y direction. Note that if the - // z-value of the midpoint is either - // 0 or 1, then the z-values of all - // vertices of the quad is like that - if ((middle(2) == 0) || (middle(2) == 1)) - return middle; - - double x=middle(0), - y=middle(1); - - if (y -class TestCases : public MultipleParameterLoop::UserClass{ - public: - TestCases (); - virtual ~TestCases (); - - virtual void create_new (const unsigned int run_no); - virtual void declare_parameters (ParameterHandler &prm); - virtual void run (ParameterHandler &prm); - - private: - Triangulation *tria; - DoFHandler *dof; -}; - - - -template -TestCases::TestCases () : - tria(0), dof(0) {}; - - -template -TestCases::~TestCases () -{ - if (dof) delete dof; - if (tria) delete tria; -}; - - - -template -void TestCases::create_new (const unsigned int) { - if (dof != 0) delete dof; - if (tria != 0) delete tria; - - tria = new Triangulation(); - GridGenerator::hyper_cube(*tria); - - dof = new DoFHandler (*tria); -}; - - - -template -void TestCases::declare_parameters (ParameterHandler &prm) { - if (dim>=2) - prm.declare_entry ("Test run", "zoom in", - Patterns::Selection("zoom in|ball|curved line|random")); - else - prm.declare_entry ("Test run", "zoom in", - Patterns::Selection("zoom in|random")); - prm.declare_entry ("Grid file", "grid.1"); - prm.declare_entry ("Sparsity file", "sparsity.1"); - prm.declare_entry ("Condensed sparsity file", "sparsity.c.1"); -}; - - - -template -void TestCases::run (ParameterHandler &prm) { - cout << "Dimension = " << dim - << ", Test case = " << prm.get ("Test run") << endl - << endl; - - string test = prm.get ("Test run"); - unsigned int test_case = 1; - if (test=="zoom in") test_case = 1; - else - if (test=="ball") test_case = 2; - else - if (test=="curved line") test_case = 3; - else - if (test=="random") test_case = 4; - else - cerr << "This test seems not to be implemented!" << endl; - - - cout << " Making grid..." << endl; - Boundary *boundary = 0; - - switch (test_case) - { - case 1: - { - // refine first cell - tria->begin_active()->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - // refine first active cell - // on coarsest level - tria->begin_active()->set_refine_flag (); - tria->execute_coarsening_and_refinement (); - - Triangulation::active_cell_iterator cell; - for (int i=0; i<17; ++i) - { - // refine the presently - // second last cell 17 - // times - cell = tria->last_active(tria->n_levels()-1); - --cell; - cell->set_refine_flag (); - tria->execute_coarsening_and_refinement (); - }; - - break; - } - - case 2: - case 3: - { - if (dim==3) - { - tria->begin_active()->face(2)->set_boundary_indicator(1); - tria->begin_active()->face(4)->set_boundary_indicator(1); - }; - - // set the boundary function - boundary = (test_case==2 ? - static_cast*>(new Ball()) : - static_cast*>(new CurvedLine())); - tria->set_boundary (0, *boundary); - tria->set_boundary (1, *boundary); - - // refine once - tria->begin_active()->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - - Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<6-dim; ++i) - { - cell = tria->begin_active(); - endc = tria->end(); - - // refine all - // boundary cells - for (; cell!=endc; ++cell) - if (cell->at_boundary()) - cell->set_refine_flag(); - - tria->execute_coarsening_and_refinement(); - }; - - break; - } - - case 4: - { - // refine once - tria->begin_active()->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - - Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<(dim==2 ? 12 : (dim==3 ? 7 : 20)); ++i) - { - int n_levels = tria->n_levels(); - cell = tria->begin_active(); - endc = tria->end(); - - for (; cell!=endc; ++cell) - { - double r = rand()*1.0/RAND_MAX, - weight = 1.* - (cell->level()*cell->level()) / - (n_levels*n_levels); - - if (r <= 0.5*weight) - cell->set_refine_flag (); - }; - - tria->execute_coarsening_and_refinement (); - }; - break; - } - }; - - // output the grid - string file_prefix ("results/"); - file_prefix += ('0'+dim); - file_prefix += "d."; - - cout << " Writing grid..." << endl; - ofstream out((file_prefix + prm.get("Grid file")).c_str()); - GridOut().write_gnuplot (*tria, out); - - - - - cout << " Distributing degrees of freedom..." << endl; - FEQ1 fe; - dof->distribute_dofs (fe); - - cout << " Renumbering degrees of freedom..." << endl; - DoFRenumbering::Cuthill_McKee (*dof); - - SparsityPattern sparsity (dof->n_dofs(), - dof->max_couplings_between_dofs()); - - - DoFTools::make_sparsity_pattern (*dof, sparsity); - int unconstrained_bandwidth = sparsity.bandwidth(); - - cout << " Writing sparsity pattern..." << endl; - ofstream sparsity_out ((file_prefix + prm.get("Sparsity file")).c_str()); - sparsity.print_gnuplot (sparsity_out); - - - - // computing constraints - cout << " Computing constraints..." << endl; - ConstraintMatrix constraints; - DoFTools::make_hanging_node_constraints (*dof, constraints); - constraints.close (); - constraints.condense (sparsity); - - cout << " Writing condensed sparsity pattern..." << endl; - ofstream c_sparsity_out ((file_prefix + - prm.get("Condensed sparsity file")).c_str()); - sparsity.print_gnuplot (c_sparsity_out); - - - cout << endl - << " Total number of cells = " << tria->n_cells() << endl - << " Total number of active cells = " << tria->n_active_cells() << endl - << " Number of DoFs = " << dof->n_dofs() << endl - << " Number of constraints = " << constraints.n_constraints() << endl - << " Unconstrained matrix bandwidth= " << unconstrained_bandwidth << endl - << " Constrained matrix bandwidth = " << sparsity.bandwidth() - << endl << endl; - - // release the lock that dof has to the - // finite element object - dof->clear (); - tria->set_boundary (0); - tria->set_boundary (1); - if (boundary) - delete boundary; -}; - - - -int main (int argc, char **argv) { - if (argc!=3) - { - cerr << "Usage: dof_test dimension parameterfile" << endl << endl; - return 1; - }; - - unsigned int dim; - if (argv[1][0] == '2') - dim = 2; - else - dim = 3; - - switch (dim) - { - case 2: - { - TestCases<2> tests; - MultipleParameterLoop input_data; - - tests.declare_parameters(input_data); - input_data.read_input (argv[2]); - input_data.loop (tests); - - break; - }; - - case 3: - { - TestCases<3> tests; - MultipleParameterLoop input_data; - - tests.declare_parameters(input_data); - input_data.read_input (argv[2]); - input_data.loop (tests); - - break; - }; - }; - - return 0; -}; - diff --git a/tests/big-tests/dof/dof_test.prm b/tests/big-tests/dof/dof_test.prm deleted file mode 100644 index 84697993f4..0000000000 --- a/tests/big-tests/dof/dof_test.prm +++ /dev/null @@ -1,4 +0,0 @@ -set Test run = { zoom in | ball | curved line | random } -set Grid file = {{ zoom_in | ball | curved_line | random }}.grid -set Sparsity file = {{ zoom_in | ball | curved_line | random }}.sparsity -set Condensed sparsity file = {{ zoom_in | ball | curved_line | random }}.sparsity.c \ No newline at end of file diff --git a/tests/big-tests/dof/make_ps b/tests/big-tests/dof/make_ps deleted file mode 100644 index 281a20bc8b..0000000000 --- a/tests/big-tests/dof/make_ps +++ /dev/null @@ -1,122 +0,0 @@ -set size 0.721,1 -set data style lines -set noxtics -set noytics -set noztics -set noxzeroaxis -set noyzeroaxis -set nokey -set term postscript eps - -!echo " Making <2d.zoom_in.grid.eps>" -set output "results/2d.zoom_in.grid.eps" -plot "results/2d.zoom_in.grid" - -!echo " Making <2d.ball.grid.eps>" -set output "results/2d.ball.grid.eps" -plot "results/2d.ball.grid" - -!echo " Making <2d.curved_line.grid.eps>" -set output "results/2d.curved_line.grid.eps" -plot "results/2d.curved_line.grid" - -!echo " Making <2d.random.grid.eps>" -set output "results/2d.random.grid.eps" -plot "results/2d.random.grid" - - - - -!echo " Making <3d.zoom_in.grid.eps>" -set output "results/3d.zoom_in.grid.eps" -splot "results/3d.zoom_in.grid" - -!echo " Making <3d.ball.grid.eps>" -set output "results/3d.ball.grid.eps" -splot "results/3d.ball.grid" - -!echo " Making <3d.curved_line.grid.eps>" -set output "results/3d.curved_line.grid.eps" -splot "results/3d.curved_line.grid" - -!echo " Making <3d.random.grid.eps>" -set output "results/3d.random.grid.eps" -splot "results/3d.random.grid" - - - - -set data style dots - -!echo " Making <2d.zoom_in.sparsity.eps>" -set output "results/2d.zoom_in.sparsity.eps" -plot "results/2d.zoom_in.sparsity" - -!echo " Making <2d.zoom_in.sparsity.c.eps>" -set output "results/2d.zoom_in.sparsity.c.eps" -plot "results/2d.zoom_in.sparsity.c" - - -!echo " Making <2d.ball.sparsity.eps>" -set output "results/2d.ball.sparsity.eps" -plot "results/2d.ball.sparsity" - -!echo " Making <2d.ball.sparsity.c.eps>" -set output "results/2d.ball.sparsity.c.eps" -plot "results/2d.ball.sparsity.c" - - -!echo " Making <2d.curved_line.sparsity.eps>" -set output "results/2d.curved_line.sparsity.eps" -plot "results/2d.curved_line.sparsity" - -!echo " Making <2d.curved_line.sparsity.c.eps>" -set output "results/2d.curved_line.sparsity.c.eps" -plot "results/2d.curved_line.sparsity.c" - - -!echo " Making <2d.random.sparsity.eps>" -set output "results/2d.random.sparsity.eps" -plot "results/2d.random.sparsity" - -!echo " Making <2d.random.sparsity.c.eps>" -set output "results/2d.random.sparsity.c.eps" -plot "results/2d.random.sparsity.c" - - - -!echo " Making <3d.zoom_in.sparsity.eps>" -set output "results/3d.zoom_in.sparsity.eps" -plot "results/3d.zoom_in.sparsity" - -!echo " Making <3d.zoom_in.sparsity.c.eps>" -set output "results/3d.zoom_in.sparsity.c.eps" -plot "results/3d.zoom_in.sparsity.c" - - -!echo " Making <3d.ball.sparsity.eps>" -set output "results/3d.ball.sparsity.eps" -plot "results/3d.ball.sparsity" - -!echo " Making <3d.ball.sparsity.c.eps>" -set output "results/3d.ball.sparsity.c.eps" -plot "results/3d.ball.sparsity.c" - - -!echo " Making <3d.curved_line.sparsity.eps>" -set output "results/3d.curved_line.sparsity.eps" -plot "results/3d.curved_line.sparsity" - -!echo " Making <3d.curved_line.sparsity.c.eps>" -set output "results/3d.curved_line.sparsity.c.eps" -plot "results/3d.curved_line.sparsity.c" - - -!echo " Making <3d.random.sparsity.eps>" -set output "results/3d.random.sparsity.eps" -plot "results/3d.random.sparsity" - -!echo " Making <3d.random.sparsity.c.eps>" -set output "results/3d.random.sparsity.c.eps" -plot "results/3d.random.sparsity.c" - diff --git a/tests/big-tests/dof/results/.cvsignore b/tests/big-tests/dof/results/.cvsignore deleted file mode 100644 index 480cb8565d..0000000000 --- a/tests/big-tests/dof/results/.cvsignore +++ /dev/null @@ -1 +0,0 @@ -?d.* diff --git a/tests/big-tests/error-estimation/.cvsignore b/tests/big-tests/error-estimation/.cvsignore deleted file mode 100644 index e3db7764a9..0000000000 --- a/tests/big-tests/error-estimation/.cvsignore +++ /dev/null @@ -1,4 +0,0 @@ -error-estimation -Makefile.dep -*.go -*.o diff --git a/tests/big-tests/error-estimation/Makefile b/tests/big-tests/error-estimation/Makefile deleted file mode 100644 index a3cf5fb381..0000000000 --- a/tests/big-tests/error-estimation/Makefile +++ /dev/null @@ -1,168 +0,0 @@ -# $Id$ -# Copyright by the deal.II authors 1998, 1999, 2000 - -############################################################ -# $D contains the root of the deal distribution and should -# be supplied as a command line argument D=../../... -############################################################ - - -# 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 = error-estimation - -# 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 - - -# 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 */*inp */*gnuplot */*eps - - - - -# -# -# 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).go $(libs.o) - flags = $(CXXFLAGS.o) -endif - - -# If in multithread mode, add the ACE library to the libraries which -# we need to link with: -ifneq ($(with-multithreading),no) - libraries += $(lib-ACE) -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) $(LDFLAGS) -o $@ $^ $(LIBS) - - -# 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) ee.gauss.prm - ./$(target) ee.singular.prm - ./$(target) ee.kink.prm - gnuplot make_ps - - - -# 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/tests/big-tests/error-estimation/ee.gauss.prm b/tests/big-tests/error-estimation/ee.gauss.prm deleted file mode 100644 index 560bd59cf3..0000000000 --- a/tests/big-tests/error-estimation/ee.gauss.prm +++ /dev/null @@ -1,8 +0,0 @@ -set Test case = Gauss shape -set Initial refinement = 2 -set Refinement criterion = { global | true error | estimated error } -set Refinement fraction = 0.3 -set Coarsening fraction = 0.03 -set Maximum cells = 10000 -set Output base filename = data-gauss/ -set Output format = ucd diff --git a/tests/big-tests/error-estimation/ee.kink.prm b/tests/big-tests/error-estimation/ee.kink.prm deleted file mode 100644 index 6bf0b8e88c..0000000000 --- a/tests/big-tests/error-estimation/ee.kink.prm +++ /dev/null @@ -1,8 +0,0 @@ -set Test case = Kink -set Initial refinement = 1 -set Refinement criterion = { global | estimated error } -set Refinement fraction = 0.1 -set Coarsening fraction = 0.02 -set Maximum cells = 100000 -set Output base filename = data-kink/ -set Output format = ucd diff --git a/tests/big-tests/error-estimation/ee.singular.prm b/tests/big-tests/error-estimation/ee.singular.prm deleted file mode 100644 index eb6988553a..0000000000 --- a/tests/big-tests/error-estimation/ee.singular.prm +++ /dev/null @@ -1,8 +0,0 @@ -set Test case = Singular -set Initial refinement = 1 -set Refinement criterion = { global | estimated error } -set Refinement fraction = 0.1 -set Coarsening fraction = 0.02 -set Maximum cells = 100000 -set Output base filename = data-singular/ -set Output format = ucd diff --git a/tests/big-tests/error-estimation/error-estimation.cc b/tests/big-tests/error-estimation/error-estimation.cc deleted file mode 100644 index 0a95248026..0000000000 --- a/tests/big-tests/error-estimation/error-estimation.cc +++ /dev/null @@ -1,759 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "../problem_base.h" -#include -#include -#include -#include - -#include -#include -#include -#include -#include - - - - -template -class PoissonEquation : public Equation { - public: - PoissonEquation (const Function &rhs) : - Equation(1), - use_coefficient(false), - right_hand_side (rhs), - coefficient (default_coefficient) {}; - - PoissonEquation (const Function &rhs, - const Function &coefficient ) : - Equation(1), - use_coefficient(true), - right_hand_side (rhs), - coefficient (coefficient) {}; - - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - protected: - const bool use_coefficient; - const Function &right_hand_side; - const Function &coefficient; - - static const ConstantFunction default_coefficient; -}; - - -const ConstantFunction<2> PoissonEquation<2>::default_coefficient(1); - - - - - -template -class PoissonProblem : public ProblemBase, public MultipleParameterLoop::UserClass { - public: - enum RefineMode { - global, true_error, error_estimator - }; - - PoissonProblem (); - ~PoissonProblem (); - - void clear (); - void create_new (const unsigned int); - void declare_parameters (ParameterHandler &prm); - void run (ParameterHandler &prm); - void print_history (const ParameterHandler &prm, - const RefineMode refine_mode) const; - - protected: - Triangulation *tria; - DoFHandler *dof; - - Function *rhs; - Function *solution_function; - Function *coefficient; - - Boundary *boundary; - - vector l2_error, linfty_error; - vector h1_error, estimated_error; - vector n_dofs; -}; - - - - - -template -class Solution { - public: - - class GaussShape : public Function { - public: - virtual double value (const Point &p, - const unsigned int component) const; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component) const; - }; - - class Singular : public Function { - public: - virtual double value (const Point &p, - const unsigned int component) const; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component) const; - }; - - class Kink : public Function { - public: - class Coefficient : public Function { - public: - virtual double value (const Point &p, - const unsigned int component) const; - }; - - virtual double value (const Point &p, - const unsigned int component) const; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component) const; - }; -}; - - - - -template -class RHS { - public: - - /** - * Right hand side constructed such that - * the exact solution is - * $x*y*exp(-(x**2+y**2)*40)$. - */ - class GaussShape : public Function { - public: - virtual double value (const Point &p, - const unsigned int component) const; - }; - - /** - * Right hand side constructed such that - * the exact solution is - * $r^{2/3}$. - */ - class Singular : public Function { - public: - virtual double value (const Point &p, - const unsigned int component) const; - }; - - /** - * Right hand side constructed such that - * the exact solution is - * $(1+4\theta(f))*f$ with - * $f=y-x**2$. - */ - class Kink : public Function { - public: - virtual double value (const Point &p, - const unsigned int component) const; - }; -}; - - - - -template <> -double Solution<2>::GaussShape::value (const Point<2> &p, - const unsigned int) const { - return p(0)*p(1)*exp(-40*p.square()); -}; - - -template <> -Tensor<1,2> Solution<2>::GaussShape::gradient (const Point<2> &p, - const unsigned int) const { - return Point<2> ((1-80.*p(0)*p(0))*p(1)*exp(-40*p.square()), - (1-80.*p(1)*p(1))*p(0)*exp(-40*p.square())); -}; - - - -template <> -double Solution<2>::Singular::value (const Point<2> &p, - const unsigned int) const { - return pow(p.square(), 1./3.); -}; - - -template <> -Tensor<1,2> Solution<2>::Singular::gradient (const Point<2> &p, - const unsigned int) const { - return 2./3.*pow(p.square(), -2./3.) * p; -}; - - - - -inline double theta(const double x) { - return (x>0 ? 1 : 0); -}; - - - -template <> -double Solution<2>::Kink::value (const Point<2> &p, - const unsigned int) const { - const double s = p(1)-p(0)*p(0); - return (1+4*theta(s))*s; -}; - - -template <> -Tensor<1,2> Solution<2>::Kink::gradient (const Point<2> &p, - const unsigned int) const { - const double s = p(1)-p(0)*p(0); - return (1+4*theta(s))*Point<2>(-2*p(0),1); -}; - - -template <> -double Solution<2>::Kink::Coefficient::value (const Point<2> &p, - const unsigned int) const { - const double s = p(1)-p(0)*p(0); - return 1./(1.+4.*theta(s)); -}; - - - -template <> -double RHS<2>::GaussShape::value (const Point<2> &p, - const unsigned int) const { - return (480.-6400.*p.square())*p(0)*p(1)*exp(-40.*p.square()); -}; - - -template <> -double RHS<2>::Singular::value (const Point<2> &p, - const unsigned int) const { - return -4./9. * pow(p.square(), -2./3.); -}; - - -template <> -double RHS<2>::Kink::value (const Point<2> &, - const unsigned int) const { - return 2; -}; - - - - - - - - -template <> -void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues<2> &fe_values, - const DoFHandler<2>::cell_iterator &) const { - for (unsigned int point=0; point -void PoissonEquation::assemble (FullMatrix &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, typename Equation::ExcPureVirtualFunctionCalled()); -}; - - - -template -void PoissonEquation::assemble (Vector &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, typename Equation::ExcPureVirtualFunctionCalled()); -}; - - - - - - - - - -template -PoissonProblem::PoissonProblem () : - tria(0), dof(0), rhs(0), - solution_function(0), coefficient(0), - boundary(0) {}; - - - -template -PoissonProblem::~PoissonProblem () -{ - clear (); -}; - - - -template -void PoissonProblem::clear () { - if (dof != 0) { delete dof; dof = 0; }; - if (tria != 0) { delete tria; tria = 0; }; - if (rhs != 0) { delete rhs; rhs = 0; }; - if (solution_function != 0) { delete solution_function; solution_function = 0; }; - if (coefficient != 0) { delete coefficient; coefficient = 0; }; - if (boundary != 0) { delete boundary; boundary = 0; }; - - // make it known to the underlying - // ProblemBase that tria and dof - // are already deleted - set_tria_and_dof (tria, dof); - - l2_error.clear (); - linfty_error.clear (); - h1_error.clear (); - estimated_error.clear(); - n_dofs.clear (); - - ProblemBase::clear (); -}; - - - - -template -void PoissonProblem::create_new (const unsigned int) { - clear (); - - tria = new Triangulation(); - dof = new DoFHandler (*tria); - set_tria_and_dof (tria, dof); - boundary = new HyperBallBoundary (); -}; - - - -template -void PoissonProblem::declare_parameters (ParameterHandler &prm) { - prm.declare_entry ("Test case", "Gauss shape", - Patterns::Selection("Gauss shape|Singular|Kink")); - prm.declare_entry ("Initial refinement", "2", - Patterns::Integer()); - prm.declare_entry ("Refinement criterion", "estimated error", - Patterns::Selection("global|true error|estimated error")); - prm.declare_entry ("Refinement fraction", "0.3", - Patterns::Double()); - prm.declare_entry ("Coarsening fraction", "0.1", - Patterns::Double()); - prm.declare_entry ("Maximum cells", "3000", - Patterns::Integer()); - prm.declare_entry ("Output base filename", ""); - prm.declare_entry ("Output format", "ucd", - Patterns::Selection("ucd|gnuplot")); -}; - - - - -template -void PoissonProblem::run (ParameterHandler &prm) { - cout << "=======================================" - << "=======================================" << endl - << "===== Test case: " << prm.get ("Test case") << endl - << "===== Doing computation with refinement criterion: "; - RefineMode refine_mode; - if (prm.get("Refinement criterion")=="global") - refine_mode = global; - else - if (prm.get("Refinement criterion")=="true error") - refine_mode = true_error; - else - if (prm.get("Refinement criterion")=="estimated error") - refine_mode = error_estimator; - else - return; - - switch (refine_mode) - { - case global: - cout << "global"; - break; - case true_error: - cout << "true error"; - break; - case error_estimator: - cout << "error estimator"; - break; - }; - - cout << endl - << "=======================================" - << "=======================================" << endl; - cout << "Making initial grid... " << endl; - const unsigned int start_level(prm.get_integer("Initial refinement")); - tria->set_boundary (0, *boundary); - GridGenerator::hyper_ball (*tria); - tria->refine_global (start_level); - - if (prm.get("Test case")=="Gauss shape") - rhs = new RHS::GaussShape(); - else - if (prm.get("Test case")=="Singular") - rhs = new RHS::Singular(); - else - if (prm.get("Test case")=="Kink") - rhs = new RHS::Kink(); - - if (prm.get("Test case")=="Gauss shape") - solution_function = new Solution::GaussShape (); - else - if (prm.get("Test case")=="Singular") - solution_function = new Solution::Singular (); - else - if (prm.get("Test case")=="Kink") - solution_function = new Solution::Kink (); - - - FEQ1 fe; - QGauss3 quadrature; - PoissonEquation *equation; - - static Solution::Kink::Coefficient kink_coefficient; - if (prm.get("Test case")=="Kink") - equation = new PoissonEquation(*rhs, kink_coefficient); - else - equation = new PoissonEquation(*rhs); - - SolutionTransfer solution_transfer (*dof_handler); - - unsigned int refine_step = 0; - const unsigned int max_cells = prm.get_integer("Maximum cells"); - while (tria->n_active_cells() < max_cells) - { - Vector old_solution = solution; - cout << "Refinement step " << refine_step - << ", using " << tria->n_active_cells() << " active cells on " - << tria->n_levels() << " levels." - << endl; - cout << " Distributing dofs... "; - dof->distribute_dofs (fe); - cout << dof->n_dofs() << " degrees of freedom." << endl; - n_dofs.push_back (dof->n_dofs()); - - cout << " Assembling matrices..." << endl; - UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | - update_gradients | update_JxW_values); - - ProblemBase::FunctionMap dirichlet_bc; - dirichlet_bc[0] = solution_function; - assemble (*equation, quadrature, update_flags, dirichlet_bc); - - // if we have an old solution lying - // around, use it to preset the solution - // vector. this reduced the quired - // number of iterations by about - // 10 per cent - if (refine_step != 0) - { - solution.reinit (dof_handler->n_dofs()); - solution_transfer.interpolate (old_solution, solution); - - // if you don't want to preset - // the solution vector, - // uncomment the following - // line and comment out the - // preceding one -// solution.reinit (dof_handler->n_dofs()); - - solution_transfer.clear (); - }; - - cout << " Solving..." << endl; - - solve (); - - - Vector l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell; - Vector estimated_error_per_cell; - QGauss3 q; - - cout << " Calculating L2 error... "; - VectorTools::integrate_difference (*dof_handler, - solution, *solution_function, - l2_error_per_cell, q, - L2_norm); - cout << l2_error_per_cell.l2_norm() << endl; - l2_error.push_back (l2_error_per_cell.l2_norm()); - - cout << " Calculating L-infinity error... "; - VectorTools::integrate_difference (*dof_handler, - solution, *solution_function, - linfty_error_per_cell, q, - Linfty_norm); - cout << linfty_error_per_cell.linfty_norm() << endl; - linfty_error.push_back (linfty_error_per_cell.linfty_norm()); - - cout << " Calculating H1 error... "; - VectorTools::integrate_difference (*dof_handler, - solution, *solution_function, - h1_error_per_cell, q, - H1_norm); - cout << h1_error_per_cell.l2_norm() << endl; - h1_error.push_back (h1_error_per_cell.l2_norm()); - - cout << " Estimating H1 error... "; - - QSimpson eq; - KellyErrorEstimator::estimate (*dof, eq, - FunctionMap::type(), - solution, - estimated_error_per_cell, - vector(), // all components - ((prm.get("Test case")=="Kink") ? - &kink_coefficient : 0 )); - cout << estimated_error_per_cell.l2_norm() << endl; - estimated_error.push_back (estimated_error_per_cell.l2_norm()); - - Vector l2_error_per_dof(dof->n_dofs()), linfty_error_per_dof(dof->n_dofs()); - Vector h1_error_per_dof(dof->n_dofs()), estimated_error_per_dof(dof->n_dofs()); - Vector error_ratio (dof->n_dofs()); - DoFTools::distribute_cell_to_dof_vector (*dof, l2_error_per_cell, l2_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, linfty_error_per_cell, - linfty_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, h1_error_per_cell, h1_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, estimated_error_per_cell, - estimated_error_per_dof); - error_ratio.ratio (h1_error_per_dof, estimated_error_per_dof); - - DataOut out; - fill_data (out); - out.add_data_vector (l2_error_per_dof, "L2_Error"); - out.add_data_vector (linfty_error_per_dof, "Linfty_Error"); - out.add_data_vector (h1_error_per_dof, "H1_Error"); - out.add_data_vector (estimated_error_per_dof, "Estimated_Error"); - out.add_data_vector (error_ratio, "Ratio_True_to_Estimated_Error"); - out.build_patches (); - string filename = prm.get ("Output base filename"); - switch (refine_mode) - { - case global: - filename += "global."; - break; - case true_error: - filename += "true_error."; - break; - case error_estimator: - filename += "estimated_error."; - break; - }; - filename += ('0'+(start_level+refine_step)/10); - filename += ('0'+(start_level+refine_step)%10); - - if (prm.get("Output format")=="ucd") - filename += ".inp"; - else - if (prm.get("Output format")=="gnuplot") - filename += ".gnuplot"; - - cout << " Writing error plots to <" << filename << ">..." << endl; - ofstream outfile(filename.c_str()); - if (prm.get("Output format")=="ucd") - out.write_ucd (outfile); - else - if (prm.get("Output format")=="gnuplot") - out.write_gnuplot (outfile); - - outfile.close(); - - cout << " Refining triangulation..."; - switch (refine_mode) - { - case global: - tria->set_all_refine_flags (); - break; - case true_error: - GridRefinement:: - refine_and_coarsen_fixed_number (*tria, - h1_error_per_cell, - prm.get_double("Refinement fraction"), - prm.get_double("Coarsening fraction")); - break; - case error_estimator: - GridRefinement:: - refine_and_coarsen_fixed_number (*tria, - estimated_error_per_cell, - prm.get_double("Refinement fraction"), - prm.get_double("Coarsening fraction")); - break; - }; - - tria->prepare_coarsening_and_refinement (); - solution_transfer.prepare_for_coarsening_and_refinement (solution); - tria->execute_coarsening_and_refinement (); - - cout << endl << endl; - ++refine_step; - }; - - string filename = prm.get ("Output base filename"); - switch (refine_mode) - { - case global: - filename += "global."; - break; - case true_error: - filename += "true_error."; - break; - case error_estimator: - filename += "estimated_error."; - break; - }; - - cout << endl; - - filename += "finest_mesh.gnuplot"; - cout << " Writing finest grid to <" << filename << ">... " << endl; - ofstream finest_mesh (filename.c_str()); - GridOut().write_gnuplot (*tria, finest_mesh); - finest_mesh.close(); - - print_history (prm, refine_mode); - cout << endl << endl << endl; - - dof->clear (); - delete equation; -}; - - -template -void PoissonProblem::print_history (const ParameterHandler &prm, - const RefineMode refine_mode) const { - string filename(prm.get("Output base filename")); - filename += "history."; - switch (refine_mode) - { - case global: - filename += "global."; - break; - case true_error: - filename += "true_error."; - break; - case error_estimator: - filename += "estimated_error."; - break; - }; - filename += "gnuplot"; - - cout << endl << "Printing convergence history to <" << filename << ">..." - << endl; - ofstream out(filename.c_str()); - out << "# n_dofs l2_error linfty_error " - << "h1_error estimated_error" - << endl; - for (unsigned int i=0; ih/2:" << endl; - cout << " L2 error : " << 1./average_l2 << endl - << " Linfty error : " << 1./average_linfty << endl - << " H1 error : " << 1./average_h1 << endl - << " Estimated error : " << 1./average_est << endl; -}; - - - - -int main (int argc, char **argv) { - if (argc!=2) - { - cout << "Usage: error-estimation parameterfile" << endl << endl; - return 1; - }; - - PoissonProblem<2> poisson; - MultipleParameterLoop input_data; - - poisson.declare_parameters(input_data); - input_data.read_input (argv[1]); - input_data.loop (poisson); - - return 0; -}; - - - diff --git a/tests/big-tests/error-estimation/make_ps b/tests/big-tests/error-estimation/make_ps deleted file mode 100644 index 5c06b6c8f2..0000000000 --- a/tests/big-tests/error-estimation/make_ps +++ /dev/null @@ -1,94 +0,0 @@ -set xlabel "Number of degrees of freedom" -set ylabel "Error" -set data style linespoints -set logscale xy - -set term postscript eps - - -set output "data-gauss/history.global.eps" - -plot "data-gauss/history.global.gnuplot" using 1:2 title "L2 error","data-gauss/history.global.gnuplot" using 1:3 title "Linfty error","data-gauss/history.global.gnuplot" using 1:4 title "H1 error","data-gauss/history.global.gnuplot" using 1:5 title "Estimated H1 error" - - -set output "data-gauss/history.true_error.eps" - -plot "data-gauss/history.true_error.gnuplot" using 1:2 title "L2 error","data-gauss/history.true_error.gnuplot" using 1:3 title "Linfty error","data-gauss/history.true_error.gnuplot" using 1:4 title "H1 error","data-gauss/history.true_error.gnuplot" using 1:5 title "Estimated H1 error" - - -set output "data-gauss/history.estimated_error.eps" - -plot "data-gauss/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-gauss/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-gauss/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-gauss/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" - - -set output "data-gauss/history.compare.eps" -plot "data-gauss/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-gauss/history.true_error.gnuplot" using 1:2 title "ref. by true error -- L2 error", "data-gauss/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", 0.1/sqrt(x) title "O(h)", "data-gauss/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-gauss/history.true_error.gnuplot" using 1:4 title "ref. by true error -- H1 error", "data-gauss/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error", 0.04/x title "O(h^2)" - - - - - -set output "data-singular/history.global.eps" - -plot "data-singular/history.global.gnuplot" using 1:2 title "L2 error","data-singular/history.global.gnuplot" using 1:3 title "Linfty error","data-singular/history.global.gnuplot" using 1:4 title "H1 error","data-singular/history.global.gnuplot" using 1:5 title "Estimated H1 error" - - -set output "data-singular/history.estimated_error.eps" - -plot "data-singular/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-singular/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-singular/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-singular/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" - - -set output "data-singular/history.compare.eps" -plot "data-singular/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-singular/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", 1.1/x**0.33 title "O(h^2/3)", 2./sqrt(x) title "O(h)", "data-singular/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-singular/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error", 0.2/x**0.4 title "O(h^0.8)", 4./x title "O(h^2)" - - - - - - - -set output "data-kink/history.global.eps" - -plot "data-kink/history.global.gnuplot" using 1:2 title "L2 error","data-kink/history.global.gnuplot" using 1:3 title "Linfty error","data-kink/history.global.gnuplot" using 1:4 title "H1 error","data-kink/history.global.gnuplot" using 1:5 title "Estimated H1 error" - - -set output "data-kink/history.estimated_error.eps" - -plot "data-kink/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-kink/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-kink/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-kink/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" - - -set output "data-kink/history.compare.eps" -plot "data-kink/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-kink/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", 5/x**0.25 title "O(h^1/2)", 20/x**0.5 title "O(h)", "data-kink/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-kink/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error", 1.5/sqrt(x) title "O(h)", 20/x**0.95 title "O(h^1.8)" - - - - -set parametric -set data style lines -set nologscale xy -set size 0.7,1 - -set output "data-gauss/finest_mesh.global.eps" -plot "data-gauss/global.finest_mesh.gnuplot" title "Finest mesh" - -set output "data-gauss/finest_mesh.true_error.eps" -plot "data-gauss/true_error.finest_mesh.gnuplot" title "Finest mesh" - -set output "data-gauss/finest_mesh.estimated_error.eps" -plot "data-gauss/estimated_error.finest_mesh.gnuplot" title "Finest mesh" - - - -set output "data-singular/finest_mesh.global.eps" -plot "data-singular/global.finest_mesh.gnuplot" title "Finest mesh" - -set output "data-singular/finest_mesh.estimated_error.eps" -plot "data-singular/estimated_error.finest_mesh.gnuplot" title "Finest mesh" - - - -set output "data-kink/finest_mesh.global.eps" -plot "data-kink/global.finest_mesh.gnuplot" title "Finest mesh" - -set output "data-kink/finest_mesh.estimated_error.eps" -plot "data-kink/estimated_error.finest_mesh.gnuplot" title "Finest mesh" diff --git a/tests/big-tests/error-estimation/strip_comments b/tests/big-tests/error-estimation/strip_comments deleted file mode 100755 index 779b6b16c7..0000000000 --- a/tests/big-tests/error-estimation/strip_comments +++ /dev/null @@ -1 +0,0 @@ -perl -pi -e 's/^#.*$\\n//g' data-*/*.inp diff --git a/tests/big-tests/grid/.cvsignore b/tests/big-tests/grid/.cvsignore deleted file mode 100644 index d91582a08e..0000000000 --- a/tests/big-tests/grid/.cvsignore +++ /dev/null @@ -1,4 +0,0 @@ -grid_test -Makefile.dep -*.go -*.o diff --git a/tests/big-tests/grid/Makefile b/tests/big-tests/grid/Makefile deleted file mode 100644 index eb60e99baa..0000000000 --- a/tests/big-tests/grid/Makefile +++ /dev/null @@ -1,175 +0,0 @@ -# $Id$ -# Copyright by the deal.II authors 1998, 1999, 2000 - -############################################################ -# $D contains the root of the deal distribution and should -# be supplied as a command line argument D=../../... -############################################################ - - -# 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 = grid_test - -# 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 - - - -# 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-deal2-3d.g) \ - $(lib-lac.g) \ - $(lib-base.g) -libs.o = $(lib-deal2-2d.o) \ - $(lib-deal2-3d.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).go $(libs.o) - flags = $(CXXFLAGS.o) -endif - - -# If in multithread mode, add the ACE library to the libraries which -# we need to link with: -ifneq ($(with-multithreading),no) - libraries += $(lib-ACE) -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) $(LDFLAGS) -o $@ $^ $(LIBS) - - -# 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) 2 1 - ./$(target) 2 2 - ./$(target) 2 3 - ./$(target) 2 4 - ./$(target) 3 1 - ./$(target) 3 2 - ./$(target) 3 3 - ./$(target) 3 4 - - - -# 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/tests/big-tests/grid/grid_test.cc b/tests/big-tests/grid/grid_test.cc deleted file mode 100644 index c5d4c0800a..0000000000 --- a/tests/big-tests/grid/grid_test.cc +++ /dev/null @@ -1,329 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - - -// 1: continuous refinement of the unit square always in the middle -// 2: refinement of the circle at the boundary -// 2: refinement of a wiggled area at the boundary -// 4: random refinement - - - - - -template -class Ball : - public StraightBoundary { - public: - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const { - Point middle = StraightBoundary::get_new_point_on_line(line); - - for (int i=0; i - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { - Point middle = StraightBoundary::get_new_point_on_quad(quad); - - for (int i=0; i -class CurvedLine : - public StraightBoundary { - public: - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - virtual Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; -}; - - - -template -Point -CurvedLine::get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - Point middle = StraightBoundary::get_new_point_on_line (line); - - // if the line is at the top of bottom - // face: do a special treatment on - // this line. Note that if the - // z-value of the midpoint is either - // 0 or 1, then the z-values of all - // vertices of the line is like that - if (dim>=3) - if (((middle(2) == 0) || (middle(2) == 1)) - // find out, if the line is in the - // interior of the top or bottom face - // of the domain, or at the edge. - // lines at the edge need to undergo - // the usual treatment, while for - // interior lines taking the midpoint - // is sufficient - // - // note: the trick with the boundary - // id was invented after the above was - // written, so we are not very strict - // here with using these flags - && (line->boundary_indicator() == 1)) - return middle; - - - double x=middle(0), - y=middle(1); - - if (y -Point -CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const -{ - Point middle = StraightBoundary::get_new_point_on_quad (quad); - - // if the face is at the top of bottom - // face: do not move the midpoint in - // x/y direction. Note that if the - // z-value of the midpoint is either - // 0 or 1, then the z-values of all - // vertices of the quad is like that - if ((middle(2) == 0) || (middle(2) == 1)) - return middle; - - double x=middle(0), - y=middle(1); - - if (y -void test (const int test_case) { - cout << "Running testcase " << test_case - << " in " << dim << " dimensions." << endl; - Triangulation tria; - GridGenerator::hyper_cube(tria); - - if ((dim==1) && ((test_case==2) || (test_case==3))) - { - cout << "Impossible for this dimension." << endl; - return; - }; - - - switch (test_case) - { - case 1: - { - // we want to log the - // refinement history -// ofstream history ("mesh.history"); - - // refine first cell - tria.begin_active()->set_refine_flag(); -// tria.save_refine_flags (history); - tria.execute_coarsening_and_refinement (); - - // refine first active cell - // on coarsest level - tria.begin_active()->set_refine_flag (); -// tria.save_refine_flags (history); - tria.execute_coarsening_and_refinement (); - - Triangulation::active_cell_iterator cell; - for (int i=0; i<17; ++i) - { - // refine the presently - // second last cell 17 - // times - cell = tria.last_active(tria.n_levels()-1); - --cell; - cell->set_refine_flag (); -// tria.save_refine_flags (history); - tria.execute_coarsening_and_refinement (); - }; - -// tria.refine_global (5); - - break; - } - - case 2: - case 3: - { - if (dim==3) - { - tria.begin_active()->face(2)->set_boundary_indicator(1); - tria.begin_active()->face(4)->set_boundary_indicator(1); - }; - - - // set the boundary function - Ball ball; - CurvedLine curved_line; - if (test_case==2) - { - tria.set_boundary (0, ball); - tria.set_boundary (1, ball); - } else { - tria.set_boundary (0, curved_line); - tria.set_boundary (1, curved_line); - }; - - // refine once - tria.begin_active()->set_refine_flag(); - tria.execute_coarsening_and_refinement (); - - Triangulation::active_cell_iterator cell, endc; - const unsigned int steps[4] = { 0, 10, 7, 2 }; - for (unsigned int i=0; iat_boundary()) - cell->set_refine_flag(); - - tria.execute_coarsening_and_refinement(); - }; - - tria.set_boundary (0); - tria.set_boundary (1); - - break; - } - - case 4: - { - // refine once - tria.begin_active()->set_refine_flag(); - tria.execute_coarsening_and_refinement (); - - Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<(dim==2 ? 13 : (dim==3 ? 7 : 30)); ++i) - { - int n_levels = tria.n_levels(); - cell = tria.begin_active(); - endc = tria.end(); - - for (; cell!=endc; ++cell) - { - double r = rand()*1.0/RAND_MAX, - weight = 1.* - (cell->level()*cell->level()) / - (n_levels*n_levels); - - if (r <= 0.5*weight) - cell->set_refine_flag (); - }; - - tria.execute_coarsening_and_refinement (); - }; - break; - } - }; - - - - // output the grid - string filename("results/"); - filename += ('0'+dim); - filename += "d."; - filename += ('0'+test_case); - filename += ".eps"; - - ofstream out(filename.c_str()); - GridOut grid_out; - GridOut::EpsFlags<3> eps_flags; - eps_flags.azimut_angle += 20; - eps_flags.turn_angle += 20; - grid_out.set_flags (eps_flags); - grid_out.write_eps (tria, out); - - cout << " Total number of cells = " << tria.n_cells() << endl - << " Total number of active cells = " << tria.n_active_cells() << endl; -}; - - - -int main (int argc, char **argv) { - if (argc!=3) - { - cout << "Usage: grid_test dimension testcase" << endl << endl - << "Dimension: 2 or 3" << endl << endl - << "Testcases:" << endl - << " 1: continuous refinement of the unit square/cube always in the middle" << endl - << " 2: refinement of the circle/sphere at the boundary" << endl - << " 3: refinement of a wiggled area at the boundary" << endl - << " 4: random refinement" << endl << endl; - return 1; - }; - - if (argv[1][0] == '2') - test<2> (argv[2][0]-'0'); - else - test<3> (argv[2][0]-'0'); - - return 0; -}; diff --git a/tests/big-tests/grid/make_ps b/tests/big-tests/grid/make_ps deleted file mode 100644 index 21782a3337..0000000000 --- a/tests/big-tests/grid/make_ps +++ /dev/null @@ -1,43 +0,0 @@ -set size 0.721,1 -set data style lines -set noxtics -set noytics -set noztics -set noxzeroaxis -set noyzeroaxis -#set nozzeroaxis -set nokey -set term postscript eps - -!echo " Making " -set output "results/2d.1.eps" -plot "results/2d.1" - -!echo " Making " -set output "results/2d.2.eps" -plot "results/2d.2" - -!echo " Making " -set output "results/2d.3.eps" -plot "results/2d.3" - -!echo " Making " -set output "results/2d.4.eps" -plot "results/2d.4" - - -!echo " Making " -set output "results/3d.1.eps" -splot "results/3d.1" - -!echo " Making " -set output "results/3d.2.eps" -splot "results/3d.2" - -!echo " Making " -set output "results/3d.3.eps" -splot "results/3d.3" - -!echo " Making " -set output "results/3d.4.eps" -splot "results/3d.4" diff --git a/tests/big-tests/grid/results/.cvsignore b/tests/big-tests/grid/results/.cvsignore deleted file mode 100644 index 480cb8565d..0000000000 --- a/tests/big-tests/grid/results/.cvsignore +++ /dev/null @@ -1 +0,0 @@ -?d.* diff --git a/tests/big-tests/multigrid/.cvsignore b/tests/big-tests/multigrid/.cvsignore deleted file mode 100644 index ab98be5ee1..0000000000 --- a/tests/big-tests/multigrid/.cvsignore +++ /dev/null @@ -1,4 +0,0 @@ -multigrid -Makefile.dep -*.go -*.o diff --git a/tests/big-tests/multigrid/Makefile b/tests/big-tests/multigrid/Makefile deleted file mode 100644 index d7a6c9c035..0000000000 --- a/tests/big-tests/multigrid/Makefile +++ /dev/null @@ -1,165 +0,0 @@ -# $Id$ -# Copyright by the deal.II authors 1998, 1999, 2000 - -############################################################ -# $D contains the root of the deal distribution and should -# be supplied as a command line argument D=../../... -############################################################ - - -# 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 = multigrid - -# 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 - - -# 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).go $(libs.o) - flags = $(CXXFLAGS.o) -endif - - -# If in multithread mode, add the ACE library to the libraries which -# we need to link with: -ifneq ($(with-multithreading),no) - libraries += $(lib-ACE) -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) $(LDFLAGS) -o $@ $^ $(LIBS) - - -# 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) - gnuplot make_ps - - -# 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/tests/big-tests/multigrid/make_ps b/tests/big-tests/multigrid/make_ps deleted file mode 100644 index 76c13a8624..0000000000 --- a/tests/big-tests/multigrid/make_ps +++ /dev/null @@ -1,52 +0,0 @@ -set term postscript eps -set xlabel "Number of degrees of freedom" -set data style linespoints -set logscale xy - - - -set ylabel "Error" - -set output "criss-cross.eps" - -plot "criss_cross.history" using 1:2 title "L1 error","criss_cross.history" using 1:3 title "L2 error","criss_cross.history" using 1:4 title "Linfty error","criss_cross.history" using 1:5 title "H1 seminorm error","criss_cross.history" using 1:6 title "H1 error" - - - -set output "linear.eps" - -plot "linear.history" using 1:2 title "L1 error","linear.history" using 1:3 title "L2 error","linear.history" using 1:4 title "Linfty error","linear.history" using 1:5 title "H1 seminorm error","linear.history" using 1:6 title "H1 error" - - - -set output "quadratic.eps" - -plot "quadratic.history" using 1:2 title "L1 error","quadratic.history" using 1:3 title "L2 error","quadratic.history" using 1:4 title "Linfty error","quadratic.history" using 1:5 title "H1 seminorm error","quadratic.history" using 1:6 title "H1 error" - - - -set output "cubic.eps" - -plot "cubic.history" using 1:2 title "L1 error","cubic.history" using 1:3 title "L2 error","cubic.history" using 1:4 title "Linfty error","cubic.history" using 1:5 title "H1 seminorm error","cubic.history" using 1:6 title "H1 error" - - - -set output "quartic.eps" - -plot "quartic.history" using 1:2 title "L1 error","quartic.history" using 1:3 title "L2 error","quartic.history" using 1:4 title "Linfty error","quartic.history" using 1:5 title "H1 seminorm error","quartic.history" using 1:6 title "H1 error" - - - -set output "l2error.eps" -set ylabel "L2-error" - -plot "criss_cross.history" using 1:3 title "Criss-cross elements", "linear.history" using 1:3 title "Linear elements", "quadratic.history" using 1:3 title "Quadratic elements", "cubic.history" using 1:3 title "Cubic elements", "quartic.history" using 1:3 title "Quartic elements" - - - -set output "h1error.eps" -set ylabel "H1-error" - -plot "criss_cross.history" using 1:6 title "Criss-cross elements", "linear.history" using 1:6 title "Linear elements", "quadratic.history" using 1:6 title "Quadratic elements", "cubic.history" using 1:6 title "Cubic elements", "quartic.history" using 1:6 title "Quartic elements" - - diff --git a/tests/big-tests/multigrid/multigrid.cc b/tests/big-tests/multigrid/multigrid.cc deleted file mode 100644 index e28b5841d2..0000000000 --- a/tests/big-tests/multigrid/multigrid.cc +++ /dev/null @@ -1,515 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - - -#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 -#include -#include - -#include - -#include - - - -template -class LaplaceProblem -{ - public: - LaplaceProblem (); - ~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; - MGDoFHandler mg_dof_handler; - - FEQ1 fe; - - ConstraintMatrix hanging_node_constraints; - - SparsityPattern global_sparsity_pattern; - SparseMatrix global_system_matrix; - - MGLevelObject level_sparsity_patterns; - MGLevelObject > level_system_matrices; - - 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, - ExcDimensionMismatch (values.size(), n_points)); - - Assert (component == 0, - ExcIndexRange (component, 0, 1)); - - for (unsigned int i=0; i > >matrices; - public: - MGSmootherLAC(MGLevelObject >&); - - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; -}; - - -MGSmootherLAC::MGSmootherLAC(MGLevelObject >& matrix) - : - matrices(&matrix) -{} - - -void -MGSmootherLAC::smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const -{ - SolverControl control(2,1.e-300,false,false); - PrimitiveVectorMemory<> mem; - SolverRichardson<> rich(control, mem); - PreconditionSSOR<> prec; - prec.initialize((*matrices)[level], 1.); - - rich.solve((*matrices)[level], u, rhs, prec); -} - - - -template -LaplaceProblem::LaplaceProblem () : - mg_dof_handler (triangulation) -{}; - - - -template -LaplaceProblem::~LaplaceProblem () -{ - mg_dof_handler.clear (); -}; - - - -template -void LaplaceProblem::setup_system () -{ - mg_dof_handler.distribute_dofs (fe); - - hanging_node_constraints.clear (); - DoFTools::make_hanging_node_constraints (mg_dof_handler, - hanging_node_constraints); - hanging_node_constraints.close (); - global_sparsity_pattern.reinit (mg_dof_handler.DoFHandler::n_dofs(), - mg_dof_handler.DoFHandler::n_dofs(), - mg_dof_handler.max_couplings_between_dofs()); - DoFTools::make_sparsity_pattern (mg_dof_handler, global_sparsity_pattern); - hanging_node_constraints.condense (global_sparsity_pattern); - global_sparsity_pattern.compress(); - - global_system_matrix.reinit (global_sparsity_pattern); - - solution.reinit (mg_dof_handler.DoFHandler::n_dofs()); - system_rhs.reinit (mg_dof_handler.DoFHandler::n_dofs()); - - - const unsigned int n_levels = triangulation.n_levels(); - level_system_matrices.resize (0, n_levels); - level_sparsity_patterns.resize (0, n_levels); - - for (unsigned int level=0; level -void LaplaceProblem::assemble_system () -{ - const Coefficient coefficient; - - 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); - - vector local_dof_indices (dofs_per_cell); - - // FIX - vector coefficient_values (n_q_points, 1.0); - - // not only active cells - MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), - endc = mg_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(); - - // FIX -// coefficient.value_list (q_points, coefficient_values); - - for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); - const unsigned int level = cell->level(); - for (unsigned int i=0; iactive()) - { - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i boundary_values; -// VectorTools::interpolate_boundary_values (mg_dof_handler, -// 0, -// ZeroFunction(), -// boundary_values); -// MatrixTools::apply_boundary_values (boundary_values, -// global_system_matrix, -// solution, -// system_rhs); -}; - - - -template -void LaplaceProblem::solve () -{ - - { - SolverControl solver_control (1000, 1e-12); - PrimitiveVectorMemory<> vector_memory; - SolverCG<> cg (solver_control, vector_memory); - - SolverControl coarse_grid_solver_control (1000, 1e-12); - PrimitiveVectorMemory<> coarse_grid_vector_memory; - - SolverCG<> coarse_grid_cg (coarse_grid_solver_control, - coarse_grid_vector_memory); - -// PreconditionRelaxation<> -// coarse_grid_solver_preconditioner(level_system_matrices[level_system_matrices.get_minlevel()], -// &SparseMatrix::template precondition_SSOR, -// 1.2); - PreconditionIdentity coarse_grid_solver_preconditioner; - - MGCoarseGridLACIteration, SparseMatrix, PreconditionIdentity> - coarse_grid_solver (coarse_grid_cg, - level_system_matrices[level_system_matrices.get_minlevel()], - coarse_grid_solver_preconditioner); - - MGSmootherLAC smoother (level_system_matrices); - MGTransferPrebuilt grid_transfer; - grid_transfer.build_matrices (mg_dof_handler); - - Multigrid<2> multigrid (mg_dof_handler, - hanging_node_constraints, - level_sparsity_patterns, - level_system_matrices, - grid_transfer); - - PreconditionMG > - mg_precondition (multigrid, smoother, smoother, coarse_grid_solver); - - solution.clear (); - cg.solve (global_system_matrix, solution, system_rhs, - mg_precondition); - - cout << " MG Outer iterations: " << solver_control.last_step() - << endl; - - cout << " MG Total inner iterations: " << coarse_grid_solver_control.last_step() - << endl; - }; - - { - SolverControl solver_control (1000, 1e-12); - PrimitiveVectorMemory<> vector_memory; - SolverCG<> cg (solver_control, vector_memory); - - PreconditionSSOR<> preconditioner; - preconditioner.initialize(global_system_matrix, 1.2); - - solution.clear (); - cg.solve (global_system_matrix, solution, system_rhs, - preconditioner); - - cout << " CG Outer iterations: " << solver_control.last_step() - << endl; - }; - - hanging_node_constraints.distribute (solution); -}; - - -template -void LaplaceProblem::refine_grid () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - FunctionMap::type neumann_boundary; - KellyErrorEstimator::estimate (mg_dof_handler, - QGauss3(), - 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 (); -}; - - - -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_cube (triangulation); - triangulation.refine_global (1); - } - else - { - refine_grid (); - }; - - - cout << " Number of active cells: " - << triangulation.n_active_cells() - << endl; - - setup_system (); - - cout << " Number of degrees of freedom: " - << mg_dof_handler.DoFHandler::n_dofs() - << endl; - - assemble_system (); - solve (); - output_results (cycle); - - typename DataOut::EpsFlags eps_flags; - eps_flags.z_scaling = 4; - - DataOut data_out; - data_out.set_flags (eps_flags); - - data_out.attach_dof_handler (mg_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 () -{ - try - { - deallog.depth_console (0); - - LaplaceProblem<2> laplace_problem_2d; - 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/tests/big-tests/nonlinear/fixed-point-iteration/.cvsignore b/tests/big-tests/nonlinear/fixed-point-iteration/.cvsignore deleted file mode 100644 index 68fd34ed36..0000000000 --- a/tests/big-tests/nonlinear/fixed-point-iteration/.cvsignore +++ /dev/null @@ -1,3 +0,0 @@ -Makefile.dep -nonlinear -Obj.g.a diff --git a/tests/big-tests/nonlinear/fixed-point-iteration/Makefile b/tests/big-tests/nonlinear/fixed-point-iteration/Makefile deleted file mode 100644 index 5322c80ace..0000000000 --- a/tests/big-tests/nonlinear/fixed-point-iteration/Makefile +++ /dev/null @@ -1,138 +0,0 @@ -# $Id$ -# Copyright by the deal.II authors 1998, 1999, 2000 - -############################################################ -# $D contains the root of the deal distribution and should -# be supplied as a command line argument D=../../... -############################################################ - - -# Template for makefiles for the examples subdirectory. In principle, -# everything should be done automatically if you set the target file -# here correctly: -target = nonlinear - -# 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 - -# If you want your program to be linked with extra object or library -# files, specify them here: -user-libs = - -# To run the program, use "make run"; to give parameters to the program, -# give the parameters to the following variable: -run-parameters = $(target).prm - -# To execute additional action apart from running the program, fill -# in this list: -additional-run-action = - -# To specify which files are to be deleted by "make clean" (apart from -# the usual ones: object files, executables, backups, etc), fill in the -# following list -delete-files = *gnuplot *.eps - - - - -############################################################################### -# Internals - - -include $D/common/Make.global_options - - - -# get lists of files we need -cc-files = $(filter-out *%, $(shell echo *.cc)) -o-files = $(cc-files:.cc=.o) -go-files = $(cc-files:.cc=.go) -h-files = $(filter-out *%, $(shell echo *.h)) -lib-h-files = $(filter-out *%, $(shell echo $D/base/include/base/*.h \ - $D/lac/include/lac/*.h \ - $D/deal.II/include/*/*.h)) - -# list of libraries needed to link with -libs.g = ./Obj.g.a \ - $(lib-deal2-2d.g) \ - $(lib-lac.g) \ - $(lib-base.g) -libs = ./Obj.a \ - $(lib-deal2-2d.o) \ - $(lib-lac.o) \ - $(lib-base.o) - - -# check whether we use debug mode or not -ifeq ($(debug-mode),on) -libraries = $(libs.g) -flags = $(CXXFLAGS.g) -endif - -ifeq ($(debug-mode),off) -libraries = $(libs) -flags = $(CXXFLAGS) -endif - -# If in multithread mode, add the ACE library to the libraries which -# we need to link with: -ifneq ($(with-multithreading),no) - libraries += $(lib-ACE) -endif - - -# make rule for the target -$(target) : $(libraries) - @echo ============================ Linking $@ - @$(CXX) $(LDFLAGS) -o $@ $^ $(LIBS) $(user-libs) - -# rule how to run the program -run: $(target) - ./$(target) $(run-parameters) - - -# rule to make object files -%.go : %.cc - @echo ============================ Compiling with debugging information: $< - @echo $(CXX) ... -c $< -o $@ - @$(CXX) $(CXXFLAGS.g) -c $< -o $@ -%.o : %.cc - @echo ============================ Compiling with optimization: $< - @echo $(CXX) ... -c $< -o $@ - @$(CXX) $(CXXFLAGS.o) -c $< -o $@ - - -# rules which files the libraries depend upon -Obj.a: ./Obj.a($(o-files)) -Obj.g.a: ./Obj.g.a($(go-files)) - - -clean: - -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files) - - - -.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. -# -#use perl to generate rules for the .go files as well -#as to make rules not for tria.o and the like, but -#rather for libnumerics.a(tria.o) -Makefile.dep: $(cc-files) $(h-files) $(lib-h-files) - @echo ============================ Remaking Makefile - @$(PERL) $D/common/scripts/make_dependencies.pl ./Obj $(INCLUDE) $(cc-files) \ - > Makefile.dep - - -include Makefile.dep - diff --git a/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc b/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc deleted file mode 100644 index c599249d46..0000000000 --- a/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc +++ /dev/null @@ -1,258 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "../../problem_base.h" -#include -#include - -#include -#include -#include -#include - - - - - -template -class RightHandSide : public Function -{ - public: - virtual double value (const Point &p, - const unsigned int) const - { - double x = 80; - for (unsigned int d=0; d -class PoissonEquation : public Equation { - public: - PoissonEquation (const Function &rhs, - const Vector &last_solution) : - Equation(1), - right_hand_side (rhs), - last_solution(last_solution) {}; - - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - protected: - const Function &right_hand_side; - const Vector &last_solution; -}; - - - - - - -template -class NonlinearProblem : public ProblemBase { - public: - NonlinearProblem (); - void run (); - - protected: - Triangulation *tria; - DoFHandler *dof; - - Vector last_solution; -}; - - - - -template -void PoissonEquation::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &) const { - const vector > >&gradients = fe_values.get_shape_grads (); - const FullMatrix &values = fe_values.get_shape_values (); - vector rhs_values (fe_values.n_quadrature_points); - const vector &weights = fe_values.get_JxW_values (); - - vector > last_solution_grads(fe_values.n_quadrature_points); - fe_values.get_function_grads (last_solution, last_solution_grads); - - - right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values); - - for (unsigned int point=0; point -void PoissonEquation::assemble (FullMatrix &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, typename Equation::ExcPureVirtualFunctionCalled()); -}; - - - -template -void PoissonEquation::assemble (Vector &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, typename Equation::ExcPureVirtualFunctionCalled()); -}; - - - - -template -NonlinearProblem::NonlinearProblem () : - tria(0), dof(0) {}; - - - -template -void NonlinearProblem::run () { - - // first reset everything to a virgin state - clear (); - - tria = new Triangulation(); - dof = new DoFHandler (*tria); - set_tria_and_dof (tria, dof); - - - RightHandSide rhs; - ZeroFunction boundary_values; - StraightBoundary boundary; - - FEQ1 fe; - PoissonEquation equation (rhs, last_solution); - QGauss2 quadrature; - - FunctionMap::type dirichlet_bc; - dirichlet_bc[0] = &boundary_values; - - - GridGenerator::hyper_cube (*tria); - tria->refine_global (4); - - for (unsigned int refinement_step=0; refinement_step<5; ++refinement_step) - { - cout << "Refinement step " << refinement_step << endl - << " Grid has " << tria->n_active_cells() << " active cells." << endl; - - cout << " Distributing dofs... "; - dof->distribute_dofs (fe); - cout << dof->n_dofs() << " degrees of freedom." << endl; - - // set the starting values for the iteration - // to a constant value of 1 - last_solution.reinit (dof->n_dofs()); - for (unsigned int i=0; in_dofs(); ++i) - last_solution(i) = 1; - - - // here comes the fixed point iteration - for (unsigned int nonlinear_step=0; nonlinear_step<10; ++nonlinear_step) - { - cout << " Nonlinear step " << nonlinear_step << endl; - cout << " Assembling matrices..." << endl; - assemble (equation, quadrature, - UpdateFlags(update_values | update_gradients | - update_JxW_values | update_q_points), - dirichlet_bc); - - cout << " Solving..." << endl; - solve (); - - if (nonlinear_step % 2 == 0) - { - string filename = "nonlinear."; - filename += ('0' + refinement_step); - filename += '.'; - filename += ('0' + (nonlinear_step/2)); - filename += ".gnuplot"; - cout << " Writing to file <" << filename << ">..." << endl; - - DataOut out; - ofstream gnuplot(filename.c_str()); - fill_data (out); - out.build_patches (); - out.write_gnuplot (gnuplot); - gnuplot.close (); - }; - - last_solution = solution; - }; - - Vector error_indicator; - KellyErrorEstimator ee; - QSimpson eq; - ee.estimate (*dof, eq, - FunctionMap::type(), - solution, - error_indicator); - GridRefinement::refine_and_coarsen_fixed_number (*tria, error_indicator, 0.3, 0); - tria->execute_coarsening_and_refinement (); - }; - - - delete dof; - delete tria; - - cout << endl; -}; - - - - -int main () -{ - deallog.depth_console (0); - - NonlinearProblem<2> problem; - problem.run (); -}; diff --git a/tests/big-tests/poisson/.cvsignore b/tests/big-tests/poisson/.cvsignore deleted file mode 100644 index 572d49119c..0000000000 --- a/tests/big-tests/poisson/.cvsignore +++ /dev/null @@ -1,6 +0,0 @@ -poisson -Makefile.dep -*.go -Makefile.dep -*.go -*.o diff --git a/tests/big-tests/poisson/Makefile b/tests/big-tests/poisson/Makefile deleted file mode 100644 index 4613ca74f7..0000000000 --- a/tests/big-tests/poisson/Makefile +++ /dev/null @@ -1,136 +0,0 @@ -# $Id$ -# Copyright by the deal.II authors 1998, 1999, 2000 - -############################################################ -# $D contains the root of the deal distribution and should -# be supplied as a command line argument D=../../... -############################################################ - -# Template for makefiles for the examples subdirectory. In principle, -# everything should be done automatically if you set the target file -# here correctly: -target = poisson - -# 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 - -# If you want your program to be linked with extra object or library -# files, specify them here: -user-libs = - -# To run the program, use "make run"; to give parameters to the program, -# give the parameters to the following variable: -run-parameters = $(target).prm - -# To execute additional action apart from running the program, fill -# in this list: -additional-run-action = cd results ; gnuplot make_ps - -# To specify which files are to be deleted by "make clean" (apart from -# the usual ones: object files, executables, backups, etc), fill in the -# following list -delete-files = results/*gnuplot results/*.eps - - - - -############################################################################### -# Internals - -include $D/common/Make.global_options - - - -# get lists of files we need -cc-files = $(filter-out *%, $(shell echo *.cc)) -o-files = $(cc-files:.cc=.o) -go-files = $(cc-files:.cc=.go) -h-files = $(filter-out *%, $(shell echo *.h)) -lib-h-files = $(filter-out *%, $(shell echo $D/include/*/*.h)) - -# list of libraries needed to link with -libs.g = ./Obj.g.a \ - $(lib-deal2-2d.g) \ - $(lib-lac.g) \ - $(lib-base.g) -libs = ./Obj.a \ - $(lib-deal2-2d.o) \ - $(lib-lac.o) \ - $(lib-base.o) - - - -# check whether we use debug mode or not -ifeq ($(debug-mode),on) -libraries = $(libs.g) -flags = $(CXXFLAGS.g) -endif - -ifeq ($(debug-mode),off) -libraries = $(libs) -flags = $(CXXFLAGS) -endif - -# If in multithread mode, add the ACE library to the libraries which -# we need to link with: -ifneq ($(with-multithreading),no) - libraries += $(lib-ACE) -endif - - -# make rule for the target -$(target) : $(libraries) - @echo ============================ Linking $@ - @$(CXX) $(LDFLAGS) -o $@ $^ $(LIBS) $(user-libs) - -# rule how to run the program -run: $(target) - ./$(target) $(run-parameters) - $(additional-run-action) - - -# rule to make object files -%.go : %.cc - @echo ============================ Compiling with debugging information: $< - @echo $(CXX) ... -c $< -o $@ - @$(CXX) $(CXXFLAGS.g) -c $< -o $@ -%.o : %.cc - @echo ============================ Compiling with optimization: $< - @echo $(CXX) ... -c $< -o $@ - @$(CXX) $(CXXFLAGS.o) -c $< -o $@ - - -# rules which files the libraries depend upon -Obj.a: ./Obj.a($(o-files)) -Obj.g.a: ./Obj.g.a($(go-files)) - - -clean: - -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files) - - - -.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. -# -#use perl to generate rules for the .go files as well -#as to make rules not for tria.o and the like, but -#rather for libnumerics.a(tria.o) -Makefile.dep: $(cc-files) $(h-files) $(lib-h-files) - @echo ============================ Remaking Makefile - @perl $D/common/scripts/make_dependencies.pl ./Obj $(INCLUDE) $(cc-files) \ - > Makefile.dep - - -include Makefile.dep - diff --git a/tests/big-tests/poisson/equation.cc b/tests/big-tests/poisson/equation.cc deleted file mode 100644 index 08eda93399..0000000000 --- a/tests/big-tests/poisson/equation.cc +++ /dev/null @@ -1,86 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - -#include "poisson.h" -#include - - - -#if deal_II_dimension == 1 - -template <> -void PoissonEquation<1>::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues<1> &fe_values, - const DoFHandler<1>::cell_iterator &) const { - for (unsigned int point=0; point= 2 - -template -void PoissonEquation::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &) const { - const vector > >&gradients = fe_values.get_shape_grads (); - const FullMatrix &values = fe_values.get_shape_values (); - vector rhs_values (fe_values.n_quadrature_points); - const vector &weights = fe_values.get_JxW_values (); - - right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values); - - for (unsigned int point=0; point -void PoissonEquation::assemble (FullMatrix &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, typename Equation::ExcPureVirtualFunctionCalled()); -}; - - - -template -void PoissonEquation::assemble (Vector &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, typename Equation::ExcPureVirtualFunctionCalled()); -}; - - - - - - -template class PoissonEquation<2>; diff --git a/tests/big-tests/poisson/poisson.cc b/tests/big-tests/poisson/poisson.cc deleted file mode 100644 index 2e647e228a..0000000000 --- a/tests/big-tests/poisson/poisson.cc +++ /dev/null @@ -1,31 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - - -#include "poisson.h" -#include -#include - - - -int main (int argc, char **argv) { - if (argc!=2) - { - cout << "Usage: poisson parameterfile" << endl << endl; - return 1; - }; - - // no additional output to console - deallog.depth_console (0); - - PoissonProblem<2> poisson; - MultipleParameterLoop input_data; - - poisson.declare_parameters(input_data); - input_data.read_input (argv[1]); - input_data.loop (poisson); - - return 0; -}; - diff --git a/tests/big-tests/poisson/poisson.h b/tests/big-tests/poisson/poisson.h deleted file mode 100644 index 021381e394..0000000000 --- a/tests/big-tests/poisson/poisson.h +++ /dev/null @@ -1,100 +0,0 @@ -/*---------------------------- poisson.h ---------------------------*/ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ -#ifndef __poisson_H -#define __poisson_H -/*---------------------------- poisson.h ---------------------------*/ - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "../problem_base.h" -#include - - -#include -#include -#include -#include - - - - - - - -template -class PoissonEquation : public Equation { - public: - PoissonEquation (const Function &rhs) : - Equation(1), - right_hand_side (rhs) {}; - - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - protected: - const Function &right_hand_side; -}; - - - - - - -template -class PoissonProblem : public ProblemBase, - public MultipleParameterLoop::UserClass { - public: - PoissonProblem (); - virtual ~PoissonProblem(); - - void clear (); - - virtual void create_new (const unsigned int run_no); - virtual void declare_parameters (ParameterHandler &prm); - virtual void run (ParameterHandler &prm); - - - bool make_grid (ParameterHandler &prm); - void make_zoom_in_grid (); - void make_random_grid (); - - bool set_right_hand_side (ParameterHandler &prm); - bool set_boundary_values (ParameterHandler &prm); - - protected: - Triangulation *tria; - DoFHandler *dof; - - Function *rhs; - Function *boundary_values; - - Boundary *boundary; -}; - - - - - -/*---------------------------- poisson.h ---------------------------*/ -/* end of #ifndef __poisson_H */ -#endif -/*---------------------------- poisson.h ---------------------------*/ diff --git a/tests/big-tests/poisson/poisson.prm b/tests/big-tests/poisson/poisson.prm deleted file mode 100644 index 434c6ff65f..0000000000 --- a/tests/big-tests/poisson/poisson.prm +++ /dev/null @@ -1,5 +0,0 @@ -set Test run = { zoom in | ball | curved line | random | jump | L-region | slit domain} -set Global refinement = {{ 2 | 5 | 6 | 0 | 3 | 5 | 5 }} -set Right hand side = {{ zero | zero | trigpoly | constant | zero | zero | poly }} -set Boundary values = {{ sine | sine | zero | zero | jump | sine | sine }} -set Output file = results/{{ zoom_in | ball | curved_line | random | jump | L-region | slit_domain }}.gnuplot diff --git a/tests/big-tests/poisson/problem.cc b/tests/big-tests/poisson/problem.cc deleted file mode 100644 index 5135bb111b..0000000000 --- a/tests/big-tests/poisson/problem.cc +++ /dev/null @@ -1,622 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - - -#include "poisson.h" -#include -#include -#include -#include - - -template -class BoundaryValuesSine : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); - - double x = 1; - - for (unsigned int i=0; i &p, - Vector &values) const { - Assert (values.size()==1, ExcDimensionMismatch (values.size(), 1)); - - double x = 1; - - for (unsigned int i=0; i > &points, - vector &values, - const unsigned int component) const { - Assert (values.size() == points.size(), - ExcDimensionMismatch(values.size(), points.size())); - for (unsigned int i=0; i::value (points[i], component); - }; -}; - - - -template -class BoundaryValuesJump : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); - switch (dim) - { - case 1: - return 0; - default: - if (p(0) == p(1)) - return 0.5; - else - return (p(0)>p(1) ? 0. : 1.); - }; - }; -}; - - - - -template -class RHSTrigPoly : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int) const; -}; - - - -/** - Right hand side constructed such that the exact solution is - $x(1-x)$ in 1d, $x(1-x)*y(1-y)$ in 2d, etc. - */ -template -class RHSPoly : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int) const; -}; - - - - - - - - - -template -class CurvedLine : - public StraightBoundary { - public: - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - virtual Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; -}; - - - -template -Point -CurvedLine::get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - Point middle = StraightBoundary::get_new_point_on_line (line); - - // if the line is at the top of bottom - // face: do a special treatment on - // this line. Note that if the - // z-value of the midpoint is either - // 0 or 1, then the z-values of all - // vertices of the line is like that - if (dim>=3) - if (((middle(2) == 0) || (middle(2) == 1)) - // find out, if the line is in the - // interior of the top or bottom face - // of the domain, or at the edge. - // lines at the edge need to undergo - // the usual treatment, while for - // interior lines taking the midpoint - // is sufficient - // - // note: the trick with the boundary - // id was invented after the above was - // written, so we are not very strict - // here with using these flags - && (line->boundary_indicator() == 1)) - return middle; - - - double x=middle(0), - y=middle(1); - - if (y -Point -CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const -{ - Point middle = StraightBoundary::get_new_point_on_quad (quad); - - // if the face is at the top of bottom - // face: do not move the midpoint in - // x/y direction. Note that if the - // z-value of the midpoint is either - // 0 or 1, then the z-values of all - // vertices of the quad is like that - if ((middle(2) == 0) || (middle(2) == 1)) - return middle; - - double x=middle(0), - y=middle(1); - - if (y -double RHSTrigPoly::value (const Point &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); - - const double pi = 3.1415926536; - switch (dim) - { - case 1: - return p(0)*p(0)*cos(2*pi*p(0)); - case 2: - return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) + - 2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) + - 5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) - - 2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi); - default: - return 0; - }; -}; - - - -template -double RHSPoly::value (const Point &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); - - double ret_val = 0; - for (unsigned int i=0; i -PoissonProblem::PoissonProblem () : - tria(0), dof(0), rhs(0), boundary_values(0), boundary(0) {}; - - - -template -PoissonProblem::~PoissonProblem () -{ - clear (); -}; - - - -template -void PoissonProblem::clear () { - if (dof != 0) { - delete dof; - dof = 0; - }; - - if (boundary != 0) - { - tria->set_boundary (0); - delete boundary; - boundary = 0; - }; - - if (tria != 0) { - delete tria; - tria = 0; - }; - - // make it known to the underlying - // ProblemBase that tria and dof - // are already deleted - set_tria_and_dof (tria, dof); - - - if (rhs != 0) - { - delete rhs; - rhs = 0; - }; - - if (boundary_values != 0) - { - delete boundary_values; - boundary_values = 0; - }; - - ProblemBase::clear (); -}; - - - - -template -void PoissonProblem::create_new (const unsigned int) { - clear (); - - tria = new Triangulation(); - dof = new DoFHandler (*tria); - set_tria_and_dof (tria, dof); -}; - - - - -template -void PoissonProblem::declare_parameters (ParameterHandler &prm) { - if (dim>=2) - prm.declare_entry ("Test run", "zoom in", - Patterns::Selection("tensor|zoom in|ball|curved line|" - "random|jump|L-region|slit domain")); - else - prm.declare_entry ("Test run", "zoom in", - Patterns::Selection("tensor|zoom in|random")); - - prm.declare_entry ("Global refinement", "0", - Patterns::Integer()); - prm.declare_entry ("Right hand side", "zero", - Patterns::Selection("zero|constant|trigpoly|poly")); - prm.declare_entry ("Boundary values", "zero", - Patterns::Selection("zero|sine|jump")); - prm.declare_entry ("Output file", "gnuplot.1"); -}; - - - - -template -bool PoissonProblem::make_grid (ParameterHandler &prm) { - string test = prm.get ("Test run"); - unsigned int test_case; - if (test=="zoom in") test_case = 1; - else - if (test=="ball") test_case = 2; - else - if (test=="curved line") test_case = 3; - else - if (test=="random") test_case = 4; - else - if (test=="tensor") test_case = 5; - else - if (test=="jump") test_case = 6; - else - if (test=="L-region") test_case = 7; - else - if (test=="slit domain") test_case = 8; - else - { - cerr << "This test seems not to be implemented!" << endl; - return false; - }; - - switch (test_case) - { - case 1: - boundary = new StraightBoundary(); - tria->set_boundary (0, *boundary); - make_zoom_in_grid (); - break; - case 2: - // make ball grid around origin with - // unit radius - { - static const Point origin; - boundary = new HyperBallBoundary(origin, 1.); - GridGenerator::hyper_ball (*tria, origin, 1.); - tria->set_boundary (0, *boundary); - break; - }; - case 3: - // set the boundary function - { - boundary = new CurvedLine(); - GridGenerator::hyper_cube (*tria); - tria->set_boundary (0, *boundary); - break; - }; - case 4: - boundary = new StraightBoundary(); - tria->set_boundary (0, *boundary); - make_random_grid (); - break; - case 5: - boundary = new StraightBoundary(); - tria->set_boundary (0, *boundary); - GridGenerator::hyper_cube (*tria); - break; - case 6: - boundary = new StraightBoundary(); - tria->set_boundary (0, *boundary); - GridGenerator::hyper_cube (*tria); - tria->refine_global (1); - for (unsigned int i=0; i<5; ++i) - { - tria->begin_active(tria->n_levels()-1)->set_refine_flag(); - (--(tria->last_active()))->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - }; - break; - case 7: - boundary = new StraightBoundary(); - tria->set_boundary (0, *boundary); - GridGenerator::hyper_L (*tria); - break; - case 8: - boundary = new StraightBoundary(); - tria->set_boundary (0, *boundary); - GridGenerator::hyper_cube_slit (*tria); - break; - default: - return false; - }; - - int refine_global = prm.get_integer ("Global refinement"); - if ((refine_global < 0) || (refine_global>10)) - return false; - else - tria->refine_global (refine_global); - - return true; -}; - - - - -template -void PoissonProblem::make_zoom_in_grid () { - GridGenerator::hyper_cube (*tria); - - // refine first cell - tria->begin_active()->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - // refine first active cell - // on coarsest level - tria->begin_active()->set_refine_flag (); - tria->execute_coarsening_and_refinement (); - - Triangulation::active_cell_iterator cell; - for (int i=0; i<(dim==3 ? 5 : 17); ++i) - { - // refine the presently - // second last cell several - // times - cell = tria->last_active(tria->n_levels()-1); - --cell; - cell->set_refine_flag (); - tria->execute_coarsening_and_refinement (); - }; -}; - - - - -template -void PoissonProblem::make_random_grid () { - GridGenerator::hyper_cube (*tria); - tria->refine_global (1); - - Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<(dim==3 ? 7 : 12); ++i) - { - int n_levels = tria->n_levels(); - cell = tria->begin_active(); - endc = tria->end(); - - for (; cell!=endc; ++cell) - { - double r = rand()*1.0/RAND_MAX, - weight = 1.* - (cell->level()*cell->level()) / - (n_levels*n_levels); - - if (r <= 0.5*weight) - cell->set_refine_flag (); - }; - - tria->execute_coarsening_and_refinement (); - }; -}; - - - - -template -bool PoissonProblem::set_right_hand_side (ParameterHandler &prm) { - string rhs_name = prm.get ("Right hand side"); - - if (rhs_name == "zero") - rhs = new ZeroFunction(); - else - if (rhs_name == "constant") - rhs = new ConstantFunction(1.); - else - if (rhs_name == "trigpoly") - rhs = new RHSTrigPoly(); - else - if (rhs_name == "poly") - rhs = new RHSPoly (); - else - return false; - - if (rhs != 0) - return true; - else - return false; -}; - - - -template -bool PoissonProblem::set_boundary_values (ParameterHandler &prm) { - string bv_name = prm.get ("Boundary values"); - - if (bv_name == "zero") - boundary_values = new ZeroFunction (); - else - if (bv_name == "sine") - boundary_values = new BoundaryValuesSine (); - else - if (bv_name == "jump") - boundary_values = new BoundaryValuesJump (); - else - { - cout << "Unknown boundary value function " << bv_name << endl; - return false; - }; - - if (boundary_values != 0) - return true; - else - return false; -}; - - - - -template -void PoissonProblem::run (ParameterHandler &prm) { - cout << "Test case = " << prm.get ("Test run") - << endl; - - cout << " Making grid... "; - if (!make_grid (prm)) - return; - cout << tria->n_active_cells() << " active cells." << endl; - - if (!set_right_hand_side (prm)) - return; - - if (!set_boundary_values (prm)) - return; - - FEQ1 fe; - PoissonEquation equation (*rhs); - QGauss2 quadrature; - - cout << " Distributing dofs... "; - dof->distribute_dofs (fe); - cout << dof->n_dofs() << " degrees of freedom." << endl; - - cout << " Assembling matrices..." << endl; - ProblemBase::FunctionMap dirichlet_bc; - dirichlet_bc[0] = boundary_values; - assemble (equation, quadrature, - UpdateFlags(update_values | update_gradients | - update_JxW_values | update_q_points), - dirichlet_bc); - - cout << " Solving..." << endl; - solve (); - - cout << " Writing to file <" << prm.get("Output file") << ">..." - << endl; - - DataOut out; - string o_filename = prm.get ("Output file"); - ofstream gnuplot(o_filename.c_str()); - out.attach_dof_handler (*dof_handler); - out.add_data_vector (solution, "solution"); - out.build_patches (); - out.write_gnuplot (gnuplot); - gnuplot.close (); - - // release the lock of the DoF object to - // the FE object - dof->clear (); - - cout << endl; -}; - - - - - -template class PoissonProblem<2>; diff --git a/tests/big-tests/poisson/results/.cvsignore b/tests/big-tests/poisson/results/.cvsignore deleted file mode 100644 index acabd78095..0000000000 --- a/tests/big-tests/poisson/results/.cvsignore +++ /dev/null @@ -1,2 +0,0 @@ -*.gnuplot -*.eps diff --git a/tests/big-tests/problem_base.h b/tests/big-tests/problem_base.h deleted file mode 100644 index b32d72d46d..0000000000 --- a/tests/big-tests/problem_base.h +++ /dev/null @@ -1,857 +0,0 @@ -/*---------------------------- problem_base.h ---------------------------*/ -/* Previously deal.II/include/numerics/base.h and deal.II/source/numerics/base.cc, - but deprecated before deal.II version 3.0 and removed after version 3.1 - - Kept here to ensure that the example programs in big-tests still - compile. It is put into this directory to allow simple access from - all the subdirectories. -*/ - -#ifndef __deal2__base_h -#define __deal2__base_h - - - - -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -/** - * The use of this class is now deprecated! - * - * This is the base class for equation objects. Equations objects describe the - * finite element discretisation of one or more equations. - * - * Equation objects need only provide functions which set up the cell - * matrices and the cell right hand side. These are then automatically inserted - * into the global matrices and vectors. - * - * @author Wolfgang Bangerth, 1998 - */ -template -class Equation -{ - public: - /** - * Constructor. You have to pass the number - * of equations you want to discretize, which - * equals the number of solution functions. - */ - Equation (const unsigned int n_equations); - - /** - * Virtual function which assembles the - * cell matrix and the right hand side - * on a given cell. - * - * This function assumes the cell matrix - * and right hand side to have the right - * size and to be empty. Functions of - * derived classes should check for - * this. - * For that purpose, the two exceptions - * @p{ExcWrongSize} and @p{ExcObjectNotEmpty} - * are declared. - */ - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::cell_iterator &cell) const; - - /** - * Virtual function which only assembles - * the cell matrix on a given cell. - * - * This function assumes the cell matrix - * and right hand side to have the right - * size and to be empty. Functions of - * derived classes should check for - * this. - * For that purpose, the two exceptions - * @p{ExcWrongSize} and @p{ExcObjectNotEmpty} - * are declared. - */ - virtual void assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const typename DoFHandler::cell_iterator &cell) const; - - /** - * Virtual function which only assembles - * the right hand side on a given cell. - * - * This function assumes the cell matrix - * and right hand side to have the right - * size and to be empty. Functions of - * derived classes should check for - * this. - * For that purpose, the two exceptions - * @p{ExcWrongSize} and @p{ExcObjectNotEmpty} - * are declared. - */ - virtual void assemble (Vector &rhs, - const FEValues &fe_values, - const typename DoFHandler::cell_iterator &cell) const; - - /** - * Return number of equations for this - * equation object. This equals the number - * of solution functions. - */ - unsigned int n_equations () const; - - /** - * Exception - */ - DeclException0 (ExcPureVirtualFunctionCalled); - /** - * Exception - */ - DeclException2 (ExcWrongSize, - int, int, - << "Object has wrong size " << arg1 - << ", but should have " << arg2 << "."); - /** - * Exception - */ - DeclException0 (ExcObjectNotEmpty); - - protected: - /** - * Store the number of solution functions, - * which is the same as the number of - * equations. - */ - const unsigned int n_eq; -}; - - - -/** - * The use of this class is now deprecated! - * - * An @p{Assembler} is a specialized version of a @p{DoFCellAccessor} which adds - * functionality to assemble global matrices and vectors from cell base ones. - * - * @author Wolfgang Bangerth, 1998 - */ -template -class Assembler : public DoFCellAccessor -{ - public: - - /** - * Structure to be passed upon - * construction of an assembler object - * through the iterator object. See - * @ref{TriaRawIterator} for a discussion - * of this mechanism. - */ - struct AssemblerData { - /** - * Constructor. - */ - AssemblerData (const DoFHandler &dof, - const bool assemble_matrix, - const bool assemble_rhs, - SparseMatrix &matrix, - Vector &rhs_vector, - const Quadrature &quadrature, - const UpdateFlags &update_flags); - - /** - * Pointer to the dof handler object - * to be used to iterate on. - */ - const DoFHandler &dof; - - /** - * Flags to assemble the matrix. - */ - const bool assemble_matrix; - - /** - * Flags whether to assemble the right hand sides. - */ - const bool assemble_rhs; - - /** - * Pointer to the matrix to be assembled - * by this object. Elements are summed - * up by the assembler, so you may want - * to clear this object (set all entries - * to zero) before use. - */ - SparseMatrix &matrix; - - /** - * Pointer to the vector to be assembled - * by this object. Elements are summed - * up by the assembler, so you may want - * to clear this object (set all entries - * to zero) before use. - */ - Vector &rhs_vector; - - /** - * Pointer to a quadrature object to be - * used for this assemblage process. - */ - const Quadrature &quadrature; - - /** - * Store which of the fields of the - * FEValues object need to be reinitialized - * on each cell. - */ - const UpdateFlags update_flags; - }; - - - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. - */ - typedef AssemblerData AccessorData; - - /** - * Default constructor, unused thus not - * implemented. - */ - Assembler (); - - /** - * Constructor. The @p{local_data} - * argument is assumed to be a pointer - * to an @p{AssemblerData} object. The data - * is copied, so the object need not live - * longer than the constructor call. - */ - Assembler (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data); - - /** - * Assemble on the present cell using - * the given equation objectand the data - * passed to the constructor. The elements - * of the local matrix and right hand side - * are added to the global matrix and - * vector so you may want to clear the - * matrix before use. - */ - void assemble (const Equation &); - - /** - * Exception. - */ - DeclException0 (ExcNoAssemblingRequired); - /** - * Exception. - */ - DeclException0 (ExcInvalidData); - /** - * Exception. - */ - /** - * Exception. - */ - private: - /** - * Store a local cell matrix. - */ - FullMatrix cell_matrix; - - /** - * Right hand side local to cell. - */ - Vector cell_vector; - - /** - * Store whether to assemble the - * global matrix. - */ - bool assemble_matrix; - - /** - * Store whether to assemble the - * right hand side. - */ - bool assemble_rhs; - - /** - * Pointer to the matrix to be assembled - * by this object. - */ - SparseMatrix &matrix; - - /** - * Pointer to the vector to be assembled - * by this object. - */ - Vector &rhs_vector; - - /** - * The finite element evaluated at the - * quadrature points. - */ - FEValues fe_values; -}; - - -#include -#include -#include -#include -#include -#include - -template class DataOut; -template class Assembler; -template class Equation; - -/** - * The use of this class is now deprecated! - * - * Base class for user problems. This class stores the system matrix and right - * hand side vectors as well as a solution vector. It initiates the assemblage - * process of matrix and vectors and so on. - * - * This class is not extremely versatile as could certainly be. For example - * it presently only supports sparse matrices and has no multigrid features. - * However, all these things depend strongly on the problem and it seems - * best to implement many of these things yourself. Thus, this class is more - * a display of concept how to work with deal.II. - * - * - * @sect3{Assemblage} - * - * The @p{assemble} member function does the assemblage of the system matrix and - * the given number of right hand sides. It does the following steps: - * @begin{itemize} - * @item Initialize solution vector with zero entries. - * @item Create sparsity pattern of the system matrix and condense it with - * the constraints induced by hanging nodes. - * @item Initialize an assembler object. - * @item Loop over all cells and assemble matrix and vectors using the given - * quadrature formula and the equation object which contains the weak - * formulation of the equation. - * @item Apply Dirichlet boundary conditions. See the section on boundary - * conditions for more details. - * @item Condense the system matrix and right hand side with the constraints - * induced by hanging nodes. - * @end{itemize} - * - * The @p{assemble} function needs an object describing the boundary of the domain, - * since for higher order finite elements, we may be tempted to use curved faces - * of cells for better approximation of the boundary. In this case, the - * transformation from the unit cell to the real cell requires knowledge of - * the exact boundary of the domain. - * - * - * @sect3{Solving} - * - * Calling the @p{solve} function with a solver object, the system of equations - * which results after having called the @p{assemble} function is solved. After - * this, the solution vector is distributed again, i.e. the constrained nodes - * are given their correct values. - * - * - * @sect3{Boundary conditions} - * - * During assemblage of matrices and right hand side, use is made of dirichlet - * boundary conditions (in short: bc) specified to the @p{assemble} function. You - * can specify a list of pairs of boundary indicators (of type @p{unsigned char}; - * see the section in the documentation of the @ref{Triangulation} class for more - * details) and the according functions denoting the dirichlet boundary values - * of the nodes on boundary faces with this boundary indicator. - * - * To actually apply the boundary conditions, use is made of the - * @p{MatrixTools::apply_boundary_values} function and by interpolation of - * the @p{boundary_values} using the @p{MatrixTool::interpolate_boundary_values} - * function. See there for more information. - * - * - * @author Wolfgang Bangerth, 1998 - */ -template -class ProblemBase -{ - public: - /** - * Declare a data type which denotes a - * mapping between a boundary indicator - * and the function denoting the boundary - * values on this part of the boundary. - * Only one boundary function may be given - * for each boundary indicator, which is - * guaranteed by the @p{map} data type. - * - * See the general documentation of this - * class for more detail. - */ - typedef typename std::map*> FunctionMap; - /** - * Typdedef an iterator which assembles - * matrices and vectors. - */ - typedef TriaActiveIterator > active_assemble_iterator; - - /** - * Constructor. - */ - ProblemBase (); - - /** - * Destructor. Declare this only to have - * a virtual destructor, which is safer - * as we have virtual functions. - * It actually does nothing spectacular. - */ - virtual ~ProblemBase (); - - /** - * Use this triangulation and - * degree of freedom object during the - * lifetime of this object. The dof - * object must refer to the given - * triangulation. - */ - void set_tria_and_dof (Triangulation *tria, - DoFHandler *dof_handler); - - /** - * Reset all fields to a state as if we - * were right after calling the - * constructor. This is useful if you - * want to use an object derived from - * this base class for multiple - * successive calculations. In special, - * all aquired memory should be freed - * until it is needed again. - */ - void clear (); - - /** - * Initiate the process of assemblage of - * vectors and system matrix. Use the - * given equation object and the given - * quadrature formula. Also use the list - * of dirichlet boundary value functions - * (by default, no dirichlet bc are assumed - * which means that all bc are included - * into the weak formulation). - * - * For what exactly happens here, refer to - * the general doc of this class. - */ - virtual void assemble (const Equation &equation, - const Quadrature &q, - const UpdateFlags update_flags, - const FunctionMap &dirichlet_bc = FunctionMap()); - - /** - * Solve the system of equations. This uses - * a simple CG method. - */ - virtual void solve (); - - /** - * Initialize the @p{DataOut} object with - * the grid and DoF handler used for this - * computation, as well as with the - * solution. Overload this function if - * you have multiple data sets to be - * written, or alternativelt call this - * function and attach the additional - * vectors directly to the @p{DataOut} - * object. - * - * The solution name are - * derived by calling the virtual - * function @p{get_solution_name}. - */ - virtual void fill_data (DataOut &) const; - - - /** - * Return the name of the solution as a - * @p{string}. The default implementation - * returns @p{"solution"}. - * Overload this function, if you - * want anything else. - */ - virtual std::string get_solution_name () const; - - /** - * Exception - */ - DeclException0 (ExcDofAndTriaDontMatch); - /** - * Exception - */ - DeclException0 (ExcNoTriaSelected); - - protected: - /** - * Pointer to the triangulation to work on. - */ - Triangulation *tria; - - /** - * Pointer to the degree of freedom handler - * to work with. - */ - DoFHandler *dof_handler; - - /** - * Sparsity pattern of the system matrix. - */ - SparsityPattern system_sparsity; - - /** - * System matrix. - */ - SparseMatrix system_matrix; - - /** - * Vector storing the right hand side. - */ - Vector right_hand_side; - - /** - * Solution vector. - */ - Vector solution; - - /** - * List of constraints introduced by - * hanging nodes. - */ - ConstraintMatrix constraints; - - friend class Assembler; -}; - - -#include -#include -#include -#include -#include -#include -#include -#include - - -// if necessary try to work around a bug in the IBM xlC compiler -#ifdef XLC_WORK_AROUND_STD_BUG -using namespace std; -#endif - - -//TODO: purge this variable -static const MappingQ1<2> mapping; - -template -Assembler::AssemblerData::AssemblerData (const DoFHandler &dof, - const bool assemble_matrix, - const bool assemble_rhs, - SparseMatrix &matrix, - Vector &rhs_vector, - const Quadrature &quadrature, - const UpdateFlags &update_flags) : - dof(dof), - assemble_matrix(assemble_matrix), - assemble_rhs(assemble_rhs), - matrix(matrix), - rhs_vector(rhs_vector), - quadrature(quadrature), - update_flags(update_flags) -{}; - - -template -Assembler::Assembler (Triangulation *tria, - const int level, - const int index, - const AssemblerData *local_data) : - DoFCellAccessor (tria,level,index, &local_data->dof), - cell_matrix (dof_handler->get_fe().dofs_per_cell), - cell_vector (Vector(dof_handler->get_fe().dofs_per_cell)), - assemble_matrix (local_data->assemble_matrix), - assemble_rhs (local_data->assemble_rhs), - matrix(local_data->matrix), - rhs_vector(local_data->rhs_vector), - fe_values (mapping, dof_handler->get_fe(), - local_data->quadrature, - local_data->update_flags) -{ - Assert (!assemble_matrix || (matrix.m() == dof_handler->n_dofs()), - ExcInvalidData()); - Assert (!assemble_matrix || (matrix.n() == dof_handler->n_dofs()), - ExcInvalidData()); - Assert (!assemble_rhs || (rhs_vector.size()==dof_handler->n_dofs()), - ExcInvalidData()); -}; - - -template -void Assembler::assemble (const Equation &equation) { - // re-init fe values for this cell - fe_values.reinit (DoFHandler::cell_iterator (*this)); - const unsigned int n_dofs = dof_handler->get_fe().dofs_per_cell; - - if (assemble_matrix) - cell_matrix.clear (); - if (assemble_rhs) - cell_vector.clear (); - - -// fill cell matrix and vector if required - DoFHandler::cell_iterator this_cell (*this); - if (assemble_matrix && assemble_rhs) - equation.assemble (cell_matrix, cell_vector, fe_values, this_cell); - else - if (assemble_matrix) - equation.assemble (cell_matrix, fe_values, this_cell); - else - if (assemble_rhs) - equation.assemble (cell_vector, fe_values, this_cell); - else - Assert (false, ExcNoAssemblingRequired()); - - -// get indices of dofs - std::vector dofs (n_dofs); - get_dof_indices (dofs); - - // one could use the - // @p{distribute_local_to_global} functions - // here, but they would require getting the - // dof indices twice, so we leave it the - // way it was originally programmed. - - // distribute cell matrix - if (assemble_matrix) - for (unsigned int i=0; i -#include - -template -Equation::Equation (const unsigned int n_equations) : - n_eq(n_equations) {}; - - -template -void Equation::assemble (FullMatrix &, - Vector &, - const FEValues &, - const typename DoFHandler::cell_iterator &) const -{ - Assert (false, ExcPureVirtualFunctionCalled()); -}; - - -template -void Equation::assemble (FullMatrix &, - const FEValues &, - const typename DoFHandler::cell_iterator &) const -{ - Assert (false, ExcPureVirtualFunctionCalled()); -}; - - -template -void Equation::assemble (Vector &, - const FEValues &, - const typename DoFHandler::cell_iterator &) const -{ - Assert (false, ExcPureVirtualFunctionCalled()); -}; - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - - -template -ProblemBase::ProblemBase () : - tria(0), - dof_handler(0), - system_sparsity(), // dummy initialisation, is later reinit'd - system_matrix() // dummy initialisation, is later reinit'd -{}; - - -template -void ProblemBase::set_tria_and_dof (Triangulation *t, - DoFHandler *d) -{ - tria = t; - dof_handler = d; - - // allow a user to reset both tria and - // dof to null pointers, but don't allow - // something other - if ((tria != 0) || (dof_handler != 0)) - { - Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); - Assert (tria == &dof_handler->get_tria(), ExcDofAndTriaDontMatch()); - }; -}; - - -template -void ProblemBase::clear () { - if (tria) { delete tria; tria = 0; }; - if (dof_handler) { delete dof_handler; dof_handler = 0; }; - system_sparsity.reinit (0,0,1); - system_matrix.clear (); - right_hand_side.reinit (0); - solution.reinit (0); - constraints.clear (); -}; - - -template -ProblemBase::~ProblemBase () {}; - - -template -void ProblemBase::assemble (const Equation &equation, - const Quadrature &quadrature, - const UpdateFlags update_flags, - const FunctionMap &dirichlet_bc) -{ - Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); - - system_sparsity.reinit (dof_handler->n_dofs(), - dof_handler->n_dofs(), - dof_handler->max_couplings_between_dofs()); - right_hand_side.reinit (dof_handler->n_dofs()); - - // make up sparsity pattern and - // compress with constraints - constraints.clear (); - DoFTools::make_hanging_node_constraints (*dof_handler, constraints); - constraints.close (); - DoFTools::make_sparsity_pattern (*dof_handler, system_sparsity); - constraints.condense (system_sparsity); - - // reinite system matrix - system_matrix.reinit (system_sparsity); - // reinit solution vector, preset - // with zeroes. - solution.reinit (dof_handler->n_dofs()); - - // create assembler - Assembler::AssemblerData data (*dof_handler, - true, true, //assemble matrix and rhs - system_matrix, - right_hand_side, - quadrature, - update_flags); - active_assemble_iterator assembler (tria, - tria->begin_active()->level(), - tria->begin_active()->index(), - &data); - // loop over all cells, fill matrix and rhs - do - { - assembler->assemble (equation); - } - while ((++assembler).state() == valid); - - // condense system matrix in-place - constraints.condense (system_matrix); - - // condense right hand side in-place - constraints.condense (right_hand_side); - - // apply Dirichlet bc as described - // in the docs - std::map boundary_value_list; - - for (typename FunctionMap::const_iterator dirichlet = dirichlet_bc.begin() ; - dirichlet != dirichlet_bc.end() ; ++dirichlet) - VectorTools::interpolate_boundary_values (*dof_handler, - dirichlet->first, - *(dirichlet->second), - boundary_value_list); - MatrixTools::apply_boundary_values (boundary_value_list, - system_matrix, solution, - right_hand_side, - true); -}; - - -template -void ProblemBase::solve () -{ - Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); - - SolverControl control(4000, 1e-16); - PrimitiveVectorMemory<> memory; - SolverCG<> cg(control,memory); - - PreconditionSSOR<> prec; - prec.initialize (system_matrix, 1.2); - - // solve - cg.solve (system_matrix, solution, right_hand_side, prec); - // distribute solution - constraints.distribute (solution); -}; - - -template -void ProblemBase::fill_data (DataOut &out) const -{ - Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); - - out.clear (); - out.attach_dof_handler (*dof_handler); - - out.add_data_vector (solution, get_solution_name()); -}; - - -template -std::string ProblemBase::get_solution_name () const -{ - return "solution"; -}; - - -#endif - -- 2.39.5