+++ /dev/null
-# 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
-
-
+++ /dev/null
-# $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========= $(<F)
- @$(CXX) $(CXXFLAGS.g) -c $< -o $@
-%.o : %.cc
- @echo ==============optimized===== $(<F)
- @$(CXX) $(CXXFLAGS) -c $< -o $@
-
-
-# The following statement tells make that the rules `run' and `clean'
-# are not expected to produce files of the same name as Makefile rules
-# usually do.
-.PHONY: run clean
-
-
-# Finally there is a rule which you normally need not care much about:
-# since the executable depends on some include files from the library,
-# besides the C++ application file of course, it is necessary to
-# re-generate the executable when one of the files it depends on has
-# changed. The following rule to created a dependency file
-# `Makefile.dep', which `make' uses to determine when to regenerate
-# the executable. 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.
-#
-# The dependency file is created using a perl script. Since the
-# script prefixes the output names by `lib/o' or `lib/go' (it was
-# written for the sublibraries' Makefile), we have to strip that again
-# since object files are placed in the present directory for this
-# application. All these things are made in the next rule:
-Makefile.dep: $(target).cc Makefile \
- $(shell echo $(include-path-base)/base/*.h \
- $(include-path-lac)/lac/*.h \
- $(include-path-deal2)/*/*.h)
- @echo ============================ Remaking Makefile
- @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
- | perl -pi -e 's!lib/g?o/!!g;' \
- > Makefile.dep
-
-# To make the dependencies known to `make', we finally have to include
-# them:
-include Makefile.dep
-
-
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-#include <base/logstream.h>
-#include <grid/tria.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary_lib.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-#include <grid/grid_generator.h>
-#include <base/function.h>
-#include <numerics/data_out.h>
-#include <fe/fe_lib.lagrange.h>
-#include <base/quadrature_lib.h>
-#include "../problem_base.h"
-#include <numerics/vectors.h>
-#include <lac/vector.h>
-
-#include <map>
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <cstdlib>
-
-
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation (const Function<dim> &rhs) :
- Equation<dim>(1),
- right_hand_side (rhs) {};
-
- virtual void assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (FullMatrix<double> &cell_matrix,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- protected:
- const Function<dim> &right_hand_side;
-};
-
-
-
-
-
-
-template <int dim>
-class PoissonProblem : public ProblemBase<dim> {
- public:
- PoissonProblem (unsigned int order);
- ~PoissonProblem ();
-
- void clear ();
- void create_new ();
- int run (unsigned int level);
- void print_history (string filename) const;
-
- protected:
- Triangulation<dim> *tria;
- DoFHandler<dim> *dof;
-
- Function<dim> *rhs;
- Function<dim> *boundary_values;
-
- vector<double> l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error;
- vector<int> n_dofs;
-
- unsigned int order;
-};
-
-
-
-
-
-/**
- Right hand side constructed such that the exact solution is
- $sin(2 pi x) + sin(2 pi y)$
- */
-template <int dim>
-class RHSPoly : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
-};
-
-
-
-template <int dim>
-class Solution : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- /**
- * Return the gradient of the function
- * at the given point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &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<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<2> &fe_values,
- const DoFHandler<2>::cell_iterator &) const {
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point)) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side.value(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, ExcPureFunctionCalled());
-};
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (Vector<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, ExcPureFunctionCalled());
-};
-
-
-
-
-
-
-
-
-
-template <int dim>
-PoissonProblem<dim>::PoissonProblem (unsigned int order) :
- tria(0), dof(0), rhs(0),
- boundary_values(0), order(order) {};
-
-
-template <int dim>
-PoissonProblem<dim>::~PoissonProblem ()
-{
- clear ();
-};
-
-
-
-template <int dim>
-void PoissonProblem<dim>::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<dim>::clear ();
-};
-
-
-
-
-template <int dim>
-void PoissonProblem<dim>::create_new () {
- clear ();
-
- tria = new Triangulation<dim>();
- dof = new DoFHandler<dim> (*tria);
- set_tria_and_dof (tria, dof);
-};
-
-
-
-
-template <int dim>
-int PoissonProblem<dim>::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<dim> 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<dim>();
- boundary_values = new Solution<dim> ();
-
-
- FiniteElement<dim> *fe;
- PoissonEquation<dim> equation (*rhs);
- Quadrature<dim> *quadrature;
- Quadrature<dim-1> *boundary_quadrature;
- switch (order) {
- case 1:
- fe = new FEQ1<dim>();
- quadrature = new QGauss3<dim>();
- boundary_quadrature = new QGauss2<dim-1>();
- break;
- case 2:
- fe = new FEQ2<dim>();
- quadrature = new QGauss4<dim>();
- boundary_quadrature = new QGauss3<dim-1>();
- break;
- case 3:
- fe = new FEQ3<dim>();
- quadrature = new QGauss5<dim>();
- boundary_quadrature = new QGauss4<dim-1>();
- break;
- case 4:
- fe = new FEQ4<dim>();
- quadrature = new QGauss6<dim>();
- boundary_quadrature = new QGauss5<dim-1>();
- 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<dim>::FunctionMap dirichlet_bc;
- dirichlet_bc[0] = boundary_values;
- assemble (equation, *quadrature, update_flags, dirichlet_bc);
-
- cout << " Solving..." << endl;
- solve ();
-
- Solution<dim> sol;
- Vector<float> l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell;
- Vector<float> 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<double> l1_error_per_dof(dof->n_dofs());
- Vector<double> l2_error_per_dof(dof->n_dofs());
- Vector<double> linfty_error_per_dof(dof->n_dofs());
- Vector<double> h1_seminorm_error_per_dof(dof->n_dofs());
- Vector<double> 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<double> projected_solution;
-// ConstraintMatrix constraints;
-// constraints.close ();
-// VectorTools::project (*dof, constraints, *fe,
-// StraightBoundary<dim>(), *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<dim> 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 <int dim>
-void PoissonProblem<dim>::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; i<n_dofs.size(); ++i)
- out << n_dofs[i]
- << " "
- << l1_error[i] << " "
- << l2_error[i] << " "
- << linfty_error[i] << " "
- << h1_seminorm_error[i] << " "
- << h1_error[i] << endl;
-
- double average_l1=0,
- average_l2=0,
- average_linfty=0,
- average_h1_semi=0,
- average_h1=0;
- for (unsigned int i=1; i<n_dofs.size(); ++i)
- {
- average_l1 += l1_error[i]/l1_error[i-1];
- average_l2 += l2_error[i]/l2_error[i-1];
- average_linfty += linfty_error[i]/linfty_error[i-1];
- average_h1_semi += h1_seminorm_error[i]/h1_seminorm_error[i-1];
- average_h1 += h1_error[i]/h1_error[i-1];
- };
-
- average_l1 /= (l1_error.size()-1);
- average_l2 /= (l1_error.size()-1);
- average_linfty /= (l1_error.size()-1);
- average_h1_semi /= (l1_error.size()-1);
- average_h1 /= (l1_error.size()-1);
-
- cout << "==========================================================\n";
- cout << "Average error reduction rates for h->h/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;
-};
+++ /dev/null
-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"
-
-
+++ /dev/null
-# $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========= $(<F)
- @$(CXX) $(CXXFLAGS.g) -c $< -o $@
-%.o : %.cc
- @echo ==============optimized===== $(<F)
- @$(CXX) $(CXXFLAGS) -c $< -o $@
-
-
-# The following statement tells make that the rules `run' and `clean'
-# are not expected to produce files of the same name as Makefile rules
-# usually do.
-.PHONY: run clean
-
-
-# Finally there is a rule which you normally need not care much about:
-# since the executable depends on some include files from the library,
-# besides the C++ application file of course, it is necessary to
-# re-generate the executable when one of the files it depends on has
-# changed. The following rule to created a dependency file
-# `Makefile.dep', which `make' uses to determine when to regenerate
-# the executable. 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.
-#
-# The dependency file is created using a perl script. Since the
-# script prefixes the output names by `lib/o' or `lib/go' (it was
-# written for the sublibraries' Makefile), we have to strip that again
-# since object files are placed in the present directory for this
-# application. All these things are made in the next rule:
-Makefile.dep: $(target).cc Makefile \
- $(shell echo $(include-path-base)/base/*.h \
- $(include-path-lac)/lac/*.h \
- $(include-path-deal2)/*/*.h)
- @echo ============================ Remaking Makefile
- @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
- | perl -pi -e 's!lib/g?o/!!g;' \
- > Makefile.dep
-
-# To make the dependencies known to `make', we finally have to include
-# them:
-include Makefile.dep
-
-
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-
-#include <grid/grid_out.h>
-#include <dofs/dof_tools.h>
-#include <dofs/dof_handler.h>
-#include <grid/tria.h>
-#include <fe/fe_lib.lagrange.h>
-#include <grid/tria_boundary.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_accessor.h>
-#include <grid/grid_generator.h>
-#include <lac/sparsity_pattern.h>
-#include <base/parameter_handler.h>
-#include <dofs/dof_constraints.h>
-#include <numerics/dof_renumbering.h>
-
-#include <fstream>
-#include <cmath>
-#include <cstdlib>
-
-
-
-// 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 <int dim>
-class Ball :
- public StraightBoundary<dim> {
- public:
- virtual Point<dim>
- get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const {
- Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line(line);
-
- for (int i=0; i<dim; ++i)
- middle(i) -= .5;
- middle *= sqrt(dim) / (sqrt(middle.square())*2);
- for (int i=0; i<dim; ++i)
- middle(i) += .5;
-
- return middle;
- };
-
-
- virtual Point<dim>
- get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const {
- Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad(quad);
-
- for (int i=0; i<dim; ++i)
- middle(i) -= .5;
- middle *= sqrt(dim) / (sqrt(middle.square())*2);
- for (int i=0; i<dim; ++i)
- middle(i) += .5;
-
- return middle;
- };
-};
-
-
-template <int dim>
-class CurvedLine :
- public StraightBoundary<dim> {
- public:
- virtual Point<dim>
- get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
-
- virtual Point<dim>
- get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
-};
-
-
-
-template <int dim>
-Point<dim>
-CurvedLine<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-{
- Point<dim> middle = StraightBoundary<dim>::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<x)
- if (y<1-x)
- middle(1) = 0.04*sin(6*3.141592*middle(0));
- else
- middle(0) = 1+0.04*sin(6*3.141592*middle(1));
-
- else
- if (y<1-x)
- middle(0) = 0.04*sin(6*3.141592*middle(1));
- else
- middle(1) = 1+0.04*sin(6*3.141592*middle(0));
-
- return middle;
-};
-
-
-
-template <int dim>
-Point<dim>
-CurvedLine<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-{
- Point<dim> middle = StraightBoundary<dim>::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<x)
- if (y<1-x)
- middle(1) = 0.04*sin(6*3.141592*middle(0));
- else
- middle(0) = 1+0.04*sin(6*3.141592*middle(1));
-
- else
- if (y<1-x)
- middle(0) = 0.04*sin(6*3.141592*middle(1));
- else
- middle(1) = 1+0.04*sin(6*3.141592*middle(0));
-
- return middle;
-};
-
-
-
-
-template <int dim>
-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<dim> *tria;
- DoFHandler<dim> *dof;
-};
-
-
-
-template <int dim>
-TestCases<dim>::TestCases () :
- tria(0), dof(0) {};
-
-
-template <int dim>
-TestCases<dim>::~TestCases ()
-{
- if (dof) delete dof;
- if (tria) delete tria;
-};
-
-
-
-template <int dim>
-void TestCases<dim>::create_new (const unsigned int) {
- if (dof != 0) delete dof;
- if (tria != 0) delete tria;
-
- tria = new Triangulation<dim>();
- GridGenerator::hyper_cube(*tria);
-
- dof = new DoFHandler<dim> (*tria);
-};
-
-
-
-template <int dim>
-void TestCases<dim>::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 <int dim>
-void TestCases<dim>::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<dim> *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<dim>::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<Boundary<dim>*>(new Ball<dim>()) :
- static_cast<Boundary<dim>*>(new CurvedLine<dim>()));
- tria->set_boundary (0, *boundary);
- tria->set_boundary (1, *boundary);
-
- // refine once
- tria->begin_active()->set_refine_flag();
- tria->execute_coarsening_and_refinement ();
-
- Triangulation<dim>::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<dim>::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<dim> 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;
-};
-
+++ /dev/null
-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
+++ /dev/null
-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"
-
+++ /dev/null
-# $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========= $(<F)
- @$(CXX) $(CXXFLAGS.g) -c $< -o $@
-%.o : %.cc
- @echo ==============optimized===== $(<F)
- @$(CXX) $(CXXFLAGS) -c $< -o $@
-
-
-# The following statement tells make that the rules `run' and `clean'
-# are not expected to produce files of the same name as Makefile rules
-# usually do.
-.PHONY: run clean
-
-
-# Finally there is a rule which you normally need not care much about:
-# since the executable depends on some include files from the library,
-# besides the C++ application file of course, it is necessary to
-# re-generate the executable when one of the files it depends on has
-# changed. The following rule to created a dependency file
-# `Makefile.dep', which `make' uses to determine when to regenerate
-# the executable. 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.
-#
-# The dependency file is created using a perl script. Since the
-# script prefixes the output names by `lib/o' or `lib/go' (it was
-# written for the sublibraries' Makefile), we have to strip that again
-# since object files are placed in the present directory for this
-# application. All these things are made in the next rule:
-Makefile.dep: $(target).cc Makefile \
- $(shell echo $(include-path-base)/base/*.h \
- $(include-path-lac)/lac/*.h \
- $(include-path-deal2)/*/*.h)
- @echo ============================ Remaking Makefile
- @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
- | perl -pi -e 's!lib/g?o/!!g;' \
- > Makefile.dep
-
-# To make the dependencies known to `make', we finally have to include
-# them:
-include Makefile.dep
-
-
+++ /dev/null
-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
+++ /dev/null
-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
+++ /dev/null
-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
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-#include <base/function.h>
-#include <base/parameter_handler.h>
-#include <base/quadrature_lib.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_out.h>
-#include <grid/tria.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_boundary_lib.h>
-#include <grid/tria_iterator.h>
-#include <grid/grid_refinement.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_tools.h>
-#include <fe/fe_lib.lagrange.h>
-#include <numerics/data_out.h>
-#include "../problem_base.h"
-#include <numerics/vectors.h>
-#include <numerics/error_estimator.h>
-#include <numerics/solution_transfer.h>
-#include <lac/vector.h>
-
-#include <map>
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <cstdlib>
-
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation (const Function<dim> &rhs) :
- Equation<dim>(1),
- use_coefficient(false),
- right_hand_side (rhs),
- coefficient (default_coefficient) {};
-
- PoissonEquation (const Function<dim> &rhs,
- const Function<dim> &coefficient ) :
- Equation<dim>(1),
- use_coefficient(true),
- right_hand_side (rhs),
- coefficient (coefficient) {};
-
- virtual void assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (FullMatrix<double> &cell_matrix,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- protected:
- const bool use_coefficient;
- const Function<dim> &right_hand_side;
- const Function<dim> &coefficient;
-
- static const ConstantFunction<dim> default_coefficient;
-};
-
-
-const ConstantFunction<2> PoissonEquation<2>::default_coefficient(1);
-
-
-
-
-
-template <int dim>
-class PoissonProblem : public ProblemBase<dim>, 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<dim> *tria;
- DoFHandler<dim> *dof;
-
- Function<dim> *rhs;
- Function<dim> *solution_function;
- Function<dim> *coefficient;
-
- Boundary<dim> *boundary;
-
- vector<double> l2_error, linfty_error;
- vector<double> h1_error, estimated_error;
- vector<int> n_dofs;
-};
-
-
-
-
-
-template <int dim>
-class Solution {
- public:
-
- class GaussShape : public Function<dim> {
- public:
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component) const;
- };
-
- class Singular : public Function<dim> {
- public:
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component) const;
- };
-
- class Kink : public Function<dim> {
- public:
- class Coefficient : public Function<dim> {
- public:
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- };
-
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component) const;
- };
-};
-
-
-
-
-template <int dim>
-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<dim> {
- public:
- virtual double value (const Point<dim> &p,
- const unsigned int component) const;
- };
-
- /**
- * Right hand side constructed such that
- * the exact solution is
- * $r^{2/3}$.
- */
- class Singular : public Function<dim> {
- public:
- virtual double value (const Point<dim> &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<dim> {
- public:
- virtual double value (const Point<dim> &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<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<2> &fe_values,
- const DoFHandler<2>::cell_iterator &) const {
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- {
- const double c = (use_coefficient ?
- coefficient.value(fe_values.quadrature_point(point)) :
- 1);
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point)) *
- fe_values.JxW(point) *
- c;
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side.value(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
- };
-};
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, typename Equation<dim>::ExcPureVirtualFunctionCalled());
-};
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (Vector<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, typename Equation<dim>::ExcPureVirtualFunctionCalled());
-};
-
-
-
-
-
-
-
-
-
-template <int dim>
-PoissonProblem<dim>::PoissonProblem () :
- tria(0), dof(0), rhs(0),
- solution_function(0), coefficient(0),
- boundary(0) {};
-
-
-
-template <int dim>
-PoissonProblem<dim>::~PoissonProblem ()
-{
- clear ();
-};
-
-
-
-template <int dim>
-void PoissonProblem<dim>::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<dim>::clear ();
-};
-
-
-
-
-template <int dim>
-void PoissonProblem<dim>::create_new (const unsigned int) {
- clear ();
-
- tria = new Triangulation<dim>();
- dof = new DoFHandler<dim> (*tria);
- set_tria_and_dof (tria, dof);
- boundary = new HyperBallBoundary<dim> ();
-};
-
-
-
-template <int dim>
-void PoissonProblem<dim>::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 <int dim>
-void PoissonProblem<dim>::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<dim>::GaussShape();
- else
- if (prm.get("Test case")=="Singular")
- rhs = new RHS<dim>::Singular();
- else
- if (prm.get("Test case")=="Kink")
- rhs = new RHS<dim>::Kink();
-
- if (prm.get("Test case")=="Gauss shape")
- solution_function = new Solution<dim>::GaussShape ();
- else
- if (prm.get("Test case")=="Singular")
- solution_function = new Solution<dim>::Singular ();
- else
- if (prm.get("Test case")=="Kink")
- solution_function = new Solution<dim>::Kink ();
-
-
- FEQ1<dim> fe;
- QGauss3<dim> quadrature;
- PoissonEquation<dim> *equation;
-
- static Solution<dim>::Kink::Coefficient kink_coefficient;
- if (prm.get("Test case")=="Kink")
- equation = new PoissonEquation<dim>(*rhs, kink_coefficient);
- else
- equation = new PoissonEquation<dim>(*rhs);
-
- SolutionTransfer<dim,double> 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<double> 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<dim>::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<float> l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell;
- Vector<float> estimated_error_per_cell;
- QGauss3<dim> 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<dim-1> eq;
- KellyErrorEstimator<dim>::estimate (*dof, eq,
- FunctionMap<dim>::type(),
- solution,
- estimated_error_per_cell,
- vector<bool>(), // 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<double> l2_error_per_dof(dof->n_dofs()), linfty_error_per_dof(dof->n_dofs());
- Vector<double> h1_error_per_dof(dof->n_dofs()), estimated_error_per_dof(dof->n_dofs());
- Vector<double> 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<dim> 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 <int dim>
-void PoissonProblem<dim>::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; i<n_dofs.size(); ++i)
- out << n_dofs[i]
- << " "
- << l2_error[i] << " "
- << linfty_error[i] << " "
- << h1_error[i] << " "
- << estimated_error[i] << " "
- << endl;
-
- double average_l2=0,
- average_linfty=0,
- average_h1=0,
- average_est=0;
-
- for (unsigned int i=1; i<n_dofs.size(); ++i)
- {
- average_l2 += l2_error[i]/l2_error[i-1];
- average_linfty += linfty_error[i]/linfty_error[i-1];
- average_h1 += h1_error[i]/h1_error[i-1];
- average_est += estimated_error[i]/estimated_error[i-1];
- };
-
- average_l2 /= (l2_error.size()-1);
- average_linfty /= (l2_error.size()-1);
- average_h1 /= (l2_error.size()-1);
- average_est /= (l2_error.size()-1);
-
- cout << "Average error reduction rates for h->h/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;
-};
-
-
-
+++ /dev/null
-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"
+++ /dev/null
-perl -pi -e 's/^#.*$\\n//g' data-*/*.inp
+++ /dev/null
-# $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========= $(<F)
- @$(CXX) $(CXXFLAGS.g) -c $< -o $@
-%.o : %.cc
- @echo ==============optimized===== $(<F)
- @$(CXX) $(CXXFLAGS) -c $< -o $@
-
-
-# The following statement tells make that the rules `run' and `clean'
-# are not expected to produce files of the same name as Makefile rules
-# usually do.
-.PHONY: run clean
-
-
-# Finally there is a rule which you normally need not care much about:
-# since the executable depends on some include files from the library,
-# besides the C++ application file of course, it is necessary to
-# re-generate the executable when one of the files it depends on has
-# changed. The following rule to created a dependency file
-# `Makefile.dep', which `make' uses to determine when to regenerate
-# the executable. 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.
-#
-# The dependency file is created using a perl script. Since the
-# script prefixes the output names by `lib/o' or `lib/go' (it was
-# written for the sublibraries' Makefile), we have to strip that again
-# since object files are placed in the present directory for this
-# application. All these things are made in the next rule:
-Makefile.dep: $(target).cc Makefile \
- $(shell echo $(include-path-base)/base/*.h \
- $(include-path-lac)/lac/*.h \
- $(include-path-deal2)/*/*.h)
- @echo ============================ Remaking Makefile
- @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
- | perl -pi -e 's!lib/g?o/!!g;' \
- > Makefile.dep
-
-# To make the dependencies known to `make', we finally have to include
-# them:
-include Makefile.dep
-
-
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-#include <grid/tria_boundary.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_out.h>
-#include <fstream>
-#include <string>
-#include <cmath>
-#include <cstdlib>
-
-
-
-// 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 <int dim>
-class Ball :
- public StraightBoundary<dim> {
- public:
- virtual Point<dim>
- get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const {
- Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line(line);
-
- for (int i=0; i<dim; ++i)
- middle(i) -= .5;
- middle *= sqrt(dim) / (sqrt(middle.square())*2);
- for (int i=0; i<dim; ++i)
- middle(i) += .5;
-
- return middle;
- };
-
-
- virtual Point<dim>
- get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const {
- Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad(quad);
-
- for (int i=0; i<dim; ++i)
- middle(i) -= .5;
- middle *= sqrt(dim) / (sqrt(middle.square())*2);
- for (int i=0; i<dim; ++i)
- middle(i) += .5;
-
- return middle;
- };
-};
-
-
-template <int dim>
-class CurvedLine :
- public StraightBoundary<dim> {
- public:
- virtual Point<dim>
- get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
-
- virtual Point<dim>
- get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
-};
-
-
-
-template <int dim>
-Point<dim>
-CurvedLine<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-{
- Point<dim> middle = StraightBoundary<dim>::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<x)
- if (y<1-x)
- middle(1) = 0.04*sin(6*3.141592*middle(0));
- else
- middle(0) = 1+0.04*sin(6*3.141592*middle(1));
-
- else
- if (y<1-x)
- middle(0) = 0.04*sin(6*3.141592*middle(1));
- else
- middle(1) = 1+0.04*sin(6*3.141592*middle(0));
-
- return middle;
-};
-
-
-
-template <int dim>
-Point<dim>
-CurvedLine<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-{
- Point<dim> middle = StraightBoundary<dim>::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<x)
- if (y<1-x)
- middle(1) = 0.04*sin(6*3.141592*middle(0));
- else
- middle(0) = 1+0.04*sin(6*3.141592*middle(1));
-
- else
- if (y<1-x)
- middle(0) = 0.04*sin(6*3.141592*middle(1));
- else
- middle(1) = 1+0.04*sin(6*3.141592*middle(0));
-
- return middle;
-};
-
-
-
-template <int dim>
-void test (const int test_case) {
- cout << "Running testcase " << test_case
- << " in " << dim << " dimensions." << endl;
- Triangulation<dim> 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<dim>::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<dim> ball;
- CurvedLine<dim> 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<dim>::active_cell_iterator cell, endc;
- const unsigned int steps[4] = { 0, 10, 7, 2 };
- for (unsigned int i=0; i<steps[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();
- };
-
- 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<dim>::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;
-};
+++ /dev/null
-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 <results/2d.1.eps>"
-set output "results/2d.1.eps"
-plot "results/2d.1"
-
-!echo " Making <results/2d.2.eps>"
-set output "results/2d.2.eps"
-plot "results/2d.2"
-
-!echo " Making <results/2d.3.eps>"
-set output "results/2d.3.eps"
-plot "results/2d.3"
-
-!echo " Making <results/2d.4.eps>"
-set output "results/2d.4.eps"
-plot "results/2d.4"
-
-
-!echo " Making <results/3d.1.eps>"
-set output "results/3d.1.eps"
-splot "results/3d.1"
-
-!echo " Making <results/3d.2.eps>"
-set output "results/3d.2.eps"
-splot "results/3d.2"
-
-!echo " Making <results/3d.3.eps>"
-set output "results/3d.3.eps"
-splot "results/3d.3"
-
-!echo " Making <results/3d.4.eps>"
-set output "results/3d.4.eps"
-splot "results/3d.4"
+++ /dev/null
-# $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========= $(<F)
- @$(CXX) $(CXXFLAGS.g) -c $< -o $@
-%.o : %.cc
- @echo ==============optimized===== $(<F)
- @$(CXX) $(CXXFLAGS) -c $< -o $@
-
-
-# The following statement tells make that the rules `run' and `clean'
-# are not expected to produce files of the same name as Makefile rules
-# usually do.
-.PHONY: run clean
-
-
-# Finally there is a rule which you normally need not care much about:
-# since the executable depends on some include files from the library,
-# besides the C++ application file of course, it is necessary to
-# re-generate the executable when one of the files it depends on has
-# changed. The following rule to created a dependency file
-# `Makefile.dep', which `make' uses to determine when to regenerate
-# the executable. 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.
-#
-# The dependency file is created using a perl script. Since the
-# script prefixes the output names by `lib/o' or `lib/go' (it was
-# written for the sublibraries' Makefile), we have to strip that again
-# since object files are placed in the present directory for this
-# application. All these things are made in the next rule:
-Makefile.dep: $(target).cc Makefile \
- $(shell echo $(include-path-base)/base/*.h \
- $(include-path-lac)/lac/*.h \
- $(include-path-deal2)/*/*.h)
- @echo ============================ Remaking Makefile
- @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
- | perl -pi -e 's!lib/g?o/!!g;' \
- > Makefile.dep
-
-# To make the dependencies known to `make', we finally have to include
-# them:
-include Makefile.dep
-
-
+++ /dev/null
-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"
-
-
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-
-#include <base/quadrature_lib.h>
-#include <base/function.h>
-#include <base/logstream.h>
-#include <lac/vector.h>
-#include <lac/full_matrix.h>
-#include <lac/sparse_matrix.h>
-#include <lac/solver_cg.h>
-#include <lac/vector_memory.h>
-#include <lac/precondition.h>
-#include <grid/tria.h>
-#include <grid/grid_refinement.h>
-#include <dofs/dof_handler.h>
-#include <grid/grid_generator.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary_lib.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_tools.h>
-#include <fe/fe_values.h>
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/data_out.h>
-#include <fe/fe_lib.lagrange.h>
-#include <grid/grid_out.h>
-#include <dofs/dof_constraints.h>
-#include <numerics/error_estimator.h>
-
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
-#include <multigrid/mg_dof_tools.h>
-#include <multigrid/mg_base.h>
-#include <multigrid/mg_smoother.h>
-#include <multigrid/multigrid.h>
-
-#include <lac/solver_richardson.h>
-
-#include <fstream>
-
-
-
-template <int dim>
-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<dim> triangulation;
- MGDoFHandler<dim> mg_dof_handler;
-
- FEQ1<dim> fe;
-
- ConstraintMatrix hanging_node_constraints;
-
- SparsityPattern global_sparsity_pattern;
- SparseMatrix<double> global_system_matrix;
-
- MGLevelObject<SparsityPattern> level_sparsity_patterns;
- MGLevelObject<SparseMatrix<double> > level_system_matrices;
-
- Vector<double> solution;
- Vector<double> system_rhs;
-};
-
-
-
-template <int dim>
-class Coefficient : public Function<dim>
-{
- public:
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- virtual void value_list (const vector<Point<dim> > &points,
- vector<double> &values,
- const unsigned int component = 0) const;
-};
-
-
-
-template <int dim>
-double Coefficient<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- if (p.square() < 0.5*0.5)
- return 20;
- else
- return 1;
-};
-
-
-
-template <int dim>
-void Coefficient<dim>::value_list (const vector<Point<dim> > &points,
- vector<double> &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<n_points; ++i)
- {
- if (points[i].square() < 0.5*0.5)
- values[i] = 20;
- else
- values[i] = 1;
- };
-};
-
-
-
-
-class MGSmootherLAC : public MGSmootherBase
-{
- private:
- SmartPointer<MGLevelObject<SparseMatrix<double> > >matrices;
- public:
- MGSmootherLAC(MGLevelObject<SparseMatrix<double> >&);
-
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const;
-};
-
-
-MGSmootherLAC::MGSmootherLAC(MGLevelObject<SparseMatrix<double> >& matrix)
- :
- matrices(&matrix)
-{}
-
-
-void
-MGSmootherLAC::smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &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 <int dim>
-LaplaceProblem<dim>::LaplaceProblem () :
- mg_dof_handler (triangulation)
-{};
-
-
-
-template <int dim>
-LaplaceProblem<dim>::~LaplaceProblem ()
-{
- mg_dof_handler.clear ();
-};
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::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<dim>::n_dofs(),
- mg_dof_handler.DoFHandler<dim>::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<dim>::n_dofs());
- system_rhs.reinit (mg_dof_handler.DoFHandler<dim>::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<n_levels; ++level)
- {
- level_sparsity_patterns[level].reinit (mg_dof_handler.n_dofs(level),
- mg_dof_handler.n_dofs(level),
- mg_dof_handler.max_couplings_between_dofs()); //xxx
- MGDoFTools::make_sparsity_pattern (mg_dof_handler,
- level_sparsity_patterns[level],
- level);
- level_sparsity_patterns[level].compress();
-
- level_system_matrices[level].reinit (level_sparsity_patterns[level]);
- };
-};
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
- const Coefficient<dim> coefficient;
-
- QGauss2<dim> quadrature_formula;
-
- FEValues<dim> 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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs (dofs_per_cell);
-
- vector<unsigned int> local_dof_indices (dofs_per_cell);
-
- // FIX
- vector<double> coefficient_values (n_q_points, 1.0);
-
- // not only active cells
- MGDoFHandler<dim>::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<double>
- & shape_values = fe_values.get_shape_values();
- const vector<vector<Tensor<1,dim> > >
- & shape_grads = fe_values.get_shape_grads();
- const vector<double>
- & JxW_values = fe_values.get_JxW_values();
- const vector<Point<dim> >
- & q_points = fe_values.get_quadrature_points();
-
- // FIX
-// coefficient.value_list (q_points, coefficient_values);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (coefficient_values[q_point] *
- (shape_grads[i][q_point] *
- shape_grads[j][q_point] +
- shape_values(i,q_point) *
- shape_values(j,q_point) ) *
- JxW_values[q_point]);
-
- cell_rhs(i) += (shape_values (i,q_point) *
- sin(4*sqrt(q_points[q_point].square())) *
- fe_values.JxW (q_point));
- };
-
-
- cell->get_mg_dof_indices (local_dof_indices);
- const unsigned int level = cell->level();
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- level_system_matrices[level].add (local_dof_indices[i],
- local_dof_indices[j],
- cell_matrix(i,j));
-
- // if active, then also into
- // global matrix
- if (cell->active())
- {
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- global_system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- cell_matrix(i,j));
-
- system_rhs(local_dof_indices[i]) += cell_rhs(i);
- };
- };
- };
-
- hanging_node_constraints.condense (global_system_matrix);
- hanging_node_constraints.condense (system_rhs);
-
-// map<unsigned int,double> boundary_values;
-// VectorTools::interpolate_boundary_values (mg_dof_handler,
-// 0,
-// ZeroFunction<dim>(),
-// boundary_values);
-// MatrixTools<dim>::apply_boundary_values (boundary_values,
-// global_system_matrix,
-// solution,
-// system_rhs);
-};
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::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<double>::template precondition_SSOR<double>,
-// 1.2);
- PreconditionIdentity coarse_grid_solver_preconditioner;
-
- MGCoarseGridLACIteration<SolverCG<>, SparseMatrix<double>, 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<Multigrid<2> >
- 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 <int dim>
-void LaplaceProblem<dim>::refine_grid ()
-{
- Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
- FunctionMap<dim>::type neumann_boundary;
- KellyErrorEstimator<dim>::estimate (mg_dof_handler,
- QGauss3<dim-1>(),
- 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 <int dim>
-void LaplaceProblem<dim>::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 <int dim>
-void LaplaceProblem<dim>::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<dim>::n_dofs()
- << endl;
-
- assemble_system ();
- solve ();
- output_results (cycle);
-
- typename DataOut<dim>::EpsFlags eps_flags;
- eps_flags.z_scaling = 4;
-
- DataOut<dim> 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;
-};
+++ /dev/null
-# $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
-
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-
-
-#include <lac/vector.h>
-#include <grid/tria.h>
-#include <dofs/dof_handler.h>
-#include <grid/tria_accessor.h>
-#include <dofs/dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_refinement.h>
-#include <numerics/data_out.h>
-#include <base/function.h>
-#include <fe/fe_lib.lagrange.h>
-#include <base/quadrature_lib.h>
-#include "../../problem_base.h"
-#include <numerics/error_estimator.h>
-#include <base/logstream.h>
-
-#include <map>
-#include <fstream>
-#include <cmath>
-#include <cstdlib>
-
-
-
-
-
-template <int dim>
-class RightHandSide : public Function<dim>
-{
- public:
- virtual double value (const Point<dim> &p,
- const unsigned int) const
- {
- double x = 80;
- for (unsigned int d=0; d<dim; ++d)
- if (p(d) < 0.5)
- x *= -p(d);
- else
- x *= (1-p(d));
-
- return x;
- };
-};
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation (const Function<dim> &rhs,
- const Vector<double> &last_solution) :
- Equation<dim>(1),
- right_hand_side (rhs),
- last_solution(last_solution) {};
-
- virtual void assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (FullMatrix<double> &cell_matrix,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- protected:
- const Function<dim> &right_hand_side;
- const Vector<double> &last_solution;
-};
-
-
-
-
-
-
-template <int dim>
-class NonlinearProblem : public ProblemBase<dim> {
- public:
- NonlinearProblem ();
- void run ();
-
- protected:
- Triangulation<dim> *tria;
- DoFHandler<dim> *dof;
-
- Vector<double> last_solution;
-};
-
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &) const {
- const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
- const FullMatrix<double> &values = fe_values.get_shape_values ();
- vector<double> rhs_values (fe_values.n_quadrature_points);
- const vector<double> &weights = fe_values.get_JxW_values ();
-
- vector<Tensor<1,dim> > 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<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- cell_matrix(i,j) += (gradients[i][point] *
- gradients[j][point]) *
- weights[point] /
- sqrt(1+last_solution_grads[i]*last_solution_grads[i]);
- rhs(i) += values(i,point) *
- rhs_values[point] *
- weights[point];
- };
-};
-
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, typename Equation<dim>::ExcPureVirtualFunctionCalled());
-};
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (Vector<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, typename Equation<dim>::ExcPureVirtualFunctionCalled());
-};
-
-
-
-
-template <int dim>
-NonlinearProblem<dim>::NonlinearProblem () :
- tria(0), dof(0) {};
-
-
-
-template <int dim>
-void NonlinearProblem<dim>::run () {
-
- // first reset everything to a virgin state
- clear ();
-
- tria = new Triangulation<dim>();
- dof = new DoFHandler<dim> (*tria);
- set_tria_and_dof (tria, dof);
-
-
- RightHandSide<dim> rhs;
- ZeroFunction<dim> boundary_values;
- StraightBoundary<dim> boundary;
-
- FEQ1<dim> fe;
- PoissonEquation<dim> equation (rhs, last_solution);
- QGauss2<dim> quadrature;
-
- FunctionMap<dim>::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; i<dof->n_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<dim> out;
- ofstream gnuplot(filename.c_str());
- fill_data (out);
- out.build_patches ();
- out.write_gnuplot (gnuplot);
- gnuplot.close ();
- };
-
- last_solution = solution;
- };
-
- Vector<float> error_indicator;
- KellyErrorEstimator<dim> ee;
- QSimpson<dim-1> eq;
- ee.estimate (*dof, eq,
- FunctionMap<dim>::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 ();
-};
+++ /dev/null
-# $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
-
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-#include "poisson.h"
-#include <lac/vector.h>
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-void PoissonEquation<1>::assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<1> &fe_values,
- const DoFHandler<1>::cell_iterator &) const {
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.total_dofs; ++i)
- {
- for (unsigned int j=0; j<fe_values.total_dofs; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point)) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
-
-#endif
-
-
-
-//#if deal_II_dimension >= 2
-
-template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &) const {
- const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
- const FullMatrix<double> &values = fe_values.get_shape_values ();
- vector<double> rhs_values (fe_values.n_quadrature_points);
- const vector<double> &weights = fe_values.get_JxW_values ();
-
- right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values);
-
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
- cell_matrix(i,j) += (gradients[i][point] *
- gradients[j][point]) *
- weights[point];
- rhs(i) += values(i,point) *
- rhs_values[point] *
- weights[point];
- };
-};
-
-//#endif
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (FullMatrix<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, typename Equation<dim>::ExcPureVirtualFunctionCalled());
-};
-
-
-
-template <int dim>
-void PoissonEquation<dim>::assemble (Vector<double> &,
- const FEValues<dim> &,
- const DoFHandler<dim>::cell_iterator &) const {
- Assert (false, typename Equation<dim>::ExcPureVirtualFunctionCalled());
-};
-
-
-
-
-
-
-template class PoissonEquation<2>;
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-
-#include "poisson.h"
-#include <lac/vector.h>
-#include <base/logstream.h>
-
-
-
-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;
-};
-
+++ /dev/null
-/*---------------------------- poisson.h ---------------------------*/
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-#ifndef __poisson_H
-#define __poisson_H
-/*---------------------------- poisson.h ---------------------------*/
-
-
-#include <grid/tria.h>
-#include <dofs/dof_handler.h>
-#include <grid/tria_accessor.h>
-#include <dofs/dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
-#include <dofs/dof_constraints.h>
-#include <numerics/data_out.h>
-#include <base/function.h>
-#include <base/parameter_handler.h>
-#include <fe/fe_lib.lagrange.h>
-#include <base/quadrature_lib.h>
-#include "../problem_base.h"
-#include <lac/sparse_matrix.h>
-
-
-#include <map>
-#include <fstream>
-#include <cmath>
-#include <cstdlib>
-
-
-
-
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation (const Function<dim> &rhs) :
- Equation<dim>(1),
- right_hand_side (rhs) {};
-
- virtual void assemble (FullMatrix<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (FullMatrix<double> &cell_matrix,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- virtual void assemble (Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &cell) const;
- protected:
- const Function<dim> &right_hand_side;
-};
-
-
-
-
-
-
-template <int dim>
-class PoissonProblem : public ProblemBase<dim>,
- 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<dim> *tria;
- DoFHandler<dim> *dof;
-
- Function<dim> *rhs;
- Function<dim> *boundary_values;
-
- Boundary<dim> *boundary;
-};
-
-
-
-
-
-/*---------------------------- poisson.h ---------------------------*/
-/* end of #ifndef __poisson_H */
-#endif
-/*---------------------------- poisson.h ---------------------------*/
+++ /dev/null
-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
+++ /dev/null
-/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1998 */
-
-
-
-#include "poisson.h"
-#include <lac/vector.h>
-#include <grid/grid_generator.h>
-#include <grid/tria_boundary_lib.h>
-#include <numerics/data_out.h>
-
-
-template <int dim>
-class BoundaryValuesSine : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
-
- double x = 1;
-
- for (unsigned int i=0; i<dim; ++i)
- x *= cos(2*3.1415926536*p(i));
- return x;
- };
-
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual void value (const Point<dim> &p,
- Vector<double> &values) const {
- Assert (values.size()==1, ExcDimensionMismatch (values.size(), 1));
-
- double x = 1;
-
- for (unsigned int i=0; i<dim; ++i)
- x *= cos(2*3.1415926536*p(i));
-
- values(0) = x;
- };
-
-
- /**
- * Set #values# to the point values
- * of the function at the #points#.
- * It is assumed that #values# be
- * empty.
- */
- virtual void value_list (const vector<Point<dim> > &points,
- vector<double> &values,
- const unsigned int component) const {
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
- for (unsigned int i=0; i<points.size(); ++i)
- values[i] = BoundaryValuesSine<dim>::value (points[i], component);
- };
-};
-
-
-
-template <int dim>
-class BoundaryValuesJump : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &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 <int dim>
-class RHSTrigPoly : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &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 <int dim>
-class RHSPoly : public Function<dim> {
- public:
- /**
- * Return the value of the function
- * at the given point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int) const;
-};
-
-
-
-
-
-
-
-
-
-template <int dim>
-class CurvedLine :
- public StraightBoundary<dim> {
- public:
- virtual Point<dim>
- get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
-
- virtual Point<dim>
- get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
-};
-
-
-
-template <int dim>
-Point<dim>
-CurvedLine<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-{
- Point<dim> middle = StraightBoundary<dim>::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<x)
- if (y<1-x)
- middle(1) = 0.04*sin(6*3.141592*middle(0));
- else
- middle(0) = 1+0.04*sin(6*3.141592*middle(1));
-
- else
- if (y<1-x)
- middle(0) = 0.04*sin(6*3.141592*middle(1));
- else
- middle(1) = 1+0.04*sin(6*3.141592*middle(0));
-
- return middle;
-};
-
-
-
-template <int dim>
-Point<dim>
-CurvedLine<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-{
- Point<dim> middle = StraightBoundary<dim>::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<x)
- if (y<1-x)
- middle(1) = 0.04*sin(6*3.141592*middle(0));
- else
- middle(0) = 1+0.04*sin(6*3.141592*middle(1));
-
- else
- if (y<1-x)
- middle(0) = 0.04*sin(6*3.141592*middle(1));
- else
- middle(1) = 1+0.04*sin(6*3.141592*middle(0));
-
- return middle;
-};
-
-
-
-
-template <int dim>
-double RHSTrigPoly<dim>::value (const Point<dim> &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 <int dim>
-double RHSPoly<dim>::value (const Point<dim> &p,
- const unsigned int component) const {
- Assert (component==0, ExcIndexRange (component, 0, 1));
-
- double ret_val = 0;
- for (unsigned int i=0; i<dim; ++i)
- ret_val += 2*p(i)*(1.-p(i));
- return ret_val;
-};
-
-
-
-
-
-template <int dim>
-PoissonProblem<dim>::PoissonProblem () :
- tria(0), dof(0), rhs(0), boundary_values(0), boundary(0) {};
-
-
-
-template <int dim>
-PoissonProblem<dim>::~PoissonProblem ()
-{
- clear ();
-};
-
-
-
-template <int dim>
-void PoissonProblem<dim>::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<dim>::clear ();
-};
-
-
-
-
-template <int dim>
-void PoissonProblem<dim>::create_new (const unsigned int) {
- clear ();
-
- tria = new Triangulation<dim>();
- dof = new DoFHandler<dim> (*tria);
- set_tria_and_dof (tria, dof);
-};
-
-
-
-
-template <int dim>
-void PoissonProblem<dim>::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 <int dim>
-bool PoissonProblem<dim>::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<dim>();
- tria->set_boundary (0, *boundary);
- make_zoom_in_grid ();
- break;
- case 2:
- // make ball grid around origin with
- // unit radius
- {
- static const Point<dim> origin;
- boundary = new HyperBallBoundary<dim>(origin, 1.);
- GridGenerator::hyper_ball (*tria, origin, 1.);
- tria->set_boundary (0, *boundary);
- break;
- };
- case 3:
- // set the boundary function
- {
- boundary = new CurvedLine<dim>();
- GridGenerator::hyper_cube (*tria);
- tria->set_boundary (0, *boundary);
- break;
- };
- case 4:
- boundary = new StraightBoundary<dim>();
- tria->set_boundary (0, *boundary);
- make_random_grid ();
- break;
- case 5:
- boundary = new StraightBoundary<dim>();
- tria->set_boundary (0, *boundary);
- GridGenerator::hyper_cube (*tria);
- break;
- case 6:
- boundary = new StraightBoundary<dim>();
- 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<dim>();
- tria->set_boundary (0, *boundary);
- GridGenerator::hyper_L (*tria);
- break;
- case 8:
- boundary = new StraightBoundary<dim>();
- 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 <int dim>
-void PoissonProblem<dim>::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<dim>::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 <int dim>
-void PoissonProblem<dim>::make_random_grid () {
- GridGenerator::hyper_cube (*tria);
- tria->refine_global (1);
-
- Triangulation<dim>::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 <int dim>
-bool PoissonProblem<dim>::set_right_hand_side (ParameterHandler &prm) {
- string rhs_name = prm.get ("Right hand side");
-
- if (rhs_name == "zero")
- rhs = new ZeroFunction<dim>();
- else
- if (rhs_name == "constant")
- rhs = new ConstantFunction<dim>(1.);
- else
- if (rhs_name == "trigpoly")
- rhs = new RHSTrigPoly<dim>();
- else
- if (rhs_name == "poly")
- rhs = new RHSPoly<dim> ();
- else
- return false;
-
- if (rhs != 0)
- return true;
- else
- return false;
-};
-
-
-
-template <int dim>
-bool PoissonProblem<dim>::set_boundary_values (ParameterHandler &prm) {
- string bv_name = prm.get ("Boundary values");
-
- if (bv_name == "zero")
- boundary_values = new ZeroFunction<dim> ();
- else
- if (bv_name == "sine")
- boundary_values = new BoundaryValuesSine<dim> ();
- else
- if (bv_name == "jump")
- boundary_values = new BoundaryValuesJump<dim> ();
- else
- {
- cout << "Unknown boundary value function " << bv_name << endl;
- return false;
- };
-
- if (boundary_values != 0)
- return true;
- else
- return false;
-};
-
-
-
-
-template <int dim>
-void PoissonProblem<dim>::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<dim> fe;
- PoissonEquation<dim> equation (*rhs);
- QGauss2<dim> quadrature;
-
- cout << " Distributing dofs... ";
- dof->distribute_dofs (fe);
- cout << dof->n_dofs() << " degrees of freedom." << endl;
-
- cout << " Assembling matrices..." << endl;
- ProblemBase<dim>::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<dim> 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>;
+++ /dev/null
-/*---------------------------- 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 <base/exceptions.h>
-#include <lac/full_matrix.h>
-#include <lac/vector.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <fe/fe_values.h>
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <vector>
-
-
-/**
- * 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 <int dim>
-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<double> &cell_matrix,
- Vector<double> &rhs,
- const FEValues<dim> &fe_values,
- const typename DoFHandler<dim>::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<double> &cell_matrix,
- const FEValues<dim> &fe_values,
- const typename DoFHandler<dim>::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<double> &rhs,
- const FEValues<dim> &fe_values,
- const typename DoFHandler<dim>::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 <int dim>
-class Assembler : public DoFCellAccessor<dim>
-{
- 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<dim> &dof,
- const bool assemble_matrix,
- const bool assemble_rhs,
- SparseMatrix<double> &matrix,
- Vector<double> &rhs_vector,
- const Quadrature<dim> &quadrature,
- const UpdateFlags &update_flags);
-
- /**
- * Pointer to the dof handler object
- * to be used to iterate on.
- */
- const DoFHandler<dim> &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<double> &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<double> &rhs_vector;
-
- /**
- * Pointer to a quadrature object to be
- * used for this assemblage process.
- */
- const Quadrature<dim> &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<dim> *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<dim> &);
-
- /**
- * Exception.
- */
- DeclException0 (ExcNoAssemblingRequired);
- /**
- * Exception.
- */
- DeclException0 (ExcInvalidData);
- /**
- * Exception.
- */
- /**
- * Exception.
- */
- private:
- /**
- * Store a local cell matrix.
- */
- FullMatrix<double> cell_matrix;
-
- /**
- * Right hand side local to cell.
- */
- Vector<double> 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<double> &matrix;
-
- /**
- * Pointer to the vector to be assembled
- * by this object.
- */
- Vector<double> &rhs_vector;
-
- /**
- * The finite element evaluated at the
- * quadrature points.
- */
- FEValues<dim> fe_values;
-};
-
-
-#include <lac/sparse_matrix.h>
-#include <base/exceptions.h>
-#include <dofs/dof_constraints.h>
-#include <fe/fe_update_flags.h>
-#include <map>
-#include <string>
-
-template <int dim> class DataOut;
-template <int dim> class Assembler;
-template <int dim> 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 <int dim>
-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<unsigned char,const Function<dim>*> FunctionMap;
- /**
- * Typdedef an iterator which assembles
- * matrices and vectors.
- */
- typedef TriaActiveIterator<dim, Assembler<dim> > 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<dim> *tria,
- DoFHandler<dim> *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<dim> &equation,
- const Quadrature<dim> &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<dim> &) 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<dim> *tria;
-
- /**
- * Pointer to the degree of freedom handler
- * to work with.
- */
- DoFHandler<dim> *dof_handler;
-
- /**
- * Sparsity pattern of the system matrix.
- */
- SparsityPattern system_sparsity;
-
- /**
- * System matrix.
- */
- SparseMatrix<double> system_matrix;
-
- /**
- * Vector storing the right hand side.
- */
- Vector<double> right_hand_side;
-
- /**
- * Solution vector.
- */
- Vector<double> solution;
-
- /**
- * List of constraints introduced by
- * hanging nodes.
- */
- ConstraintMatrix constraints;
-
- friend class Assembler<dim>;
-};
-
-
-#include <grid/tria_iterator.h>
-#include <grid/tria_iterator.templates.h>
-#include <fe/fe.h>
-#include <lac/full_matrix.h>
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-#include <base/quadrature.h>
-#include <fe/mapping_q1.h>
-
-
-// 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 <int dim>
-Assembler<dim>::AssemblerData::AssemblerData (const DoFHandler<dim> &dof,
- const bool assemble_matrix,
- const bool assemble_rhs,
- SparseMatrix<double> &matrix,
- Vector<double> &rhs_vector,
- const Quadrature<dim> &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 <int dim>
-Assembler<dim>::Assembler (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AssemblerData *local_data) :
- DoFCellAccessor<dim> (tria,level,index, &local_data->dof),
- cell_matrix (dof_handler->get_fe().dofs_per_cell),
- cell_vector (Vector<double>(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 <int dim>
-void Assembler<dim>::assemble (const Equation<dim> &equation) {
- // re-init fe values for this cell
- fe_values.reinit (DoFHandler<dim>::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<dim>::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<unsigned int> 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<n_dofs; ++i)
- for (unsigned int j=0; j<n_dofs; ++j)
- matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
-
- // distribute cell vector
- if (assemble_rhs)
- for (unsigned int j=0; j<n_dofs; ++j)
- rhs_vector(dofs[j]) += cell_vector(j);
-};
-
-
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-
-template <int dim>
-Equation<dim>::Equation (const unsigned int n_equations) :
- n_eq(n_equations) {};
-
-
-template <int dim>
-void Equation<dim>::assemble (FullMatrix<double> &,
- Vector<double> &,
- const FEValues<dim> &,
- const typename DoFHandler<dim>::cell_iterator &) const
-{
- Assert (false, ExcPureVirtualFunctionCalled());
-};
-
-
-template <int dim>
-void Equation<dim>::assemble (FullMatrix<double> &,
- const FEValues<dim> &,
- const typename DoFHandler<dim>::cell_iterator &) const
-{
- Assert (false, ExcPureVirtualFunctionCalled());
-};
-
-
-template <int dim>
-void Equation<dim>::assemble (Vector<double> &,
- const FEValues<dim> &,
- const typename DoFHandler<dim>::cell_iterator &) const
-{
- Assert (false, ExcPureVirtualFunctionCalled());
-};
-
-#include <dofs/dof_tools.h>
-#include <base/function.h>
-#include <fe/fe.h>
-#include <base/quadrature.h>
-#include <lac/vector.h>
-#include <lac/precondition.h>
-#include <lac/solver_cg.h>
-#include <lac/vector_memory.h>
-
-#include <map>
-#include <numeric>
-#include <algorithm>
-#include <cmath>
-
-
-template <int dim>
-ProblemBase<dim>::ProblemBase () :
- tria(0),
- dof_handler(0),
- system_sparsity(), // dummy initialisation, is later reinit'd
- system_matrix() // dummy initialisation, is later reinit'd
-{};
-
-
-template <int dim>
-void ProblemBase<dim>::set_tria_and_dof (Triangulation<dim> *t,
- DoFHandler<dim> *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 <int dim>
-void ProblemBase<dim>::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 <int dim>
-ProblemBase<dim>::~ProblemBase () {};
-
-
-template <int dim>
-void ProblemBase<dim>::assemble (const Equation<dim> &equation,
- const Quadrature<dim> &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<dim>::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<unsigned int, double> 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<dim>::apply_boundary_values (boundary_value_list,
- system_matrix, solution,
- right_hand_side,
- true);
-};
-
-
-template <int dim>
-void ProblemBase<dim>::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 <int dim>
-void ProblemBase<dim>::fill_data (DataOut<dim> &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 <int dim>
-std::string ProblemBase<dim>::get_solution_name () const
-{
- return "solution";
-};
-
-
-#endif
-