--- /dev/null
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html>
+<head>
+ <title>DEAL tutorial: the Laplace problem</title>
+ <meta name="keyword" content="DEAL,DEAL tutorial">
+ <link href="../dealtut.css" rel="StyleSheet" media="screen" type="text/css" title="DEAL tutorial">
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="DEAL,DEAL tutorial">
+</head>
+<body>
+<h1>Assembling the problem</h1>
+
+<h2>What's to be done</h2>
+In order to assemble the matrices we basically need to:
+<ol>
+<li><a href=#matrix">Generate the matrices</a>, i.e. call the <em>DEAL</em> functions that reserve
+storage space for us.
+</li>
+<li>
+<a href="#calcfe">Calculate the finite element trial functions</a>
+</li>
+<li>
+Traverse all existing cells and <a href="#integrate">integrate the problem</a>
+using the discretized laplace operator
+</li>
+<li>
+Traverse all the cell faces and <a href="#boundary">set the appropriate boundary conditions</a>
+</li>
+<li>
+Insert the local matrices we have used into the global matrix using the
+appropriate <em>DEAL</em> functions
+</li>
+</ol>
+
+<h2>...and how to do it</h2>
+
+<h3>Function parameters</h3>
+<pre>
+<code>
+void
+Laplace::assemble_primal(const Function<2>& exact, const Function<2>&)
+{
+</code>
+</pre>
+
+<h3><a name="matrix">Generating the matrix structures<a></h3>
+
+<p>
+First we generate an n times n square matrix where n is the number
+of the degrees of freedom, i.e. the number of points of our discretization.
+The parameter <em>max_couplings_between_dofs()</em> returns the maximum
+number of couplings between degrees of freedom and allows <em>DEAL</em>
+to generate the matrix structure more efficiently, for most of its
+elements are zero.
+</p>
+<p>
+Afterwards the <em>hanging nodes</em> are copied into the matrix, i.e.
+the matrix is generated.
+</p>
+<pre>
+<code>
+ matrix_structure.reinit(dof_primal.n_dofs(),dof_primal.n_dofs(),
+ dof_primal.max_couplings_between_dofs());
+ dof_primal.make_sparsity_pattern(matrix_structure);
+ hanging_nodes.clear();
+ dof_primal.make_constraint_matrix(hanging_nodes);
+ hanging_nodes.condense(matrix_structure);
+</code>
+</pre>
+<p>
+The problem is of the form <tt>Au=f</tt>:
+</p>
+<pre>
+<code>
+ A.reinit(matrix_structure);
+ f.reinit(dof_primal.n_dofs());
+</code>
+</pre>
+
+<h3><a name+"calcfe">Calculatinginite element trial functions</a></h3>
+<p>
+The two lines below calculate trial functions for the finite elements and
+for their faces using Gaussian quadrature.
+</p>
+<pre>
+<code>
+ FEValues<2> fevalues(fe_primal, qc_primal, UpdateFlags(update_gradients |
+ update_JxW_values));
+ FEFaceValues<2> ffvalues(fe_primal, qf_primal,
+ UpdateFlags(update_JxW_values | update_q_points));
+</code>
+</pre>
+
+<h3><a name="integrate">Integrating the problem</a></h3>
+<p>
+Integration is done locally. Therefore we need appropriate definitions for
+<ul>
+<li>
+an index vector that will allow us to reassemble the global matrix later on
+</li>
+<li>a vector of doubles with the dimension of the total number of degrees of freedom</li>
+<li>and a square matrix of doubles with the same dimension</li>
+</ul>
+</p>
+<pre>
+<code>
+ vector<int> indices(fe_primal.total_dofs);
+ dVector elvec(fe_primal.total_dofs);
+
+ dFMatrix elmat(fe_primal.total_dofs);
+</code>
+</pre>
+<p>
+Next we traverse all the cells and integrate the Laplace problem using the
+discretized Laplace operator. <tt>qc_primal</tt> is a Gaussian quadrature.
+</p>
+<pre>
+<code>
+ for (DoFHandler<2>::active_cell_iterator c = dof_primal.begin_active()
+ ; c != dof_primal.end() ; ++c)
+ {
+ fevalues.reinit(c, stb);
+ elmat.clear();
+ elvec.clear();
+ c->get_dof_indices(indices);
+
+ for (unsigned k=0;k<qc_primal.n_quadrature_points;++k)
+ {
+ for (unsigned i=0;i<fe_primal.total_dofs;++i)
+ {
+ const Point<2> dv = fevalues.shape_grad(i,k);
+
+ for (unsigned j=0;j<fe_primal.total_dofs;++j)
+ {
+ const Point<2> du = fevalues.shape_grad(j,k);
+
+ elmat(i,j) += fevalues.JxW(k)
+ * du * dv
+ ;
+
+ }
+ }
+ }
+</code>
+</pre>
+
+<h3><a name="boundary">Setting boundary conditions</a></h3>
+
+<p>
+There are two <em>DEAL</em> functions relevant for us at the moment:
+<pre>
+<code>
+static_void interpolate_boundary_values(...)
+</code>
+</pre>
+which does exactly what it says. This function returns a list of pairs
+of boundary indicators and the according functions denoting the respective
+Dirichlet boundary values.
+</p>
+<p>
+This output is used by
+<pre>
+<code>
+static void apply_boundary_values(...)
+</code>
+</pre>
+that inserts the proper boundary conditions into the equation system:
+</p>
+<P>
+<pre><code>
+ map<int,double> boundary_values;
+ DoFHandler<2>::FunctionMap dirichlet_bc;
+ BoundaryFct bfkt;
+ dirichlet_bc[0]=&bfkt;
+ VectorTools<2>::interpolate_boundary_values(dof_primal,dirichlet_bc,fe_primal,boundary,boundary_values);
+ u.reinit(f);
+ MatrixTools<2>::apply_boundary_values(boundary_values,A,u,f);
+</code></pre>
+First, we need a few definitions:
+<ul>
+<li>
+<code>boundary_values</code> maps boundary values computed by <code>interpolate_boundary_values</code> to boundary indicators,i.e. to boundaries.
+</li>
+<li><code>dirichlet_bc</code> maps boundary functions, supplied by us, to boundary indicators. The boundary functions compute the boundary values.
+</li>
+<li><code>bfkt</bfkt> is a function returning <code>sin(x)*sin(y)
+</code>, thereby supplying boundary values.
+</ul>
+This may seem a bit confusing. What actually happens is the following:
+<ol>
+<li><code>interpolate_boundary_values</code> takes the boundary functions
+<code>bfkt</code>, its relation to boundaries <code>dirichlet_bc</code> and
+the triangulation <code>dof_primal, fe_primal</code> and returns a
+mapping <code>boundary_values</code> that maps values instead of functions
+to our boundaries. The function looks at <em>all</em> the boundaries. All we
+ever need to do is specify the initial triangulation.
+</li>
+<li><code>apply_boundary_values</code> subsequently takes that mapping and
+our equation system <code>Au=f</code> and inserts the boundary values into
+the equation system which can then be solved.
+</li>
+</ol>
+</p>
+
+<hr>
+<p>
+<a href="../index.html">Back to the tutorial index</a>
+</p>
+
+<hr>
+<address><a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
+<p>
+Last modified: Fri Feb 12, 1999
+</p>
+</body>
+</html>
--- /dev/null
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+root = ../../../../
+
+vpath %.a $(root)/deal.II/lib
+vpath %.a $(root)/lac/lib
+vpath %.a $(root)/base/lib
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target = laplace
+
+# 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 = 60 4
+
+# 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
+
+deal_II_dimension=2
+
+
+###############################################################################
+# Internals
+
+INCLUDE=-I$(root)/deal.II/include -I$(root)/lac/include -I$(root)/base/include
+
+CXXFLAGS.g= -DDEBUG -g -Wall -W -pedantic -Wconversion \
+ -Winline -Woverloaded-virtual \
+ $(INCLUDE) -Ddeal_II_dimension=$(deal_II_dimension)
+CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \
+ -felide-constructors -fnonnull-objects \
+ $(INCLUDE) \
+ -Ddeal_II_dimension=$(deal_II_dimension)
+
+ifeq ($(shell uname),Linux)
+CXX = g++
+endif
+
+ifeq ($(shell uname),SunOS)
+CXX = g++
+endif
+
+
+%.go : %.cc #Makefile
+ @echo ============================ Compiling with debugging information: $<
+ @echo $(CXX) ... -c $< -o $@
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc #Makefile
+ @echo ============================ Compiling with optimization: $<
+ @echo $(CXX) ... -c $< -o $@
+ @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# get lists of files we need
+cc-files = $(wildcard *.cc)
+o-files = $(cc-files:.cc=.o)
+go-files = $(cc-files:.cc=.go)
+h-files = $(wildcard *.h)
+lib-h-files = $(wildcard $(root)/deal.II/include/*/*.h)
+
+# list of libraries needed to link with
+libs = ./Obj.a $(wildcard $(root)/deal.II/lib/lib*2d.a) $(root)/lac/lib/liblac.a
+libs.g = ./Obj.g.a $(wildcard $(root)/deal.II/lib/lib*2d.g.a) $(root)/lac/lib/liblac.a $(root)/base/lib/libbase.g.a
+
+#$(root)/deal.II/lib/deal_II_2d.g -llac.g -lbase.g#-lbasic.g -lfe.g -lgrid.g -lbasic.g -llac.g
+
+
+# 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
+
+
+
+# make rule for the target
+$(target) : $(libraries) $(user-libs)
+ @echo ============================ Linking $@
+ $(CXX) $(flags) -o $@ $^
+
+# 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) -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 $(root)/deal.II/Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+ > Makefile.dep
+
+
+include Makefile.dep
+
--- /dev/null
+// $Id$
+
+// JS.
+const char* funcversion = "Functions: $Revision$";
+
+#include "functions.h"
+
+#include <cmath>
+
+double
+WeightFunction::operator() (const Point<2>& p) const
+{
+ // double r = p(0)*p(0) + p(1) * p(1);
+ //if (r>=.8) return 0.;
+ //return 1.-r*(2.-r);
+ return 1.;
+
+}
+
+double
+BoundaryFct::operator ()(const Point<2> &p) const
+{
+ return sin(4*M_PI*p(0))*sin(4*M_PI*p(1));
+}
--- /dev/null
+// $Id$
+
+// JS.Wird das File ueberhaupt gebraucht ?
+
+#include <base/function.h>
+
+class WeightFunction
+ : public Function<2>
+{
+public:
+ WeightFunction()
+ {}
+ virtual double operator()(const Point<2> &p) const;
+};
+
+class BoundaryFct
+ : public Function<2>
+{
+ public:
+ virtual double operator()(const Point<2> &p) const;
+};
+
--- /dev/null
+// $Id$
+
+const char* laplaceversion = "Laplace: $Revision$";
+
+#include "laplace.h"
+#include "functions.h"
+
+#include <lac/solver_cg.h>
+#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary.h>
+#include <grid/dof_constraints.h>
+#include <basic/data_io.h>
+#include <base/function.h>
+#include <base/parameter_handler.h>
+#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_lib.criss_cross.h>
+#include <fe/fe_values.h>
+#include <base/quadrature_lib.h>
+#include <numerics/matrices.h>
+#include <numerics/vectors.h>
+
+#include <base/logstream.h>
+
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+#define PRIMEL FELinear<2>
+#define DUEL FEQuadraticSub<2>
+
+// Finite Elements
+
+static PRIMEL fe_primal;
+// JS.static DUEL fe_dual;
+
+// Quadrature formulae
+
+static QGauss2<2> qc_primal;
+static QGauss2<1> qf_primal;
+// static QGauss3<2> qc_dual;
+// static QGauss3<1> qf_dual;
+static QGauss5<2> qc_integrate;
+static QGauss5<1> qf_integrate;
+
+StraightBoundary<2> stb;
+
+// JS.ist im Moment noch PureTransportSolution...
+Laplace::Laplace()
+ : dof_primal(&tr)
+{
+ // JS.Triangulation generieren. Zellränder werden numeriert.
+ tr.create_hypercube(-1.,1.);
+
+ // JS.Freiheitsgrade verteilen. D.h. Zellen numerieren.
+ dof_primal.distribute_dofs(fe_primal);
+ }
+
+Laplace::~Laplace()
+{}
+
+
+// JS.Gitter verfeinern.
+void
+Laplace::remesh(unsigned int steps)
+{
+ if (tr.n_levels() <= 1)
+ {
+ tr.refine_global(1); //JS.Lokal ist execute_coarsening_etc...
+ }
+
+ if (steps)
+ tr.refine_global(steps);
+ else
+ tr.execute_coarsening_and_refinement();
+
+ // JS. Freiheitsgrade neu verteilen.
+ dof_primal.distribute_dofs(fe_primal);
+ // JS. und dem Problem angemessener nochmal numerieren.
+ dof_primal.renumber_dofs(Cuthill_McKee);
+ deallog << "Cells " << tr.n_active_cells()
+ << " PrimalDoFs " << dof_primal.n_dofs()
+ << endl;
+}
+
+// JS.Primales Problem zusammenstellen.
+void
+Laplace::assemble_primal(const Function<2>&,const Function<2>&)
+{
+ deallog << "Assembling primal problem" << endl;
+ // JS. Platz für neue Matrix mit (?) (2x = quadratisch) Anzahl der Zellen,
+ // Anzahl der Kopplungen (Matrix dünn besetzt, für effizientes Speichern)
+ matrix_structure.reinit(dof_primal.n_dofs(),dof_primal.n_dofs(),
+ dof_primal.max_couplings_between_dofs());
+ // JS.Hängende Noden in die Matrix einbauen; d.h. Matrix generieren.
+ dof_primal.make_sparsity_pattern(matrix_structure);
+ hanging_nodes.clear();
+ dof_primal.make_constraint_matrix(hanging_nodes);
+ hanging_nodes.condense(matrix_structure);
+
+ //JS.Problem der Form Au=f.
+ A.reinit(matrix_structure);
+ f.reinit(dof_primal.n_dofs());
+
+ // JS.Ansatzfunktionen auf Zellrändern im Voraus berechnen aus
+ // Effizienzgründen.
+ FEValues<2> fevalues(fe_primal, qc_primal, UpdateFlags(update_gradients |
+ update_JxW_values));
+ FEFaceValues<2> ffvalues(fe_primal, qf_primal,
+ UpdateFlags(update_JxW_values | update_q_points));
+ //JS.Ab hier lokales Problem el... = Finites Element...
+ //JS. Index für eine Zelle, für späteren Einbau in globale Matrix.
+ vector<int> indices(fe_primal.total_dofs);
+ dVector elvec(fe_primal.total_dofs);
+
+ dFMatrix elmat(fe_primal.total_dofs);
+
+ // JS.Einmal alle Zellen durchlaufen
+ for (DoFHandler<2>::active_cell_iterator c = dof_primal.begin_active()
+ ; c != dof_primal.end() ; ++c)
+ {
+ fevalues.reinit(c, stb);
+ elmat.clear();
+ elvec.clear();
+ c->get_dof_indices(indices);
+
+ // JS.Integration des Problems. Diese Schleifenfolge für Effizienz.
+ for (unsigned k=0;k<qc_primal.n_quadrature_points;++k)
+ {
+ for (unsigned i=0;i<fe_primal.total_dofs;++i)
+ {
+ const Point<2> dv = fevalues.shape_grad(i,k);
+
+
+ for (unsigned j=0;j<fe_primal.total_dofs;++j)
+ {
+ const Point<2> du = fevalues.shape_grad(j,k);
+
+ elmat(i,j) += fevalues.JxW(k)
+ * du * dv
+ ;
+
+ }
+ }
+ }
+ // JS.Lokale Matrix in globale einbauen.
+ for (unsigned i=0;i<fe_primal.total_dofs;++i)
+ {
+ // f(indices[i]) += elvec(i);
+ f(indices[i]) = 0;
+
+ for (unsigned j=0;j<fe_primal.total_dofs;++j)
+ {
+ A.add(indices[i], indices[j], elmat(i,j));
+ }
+ }
+ }
+
+ // JS. Randwerte.
+ map<int,double> boundary_values;
+ DoFHandler<2>::FunctionMap dirichlet_bc;
+ BoundaryFct bfkt;
+ dirichlet_bc[0]=&bfkt;
+ VectorTools<2>::interpolate_boundary_values(dof_primal,dirichlet_bc,
+ fe_primal,boundary,
+ boundary_values);
+ u.reinit(f);
+ MatrixTools<2>::apply_boundary_values(boundary_values,A,u,f);
+
+ cout << "u: " << u.l2_norm() << endl;
+ cout << "f: " << f.l2_norm() << endl;
+
+ // JS.Hängende Noden einbauen.
+ // hanging_nodes.condense(A);
+ //hanging_nodes.condense(f);
+}
+
+// JS. Primales Problem lösen.
+void
+Laplace::solve_primal()
+{
+ deallog.push("Solve");
+
+ // JS.Empfindlichkeit des Lösers einstellen.
+ SolverControl control(1000, 1.e-10);
+ // JS. Löser definieren. modifiziertes cg-Verfahren,
+ SolverCG<AdvMatrix, dVector> solver(control, mem);
+
+ // JS.???
+ cout << "f:L2-norm=" << f.l2_norm() << endl;
+ // u.reinit(f);
+
+ // JS.lösen.
+ solver.solve(A,u,f);
+ cout << "u:L2-norm=" << u.l2_norm() << endl;
+
+ // JS.???
+ hanging_nodes.distribute(u);
+
+ deallog.pop();
+}
+
+
+// JS. Datenausgabe im Gnuplot-Format.
+void Laplace::write_data(const char* name)
+{
+ deallog << "Writing gnuplot" << endl;
+
+ DataOut<2> out;
+ char fname[100];
+ sprintf(fname,"P_%s",name);
+
+ {
+ ofstream gnuplot(fname);
+
+ out.clear_data_vectors();
+ out.attach_dof_handler(dof_primal);
+ out.add_data_vector(u,"solution","kg");
+
+ out.write_gnuplot (gnuplot, 1);
+ gnuplot.close ();
+ }
+
+}
+
+
+// JS. Ergebnis zurückgeben. Wie funktioniert das ?
+double
+Laplace::result(const Function<2>& interior, const Function<2>& boundary)
+{
+ double erg = 0., ergex = 0.;
+ FEValues<2> fevalues(fe_primal, qc_integrate,
+ UpdateFlags(update_q_points | update_JxW_values));
+ FEFaceValues<2> ffvalues(fe_primal, qf_integrate,
+ UpdateFlags(update_q_points | update_JxW_values));
+ vector<double> uh(qc_integrate.n_quadrature_points);
+ vector<double> uf(qf_integrate.n_quadrature_points);
+
+ // JS.Alle Zellen durch.
+ for (DoFHandler<2>::active_cell_iterator c = dof_primal.begin_active()
+ ; c != dof_primal.end() ; ++c)
+ {
+ double s = 0.;
+ fevalues.reinit(c, stb);
+
+ fevalues.get_function_values(u, uh);
+
+ for (unsigned int k=0;k<qc_integrate.n_quadrature_points;++k)
+ {
+ s += fevalues.JxW(k)
+ * uh[k] * interior(fevalues.get_quadrature_points()[k]);
+ ergex += fevalues.JxW(k)
+ // * exact(fevalues.get_quadrature_points()[k])
+ * interior(fevalues.get_quadrature_points()[k]);
+ }
+ // JS.Alle Zellränder.
+ for (unsigned fi=0;fi<GeometryInfo<2>::faces_per_cell;++fi)
+ {
+ DoFHandler<2>::face_iterator f = c->face(fi);
+ unsigned char bi = f->boundary_indicator();
+ if (bi == 0xFF) continue;
+ ffvalues.reinit(c, fi, stb);
+ ffvalues.get_function_values(u, uf);
+
+ // JS.??? Integrationspunkte ???
+ for (unsigned k=0;k<qf_primal.n_quadrature_points;++k)
+ {
+ s += ffvalues.JxW(k)
+ * uf[k] * boundary(ffvalues.get_quadrature_points()[k]);
+ }
+ }
+ erg += s;
+ }
+ deallog << "Results " << setw(18) << setprecision(15) << erg << " " << ergex << endl;
+ return erg;
+}
+
+// JS. Gitter anpassen
+void
+Laplace::adapt()
+{
+ tr.refine_and_coarsen_fixed_fraction(f, .5, 0);
+}
+
+
+
--- /dev/null
+// $Id$
+
+#include <base/exceptions.h>
+#include <base/point.h>
+#include <grid/tria.h>
+#include <grid/dof.h>
+#include <grid/dof_constraints.h>
+#include <grid/tria_boundary.h>
+#include <lac/vector_memory.h>
+#include <lac/dvector.h>
+#include <lac/dsmatrix.h>
+
+class AdvMatrix :
+ public dSMatrix
+{
+public:
+ void precondition(dVector& dst, const dVector& src) const
+ {
+ dSMatrix::precondition_SSOR(dst, src);
+ }
+};
+
+
+class Laplace
+{
+ Function<2>& exact;
+protected:
+ Point<2> direction;
+ Triangulation<2> tr;
+ DoFHandler<2> dof_primal;
+// DoFHandler<2> dof_dual;
+
+ dSMatrixStruct matrix_structure;
+ AdvMatrix A;
+
+ dVector u;
+ dVector z;
+ dVector f;
+
+ PrimitiveVectorMemory<dVector> mem;
+
+ ConstraintMatrix hanging_nodes;
+
+ StraightBoundary<2> boundary;
+
+public:
+ Laplace();
+ ~Laplace();
+
+ void remesh(unsigned int global_refine = 0);
+ void assemble_primal(const Function<2>& boundary, const Function<2>& rhs);
+ void solve_primal();
+
+ double result(const Function<2>& interior, const Function<2>& boundary);
+
+ void adapt();
+
+ void write_data(const char* name);
+
+ void fill_vector(dVector& v, const Function<2>& f) const;
+};
+
--- /dev/null
+// $Id$
+
+#include "laplace.h"
+#include "functions.h"
+
+#include <iostream.h>
+#include <fstream.h>
+#include <stdlib.h>
+
+#include <base/logstream.h>
+#include <base/jobidentifier.h>
+
+ZeroFunction<2> zero;
+
+char fname[30];
+
+main(int argc, char** argv)
+{
+ // JS.Logfile erzeugen, als Stream konzipiert
+
+ ofstream logfile("T");
+ deallog.attach(logfile);
+
+ if (argc==1)
+ cerr << "Usage: " << argv[0] << "firstgrid\nUsing 3"
+ << endl;
+
+ int firstgrid = 3;
+
+ if (argc>=2) firstgrid = atoi(argv[1]);
+
+ deallog << "Firstgrid " << firstgrid << endl;
+
+ // JS.Benötigte Funktionen zur Lösung des Problems
+ WeightFunction weight;
+ BoundaryFct boundary;
+ Laplace lap;
+
+ // JS.Logstream ist ein stack, ab hier wird "Adaptive:" vor jede Zeile
+ // gestellt
+ deallog.push("Adaptive");
+ deallog.depth_console(2);
+
+ for (unsigned step = 0; step < 3 ; ++step)
+ {
+ deallog << "Step " << step << endl;
+ // JS.Beim ersten Mal ein verfeinertes Grid erzeugen, firstgrid=Anzahl
+ // der Verfeinerungen
+ if (!step)
+ lap.remesh(firstgrid);
+ else
+ {
+//JS. lap.adapt();
+ lap.remesh(1);
+ }
+
+ deallog.push("Primal");
+
+ // JS.exakte Lösung mit 0 auf der rechten Seite
+ lap.assemble_primal(boundary,zero);
+ lap.solve_primal();
+
+ // JS.ab hier kein "Adaptive:" mehr
+ deallog.pop();
+
+ sprintf(fname,"T%02d",step);
+ // JS.Daten zurückschreiben. Aber welche ?
+ lap.result(weight,boundary);
+ lap.write_data(fname);
+ }
+ // JS.Logfile sauber abschließen und Schluß
+ deallog.pop();
+ deallog.detach();
+}
<head>
<title>The Laplace Problem</title>
<link href="../dealtut.css" rel="StyleSheet" title="DEAL Tutorial">
- <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de
+i <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de
>">
<meta name="keywords" content="DEAL,DEAL tutorial">
</head>
where a triangulation is generated and degrees of freedom are discussed</p>
</li>
<li>
-<strong><a href="assemble.html">Assembling the problem</a></strong>
+<strong><a href="assemble.html">Assembling the problem matrix</a></strong>
<p>
-where the matrices describing the problem and the boundary conditions are
-built</p>
+where the matrices describing the problem is assembled
+and the boundary conditions are set
+</p>
</li>
<li>
<strong><a href="solution.html">Solving the problem</a></strong>
<hr>
<p>
-<a href="../../index.html">Back to the tutorial index</a></p>
+<a href="../index.html">Back to the tutorial index</a>
+</p>
<hr>
<address>
<a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
<!-- Created: Tue Jan 5 12:50:29 MET 1999 -->
-<!-- hhmts start -->
<p>
-Last modified: Tue 5 Jan 1999 <!-- hhmts end-->
+Last modified: Mon 15 Feb 1999
</p>
</body>
</html>
--- /dev/null
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html>
+<head>
+ <title>DEAL tutorial: the Laplace problem</title>
+ <meta name="keyword" content="DEAL,DEAL tutorial">
+ <link href="../dealtut.css" rel="StyleSheet" media="screen" type="text/css" title="DEAL tutorial">
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="DEAL,DEAL tutorial">
+</head>
+<body>
+<h1>The class Laplace</h1>
+<p>
+The class <em>Laplace</em> contains most of the code we need to actually solve the problem at
+hand. The purpose of this chapter is to establish the elements that are needed for this class
+or any one like it.
+</p>
+
+<p>
+Let's have a look at the class definition:
+<pre>
+<code>
+class Laplace
+{
+ Function<2>& exact;
+protected:
+ Triangulation<2> tr;
+ DoFHandler<2> dof_primal;
+
+ dSMatrixStruct matrix_structure;
+ LapMatrix A;
+
+ ConstraintMatrix hanging_nodes;
+
+<code>
+</pre>
+<p>
+These few lines define several important elements: The right hand side of the equation
+<em>exact</em>, the triangulation <em>tr</em>, i.e. the grid, and a handler for the degrees
+of freedom for the finite elements <em>dof_primal</em>, all for the two-dimensional case.
+In addition three matrices are defined (the matrix <em>A</em> defining our problem). Note that
+in order to solve any problem at all with DEAL the definitions above are paramount.
+</p>
+
+<p>
+The next bit are the functions actually called by our program, the working part, so to say.
+The constructor has the task of generating a triangulation, too.
+</p>
+<pre>
+<code>
+public:
+ Laplace(Function<2>& solution);
+ ~Laplace();
+</code>
+</pre>
+<p>
+The next few functions refine the grid - non-adaptively - assemble the primal problem and
+call the appropriate solver.
+</p>
+<pre>
+<code>
+ void remesh(unsigned int global_refine = 0);
+ void assemble_primal(const Function<2>& boundary, const Function<2>& rhs);
+ void solve_primal();
+</code>
+</pre>
+<!--
+ double result(const Function<2>& interior, const Function<2>& boundary);
+ double estimate(const Function<2>& boundary, const Function<2>& rhs);
+
+ void adapt();
+
+ void write_data(const char* name);
+
+ void fill_vector(dVector& v, const Function<2>& f) const;
+-->
+};
+</code>
+</pre>
+
+<hr>
+<p>
+<a href="../index.html">Back to the tutorial index</a>
+</p>
+
+<hr>
+<p>
+<address><a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
+</p>
+
+<p>
+Last modified: Fri Feb 12, 1999
+</p>
+</body>
+</html>
-<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
- "http://www.w3.org/TR/REC-html40/strict.dtd">
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
<html>
<head>
<title>DEAL tutorial: the Laplace problem</title>
<li> <a href="#program">Inclusion of the necessary headers and variable
definitions.</a>
</li>
- <li>
- <a href="#Data">Creation of a logfile and data logging.</a>
- </li>
<li>
<a href="#Instantiation">Instantiation of the C++ class <em>Laplace</em>
dealing with the problem</a>.
</li>
- <li>Calling the <em>Laplace</em> methods for assemblage and solution of the
- problem.
+ <li><a href="#laplace">Calling the <em>Laplace</em> methods for assemblage and solution of the
+ problem.</a>
</li>
<li>
- Writing out of the output data and termination.
+ <a href="#output">Writing out of the output data and termination.</a>
</li>
</ol>
</code>
</pre>
-<p>
-Next, we need some standard include files
-</p>
-
-<pre>
-<code>
-#include <iostream.h&rt;
-#include <fstream.h&rt;
-#include <stdlib.h&rt;
-</code>
-</pre>
-
-<p>
-and some DEAL specific include files that enable us to do runtime logging and
-nicely identify our program:</p>
-<pre>
-<code>
-#include <base/logstream.h&rt;
-#include <base/jobidentifier.h&rt;
-</code>
-</pre>
-
-<p>
-and some more stuff which will be discussed when appropriate:
-</p>
-
-<pre>
-<code>
-ZeroFunction<2> zero;
-char fname[30];
-extern const char* funcversion;
-extern const char* laplaceversion;
-char versioninfo[500];
-
-const char* JobIdentifier::program_id()
-{
- sprintf(versioninfo,"%s\n%s", funcversion, laplaceversion);
- return versioninfo;
-}
-</code>
-</pre>
<p>
Our main program starts here. It will take the initial number of global grid
refinements as its argument. That means that the number 3 as argument will
{
</code>
</pre>
-
-<h4><a name="Data">Data logging</a></h4>
-<p>
-We want to log info about the program run to a file. Logfiles are streams,
-which behave like stacks in several ways. How to make use of the stack concept
-will be apparent later.
-</p>
-<pre>
-<code>
-ofstream logfile("T");
-deallog.attach(logfile);
-</code>
-</pre>
-
<h4><a name="Instantiation">Instantiation of needed classes</a></h4>
<p>
We need a solution function and an instance of the class dealing with the
<p>
In addition we need to do some refinement (the command line argument was
-previously stored in the variable <var>firstgrid</var>.Please note that every
-log entry will now be preceeded by "Adaptive: " and that the stack height is
-increased by one. Only entries with a stack height of two or less will we
-written to the console.
+previously stored in the variable <var>firstgrid</var>.
</p>
<pre><code>
-deallog.push("Adaptive");
-deallog.depth_console(2);
-
for (unsigned step = 0; step < 3 ; ++step)
{
- deallog << "Step " << step << endl;
+
if (!step)
lap.remesh(firstgrid);
else
}
</code></pre>
+<h4><a name="laplace">Problem assemblage and solution</a></h4>
<p>
-Now we have our log entries preceded by "Primal" since we solve the primal problem,
-our class assembles and solves the problem; the solution is exact (as defined above)
+Our class assembles and solves the primal problem; the solution is exact (as defined above)
and the right hand side of the equation is zero (as defined above). If the right
-hand side were not zero we would solve the Poisson equation instead. Afterwards
-"Primal" is removed from all the following log entries.
+hand side were not zero we would solve the Poisson equation instead.
</p>
<pre><code>
- deallog.push("Primal");
lap.assemble_primal(exact, zero);
lap.solve_primal();
- deallog.pop();
</code></pre>
+
+<h4><a name="output">Data output</a></h4>
<p>
-Finally the solution is written to a file and the logfile is closed.
-</p>
+Finally the solution is written to a file.</p>
<pre>
<code>
sprintf(fname,"T%02d",step);
lap.write_data(fname);
}
- deallog.pop();
- deallog.detach();
}
</code></pre>
<hr>
<p>
-<a href="../../index.html">Back to the tutorial index</a>
+<a href="../index.html">Back to the tutorial index</a>
</p>
<hr>
<address>
--- /dev/null
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html>
+<head>
+ <title>DEAL tutorial: the Laplace problem</title>
+ <meta name="keyword" content="DEAL,DEAL tutorial">
+ <link href="../dealtut.css" rel="StyleSheet" media="screen" type="text/css" title="DEAL tutorial">
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="DEAL,DEAL tutorial">
+</head>
+<body>
+<h1>Solving the problem</h1>
+
+<hr>
+<p>
+<a href="../index.html">Back to the tutorial index</a>
+</p>
+
+<hr>
+<p>
+<address><a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
+</p>
+
+<p>
+Last modified: Fri Feb 12, 1999
+</p>
+</body>
+</html>
--- /dev/null
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html>
+<head>
+ <title>DEAL tutorial: the Laplace problem</title>
+ <meta name="keyword" content="DEAL,DEAL tutorial">
+ <link href="../dealtut.css" rel="StyleSheet" media="screen" type="text/css" title="DEAL tutorial">
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="DEAL,DEAL tutorial">
+</head>
+ <body>
+ <h1>Program structure</h1>
+ <p>
+ In order to solve our differential equation we have to deal with the
+ the following tasks:
+ </p>
+ <ul>
+ <li> <a href="triangulation.html">Grid generation</a>
+ </li>
+ <li> <a href="assemble.html">Assembly of the problem matrix</a>
+ </li>
+ <li> <a href="assemble.html#boundary">Description of the boundary conditions</a>
+ </li>
+ <li> <a href="solution.html">Solution of the equation system</a>
+ </li>
+ </ul>
+ <p>
+ All these tasks are problem specific, therefore we will put them into
+ a class called <a href="laplace.html"><em>Laplace</em></a>.
+ The <a href="main.html">main program</a>
+ will - for now - do nothing but call the methods of this class in their
+ proper order.
+ </p>
+ <p>
+ We end up with a program consisting of two major components:
+ <ul>
+ <li>a <em>main program</em> that starts the process of problem
+ generation and solution by calling methods of
+ </li>
+ <li>the <em>class Laplace</em> that contains all the parts specific
+ to the problem at hand
+ </li>
+ </ul>
+
+
+<hr>
+<p>
+<a href="../index.html">Back to the tutorial index</a>
+</p>
+
+<hr>
+<address><a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
+
+<p>
+Last modified: Mon Feb 15, 1999
+</p>
+</body>
+</html>
--- /dev/null
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html>
+<head>
+ <title>DEAL tutorial: the Laplace problem</title>
+ <meta name="keyword" content="DEAL,DEAL tutorial">
+ <link href="../dealtut.css" rel="StyleSheet" media="screen" type="text/css" title="DEAL tutorial">
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="DEAL,DEAL tutorial">
+</head>
+<body>
+<h1>Generating a triangulation</h1>
+
+The initial triangulation is generated by the constructor of the class <em>Laplace</em>.
+It is a hypercube with the coordinates along each axis ranging from -1 to 1. This first
+triangulation consists of only one finite element.
+<pre>
+<code>
+Laplace::Laplace(Function<2>& solution)
+{
+ tr.create_hypercube(-1.,1.);
+</code>
+</pre>
+<p>
+Right afterwards an iterator for the cells is defined (in effect a pointer to the first cell
+in this case) and all the cell's faces are numbered.
+</p>
+<pre>
+<code>
+ Triangulation<2>::cell_iterator c = tr.begin();
+ for (unsigned int i=0;i<GeometryInfo<2>::faces_per_cell;++i)
+ {
+ Triangulation<2>::face_iterator f = c->face(i);
+ f->set_boundary_indicator(i+1);
+ }
+</code>
+</pre>
+<p>
+Afterwards all the degrees of freedom, i.e. all the cell vertices, are renumbered.
+This step is necessary whenever the triangulation has changed, e.g. after each refinement.
+</p>
+<pre>
+<code>
+ dof_primal.distribute_dofs(fe_primal);
+ }
+</code>
+</pre>
+
+<hr>
+<p>
+<a href="../index.html">Back to the tutorial index</a>
+</p>
+
+<hr>
+<p>
+<address><a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
+</p>
+
+<p>
+Last modified: Fri Feb 12, 1999
+</p>
+</body>
+</html>
--- /dev/null
+<!-- DEAL Documentation Style Sheet.
+ Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de> 1999
+-->
+BODY {background: white}
+
--- /dev/null
+<!-- DEAL Tutorial Style Sheet
+ Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de 1999
+-->
+@import url(dealdoc.css);
--- /dev/null
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html>
+ <head>
+ <link href="dealdoc.css" rel="StyleSheet" type="text/css" title="DEAL Documentation" media="screen">
+ <title>The DEAL Tutorial</title>
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="DEAL,DEAL tutorial"></head>
+ <body>
+
+ <h1>The <acronym>DEAL</acronym> Tutorial</h1>
+
+ <h2>or <acronym>DEAL</acronym> in several easy and some more hard lessons</h2>
+
+ <h3>Prerequisites</h3>
+ <p>
+ This tutorial is written in <a href="http://www.w3.org/TR/REC-html40">HTML
+ 4.0</a> which your browser may or may not support. If you run into problems
+ you are advised to try a recent browser or one that is able to parse Document
+ Type Definitions. The reason for this is the math support of HTML 4.0. Most
+ people should, however, not run into any
+ problems.</p>
+ <p>
+ You yourself should have some background knowledge on numerical methods, finite
+ elements and C++.
+ </p>
+
+ <h3>Usage of this document</h3>
+ <p>
+ This tutorial is split into several parts, starting with the stationary
+ Laplace problem and advancing to more complicated problems. If you have never
+ used <acronym>DEAL</acronym> we recommend that you begin at the beginning and
+ advance through the chapters in their proper order. This tutorial contains the
+ complete code for every example program and can be used as a reference, too,
+ if you are looking for code to solve a specific problem with
+ <acronym>DEAL</acronym>.
+ </p>
+
+ <h3>On <acronym>DEAL</acronym></h3>
+ <p>
+ <acronym>DEAL</acronym> is short for Differential Equations Analysis Library,
+ which says it all, really. Its purpose it the solution of partial differential
+ equations on unstructured grids with error control. You can find more
+ information on <acronym>DEAL</acronym> and some examples at
+ <a href="http://gaia.iwr.uni-heidelberg.de/DEAL/">
+ http://gaia.iwr.uni-heidelberg.de/DEAL/</a>. <acronym>DEAL</acronym> was written
+ by the numerics group at the IWR, University of Heidelberg.
+ </p>
+
+ <h3><a href="01.laplace/index.html">Chapter 1: The Laplace
+ problem</a></h3>
+ <p>
+ We will solve the stationary Laplace problem on a square grid with zero
+ boundary condition on three of the edges and a constant boundary condition on
+ one edge.
+ </p>
+
+ <h3>Technical terms</h3>
+ <p>
+ There is a <a href="glossary.html">glossary</a> available to explain technical terms.
+ </p>
+
+ <hr>
+
+ <address>
+ <a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
+ <p>
+ Last modified: Tue 5 Jan 1999
+ </p>
+ </body>
+</html>
+
--- /dev/null
+#!/bin/sh
+echo Validating $1
+cat $1 | sed -e 's/http:\/\/www.w3.org\/TR\/REC-html40\/strict.dtd//' -e 's/\"\"//' | nsgmls -s -m /home/people/schrage/lib/sgml/pubtext/HTML4.soc