774 * <a name=
"step_3-Step3make_grid"></a>
780 * which we would like to do our computation and number each vertex with a
781 * degree of freedom. We have seen these two steps in @ref step_1 "step-1" and @ref step_2 "step-2"
782 * before, respectively.
786 * This function does the first part, creating the mesh. We create the grid
787 * and refine all cells five times. Since the initial grid (which is the
788 * square @f$[-1,1] \times [-1,1]@f$) consists of only one cell, the final grid
789 * has 32 times 32 cells, for a total of 1024.
793 * Unsure that 1024 is the correct number? We can check that by outputting the
794 * number of cells using the <code>n_active_cells()</code> function on the
798 * void Step3::make_grid()
800 * GridGenerator::hyper_cube(triangulation, -1, 1);
801 * triangulation.refine_global(5);
803 * std::cout << "Number of active cells: " << triangulation.n_active_cells()
809 * @note We call the Triangulation::n_active_cells() function, rather than
810 * Triangulation::n_cells(). Here, <i>active</i> means the cells that aren't
827 * <a name=
"step_3-Step3setup_system"></a>
837 *
with each vertex.
While we're at generating output, let us also take a look
838 * at how many degrees of freedom are generated:
841 * void Step3::setup_system()
843 * dof_handler.distribute_dofs(fe);
844 * std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
848 * There should be one DoF for each vertex. Since we have a 32 times 32
849 * grid, the number of DoFs should be 33 times 33, or 1089.
853 * As we have seen in the previous example, we set up a sparsity pattern by
854 * first creating a temporary structure, tagging those entries that might be
855 * nonzero, and then copying the data over to the SparsityPattern object
856 * that can then be used by the system matrix.
859 * DynamicSparsityPattern dsp(dof_handler.n_dofs());
860 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
861 * sparsity_pattern.copy_from(dsp);
865 * Note that the SparsityPattern object does not hold the values of the
866 * matrix, it only stores the places where entries are. The entries
867 * themselves are stored in objects of type SparseMatrix, of which our
868 * variable system_matrix is one.
872 * The distinction between sparsity pattern and matrix was made to allow
873 * several matrices to use the same sparsity pattern. This may not seem
874 * relevant here, but when you consider the size which matrices can have,
875 * and that it may take some time to build the sparsity pattern, this
876 * becomes important in large-scale problems if you have to store several
877 * matrices in your program.
880 * system_matrix.reinit(sparsity_pattern);
884 * The last thing to do in this function is to set the sizes of the right
885 * hand side vector and the solution vector to the right values:
888 * solution.reinit(dof_handler.n_dofs());
889 * system_rhs.reinit(dof_handler.n_dofs());
895 * <a name="step_3-Step3assemble_system"></a>
896 * <h4>Step3::assemble_system</h4>
903 * The next step is to compute the entries of the matrix and right hand side
904 * that form the linear system from which we compute the solution. This is the
905 * central function of each finite element program and we have discussed the
906 * primary steps in the introduction already.
910 * The general approach to assemble matrices and vectors is to loop over all
911 * cells, and on each cell compute the contribution of that cell to the global
912 * matrix and right hand side by quadrature. The point to realize now is that
913 * we need the values of the shape functions at the locations of quadrature
914 * points on the real cell. However, both the finite element shape functions
915 * as well as the quadrature points are only defined on the reference
916 * cell. They are therefore of little help to us, and we will in fact hardly
917 * ever query information about finite element shape functions or quadrature
918 * points from these objects directly.
922 * Rather, what is required is a way to map this data from the reference cell
923 * to the real cell. Classes that can do that are derived from the Mapping
924 * class, though one again often does not have to deal with them directly:
925 * many functions in the library can take a mapping object as argument, but
926 * when it is omitted they simply resort to the standard bilinear Q1
927 * mapping. We will go this route, and not bother with it for the moment (we
928 * come back to this in @ref step_10 "step-10", @ref step_11 "step-11", and @ref step_12 "step-12").
932 * So what we now have is a collection of three classes to deal with: finite
938 * quadrature points
on a
real cell.
950 *
Ok,
let's start: we need a quadrature formula for the evaluation of the
1068 *
const unsigned
int dofs_per_cell = fe.n_dofs_per_cell();
1079 * allocate these objects (these being local objects, all degrees of freedom
1080 * are coupling with all others, and we should use a full matrix object
1081 * rather than a sparse one for the local operations; everything will be
1082 * transferred to a global sparse matrix later on):
1085 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1086 * Vector<double> cell_rhs(dofs_per_cell);
1090 * When assembling the contributions of each cell, we do this with the local
1091 * numbering of the degrees of freedom (i.e. the number running from zero
1092 * through dofs_per_cell-1). However, when we transfer the result into the
1093 * global matrix, we have to know the global numbers of the degrees of
1094 * freedom. When we query them, we need a scratch (temporary) array for
1095 * these numbers (see the discussion at the end of the introduction for
1096 * the type, types::global_dof_index, used here):
1099 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1103 * Now for the loop over all cells. We have seen before how this works for a
1104 * triangulation. A DoFHandler has cell iterators that are exactly analogous
1105 * to those of a Triangulation, but with extra information about the degrees
1114 * examining the cells without modifying them, so it's
good practice to
1118 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1130 *
fe_values.
reinit(cell);
1134 *
Next, reset
the local cell
's contributions to global matrix and
1135 * global right hand side to zero, before we fill them:
1143 * Now it is time to start integration over the cell, which we
1144 * do by looping over all quadrature points, which we will
1145 * number by q_index.
1148 * for (const unsigned int q_index : fe_values.quadrature_point_indices())
1152 * First assemble the matrix: For the Laplace problem, the
1153 * matrix on each cell is the integral over the gradients of
1154 * shape function i and j. Since we do not integrate, but
1155 * rather use quadrature, this is the sum over all
1156 * quadrature points of the integrands times the determinant
1157 * of the Jacobian matrix at the quadrature point times the
1158 * weight of this quadrature point. You can get the gradient
1159 * of shape function @f$i@f$ at quadrature point with number q_index by
1160 * using <code>fe_values.shape_grad(i,q_index)</code>; this
1161 * gradient is a 2-dimensional vector (in fact it is of type
1162 * Tensor@<1,dim@>, with here dim=2) and the product of two
1163 * such vectors is the scalar product, i.e. the product of
1164 * the two shape_grad function calls is the dot
1165 * product. This is in turn multiplied by the Jacobian
1166 * determinant and the quadrature point weight (that one
1167 * gets together by the call to FEValues::JxW() ). Finally,
1168 * this is repeated for all shape functions @f$i@f$ and @f$j@f$:
1171 * for (const unsigned int i : fe_values.dof_indices())
1172 * for (const unsigned int j : fe_values.dof_indices())
1173 * cell_matrix(i, j) +=
1174 * (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
1175 * fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
1176 * fe_values.JxW(q_index)); // dx
1180 * We then do the same thing for the right hand side. Here,
1181 * the integral is over the shape function i times the right
1182 * hand side function, which we choose to be the function
1183 * with constant value one (more interesting examples will
1184 * be considered in the following programs).
1187 * for (const unsigned int i : fe_values.dof_indices())
1188 * cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1190 * fe_values.JxW(q_index)); // dx
1194 * Now that we have the contribution of this cell, we have to transfer
1195 * it to the global matrix and right hand side. To this end, we first
1196 * have to find out which global numbers the degrees of freedom on this
1200 *
cell->get_dof_indices(local_dof_indices);
1206 *
obtained using local_dof_indices[i]:
1209 *
for (
const unsigned int i : fe_values.dof_indices())
1211 *
system_matrix.add(local_dof_indices[i],
1212 *
local_dof_indices[
j],
1220 *
for (
const unsigned int i : fe_values.dof_indices())
1229 * Laplace
's equation without Dirichlet boundary values is not even uniquely
1230 * solvable, since you can add an arbitrary constant to the discrete
1231 * solution). We therefore have to do something about the situation.
1235 * For this, we first obtain a list of the degrees of freedom on the
1236 * boundary and the value the shape function shall have there. For
1237 * simplicity, we only interpolate the boundary value function, rather than
1238 * projecting it onto the boundary. There is a function in the library which
1239 * does exactly this: VectorTools::interpolate_boundary_values(). Its
1240 * parameters are (omitting parameters for which default values exist and
1282 *
VectorTools::interpolate_boundary_values(dof_handler,
1283 *
types::boundary_id(0),
1329 * yet, so we tell it to use the identity operation as preconditioner. Later
1330 * tutorial programs will spend significant amount of time and space on
1331 * constructing better preconditioners.
1335 * At the end of this process, the `solution` variable contains the
1336 * nodal values of the solution function. At the end of the function, we
1337 * output how many Conjugate Gradients iterations it took to solve the
1341 * void Step3::solve()
1343 * SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
1344 * SolverCG<Vector<double>> solver(solver_control);
1345 * solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
1347 * std::cout << solver_control.last_step()
1348 * << " CG iterations needed to obtain convergence." << std::endl;
1355 * <a name="step_3-Step3output_results"></a>
1356 * <h4>Step3::output_results</h4>
1360 * The last part of a typical finite element program is to output the results
1361 * and maybe do some postprocessing (for example compute the maximal stress
1362 * values at the boundary, or the average flux across the outflow, etc). We
1363 * have no such postprocessing here, but we would like to write the solution
1367 * void Step3::output_results() const
1371 * To write the output to a file, we need an object which knows about output
1372 * formats and the like. This is the DataOut class, and we need an object of
1376 * DataOut<2> data_out;
1379 * Now we have to tell it where to take the values from which it shall
1380 * write. We tell it which DoFHandler object to use, and the solution vector
1381 * (and the name by which the solution variable shall appear in the output
1382 * file). If we had more than one vector which we would like to look at in
1383 * the output (for example right hand sides, errors per cell, etc) we would
1387 * data_out.attach_dof_handler(dof_handler);
1388 * data_out.add_data_vector(solution, "solution");
1391 * After the DataOut object knows which data it is to work on, we have to
1392 * tell it to process them into something the back ends can handle. The
1393 * reason is that we have separated the frontend (which knows about how to
1394 * treat DoFHandler objects and data vectors) from the back end (which knows
1395 * many different output formats) and use an intermediate data format to
1396 * transfer data from the front- to the backend. The data is transformed
1397 * into this intermediate format by the following function:
1400 * data_out.build_patches();
1404 * Now we have everything in place for the actual output. Just open a file
1405 * and write the data into it, using VTK format (there are many other
1406 * functions in the DataOut class we are using here that can write the
1407 * data in postscript, AVS, GMV, Gnuplot, or some other file
1411 * const std::string filename = "solution.vtk";
1412 * std::ofstream output(filename);
1413 * data_out.write_vtk(output);
1414 * std::cout << "Output written to " << filename << std::endl;
1421 * <a name="step_3-Step3run"></a>
1422 * <h4>Step3::run</h4>
1426 * Finally, the last function of this class is the main function which calls
1427 * all the other functions of the <code>Step3</code> class. The order in which
1428 * this is done resembles the order in which most finite element programs
1429 * work. Since the names are mostly self-explanatory, there is not much to
1437 * assemble_system();
1446 * <a name="step_3-Thecodemaincodefunction"></a>
1447 * <h3>The <code>main</code> function</h3>
1451 * This is the main function of the program. Since the concept of a
1452 * main function is mostly a remnant from the pre-object oriented era
1453 * before C++ programming, it often does not do much more than
1454 * creating an object of the top-level class and calling its principle
1460 * Step3 laplace_problem;
1461 * laplace_problem.run();
1466<a name="step_3-Results"></a><h1>Results</h1>
1469The output of the program looks as follows:
1471Number of active cells: 1024
1472Number of degrees of freedom: 1089
147336 CG iterations needed to obtain convergence.
1474Output written to solution.vtk
1477The last line is the output we generated at the bottom of the
1478`output_results()` function: The program generated the file
1479<code>solution.vtk</code>, which is in the VTK format that is widely
1480used by many visualization programs today -- including the two
1481heavy-weights <a href="https://www.llnl.gov/visit">VisIt</a> and
1482<a href="https://www.paraview.org">Paraview</a> that are the most
1483commonly used programs for this purpose today.
1485Using VisIt, it is not very difficult to generate a picture of the
1487<table width="60%" align="center">
1490 <img src="https://www.dealii.org/images/steps/developer/step-3.solution-3.png" alt="Visualization of the solution of step-3">
1494It shows both the solution and the mesh, elevated above the @f$x@f$-@f$y@f$ plane
1495based on the value of the solution at each point. Of course the solution
1496here is not particularly exciting, but that is a result of both what the
1497Laplace equation represents and the right hand side @f$f(\mathbf x)=1@f$ we
1498have chosen for this program: The Laplace equation describes (among many
1499other uses) the vertical deformation of a membrane subject to an external
1500(also vertical) force. In the current example, the membrane's
borders
1507See
also <a
href=
"https://www.math.colostate.edu/~bangerth/videos.676.11.html">
video lecture 11</a>, <a
href=
"https://www.math.colostate.edu/~bangerth/videos.676.32.html">
video lecture 32</a>.
1511<a name=
"step-3-extensions"></a>
1512<a name=
"step_3-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
1531 non-zero constant boundary values
using
1544 interpolate boundary values
to zero
on all boundary components
with
1582 if (face->at_boundary())
1583 if (
std::fabs(face->center()(1) - (-1.0)) < 1e-12 ||
1584 std::fabs(face->center()(1) - (1.0)) < 1e-12)
1585 face->set_boundary_id(1);
1596 VectorTools::interpolate_boundary_values(dof_handler,
1597 types::boundary_id(1),
1636--------------------------------------------------------
1647 the active_fe_index?
1661 have to add `
#include <deal.II/fe/fe_simplex_p.h>` at the top of the file.)
1670 <
img src=
"https://www.dealii.org/images/steps/developer/step-3.solution-triangles.png" alt=
"Visualization of the solution of step-3 using triangles">
1682 std::cout <<
"Solution at (1/3,1/3): "
1689 <table
align=
"center" class=
"doxtable">
1721 QGauss<2>(fe.degree + 1),
1730 <table
align=
"center" class=
"doxtable">
1765n_refinement_steps = 5;
1786const std::string
filename_h5 =
"solution_" + std::to_string(n_refinement_steps) +
".h5";
1810#include <deal.II/base/hdf5.h>
1821data_file.write_dataset(
"point_value", point_value);
1826data_file.write_dataset(
"mean_value",mean_value);
1838<a
href=
"https://en.wikipedia.org/wiki/R_(programming_language)">R
1842<a
href=
"https://datacarpentry.org/R-ecology-lesson/index.html">
here</a>.
1893# Counting in R starts at 1 instead of 0, so we need to increment all
1894# vertex indices by one:
1897# Store the x and y positions of each vertex in one big vector in a
1898# cell by cell fashion (every 4 entries belong to one cell):
1902# Construct a vector that stores the matching cell by cell grouping
1903# (1,1,1,1,2,2,2,2,...):
1925<table width=
"60%" align=
"center">
1928 <
img src=
"https://www.dealii.org/images/steps/developer/step-3.extensions.grid_5.png" alt=
"Grid after 5 refinement steps of step-3">
1952<table width=
"60%" align=
"center">
1955 <
img src=
"https://www.dealii.org/images/steps/developer/step-3.extensions.pseudocolor_5.png" alt=
"Solution after 5 refinement steps of step-3">
1985# Calculate the error w.r.t. our maximum refinement step
1990dofs <- dofs[1:
n_ref-1]
1997pdf (
paste(
"convergence.pdf",sep=
""),width = 5,height = 4.2)
2017 <
td><
img src=
"https://www.dealii.org/images/steps/developer/step-3.extensions.convergence_mean.png" alt=
""></
td>
2018 <
td><
img src=
"https://www.dealii.org/images/steps/developer/step-3.extensions.convergence_point.png" alt=
""></
td>
2052cells : <
HDF5 dataset "cells": shape (1024, 4), type
"<u4">
2053mean_value : <
HDF5 dataset "mean_value": shape (1,), type
"<f8">
2056solution : <
HDF5 dataset "solution": shape (1089, 1), type
"<f8">
2065nodes =
np.array(file[
"/nodes"])
2066cells =
np.array(file[
"/cells"])
2067solution =
np.array(file[
"/solution"])
20954091 0.93750 1.00000 1022
20964092 0.96875 0.96875 1023
20974093 1.00000 0.96875 1023
20984094 1.00000 1.00000 1023
20994095 0.96875 1.00000 1023
21014096 rows
× 3 columns
2107ax.set_aspect(
"equal",
"box")
2108ax.set_title(
"grid at refinement level #%d" %
refinement)
2110for cell_id, cell in
meshdata.groupby([
"ids"]):
2111 cell =
pd.concat([cell, cell.head(1)])
2112 plt.plot(cell[
"x"], cell[
"y"], c=
"k")
2117ax.set_aspect(
"equal",
"box")
2118ax.set_title(
"grid at refinement level #%d" %
refinement)
2119for cell_id, cell in
meshdata.groupby([
"ids"]):
2129xg =
np.linspace(
x.min(),
x.max(),
nx)
2130yg =
np.linspace(
y.min(),
y.max(),
nx)
2138c = ax.pcolor(xgrid, ygrid, solution_grid, cmap=cm.viridis)
2139fig.colorbar(c, ax=ax)
2144To check the convergence of `mean_value` and `point_value`
2145we loop over data of all refinements and store into vectors <code>mean_values</code> and <code>mean_values</code>
2147mean_values = np.zeros((8,))
2148point_values = np.zeros((8,))
2149dofs = np.zeros((8,))
2151for refinement in range(1, 9):
2152 filename = "solution_%d.h5" % refinement
2153 file = h5py.File(filename, "r")
2154 mean_values[refinement - 1] = np.array(file["/mean_value"])[0]
2155 point_values[refinement - 1] = np.array(file["/point_value"])[0]
2156 dofs[refinement - 1] = np.array(file["/solution"]).shape[0]
2159Following is the plot of <code>mean_values</code> on `log-log` scale
2161mean_error = np.abs(mean_values[1:] - mean_values[:-1])
2162plt.loglog(dofs[:-1], mean_error)
2165plt.ylabel("mean value error")
2169Following is the plot of <code>point_values</code> on `log-log` scale
2171point_error = np.abs(point_values[1:] - point_values[:-1])
2172plt.loglog(dofs[:-1], point_error)
2175plt.ylabel("point value error")
2179A Python package which mimics the `R` ggplot2 (which is based on specifying the grammar of the graphics) is
2183from plotnine import (
2195Then plot the mesh <code>meshdata</code> dataframe
2198 ggplot(meshdata, aes(x="x", y="y", group="ids"))
2199 + geom_polygon(fill="white", colour="black")
2200 + ggtitle("grid at refinement level #%d" % refinement)
2205Collect the solution into a dataframe
2207colordata = pd.DataFrame({"x": x, "y": y, "solution": solution.flatten()})
2212 ggplot(colordata, aes(x="x", y="y", fill="solution"))
2213 + geom_raster(interpolate=True)
2214 + ggtitle("solution at refinement level #%d" % refinement)
2220Collect the convergence data into a dataframe
2222convdata = pd.DataFrame(
2223 {"dofs": dofs[:-1], "mean_value": mean_error, "point_value": point_error}
2227Following is the plot of <code>mean_values</code> on `log-log` scale
2230 ggplot(convdata, mapping=aes(x="dofs", y="mean_value"))
2232 + labs(x="#DoFs", y="mean value error")
2237plot.save("mean_error.pdf", dpi=60)
2241Following is the plot of <code>point_values</code> on `log-log` scale
2244 ggplot(convdata, mapping=aes(x="dofs", y="point_value"))
2246 + labs(x="#DoFs", y="point value error")
2251plot.save("point_error.pdf", dpi=60)
2256<a name="step_3-PlainProg"></a>
2257<h1> The plain program</h1>
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
void write_dataset(const std::string &name, const Container &data) const
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
std::vector< index_type > data
void convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim > > &vertices)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
@ matrix
Contents is actually a matrix.
@ general
No special properties.
constexpr types::blas_int zero
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation