laplace.cc File Reference

#include <amandus/apps.h>
#include <amandus/laplace/matrix.h>
#include <amandus/laplace/matrix_factor.h>
#include <amandus/laplace/noforce.h>
#include <amandus/laplace/polynomial.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/numerics/dof_output_operator.h>
#include <deal.II/numerics/dof_output_operator.templates.h>
#include <deal.II/base/function.h>
#include <deal.II/algorithms/newton.h>

Classes
class	Startup< dim >

Functions
int	main ()

Detailed Description

• About this example

This example is the easiest program you can write in Amandus. It shows how Amandus work in practice and it gives a brief description of the functions used. It doesn't discuss in detail any individual function but it wants to give a big picture of how things work together. If you wants to learn deeply how a single function work, you can search inside the manual. For informations about deal.II functions and classes, use deal.II manual . As concerns the Amandus manual, it is not complete and it is still in development.

How to compile and run

Go to example/laplace directory and type:

make laplace_laplace

Now you can see the executable. To execute it and produce the output, type:

./laplace

Introduction

We solve a simple version of Poisson's equation with nonzero right hand side:

\begin{align*} - \Delta u & = f \qquad\qquad & \text{in}\ \Omega \end{align*}

We solve this equation on the unit square, \(\Omega=[0,1]^2\), for which we already know how to generate a mesh from dealii documentation. In this program, we also only consider the particular case \(f(\mathbf x)=1\).

We choose the following Dirichlet boundary conditions:

\begin{align*} u(x) &= x \qquad\qquad & \text{on}\ \partial\Omega. \end{align*}

From the basics of the finite element method, you remember the steps we need to take to approximate the solution \(u\) by a finite dimensional approximation. Multiplying from the left by the test function \(\varphi \), integrating over the domain \(\Omega\) and integrating by parts, we get the following weak formulation:

\begin{align*} \int_\Omega \nabla\varphi \cdot \nabla u - \int_{\partial\Omega} \varphi \mathbf{n}\cdot \nabla u = \int_\Omega \varphi f. \end{align*}

We build a discreate approximation of u:

\begin{align*} u_h(\mathbf x)=\sum_j U_j \varphi_j(\mathbf x) \end{align*}

where \(\varphi_i(\mathbf x)\) are the finite element shape functions we will use.

The weak form of the discrete problem is : Find a function \(u_h\), i.e. find the expansion coefficients \(U_i\) mentioned above, so that

\begin{align*} (\nabla\varphi_i, \nabla u_h) = (\varphi_i, f), \qquad\qquad i=0\ldots N-1. \end{align*}

This equation can be rewritten as a linear system by inserting the representation \(u_h(\mathbf x)=\sum_j U_j \varphi_j(\mathbf x)\): Find a vector \(U\) so that

\begin{align*} A U = F, \end{align*}

where the matrix \(A\) and the right hand side \(F\) are defined as

\begin{align*} A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j), \\ F_i &= (\varphi_i, f). \end{align*}

If you need more details about how to construct a linear system starting from the Poisson's equation, you might be interested in step-3 of the deal.II tutorial.

Now we can go forward and describe how this quantities can be computed with Amandus.

Function Documentation

int main

(

)

We fix the dimension of the domain

const unsigned int d=2;

Construction of a ofstream and the attachment of the stream to the open file deallog

std::ofstream logfile("deallog");
deallog.attach(logfile);

Now we define an object for a triangulation and we fill it with a single cell of a square domain. The triangulation is refined 3 times, to yield \(4^3=64\) cells in total.

Triangulation<d> tr;
GridGenerator::hyper_cube (tr, -1, 1);
tr.refine_global(3);

We construct the tensor product polynomials of degree p. We use continuous polynomials.

const unsigned int degree = 2;
FE_Q<d> fe(degree);

Creation of a Matrix object. The Matrix class permits to construct the matrix A using LocalIntegrators functions. It is essential to implement a class matrix for every kind of problem you want to solve. For more details, read the documentation for the matrix class of the Laplace problem.

LaplaceIntegrators::Matrix<d> matrix_integrator;

Todo:

The code does't compile. Correct it and finish the code description.

Here is the call graph for this function: