// @sect3{Include files}
- // The first few (many?) include
- // files have already been used in
- // the previous example, so we will
- // not explain their meaning here
- // again.
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <iostream>
#include <list>
- // This is new, however: in the previous
- // example we got some unwanted output from
- // the linear solvers. If we want to suppress
- // it, we have to include this file and add a
- // single line somewhere to the program (see
- // the main() function below for that):
-#include <deal.II/base/logstream.h>
-
- // The final step, as in previous
- // programs, is to import all the
- // deal.II class and function names
- // into the global namespace:
using namespace dealii;
- // @sect3{The <code>Step4</code> class template}
-
- // This is again the same
- // <code>Step4</code> class as in the
- // previous example. The only
- // difference is that we have now
- // declared it as a class with a
- // template parameter, and the
- // template parameter is of course
- // the spatial dimension in which we
- // would like to solve the Laplace
- // equation. Of course, several of
- // the member variables depend on
- // this dimension as well, in
- // particular the Triangulation
- // class, which has to represent
- // quadrilaterals or hexahedra,
- // respectively. Apart from this,
- // everything is as before.
+ // @sect3{The <code>Step41</code> class template}
+
+ // This class supply all function and variables
+ // to an obstacle problem. The projection_active_set
+ // function and the ConstaintMatrix are important
+ // for the handling of the active set as we see
+ // later.
+
template <int dim>
-class Step4
+class Step41
{
public:
- Step4 ();
+ Step41 ();
void run ();
private:
// @sect3{Right hand side and boundary values}
- // In the following, we declare two more
- // classes denoting the right hand side and
- // the non-homogeneous Dirichlet boundary
- // values. Both are functions of a
- // dim-dimensional space variable, so we
- // declare them as templates as well.
- //
- // Each of these classes is derived from a
- // common, abstract base class Function,
- // which declares the common interface which
- // all functions have to follow. In
- // particular, concrete classes have to
- // overload the <code>value</code> function,
- // which takes a point in dim-dimensional
- // space as parameters and shall return the
- // value at that point as a
- // <code>double</code> variable.
- //
- // The <code>value</code> function takes a
- // second argument, which we have here named
- // <code>component</code>: This is only meant
- // for vector valued functions, where you may
- // want to access a certain component of the
- // vector at the point
- // <code>p</code>. However, our functions are
- // scalar, so we need not worry about this
- // parameter and we will not use it in the
- // implementation of the functions. Inside
- // the library's header files, the Function
- // base class's declaration of the
- // <code>value</code> function has a default
- // value of zero for the component, so we
- // will access the <code>value</code>
- // function of the right hand side with only
- // one parameter, namely the point where we
- // want to evaluate the function. A value for
- // the component can then simply be omitted
- // for scalar functions.
- //
- // Note that the C++ language forces
- // us to declare and define a
- // constructor to the following
- // classes even though they are
- // empty. This is due to the fact
- // that the base class has no default
- // constructor (i.e. one without
- // arguments), even though it has a
- // constructor which has default
- // values for all arguments.
template <int dim>
class RightHandSide : public Function<dim>
{
const unsigned int component = 0) const;
};
-
-
template <int dim>
class BoundaryValues : public Function<dim>
{
// For this example, we choose as right hand
- // side function to function $4(x^4+y^4)$ in
- // 2D, or $4(x^4+y^4+z^4)$ in 3D. We could
- // write this distinction using an
- // if-statement on the space dimension, but
- // here is a simple way that also allows us
- // to use the same function in 1D (or in 4D,
- // if you should desire to do so), by using a
- // short loop. Fortunately, the compiler
- // knows the size of the loop at compile time
- // (remember that at the time when you define
- // the template, the compiler doesn't know
- // the value of <code>dim</code>, but when it later
- // encounters a statement or declaration
- // <code>RightHandSide@<2@></code>, it will take the
- // template, replace all occurrences of dim
- // by 2 and compile the resulting function);
- // in other words, at the time of compiling
- // this function, the number of times the
- // body will be executed is known, and the
- // compiler can optimize away the overhead
- // needed for the loop and the result will be
- // as fast as if we had used the formulas
- // above right away.
- //
- // The last thing to note is that a
- // <code>Point@<dim@></code> denotes a point in
- // dim-dimensionsal space, and its individual
- // components (i.e. $x$, $y$,
- // ... coordinates) can be accessed using the
- // () operator (in fact, the [] operator will
- // work just as well) with indices starting
- // at zero as usual in C and C++.
+ // side function a constant force density
+ // like the gravitation attraction.
template <int dim>
double RightHandSide<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
}
- // As boundary values, we choose x*x+y*y in
- // 2D, and x*x+y*y+z*z in 3D. This happens to
- // be equal to the square of the vector from
- // the origin to the point at which we would
- // like to evaluate the function,
- // irrespective of the dimension. So that is
- // what we return:
+ // As boundary values, we choose the zero.
template <int dim>
double BoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
return return_value;
}
+
+ // The obstacle function describes a cascaded
+ // barrier. So if the gravitation attraction
+ // pulls the membrane down it blows over the
+ // steps.
template <int dim>
double Obstacle<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
- // @sect3{Implementation of the <code>Step4</code> class}
-
- // Next for the implementation of the class
- // template that makes use of the functions
- // above. As before, we will write everything
- // as templates that have a formal parameter
- // <code>dim</code> that we assume unknown at
- // the time we define the template
- // functions. Only later, the compiler will
- // find a declaration of
- // <code>Step4@<2@></code> (in the
- // <code>main</code> function, actually) and
- // compile the entire class with
- // <code>dim</code> replaced by 2, a process
- // referred to as `instantiation of a
- // template'. When doing so, it will also
- // replace instances of
- // <code>RightHandSide@<dim@></code> by
- // <code>RightHandSide@<2@></code> and
- // instantiate the latter class from the
- // class template.
- //
- // In fact, the compiler will also find a
- // declaration
- // <code>Step4@<3@></code> in
- // <code>main()</code>. This will cause it to
- // again go back to the general
- // <code>Step4@<dim@></code>
- // template, replace all occurrences of
- // <code>dim</code>, this time by 3, and
- // compile the class a second time. Note that
- // the two instantiations
- // <code>Step4@<2@></code> and
- // <code>Step4@<3@></code> are
- // completely independent classes; their only
- // common feature is that they are both
- // instantiated from the same general
- // template, but they are not convertible
- // into each other, for example, and share no
- // code (both instantiations are compiled
- // completely independently).
-
-
- // @sect4{Step4::Step4}
-
- // After this introduction, here is the
- // constructor of the <code>Step4</code>
- // class. It specifies the desired polynomial
- // degree of the finite elements and
- // associates the DoFHandler to the
- // triangulation just as in the previous
- // example program, step-3:
+ // @sect3{Implementation of the <code>Step41</code> class}
+
+
+ // @sect4{Step41::Step41}
+
template <int dim>
-Step4<dim>::Step4 ()
+Step41<dim>::Step41 ()
:
fe (1),
dof_handler (triangulation)
{}
- // @sect4{Step4::make_grid}
-
- // Grid creation is something inherently
- // dimension dependent. However, as long as
- // the domains are sufficiently similar in 2D
- // or 3D, the library can abstract for
- // you. In our case, we would like to again
- // solve on the square $[-1,1]\times [-1,1]$
- // in 2D, or on the cube $[-1,1] \times
- // [-1,1] \times [-1,1]$ in 3D; both can be
- // termed GridGenerator::hyper_cube(), so we may
- // use the same function in whatever
- // dimension we are. Of course, the functions
- // that create a hypercube in two and three
- // dimensions are very much different, but
- // that is something you need not care
- // about. Let the library handle the
- // difficult things.
+ // @sect4{Step41::make_grid}
+
+ // We solve our obstacle problem on the square
+ // $[-1,1]\times [-1,1]$ in 2D.
template <int dim>
-void Step4<dim>::make_grid ()
+void Step41<dim>::make_grid ()
{
GridGenerator::hyper_cube (triangulation, -1, 1);
n_refinements = 6;
<< std::endl;
}
- // @sect4{Step4::setup_system}
-
- // This function looks
- // exactly like in the previous example,
- // although it performs actions that in their
- // details are quite different if
- // <code>dim</code> happens to be 3. The only
- // significant difference from a user's
- // perspective is the number of cells
- // resulting, which is much higher in three
- // than in two space dimensions!
+ // @sect4{Step41::setup_system}
+
template <int dim>
-void Step4<dim>::setup_system ()
+void Step41<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
}
- // @sect4{Step4::assemble_system}
+ // @sect4{Step41::assemble_system}
- // Unlike in the previous example, we
- // would now like to use a
- // non-constant right hand side
- // function and non-zero boundary
- // values. Both are tasks that are
- // readily achieved with a only a few
- // new lines of code in the
- // assemblage of the matrix and right
- // hand side.
- //
- // More interesting, though, is the
- // way we assemble matrix and right
- // hand side vector dimension
- // independently: there is simply no
- // difference to the
- // two-dimensional case. Since the
- // important objects used in this
- // function (quadrature formula,
- // FEValues) depend on the dimension
- // by way of a template parameter as
- // well, they can take care of
- // setting up properly everything for
- // the dimension for which this
- // function is compiled. By declaring
- // all classes which might depend on
- // the dimension using a template
- // parameter, the library can make
- // nearly all work for you and you
- // don't have to care about most
- // things.
+
+ // At once with assembling the system matrix and
+ // right-hand-side we apply the constraints
+ // to our system. The constraint consists not
+ // only of the zero Dirichlet boundary values,
+ // in addition they contain the obstacle values.
+ // The projection_active_set function are used
+ // to fill the ConstraintMatrix.
template <int dim>
-void Step4<dim>::assemble_system ()
+void Step41<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(2);
- // We wanted to have a non-constant right
- // hand side, so we use an object of the
- // class declared above to generate the
- // necessary data. Since this right hand
- // side object is only used locally in the
- // present function, we declare it here as
- // a local variable:
const RightHandSide<dim> right_hand_side;
- // Compared to the previous example, in
- // order to evaluate the non-constant right
- // hand side function we now also need the
- // quadrature points on the cell we are
- // presently on (previously, we only
- // required values and gradients of the
- // shape function from the
- // FEValues object, as well as
- // the quadrature weights,
- // FEValues::JxW() ). We can tell the
- // FEValues object to do for
- // us by also giving it the
- // #update_quadrature_points
- // flag:
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
- // We then again define a few
- // abbreviations. The values of these
- // variables of course depend on the
- // dimension which we are presently
- // using. However, the FE and Quadrature
- // classes do all the necessary work for
- // you and you don't have to care about the
- // dimension dependent parts:
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- // Next, we again have to loop over all
- // cells and assemble local contributions.
- // Note, that a cell is a quadrilateral in
- // two space dimensions, but a hexahedron
- // in 3D. In fact, the
- // <code>active_cell_iterator</code> data
- // type is something different, depending
- // on the dimension we are in, but to the
- // outside world they look alike and you
- // will probably never see a difference
- // although the classes that this typedef
- // stands for are in fact completely
- // unrelated:
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
cell_matrix = 0;
cell_rhs = 0;
- // Now we have to assemble the
- // local matrix and right hand
- // side. This is done exactly
- // like in the previous
- // example, but now we revert
- // the order of the loops
- // (which we can safely do
- // since they are independent
- // of each other) and merge the
- // loops for the local matrix
- // and the local vector as far
- // as possible to make
- // things a bit faster.
- //
- // Assembling the right hand side
- // presents the only significant
- // difference to how we did things in
- // step-3: Instead of using a constant
- // right hand side with value 1, we use
- // the object representing the right
- // hand side and evaluate it at the
- // quadrature points:
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
right_hand_side.value (fe_values.quadrature_point (q_point)) *
fe_values.JxW (q_point));
}
- // As a final remark to these loops:
- // when we assemble the local
- // contributions into
- // <code>cell_matrix(i,j)</code>, we
- // have to multiply the gradients of
- // shape functions $i$ and $j$ at point
- // q_point and multiply it with the
- // scalar weights JxW. This is what
- // actually happens:
- // <code>fe_values.shape_grad(i,q_point)</code>
- // returns a <code>dim</code>
- // dimensional vector, represented by a
- // <code>Tensor@<1,dim@></code> object,
- // and the operator* that multiplies it
- // with the result of
- // <code>fe_values.shape_grad(j,q_point)</code>
- // makes sure that the <code>dim</code>
- // components of the two vectors are
- // properly contracted, and the result
- // is a scalar floating point number
- // that then is multiplied with the
- // weights. Internally, this operator*
- // makes sure that this happens
- // correctly for all <code>dim</code>
- // components of the vectors, whether
- // <code>dim</code> be 2, 3, or any
- // other space dimension; from a user's
- // perspective, this is not something
- // worth bothering with, however,
- // making things a lot simpler if one
- // wants to write code dimension
- // independently.
-
- // With the local systems assembled,
- // the transfer into the global matrix
- // and right hand side is done exactly
- // as before, but here we have again
- // merged some loops for efficiency:
- cell->get_dof_indices (local_dof_indices);
+
+ cell->get_dof_indices (local_dof_indices);
+
+ // This function apply the constraints
+ // to the system matrix and system rhs.
+ // The true parameter is set to make sure
+ // that the system rhs contains correct
+ // values in the rows with inhomogeneity
+ // constraints.
constraints.distribute_local_to_global (cell_matrix, cell_rhs,
local_dof_indices,
system_matrix, system_rhs, true);
}
}
- // @sect4{Step4::projection_active_set}
+ // @sect4{Step41::projection_active_set}
- // Projection and updating of the active set
- // for the dofs which penetrates the obstacle.
+ // Updating of the active set which means to
+ // set a inhomogeneity constraint in the
+ // ConstraintMatrix. At the same time we set
+ // the solution to the correct value - the obstacle value.
+ // To control the active set we use the vector
+ // active_set which contains a zero in a component
+ // that is not in the active set and elsewise a
+ // one. With the output file you can visualize it.
template <int dim>
-void Step4<dim>::projection_active_set ()
+void Step41<dim>::projection_active_set ()
{
const Obstacle<dim> obstacle;
std::vector<bool> vertex_touched (triangulation.n_vertices(),
double obstacle_value = obstacle.value (point);
double solution_index_x = solution (index_x);
- // to decide which dof belongs to the
- // active-set. for that we scale the
+ // To decide which dof belongs to the
+ // active-set. For that we scale the
// residual-vector with the cell-size and
// the diag-entry of the mass-matrix.
+
+ // TODO: I have to check the condition
if ((resid_vector (index_x)*std::pow (2, 2*n_refinements)*diag_mass_matrix_vector (index_x) >= solution_index_x - obstacle_value))
{
constraints.add_line (index_x);
}
std::cout<< "Number of Contact-Constaints: " << counter_contact_constraints <<std::endl;
- // to supply the boundary values of the
+ // To supply the boundary values of the
// dirichlet-boundary in constraints
VectorTools::interpolate_boundary_values (dof_handler,
0,
constraints.close ();
}
- // @sect4{Step4::solve}
+ // @sect4{Step41::solve}
- // Solving the linear system of
- // equations is something that looks
- // almost identical in most
- // programs. In particular, it is
- // dimension independent, so this
- // function is copied verbatim from the
- // previous example.
template <int dim>
-void Step4<dim>::solve ()
+void Step41<dim>::solve ()
{
ReductionControl reduction_control (100, 1e-12, 1e-2);
SolverCG<TrilinosWrappers::Vector> solver (reduction_control);
<< std::endl;
}
- // @sect4{Step4::output_results}
+ // @sect4{Step41::output_results}
- // This function also does what the
- // respective one did in step-3. No changes
- // here for dimension independence either.
- //
- // The only difference to the previous
- // example is that we want to write output in
- // VTK format, rather than for gnuplot. VTK
- // format is currently the most widely used
- // one and is supported by a number of
- // visualization programs such as Visit and
- // Paraview (for ways to obtain these
- // programs see the ReadMe file of
- // deal.II). To write data in this format, we
- // simply replace the
- // <code>data_out.write_gnuplot</code> call
- // by <code>data_out.write_vtk</code>.
- //
- // Since the program will run both 2d and 3d
- // versions of the laplace solver, we use the
- // dimension in the filename to generate
- // distinct filenames for each run (in a
- // better program, one would check whether
- // <code>dim</code> can have other values
- // than 2 or 3, but we neglect this here for
- // the sake of brevity).
+ // We use the vtk-format for the output.
+ // The file contains the displacement,
+ // the residual and active set vectors.
template <int dim>
-void Step4<dim>::output_results (const std::string& title) const
+void Step41<dim>::output_results (const std::string& title) const
{
DataOut<dim> data_out;
data_out.build_patches ();
- std::ofstream output_vtk (dim == 2 ?
- (title + ".vtk").c_str () :
- (title + ".vtk").c_str ());
+ std::ofstream output_vtk ((title + ".vtk").c_str ());
data_out.write_vtk (output_vtk);
}
- // @sect4{Step4::run}
+ // @sect4{Step41::run}
// This is the function which has the
- // top-level control over
- // everything. Apart from one line of
- // additional output, it is the same
- // as for the previous example.
+ // top-level control over everything.
+ // Here the active set method is implemented.
+
+ // TODO: I have to compare it with the algorithm
+ // in the Wohlmuth-paper
template <int dim>
-void Step4<dim>::run ()
+void Step41<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
// looks mostly like in step-3, but if you
// look at the code below, note how we first
// create a variable of type
- // <code>Step4@<2@></code> (forcing
+ // <code>Step41@<2@></code> (forcing
// the compiler to compile the class template
// with <code>dim</code> replaced by
// <code>2</code>) and run a 2d simulation,
Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
- Step4<2> laplace_problem_2d;
+ Step41<2> laplace_problem_2d;
laplace_problem_2d.run ();
return 0;