--- /dev/null
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1999
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly. We get deduce it from the files in the present
+# directory:
+target = $(basename $(shell echo step-*.cc))
+
+# 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
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+
+
+# list of libraries needed to link with
+libs = -ldeal_II_2d -llac -lbase
+libs.g = -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+ libraries = $(target).go $(libs.g)
+ flags = $(CXXFLAGS.g)
+else
+ libraries = $(target).go $(libs)
+ flags = $(CXXFLAGS.o)
+endif
+
+
+
+# make rule for the target. $^ is the object file $(target).g?o
+$(target) : $(libraries)
+ @echo ============================ Linking $@
+ @$(CXX) $(flags) -o $@ $^
+
+# rule how to run the program
+run: $(target)
+ @echo ============================ Running $<
+ @./$(target)
+
+
+# rule to make object files
+%.go : %.cc
+ @echo ============================ Compiling with debugging information: $<
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+ @echo ============================ Compiling with optimization: $<
+ @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+clean:
+ -rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps
+
+
+
+.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.
+#
+# Since the script prefixes the output names by lib/g?o, we have to
+# strip that again (the script was written for the main libraries and
+# large projects where object files are put into subdirs)
+Makefile.dep: $(target).cc Makefile \
+ $(shell echo $D/base/include/base/*.h \
+ $D/lac/include/lac/*.h \
+ $D/deal.II/include/*/*.h)
+ @echo ============================ Remaking Makefile
+ @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
+ | perl -pi -e 's!lib/g?o/!!g;' \
+ > Makefile.dep
+
+
+include Makefile.dep
+
--- /dev/null
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+
+ // As usual, the first few include
+ // files are already known, so we
+ // will not comment on them further.
+#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 <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 <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+ // In this example, we need
+ // vector-valued finite elements. The
+ // suuport for these can be found in
+ // the following include file:
+#include <fe/fe_system.h>
+ // We will compose the vector-valued
+ // finite elements from regular Q1
+ // elements which can be found here,
+ // as usual:
+#include <fe/fe_lib.lagrange.h>
+
+ // This again is C++:
+#include <fstream>
+
+
+ // The main class is, except for its
+ // name, almost unchanged with
+ // respect to the step-6 example. The
+ // only change is the use of a
+ // different class for the ``fe''
+ // variable.
+template <int dim>
+class ElasticProblem
+{
+ public:
+ ElasticProblem ();
+ ~ElasticProblem ();
+ 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;
+ DoFHandler<dim> dof_handler;
+
+ // Instead of a concrete finite
+ // element class such as
+ // ``FEQ1'', we now use a more
+ // generic one, ``FESystem''. In
+ // fact, it is not a finite
+ // element itself, but rather a
+ // class that can be used to
+ // stack several usual elements
+ // together to form one
+ // vector-valued finite
+ // element. In our case, we will
+ // compose the vector-valued
+ // element of ``FEQ1'' objects,
+ // as shown below in the
+ // constructor of this class.
+ FESystem<dim> fe;
+
+ ConstraintMatrix hanging_node_constraints;
+
+ SparseMatrixStruct sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+ // Before going over to the
+ // implementation of the main class,
+ // we declare and define the class
+ // which describes the right hand
+ // side. This time, the right hand
+ // side is vector-valued, as is the
+ // solution, so we will describe the
+ // new elements in some more detail.
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+ public:
+ // The first thing is that
+ // vector-valued functions have a
+ // constructor, since they need
+ // to pass down to the base class
+ // of how many components the
+ // function consists. The default
+ // value in the constructor of
+ // the base class is one, so we
+ // need not define a constructor
+ // for the usual scalar function.
+ RightHandSide ();
+
+ // The next function is a
+ // replacement for the ``value''
+ // function of the previous
+ // examples. There, a second
+ // parameter ``component'' was
+ // given, which denoted which
+ // component was requested. Here,
+ // we implement a function that
+ // returns the whole vector of
+ // values at the given place at
+ // once.
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &values) const;
+
+ // Then, in analogy to the
+ // ``value_list'' function, there
+ // is a function
+ // ``vector_value_list'', which
+ // returns the values of the
+ // vector-valued function at
+ // several points at once:
+ virtual void vector_value_list (const vector<Point<dim> > &points,
+ vector<Vector<double> > &value_list) const;
+};
+
+
+ // This is the constructor of the
+ // right hand side class. As said
+ // above, it only passes down to the
+ // base class the number of
+ // components, which is ``dim'' in
+ // the present case. Note that
+ // although the implementation is
+ // very short here, we do not move it
+ // into the class declaration, since
+ // our style guides require that
+ // inside the class declaration only
+ // declarations have to happen and
+ // that definitions are always to be
+ // found outside.
+template <int dim>
+RightHandSide<dim>::RightHandSide () :
+ Function<dim> (dim)
+{};
+
+
+ // This is the function that returns
+ // the whole vector of values at the
+ // point ``p'' at once:
+template <int dim>
+inline
+void RightHandSide<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+{
+ // To prevent cases where the
+ // return value has not previously
+ // been set to the right size
+ // (which is kind of a convention
+ // in the deal.II library), we test
+ // for this case and otherwise
+ // throw an exception:
+ Assert (values.size() == dim,
+ ExcVectorHasWrongSize (values.size(), dim));
+ // Likewise, if by some accident
+ // someone tried to compile and run
+ // the program in only one space
+ // dimension (in which the elastic
+ // equations do not make much sense
+ // since they reduce to the
+ // ordinary Laplace equation), we
+ // terminate the program if the
+ // dimension is not as expected.
+ Assert (dim >= 2, ExcInternalError());
+
+ // The rest of the function is as
+ // would probably be expected given
+ // the form of the right hand side
+ // function. First we define the
+ // centers of the two points around
+ // which are the sources of
+ // x-displacement, i.e. (0.5,0) and
+ // (-0.5,0). Note that upon
+ // construction of the ``Point''
+ // objects, all components are set
+ // to zero.
+ Point<dim> point_1, point_2;
+ point_1(0) = 0.5;
+ point_2(0) = -0.5;
+
+ // If now the point ``p'' is in the
+ // circle of radius 0.2 around one
+ // of these points, then set the
+ // force in x-direction to one,
+ // otherwise to zero:
+ if (((p-point_1).square() < 0.2*0.2) ||
+ ((p-point_2).square() < 0.2*0.2))
+ values(0) = 1;
+ else
+ values(0) = 0;
+
+ // Likewise, if ``p'' is in the
+ // vicinity of the origin, then set
+ // the y-force to 1, otherwise to
+ // zero:
+ if (p.square() < 0.2*0.2)
+ values(1) = 1;
+ else
+ values(1) = 0;
+};
+
+
+
+ // Now, this is the function of the
+ // right hand side class that returns
+ // the values at several points at
+ // once.
+template <int dim>
+void RightHandSide<dim>::vector_value_list (const vector<Point<dim> > &points,
+ vector<Vector<double> > &value_list) const
+{
+ // First we define an abbreviation
+ // for the number of points which
+ // we shall work on:
+ const unsigned int n_points = points.size();
+
+ // Then we check whether the number
+ // of output slots has been set
+ // correctly, i.e. to the number of
+ // input points:
+ Assert (value_list.size() == n_points,
+ ExcVectorHasWrongSize (value_list.size(), n_points));
+
+ // Finally we treat each of the
+ // points. In one of the previous
+ // examples, we have explained why
+ // the
+ // ``value_list''/``vector_value_list''
+ // function had been introduced: to
+ // prevent us from calling virtual
+ // functions too frequently. On the
+ // other hand, we now need to
+ // implement the same function
+ // twice, which can lead to
+ // confusion if one function is
+ // changed but the other is
+ // not. However, we can prevent
+ // this situation using the
+ // following construct:
+ for (unsigned int p=0; p<n_points; ++p)
+ RightHandSide<dim>::vector_value (points[p],
+ value_list[p]);
+ // It calls the ``vector_value''
+ // function defined above for each
+ // point, and thus preempts all
+ // chances for inconsistency. It is
+ // important to note how the
+ // function was called: using the
+ // full class qualification using
+ // ``RightHandSide::'', since this
+ // calls the function directly and
+ // not using the virtual function
+ // table. The call is thus as fast
+ // as a call to any non-virtual
+ // function. In addition, we have
+ // declared the ``vector_value''
+ // function ``inline'', i.e. the
+ // compiler can remove the function
+ // call altogether and the
+ // resulting code can in principle
+ // be as fast as if we had
+ // duplicated the code.
+};
+
+
+
+
+template <int dim>
+ElasticProblem<dim>::ElasticProblem () :
+ dof_handler (triangulation),
+ // As said before, we
+ // would like to
+ // construct one
+ // vector-valued
+ // finite element as
+ // outer product of
+ // several scala
+ // finite
+ // elements. Of
+ // course, the number
+ // of scalar finite
+ // element we would
+ // like to stack
+ // together equals
+ // the number of
+ // components the
+ // solution function
+ // has, which is
+ // ``dim'' since we
+ // consider
+ // displacement in
+ // each space
+ // direction. The
+ // ``FESystem'' class
+ // can handle this:
+ // we pass it the
+ // finite element of
+ // which we would
+ // like to compose
+ // the system of, and
+ // how often it shall
+ // be repeated:
+ fe (FEQ1<dim>(), dim)
+ // In fact, the ``FESystem'' class
+ // has several more constructors
+ // which can perform more complex
+ // operations that just stacking
+ // together several scalar finite
+ // elements of the same type into
+ // one; we will get to know these
+ // possibilities in later examples.
+ //
+ // It should be noted that the
+ // ``FESystem'' object thus created
+ // does not actually use the finite
+ // element which we have passed to it
+ // as first parameter. We could thus
+ // use an anonymous object created
+ // in-place. The ``FESystem''
+ // constructor only needs the
+ // parameter to deduce the type of
+ // the finite element from this and
+ // then creates objects of the
+ // underlying finite element type
+ // itself.
+{};
+
+
+
+template <int dim>
+ElasticProblem<dim>::~ElasticProblem ()
+{
+ dof_handler.clear ();
+};
+
+
+ // Setting up the system of equations
+ // is equal to the function used in
+ // the step-6 example. The
+ // ``DoFHandler'' class and all other
+ // classes used take care of the
+ // vector-valuedness of the finite
+ // element themselves (in fact, the
+ // do not do so, since they only take
+ // care how many degrees of freedom
+ // there are per vertex, line and
+ // cell, and they do not askwhat they
+ // represent, i.e. whether the finite
+ // element under consideration is
+ // vector-valued or whether it is,
+ // for example, a scalar Hermite
+ // element with several degrees of
+ // freedom on each vertex).
+template <int dim>
+void ElasticProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ // When making the sparsity
+ // pattern, there is some potential
+ // for optimization if not all
+ // components couple to all
+ // others. However, this is not the
+ // case for the elastic equations,
+ // so we use the standard call:
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+ hanging_node_constraints.condense (sparsity_pattern);
+
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+};
+
+
+ // The big changes in this program
+ // are in the creation of matrix and
+ // right hand side, since they are
+ // problem-dependent. We will go
+ // through that process step-by-step,
+ // since it is a bit more complicated
+ // than in previous examples.
+template <int dim>
+void ElasticProblem<dim>::assemble_system ()
+{
+ // First thing: the quadrature
+ // formula does not need
+ // modification since we still deal
+ // with bilinear functions.
+ QGauss2<dim> quadrature_formula;
+ // Also, the ``FEValues'' objects
+ // takes care of everything for us
+ // (or better: it does not really
+ // so; as in the comment in the
+ // function setting up the system,
+ // here as well the ``FEValues''
+ // object computes the same data on
+ // each cell, but it has some
+ // functionality to access data
+ // stored inside the finite element
+ // where they are precomputed upon
+ // construction).
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
+
+ // The number of degrees of freedom
+ // per cell we now obviously ask
+ // from the composed finite element
+ // rather than from the underlying
+ // scalar Q1 element. Here, it is
+ // ``dim'' times the number of
+ // degrees of freedom per cell of
+ // the Q1 element, but this is not
+ // something we need to care about.
+ 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<int> local_dof_indices (dofs_per_cell);
+
+ // As was shown in previous
+ // examples as well, we need a
+ // place where to store the values
+ // of the coefficients at all the
+ // quadrature points on a cell. In
+ // the present situation, we have
+ // two coefficients, lambda and mu.
+ vector<double> lambda_values (n_q_points);
+ vector<double> mu_values (n_q_points);
+
+ // Well, we could as well have
+ // omitted the above two arrays
+ // since we will use constant
+ // coefficients for both lambda and
+ // mu, which can be declared like
+ // this. They both represent
+ // functions always returning the
+ // constant value 1.0. Although we
+ // could omit the respective
+ // factors in the assemblage of the
+ // matrix, we use them here for
+ // purpose of demonstration.
+ ConstantFunction<dim> lambda(1.), mu(1.);
+
+ // Then again, we need to have the
+ // same for the right hand
+ // side. This is exactly as before
+ // in previous examples. However,
+ // we now have a vector-valued
+ // right hand side, which is why
+ // the data type of the
+ // ``rhs_values'' array is
+ // changed. We initialize it by
+ // ``n_q_points'' elements, each of
+ // which is a ``Vector<double>''
+ // with ``dim'' elements.
+ RightHandSide<dim> right_hand_side;
+ vector<Vector<double> > rhs_values (n_q_points,
+ Vector<double>(dim));
+
+
+ // Now we can begin with the loop
+ // over all cells:
+ DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = 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();
+
+ lambda.value_list (q_points, lambda_values);
+ mu.value_list (q_points, mu_values);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int component_i = fe.system_to_component_index(i).first;
+
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const unsigned int component_j = fe.system_to_component_index(j).first;
+
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ {
+ cell_matrix(i,j)
+ +=
+ // (lambda d_i u_i, d_j v_j)
+ (
+ (shape_grads[i][q_point][component_i] *
+ shape_grads[j][q_point][component_j] *
+ lambda_values[q_point])
+ + // (mu d_i u_j, d_i v_j)
+ (shape_grads[i][q_point][component_j] *
+ shape_grads[j][q_point][component_i] *
+ mu_values[q_point])
+ + // (mu d_i v_j, d_i v_j)
+ ((component_i == component_j) ?
+ (shape_grads[i][q_point] *
+ shape_grads[j][q_point] *
+ mu_values[q_point]) :
+ 0)
+ )
+ *
+ JxW_values[q_point];
+ };
+ };
+ };
+
+
+ right_hand_side.vector_value_list (q_points, rhs_values);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int component_i = fe.system_to_component_index(i).first;
+
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ cell_rhs(i) += shape_values(i,q_point) *
+ rhs_values[q_point](component_i) *
+ JxW_values[q_point];
+ };
+
+ // The transfer from local
+ // degrees of freedom into the
+ // global matrix and right hand
+ // side vector does not depend
+ // on the equation under
+ // consideration, and is thus
+ // the same as in all previous
+ // examples.
+ 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)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ };
+ };
+
+ // The interpolation of the
+ // boundary values needs a small
+ // modification: since the solution
+ // function is vector-valued, so
+ // needs to be the boundary
+ // values. The ``ZeroFunction''
+ // constructor accepts a parameter
+ // that tells it that it shall
+ // represent a vector valued,
+ // constant zero function with that
+ // many components. By default,
+ // this parameter is equal to one,
+ // in which case the
+ // ``ZeroFunction'' object would
+ // represent a scalar
+ // function. Since the solution
+ // vector has ``dim'' components,
+ // we need to pass ``dim'' as
+ // number of components to the zero
+ // function as well.
+ map<int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(dim),
+ boundary_values);
+ MatrixTools<dim>::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+
+ hanging_node_constraints.condense (system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+};
+
+
+
+ // The solver does not care about
+ // where the system of equations
+ // comes, as long as it stays
+ // positive definite and symmetric
+ // (which are the requirements for
+ // the use of the CG solver), which
+ // the system is. Therefore, we need
+ // not change anything.
+template <int dim>
+void ElasticProblem<dim>::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+
+ PreconditionRelaxation<>
+ preconditioner(system_matrix,
+ &SparseMatrix<double>::template precondition_SSOR<double>,
+ 1.2);
+
+ cg.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+};
+
+
+
+ // The function that does the
+ // refinement of the grid is the same
+ // as in the step-6 example. The
+ // quadrature formula is adapted to
+ // the linear elements again. Note
+ // that the error estimator by
+ // default adds up the estimated
+ // obtained from all components of
+ // the finite element solution, that
+ // is it uses the displacement in all
+ // directions with the same
+ // weight. If we would like the grid
+ // to be adapted to the
+ // x-displacement only, we could pass
+ // the function an additional
+ // parameter which tells it to do so
+ // and do not consider the
+ // displacements in all other
+ // directions for the error
+ // indicators.
+template <int dim>
+void ElasticProblem<dim>::refine_grid ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss2<dim-1>(),
+ neumann_boundary,
+ solution,
+ estimated_error_per_cell);
+
+ triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
+ 0.3, 0.03);
+
+ triangulation.execute_coarsening_and_refinement ();
+};
+
+
+ // The output happens mostly as has
+ // been shown in previous examples
+ // already. The only difference is
+ // not that the solution function is
+ // vector values. The ``DataOut''
+ // class takes care of this
+ // automatically, but we have to give
+ // each component of the solution
+ // vector a different name.
+template <int dim>
+void ElasticProblem<dim>::output_results (const unsigned int cycle) const
+{
+ string filename = "solution-";
+ filename += ('0' + cycle);
+ Assert (cycle < 10, ExcInternalError());
+
+ filename += ".gmv";
+ ofstream output (filename.c_str());
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+
+
+
+ // As said above, we need a
+ // different name for each
+ // component of the solution
+ // function. To pass one name for
+ // each component, a vector of
+ // strings is used. Since the
+ // number of components is the same
+ // as the number of dimensions we
+ // are working in, the following
+ // ``switch'' statement is used.
+ //
+ // We note that some graphics
+ // programs have restriction as to
+ // what characters are allowed in
+ // the names of variables. The
+ // library therefore supports only
+ // the minimal subset of these
+ // characters that is supported by
+ // all programs. Basically, these
+ // are letters, numbers,
+ // underscores, and some other
+ // characters, but in particular no
+ // whitespace and minus/hyphen. The
+ // library will throw an exception
+ // otherwise, at least if in debug
+ // mode.
+ vector<string> solution_names;
+ switch (dim)
+ {
+ case 1:
+ solution_names.push_back ("displacement");
+ break;
+ case 2:
+ solution_names.push_back ("x_displacement");
+ solution_names.push_back ("y_displacement");
+ break;
+ case 3:
+ solution_names.push_back ("x_displacement");
+ solution_names.push_back ("y_displacement");
+ solution_names.push_back ("z_displacement");
+ break;
+ // It is good style to
+ // let the program die if
+ // we run upon a case
+ // which we did not
+ // consider. Remember
+ // that the ``Assert''
+ // macro throws an
+ // exception if the
+ // condition in the first
+ // parameter is not
+ // satisfied. Of course,
+ // the condition
+ // ``false'' can never be
+ // satisfied, so the
+ // program will always
+ // abort whenever it gets
+ // to this statement:
+ default:
+ Assert (false, ExcInternalError());
+ };
+
+ // After setting up the names for
+ // the different components of the
+ // solution vector, we can add the
+ // solution vector to the list of
+ // data vectors scheduled for
+ // output. Note that the following
+ // function takes a vector of
+ // strings as second argument,
+ // whereas the one which we have
+ // used in all previous examples
+ // accepted a string there. In
+ // fact, the latter function is
+ // only a shortcut for the function
+ // which we call here: it puts the
+ // single string that is passed to
+ // it into a vector of strings with
+ // only one element and forwards
+ // that to the other function.
+ data_out.add_data_vector (solution, solution_names);
+ data_out.build_patches ();
+ data_out.write_gmv (output);
+};
+
+
+
+template <int dim>
+void ElasticProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<8; ++cycle)
+ {
+ cout << "Cycle " << cycle << ':' << endl;
+
+ if (cycle == 0)
+ {
+ // As in previous examples,
+ // we use the unit square
+ // (or cube) as domain.
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ // This time, we have to
+ // refine the coarse grid
+ // twice before we first
+ // solve on it. The reason
+ // is the following: we use
+ // the ``Gauss2''
+ // quadrature formula for
+ // integration of the right
+ // hand side; that means
+ // that there are four
+ // quadrature points on
+ // each cell (in 2D). If we
+ // only refine the initial
+ // grid once globally, then
+ // there will be only four
+ // quadrature points in
+ // each direction on the
+ // domain. However, the
+ // right hand side function
+ // was chosen to be rather
+ // localized and in that
+ // case all quadrature
+ // points lie outside the
+ // support of the right
+ // hand side function. The
+ // right hand side vector
+ // will then contain only
+ // zeroes and the solution
+ // of the system of
+ // equations is the zero
+ // vector, i.e. a finite
+ // element function that it
+ // zero everywhere. We
+ // should not be surprised
+ // about such things
+ // happening, since we have
+ // chosen an initial grid
+ // that is totally
+ // unsuitable for the
+ // problem at hand.
+ //
+ // The unfortunate thing is
+ // that if the discrete
+ // solution is constant,
+ // then the error
+ // indicators computed by
+ // the
+ // ``KellyErrorEstimator''
+ // class are zero for each
+ // cell as well, and the
+ // call to
+ // ``refine_and_coarsen_fixed_number''
+ // of the ``triangulation''
+ // object will not flag any
+ // cells for refinement
+ // (why should it if the
+ // indicated error is zero
+ // for each cell?). The
+ // grid in the next
+ // iteration will therefore
+ // consist of four cells
+ // only as well, and the
+ // same problem occurs
+ // again.
+ //
+ // The conclusion needs to
+ // be: while of course we
+ // will not choose the
+ // initial grid to be
+ // well-suited for the
+ // accurate solution of the
+ // problem, we must at
+ // least choose it such
+ // that it has the chance
+ // to capture the most
+ // striking features of the
+ // solution. In this case,
+ // it needs to be able to
+ // see the right hand
+ // side. Thus, we refine
+ // twice globally.
+ triangulation.refine_global (2);
+ }
+ else
+ refine_grid ();
+
+ cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << endl;
+
+ setup_system ();
+
+ cout << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << endl;
+
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ };
+};
+
+
+ // The main function is again exactly
+ // like in step-6 (apart from the
+ // changed class names, of course).
+int main ()
+{
+ try
+ {
+ deallog.depth_console (0);
+
+ ElasticProblem<2> elastic_problem_2d;
+ elastic_problem_2d.run ();
+ }
+ catch (exception &exc)
+ {
+ cerr << endl << endl
+ << "----------------------------------------------------"
+ << endl;
+ cerr << "Exception on processing: " << 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
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+<p>
+In real life, most partial differential equations are really systems
+of equations. Accordingly, the solutions are usually
+vector-valued. The deal.II library supports such problems we will show
+that that is mostly rather simple. The only more complicated problems
+are in assembling matrix and right hand side, but these are easily
+understood as well.
+</p>
+
+<p>
+In the example, we will want to solve the elastic equations. They are
+an extension to Laplace's equation with a vector-valued solution that
+describes the displacement in each space direction of a rigid body
+which is subject to a force. Of course, the force is also
+vector-valued, meaning that in each point it has a direction and an
+absolute value. The elastic equations are the following:
+</p>
+<P ALIGN=CENTER><MATH CLASS=EQUATION>
+ -
+ &pd;<SUB>j</SUB> (c<SUB>ijkl</SUB> &pd;<SUB>k</SUB> u<SUB>l</SUB>)
+ =
+ 0,
+ &quad;
+ i=1&ldots;d,
+</MATH></P>
+<p>
+where the values c<SUB>ijkl</SUB> are the stiffness coefficients and
+will usually depend on the space coordinates. In
+many cases, one knows that the material under consideration is
+isotropic, in which case by introduction of the two coefficients
+lambda and mu the coefficient tensor reduces to
+</p>
+<P ALIGN=CENTER><MATH CLASS=EQUATION>
+ c<SUB>ijkl</SUB>
+ =
+ λδ<SUB>ij</SUB> δ<SUB>kl</SUB> +
+ μ(δ<SUB>ik</SUB> δ<SUB>jl</SUB> +
+ δ<SUB>il</SUB> δ<SUB>jk</SUB>).
+</MATH></P>
+
+<p>
+The elastic equations can then be rewritten in much simpler a form:
+</p>
+<P ALIGN=CENTER><MATH CLASS=EQUATION>
+ -
+ ∇(λ+μ) (∇u)
+ -
+ (∇μ∇) u
+ =
+ 0
+</MATH></P>
+<p>
+This is the form that we will try to solve. The program that does so
+is based on the step-6 program.
+</p>
--- /dev/null
+<a name="Results"></a>
+<h1>Results</h1>
+
+<p>
+There is not much to be said about the results of this program, apart
+from that they look nice. All images were made using GMV from the
+output files that the program wrote to disk. The first picture shows
+the displacement as a vector field, where one vector is shown at each
+vertex of the grid:
+</p>
+
+<p ALIGN=CENTER>
+<IMG SRC="step-8.data/displacement-vectors.jpg" ALT="displacement-vectors" WIDTH="300">
+</p>
+
+<p>
+You can clearly see the sources of x-displacement around x=0.5 and
+x=-0.5, and of y-displacement at the origin. The next image shows the
+final grid after eight steps of refinement:
+</p>
+
+<p ALIGN=CENTER>
+<IMG SRC="step-8.data/final-grid.jpg" ALT="final-grid" WIDTH="300">
+</p>
+
+<p>
+Finally, the x-displacement and y-displacement are displayed separately:
+</p>
+
+<p>
+<TABLE WIDTH="100%">
+<tr>
+<td>
+<IMG SRC="step-8.data/displacement-x.jpg" ALT="displacement-x" WIDTH="300">
+</td>
+<td>
+<IMG SRC="step-8.data/displacement-y.jpg" ALT="displacement-y" WIDTH="300">
+</td>
+</tr>
+</table>
+</p>
+
+<p>
+It should be noted that intuitively one would have expected the
+solution to be symmetric about the x- and y-axes since the x- and
+y-forces are symmetric with respect to these axes. However, the force
+considered as a vector is not symmetric and so not is the solution.
+</p>
+
--- /dev/null
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1999
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly. We get deduce it from the files in the present
+# directory:
+target = $(basename $(shell echo step-*.cc))
+
+# 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
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+
+
+# list of libraries needed to link with
+libs = -ldeal_II_2d -llac -lbase
+libs.g = -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+ libraries = $(target).go $(libs.g)
+ flags = $(CXXFLAGS.g)
+else
+ libraries = $(target).go $(libs)
+ flags = $(CXXFLAGS.o)
+endif
+
+
+
+# make rule for the target. $^ is the object file $(target).g?o
+$(target) : $(libraries)
+ @echo ============================ Linking $@
+ @$(CXX) $(flags) -o $@ $^
+
+# rule how to run the program
+run: $(target)
+ @echo ============================ Running $<
+ @./$(target)
+
+
+# rule to make object files
+%.go : %.cc
+ @echo ============================ Compiling with debugging information: $<
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+ @echo ============================ Compiling with optimization: $<
+ @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+clean:
+ -rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps
+
+
+
+.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.
+#
+# Since the script prefixes the output names by lib/g?o, we have to
+# strip that again (the script was written for the main libraries and
+# large projects where object files are put into subdirs)
+Makefile.dep: $(target).cc Makefile \
+ $(shell echo $D/base/include/base/*.h \
+ $D/lac/include/lac/*.h \
+ $D/deal.II/include/*/*.h)
+ @echo ============================ Remaking Makefile
+ @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \
+ | perl -pi -e 's!lib/g?o/!!g;' \
+ > Makefile.dep
+
+
+include Makefile.dep
+
--- /dev/null
+/* $Id$ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+
+ // As usual, the first few include
+ // files are already known, so we
+ // will not comment on them further.
+#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 <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 <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+ // In this example, we need
+ // vector-valued finite elements. The
+ // suuport for these can be found in
+ // the following include file:
+#include <fe/fe_system.h>
+ // We will compose the vector-valued
+ // finite elements from regular Q1
+ // elements which can be found here,
+ // as usual:
+#include <fe/fe_lib.lagrange.h>
+
+ // This again is C++:
+#include <fstream>
+
+
+ // The main class is, except for its
+ // name, almost unchanged with
+ // respect to the step-6 example. The
+ // only change is the use of a
+ // different class for the ``fe''
+ // variable.
+template <int dim>
+class ElasticProblem
+{
+ public:
+ ElasticProblem ();
+ ~ElasticProblem ();
+ 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;
+ DoFHandler<dim> dof_handler;
+
+ // Instead of a concrete finite
+ // element class such as
+ // ``FEQ1'', we now use a more
+ // generic one, ``FESystem''. In
+ // fact, it is not a finite
+ // element itself, but rather a
+ // class that can be used to
+ // stack several usual elements
+ // together to form one
+ // vector-valued finite
+ // element. In our case, we will
+ // compose the vector-valued
+ // element of ``FEQ1'' objects,
+ // as shown below in the
+ // constructor of this class.
+ FESystem<dim> fe;
+
+ ConstraintMatrix hanging_node_constraints;
+
+ SparseMatrixStruct sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+ // Before going over to the
+ // implementation of the main class,
+ // we declare and define the class
+ // which describes the right hand
+ // side. This time, the right hand
+ // side is vector-valued, as is the
+ // solution, so we will describe the
+ // new elements in some more detail.
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+ public:
+ // The first thing is that
+ // vector-valued functions have a
+ // constructor, since they need
+ // to pass down to the base class
+ // of how many components the
+ // function consists. The default
+ // value in the constructor of
+ // the base class is one, so we
+ // need not define a constructor
+ // for the usual scalar function.
+ RightHandSide ();
+
+ // The next function is a
+ // replacement for the ``value''
+ // function of the previous
+ // examples. There, a second
+ // parameter ``component'' was
+ // given, which denoted which
+ // component was requested. Here,
+ // we implement a function that
+ // returns the whole vector of
+ // values at the given place at
+ // once.
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &values) const;
+
+ // Then, in analogy to the
+ // ``value_list'' function, there
+ // is a function
+ // ``vector_value_list'', which
+ // returns the values of the
+ // vector-valued function at
+ // several points at once:
+ virtual void vector_value_list (const vector<Point<dim> > &points,
+ vector<Vector<double> > &value_list) const;
+};
+
+
+ // This is the constructor of the
+ // right hand side class. As said
+ // above, it only passes down to the
+ // base class the number of
+ // components, which is ``dim'' in
+ // the present case. Note that
+ // although the implementation is
+ // very short here, we do not move it
+ // into the class declaration, since
+ // our style guides require that
+ // inside the class declaration only
+ // declarations have to happen and
+ // that definitions are always to be
+ // found outside.
+template <int dim>
+RightHandSide<dim>::RightHandSide () :
+ Function<dim> (dim)
+{};
+
+
+ // This is the function that returns
+ // the whole vector of values at the
+ // point ``p'' at once:
+template <int dim>
+inline
+void RightHandSide<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+{
+ // To prevent cases where the
+ // return value has not previously
+ // been set to the right size
+ // (which is kind of a convention
+ // in the deal.II library), we test
+ // for this case and otherwise
+ // throw an exception:
+ Assert (values.size() == dim,
+ ExcVectorHasWrongSize (values.size(), dim));
+ // Likewise, if by some accident
+ // someone tried to compile and run
+ // the program in only one space
+ // dimension (in which the elastic
+ // equations do not make much sense
+ // since they reduce to the
+ // ordinary Laplace equation), we
+ // terminate the program if the
+ // dimension is not as expected.
+ Assert (dim >= 2, ExcInternalError());
+
+ // The rest of the function is as
+ // would probably be expected given
+ // the form of the right hand side
+ // function. First we define the
+ // centers of the two points around
+ // which are the sources of
+ // x-displacement, i.e. (0.5,0) and
+ // (-0.5,0). Note that upon
+ // construction of the ``Point''
+ // objects, all components are set
+ // to zero.
+ Point<dim> point_1, point_2;
+ point_1(0) = 0.5;
+ point_2(0) = -0.5;
+
+ // If now the point ``p'' is in the
+ // circle of radius 0.2 around one
+ // of these points, then set the
+ // force in x-direction to one,
+ // otherwise to zero:
+ if (((p-point_1).square() < 0.2*0.2) ||
+ ((p-point_2).square() < 0.2*0.2))
+ values(0) = 1;
+ else
+ values(0) = 0;
+
+ // Likewise, if ``p'' is in the
+ // vicinity of the origin, then set
+ // the y-force to 1, otherwise to
+ // zero:
+ if (p.square() < 0.2*0.2)
+ values(1) = 1;
+ else
+ values(1) = 0;
+};
+
+
+
+ // Now, this is the function of the
+ // right hand side class that returns
+ // the values at several points at
+ // once.
+template <int dim>
+void RightHandSide<dim>::vector_value_list (const vector<Point<dim> > &points,
+ vector<Vector<double> > &value_list) const
+{
+ // First we define an abbreviation
+ // for the number of points which
+ // we shall work on:
+ const unsigned int n_points = points.size();
+
+ // Then we check whether the number
+ // of output slots has been set
+ // correctly, i.e. to the number of
+ // input points:
+ Assert (value_list.size() == n_points,
+ ExcVectorHasWrongSize (value_list.size(), n_points));
+
+ // Finally we treat each of the
+ // points. In one of the previous
+ // examples, we have explained why
+ // the
+ // ``value_list''/``vector_value_list''
+ // function had been introduced: to
+ // prevent us from calling virtual
+ // functions too frequently. On the
+ // other hand, we now need to
+ // implement the same function
+ // twice, which can lead to
+ // confusion if one function is
+ // changed but the other is
+ // not. However, we can prevent
+ // this situation using the
+ // following construct:
+ for (unsigned int p=0; p<n_points; ++p)
+ RightHandSide<dim>::vector_value (points[p],
+ value_list[p]);
+ // It calls the ``vector_value''
+ // function defined above for each
+ // point, and thus preempts all
+ // chances for inconsistency. It is
+ // important to note how the
+ // function was called: using the
+ // full class qualification using
+ // ``RightHandSide::'', since this
+ // calls the function directly and
+ // not using the virtual function
+ // table. The call is thus as fast
+ // as a call to any non-virtual
+ // function. In addition, we have
+ // declared the ``vector_value''
+ // function ``inline'', i.e. the
+ // compiler can remove the function
+ // call altogether and the
+ // resulting code can in principle
+ // be as fast as if we had
+ // duplicated the code.
+};
+
+
+
+
+template <int dim>
+ElasticProblem<dim>::ElasticProblem () :
+ dof_handler (triangulation),
+ // As said before, we
+ // would like to
+ // construct one
+ // vector-valued
+ // finite element as
+ // outer product of
+ // several scala
+ // finite
+ // elements. Of
+ // course, the number
+ // of scalar finite
+ // element we would
+ // like to stack
+ // together equals
+ // the number of
+ // components the
+ // solution function
+ // has, which is
+ // ``dim'' since we
+ // consider
+ // displacement in
+ // each space
+ // direction. The
+ // ``FESystem'' class
+ // can handle this:
+ // we pass it the
+ // finite element of
+ // which we would
+ // like to compose
+ // the system of, and
+ // how often it shall
+ // be repeated:
+ fe (FEQ1<dim>(), dim)
+ // In fact, the ``FESystem'' class
+ // has several more constructors
+ // which can perform more complex
+ // operations that just stacking
+ // together several scalar finite
+ // elements of the same type into
+ // one; we will get to know these
+ // possibilities in later examples.
+ //
+ // It should be noted that the
+ // ``FESystem'' object thus created
+ // does not actually use the finite
+ // element which we have passed to it
+ // as first parameter. We could thus
+ // use an anonymous object created
+ // in-place. The ``FESystem''
+ // constructor only needs the
+ // parameter to deduce the type of
+ // the finite element from this and
+ // then creates objects of the
+ // underlying finite element type
+ // itself.
+{};
+
+
+
+template <int dim>
+ElasticProblem<dim>::~ElasticProblem ()
+{
+ dof_handler.clear ();
+};
+
+
+ // Setting up the system of equations
+ // is equal to the function used in
+ // the step-6 example. The
+ // ``DoFHandler'' class and all other
+ // classes used take care of the
+ // vector-valuedness of the finite
+ // element themselves (in fact, the
+ // do not do so, since they only take
+ // care how many degrees of freedom
+ // there are per vertex, line and
+ // cell, and they do not askwhat they
+ // represent, i.e. whether the finite
+ // element under consideration is
+ // vector-valued or whether it is,
+ // for example, a scalar Hermite
+ // element with several degrees of
+ // freedom on each vertex).
+template <int dim>
+void ElasticProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ // When making the sparsity
+ // pattern, there is some potential
+ // for optimization if not all
+ // components couple to all
+ // others. However, this is not the
+ // case for the elastic equations,
+ // so we use the standard call:
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+ hanging_node_constraints.condense (sparsity_pattern);
+
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+};
+
+
+ // The big changes in this program
+ // are in the creation of matrix and
+ // right hand side, since they are
+ // problem-dependent. We will go
+ // through that process step-by-step,
+ // since it is a bit more complicated
+ // than in previous examples.
+template <int dim>
+void ElasticProblem<dim>::assemble_system ()
+{
+ // First thing: the quadrature
+ // formula does not need
+ // modification since we still deal
+ // with bilinear functions.
+ QGauss2<dim> quadrature_formula;
+ // Also, the ``FEValues'' objects
+ // takes care of everything for us
+ // (or better: it does not really
+ // so; as in the comment in the
+ // function setting up the system,
+ // here as well the ``FEValues''
+ // object computes the same data on
+ // each cell, but it has some
+ // functionality to access data
+ // stored inside the finite element
+ // where they are precomputed upon
+ // construction).
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
+
+ // The number of degrees of freedom
+ // per cell we now obviously ask
+ // from the composed finite element
+ // rather than from the underlying
+ // scalar Q1 element. Here, it is
+ // ``dim'' times the number of
+ // degrees of freedom per cell of
+ // the Q1 element, but this is not
+ // something we need to care about.
+ 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<int> local_dof_indices (dofs_per_cell);
+
+ // As was shown in previous
+ // examples as well, we need a
+ // place where to store the values
+ // of the coefficients at all the
+ // quadrature points on a cell. In
+ // the present situation, we have
+ // two coefficients, lambda and mu.
+ vector<double> lambda_values (n_q_points);
+ vector<double> mu_values (n_q_points);
+
+ // Well, we could as well have
+ // omitted the above two arrays
+ // since we will use constant
+ // coefficients for both lambda and
+ // mu, which can be declared like
+ // this. They both represent
+ // functions always returning the
+ // constant value 1.0. Although we
+ // could omit the respective
+ // factors in the assemblage of the
+ // matrix, we use them here for
+ // purpose of demonstration.
+ ConstantFunction<dim> lambda(1.), mu(1.);
+
+ // Then again, we need to have the
+ // same for the right hand
+ // side. This is exactly as before
+ // in previous examples. However,
+ // we now have a vector-valued
+ // right hand side, which is why
+ // the data type of the
+ // ``rhs_values'' array is
+ // changed. We initialize it by
+ // ``n_q_points'' elements, each of
+ // which is a ``Vector<double>''
+ // with ``dim'' elements.
+ RightHandSide<dim> right_hand_side;
+ vector<Vector<double> > rhs_values (n_q_points,
+ Vector<double>(dim));
+
+
+ // Now we can begin with the loop
+ // over all cells:
+ DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = 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();
+
+ lambda.value_list (q_points, lambda_values);
+ mu.value_list (q_points, mu_values);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int component_i = fe.system_to_component_index(i).first;
+
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const unsigned int component_j = fe.system_to_component_index(j).first;
+
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ {
+ cell_matrix(i,j)
+ +=
+ // (lambda d_i u_i, d_j v_j)
+ (
+ (shape_grads[i][q_point][component_i] *
+ shape_grads[j][q_point][component_j] *
+ lambda_values[q_point])
+ + // (mu d_i u_j, d_i v_j)
+ (shape_grads[i][q_point][component_j] *
+ shape_grads[j][q_point][component_i] *
+ mu_values[q_point])
+ + // (mu d_i v_j, d_i v_j)
+ ((component_i == component_j) ?
+ (shape_grads[i][q_point] *
+ shape_grads[j][q_point] *
+ mu_values[q_point]) :
+ 0)
+ )
+ *
+ JxW_values[q_point];
+ };
+ };
+ };
+
+
+ right_hand_side.vector_value_list (q_points, rhs_values);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int component_i = fe.system_to_component_index(i).first;
+
+ for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
+ cell_rhs(i) += shape_values(i,q_point) *
+ rhs_values[q_point](component_i) *
+ JxW_values[q_point];
+ };
+
+ // The transfer from local
+ // degrees of freedom into the
+ // global matrix and right hand
+ // side vector does not depend
+ // on the equation under
+ // consideration, and is thus
+ // the same as in all previous
+ // examples.
+ 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)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ };
+ };
+
+ // The interpolation of the
+ // boundary values needs a small
+ // modification: since the solution
+ // function is vector-valued, so
+ // needs to be the boundary
+ // values. The ``ZeroFunction''
+ // constructor accepts a parameter
+ // that tells it that it shall
+ // represent a vector valued,
+ // constant zero function with that
+ // many components. By default,
+ // this parameter is equal to one,
+ // in which case the
+ // ``ZeroFunction'' object would
+ // represent a scalar
+ // function. Since the solution
+ // vector has ``dim'' components,
+ // we need to pass ``dim'' as
+ // number of components to the zero
+ // function as well.
+ map<int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(dim),
+ boundary_values);
+ MatrixTools<dim>::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+
+ hanging_node_constraints.condense (system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+};
+
+
+
+ // The solver does not care about
+ // where the system of equations
+ // comes, as long as it stays
+ // positive definite and symmetric
+ // (which are the requirements for
+ // the use of the CG solver), which
+ // the system is. Therefore, we need
+ // not change anything.
+template <int dim>
+void ElasticProblem<dim>::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+
+ PreconditionRelaxation<>
+ preconditioner(system_matrix,
+ &SparseMatrix<double>::template precondition_SSOR<double>,
+ 1.2);
+
+ cg.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+};
+
+
+
+ // The function that does the
+ // refinement of the grid is the same
+ // as in the step-6 example. The
+ // quadrature formula is adapted to
+ // the linear elements again. Note
+ // that the error estimator by
+ // default adds up the estimated
+ // obtained from all components of
+ // the finite element solution, that
+ // is it uses the displacement in all
+ // directions with the same
+ // weight. If we would like the grid
+ // to be adapted to the
+ // x-displacement only, we could pass
+ // the function an additional
+ // parameter which tells it to do so
+ // and do not consider the
+ // displacements in all other
+ // directions for the error
+ // indicators.
+template <int dim>
+void ElasticProblem<dim>::refine_grid ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss2<dim-1>(),
+ neumann_boundary,
+ solution,
+ estimated_error_per_cell);
+
+ triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
+ 0.3, 0.03);
+
+ triangulation.execute_coarsening_and_refinement ();
+};
+
+
+ // The output happens mostly as has
+ // been shown in previous examples
+ // already. The only difference is
+ // not that the solution function is
+ // vector values. The ``DataOut''
+ // class takes care of this
+ // automatically, but we have to give
+ // each component of the solution
+ // vector a different name.
+template <int dim>
+void ElasticProblem<dim>::output_results (const unsigned int cycle) const
+{
+ string filename = "solution-";
+ filename += ('0' + cycle);
+ Assert (cycle < 10, ExcInternalError());
+
+ filename += ".gmv";
+ ofstream output (filename.c_str());
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+
+
+
+ // As said above, we need a
+ // different name for each
+ // component of the solution
+ // function. To pass one name for
+ // each component, a vector of
+ // strings is used. Since the
+ // number of components is the same
+ // as the number of dimensions we
+ // are working in, the following
+ // ``switch'' statement is used.
+ //
+ // We note that some graphics
+ // programs have restriction as to
+ // what characters are allowed in
+ // the names of variables. The
+ // library therefore supports only
+ // the minimal subset of these
+ // characters that is supported by
+ // all programs. Basically, these
+ // are letters, numbers,
+ // underscores, and some other
+ // characters, but in particular no
+ // whitespace and minus/hyphen. The
+ // library will throw an exception
+ // otherwise, at least if in debug
+ // mode.
+ vector<string> solution_names;
+ switch (dim)
+ {
+ case 1:
+ solution_names.push_back ("displacement");
+ break;
+ case 2:
+ solution_names.push_back ("x_displacement");
+ solution_names.push_back ("y_displacement");
+ break;
+ case 3:
+ solution_names.push_back ("x_displacement");
+ solution_names.push_back ("y_displacement");
+ solution_names.push_back ("z_displacement");
+ break;
+ // It is good style to
+ // let the program die if
+ // we run upon a case
+ // which we did not
+ // consider. Remember
+ // that the ``Assert''
+ // macro throws an
+ // exception if the
+ // condition in the first
+ // parameter is not
+ // satisfied. Of course,
+ // the condition
+ // ``false'' can never be
+ // satisfied, so the
+ // program will always
+ // abort whenever it gets
+ // to this statement:
+ default:
+ Assert (false, ExcInternalError());
+ };
+
+ // After setting up the names for
+ // the different components of the
+ // solution vector, we can add the
+ // solution vector to the list of
+ // data vectors scheduled for
+ // output. Note that the following
+ // function takes a vector of
+ // strings as second argument,
+ // whereas the one which we have
+ // used in all previous examples
+ // accepted a string there. In
+ // fact, the latter function is
+ // only a shortcut for the function
+ // which we call here: it puts the
+ // single string that is passed to
+ // it into a vector of strings with
+ // only one element and forwards
+ // that to the other function.
+ data_out.add_data_vector (solution, solution_names);
+ data_out.build_patches ();
+ data_out.write_gmv (output);
+};
+
+
+
+template <int dim>
+void ElasticProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<8; ++cycle)
+ {
+ cout << "Cycle " << cycle << ':' << endl;
+
+ if (cycle == 0)
+ {
+ // As in previous examples,
+ // we use the unit square
+ // (or cube) as domain.
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ // This time, we have to
+ // refine the coarse grid
+ // twice before we first
+ // solve on it. The reason
+ // is the following: we use
+ // the ``Gauss2''
+ // quadrature formula for
+ // integration of the right
+ // hand side; that means
+ // that there are four
+ // quadrature points on
+ // each cell (in 2D). If we
+ // only refine the initial
+ // grid once globally, then
+ // there will be only four
+ // quadrature points in
+ // each direction on the
+ // domain. However, the
+ // right hand side function
+ // was chosen to be rather
+ // localized and in that
+ // case all quadrature
+ // points lie outside the
+ // support of the right
+ // hand side function. The
+ // right hand side vector
+ // will then contain only
+ // zeroes and the solution
+ // of the system of
+ // equations is the zero
+ // vector, i.e. a finite
+ // element function that it
+ // zero everywhere. We
+ // should not be surprised
+ // about such things
+ // happening, since we have
+ // chosen an initial grid
+ // that is totally
+ // unsuitable for the
+ // problem at hand.
+ //
+ // The unfortunate thing is
+ // that if the discrete
+ // solution is constant,
+ // then the error
+ // indicators computed by
+ // the
+ // ``KellyErrorEstimator''
+ // class are zero for each
+ // cell as well, and the
+ // call to
+ // ``refine_and_coarsen_fixed_number''
+ // of the ``triangulation''
+ // object will not flag any
+ // cells for refinement
+ // (why should it if the
+ // indicated error is zero
+ // for each cell?). The
+ // grid in the next
+ // iteration will therefore
+ // consist of four cells
+ // only as well, and the
+ // same problem occurs
+ // again.
+ //
+ // The conclusion needs to
+ // be: while of course we
+ // will not choose the
+ // initial grid to be
+ // well-suited for the
+ // accurate solution of the
+ // problem, we must at
+ // least choose it such
+ // that it has the chance
+ // to capture the most
+ // striking features of the
+ // solution. In this case,
+ // it needs to be able to
+ // see the right hand
+ // side. Thus, we refine
+ // twice globally.
+ triangulation.refine_global (2);
+ }
+ else
+ refine_grid ();
+
+ cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << endl;
+
+ setup_system ();
+
+ cout << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << endl;
+
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ };
+};
+
+
+ // The main function is again exactly
+ // like in step-6 (apart from the
+ // changed class names, of course).
+int main ()
+{
+ try
+ {
+ deallog.depth_console (0);
+
+ ElasticProblem<2> elastic_problem_2d;
+ elastic_problem_2d.run ();
+ }
+ catch (exception &exc)
+ {
+ cerr << endl << endl
+ << "----------------------------------------------------"
+ << endl;
+ cerr << "Exception on processing: " << exc.what() << endl
+ << "Aborting!" << endl
+ << "----------------------------------------------------"
+ << endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ cerr << endl << endl
+ << "----------------------------------------------------"
+ << endl;
+ cerr << "Unknown exception!" << endl
+ << "Aborting!" << endl
+ << "----------------------------------------------------"
+ << endl;
+ return 1;
+ };
+
+ return 0;
+};