From: wolf Date: Sun, 9 Jan 2000 17:18:02 +0000 (+0000) Subject: Add step-8. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f54c66a1a75e81f58273592a7636120a8ec5e9b7;p=dealii-svn.git Add step-8. git-svn-id: https://svn.dealii.org/trunk@2174 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-8/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-8/Makefile new file mode 100644 index 0000000000..8753051799 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-8/Makefile @@ -0,0 +1,94 @@ +# $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 + diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-8/step-8.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-8/step-8.cc new file mode 100644 index 0000000000..9fac1af41e --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-8/step-8.cc @@ -0,0 +1,961 @@ +/* $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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + // In this example, we need + // vector-valued finite elements. The + // suuport for these can be found in + // the following include file: +#include + // We will compose the vector-valued + // finite elements from regular Q1 + // elements which can be found here, + // as usual: +#include + + // This again is C++: +#include + + + // 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 +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 triangulation; + DoFHandler 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 fe; + + ConstraintMatrix hanging_node_constraints; + + SparseMatrixStruct sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector 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 +class RightHandSide : public Function +{ + 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 &p, + Vector &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 > &points, + vector > &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 +RightHandSide::RightHandSide () : + Function (dim) +{}; + + + // This is the function that returns + // the whole vector of values at the + // point ``p'' at once: +template +inline +void RightHandSide::vector_value (const Point &p, + Vector &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 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 +void RightHandSide::vector_value_list (const vector > &points, + vector > &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::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 +ElasticProblem::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) + // 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 +ElasticProblem::~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 +void ElasticProblem::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 +void ElasticProblem::assemble_system () +{ + // First thing: the quadrature + // formula does not need + // modification since we still deal + // with bilinear functions. + QGauss2 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 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + vector 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 lambda_values (n_q_points); + vector 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 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'' + // with ``dim'' elements. + RightHandSide right_hand_side; + vector > rhs_values (n_q_points, + Vector(dim)); + + + // Now we can begin with the loop + // over all cells: + DoFHandler::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 + & shape_values = fe_values.get_shape_values(); + const vector > > + & shape_grads = fe_values.get_shape_grads(); + const vector + & JxW_values = fe_values.get_JxW_values(); + const vector > + & 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; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(dim), + boundary_values); + MatrixTools::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 +void ElasticProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionRelaxation<> + preconditioner(system_matrix, + &SparseMatrix::template precondition_SSOR, + 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 +void ElasticProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::FunctionMap neumann_boundary; + KellyErrorEstimator::estimate (dof_handler, + QGauss2(), + 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 +void ElasticProblem::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 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 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 +void ElasticProblem::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; +}; diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-vectors.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-vectors.jpg new file mode 100644 index 0000000000..4078b9fbb4 Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-vectors.jpg differ diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-x.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-x.jpg new file mode 100644 index 0000000000..306bc68dde Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-x.jpg differ diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-y.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-y.jpg new file mode 100644 index 0000000000..3028b9f6e9 Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/displacement-y.jpg differ diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/final-grid.jpg b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/final-grid.jpg new file mode 100644 index 0000000000..1d1fe99dd8 Binary files /dev/null and b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.data/final-grid.jpg differ diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.intro b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.intro new file mode 100644 index 0000000000..6fefefb616 --- /dev/null +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.intro @@ -0,0 +1,58 @@ + +

Introduction

+ +

+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. +

+ +

+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: +

+

+ - + &pd;j (cijkl &pd;k ul) + = + 0, + &quad; + i=1&ldots;d, +

+

+where the values cijkl 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 +

+

+ cijkl + = + λδij δkl + + μ(δik δjl + + δil δjk). +

+ +

+The elastic equations can then be rewritten in much simpler a form: +

+

+ - + ∇(λ+μ) (∇u) + - + (∇μ∇) u + = + 0 +

+

+This is the form that we will try to solve. The program that does so +is based on the step-6 program. +

diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.results b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.results new file mode 100644 index 0000000000..a0a2f7357f --- /dev/null +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-8.results @@ -0,0 +1,49 @@ + +

Results

+ +

+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: +

+ +

+displacement-vectors +

+ +

+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: +

+ +

+final-grid +

+ +

+Finally, the x-displacement and y-displacement are displayed separately: +

+ +

+ + + + + +
+displacement-x + +displacement-y +
+

+ +

+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. +

+ diff --git a/deal.II/examples/step-8/Makefile b/deal.II/examples/step-8/Makefile new file mode 100644 index 0000000000..8753051799 --- /dev/null +++ b/deal.II/examples/step-8/Makefile @@ -0,0 +1,94 @@ +# $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 + diff --git a/deal.II/examples/step-8/step-8.cc b/deal.II/examples/step-8/step-8.cc new file mode 100644 index 0000000000..9fac1af41e --- /dev/null +++ b/deal.II/examples/step-8/step-8.cc @@ -0,0 +1,961 @@ +/* $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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + // In this example, we need + // vector-valued finite elements. The + // suuport for these can be found in + // the following include file: +#include + // We will compose the vector-valued + // finite elements from regular Q1 + // elements which can be found here, + // as usual: +#include + + // This again is C++: +#include + + + // 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 +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 triangulation; + DoFHandler 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 fe; + + ConstraintMatrix hanging_node_constraints; + + SparseMatrixStruct sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector 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 +class RightHandSide : public Function +{ + 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 &p, + Vector &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 > &points, + vector > &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 +RightHandSide::RightHandSide () : + Function (dim) +{}; + + + // This is the function that returns + // the whole vector of values at the + // point ``p'' at once: +template +inline +void RightHandSide::vector_value (const Point &p, + Vector &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 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 +void RightHandSide::vector_value_list (const vector > &points, + vector > &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::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 +ElasticProblem::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) + // 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 +ElasticProblem::~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 +void ElasticProblem::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 +void ElasticProblem::assemble_system () +{ + // First thing: the quadrature + // formula does not need + // modification since we still deal + // with bilinear functions. + QGauss2 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 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + vector 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 lambda_values (n_q_points); + vector 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 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'' + // with ``dim'' elements. + RightHandSide right_hand_side; + vector > rhs_values (n_q_points, + Vector(dim)); + + + // Now we can begin with the loop + // over all cells: + DoFHandler::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 + & shape_values = fe_values.get_shape_values(); + const vector > > + & shape_grads = fe_values.get_shape_grads(); + const vector + & JxW_values = fe_values.get_JxW_values(); + const vector > + & 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; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(dim), + boundary_values); + MatrixTools::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 +void ElasticProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionRelaxation<> + preconditioner(system_matrix, + &SparseMatrix::template precondition_SSOR, + 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 +void ElasticProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::FunctionMap neumann_boundary; + KellyErrorEstimator::estimate (dof_handler, + QGauss2(), + 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 +void ElasticProblem::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 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 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 +void ElasticProblem::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; +};